Open discussion regarding features, bugs, issues, vendors, etc.
by traptij » Mon May 22, 2006 5:47 am
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
-
traptij
-
- Posts: 11
- Joined: Sat May 13, 2006 7:21 am
- Location: IIT Kanpur
-
by traptij » Tue May 23, 2006 12:33 am
Hello,
No, it is not working, giving the same error.
Thanks,
T.Jain
-
traptij
-
- Posts: 11
- Joined: Sat May 13, 2006 7:21 am
- Location: IIT Kanpur
-
by Julien Langou » Tue May 23, 2006 6:52 am
Hello,
you will need to post a piece of code if you want more help, all this should work then.
Julien
-
Julien Langou
-
- Posts: 835
- Joined: Thu Dec 09, 2004 12:32 pm
- Location: Denver, CO, USA
-
by traptij » Tue May 23, 2006 8:05 am
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
-
traptij
-
- Posts: 11
- Joined: Sat May 13, 2006 7:21 am
- Location: IIT Kanpur
-
by Julien Langou » Tue May 23, 2006 10:42 am
Hello
in routine calc_nonsymmetric
I think you messed up your call to dgeev
what about replacing nwork by work ?
-j
-
Julien Langou
-
- Posts: 835
- Joined: Thu Dec 09, 2004 12:32 pm
- Location: Denver, CO, USA
-
by traptij » Wed May 24, 2006 7:47 am
Hello,
Thank you, now it is working.
Thanks,
T.Jain
-
traptij
-
- Posts: 11
- Joined: Sat May 13, 2006 7:21 am
- Location: IIT Kanpur
-
Return to User Discussion
Who is online
Users browsing this forum: No registered users and 1 guest