在 Mac OS X 上的 Fortran 中链接 LAPACK

Question

我想这是一个标准的菜鸟问题，但是在花了一上午的时间搜索网络之后，我还是决定给你找麻烦。我在 Mac OS 10.9 上，我想从 Fortran 程序中调用 LAPACK 特征值例程。昨天我很高兴被介绍给 Fortran，所以请原谅任何愚蠢的错误。

这是我想要运行的最小示例：

program eigtest
    complex A(3,3)
    real eigs(3)
    A(1,1) = cmplx(1,0)
    A(1,2) = cmplx(0,2)
    A(1,3) = cmplx(3,0)
    A(2,1) = cmplx(0,-2)
    A(2,2) = cmplx(5,0)
    A(2,3) = cmplx(1,-1)
    A(3,1) = cmplx(3,0)
    A(3,2) = cmplx(1,1)
    A(3,3) = cmplx(7,0)
    call heevd(A, eigs)
    write(*,*) eigs
end

我了解到在 OS X 上，LAPACK 是 Accelerate 框架的一部分，所以我尝试了以下方法：

gfortran -o eigtest -framework accelerate eigtest.f95

但链接器抱怨：

Undefined symbols for architecture x86_64:
  "_heevd_", referenced from:
      _MAIN__ in ccleuVFO.o
ld: symbol(s) not found for architecture x86_64
collect2: ld returned 1 exit status

Answer 1

我不确定 Accelerate 框架是否对 heevd(A, eigs) 进行了很好的 fortran95 调用。下面，我总结了 zheevd lapack 子例程（以及它的所有工作区变量），以使调用保持良好和紧凑。您可以将这样的丑陋存储在可以加载到程序中的外部模块中。我认为intel MKL拥有大部分lapack95接口，os x框架可能还有其他更好的方法。下面的代码编译为：

gfortran -o eigtest -framework Accelerate eigtest.f95    

program eigtest
    complex(kind=8),allocatable :: A(:,:), eigs(:), vecs(:,:)
    integer                     :: ierr

    allocate(A(3,3),stat=ierr)
    if (ierr /= 0) STOP "*** not enough memory ***"

    A(1,1) = cmplx(1,0)
    A(1,2) = cmplx(0,2)
    A(1,3) = cmplx(3,0)
    A(2,1) = cmplx(0,-2)
    A(2,2) = cmplx(5,0)
    A(2,3) = cmplx(1,-1)
    A(3,1) = cmplx(3,0)
    A(3,2) = cmplx(1,1)
    A(3,3) = cmplx(7,0)
    !call heevd(A, eigs)
    call wrapped_zheevd(A, eigs,vecs)
    write(*,*) eigs

    contains

subroutine wrapped_zheevd (matin, &
                         zvals,zvecs )
    integer                                  :: ndim
    complex(kind=8),intent(in),  allocatable :: matin(:,:)
    complex(kind=8),intent(out), allocatable :: zvals(:),zvecs(:,:)
    character*1                              :: jobz='V',uplo='U'
    integer                                  :: info,lda,liwork,lrwork,lwork,n
    integer,                     allocatable :: iwork(:)
    real(kind=8),                allocatable :: rwork(:), w(:)
    complex(kind=8),             allocatable :: A(:,:),   work(:)
    integer                                  :: ierr

    ndim=size(matin(1,:))

    if (allocated(zvecs)) deallocate(zvecs)
    if (allocated(zvals)) deallocate(zvals)
    allocate(zvecs(ndim,ndim),zvals(ndim),stat=ierr)
    if (ierr /= 0) STOP "*** not enough memory ***"

    n=ndim
    lda=n

    lwork  = 2*n+n**2
    lrwork = 1+5*n+2*n**2
    liwork = 3+5*n

    allocate(a(ndim,ndim),w(ndim),work(lwork),&
             rwork(lrwork),iwork(liwork),stat=ierr)
    if (ierr /= 0) STOP "*** not enough memory ***"

    a=matin

    call zheevd(jobz,uplo,n,a,lda,w,work,lwork,rwork,lrwork,iwork,liwork,info)

    zvals=w
    zvecs=a

    deallocate(a,w,rwork,iwork,work)

end subroutine

end

Answer 2

此链接表明缺少的方法是英特尔库调用：

http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htm#GUID-9AD3B5B7-DC35-4DF7-A126-9A8730FE98CA.htm

我想知道它的等价物在 Mac 上的什么地方？你有 gfortran 吗？

http://gcc.gnu.org/wiki/GFortranBinaries

底部消息是否表明没有 64 位版本可用？

抱歉 - 问题多于答案。