module m_dgeevx
!$$$ subprogram documentation block
! . . . .
! subprogram: module m_dgeevx
! prgmmr: j guo <jguo@nasa.gov>
! org: NASA/GSFC, Global Modeling and Assimilation Office, 900.3
! date: 2010-03-24
!
! abstract: an alternative interface of LAPACK DGEEV()
!
! program history log:
! 2010-03-24 j guo - added this document block
!
! input argument list: see Fortran 90 style document below
!
! output argument list: see Fortran 90 style document below
!
! attributes:
! language: Fortran 90 and/or above
! machine:
!
!$$$ end subprogram documentation block
! module interface:
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS !
!BOP -------------------------------------------------------------------
!
! !MODULE: m_dgeevx - an alternative interface of LAPACK DGEEV()
!
! !DESCRIPTION:
!
! !INTERFACE:
implicit none
private ! except
public :: dgeevx ! alternative DGEEV()
! !REVISION HISTORY:
! 29May08 - Jing Guo <guo@gmao.gsfc.nasa.gov>
! - initial prototype/prolog/code
!EOP ___________________________________________________________________
character(len=*),parameter :: myname='m_dgeevx'
contains
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NASA/GSFC, Global Modeling and Assimilation Office, 900.3, GEOS/DAS !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: dgeevx - an alternative interface of DGEEV
!
! !DESCRIPTION:
!
! !INTERFACE:
subroutine dgeevx(nsig,qmat,ldqmat,www,wwwd,zzz,ldzzz,zzzd,ldzzzd,info,mype)
use kinds,only : i_kind
use kinds,only : r_double,r_kind
use constants,only : zero
implicit none
integer(i_kind) ,intent(in ) :: nsig ! dimension matrix A (n-by-n (n=nsig))
integer(i_kind) ,intent(in ) :: ldqmat ! leading-dimension of qmat
real(r_kind) ,dimension(ldqmat ,nsig),intent(in ) :: qmat ! matrix A
real(r_double),dimension(2,nsig) ,intent( out) :: www ! right-eigen-values
real(r_double),dimension(2,nsig) ,intent( out) :: wwwd ! left-eigen-values
integer(i_kind) ,intent(in ) :: ldzzz ! leading-dimension of zzz
real(r_double),dimension(2,ldzzz ,nsig),intent( out) :: zzz ! right-eigen-vectors
integer(i_kind) ,intent(in ) :: ldzzzd ! leading-dimension of zzzd
real(r_double),dimension(2,ldzzzd,nsig),intent( out) :: zzzd ! left-eigen-vectors
integer(i_kind) ,intent( out) :: info ! return status
integer(i_kind) ,intent(in ) :: mype ! PE rank for diagnostics, -1 to turnoff
! !REVISION HISTORY:
! 29May08 - Jing Guo <guo@gmao.gsfc.nasa.gov>
! - initial prototype/prolog/code
!EOP ___________________________________________________________________
character(len=*),parameter :: myname_=myname//'::dgeevx'
#ifdef ibm_sp
integer(i_kind) :: iopt,i,j,k
real(r_double) :: aaa(nsig,nsig)
logical :: select(nsig) ! an argument of dgeev(), but not used.
integer(i_kind) :: naux
real(r_kind) :: aux,factor
! Use IBM ESSL routine dgeev()
! Below is a copy of the original code in mod_vtrans.f90_1.2.
!<<<<
! next get eigenvalues, eigenvectors and compare to singular values, vectors.
! use essl routine dgeev
iopt=1 ! eigenvalues and eigenvectors are computed
do j=1,nsig
do i=1,nsig
aaa(i,j)=qmat(i,j)
end do
end do
naux=0
call dgeev(iopt,aaa,nsig,www,zzz,nsig,select,nsig,aux,naux)
! sort from largest to smallest eigenvalue
do j=1,nsig-1
do i=j+1,nsig
if(www(1,i)>www(1,j)) then
factor=www(1,j)
www(1,j)=www(1,i)
www(1,i)=factor
factor=www(2,j)
www(2,j)=www(2,i)
www(2,i)=factor
do k=1,nsig
factor=zzz(1,k,j)
zzz(1,k,j)=zzz(1,k,i)
zzz(1,k,i)=factor
factor=zzz(2,k,j)
zzz(2,k,j)=zzz(2,k,i)
zzz(2,k,i)=factor
end do
end if
end do
end do
! checks and print out eigenvalues (removed)
!
iopt=1 ! eigenvalues and dual eigenvectors are computed next
do j=1,nsig
do i=1,nsig
aaa(i,j)=qmat(j,i) ! to get dual vectors, use transpose of qmat
end do
end do
naux=0
call dgeev(iopt,aaa,nsig,wwwd,zzzd,nsig,select,nsig,aux,naux)
! sort from largest to smallest eigenvalue
do j=1,nsig-1
do i=j+1,nsig
if(wwwd(1,i)>wwwd(1,j)) then
factor=wwwd(1,j)
wwwd(1,j)=wwwd(1,i)
wwwd(1,i)=factor
factor=wwwd(2,j)
wwwd(2,j)=wwwd(2,i)
wwwd(2,i)=factor
do k=1,nsig
factor=zzzd(1,k,j)
zzzd(1,k,j)=zzzd(1,k,i)
zzzd(1,k,i)=factor
factor=zzzd(2,k,j)
zzzd(2,k,j)=zzzd(2,k,i)
zzzd(2,k,i)=factor
end do
end if
end do
end do
info=0
#else
! Use standard LAPACK routine dgeev()
!----------------------------------------
! Modification for using GSI dgeev -RY
!----------------------------------------
character*1,parameter:: jobvr='V'
character*1,parameter:: jobvl='V'
! double precision
real(r_double) :: aaa(nsig,nsig)
real(r_double) :: wr(nsig),wi(nsig)
real(r_double) :: vl(nsig,nsig),vr(nsig,nsig)
real(r_double) :: work(4*nsig+1)
integer(i_kind) :: lwork,i,j,k
real(r_kind) :: factor,factor2
! size of work(:), used by dgeev subroutine
lwork=size(work)
!------------------------------------
! ryang:
! use SGI scslib dgeev subroutine
do j=1,nsig
do i=1,nsig
aaa(i,j)=qmat(i,j)
end do
end do
if (mype ==0) write (6,*) myname, ' before CALL DGEEV'
call dgeev(jobvl,jobvr,nsig,aaa,nsig,wr,wi,vl,nsig,vr,nsig,work,lwork,info)
if (mype ==0) write (6,*) myname, ' after CALL DGEEV', 'status: info = ', info
! use Dave's array names
do j=1,nsig
if (wi(j) /= zero) write (6,*) &
'wrong eigen computation: create_vtrans'
www(1,j)=wr(j)
www(2,j)=wi(j)
!-----------------------------------------------------------
!ADD explanation:!!!
!
! vl is zzzd
!the following line is not generic, only hold when wi (j)=0.0
!eigenvalues of A^transpose are the same as A
!----------------------------------------------
do k=1,nsig
zzz (1,k,j)=vr(k,j)
zzzd(1,k,j)=vl(k,j)
zzz (2,k,j)=zero
zzzd(2,k,j)=zero
enddo
enddo
! back to Dave's code
!sort from largest to smallest eigenvalues
! sort from largest to smallest eigenvalue
do j=1,nsig-1
do i=j+1,nsig
if(www(1,i)>www(1,j)) then
factor=www(1,j)
www(1,j)=www(1,i)
www(1,i)=factor
factor=www(2,j)
www(2,j)=www(2,i)
www(2,i)=factor
do k=1,nsig
factor =zzz (1,k,j)
factor2=zzzd(1,k,j)
zzz (1,k,j)=zzz (1,k,i)
zzzd(1,k,j)=zzzd(1,k,i)
zzz (1,k,i)=factor
zzzd(1,k,i)=factor2
factor =zzz (2,k,j)
factor2=zzzd(2,k,j)
zzz (2,k,j)=zzz (2,k,i)
zzzd(2,k,j)=zzzd(2,k,i)
zzz (2,k,i)=factor
zzzd(2,k,i)=factor2
end do
end if
end do
end do
! checks and print out eigenvalues (removed)
!
!**********************************************************
!*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! do not need to do the transpose
! check the eigenvalues
! using the component form
wwwd=www
if (mype==0) write (6,*) '****************************'
if (mype==0) write (6,*) 'SECOND TIME CALL DGEEV'
#endif
end subroutine dgeevx
end module m_dgeevx