# [FORTRAN] Help figuring out zheev

• Fortran
• Matterwave
In summary: LWORK.*> \endverbatim* In summary, the zheev subroutine computes all eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix A.
Matterwave
Gold Member
Hello guys, I am in need of a subroutine that diagonalizes a nxn Hermitian matrix (really I'm just looking for the eigenvalues and eigenvectors of course). Looking online I found that LAPACK has a zheev subroutine that presumably does just this. I am a little confused on how to use this subroutine though, I'd appreciate it if you guys could take a look and help me with some of my questions:

Code:
*  Definition:
*  ===========
*
*       SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
*                         INFO )
*
*       .. Scalar Arguments ..
*       CHARACTER          JOBZ, UPLO
*       INTEGER            INFO, LDA, LWORK, N
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   RWORK( * ), W( * )
*       COMPLEX*16         A( LDA, * ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZHEEV computes all eigenvalues and, optionally, eigenvectors of a
*> complex Hermitian matrix A.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] JOBZ
*> \verbatim
*>          JOBZ is CHARACTER*1
*>          = 'N':  Compute eigenvalues only;
*>          = 'V':  Compute eigenvalues and eigenvectors.
*> \endverbatim
*>
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*>          = 'U':  Upper triangle of A is stored;
*>          = 'L':  Lower triangle of A is stored.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix A.  N >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is COMPLEX*16 array, dimension (LDA, N)
*>          On entry, the Hermitian matrix A.  If UPLO = 'U', the
*>          leading N-by-N upper triangular part of A contains the
*>          upper triangular part of the matrix A.  If UPLO = 'L',
*>          the leading N-by-N lower triangular part of A contains
*>          the lower triangular part of the matrix A.
*>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
*>          orthonormal eigenvectors of the matrix A.
*>          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
*>          or the upper triangle (if UPLO='U') of A, including the
*>          diagonal, is destroyed.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max(1,N).
*> \endverbatim
*>
*> \param[out] W
*> \verbatim
*>          W is DOUBLE PRECISION array, dimension (N)
*>          If INFO = 0, the eigenvalues in ascending order.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The length of the array WORK.  LWORK >= max(1,2*N-1).
*>          For optimal efficiency, LWORK >= (NB+1)*N,
*>          where NB is the blocksize for ZHETRD returned by ILAENV.
*>
*>          If LWORK = -1, then a workspace query is assumed; the routine
*>          only calculates the optimal size of the WORK array, returns
*>          this value as the first entry of the WORK array, and no error
*>          message related to LWORK is issued by XERBLA.
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*>          RWORK is DOUBLE PRECISION array, dimension (max(1, 3*N-2))
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>          = 0:  successful exit
*>          < 0:  if INFO = -i, the i-th argument had an illegal value
*>          > 0:  if INFO = i, the algorithm failed to converge; i
*>                off-diagonal elements of an intermediate tridiagonal
*>                form did not converge to zero.
*> \endverbatim

So, the above are the comments on what all the inputs and outputs of this subroutine are. There's a few things I'm not sure of.

1) It looks to me like this subroutine will take a nxn matrix A and then replace that matrix with a matrix of its eigenvectors, is that correct?

In this case, I should first create a dummy array for A so that my matrix is not destroyed in the process of calling this subroutine correct?

2) What is this LDA? I can't make heads or tails of it other than it should be "N" since Hermitian matrices are square...I don't understand why someone would write this redundant parameter if that's the case though.

3) I have no idea what's the use of WORK, LWORK, and RWORK. Can anyone make heads or tails of this?

Next, take a look at how this subroutine is defined:

Code:
    SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
$INFO ) * * -- LAPACK driver routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2011 * * .. Scalar Arguments .. CHARACTER JOBZ, UPLO INTEGER INFO, LDA, LWORK, N * .. * .. Array Arguments .. DOUBLE PRECISION RWORK( * ), W( * ) COMPLEX*16 A( LDA, * ), WORK( * ) I am completely confused on how the arrays are defined. I am used to definitions of arrays like: complex,dimension(N,N)::A How come his subroutine's array definitions don't seem to require a "dimension" declaration? How does the compiler know how much memory to allocate to these arrays? Not to mention the syntax lacks a "::" now. Zheev should be written in fortran, it's extension is .f. I am not familiar with this syntax at all... Last edited: Matterwave said: Hello guys, I am in need of a subroutine that diagonalizes a nxn Hermitian matrix (really I'm just looking for the eigenvalues and eigenvectors of course). Looking online I found that LAPACK has a zheev subroutine that presumably does just this. I am a little confused on how to use this subroutine though, I'd appreciate it if you guys could take a look and help me with some of my questions: Code: * Definition: * =========== * * SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, * INFO ) * * .. Scalar Arguments .. * CHARACTER JOBZ, UPLO * INTEGER INFO, LDA, LWORK, N * .. * .. Array Arguments .. * DOUBLE PRECISION RWORK( * ), W( * ) * COMPLEX*16 A( LDA, * ), WORK( * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> ZHEEV computes all eigenvalues and, optionally, eigenvectors of a *> complex Hermitian matrix A. *> \endverbatim * * Arguments: * ========== * *> \param[in] JOBZ *> \verbatim *> JOBZ is CHARACTER*1 *> = 'N': Compute eigenvalues only; *> = 'V': Compute eigenvalues and eigenvectors. *> \endverbatim *> *> \param[in] UPLO *> \verbatim *> UPLO is CHARACTER*1 *> = 'U': Upper triangle of A is stored; *> = 'L': Lower triangle of A is stored. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER *> The order of the matrix A. N >= 0. *> \endverbatim *> *> \param[in,out] A *> \verbatim *> A is COMPLEX*16 array, dimension (LDA, N) *> On entry, the Hermitian matrix A. If UPLO = 'U', the *> leading N-by-N upper triangular part of A contains the *> upper triangular part of the matrix A. If UPLO = 'L', *> the leading N-by-N lower triangular part of A contains *> the lower triangular part of the matrix A. *> On exit, if JOBZ = 'V', then if INFO = 0, A contains the *> orthonormal eigenvectors of the matrix A. *> If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') *> or the upper triangle (if UPLO='U') of A, including the *> diagonal, is destroyed. *> \endverbatim *> *> \param[in] LDA *> \verbatim *> LDA is INTEGER *> The leading dimension of the array A. LDA >= max(1,N). *> \endverbatim *> *> \param[out] W *> \verbatim *> W is DOUBLE PRECISION array, dimension (N) *> If INFO = 0, the eigenvalues in ascending order. *> \endverbatim *> *> \param[out] WORK *> \verbatim *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. *> \endverbatim *> *> \param[in] LWORK *> \verbatim *> LWORK is INTEGER *> The length of the array WORK. LWORK >= max(1,2*N-1). *> For optimal efficiency, LWORK >= (NB+1)*N, *> where NB is the blocksize for ZHETRD returned by ILAENV. *> *> If LWORK = -1, then a workspace query is assumed; the routine *> only calculates the optimal size of the WORK array, returns *> this value as the first entry of the WORK array, and no error *> message related to LWORK is issued by XERBLA. *> \endverbatim *> *> \param[out] RWORK *> \verbatim *> RWORK is DOUBLE PRECISION array, dimension (max(1, 3*N-2)) *> \endverbatim *> *> \param[out] INFO *> \verbatim *> INFO is INTEGER *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value *> > 0: if INFO = i, the algorithm failed to converge; i *> off-diagonal elements of an intermediate tridiagonal *> form did not converge to zero. *> \endverbatim So, the above are the comments on what all the inputs and outputs of this subroutine are. There's a few things I'm not sure of. 1) It looks to me like this subroutine will take a nxn matrix A and then replace that matrix with a matrix of its eigenvectors, is that correct? You can use it for an N x N matrix, but the subroutine is more general than that, in that it can work with a rectangular matrix with a different number of rows than columns. Matterwave said: In this case, I should first create a dummy array for A so that my matrix is not destroyed in the process of calling this subroutine correct? The A parameter in the ZHEEV subroutine is both in and out, which means that this parameter is modified. It would be a good idea to store the matrix in another variable. Matterwave said: 2) What is this LDA? I can't make heads or tails of it other than it should be "N" since Hermitian matrices are square...I don't understand why someone would write this redundant parameter if that's the case though. LDA stands for "leading dimension of A". IOW, this is the number of rows in your matrix. Matterwave said: 3) I have no idea what's the use of WORK, LWORK, and RWORK. Can anyone make heads or tails of this? The documentation isn't very clear on what these are for. I suggest you run the routine, supplying all parameters as indicated by the documentation, and see what you get. Be sure to supply initialized variables of the appropriate size of any parameters marked "in" and supply uninitialized parameters for those marked "out". Matterwave said: Next, take a look at how this subroutine is defined: Code:  SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,$                  INFO )
*
*  -- LAPACK driver routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
CHARACTER          JOBZ, UPLO
INTEGER            INFO, LDA, LWORK, N
*     ..
*     .. Array Arguments ..
DOUBLE PRECISION   RWORK( * ), W( * )
COMPLEX*16         A( LDA, * ), WORK( * )

I am completely confused on how the arrays are defined. I am used to definitions of arrays like:

complex,dimension(N,N)::A

How come his subroutine's array definitions don't seem to require a "dimension" declaration? How does the compiler know how much memory to allocate to these arrays?

Not to mention the syntax lacks a "::" now. Zheev should be written in fortran, it's extension is .f. I am not familiar with this syntax at all...

The LAPACK routines are written in an older form of Fortran. Possibly it's Fortran 77, in which arrays are declared without a DIMENSION keyword, and also without the :: thing.

Hmmm, ok, I can try it out I guess. It's somewhat hard since I'm working on only a snippet of the full code, and to run the full code will take a lot more modifications. I guess I will run a dummy program to test out this code and see what it does.

I don't understand how a Hermitian matrix can be a rectangular matrix though? Does that even make sense?

Also, perhaps you can take a look at my code which calls this function and see if you notice any big errors?

Code:
subroutine exphcalsp(hnow,drdis,exphel)
use params
use helectron
implicit real(kind=reel8) (a-h,o-z)
complex(kind=reel8),dimension(nflavor2,nflavor2),intent(IN) :: hnow
complex(kind=reel8),dimension(nflavor2,nflavor2),intent(OUT)::exphel
complex(kind=reel8),dimension(nflavor2,nflavor2) :: eigvec
real(kind=reel8),dimension(nflavor2) :: eigval
complex(kind=reel8),dimension(nflavor2) :: ceigval
real(kind=reel8) :: drdis
real(kind=reel8) :: xa,xb,xc,xl,xldr,xli,cldr,sldr,xbs,xas
complex*16,dimension(:),allocatable::work
double precision,dimension(:),allocatable::rwork
integer :: info

if(nflavor.eq.2) then

allocate(work(7))
allocate(rwork(10))
eigvec(:,:)=hnow(:,:)

call zheev(v,u,4,eigvec,4,eigval,work,7,rwork,info)
ceigval(:)=exp(-ci*eigval(:)*drdis)
exphel(:,:)=(0.d0,0.d0)
do i=1,nflavor2
do j=1,nflavor2
do k=1,nflavor2
exphel(i,j)= exphel(i,j)+ eigvec(i,k)*ceigval(k)*conjg(eigvec(j,k))
enddo
enddo
enddo

elseif(nflavor.eq.3)then

allocate(work(10))
allocate(rwork(16))
eigvec(:,:)=hnow(:,:)

call zheev(v,u,6,eigvec,6,eigval,work,10,rwork,info)
ceigval(:)=exp(-ci*eigval(:)*drdis)
exphel(:,:)=(0.d0,0.d0)
do i=1,nflavor2
do j=1,nflavor2
do k=1,nflavor2
exphel(i,j)= exphel(i,j)+ eigvec(i,k)*ceigval(k)*conjg(eigvec(j,k))
enddo
enddo
enddo

else
stop 'nflavor > 3 not implemented yet'
endif
return
end subroutine exphcalsp

end module spincohere

Matterwave said:
2) What is this LDA? I can't make heads or tails of it other than it should be "N" since Hermitian matrices are square...I don't understand why someone would write this redundant parameter if that's the case though.
As you mention yourself, indeed a Hermitian matrix should be square. But the memory storage of that matrix doesn't need to be the same size. This is a remnant from old Fortran programming constraints, where the size of an array would need to be fixed at compile time. You could then declare the matrix A to be of dimensions (Nmax,Nmax), but only at runtime decide the actual size the matrix N ≤ Nmax. The spacing between the first element of each column is then LDA = Nmax.

Matterwave said:
3) I have no idea what's the use of WORK, LWORK, and RWORK. Can anyone make heads or tails of this?
This comes again from the constraints I mentionned. Before Fortran 90, a subroutine couldn't allocate any additional memory it needed. So the subroutine is asking you to provide arrays WORK and RWORK, of big enough size, which it will use as temporary (scratch) memory.

Matterwave said:
Next, take a look at how this subroutine is defined:

Code:
    SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
\$                  INFO )
*
*  -- LAPACK driver routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
CHARACTER          JOBZ, UPLO
INTEGER            INFO, LDA, LWORK, N
*     ..
*     .. Array Arguments ..
DOUBLE PRECISION   RWORK( * ), W( * )
COMPLEX*16         A( LDA, * ), WORK( * )

I am completely confused on how the arrays are defined. I am used to definitions of arrays like:

complex,dimension(N,N)::A

How come his subroutine's array definitions don't seem to require a "dimension" declaration? How does the compiler know how much memory to allocate to these arrays?
The point is that they are not allocated in the subroutine. These arrays will be whatever size they are declared as in the portion of the program calling this subroutine. The subroutine doesn't known about the size, except through the other arguments. Notice that it still needs to know exactly the number of elements between each column, so the leading dimension of A has to be explicit: A(LDA,*).

Matterwave said:
Not to mention the syntax lacks a "::" now. Zheev should be written in fortran, it's extension is .f. I am not familiar with this syntax at all...
It's written in Fortran 77, not this newfangled "Fortran 90" or whatever it is kids are using nowadays. In my days...

Matterwave said:
Code:
      if(nflavor.eq.2) then

allocate(work(7))
allocate(rwork(10))
eigvec(:,:)=hnow(:,:)

call zheev(v,u,4,eigvec,4,eigval,work,7,rwork,info)
Are you using explicit values for the size of work and rwork, based a conditional test of nflavor? If so, please change that. Calculate LWORK and the size of RWORK as variables, and use them. Also N is the order of the matrix. If eigvec if of size(nflavor,nflavor), then N and LDA = nflavor.

Code:
implicit real(kind=reel8) (a-h,o-z)
and replace it by
Code:
implicit none
No serious modern program should have any variable that is not explicitly declared.

DrClaude said:
Are you using explicit values for the size of work and rwork, based a conditional test of nflavor? If so, please change that. Calculate LWORK and the size of RWORK as variables, and use them. Also N is the order of the matrix. If eigvec if of size(nflavor,nflavor), then N and LDA = nflavor.

They are of size nflavor2=2*nflavor, so I put that explicitly in as integers. I can replace the 4's and 6's with the nflavor2 parameter, but it shouldn't make a difference right?

I put the values of Lwork I found as the documentation stated LWORK>=blah. I just took the smallest value LWORK=blah, with N my dimension size, and put it in. I don't understand how it's calculated?

Are you saying I should declare an integer LWORK = 2*nflavor2-1 and put that in there instead of 7 and 10 respectively?

If that's the case, I thought that array allocations must be given parameters instead of variables so that the compiler knows how much memory to allocate it?
Code:
implicit real(kind=reel8) (a-h,o-z)
and replace it by
Code:
implicit none
No serious modern program should have any variable that is not explicitly declared.

The code is a very large code spanning 6-10 modules, all of which basically use this convention. Since I'm not the one who wrote this code (I'm just modifying a portion of it to do what I want it to do) I'm too afraid to change this kind of stuff lol.

Last edited:
Matterwave said:
They are of size nflavor2=2*nflavor, so I put that explicitly in as integers.
Sorry, I didn't notice that there were two variables, nflavor and nflavor2.

Matterwave said:
I can replace the 4's and 6's with the nflavor2 parameter, but it shouldn't make a difference right?
Yes, but it's bad practice. First, it forces you to have separate code for each case. Second, it makes the program much less readable.

Matterwave said:
I put the values of Lwork I found as the documentation stated LWORK>=blah. I just took the smallest value LWORK=blah, with N my dimension size, and put it in. I don't understand how it's calculated?
You have have to understand the inner workings of zheev to understand how it's calculated. It is indeed standard practice to calculate LWORK to its minimum value, which is often exactly what is needed.

Matterwave said:
Are you saying I should declare an integer LWORK = 2*nflavor2-1 and put that in there instead of 7 and 10 respectively?

Matterwave said:
If that's the case, I thought that array allocations must be given parameters instead of variables so that the compiler knows how much memory to allocate it?
When you are using allocate, the allocation is in any case performed at runtime. It is only when an array is declared with a fixed size that the size must be known at compile time.

Matterwave said:
The code is a very large code spanning 6-10 modules, all of which basically use this convention. Since I'm not the one who wrote this code (I'm just modifying a portion of it to do what I want it to do) I'm too afraid to change this kind of stuff lol.
<feigned indignation> Then I won't be helping you anymore! I refuse to debug any code with implicit declarations. </feigned indignation>

In Pf, we usually teach how to fish, not give out fish, but since you're a SA yourself, I guess it's ok if I give you an example of how I would use zheev.

Code:
subroutine exphcalsp(hnow,drdis,exphel)
use params
use helectron
implicit real(kind=reel8) (a-h,o-z)
complex(kind=reel8),dimension(nflavor2,nflavor2),intent(IN) :: hnow
complex(kind=reel8),dimension(nflavor2,nflavor2),intent(OUT)::exphel
complex(kind=reel8),dimension(nflavor2,nflavor2) :: eigvec
real(kind=reel8),dimension(nflavor2) :: eigval
complex(kind=reel8),dimension(nflavor2) :: ceigval
real(kind=reel8) :: drdis  ! Shouldn't this have intent(in)?
real(kind=reel8) :: xa,xb,xc,xl,xldr,xli,cldr,sldr,xbs,xas ! Not used !
complex(kind=reel8), dimension(:), allocatable :: work  ! Declarations should be consistent
real(kind=reel8), allocatable :: rwork ! Declarations should be consistent
integer :: lwork, info

! Allocate temporary work space
lwork = 2 * nflavor2 + 1
allocate(work(lwork))
allocate(rwork(3 * nflavor - 2))

! Assign matrix
eigvec = hnow  ! Note that the (:,:) are not necessary

! Calculate eigenvalues and eigenvectors
call zheev('V', 'U' , nflavor2, eigvec, nflavor2, eigval, work, lwork, rwork, info)
! Note that 'U' and 'V' are character constants, not variables

! Free temporary memory
deallocate(work, rwork)

! I don't know what this part does, so I don't know what comment to put here
ceigval = exp(-ci * eigval * drdis) ! Note that (:) is not necessary
exphel = (0.d0,0.d0) ! Note that (:,:) is not necessary
do k=1,nflavor2
do j=1,nflavor2
do i=1,nflavor2
exphel(i,j) = exphel(i,j) + eigvec(i,k) * ceigval(k) * conjg(eigvec(j,k))
enddo
enddo
enddo
! Note that I inversed the order of the loops.  Multi-dim. arrays in Fortran are stored in a
! row-wise order, such that the first index should be the one changing the fastest

end module spincohere

Interesting edits!

As I said, I did not write this code, I am merely modifying it, so a lot of those points should be directed at the original writer of this code lol!

For example, everywhere I see array assignments I see the (:,:) in this code, so I just stuck with it. Also, most of those declarations at the top (e.g. of the variables that aren't used) were not made by me lol. When reading this code, I am often very confused by unused declared variables!

A couple of follow up questions:
1) I was told in the documentation of zheev what "type" of arrays rwork, etc. have to be, which is why I declared them with the type that was asked of me (double precision, etc). The code I am running uses higher (I think) precision floats (kind=reel8 is selected_real_kind(15,9)), is it legit for me to declare the arrays zheev will be working on using the precision in my code, rather than the precision used in zheev?

2) In my code, the nflavor parameter is a parameter and will never change in anyone running of the code. Which is why I have allocate statements within if-statements and never any deallocate statements. However, this subroutine exphcal will be called many (thousands/millions) times (but never changing nflavor), is it essential that I deallocate the arrays after zheev is done with them?

Here's my modified code, what do you think?

Code:
    subroutine exphcalsp(hnow,drdis,exphel)
use params
use helectron
implicit real(kind=reel8) (a-h,o-z)
complex(kind=reel8),dimension(nflavor2,nflavor2),intent(IN) :: hnow
complex(kind=reel8),dimension(nflavor2,nflavor2),intent(OUT):: exphel
complex(kind=reel8),dimension(nflavor2,nflavor2) :: eigvec
real(kind=reel8),dimension(nflavor2) :: eigval
complex(kind=reel8),dimension(nflavor2) :: ceigval
real(kind=reel8) :: drdis
real(kind=reel8) :: xa,xb,xc,xl,xldr,xli,cldr,sldr,xbs,xas
complex(kind=reel8),dimension(:),allocatable :: work
real(kind=reel8),dimension(:),allocatable :: rwork
integer :: lwork
integer :: info

lwork=2*nflavor2-1

if(nflavor.eq.1) then
xa=real(hnow(1,1),kind=reel8)
xb=real(hnow(1,2),kind=reel8)
xc=dimag(hnow(1,2))

xl = sqrt(xa**2+xb**2+xc**2)
xldr=xl*drdis
xli=1.d0/xl

cldr=cos(xldr)
sldr=sin(xldr)
xbs=xli*sldr
xas=xa*xbs

exphel(1,1)=cmplx(cldr,-xas,kind=reel8)
exphel(2,2)=cmplx(cldr, xas,kind=reel8)
exphel(1,2)=xbs*cmplx(xc,-xb,kind=reel8)
exphel(2,1)=xbs*cmplx(-xc,-xb,kind=reel8)

if(nflavor.eq.2) then

allocate(work(lwork))
allocate(rwork(3*nflavor2-2))
eigvec(:,:)=hnow(:,:)

call zheev(v,u,nflavor2,eigvec,nflavor2,eigval,work,lwork,rwork,info)
ceigval(:)=exp(-ci*eigval(:)*drdis)
exphel(:,:)=(0.d0,0.d0)
do i=1,nflavor2
do j=1,nflavor2
do k=1,nflavor2
exphel(i,j)= exphel(i,j)+ eigvec(i,k)*ceigval(k)*conjg(eigvec(j,k))
enddo
enddo
enddo

elseif(nflavor.eq.3)then

allocate(work(lwork))
allocate(rwork(3*nflavor2-2))
eigvec(:,:)=hnow(:,:)

call zheev(v,u,nflavor2,eigvec,nflavor2,eigval,work,lwork,rwork,info)
ceigval(:)=exp(-ci*eigval(:)*drdis)
exphel(:,:)=(0.d0,0.d0)
do i=1,nflavor2
do j=1,nflavor2
do k=1,nflavor2
exphel(i,j)= exphel(i,j)+ eigvec(i,k)*ceigval(k)*conjg(eigvec(j,k))
enddo
enddo
enddo

else
stop 'nflavor > 3 not implemented yet'
endif
return
end subroutine exphcalsp

EDIT: In fact, I think I can just have 1 set of allocate statements outside the if statements here if nflavor never changes?

Also, I'm not sure what you mean with the comment on changing the order of the loops. Does that really matter since eventually all the columns and rows of the matrices will be filled? o.o

Last edited:
What's with the "reel8"? Is this a user-defined variable that is an alias for REAL*8?

It is user defined integer parameter, in the module params. I...don't really know much about it other than it is selecting a set precision:

Code:
integer,parameter :: reel4 = selected_real_kind (6,5)
integer,parameter :: reel8 = selected_real_kind (15,9)
!   for using quadruple precision for distances
!     integer,parameter :: reel16 = selected_real_kind (30,9)
!   test using double precision only
integer,parameter :: reel16 = selected_real_kind (15,9)

I am not familiar with the "selected_real_kind(:,:)" statement...

Matterwave said:
I am not familiar with the "selected_real_kind(:,:)" statement...
Actually it's a built-in function, not a statement.

The idea is that it returns the "best" type of variable to handle a given range of values. selected_real_kind (15,9) means you want at least 15 significant figures of accuracy, and exponents up to at least 109.

Most compilers have now unofficially standardized on kind=4 for 32 bit (4 byte) IEEE floating point and kind=8 for 64 bit (8 byte), but historically there were machines with different hardware designs. For example the old CDC scientific computers used 60-bit float point arithmetic (and nothing smaller) and on Cray supercomputers the only size was 64-bit.

I suppose somebody thought reel8 would not clash with the "reserved word" REAL, but actually in fortran there are no reserved words. If you want to write code like
Code:
real integer,character,complex
integer real
...
real = sqrt(character) + integer**complex
you can, but it's not recommended!

AlephZero said:
Actually it's a built-in function, not a statement.

The idea is that it returns the "best" type of variable to handle a given range of values. selected_real_kind (15,9) means you want at least 15 significant figures of accuracy, and exponents up to at least 109.

Most compilers have now unofficially standardized on kind=4 for 32 bit (4 byte) IEEE floating point and kind=8 for 64 bit (8 byte), but historically there were machines with different hardware designs. For example the old CDC scientific computers used 60-bit float point arithmetic (and nothing smaller) and on Cray supercomputers the only size was 64-bit.

I suppose somebody thought reel8 would not clash with the "reserved word" REAL, but actually in fortran there are no reserved words. If you want to write code like
Code:
real integer,character,complex
integer real
...
real = sqrt(character) + integer**complex
you can, but it's not recommended!

Interesting! How does this type of selection of the accuracy get stored into an integer parameter tho? I thought integers should just be one single integer?

Matterwave said:
For example, everywhere I see array assignments I see the (:,:) in this code, so I just stuck with it.
That one I hesitated with. I do admit that adding (:) and (:,:) does make it clearer that they are arrays assignments.

Matterwave said:
1) I was told in the documentation of zheev what "type" of arrays rwork, etc. have to be, which is why I declared them with the type that was asked of me (double precision, etc). The code I am running uses higher (I think) precision floats (kind=reel8 is selected_real_kind(15,9)), is it legit for me to declare the arrays zheev will be working on using the precision in my code, rather than the precision used in zheev?
The use of "double precision" or "complex*16" is deprecated in Fortran 90/95, and the new declaration format should be used. As for zheev, in any case everything must be consistant: if declaring A complex(kind=reel8) is correct for zheev, which expects complex*16, then using complex(kind=reel8) for WORK and real(kind=reel8) for RWORK is correct.

Matterwave said:
2) In my code, the nflavor parameter is a parameter and will never change in anyone running of the code. Which is why I have allocate statements within if-statements and never any deallocate statements. However, this subroutine exphcal will be called many (thousands/millions) times (but never changing nflavor), is it essential that I deallocate the arrays after zheev is done with them?
If the arrays are allocated in the subroutine, they will be deallocated automatically on exit anyway. If the size of the matrix A will not change, I would put the work arrays in a module so that there is no allocation/deallocation at each call, which is a waste of time. Here is what it would look like:

Code:
  module workspace
integer :: lwork, lrwork, size, info
real(kind=reel8), dimension (:), allocatable :: rwork
complex(kind=reel8), dimension (:), allocatable :: work
logical :: isInitialized = .false.
end module workspace

subroutine exphcalsp(hnow,drdis,exphel)
use params
use helectron
use workspace
implicit real(kind=reel8) (a-h,o-z)
complex(kind=reel8),dimension(nflavor2,nflavor2),intent(IN) :: hnow
complex(kind=reel8),dimension(nflavor2,nflavor2),intent(OUT)::exphel
complex(kind=reel8),dimension(nflavor2,nflavor2) :: eigvec
real(kind=reel8),dimension(nflavor2) :: eigval
complex(kind=reel8),dimension(nflavor2) :: ceigval
real(kind=reel8) :: drdis
real(kind=reel8) :: xa,xb,xc,xl,xldr,xli,cldr,sldr,xbs,xas

! Verify that the matrix size has not changed
if (isInitialized .and. nflavor2 /= size) then
deallocate(work, rwork)
isInitialized = .false.
end if

! Initialize work arrays if needed
if (.not.isInitialized) then
lwork = 2 * nflavor2 - 1
lrwork = 3 * nflavor2 - 2
size = nflavor2
allocate(work(lwork), rwork(lrwork))
isInitialized = .true.
end if

! Rest of code...

end subroutine exphcalsp
Size work and rwork are part of a module, they will not be deallocated upon exit, and the same memory space will be resued with little overhead. Note that I also included a failsafe in case the code ever gets modified and different values of nflavor2 are used in subsequent calls.

Matterwave said:
Here's my modified code, what do you think?
There are some "end if" or "else" missing. Every piece of code that is repeated can lead to debugging headaches, if one part that is repeated gets modified for one case but not others. The following arrangement:
Code:
 if (case == 1) then
! case 1 specific code
else if (case == 2) then
! case 2 specific code
else
! generic code for other cases
end if

! Common code A

if (case == 1) then
! case 1 specific code
else if (case == 2) then
! case 2 specific code
else
! generic code for other cases
end if

! Common code B
is much better than
Code:
 if (case == 1) then
! case 1 specific code
! Common code A
! case 1 specific code
! Common code B
else if (case == 2) then
! case 2 specific code
! Common code A
! case 2 specific code
! Common code B
else
! generic code for other cases
! Common code A
! generic code for other cases
! Common code B
end if
because any change in code A or code B would require things to be changed in many places, in exactly the same way.

Matterwave said:
Also, I'm not sure what you mean with the comment on changing the order of the loops. Does that really matter since eventually all the columns and rows of the matrices will be filled?
Since your arrays are small, it won't change anything, but it is good practice to do it right! If the arrays are so big that they don't fit in cache, then the CPU will have to keep fetching the data from slower memory at each iteration of the loop. The cache is usually filled by blocks of memory, so addressing elements in the same order they appear in memory will reduce the number of out-of-cache calls.

Matterwave said:
Interesting! How does this type of selection of the accuracy get stored into an integer parameter tho? I thought integers should just be one single integer?
In most systems, the integer is simply the number of bytes, so different requests for selected_real_kind will return the same value. For instance:
Code:
  write(*,*) selected_real_kind(15,9)
write(*,*) selected_real_kind(14,9)
write(*,*) selected_real_kind(7,8)
write(*,*) selected_real_kind(6,9)
outputs 8, 8, 8, and 4 on my computer.

Ah ok. Thank you DrClaude! You have been most helpful. I'll take your advice into consideration and implement them into my code. :)

One more question though. You suggested I write a module workspace. The piece of code I am working on is a module in itself (module spincohere which includes all these subroutines), I assume that I can't have a module within a module, so I would have to write that piece of code outside?

Matterwave said:
One more question though. You suggested I write a module workspace. The piece of code I am working on is a module in itself (module spincohere which includes all these subroutines), I assume that I can't have a module within a module, so I would have to write that piece of code outside?
In that case, you indeed can't create a module within a module. What you can do is make these variables a part of the module, but declare them as "private".
Code:
 module foobar
! Original variables that are part of the module appear here

! all the variables that follow will only be seen within the module foobar
integer, private :: lwork, lrwork, size, info
real(kind=reel8), dimension (:), allocatable, private :: rwork
complex(kind=reel8), dimension (:), allocatable, private :: work
logical, private :: isInitialized = .false.

contains

subroutine exphcalsp(hnow,drdis,exphel)
! etc.

end module foobar

So, the reason I would make these variables private and the like, is purely for speed, and memory purposes of running the code?

Matterwave said:
So, the reason I would make these variables private and the like, is purely for speed, and memory purposes of running the code?
The reason to have them private is to insure that there is no interference with anything piece of code outside the module that uses the module. The reason to have the variables as part of the entire module, instead of simply within the subroutine, was to save time by not continuously allocating and deallocating memory. But I now realize that the method I suggested is too complicated for no good reason. The variables should be within the subroutine, since they are only used there, and declared with a SAVE attribute. Here is what this looks like:
Code:
  subroutine exphcalsp(hnow,drdis,exphel)
use params
use helectron
implicit real(kind=reel8) (a-h,o-z)
complex(kind=reel8),dimension(nflavor2,nflavor2),intent(IN) :: hnow
complex(kind=reel8),dimension(nflavor2,nflavor2),intent(OUT)::exphel
complex(kind=reel8),dimension(nflavor2,nflavor2) :: eigvec
real(kind=reel8),dimension(nflavor2) :: eigval
complex(kind=reel8),dimension(nflavor2) :: ceigval
real(kind=reel8) :: drdis
real(kind=reel8) :: xa,xb,xc,xl,xldr,xli,cldr,sldr,xbs,xas

integer :: info
integer, save :: lwork, lrwork, size
real(kind=reel8), dimension (:), allocatable, save :: rwork
complex(kind=reel8), dimension (:), allocatable, save :: work
logical :: isInitialized = .false.  ! Note: assingment gives it automatically a save attribute

! Verify that the matrix size has not changed
if (isInitialized .and. nflavor2 /= size) then
deallocate(work, rwork)
isInitialized = .false.
end if

! Initialize work arrays if needed
if (.not.isInitialized) then
lwork = 2 * nflavor2 - 1
lrwork = 3 * nflavor2 - 2
size = nflavor2
allocate(work(lwork), rwork(lrwork))
isInitialized = .true.
end if

! Rest of code...

end subroutine exphcalsp

Thank you for the help DrClaude!

## 1. What is FORTRAN?

FORTRAN (Formula Translation) is a high-level programming language used for scientific and numerical computations. It was developed in the 1950s and is still widely used in scientific and engineering fields.

## 2. What does zheev do in FORTRAN?

ZHEEV is a FORTRAN subroutine used for computing all the eigenvalues and eigenvectors of a complex Hermitian matrix. It is part of the Linear Algebra Package (LAPACK) library and is commonly used in scientific computing applications.

## 3. How do I use zheev in my FORTRAN program?

To use zheev in your FORTRAN program, you need to include the LAPACK library and call the subroutine with the appropriate arguments. The exact steps will depend on your programming environment and the specific problem you are trying to solve.

## 4. Can zheev be used for non-Hermitian matrices?

No, zheev is specifically designed for complex Hermitian matrices. For non-Hermitian matrices, you can use the similar subroutine ZGEEV, which computes all the eigenvalues and eigenvectors of a general square matrix.

## 5. Are there any alternatives to zheev in FORTRAN?

Yes, there are several other subroutines in the LAPACK library that can be used for computing eigenvalues and eigenvectors in FORTRAN. Some alternatives include ZHEGV, which computes eigenvalues and eigenvectors of a generalized Hermitian-definite eigenproblem, and ZGESVD, which computes the singular value decomposition of a matrix.

• Programming and Computer Science
Replies
5
Views
1K
• Programming and Computer Science
Replies
2
Views
1K
• Programming and Computer Science
Replies
4
Views
1K
• Programming and Computer Science
Replies
4
Views
6K
• Programming and Computer Science
Replies
1
Views
8K
• Programming and Computer Science
Replies
7
Views
4K
• Programming and Computer Science
Replies
6
Views
2K
• Programming and Computer Science
Replies
21
Views
2K
• Programming and Computer Science
Replies
1
Views
2K
• Programming and Computer Science
Replies
4
Views
2K