Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

[FORTRAN] Help figuring out zheev

Tags:
  1. Sep 7, 2014 #1

    Matterwave

    User Avatar
    Science Advisor
    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 (Text):

    *  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 (Text):
        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: Sep 7, 2014
  2. jcsd
  3. Sep 8, 2014 #2

    Mark44

    Staff: Mentor

    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.
    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.
    LDA stands for "leading dimension of A". IOW, this is the number of rows in your matrix.
    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".
    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.
     
  4. Sep 8, 2014 #3

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    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 (Text):
    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
     
  5. Sep 8, 2014 #4

    DrClaude

    User Avatar

    Staff: Mentor

    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.

    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.

    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,*).

    It's written in Fortran 77, not this newfangled "Fortran 90" or whatever it is kids are using nowadays. In my days... :grumpy:
     
  6. Sep 8, 2014 #5

    DrClaude

    User Avatar

    Staff: Mentor

    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.

    Annd please please please remove that eyesore
    Code (Text):
    implicit real(kind=reel8) (a-h,o-z)
    and replace it by
    Code (Text):
    implicit none
    No serious modern program should have any variable that is not explicitly declared.
     
  7. Sep 8, 2014 #6

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    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?


    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: Sep 8, 2014
  8. Sep 9, 2014 #7

    DrClaude

    User Avatar

    Staff: Mentor

    Sorry, I didn't notice that there were two variables, nflavor and nflavor2.

    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.

    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.

    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.

    <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 (Text):
    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
     
  9. Sep 9, 2014 #8

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    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 any one 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 (Text):
        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: Sep 9, 2014
  10. Sep 9, 2014 #9

    Mark44

    Staff: Mentor

    What's with the "reel8"? Is this a user-defined variable that is an alias for REAL*8?
     
  11. Sep 9, 2014 #10

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    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 (Text):


    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...
     
  12. Sep 9, 2014 #11

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    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 (Text):

    real integer,character,complex
    integer real
    ...
    real = sqrt(character) + integer**complex
     
    you can, but it's not recommended!
     
  13. Sep 9, 2014 #12

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    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?
     
  14. Sep 10, 2014 #13

    DrClaude

    User Avatar

    Staff: Mentor

    That one I hesitated with. I do admit that adding (:) and (:,:) does make it clearer that they are arrays assignments.

    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.

    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 (Text):

      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.

    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 (Text):

     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 (Text):

     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.

    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.
     
  15. Sep 10, 2014 #14

    DrClaude

    User Avatar

    Staff: Mentor

    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 (Text):

      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.
     
  16. Sep 10, 2014 #15

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    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?
     
  17. Sep 10, 2014 #16

    DrClaude

    User Avatar

    Staff: Mentor

    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 (Text):

     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
     
     
  18. Sep 10, 2014 #17

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    So, the reason I would make these variables private and the like, is purely for speed, and memory purposes of running the code?
     
  19. Sep 10, 2014 #18

    DrClaude

    User Avatar

    Staff: Mentor

    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 (Text):

      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
     
     
  20. Sep 11, 2014 #19

    Matterwave

    User Avatar
    Science Advisor
    Gold Member

    Thank you for the help DrClaude!
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: [FORTRAN] Help figuring out zheev
  1. Fortran Help (Replies: 2)

  2. Fortran If help (Replies: 3)

Loading...