Вопрос

I am currently trying to use the PRIMME library in a code written in Fortran90 (and a bit of Fortran03). PRIMME itself is written in C, but it comes with a Fortran77 interface that should make using it from Fortran relatively easy. I have run into a question though that is related to bringing the different Fortran styles together, to which I have not found a convincing answer yet:

PRIMME is an iterative eigensolver library that calculates the eigenvectors and eigenvalues of huge, symmetric/hermitian matrix A without ever storing the matrix itself in memory. The matrix A is only used in a matrix-matrix multiplication that is implemented outside of the PRIMME library. A function pointer to the matrix-matrix multiplication routine is then passed into the PRIMME library which will then internally call it a number of times to calculate eigenvectors and eigenvalues.

Looking at the Fortran77 binding example that comes with the library, the matrix-matrix multiplication routine should have the following signature.

subroutine MatMatMul(B,C,k,primme)

    real*8 B(*), C(*)
    integer k, primme

    ! calculate C = A * B
end

Here k is the number of columns of the matrices B and C. primme is passed as a Fortran integer, but it is actually a C pointer to a datastructure that contains the configuration of the PRIMME library, most notably the dimension of the eigenvalue problem, which is the number of rows of B and C. So within the MatMatMul method I have all the information I need about the sizes of the involved matrices.

Unfortunately the matrices B and C are passed into MatMatMul as one-dimensional assumed size arrays B(*) and C(*). While I could of course do all the index arithmetic myself, I would prefer to access B and C as two-dimensional arrays. Also it would be nice to use two-dimensional slicing on them. The rest of the Fortran codebase is written using modules/explicit interfaces and assumed shape arrays, so this would also be more consistent with the rest of the code.

Is there any way to treat B and C as two-dimensional Fortan arrays with runtime determined sizes?

Given the fact that I know the number of rows and columns at runtime, I thought it should not be a big problem. I searched the internet but did not find a solution (or any discussion of the problem actually). I've read a lot of times that what you can do with these assumed size arrays is very limited, as the compiler does not have enough information about them to do for example slicing. I don't really understand this point, as the shape of all my other allocatable arrays that I pass around as assumed shape arrays is determined at runtime, so the compiler has no information about them anyway but still allows me to do slicing with them. In view of this fact there should be absolutely no problem in treating B(*) as B(n,k).

I also talked to my collegues who have more experience with Fortran, and while they suggested some workarounds, e.g. making MatMatMul a wrapper to another method that accepts B and C as explicit shape two-dimensional arrays and putting the new method outside of a module (so that the interface is not checked by the compiler), they did not have the simple solution I was hoping for.

I do not have too much experience with Fortran yet, so I thought I should ask here to make sure I'm not missing something.

EDIT: Based on the accepted answer below, I went for the following solution:

subroutine WrapperMatMatMul(B,C,k,primme)

    real*8 B(*), C(*)
    integer k, primme

    call primme_get_member_f77(primme, PRIMMEF77_n, n)
    call MatMatMul(B,C,k,n,primme)
end

subroutine MatMatMul(B,C,k,n,primme)

    real*8 B(k,n), C(k,n)
    integer k, n, primme

    ! calculate C = A * B
end

The pointer to the wrapper is then passed into PRIMME. It's not exactly pretty, but works perfectly.

Это было полезно?

Решение

It is allowed to pass two dimensional arrays thanks to the rules of sequence association. The arrays are re-interpreted as on dimensional in the subroutine in the memory storage order (column major).

See the Fortran standard, or for more informal introduction http://michaelgoerz.net/blog/2011/05/advanced-array-passing-in-fortran/

The only place where it is a problem is generating generic interfaces, where the ranks are strictly enforced in generic disambiguation.

subroutine MatMatMul(B,C,k,primme)

    real*8 B(*), C(*)
    integer k, primme

    ! calculate C = A * B
end

real*8 A(10,20), B(20,10)

call MatMatMul(A,B,20,0)

end

or even

real*8 A(10,20), B(20,10)

call MatMatMul(A,B,20,0)

contains

subroutine MatMatMul(B,C,k,primme)

    real*8 B(*), C(*)
    integer k, primme

    ! calculate C = A * B
end

end

But I would prefer to use the Fortran 90 intrinsic function MATMUL instead. Or BLAS subroutine dgemv, if you need top performance, but some compilers call it from matmul automatically.

Лицензировано под: CC-BY-SA с атрибуция
Не связан с StackOverflow
scroll top