Linking LAPACK in Fortran on Mac OS X

Question

I imagine this to be a standard noob problem but after spending all morning searching the web, I decided to bug you anyway. I'm on Mac OS 10.9 and I'd like to call a LAPACK eigenvalue routine from a Fortran program. I had the pleasure of being introduced to Fortran just yesterday, so please excuse any foolish mistakes.

This is the minimal example I want to get to run:

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

I learned that on OS X, LAPACK is part of the Accelerate framework, so I tried things like:

gfortran -o eigtest -framework accelerate eigtest.f95

but the linker complains:

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

I'm not sure that the Accelerate framework has that nice fortran95 call to heevd(A, eigs). Below, I wrap up the zheevd lapack subroutine (with all of its workspace variables) to keep the call nice and tight. You can store such ugliness inside an external module somewhere that you can load in your programs. I think intel MKL has most of the lapack95 interfaces, and there may be other better ways with the os x framework. The code below compiles with:

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