Subversion Repositories OpenSees

Rev

Blame | Last modification | View Log | RSS feed

c\BeginDoc
c
c\Name: zneupd
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 ZNAUPD.  ZNAUPD 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 ZNAUPD.  For a brief
c  description, see definitions of IPARAM(7), MODE and WHICH in the
c  documentation of ZNAUPD.
c
c\Usage:
c  call zneupd
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*16 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*16 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 ZNAUPD.  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*16  (INPUT)
c          If IPARAM(7) = 3 then SIGMA represents the shift.
c          Not referenced if IPARAM(7) = 1 or 2.
c
c  WORKEV  Complex*16 work array of dimension 2*NCV.  (WORKSPACE)
c
c  **** The remaining arguments MUST be the same as for the   ****
c  **** call to ZNAUPD 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 ZNEUPD following the last call
c         to ZNAUPD.  These arguments MUST NOT BE MODIFIED between
c         the the last call to ZNAUPD and the call to ZNEUPD.
c
c  Three of these parameters (V, WORKL and INFO) are also output parameters:
c
c  V       Complex*16 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 ZNAUPD .
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   Double precision work array of length LWORKL.  (OUTPUT/WORKSPACE)
c          WORKL(1:ncv*ncv+2*ncv) contains information obtained in
c          znaupd.  They are not changed by zneupd.
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 zneupd.
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                     zneupd 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 ztrsen.
c                Re-enter subroutine zneupd 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 ztrevc.
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: ZNAUPD 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     zmout   ARPACK utility routine that prints matrices
c     zvout   ARPACK utility routine that prints vectors.
c     zgeqr2  LAPACK routine that computes the QR factorization of
c             a matrix.
c     zlacpy  LAPACK matrix copy routine.
c     zlahqr  LAPACK routine that computes the Schur form of a
c             upper Hessenberg matrix.
c     zlaset  LAPACK matrix initialization routine.
c     ztrevc  LAPACK routine to compute the eigenvectors of a matrix
c             in upper triangular form.
c     ztrsen  LAPACK routine that re-orders the Schur form.
c     zunm2r  LAPACK routine that applies an orthogonal matrix in
c             factored form.
c     dlamch  LAPACK routine that determines machine constants.
c     ztrmm   Level 3 BLAS matrix times an upper triangular matrix.
c     zgeru   Level 2 BLAS rank one update to a matrix.
c     zcopy   Level 1 BLAS that copies one vector to another .
c     zscal   Level 1 BLAS that scales a vector.
c     zdscal  Level 1 BLAS that scales a complex vector by a real number.
c     dznrm2  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.5   DATE OF SID: 11/05/99   RELEASE: 2
c
c\EndLib
c
c-----------------------------------------------------------------------
      subroutine zneupd (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*16    
     &           sigma
      Double precision
     &           tol
c
c     %-----------------%
c     | Array Arguments |
c     %-----------------%
c
      integer    iparam(11), ipntr(14)
      logical    select(ncv)
      Double precision
     &           rwork(ncv)
      Complex*16
     &           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*16
     &           one, zero
      parameter  (one = (1.0D+0, 0.0D+0), zero = (0.0D+0, 0.0D+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*16
     &           rnorm, temp, vl(1)
      Double precision
     &           thres, conds, sep, rtemp, eps23
      logical    reord
c
c     %----------------------%
c     | External Subroutines |
c     %----------------------%
c
      external   zcopy, zgeru, zgeqr2, zlacpy, zmout,
     &           zunm2r, ztrmm, zvout, ivout,
     &           zlahqr
c  
c     %--------------------%
c     | External Functions |
c     %--------------------%
c
      Double precision
     &           dznrm2, dlamch, dlapy2
      external   dznrm2, dlamch, dlapy2
c
      Complex*16
     &           zdotc
      external   zdotc
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 = dlamch('Epsilon-Machine')
      eps23 = eps23**(2.0D+0 / 3.0D+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 ZNEUPD.                  |
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 = dlapy2(dble(workl(ritz)),dimag(workl(ritz)))
         else if (which .eq. 'LR' .or. which .eq. 'SR') then
            thres = dble(workl(ritz))
         else if (which .eq. 'LI' .or. which .eq. 'SI') then
            thres = dimag(workl(ritz))
         end if
         if (msglvl .gt. 2) then
            call dvout(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 ( dlapy2(dble(workl(irz+j)),
     &                      dimag(workl(irz+j))) .ge. thres ) then
                  rtemp = max( eps23, dlapy2(dble(workl(irz+j)),
     &                         dimag(workl(irz+j))) )
                  if ( dlapy2(dble(workl(ibd+j)),
     &                 dimag(workl(ibd+j))) .le. tol*rtemp )
     &               select(j+1) = .true.
               end if
            else if (which .eq. 'SM') then
               if ( dlapy2(dble(workl(irz+j)),
     &                      dimag(workl(irz+j))) .le. thres )  then
                  rtemp = max( eps23, dlapy2(dble(workl(irz+j)),
     &                         dimag(workl(irz+j))) )
                  if ( dlapy2(dble(workl(ibd+j)),
     &                 dimag(workl(ibd+j))) .le. tol*rtemp )
     &               select(j+1) = .true.
               end if
            else if (which .eq. 'LR') then
               if ( dble(workl(irz+j)) .ge. thres ) then
                  rtemp = max( eps23, dlapy2(dble(workl(irz+j)),
     &                         dimag(workl(irz+j))) )
                  if ( dlapy2(dble(workl(ibd+j)),
     &                 dimag(workl(ibd+j))) .le. tol*rtemp )
     &               select(j+1) = .true.
               end if
            else if (which .eq. 'SR') then
               if ( dble(workl(irz+j)) .le. thres ) then
                  rtemp = max( eps23, dlapy2(dble(workl(irz+j)),
     &                    dimag(workl(irz+j))) )
                  if ( dlapy2(dble(workl(ibd+j)),
     &                 dimag(workl(ibd+j))) .le. tol*rtemp )
     &               select(j+1) = .true.
               end if
            else if (which .eq. 'LI') then
               if ( dimag(workl(irz+j)) .ge. thres ) then
                  rtemp = max( eps23, dlapy2(dble(workl(irz+j)),
     &                    dimag(workl(irz+j))) )
                  if ( dlapy2(dble(workl(ibd+j)),
     &                 dimag(workl(ibd+j))) .le. tol*rtemp )
     &               select(j+1) = .true.
               end if
            else if (which .eq. 'SI') then
               if ( dimag(workl(irz+j)) .le. thres ) then
                  rtemp = max( eps23, dlapy2(dble(workl(irz+j)),
     &                    dimag(workl(irz+j))) )
                  if ( dlapy2(dble(workl(ibd+j)),
     &                 dimag(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 zlahqr to compute the Schur form  |
c        | of the upper Hessenberg matrix returned by ZNAUPD.    |
c        | Make a copy of the upper Hessenberg matrix.           |
c        | Initialize the Schur vector matrix Q to the identity. |
c        %-------------------------------------------------------%
c
         call zcopy (ldh*ncv, workl(ih), 1, workl(iuptri), 1)
         call zlaset ('All', ncv, ncv, zero, one, workl(invsub), ldq)
         call zlahqr (.true., .true., ncv, 1, ncv, workl(iuptri),
     &        ldh, workl(iheig), 1, ncv, workl(invsub), ldq, ierr)
         call zcopy (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 zvout (logfil, ncv, workl(iheig), ndigit,
     &           '_neupd: Eigenvalues of H')
            call zvout (logfil, ncv, workl(ihbds), ndigit,
     &           '_neupd: Last row of the Schur vector matrix')
            if (msglvl .gt. 3) then
               call zmout (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 ztrsen ('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 zvout (logfil, ncv, workl(iheig), ndigit,
     &           '_neupd: Eigenvalues of H--reordered')
                if (msglvl .gt. 3) then
                   call zmout (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 zcopy (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 zcopy (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 zgeqr2 (ncv, nconv, workl(invsub), ldq, workev,
     &                workev(ncv+1), ierr)
c
c        %--------------------------------------------------------%
c        | * Postmultiply V by Q using zunm2r.                    |
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 zunm2r ('Right', 'Notranspose', n, ncv, nconv,
     &        workl(invsub), ldq, workev, v, ldv, workd(n+1),
     &        ierr)
         call zlacpy ('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 ( dble( workl(invsub+(j-1)*ldq+j-1) ) .lt.
     &                  dble(zero) ) then
               call zscal (nconv, -one, workl(iuptri+j-1), ldq)
               call zscal (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 ztrevc ('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           | ztrevc returns each eigenvector normalized so  |
c           | that the element of largest magnitude has      |
c           | magnitude 1.                                   |
c           %------------------------------------------------%
c
            do 40 j=1, nconv
                  rtemp = dznrm2(ncv, workl(invsub+(j-1)*ldq), 1)
                  rtemp = dble(one) / rtemp
                  call zdscal ( 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) = zdotc(j, workl(ihbds), 1,
     &                        workl(invsub+(j-1)*ldq), 1)
 40         continue
c
            if (msglvl .gt. 2) then
               call zcopy(nconv, workl(invsub+ncv-1), ldq,
     &                    workl(ihbds), 1)
               call zvout (logfil, nconv, workl(ihbds), ndigit,
     &              '_neupd: Last row of the eigenvector matrix for T')
               if (msglvl .gt. 3) then
                  call zmout (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 zcopy(nconv, workev, 1, workl(ihbds), 1)
c
c           %----------------------------------------------%
c           | The eigenvector matrix Q of T is triangular. |
c           | Form Z*Q.                                    |
c           %----------------------------------------------%
c
            call ztrmm ('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 ZNAUPD into D.    |
c        %--------------------------------------------------%
c
         call zcopy (nconv, workl(ritz), 1, d, 1)
         call zcopy (nconv, workl(ritz), 1, workl(iheig), 1)
         call zcopy (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 zscal(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 zscal(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 zvout (logfil, nconv, d, ndigit,
     &     '_neupd: Untransformed Ritz values.')
         call zvout (logfil, nconv, workl(ihbds), ndigit,
     &     '_neupd: Ritz estimates of the untransformed Ritz values.')
      else if ( msglvl .gt. 1) then
         call zvout (logfil, nconv, d, ndigit,
     &     '_neupd: Converged Ritz values.')
         call zvout (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 zgeru (n, nconv, one, resid, 1, workev, 1, z, ldz)
c
      end if
c
 9000 continue
c
      return
c    
c     %---------------%
c     | End of zneupd |
c     %---------------%
c
      end