Program f08wbfe
! F08WBF Example Program Text
! Mark 29.2 Release. NAG Copyright 2023.
! .. Use Statements ..
Use nag_library, Only: dggevx, dnrm2, m01daf, m01eaf, nag_wp, x02ajf, &
x02amf
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: nb = 64, nin = 5, nout = 6
Logical, Parameter :: verbose = .False.
! .. Local Scalars ..
Complex (Kind=nag_wp) :: eig
Real (Kind=nag_wp) :: abnrm, bbnrm, eps, jswap, rcnd, &
scal_i, scal_r, small
Integer :: i, ifail, ihi, ilo, info, j, k, lda, &
ldb, ldvr, lwork, n
Logical :: pair
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: a(:,:), alphai(:), alphar(:), &
b(:,:), beta(:), lscale(:), &
rconde(:), rcondv(:), rscale(:), &
vr(:,:), vr_row(:), work(:)
Real (Kind=nag_wp) :: dummy(1,1)
Integer, Allocatable :: iwork(:)
Logical, Allocatable :: bwork(:)
! .. Intrinsic Procedures ..
Intrinsic :: abs, all, cmplx, max, maxloc, nint, &
sqrt, sum
! .. Executable Statements ..
Write (nout,*) 'F08WBF Example Program Results'
! Skip heading in data file
Read (nin,*)
Read (nin,*) n
lda = n
ldb = n
ldvr = n
Allocate (a(lda,n),alphai(n),alphar(n),b(ldb,n),beta(n),lscale(n), &
rconde(n),rcondv(n),rscale(n),vr(ldvr,n),iwork(n+6),bwork(n))
! Use routine workspace query to get optimal workspace.
lwork = -1
! The NAG name equivalent of dggevx is f08wbf
Call dggevx('Balance','No vectors (left)','Vectors (right)', &
'Both reciprocal condition numbers',n,a,lda,b,ldb,alphar,alphai,beta, &
dummy,1,vr,ldvr,ilo,ihi,lscale,rscale,abnrm,bbnrm,rconde,rcondv,dummy, &
lwork,iwork,bwork,info)
! Make sure that there is enough workspace for block size nb.
lwork = max((nb+2*n)*n,nint(dummy(1,1)))
Allocate (work(lwork))
! Read in the matrices A and B
Read (nin,*)(a(i,1:n),i=1,n)
Read (nin,*)(b(i,1:n),i=1,n)
! Solve the generalized eigenvalue problem
! The NAG name equivalent of dggevx is f08wbf
Call dggevx('Balance','No vectors (left)','Vectors (right)', &
'Both reciprocal condition numbers',n,a,lda,b,ldb,alphar,alphai,beta, &
dummy,1,vr,ldvr,ilo,ihi,lscale,rscale,abnrm,bbnrm,rconde,rcondv,work, &
lwork,iwork,bwork,info)
If (info>0) Then
Write (nout,*)
Write (nout,99999) 'Failure in DGGEVX. INFO =', info
Else
! Compute the machine precision and the safe range parameter
! small
eps = x02ajf()
small = x02amf()
! If beta(:) > eps, Order eigenvalues by ascending real parts
If (all(abs(beta(1:n))>eps)) Then
work(1:n) = alphar(1:n)/beta(1:n)
ifail = 0
Call m01daf(work,1,n,'Ascending',iwork,ifail)
Call m01eaf(alphar,1,n,iwork,ifail)
Call m01eaf(alphai,1,n,iwork,ifail)
Call m01eaf(beta,1,n,iwork,ifail)
! Order the eigenvectors in the same way
Allocate (vr_row(n))
Do j = 1, n
vr_row(1:n) = vr(j,1:n)
Call m01eaf(vr_row,1,n,iwork,ifail)
vr(j,1:n) = vr_row(1:n)
End Do
Deallocate (vr_row)
End If
! Print out eigenvalues and vectors and associated condition
! number and bounds
pair = .False.
Do j = 1, n
! Print out information on the j-th eigenvalue
Write (nout,*)
If ((abs(alphar(j))+abs(alphai(j)))*small>=abs(beta(j))) Then
Write (nout,99998) 'Eigenvalue(', j, ')', &
' is numerically infinite or undetermined', 'ALPHAR(', j, &
') = ', alphar(j), ', ALPHAI(', j, ') = ', alphai(j), ', BETA(', &
j, ') = ', beta(j)
Else
If (.Not. pair) Then
jswap = 1.0_nag_wp
If (alphai(j)>0.0_nag_wp) Then
jswap = -jswap
End If
End If
If (alphai(j)==0.0E0_nag_wp) Then
Write (nout,99997) 'Eigenvalue(', j, ') = ', alphar(j)/beta(j)
Else
eig = cmplx(alphar(j),jswap*alphai(j),kind=nag_wp)/ &
cmplx(beta(j),kind=nag_wp)
Write (nout,99996) 'Eigenvalue(', j, ') = ', eig
End If
End If
If (verbose) Then
rcnd = rconde(j)
Write (nout,*)
Write (nout,99995) ' Reciprocal condition number = ', rcnd
End If
! Print out information on the j-th eigenvector
Write (nout,*)
Write (nout,99994) 'Eigenvector(', j, ')'
If (alphai(j)==0.0E0_nag_wp) Then
! Let largest element be positive
work(1:n) = abs(vr(1:n,j))
k = maxloc(work(1:n),1)
If (vr(k,j)<0.0_nag_wp) Then
vr(1:n,j) = -vr(1:n,j)/dnrm2(n,vr(1,j),1)
End If
Write (nout,99993)(vr(i,j),i=1,n)
Else
If (pair) Then
Write (nout,99992)(vr(i,j-1),-jswap*vr(i,j),i=1,n)
Else
! Let largest element be real (and positive).
work(1:n) = vr(1:n,j)**2 + vr(1:n,j+1)**2
k = maxloc(work(1:n),1)
scal_r = vr(k,j)/sqrt(work(k))/sqrt(sum(work(1:n)))
scal_i = -vr(k,j+1)/sqrt(work(k))/sqrt(sum(work(1:n)))
work(1:n) = vr(1:n,j)
vr(1:n,j) = scal_r*work(1:n) - scal_i*vr(1:n,j+1)
vr(1:n,j+1) = scal_r*vr(1:n,j+1) + scal_i*work(1:n)
vr(k,j+1) = 0.0_nag_wp
Write (nout,99992)(vr(i,j),jswap*vr(i,j+1),i=1,n)
End If
pair = .Not. pair
End If
If (verbose) Then
rcnd = rcondv(j)
Write (nout,*)
Write (nout,99995) ' Reciprocal condition number = ', rcnd
End If
End Do
End If
99999 Format (1X,A,I4)
99998 Format (/,1X,A,I2,2A,/,1X,2(A,I2,A,F11.5),A,I2,A,F11.5)
99997 Format (/,1X,A,I2,A,F11.5)
99996 Format (/,1X,A,I2,A,'(',F11.5,',',F11.5,')')
99995 Format (1X,A,1P,E8.1)
99994 Format (1X,A,I2,A)
99993 Format (1X,F11.5)
99992 Format (1X,'(',F11.5,',',F11.5,')')
End Program f08wbfe