我目前正在尝试在用Fortran90(以及一点Fortran03(编写的代码中使用PRIMME库。PRIMME本身是用C编写的,但它附带了一个Fortran77接口,这应该会使从Fortran中使用它相对容易。不过,我遇到了一个与将不同的Fortran风格结合在一起有关的问题,我还没有找到一个令人信服的答案:
PRIMME是一个迭代本征解算器库,用于计算大型对称/厄米矩阵A的本征向量和本征值,而无需将矩阵本身存储在内存中。矩阵A仅用于在PRIMME库之外实现的矩阵矩阵乘法。然后,矩阵矩阵乘法例程的函数指针被传递到PRIMME库,PRIMME将在内部多次调用它来计算特征向量和特征值。
看看库附带的Fortran77绑定示例,矩阵矩阵乘法例程应该具有以下签名。
subroutine MatMatMul(B,C,k,primme)
real*8 B(*), C(*)
integer k, primme
! calculate C = A * B
end
这里k
是矩阵B和C的列数。primme
是作为Fortran整数传递的,但它实际上是指向数据结构的C指针,该数据结构包含PRIMME库的配置,最显著的是特征值问题的维度,即B和C的行数。因此,在MatMatMul
方法中,我有关于所涉及矩阵大小的所有信息。
不幸的是,矩阵B和C被作为一维假定大小数组B(*)
和C(*)
传递到MatMatMul
中。虽然我当然可以自己完成所有的索引运算,但我更喜欢将B
和C
作为二维数组访问。此外,在它们上使用二维切片也很好。Fortran代码库的其余部分是使用模块/显式接口和假定形状数组编写的,因此这也将与代码的其余部分更加一致。
是否有任何方法可以将B
和C
视为具有运行时确定大小的二维Fortan数组
考虑到我知道运行时的行数和列数,我认为这应该不是什么大问题。我在网上搜索了一下,但没有找到解决方案(实际上也没有任何关于这个问题的讨论(。我读过很多次,你可以对这些假设大小的数组做的事情非常有限,因为编译器没有足够的信息来进行切片。我真的不理解这一点,因为我作为假定形状数组传递的所有其他可分配数组的形状都是在运行时确定的,所以编译器无论如何都没有关于它们的信息,但仍然允许我对它们进行切片。鉴于这一事实,将CCD_ 11视为CCD_。
我还与我的同事们进行了交谈,他们对Fortran有更多的经验,虽然他们提出了一些解决方案,例如将MatMatMul
作为另一个方法的包装器,该方法接受B
和C
作为显式形状的二维数组,并将新方法放在模块之外(这样编译器就不会检查接口(,但他们没有我所希望的简单解决方案。
我对Fortran还没有太多的经验,所以我想我应该在这里询问一下,以确保我没有遗漏什么。
编辑:根据下面接受的答案,我选择了以下解决方案:
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
然后将指向包装器的指针传递到PRIMME中。它不太漂亮,但效果很好。
由于sequence association
的规则,允许传递二维数组。数组在子程序中按内存存储顺序(主列(重新解释为一维数组。
请参阅Fortran标准,或了解更多非正式介绍http://michaelgoerz.net/blog/2011/05/advanced-array-passing-in-fortran/
唯一存在问题的地方是生成通用接口,在通用消歧中严格执行等级。
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
甚至
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
但我更喜欢使用Fortran 90的内在函数MATMUL。或者BLAS子例程dgemv
,如果您需要最高性能,但有些编译器会自动从matmul
调用它。