Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

added a routine to use DSTEGR #2

Merged
merged 1 commit into from
Apr 5, 2018
Merged

added a routine to use DSTEGR #2

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
104 changes: 100 additions & 4 deletions src/dsyevr_check.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,11 +43,13 @@ module diagmod
! utils/diag.F
!
! Written by Miro Ilias, August 2014
! Vishnu V. Krishnan, April 2018
!-------------------------------------------------------------------------------------------
integer :: N=0,NZ=0, print_level
integer :: lu=12, LUPRI=6
real*8, allocatable :: Ar(:,:), ASr(:,:),Ai(:,:),ASi(:,:), &
Aj(:,:),ASj(:,:), Ak(:,:), ASk(:,:)
Aj(:,:),ASj(:,:), Ak(:,:), ASk(:,:), &
ArTemp(:,:)
character*50 :: title_text
! hardcoded input matrix name
character*50 :: matrix_file_name = "Jz_SS_matrix.fermirp2-2"
Expand Down Expand Up @@ -120,6 +122,9 @@ subroutine test_lapack_dsyevr
integer :: IWORKDUM,IDUM
real*8, parameter :: D0=0.0D0

allocate(ArTemp(N,N), stat=ierr)
ArTemp = Ar

ABSTOL=0.0d0; JOBZ='V'; RANGE = 'A'; UPLO = 'U'
DDUM = D0; IDUM = 0; WORKDUM = D0; IWORKDUM = 0

Expand All @@ -139,7 +144,7 @@ subroutine test_lapack_dsyevr
IERR = 0
! FYI: LLWORK=-1; LIWORK=-1 on input means calculate temporary storage
M = N ; LRA = N ;LRV = N
CALL DSYEVR(JOBZ,'A','U',N,Ar,LRA,DDUM,DDUM,IDUM,IDUM,ABSTOL, &
CALL DSYEVR(JOBZ,'A','U',N,ArTemp,LRA,DDUM,DDUM,IDUM,IDUM,ABSTOL, &
M,EIG,VEC,LRV,IDUM,EIG,-1,IWORKDUM,-1,IERR)

LLWORK = NINT(EIG(1))
Expand All @@ -150,7 +155,7 @@ subroutine test_lapack_dsyevr
allocate(ISUPPZ(2*N),stat=ierr)

IERR = 0
CALL DSYEVR(JOBZ,'A','U',N,Ar,LRA,DDUM,DDUM,IDUM,IDUM,ABSTOL, &
CALL DSYEVR(JOBZ,'A','U',N,ArTemp,LRA,DDUM,DDUM,IDUM,IDUM,ABSTOL, &
M,EIG,VEC,LRV,ISUPPZ,WORK,LLWORK,IWORK,LILWORK,IERR)
if (IERR /= 0) then
print *,'The lapack_dsyerv routine ended with error ! ierr=',ierr
Expand All @@ -166,10 +171,100 @@ subroutine test_lapack_dsyevr

call eigv_check(EIG,VEC,'**** LAPACK DSYEVR ****')

deallocate(WORK,IWORK,EIG,VEC,ISUPPZ,stat=ierr)
deallocate(ArTemp,WORK,IWORK,EIG,VEC,ISUPPZ,stat=ierr)

end subroutine test_lapack_dsyevr

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
subroutine test_lapack_dstegr
!
! tests lapack's routine DSTEGR for real symmetric matrices
! after tri-diagonalising with DSYTRD
!
implicit none
integer :: ierr, MATZ, LLWORK, KLWORK, LILWORK
integer :: i
integer :: IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M, LRA, LRV
character*1 :: JOBZ, RANGE, UPLO
real*8 :: ABSTOL, VL, VU
integer, allocatable :: ISUPPZ(:), IWORK(:)
real*8, allocatable :: WORK(:), EIG(:),VEC(:,:),DIAG(:),OFFDIAG(:), &
RFLCT(:)
real*8 :: DDUM, WORKDUM
integer :: IWORKDUM,IDUM
real*8, parameter :: D0=0.0D0

allocate(ArTemp(N,N), stat=ierr)
ArTemp = Ar

ABSTOL=0.0d0; JOBZ='V'; RANGE = 'A'; UPLO = 'U'
DDUM = D0; IDUM = 0; WORKDUM = D0; IWORKDUM = 0

M = N ; LRA = N ;LRV = N

allocate(DIAG(N),OFFDIAG(N),RFLCT(N-1))

IERR = 0
CALL DSYTRD(UPLO,N,ArTemp,LRA,DIAG,OFFDIAG(:N-1),RFLCT,WORKDUM,-1,IERR)

LLWORK = NINT(WORKDUM)
allocate(WORK(LLWORK))

IERR = 0
CALL DSYTRD(UPLO,N,ArTemp,LRA,DIAG,OFFDIAG(:N-1),RFLCT,WORK,LLWORK,IERR)
if (IERR /= 0) then
print *,'The lapack_dsytrd routine ended with error ! ierr=',ierr
endif

deallocate(RFLCT,WORK)

MATZ = 1
IF(MATZ.EQ.1) THEN
JOBZ='V'
ELSEIF(MATZ.EQ.0) THEN
JOBZ='N'
ELSE
WRITE(*,*) 'Illegal value of MATZ: ',MATZ
STOP 2
ENDIF

allocate(EIG(N),VEC(N,N))

WORKDUM = D0
! 1.call to determine size
IERR = 0
! FYI: LLWORK=-1; LIWORK=-1 on input means calculate temporary storage
CALL DSTEGR(JOBZ,'A',N,DIAG,OFFDIAG,DDUM,DDUM,IDUM,IDUM,ABSTOL, &
M,EIG,VEC,LRV,IDUM,WORKDUM,-1,IWORKDUM,-1,IERR)

LLWORK = NINT(WORKDUM)
LILWORK = IWORKDUM

allocate(WORK(LLWORK), stat=ierr)
allocate(IWORK(LILWORK), stat=ierr)
allocate(ISUPPZ(2*N),stat=ierr)

IERR = 0
CALL DSTEGR(JOBZ,'A',N,DIAG,OFFDIAG,DDUM,DDUM,IDUM,IDUM,ABSTOL, &
M,EIG,VEC,LRV,ISUPPZ,WORK,LLWORK,IWORK,LILWORK,IERR)
if (IERR /= 0) then
print *,'The lapack_dstegr routine ended with error ! ierr=',ierr
endif

! ... print out eigenvalues
if (print_level >= 0) then
print *,"LAPACK DSTEGR eigenvalues:"
do i=1,N
print *,i,EIG(i)
enddo
endif

call eigv_check(EIG,VEC,'**** LAPACK DSYTRD+DSTEGR ****')

deallocate(ArTemp,DIAG,OFFDIAG,WORK,IWORK,EIG,VEC,ISUPPZ)

end subroutine test_lapack_dstegr

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
subroutine test_dirac_rs
! tests DIRAC's EISPACK routine RS real Hermitian matrix
Expand Down Expand Up @@ -293,5 +388,6 @@ Program Test_DIRAC_Diagonalization_Routines
call read_matrix
call test_dirac_rs
call test_lapack_dsyevr
call test_lapack_dstegr

End Program Test_DIRAC_Diagonalization_Routines