c\BeginDoc
c
c\Name: cneupd
c
c\Description:
c This subroutine returns the converged approximations to eigenvalues
c of A*z = lambda*B*z and (optionally):
c
c (1) The corresponding approximate eigenvectors;
c
c (2) An orthonormal basis for the associated approximate
c invariant subspace;
c
c (3) Both.
c
c There is negligible additional cost to obtain eigenvectors. An orthonormal
c basis is always computed. There is an additional storage cost of n*nev
c if both are requested (in this case a separate array Z must be supplied).
c
c The approximate eigenvalues and eigenvectors of A*z = lambda*B*z
c are derived from approximate eigenvalues and eigenvectors of
c of the linear operator OP prescribed by the MODE selection in the
c call to CNAUPD. CNAUPD must be called before this routine is called.
c These approximate eigenvalues and vectors are commonly called Ritz
c values and Ritz vectors respectively. They are referred to as such
c in the comments that follow. The computed orthonormal basis for the
c invariant subspace corresponding to these Ritz values is referred to as a
c Schur basis.
c
c The definition of OP as well as other terms and the relation of computed
c Ritz values and vectors of OP with respect to the given problem
c A*z = lambda*B*z may be found in the header of CNAUPD. For a brief
c description, see definitions of IPARAM(7), MODE and WHICH in the
c documentation of CNAUPD.
c
c\Usage:
c call cneupd
c ( RVEC, HOWMNY, SELECT, D, Z, LDZ, SIGMA, WORKEV, BMAT,
c N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR, WORKD,
c WORKL, LWORKL, RWORK, INFO )
c
c\Arguments:
c RVEC LOGICAL (INPUT)
c Specifies whether a basis for the invariant subspace corresponding
c to the converged Ritz value approximations for the eigenproblem
c A*z = lambda*B*z is computed.
c
c RVEC = .FALSE. Compute Ritz values only.
c
c RVEC = .TRUE. Compute Ritz vectors or Schur vectors.
c See Remarks below.
c
c HOWMNY Character*1 (INPUT)
c Specifies the form of the basis for the invariant subspace
c corresponding to the converged Ritz values that is to be computed.
c
c = 'A': Compute NEV Ritz vectors;
c = 'P': Compute NEV Schur vectors;
c = 'S': compute some of the Ritz vectors, specified
c by the logical array SELECT.
c
c SELECT Logical array of dimension NCV. (INPUT)
c If HOWMNY = 'S', SELECT specifies the Ritz vectors to be
c computed. To select the Ritz vector corresponding to a
c Ritz value D(j), SELECT(j) must be set to .TRUE..
c If HOWMNY = 'A' or 'P', SELECT need not be initialized
c but it is used as internal workspace.
c
c D Complex array of dimension NEV+1. (OUTPUT)
c On exit, D contains the Ritz approximations
c to the eigenvalues lambda for A*z = lambda*B*z.
c
c Z Complex N by NEV array (OUTPUT)
c On exit, if RVEC = .TRUE. and HOWMNY = 'A', then the columns of
c Z represents approximate eigenvectors (Ritz vectors) corresponding
c to the NCONV=IPARAM(5) Ritz values for eigensystem
c A*z = lambda*B*z.
c
c If RVEC = .FALSE. or HOWMNY = 'P', then Z is NOT REFERENCED.
c
c NOTE: If if RVEC = .TRUE. and a Schur basis is not required,
c the array Z may be set equal to first NEV+1 columns of the Arnoldi
c basis array V computed by CNAUPD. In this case the Arnoldi basis
c will be destroyed and overwritten with the eigenvector basis.
c
c LDZ Integer. (INPUT)
c The leading dimension of the array Z. If Ritz vectors are
c desired, then LDZ .ge. max( 1, N ) is required.
c In any case, LDZ .ge. 1 is required.
c
c SIGMA Complex (INPUT)
c If IPARAM(7) = 3 then SIGMA represents the shift.
c Not referenced if IPARAM(7) = 1 or 2.
c
c WORKEV Complex work array of dimension 2*NCV. (WORKSPACE)
c
c **** The remaining arguments MUST be the same as for the ****
c **** call to CNAUPD that was just completed. ****
c
c NOTE: The remaining arguments
c
c BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR,
c WORKD, WORKL, LWORKL, RWORK, INFO
c
c must be passed directly to CNEUPD following the last call
c to CNAUPD. These arguments MUST NOT BE MODIFIED between
c the the last call to CNAUPD and the call to CNEUPD.
c
c Three of these parameters (V, WORKL and INFO) are also output parameters:
c
c V Complex N by NCV array. (INPUT/OUTPUT)
c
c Upon INPUT: the NCV columns of V contain the Arnoldi basis
c vectors for OP as constructed by CNAUPD .
c
c Upon OUTPUT: If RVEC = .TRUE. the first NCONV=IPARAM(5) columns
c contain approximate Schur vectors that span the
c desired invariant subspace.
c
c NOTE: If the array Z has been set equal to first NEV+1 columns
c of the array V and RVEC=.TRUE. and HOWMNY= 'A', then the
c Arnoldi basis held by V has been overwritten by the desired
c Ritz vectors. If a separate array Z has been passed then
c the first NCONV=IPARAM(5) columns of V will contain approximate
c Schur vectors that span the desired invariant subspace.
c
c WORKL Real work array of length LWORKL. (OUTPUT/WORKSPACE)
c WORKL(1:ncv*ncv+2*ncv) contains information obtained in
c cnaupd. They are not changed by cneupd.
c WORKL(ncv*ncv+2*ncv+1:3*ncv*ncv+4*ncv) holds the
c untransformed Ritz values, the untransformed error estimates of
c the Ritz values, the upper triangular matrix for H, and the
c associated matrix representation of the invariant subspace for H.
c
c Note: IPNTR(9:13) contains the pointer into WORKL for addresses
c of the above information computed by cneupd.
c -------------------------------------------------------------
c IPNTR(9): pointer to the NCV RITZ values of the
c original system.
c IPNTR(10): Not used
c IPNTR(11): pointer to the NCV corresponding error estimates.
c IPNTR(12): pointer to the NCV by NCV upper triangular
c Schur matrix for H.
c IPNTR(13): pointer to the NCV by NCV matrix of eigenvectors
c of the upper Hessenberg matrix H. Only referenced by
c cneupd if RVEC = .TRUE. See Remark 2 below.
c -------------------------------------------------------------
c
c INFO Integer. (OUTPUT)
c Error flag on output.
c = 0: Normal exit.
c
c = 1: The Schur form computed by LAPACK routine csheqr
c could not be reordered by LAPACK routine ctrsen.
c Re-enter subroutine cneupd with IPARAM(5)=NCV and
c increase the size of the array D to have
c dimension at least dimension NCV and allocate at least NCV
c columns for Z. NOTE: Not necessary if Z and V share
c the same space. Please notify the authors if this error
c occurs.
c
c = -1: N must be positive.
c = -2: NEV must be positive.
c = -3: NCV-NEV >= 2 and less than or equal to N.
c = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'
c = -6: BMAT must be one of 'I' or 'G'.
c = -7: Length of private work WORKL array is not sufficient.
c = -8: Error return from LAPACK eigenvalue calculation.
c This should never happened.
c = -9: Error return from calculation of eigenvectors.
c Informational error from LAPACK routine ctrevc.
c = -10: IPARAM(7) must be 1,2,3
c = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
c = -12: HOWMNY = 'S' not yet implemented
c = -13: HOWMNY must be one of 'A' or 'P' if RVEC = .true.
c = -14: CNAUPD did not find any eigenvalues to sufficient
c accuracy.
c
c\BeginLib
c
c\References:
c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
c pp 357-385.
c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
c Restarted Arnoldi Iteration", Rice University Technical Report
c TR95-13, Department of Computational and Applied Mathematics.
c 3. B. Nour-Omid, B. N. Parlett, T. Ericsson and P. S. Jensen,
c "How to Implement the Spectral Transformation", Math Comp.,
c Vol. 48, No. 178, April, 1987 pp. 664-673.
c
c\Routines called:
c ivout ARPACK utility routine that prints integers.
c cmout ARPACK utility routine that prints matrices
c cvout ARPACK utility routine that prints vectors.
c cgeqr2 LAPACK routine that computes the QR factorization of
c a matrix.
c clacpy LAPACK matrix copy routine.
c clahqr LAPACK routine that computes the Schur form of a
c upper Hessenberg matrix.
c claset LAPACK matrix initialization routine.
c ctrevc LAPACK routine to compute the eigenvectors of a matrix
c in upper triangular form.
c ctrsen LAPACK routine that re-orders the Schur form.
c cunm2r LAPACK routine that applies an orthogonal matrix in
c factored form.
c slamch LAPACK routine that determines machine constants.
c ctrmm Level 3 BLAS matrix times an upper triangular matrix.
c cgeru Level 2 BLAS rank one update to a matrix.
c ccopy Level 1 BLAS that copies one vector to another .
c cscal Level 1 BLAS that scales a vector.
c csscal Level 1 BLAS that scales a complex vector by a real number.
c scnrm2 Level 1 BLAS that computes the norm of a complex vector.
c
c\Remarks
c
c 1. Currently only HOWMNY = 'A' and 'P' are implemented.
c
c 2. Schur vectors are an orthogonal representation for the basis of
c Ritz vectors. Thus, their numerical properties are often superior.
c If RVEC = .true. then the relationship
c A * V(:,1:IPARAM(5)) = V(:,1:IPARAM(5)) * T, and
c V(:,1:IPARAM(5))' * V(:,1:IPARAM(5)) = I are approximately satisfied.
c Here T is the leading submatrix of order IPARAM(5) of the
c upper triangular matrix stored workl(ipntr(12)).
c
c\Authors
c Danny Sorensen Phuong Vu
c Richard Lehoucq CRPC / Rice University
c Chao Yang Houston, Texas
c Dept. of Computational &
c Applied Mathematics
c Rice University
c Houston, Texas
c
c\SCCS Information: @(#)
c FILE: neupd.F SID: 2.4 DATE OF SID: 7/31/96 RELEASE: 2
c
c\EndLib
c
c-----------------------------------------------------------------------
subroutine cneupd (rvec, howmny, select, d, z, ldz, sigma,
& workev, bmat, n, which, nev, tol,
& resid, ncv, v, ldv, iparam, ipntr, workd,
& workl, lworkl, rwork, info)
c
c %----------------------------------------------------%
c | Include files for debugging and timing information |
c %----------------------------------------------------%
c
include 'debug.h'
include 'stat.h'
c
c %------------------%
c | Scalar Arguments |
c %------------------%
c
character bmat, howmny, which*2
logical rvec
integer info, ldz, ldv, lworkl, n, ncv, nev
Complex
& sigma
Real
& tol
c
c %-----------------%
c | Array Arguments |
c %-----------------%
c
integer iparam(11), ipntr(14)
logical select(ncv)
Real
& rwork(ncv)
Complex
& d(nev), resid(n), v(ldv,ncv), z(ldz, nev),
& workd(3*n), workl(lworkl), workev(2*ncv)
c
c %------------%
c | Parameters |
c %------------%
c
Complex
& one, zero
parameter (one = (1.0E+0, 0.0E+0), zero = (0.0E+0, 0.0E+0))
c
c %---------------%
c | Local Scalars |
c %---------------%
c
character type*6
integer bounds, ierr, ih, ihbds, iheig, nconv,
& invsub, iuptri, iwev, j,
& ldh, ldq, mode, msglvl, ritz, wr, k,
& irz, ibd, ktrord, outncv, iq
Complex
& rnorm, temp, vl(1)
Real
& thres, conds, sep, rtemp, eps23
logical reord
c
c %----------------------%
c | External Subroutines |
c %----------------------%
c
external ccopy, cgeru, cgeqr2, clacpy, cmout,
& cunm2r, ctrmm, cvout, ivout,
& clahqr
c
c %--------------------%
c | External Functions |
c %--------------------%
c
Real
& scnrm2, slamch, slapy2
external scnrm2, slamch, slapy2
c
Complex
& cdotc
external cdotc
c
c %-----------------------%
c | Executable Statements |
c %-----------------------%
c
c %------------------------%
c | Set default parameters |
c %------------------------%
c
msglvl = mceupd
mode = iparam(7)
nconv = iparam(5)
info = 0
c
c
c %---------------------------------%
c | Get machine dependent constant. |
c %---------------------------------%
c
eps23 = slamch('Epsilon-Machine')
eps23 = eps23**(2.0E+0 / 3.0E+0)
c
c %-------------------------------%
c | Quick return |
c | Check for incompatible input |
c %-------------------------------%
c
ierr = 0
c
if (nconv .le. 0) then
ierr = -14
else if (n .le. 0) then
ierr = -1
else if (nev .le. 0) then
ierr = -2
else if (ncv .le. nev+1 .or. ncv .gt. n) then
ierr = -3
else if (which .ne. 'LM' .and.
& which .ne. 'SM' .and.
& which .ne. 'LR' .and.
& which .ne. 'SR' .and.
& which .ne. 'LI' .and.
& which .ne. 'SI') then
ierr = -5
else if (bmat .ne. 'I' .and. bmat .ne. 'G') then
ierr = -6
else if (lworkl .lt. 3*ncv**2 + 4*ncv) then
ierr = -7
else if ( (howmny .ne. 'A' .and.
& howmny .ne. 'P' .and.
& howmny .ne. 'S') .and. rvec ) then
ierr = -13
else if (howmny .eq. 'S' ) then
ierr = -12
end if
c
if (mode .eq. 1 .or. mode .eq. 2) then
type = 'REGULR'
else if (mode .eq. 3 ) then
type = 'SHIFTI'
else
ierr = -10
end if
if (mode .eq. 1 .and. bmat .eq. 'G') ierr = -11
c
c %------------%
c | Error Exit |
c %------------%
c
if (ierr .ne. 0) then
info = ierr
go to 9000
end if
c
c %--------------------------------------------------------%
c | Pointer into WORKL for address of H, RITZ, WORKEV, Q |
c | etc... and the remaining workspace. |
c | Also update pointer to be used on output. |
c | Memory is laid out as follows: |
c | workl(1:ncv*ncv) := generated Hessenberg matrix |
c | workl(ncv*ncv+1:ncv*ncv+ncv) := ritz values |
c | workl(ncv*ncv+ncv+1:ncv*ncv+2*ncv) := error bounds |
c %--------------------------------------------------------%
c
c %-----------------------------------------------------------%
c | The following is used and set by CNEUPD. |
c | workl(ncv*ncv+2*ncv+1:ncv*ncv+3*ncv) := The untransformed |
c | Ritz values. |
c | workl(ncv*ncv+3*ncv+1:ncv*ncv+4*ncv) := The untransformed |
c | error bounds of |
c | the Ritz values |
c | workl(ncv*ncv+4*ncv+1:2*ncv*ncv+4*ncv) := Holds the upper |
c | triangular matrix |
c | for H. |
c | workl(2*ncv*ncv+4*ncv+1: 3*ncv*ncv+4*ncv) := Holds the |
c | associated matrix |
c | representation of |
c | the invariant |
c | subspace for H. |
c | GRAND total of NCV * ( 3 * NCV + 4 ) locations. |
c %-----------------------------------------------------------%
c
ih = ipntr(5)
ritz = ipntr(6)
iq = ipntr(7)
bounds = ipntr(8)
ldh = ncv
ldq = ncv
iheig = bounds + ldh
ihbds = iheig + ldh
iuptri = ihbds + ldh
invsub = iuptri + ldh*ncv
ipntr(9) = iheig
ipntr(11) = ihbds
ipntr(12) = iuptri
ipntr(13) = invsub
wr = 1
iwev = wr + ncv
c
c %-----------------------------------------%
c | irz points to the Ritz values computed |
c | by _neigh before exiting _naup2. |
c | ibd points to the Ritz estimates |
c | computed by _neigh before exiting |
c | _naup2. |
c %-----------------------------------------%
c
irz = ipntr(14) + ncv*ncv
ibd = irz + ncv
c
c %------------------------------------%
c | RNORM is B-norm of the RESID(1:N). |
c %------------------------------------%
c
rnorm = workl(ih+2)
workl(ih+2) = zero
c
if (rvec) then
c
c %-------------------------------------------%
c | Get converged Ritz value on the boundary. |
c | Note: converged Ritz values have been |
c | placed in the first NCONV locations in |
c | workl(ritz). They have been sorted |
c | (in _naup2) according to the WHICH |
c | selection criterion |
c %-------------------------------------------%
c
if (which .eq. 'LM' .or. which .eq. 'SM') then
thres = slapy2(real(workl(ritz)),aimag(workl(ritz)))
else if (which .eq. 'LR' .or. which .eq. 'SR') then
thres = real(workl(ritz))
else if (which .eq. 'LI' .or. which .eq. 'SI') then
thres = aimag(workl(ritz))
end if
if (msglvl .gt. 2) then
call svout(logfil, 1, [thres], ndigit,
& '_neupd: Threshold eigenvalue used for re-ordering')
end if
c
c %---------------------------------------------------------%
c | Check to see if all converged Ritz values appear at the |
c | at the top of the upper triangular matrix computed by |
c | _neigh in _naup2. This is done in the following way: |
c | |
c | 1) For each Ritz value from _neigh, compare it with the |
c | threshold Ritz value computed above to determine |
c | whether it is a wanted one. |
c | |
c | 2) If it is wanted, then check the corresponding Ritz |
c | estimate to see if it has converged. If it has, set |
c | correponding entry in the logical array SELECT to |
c | .TRUE.. |
c | |
c | If SELECT(j) = .TRUE. and j > NCONV, then there is a |
c | converged Ritz value that does not appear at the top of |
c | the upper triangular matrix computed by _neigh in |
c | _naup2. Reordering is needed. |
c %---------------------------------------------------------%
c
reord = .false.
ktrord = 0
do 10 j = 0, ncv-1
select(j+1) = .false.
if (which .eq. 'LM') then
if ( slapy2(real(workl(irz+j)),
& aimag(workl(irz+j))) .ge. thres ) then
rtemp = max( eps23, slapy2(real(workl(irz+j-1)),
& aimag(workl(irz+j-1))) )
if ( slapy2(real(workl(ibd+j)),
& aimag(workl(ibd+j))) .le. tol*rtemp )
& select(j+1) = .true.
end if
else if (which .eq. 'SM') then
if ( slapy2(real(workl(irz+j)),
& aimag(workl(irz+j))) .le. thres ) then
rtemp = max( eps23, slapy2(real(workl(irz+j-1)),
& aimag(workl(irz+j-1))) )
if ( slapy2(real(workl(ibd+j)),
& aimag(workl(ibd+j))) .le. tol*rtemp )
& select(j+1) = .true.
end if
else if (which .eq. 'LR') then
if ( real(workl(irz+j)) .ge. thres ) then
rtemp = max( eps23, slapy2(real(workl(irz+j-1)),
& aimag(workl(irz+j-1))) )
if ( slapy2(real(workl(ibd+j)),
& aimag(workl(ibd+j))) .le. tol*rtemp )
& select(j+1) = .true.
end if
else if (which .eq. 'SR') then
if ( real(workl(irz+j)) .le. thres ) then
rtemp = max( eps23, slapy2(real(workl(irz+j-1)),
& aimag(workl(irz+j-1))) )
if ( slapy2(real(workl(ibd+j)),
& aimag(workl(ibd+j))) .le. tol*rtemp )
& select(j+1) = .true.
end if
else if (which .eq. 'LI') then
if ( aimag(workl(irz+j)) .ge. thres ) then
rtemp = max( eps23, slapy2(real(workl(irz+j-1)),
& aimag(workl(irz+j-1))) )
if ( slapy2(real(workl(ibd+j)),
& aimag(workl(ibd+j))) .le. tol*rtemp )
& select(j+1) = .true.
end if
else if (which .eq. 'SI') then
if ( aimag(workl(irz+j)) .le. thres ) then
rtemp = max( eps23, slapy2(real(workl(irz+j-1)),
& aimag(workl(irz+j-1))) )
if ( slapy2(real(workl(ibd+j)),
& aimag(workl(ibd+j))) .le. tol*rtemp )
& select(j+1) = .true.
end if
end if
if (j+1 .gt. nconv ) reord = ( select(j+1) .or. reord )
if (select(j+1)) ktrord = ktrord + 1
10 continue
c
if (msglvl .gt. 2) then
call ivout(logfil, 1, [ktrord], ndigit,
& '_neupd: Number of specified eigenvalues')
call ivout(logfil, 1, [nconv], ndigit,
& '_neupd: Number of "converged" eigenvalues')
end if
c
c if (ktrord .gt. nconv) then
c
c %-----------------------------------%
c | More than NCONV Ritz values have |
c | "converged", and they all satisfy |
c | the WHICH selection criterion. |
c %-----------------------------------%
c
c iparam(6) = ktrord
c
c end if
c
c %-------------------------------------------------------%
c | Call LAPACK routine clahqr to compute the Schur form |
c | of the upper Hessenberg matrix returned by CNAUPD. |
c | Make a copy of the upper Hessenberg matrix. |
c | Initialize the Schur vector matrix Q to the identity. |
c %-------------------------------------------------------%
c
call ccopy (ldh*ncv, workl(ih), 1, workl(iuptri), 1)
call claset ('All', ncv, ncv, zero, one, workl(invsub), ldq)
call clahqr (.true., .true., ncv, 1, ncv, workl(iuptri),
& ldh, workl(iheig), 1, ncv, workl(invsub), ldq, ierr)
call ccopy (ncv, workl(invsub+ncv-1), ldq, workl(ihbds), 1)
c
if (ierr .ne. 0) then
info = -8
go to 9000
end if
c
if (msglvl .gt. 1) then
call cvout (logfil, ncv, workl(iheig), ndigit,
& '_neupd: Eigenvalues of H')
call cvout (logfil, ncv, workl(ihbds), ndigit,
& '_neupd: Last row of the Schur vector matrix')
if (msglvl .gt. 3) then
call cmout (logfil, ncv, ncv, workl(iuptri), ldh, ndigit,
& '_neupd: The upper triangular matrix ')
end if
end if
c
if (reord) then
c
c %-----------------------------------------------%
c | Reorder the computed upper triangular matrix. |
c %-----------------------------------------------%
c
call ctrsen ('None', 'V', select, ncv, workl(iuptri), ldh,
& workl(invsub), ldq, workl(iheig), nconv, conds, sep,
& workev, ncv, ierr)
c
if (ierr .eq. 1) then
info = 1
go to 9000
end if
c
if (msglvl .gt. 2) then
call cvout (logfil, ncv, workl(iheig), ndigit,
& '_neupd: Eigenvalues of H--reordered')
if (msglvl .gt. 3) then
call cmout (logfil, ncv, ncv, workl(iuptri), ldq,
& ndigit,
& '_neupd: Triangular matrix after re-ordering')
end if
end if
c
end if
c
c %---------------------------------------------%
c | Copy the last row of the Schur basis matrix |
c | to workl(ihbds). This vector will be used |
c | to compute the Ritz estimates of converged |
c | Ritz values. |
c %---------------------------------------------%
c
call ccopy (ncv, workl(invsub+ncv-1), ldq, workl(ihbds), 1)
c
c %--------------------------------------------%
c | Place the computed eigenvalues of H into D |
c | if a spectral transformation was not used. |
c %--------------------------------------------%
c
if (type .eq. 'REGULR') then
call ccopy (nconv, workl(iheig), 1, d, 1)
end if
c
c %----------------------------------------------------------%
c | Compute the QR factorization of the matrix representing |
c | the wanted invariant subspace located in the first NCONV |
c | columns of workl(invsub,ldq). |
c %----------------------------------------------------------%
c
call cgeqr2 (ncv, nconv, workl(invsub), ldq, workev,
& workev(ncv+1), ierr)
c
c %--------------------------------------------------------%
c | * Postmultiply V by Q using cunm2r. |
c | * Copy the first NCONV columns of VQ into Z. |
c | * Postmultiply Z by R. |
c | The N by NCONV matrix Z is now a matrix representation |
c | of the approximate invariant subspace associated with |
c | the Ritz values in workl(iheig). The first NCONV |
c | columns of V are now approximate Schur vectors |
c | associated with the upper triangular matrix of order |
c | NCONV in workl(iuptri). |
c %--------------------------------------------------------%
c
call cunm2r ('Right', 'Notranspose', n, ncv, nconv,
& workl(invsub), ldq, workev, v, ldv, workd(n+1),
& ierr)
call clacpy ('All', n, nconv, v, ldv, z, ldz)
c
do 20 j=1, nconv
c
c %---------------------------------------------------%
c | Perform both a column and row scaling if the |
c | diagonal element of workl(invsub,ldq) is negative |
c | I'm lazy and don't take advantage of the upper |
c | triangular form of workl(iuptri,ldq). |
c | Note that since Q is orthogonal, R is a diagonal |
c | matrix consisting of plus or minus ones. |
c %---------------------------------------------------%
c
if ( real( workl(invsub+(j-1)*ldq+j-1) ) .lt.
& real(zero) ) then
call cscal (nconv, -one, workl(iuptri+j-1), ldq)
call cscal (nconv, -one, workl(iuptri+(j-1)*ldq), 1)
end if
c
20 continue
c
if (howmny .eq. 'A') then
c
c %--------------------------------------------%
c | Compute the NCONV wanted eigenvectors of T |
c | located in workl(iuptri,ldq). |
c %--------------------------------------------%
c
do 30 j=1, ncv
if (j .le. nconv) then
select(j) = .true.
else
select(j) = .false.
end if
30 continue
c
call ctrevc ('Right', 'Select', select, ncv, workl(iuptri),
& ldq, vl, 1, workl(invsub), ldq, ncv, outncv, workev,
& rwork, ierr)
c
if (ierr .ne. 0) then
info = -9
go to 9000
end if
c
c %------------------------------------------------%
c | Scale the returning eigenvectors so that their |
c | Euclidean norms are all one. LAPACK subroutine |
c | ctrevc returns each eigenvector normalized so |
c | that the element of largest magnitude has |
c | magnitude 1. |
c %------------------------------------------------%
c
do 40 j=1, nconv
rtemp = scnrm2(ncv, workl(invsub+(j-1)*ldq), 1)
rtemp = real(one) / rtemp
call csscal ( ncv, rtemp,
& workl(invsub+(j-1)*ldq), 1 )
c
c %------------------------------------------%
c | Ritz estimates can be obtained by taking |
c | the inner product of the last row of the |
c | Schur basis of H with eigenvectors of T. |
c | Note that the eigenvector matrix of T is |
c | upper triangular, thus the length of the |
c | inner product can be set to j. |
c %------------------------------------------%
c
workev(j) = cdotc(j, workl(ihbds), 1,
& workl(invsub+(j-1)*ldq), 1)
40 continue
c
if (msglvl .gt. 2) then
call ccopy(nconv, workl(invsub+ncv-1), ldq,
& workl(ihbds), 1)
call cvout (logfil, nconv, workl(ihbds), ndigit,
& '_neupd: Last row of the eigenvector matrix for T')
if (msglvl .gt. 3) then
call cmout (logfil, ncv, ncv, workl(invsub), ldq,
& ndigit, '_neupd: The eigenvector matrix for T')
end if
end if
c
c %---------------------------------------%
c | Copy Ritz estimates into workl(ihbds) |
c %---------------------------------------%
c
call ccopy(nconv, workev, 1, workl(ihbds), 1)
c
c %----------------------------------------------%
c | The eigenvector matrix Q of T is triangular. |
c | Form Z*Q. |
c %----------------------------------------------%
c
call ctrmm ('Right', 'Upper', 'No transpose', 'Non-unit',
& n, nconv, one, workl(invsub), ldq, z, ldz)
c
end if
c
else
c
c %--------------------------------------------------%
c | An approximate invariant subspace is not needed. |
c | Place the Ritz values computed CNAUPD into D. |
c %--------------------------------------------------%
c
call ccopy (nconv, workl(ritz), 1, d, 1)
call ccopy (nconv, workl(ritz), 1, workl(iheig), 1)
call ccopy (nconv, workl(bounds), 1, workl(ihbds), 1)
c
end if
c
c %------------------------------------------------%
c | Transform the Ritz values and possibly vectors |
c | and corresponding error bounds of OP to those |
c | of A*x = lambda*B*x. |
c %------------------------------------------------%
c
if (type .eq. 'REGULR') then
c
if (rvec)
& call cscal(ncv, rnorm, workl(ihbds), 1)
c
else
c
c %---------------------------------------%
c | A spectral transformation was used. |
c | * Determine the Ritz estimates of the |
c | Ritz values in the original system. |
c %---------------------------------------%
c
if (rvec)
& call cscal(ncv, rnorm, workl(ihbds), 1)
c
do 50 k=1, ncv
temp = workl(iheig+k-1)
workl(ihbds+k-1) = workl(ihbds+k-1) / temp / temp
50 continue
c
end if
c
c %-----------------------------------------------------------%
c | * Transform the Ritz values back to the original system. |
c | For TYPE = 'SHIFTI' the transformation is |
c | lambda = 1/theta + sigma |
c | NOTES: |
c | *The Ritz vectors are not affected by the transformation. |
c %-----------------------------------------------------------%
c
if (type .eq. 'SHIFTI') then
do 60 k=1, nconv
d(k) = one / workl(iheig+k-1) + sigma
60 continue
end if
c
if (type .ne. 'REGULR' .and. msglvl .gt. 1) then
call cvout (logfil, nconv, d, ndigit,
& '_neupd: Untransformed Ritz values.')
call cvout (logfil, nconv, workl(ihbds), ndigit,
& '_neupd: Ritz estimates of the untransformed Ritz values.')
else if ( msglvl .gt. 1) then
call cvout (logfil, nconv, d, ndigit,
& '_neupd: Converged Ritz values.')
call cvout (logfil, nconv, workl(ihbds), ndigit,
& '_neupd: Associated Ritz estimates.')
end if
c
c %-------------------------------------------------%
c | Eigenvector Purification step. Formally perform |
c | one of inverse subspace iteration. Only used |
c | for MODE = 3. See reference 3. |
c %-------------------------------------------------%
c
if (rvec .and. howmny .eq. 'A' .and. type .eq. 'SHIFTI') then
c
c %------------------------------------------------%
c | Purify the computed Ritz vectors by adding a |
c | little bit of the residual vector: |
c | T |
c | resid(:)*( e s ) / theta |
c | NCV |
c | where H s = s theta. |
c %------------------------------------------------%
c
do 100 j=1, nconv
if (workl(iheig+j-1) .ne. zero) then
workev(j) = workl(invsub+(j-1)*ldq+ncv-1) /
& workl(iheig+j-1)
endif
100 continue
c %---------------------------------------%
c | Perform a rank one update to Z and |
c | purify all the Ritz vectors together. |
c %---------------------------------------%
c
call cgeru (n, nconv, one, resid, 1, workev, 1, z, ldz)
c
end if
c
9000 continue
c
return
c
c %---------------%
c | End of cneupd |
c %---------------%
c
end
