error in deallocating memory

[traptij]

Hi all,
I am using lapack routine DGEEV for finding the eigenvalue and eigenvector of a real non symmetric matrix which is written in fortran 77. The main program which calls DGEEV is written in fortran 90. I am using mingw g95 fortran compiler. In the main program DGEEV is first called to calculate the workspace required and then allocate the workspace. The second time DGEEV is called to calculate the eigenvalues and eigenvectors. After that a command is given to deallocate the memory for workspace. But at that line it gives the error that "could not deallocate memory" and sometimes "exception :access violation". Can anybody help me where could be the problem and how can i remove it?

Thanks,
T.Jain

[Julien Langou]

hello,
the routine dgeev in the lapack distribution provided at:
http://www.netlib.org/lapack/lapack.tgz
has a wrong workspace calculation. Try the one at:
http://www.netlib.org/lapack/patch/src/dgeev.f
There have been several fixes to the current version of LAPACK and they are all available at:
http://www.netlib.org/lapack-dev/lapack--3.0--patch--10042002.html
Let me know if this works better.
Julien

[traptij]

Hello,
No, it is not working, giving the same error.

Thanks,
T.Jain

[Julien Langou]

Hello,
you will need to post a piece of code if you want more help, all this should work then.
Julien

[traptij]

Hello,
Here is the code. The error is in eigenvector.f90 at the line "deallocate (work)". It is not able to deallocate the memory.
One more thing i could not understand is when the routine DTREVC is called in DGEEV sometimes it gives error that parameter 8 i.e.LDVL has an illegal value. But it is an input from eigenvector.f90 and not manipulated elsewhere in DGEEV then why the value of LDVL at that point is becoming zero(this i have checked).

Thanks,
T.Jain

Code: Select all

program eigenvector

!
! usage:
!    eigenvector gdxin i a gdxout eval evec
!
  use f90client

  implicit none

  character(len=255) :: iname,aname,evrname,eviname,evecrname,eveciname  ! name of parameter
  character(len=255) :: gdxin,gdxout   ! name of gdx files
  real(kind=8), allocatable :: a(:,:)  ! matrix
  integer :: n                         ! matrix size
  real(kind=8), allocatable :: wr(:) ! vector with eigenvalues
   real(kind=8), allocatable :: wi(:)   ! vector with imaginary eigenvalues
  real(kind=8), allocatable :: vr(:,:)   ! right eigenvectors
  real(kind=8), allocatable :: vl(:,:)   ! left eigenvectors
  real(kind=8), allocatable :: vrr(:,:)   ! real part of right eigenvectors
  real(kind=8), allocatable :: vri(:,:)   ! imaginary part of right eigenvectors
!  real :: vrr(n,n)
!  real :: vri(n,n)
  integer :: ap1,ap2,i,j,k                   ! gdxio pointers
  integer, allocatable :: index(:)     ! raw index elements
  integer :: ierr
  write(*,*) 'EIGENVECTOR, calculates eigenvalues/vectors of a non-symmetric matrix'

  call commandline

  call readmatrix

  call calc_nonsymmetric
   allocate(vrr(n,n),stat=ierr)
   if (ierr /= 0) call err('Could not allocate memory')
   allocate(vri(n,n),stat=ierr)
   if (ierr /= 0) call err('Could not allocate memory')

  k=0
   do i=1,n
    if   (wi(i) /= 0) then
      do j=1,n
        if (i /= n) then
        vrr(j,i)=vr(j,i+k)
        vri(j,i)=vr(j,i+k)
        k=k+1
        else
        vrr(j,i)=vr(j,n-1)
        vri(j,i)=-vr(j,n)
     endif
      enddo
    else
     do j=1,n
       vrr(j,i)=vr(j,i)
      vri(j,i)=0
      enddo
    endif
   enddo
    
  write(*,*) 'EVR:',wr
  write(*,*) 'EVI:',wi
  write(*,*) 'RREV:',vrr
  write(*,*) 'IREV:',vri

  call writeresults

  deallocate(vr,stat=ierr)
    if (ierr /= 0) call err('2.Could not deallocate memory')
  deallocate(vrr,stat=ierr)
   if (ierr /= 0) call err('3.Could not deallocate memory')
  deallocate(vri,stat=ierr)
    if (ierr /= 0) call err('4.Could not deallocate memory')
  deallocate(wr,stat=ierr)
  if (ierr /= 0) call err('5.Could not allocate memory')
  deallocate(wi,stat=ierr)
  if (ierr /= 0) call err('6.Could not allocate memory')


contains

!-----------------------------------------------------------------------------
! calculate eigenvalues for symmetric matrix
!-----------------------------------------------------------------------------

 subroutine calc_nonsymmetric
    real(kind=8),allocatable :: work(:)
    integer :: lwork
    integer :: info
    integer :: ierr
    real(kind=8) :: work1
   real(kind=8) :: nwork

    write(*,*) 'LAPACK DGEEV'
    write(*,*) 'Size : ',n

    allocate(wr(1:n),stat=ierr)
    if (ierr /= 0) call err('Could not allocate memory')
   allocate(wi(1:n),stat=ierr)
    if (ierr /= 0) call err('Could not allocate memory')
   call dgeev('N','V',n,a,n,wr,wi,vl,1,vr,n,work1,-1,info)
    if (info /= 0) then
       write(*,*) 'Workspace query call, INFO = ',info
       call err('Error in DGEEV')
    end if

    lwork = int(work1)
   write(*,*) 'Workspace (doubles) : ',lwork
    allocate(work(1:lwork),stat=ierr)
!   write(*,*)'ierr:',ierr
    if (ierr /= 0) call err('Could not allocate memory')
   allocate(vr(n,n),stat=ierr)
   if (ierr /= 0) call err('Could not allocate memory')
!    allocate(vl(n,n),stat=ierr)
!   if (ierr /= 0) call err('Could not allocate memory')
    call dgeev('N','V',n,a,n,wr,wi,vl,1,vr,n,nwork,lwork,info)
!   write(*,*)'info:',info
    if (info /= 0) then
       write(*,*) 'INFO = ',info
       call err('Error in DGEEV')
    end if
    
!   if (allocated(work)) deallocate(work)   
   deallocate(work,stat=ierr)
!   write(*,*)'ierr:',ierr
   if (ierr /= 0) call err('1.Could not deallocate memory')
   
  end subroutine calc_nonsymmetric

!-----------------------------------------------------------------------------
! write message and exit
!-----------------------------------------------------------------------------

  subroutine err(mess)

    character(len=*) mess

    write(*,*) trim(mess)
    stop

  end subroutine err

!-----------------------------------------------------------------------------
! report error in gdxio
!-----------------------------------------------------------------------------
  subroutine gdxerr(ap)

    integer :: ap
    integer :: ierr
    logical :: ok
    character(len=255) msg

    ierr = gdxgetlasterror(ap)
    ok = gdxerrorstr(ierr,msg)
    if (.not.ok) msg = 'Unknown error'
    call err(msg)

  end subroutine gdxerr

!-----------------------------------------------------------------------------
! binary search in integer array
!-----------------------------------------------------------------------------
  integer function bsearch(a,n,key)
    integer n,key
    integer a(1:n)
    integer i,j,k

    i = 1
    j = n
    do while (i <= j)
       k = (i+j)/2
       if (key.eq.a(k)) then
          bsearch = k
          return
       else if (key.lt.a(k)) then
          j = k - 1
       else
          i = k + 1
       end if
    end do
    bsearch = 0
  end function bsearch

!-----------------------------------------------------------------------------
! read matrix from GDX file
!-----------------------------------------------------------------------------
  subroutine readmatrix

    character(len=255) :: dllstr
    logical :: ok
    integer :: rc
    integer :: synr
    integer :: nrrecs
    character(len=255) :: name
    integer :: adim
    integer :: iatyp
    integer :: i,j,k
    integer :: ierr
    integer :: afdim

    integer :: aelements(1:10)
    real(kind=8) :: avals(5)

    ok = gdxdllversion(dllstr)
    if (.not.ok) call err('Could get dll version')
    write(*,*) 'GDXIO.DLL version: ',trim(dllstr)

    rc = gdxopenread(ap1,gdxin)
    if (rc /= 0) call err('Could not read file '//trim(gdxin))

!
! read set iname
!
    ok = gdxfindsymbol(ap1,iname,synr)
    if (.not.ok) call gdxerr(ap1)

    ok = gdxsymbolinfo(ap1,synr,name,adim,iatyp)
    if (.not.ok) call gdxerr(ap1)
    if (adim /= 1) call err('Dimension of set should 1')
    if (iatyp /= dt_set) call err('Set '//trim(iname)//' should be a set')

    ok = gdxdatareadrawstart(ap1,synr,n)
    if (.not.ok) call gdxerr(ap1)

    allocate(index(1:n),stat=ierr)
    if (ierr /= 0) call err('Could not allocate vector')

    do i=1,n
       ok = gdxdatareadraw(ap1,aelements,avals,afdim)
      if (.not.ok) call gdxerr(ap1)
       index(i) = aelements(1)
    end do
!
! read parameter aname
!
    ok = gdxfindsymbol(ap1,aname,synr)
    if (.not.ok) call gdxerr(ap1)

    ok = gdxsymbolinfo(ap1,synr,name,adim,iatyp)
    if (.not.ok) call gdxerr(ap1)
    if (adim /= 2) call err('Dimension of matrix should 2')
    if (iatyp /= dt_par) call err('Parameter '//trim(aname)//' should be a parameter')

    ok = gdxdatareadrawstart(ap1,synr,nrrecs)
    if (.not.ok) call gdxerr(ap1)

    allocate(a(1:n,1:n),stat=ierr)
    if (ierr /= 0) call err('Could not allocate matrix')

    a = 0

    do k=1,nrrecs
       ok = gdxdatareadraw(ap1,aelements,avals,afdim)
       if (.not.ok) call gdxerr(ap1)     
       i = bsearch(index,n,aelements(1))
       if (i == 0) call err('matrix element not in set')
       j = bsearch(index,n,aelements(2))
       if (j == 0) call err('matrix element not in set')
       a(i,j) = avals(1)
    end do

!
! check for symmetry
!

!    do i=1,n
!       do j=i+1,n
!          if (abs(a(i,j)-a(j,i))>1.0e-9) call err('Matrix is not symmetric')
!       end do
!    end do

  end subroutine readmatrix

!-----------------------------------------------------------------------------
! print usage and stop
!-----------------------------------------------------------------------------
  subroutine usage

    write(*,1000)
    stop


    1000 format( ' Usage '/,                                                     &
                 '    > eigenvector gdxin i a gdxout eval evecr eveci'/,                &
                 ' where'/,                                                      &
                 '    gdxin  : name of gdxfile with matrix'/,                    &
                 '    i      : name of set used in matrix'/,                     &
                 '    a      : name of 2 dimensional parameter inside gdxin'/,   &
                 '    gdxout : name of gdxfile for results (eigenvalues)'/,      &
                 '    eval   : name of 1 dimensional parameter inside gdxout'/,  &
                 '    evecr   : name of 2 dimensional parameter inside gdxout'/,  &
             '    eveci   : name of 2 dimensional parameter inside gdxout'/,  &
                 ' '/,                                                           &
                 ' Calculates eigenvalues/vectors of symmetric matrix a(i,j) where'/,   &
                 ' i and j are aliased sets. eval will contain the eigenvalues'/, &
                 ' and evec will contain the eigenvectors'/ &
                )

  end subroutine usage

!-----------------------------------------------------------------------------
! get command line parameters
!-----------------------------------------------------------------------------
  subroutine commandline

    integer(2) :: err

    call getarg(1,gdxin)
!    if (err == -1) call usage

    call getarg(2,iname)
!    if (err == -1) call usage

    call getarg(3,aname)
!    if (err == -1) call usage

    call getarg(4,gdxout)
!    if (err == -1) call usage

    call getarg(5,evrname)
!    if (err == -1) call usage

    call getarg(6,eviname)
!    if (err == -1) call usage

   call getarg(7,evecrname)
   call getarg(8,eveciname)

  end subroutine commandline

!-----------------------------------------------------------------------------
! write result vector to gdx file
!-----------------------------------------------------------------------------

  subroutine writeresults

    integer :: rc
    integer :: ierr
    integer :: i,j
    logical :: ok
    character(len=31) :: astrelements(1:10)
    real(kind=8) :: avals(1:5)
    character(len=31) :: s
    integer :: umap


    rc = gdxopenwrite(ap2,gdxout,'eigenvector')
    if (rc /= 0) call err('Could not write file '//trim(gdxout))

    ok = gdxdatawritestrstart(ap2,evrname,'real eigenvalues',1,dt_par,0)
    if (.not. ok) call gdxerr(ap2)

    astrelements(1:10) = ' '
    avals = 0.0d0

    do i=1,n
       ok = gdxumuelget(ap1,index(i),s,umap)
       if (.not. ok) call gdxerr(ap1)
       astrelements(1) = s
       avals(1) = wr(i)
       ok = gdxdatawritestr(ap2,astrelements,avals)
       if (.not. ok) call gdxerr(ap2)
    end do

    ok = gdxdatawritedone(ap2)
    if (.not. ok) call gdxerr(ap2)

    ok = gdxdatawritestrstart(ap2,eviname,'imaginary eigenvalues',1,dt_par,0)
    if (.not. ok) call gdxerr(ap2)

    astrelements(1:10) = ' '
    avals = 0.0d0

    do i=1,n
       ok = gdxumuelget(ap1,index(i),s,umap)
       if (.not. ok) call gdxerr(ap1)
       astrelements(1) = s
       avals(1) = wi(i)
       ok = gdxdatawritestr(ap2,astrelements,avals)
       if (.not. ok) call gdxerr(ap2)
    end do

    ok = gdxdatawritedone(ap2)
    if (.not. ok) call gdxerr(ap2)


    ok = gdxdatawritestrstart(ap2,evecrname,'real eigenvectors',2,dt_par,0)
    if (.not. ok) call gdxerr(ap2)

    astrelements(1:10) = ' '
    avals = 0.0d0

    do i=1,n
       ok = gdxumuelget(ap1,index(i),s,umap)
       if (.not. ok) call gdxerr(ap1)
       astrelements(1) = s
       do j=1,n
          ok = gdxumuelget(ap1,index(j),s,umap)
          if (.not. ok) call gdxerr(ap1)
          astrelements(2) = s
          avals(1) = vr(i,j)
          ok = gdxdatawritestr(ap2,astrelements,avals)
          if (.not. ok) call gdxerr(ap2)
       end do
    end do

    ok = gdxdatawritedone(ap2)
    if (.not. ok) call gdxerr(ap2)

   ok = gdxdatawritestrstart(ap2,eveciname,'imaginary eigenvectors',2,dt_par,0)
    if (.not. ok) call gdxerr(ap2)

    astrelements(1:10) = ' '
    avals = 0.0d0

    do i=1,n
       ok = gdxumuelget(ap1,index(i),s,umap)
       if (.not. ok) call gdxerr(ap1)
       astrelements(1) = s
       do j=1,n
          ok = gdxumuelget(ap1,index(j),s,umap)
          if (.not. ok) call gdxerr(ap1)
          astrelements(2) = s
          avals(1) = vri(i,j)
          ok = gdxdatawritestr(ap2,astrelements,avals)
          if (.not. ok) call gdxerr(ap2)
       end do
    end do

    ok = gdxdatawritedone(ap2)
    if (.not. ok) call gdxerr(ap2)


    deallocate(index,stat=ierr)
    if (ierr /= 0) call err('Could not deallocate vector')

    rc = gdxclose(ap1)
    rc = gdxclose(ap2)

  end subroutine writeresults

end program eigenvector
     SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
     $                  LDVR, WORK, LWORK, INFO )
*
*  -- LAPACK driver routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1999
*     8-15-00:  Improve consistency of WS calculations (eca)
*
*     .. Scalar Arguments ..
      CHARACTER          JOBVL, JOBVR
      INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
     $                   WI( * ), WORK( * ), WR( * )
*     ..
*
*  Purpose
*  =======
*
*  DGEEV computes for an N-by-N real nonsymmetric matrix A, the
*  eigenvalues and, optionally, the left and/or right eigenvectors.
*
*  The right eigenvector v(j) of A satisfies
*                   A * v(j) = lambda(j) * v(j)
*  where lambda(j) is its eigenvalue.
*  The left eigenvector u(j) of A satisfies
*                u(j)**H * A = lambda(j) * u(j)**H
*  where u(j)**H denotes the conjugate transpose of u(j).
*
*  The computed eigenvectors are normalized to have Euclidean norm
*  equal to 1 and largest component real.
*
*  Arguments
*  =========
*
*  JOBVL   (input) CHARACTER*1
*          = 'N': left eigenvectors of A are not computed;
*          = 'V': left eigenvectors of A are computed.
*
*  JOBVR   (input) CHARACTER*1
*          = 'N': right eigenvectors of A are not computed;
*          = 'V': right eigenvectors of A are computed.
*
*  N       (input) INTEGER
*          The order of the matrix A. N >= 0.
*
*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
*          On entry, the N-by-N matrix A.
*          On exit, A has been overwritten.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,N).
*
*  WR      (output) DOUBLE PRECISION array, dimension (N)
*  WI      (output) DOUBLE PRECISION array, dimension (N)
*          WR and WI contain the real and imaginary parts,
*          respectively, of the computed eigenvalues.  Complex
*          conjugate pairs of eigenvalues appear consecutively
*          with the eigenvalue having the positive imaginary part
*          first.
*
*  VL      (output) DOUBLE PRECISION array, dimension (LDVL,N)
*          If JOBVL = 'V', the left eigenvectors u(j) are stored one
*          after another in the columns of VL, in the same order
*          as their eigenvalues.
*          If JOBVL = 'N', VL is not referenced.
*          If the j-th eigenvalue is real, then u(j) = VL(:,j),
*          the j-th column of VL.
*          If the j-th and (j+1)-st eigenvalues form a complex
*          conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
*          u(j+1) = VL(:,j) - i*VL(:,j+1).
*
*  LDVL    (input) INTEGER
*          The leading dimension of the array VL.  LDVL >= 1; if
*          JOBVL = 'V', LDVL >= N.
*
*  VR      (output) DOUBLE PRECISION array, dimension (LDVR,N)
*          If JOBVR = 'V', the right eigenvectors v(j) are stored one
*          after another in the columns of VR, in the same order
*          as their eigenvalues.
*          If JOBVR = 'N', VR is not referenced.
*          If the j-th eigenvalue is real, then v(j) = VR(:,j),
*          the j-th column of VR.
*          If the j-th and (j+1)-st eigenvalues form a complex
*          conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
*          v(j+1) = VR(:,j) - i*VR(:,j+1).
*
*  LDVR    (input) INTEGER
*          The leading dimension of the array VR.  LDVR >= 1; if
*          JOBVR = 'V', LDVR >= N.
*
*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK.  LWORK >= max(1,3*N), and
*          if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N.  For good
*          performance, LWORK must generally be larger.
*
*          If LWORK = -1, a workspace query is assumed.  The optimal
*          size for the WORK array is calculated and stored in WORK(1),
*          and no other work except argument checking is performed.
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
*          > 0:  if INFO = i, the QR algorithm failed to compute all the
*                eigenvalues, and no eigenvectors have been computed;
*                elements i+1:N of WR and WI contain eigenvalues which
*                have converged.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            LQUERV
      PARAMETER          ( LQUERV = -1 )
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            SCALEA, WANTVL, WANTVR
      CHARACTER          SIDE
      INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
     $                   MAXB, MAXWRK, MINWRK, NOUT
      DOUBLE PRECISION   ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
     $                   SN
*     ..
*     .. Local Arrays ..
      LOGICAL            SELECT( 1 )
      DOUBLE PRECISION   DUM( 1 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY,
     $                   DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC,
     $                   XERBLA
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            IDAMAX, ILAENV
      DOUBLE PRECISION   DLAMCH, DLANGE, DLAPY2, DNRM2
      EXTERNAL           LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2,
     $                   DNRM2
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN, SQRT
*     ..
*     .. Executable Statements ..
*
*     Test the input arguments
*
      INFO = 0
      WANTVL = LSAME( JOBVL, 'V' )
      WANTVR = LSAME( JOBVR, 'V' )
      IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
         INFO = -1
      ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
         INFO = -9
      ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
         INFO = -11
      END IF
*
*     Compute workspace
*      (Note: Comments in the code beginning "Workspace:" describe the
*       minimal amount of workspace needed at that point in the code,
*       as well as the preferred amount for good performance.
*       NB refers to the optimal block size for the immediately
*       following subroutine, as returned by ILAENV.
*       HSWORK refers to the workspace preferred by DHSEQR, as
*       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
*       the worst case.)
*
      MINWRK = 1
      IF( INFO.EQ.0 ) THEN
         MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
         IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN
            MINWRK = MAX( 1, 3*N )
            MAXB = MAX( ILAENV( 8, 'DHSEQR', 'EN', N, 1, N, -1 ), 2 )
            K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'DHSEQR', 'EN', N, 1,
     $          N, -1 ) ) )
            HSWORK = MAX( K*( K+2 ), 2*N )
            MAXWRK = MAX( MAXWRK, N+1, N+HSWORK )
         ELSE
            MINWRK = MAX( 1, 4*N )
            MAXWRK = MAX( MAXWRK, 2*N+( N-1 )*
     $               ILAENV( 1, 'DORGHR', ' ', N, 1, N, -1 ) )
            MAXB = MAX( ILAENV( 8, 'DHSEQR', 'SV', N, 1, N, -1 ), 2 )
            K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'DHSEQR', 'SV', N, 1,
     $          N, -1 ) ) )
            HSWORK = MAX( K*( K+2 ), 2*N )
            MAXWRK = MAX( MAXWRK, N+1, N+HSWORK )
            MAXWRK = MAX( MAXWRK, 4*N )
         END IF
         WORK( 1 ) = MAXWRK
         IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV )
     $      INFO = -13
      END IF
*
*     Quick returns
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'DGEEV ', -INFO )
         RETURN
      END IF
      IF( LWORK.EQ.LQUERV )
     $   RETURN
      IF( N.EQ.0 )
     $   RETURN
*
*     Get machine constants
*
      EPS = DLAMCH( 'P' )
      SMLNUM = DLAMCH( 'S' )
      BIGNUM = ONE / SMLNUM
      CALL DLABAD( SMLNUM, BIGNUM )
      SMLNUM = SQRT( SMLNUM ) / EPS
      BIGNUM = ONE / SMLNUM
*
*     Scale A if max element outside range [SMLNUM,BIGNUM]
*
      ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
      SCALEA = .FALSE.
      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
         SCALEA = .TRUE.
         CSCALE = SMLNUM
      ELSE IF( ANRM.GT.BIGNUM ) THEN
         SCALEA = .TRUE.
         CSCALE = BIGNUM
      END IF
      IF( SCALEA )
     $   CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
*
*     Balance the matrix
*     (Workspace: need N)
*
      IBAL = 1
      CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
*
*     Reduce to upper Hessenberg form
*     (Workspace: need 3*N, prefer 2*N+N*NB)
*
      ITAU = IBAL + N
      IWRK = ITAU + N
      CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
     $             LWORK-IWRK+1, IERR )
*
      IF( WANTVL ) THEN
*
*        Want left eigenvectors
*        Copy Householder vectors to VL
*
         SIDE = 'L'
         CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL )
*
*        Generate orthogonal matrix in VL
*        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
*
         CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
     $                LWORK-IWRK+1, IERR )
*
*        Perform QR iteration, accumulating Schur vectors in VL
*        (Workspace: need N+1, prefer N+HSWORK (see comments) )
*
         IWRK = ITAU
         CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL,
     $                WORK( IWRK ), LWORK-IWRK+1, INFO )
*
         IF( WANTVR ) THEN
*
*           Want left and right eigenvectors
*           Copy Schur vectors to VR
*
            SIDE = 'B'
            CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
         END IF
*
      ELSE IF( WANTVR ) THEN
*
*        Want right eigenvectors
*        Copy Householder vectors to VR
*
         SIDE = 'R'
         CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR )
*
*        Generate orthogonal matrix in VR
*        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
*
         CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
     $                LWORK-IWRK+1, IERR )
*
*        Perform QR iteration, accumulating Schur vectors in VR
*        (Workspace: need N+1, prefer N+HSWORK (see comments) )
*
         IWRK = ITAU
         CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
     $                WORK( IWRK ), LWORK-IWRK+1, INFO )
*
      ELSE
*
*        Compute eigenvalues only
*        (Workspace: need N+1, prefer N+HSWORK (see comments) )
*
         IWRK = ITAU
         CALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR,
     $                WORK( IWRK ), LWORK-IWRK+1, INFO )
      END IF
*
*     If INFO > 0 from DHSEQR, then quit
*
      IF( INFO.GT.0 )
     $   GO TO 50
*
      IF( WANTVL .OR. WANTVR ) THEN
*
*        Compute left and/or right eigenvectors
*        (Workspace: need 4*N)
*
         CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
     $                N, NOUT, WORK( IWRK ), IERR )
      END IF
*
      IF( WANTVL ) THEN
*
*        Undo balancing of left eigenvectors
*        (Workspace: need N)
*
         CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL,
     $                IERR )
*
*        Normalize left eigenvectors and make largest component real
*
         DO 20 I = 1, N
            IF( WI( I ).EQ.ZERO ) THEN
               SCL = ONE / DNRM2( N, VL( 1, I ), 1 )
               CALL DSCAL( N, SCL, VL( 1, I ), 1 )
            ELSE IF( WI( I ).GT.ZERO ) THEN
               SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ),
     $               DNRM2( N, VL( 1, I+1 ), 1 ) )
               CALL DSCAL( N, SCL, VL( 1, I ), 1 )
               CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 )
               DO 10 K = 1, N
                  WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2
   10          CONTINUE
               K = IDAMAX( N, WORK( IWRK ), 1 )
               CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R )
               CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN )
               VL( K, I+1 ) = ZERO
            END IF
   20    CONTINUE
      END IF
*
      IF( WANTVR ) THEN
*
*        Undo balancing of right eigenvectors
*        (Workspace: need N)
*
         CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR,
     $                IERR )
*
*        Normalize right eigenvectors and make largest component real
*
         DO 40 I = 1, N
            IF( WI( I ).EQ.ZERO ) THEN
               SCL = ONE / DNRM2( N, VR( 1, I ), 1 )
               CALL DSCAL( N, SCL, VR( 1, I ), 1 )
            ELSE IF( WI( I ).GT.ZERO ) THEN
               SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ),
     $               DNRM2( N, VR( 1, I+1 ), 1 ) )
               CALL DSCAL( N, SCL, VR( 1, I ), 1 )
               CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 )
               DO 30 K = 1, N
                  WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2
   30          CONTINUE
               K = IDAMAX( N, WORK( IWRK ), 1 )
               CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R )
               CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN )
               VR( K, I+1 ) = ZERO
            END IF
   40    CONTINUE
      END IF
*
*     Undo scaling if necessary
*
   50 CONTINUE
      IF( SCALEA ) THEN
         CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ),
     $                MAX( N-INFO, 1 ), IERR )
         CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ),
     $                MAX( N-INFO, 1 ), IERR )
         IF( INFO.GT.0 ) THEN
            CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N,
     $                   IERR )
            CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
     $                   IERR )
         END IF
      END IF
*
      WORK( 1 ) = MAXWRK
      RETURN
*
*     End of DGEEV
*
      END


[Julien Langou]

Hello
in routine calc_nonsymmetric
I think you messed up your call to dgeev
what about replacing nwork by work ?
-j

[traptij]

Hello,
Thank you, now it is working.

Thanks,
T.Jain