On 3/12/22 4:34 AM, Neil wrote:
I think I probably need to roll up my sleeves and replace the array
elements with arrays in calls to BLAS & LAPACK routines. But I don't
see why in principle compilers can't cope with overloading of
subprograms with assumed-size array arguments, allowing either array
segments or initial array elements to be given in the calls to the
overloaded subprograms.
I did not follow exactly what your code examples do, but in general you
DO NOT want to do what you are attempting to do when performance is
important. The reason is that when you pass an array section actual
argument to an assumed size dummy argument (or to a subprogram with an
implicit interface), the compiler will invoke copy-in/copy-out
semantics. That is because the assumed size array must be contiguous.
The compiler might look ahead and skip making that copy in some cases
when it knows the actual argument is also contiguous, but in other cases
it is not possible for the compiler to know, so it will add the overhead anyway. Or at least it will add overhead to do the checking. For a
vector operation (level-1 BLAS, like ddot()), there are N memory fetches
for each argument and a comparable number of floating point operations,
so the performance is already memory-limited. If you add an extra N
memory fetches and an extra N memory stores to do the copy-in for each argument, then you have tripled the memory overhead. Then if it needs to
do the copy-out (some level-1 BLAS will need this, like daxpy()), there
will be an extra 2*N memory overhead for each such vector, multiplying
the memory overhead further.
Furthermore, there will be the odd looking feature that if you pass in noncontiguous array slices, the incx and incy values will be 1, while in
the original ddot() call they would have been >1. You should certainly
leave comments in the code explaining this quirk to other programmers
(or to future-you a month later).
This confusing feature of actual and dummy argument association is not
limited to BLAS and LAPACK, it appears everywhere with legacy f77
libraries. MPI calls are another place where the programmer must be
careful, for example. The answer to this problem is that all of these
libraries must be rewritten to use f90+ calling conventions, from top to bottom, in order to eliminate these redundant copy-in/copy-out
operations. It cannot be done simply with interface blocks. The reason
is that somewhere down the calling tree there will be the possibility of
a noncontiguous array actual argument that is associated with a
contiguous assumed size dummy argument, and the compiler must make the temporary copy for that to work.
You cannot use tricks like defining the interface block dummy arguments
to be scalars rather than vectors. That might work in some cases,
despite violating the standard, but eventually you will find a compiler,
or an optimization level of a compiler, where that fails.
It does not look like that major rewrite is ever going to happen. At
least for BLAS, LAPACK, and MPI, it has been three decades and the code
is still written and maintained using the old conventions. So I think
the answer is to learn how storage association works (yes, it is
complicated), and to continue to use the old f77 conventions for those
types of libraries when performance is important.
$.02 -Ron Shepard
--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)