docs/anasazi/mypdsaupd_8f_source.html

c @HEADER

c *****************************************************************************

c                 Anasazi: Block Eigensolvers Package

c

c Copyright 2004 NTESS and the Anasazi contributors.

c SPDX-License-Identifier: BSD-3-Clause

c *****************************************************************************

c @HEADER


c This software is a result of the research described in the report

c

c     "A comparison of algorithms for modal analysis in the absence

c     of a sparse direct method", P. Arbenz, R. Lehoucq, and U. Hetmaniuk,

c     Sandia National Laboratories, Technical report SAND2003-1028J.

c

c It is based on the Epetra, AztecOO, and ML packages defined in the Trilinos

c framework ( http://trilinos.org/ ).

c-----------------------------------------------------------------------

c\BeginDoc

c

c\Name: mypdsaupd

c

c Message Passing Layer: MPI

c

c\Description:

c

c  Reverse communication interface for the Implicitly Restarted Arnoldi

c  Iteration.  For symmetric problems this reduces to a variant of the Lanczos

c  method.  This method has been designed to compute approximations to a

c  few eigenpairs of a linear operator OP that is real and symmetric

c  with respect to a real positive semi-definite symmetric matrix B,

c  i.e.

c

c       B*OP = (OP`)*B.

c

c  Another way to express this condition is

c

c       < x,OPy > = < OPx,y >  where < z,w > = z`Bw  .

c

c  In the standard eigenproblem B is the identity matrix.

c  ( A` denotes transpose of A)

c

c  The computed approximate eigenvalues are called Ritz values and

c  the corresponding approximate eigenvectors are called Ritz vectors.

c

c  mypdsaupd  is usually called iteratively to solve one of the

c  following problems:

c

c  Mode 1:  A*x = lambda*x, A symmetric

c           ===> OP = A  and  B = I.

c

c  Mode 2:  A*x = lambda*M*x, A symmetric, M symmetric positive definite

c           ===> OP = inv[M]*A  and  B = M.

c           ===> (If M can be factored see remark 3 below)

c

c  Mode 3:  K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite

c           ===> OP = (inv[K - sigma*M])*M  and  B = M.

c           ===> Shift-and-Invert mode

c

c  Mode 4:  K*x = lambda*KG*x, K symmetric positive semi-definite,

c           KG symmetric indefinite

c           ===> OP = (inv[K - sigma*KG])*K  and  B = K.

c           ===> Buckling mode

c

c  Mode 5:  A*x = lambda*M*x, A symmetric, M symmetric positive semi-definite

c           ===> OP = inv[A - sigma*M]*[A + sigma*M]  and  B = M.

c           ===> Cayley transformed mode

c

c  NOTE: The action of w <- inv[A - sigma*M]*v or w <- inv[M]*v

c        should be accomplished either by a direct method

c        using a sparse matrix factorization and solving

c

c           [A - sigma*M]*w = v  or M*w = v,

c

c        or through an iterative method for solving these

c        systems.  If an iterative method is used, the

c        convergence test must be more stringent than

c        the accuracy requirements for the eigenvalue

c        approximations.

c

c\Usage:

c  call mypdsaupd

c     ( COMM, IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM,

c       IPNTR, WORKD, WORKL, LWORKL, INFO, VERBOSE )

c

c\Arguments

c  COMM    MPI  Communicator for the processor grid.  (INPUT)

c

c  IDO     Integer.  (INPUT/OUTPUT)

c          Reverse communication flag.  IDO must be zero on the first

c          call to pdsaupd .  IDO will be set internally to

c          indicate the type of operation to be performed.  Control is

c          then given back to the calling routine which has the

c          responsibility to carry out the requested operation and call

c          pdsaupd  with the result.  The operand is given in

c          WORKD(IPNTR(1)), the result must be put in WORKD(IPNTR(2)).

c          (If Mode = 2 see remark 5 below)

c          -------------------------------------------------------------

c          IDO =  0: first call to the reverse communication interface

c          IDO = -1: compute  Y = OP * X  where

c                    IPNTR(1) is the pointer into WORKD for X,

c                    IPNTR(2) is the pointer into WORKD for Y.

c                    This is for the initialization phase to force the

c                    starting vector into the range of OP.

c          IDO =  1: compute  Y = OP * X where

c                    IPNTR(1) is the pointer into WORKD for X,

c                    IPNTR(2) is the pointer into WORKD for Y.

c                    In mode 3,4 and 5, the vector B * X is already

c                    available in WORKD(ipntr(3)).  It does not

c                    need to be recomputed in forming OP * X.

c          IDO =  2: compute  Y = B * X  where

c                    IPNTR(1) is the pointer into WORKD for X,

c                    IPNTR(2) is the pointer into WORKD for Y.

c          IDO =  3: compute the IPARAM(8) shifts where

c                    IPNTR(11) is the pointer into WORKL for

c                    placing the shifts. See remark 6 below.

c    IDO =  4: user checks convergence. See remark 7 below.

c          IDO = 99: done

c          -------------------------------------------------------------

c

c  BMAT    Character*1.  (INPUT)

c          BMAT specifies the type of the matrix B that defines the

c          semi-inner product for the operator OP.

c          B = 'I' -> standard eigenvalue problem A*x = lambda*x

c          B = 'G' -> generalized eigenvalue problem A*x = lambda*B*x

c

c  N       Integer.  (INPUT)

c          Dimension of the eigenproblem.

c

c  WHICH   Character*2.  (INPUT)

c          Specify which of the Ritz values of OP to compute.

c

c          'LA' - compute the NEV largest (algebraic) eigenvalues.

c          'SA' - compute the NEV smallest (algebraic) eigenvalues.

c          'LM' - compute the NEV largest (in magnitude) eigenvalues.

c          'SM' - compute the NEV smallest (in magnitude) eigenvalues.

c          'BE' - compute NEV eigenvalues, half from each end of the

c                 spectrum.  When NEV is odd, compute one more from the

c                 high end than from the low end.

c           (see remark 1 below)

c

c  NEV     Integer.  (INPUT)

c          Number of eigenvalues of OP to be computed. 0 < NEV < N.

c

c  TOL     Double precision  scalar.  (INPUT)

c          Stopping criterion: the relative accuracy of the Ritz value

c          is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I)).

c          If TOL .LE. 0. is passed a default is set:

c          DEFAULT = DLAMCH ('EPS')  (machine precision as computed

c                    by the LAPACK auxiliary subroutine DLAMCH ).

c          Not referenced if IPARAM(4) is not equal to zero. See remark

c          7 below.

c

c

c  RESID   Double precision  array of length N.  (INPUT/OUTPUT)

c          On INPUT:

c          If INFO .EQ. 0, a random initial residual vector is used.

c          If INFO .NE. 0, RESID contains the initial residual vector,

c                          possibly from a previous run.

c          On OUTPUT:

c          RESID contains the final residual vector.

c

c  NCV     Integer.  (INPUT)

c          Number of columns of the matrix V (less than or equal to N).

c          This will indicate how many Lanczos vectors are generated

c          at each iteration.  After the startup phase in which NEV

c          Lanczos vectors are generated, the algorithm generates

c          NCV-NEV Lanczos vectors at each subsequent update iteration.

c          Most of the cost in generating each Lanczos vector is in the

c          matrix-vector product OP*x. (See remark 4 below).

c

c  V       Double precision  N by NCV array.  (OUTPUT)

c          The NCV columns of V contain the Lanczos basis vectors.

c

c  LDV     Integer.  (INPUT)

c          Leading dimension of V exactly as declared in the calling

c          program.

c

c  IPARAM  Integer array of length 11.  (INPUT/OUTPUT)

c          IPARAM(1) = ISHIFT: method for selecting the implicit shifts.

c          The shifts selected at each iteration are used to restart

c          the Arnoldi iteration in an implicit fashion.

c          -------------------------------------------------------------

c          ISHIFT = 0: the shifts are provided by the user via

c                      reverse communication.  The NCV eigenvalues of

c                      the current tridiagonal matrix T are returned in

c                      the part of WORKL array corresponding to RITZ.

c                      See remark 6 below.

c          ISHIFT = 1: exact shifts with respect to the reduced

c                      tridiagonal matrix T.  This is equivalent to

c                      restarting the iteration with a starting vector

c                      that is a linear combination of Ritz vectors

c                      associated with the "wanted" Ritz values.

c          -------------------------------------------------------------

c

c          IPARAM(2) = LEVEC

c          No longer referenced. See remark 2 below.

c

c          IPARAM(3) = MXITER

c          On INPUT:  maximum number of Arnoldi update iterations allowed.

c          On OUTPUT: actual number of Arnoldi update iterations taken.

c

c          IPARAM(4) = USRCHK: a non-zero value denotes that the user

c    will check convergence. See remark 7 below.

c

c          IPARAM(5) = NCONV: number of "converged" Ritz values.

c          If USRCHK is equal to zero, then, upon exit (IDO=99)

c          NCONV is the number of Ritz values that satisfy

c          the convergence criterion as determined by pdsaupd.f

c          If USRCHK is not equal to zero, then the user must set

c    NCONV to be the number of Ritz values that satisfy the

c          convergence criterion.

c

c          IPARAM(6) = IUPD

c          No longer referenced. Implicit restarting is ALWAYS used.

c

c          IPARAM(7) = MODE

c          On INPUT determines what type of eigenproblem is being solved.

c          Must be 1,2,3,4,5; See under \Description of pdsaupd  for the

c          five modes available.

c

c          IPARAM(8) = NP

c          When ido = 3 and the user provides shifts through reverse

c          communication (IPARAM(1)=0), pdsaupd  returns NP, the number

c          of shifts the user is to provide. 0 < NP <=NCV-NEV. See Remark

c          6 below.

c

c          IPARAM(9) = NUMOP, IPARAM(10) = NUMOPB, IPARAM(11) = NUMREO,

c          OUTPUT: NUMOP  = total number of OP*x operations,

c                  NUMOPB = total number of B*x operations if BMAT='G',

c                  NUMREO = total number of steps of re-orthogonalization.

c

c  IPNTR   Integer array of length 11.  (OUTPUT)

c          Pointer to mark the starting locations in the WORKD and WORKL

c          arrays for matrices/vectors used by the Lanczos iteration.

c          -------------------------------------------------------------

c          IPNTR(1): pointer to the current operand vector X in WORKD.

c          IPNTR(2): pointer to the current result vector Y in WORKD.

c          IPNTR(3): pointer to the vector B * X in WORKD when used in

c                    the shift-and-invert mode.

c          IPNTR(4): pointer to the next available location in WORKL

c                    that is untouched by the program.

c          IPNTR(5): pointer to the NCV by 2 tridiagonal matrix T in WORKL.

c          IPNTR(6): pointer to the NCV RITZ values array in WORKL.

c          IPNTR(7): pointer to the Ritz estimates in array WORKL associated

c                    with the Ritz values located in RITZ in WORKL.

c          IPNTR(11): pointer to the NP shifts in WORKL. See Remark 6 below.

c

c          Note: IPNTR(8:10) is only referenced by pdseupd . See Remark 2.

c          IPNTR(8): pointer to the NCV RITZ values of the original system.

c          IPNTR(9): pointer to the NCV corresponding error bounds.

c          IPNTR(10): pointer to the NCV by NCV matrix of eigenvectors

c                     of the tridiagonal matrix T. Only referenced by

c                     pdseupd  if RVEC = .TRUE. See Remarks.

c          -------------------------------------------------------------

c

c  WORKD   Double precision  work array of length 3*N.  (REVERSE COMMUNICATION)

c          Distributed array to be used in the basic Arnoldi iteration

c          for reverse communication.  The user should not use WORKD

c          as temporary workspace during the iteration. Upon termination

c          WORKD(1:N) contains B*RESID(1:N). If the Ritz vectors are desired

c          subroutine pdseupd  uses this output.

c          See Data Distribution Note below.

c

c  WORKL   Double precision  work array of length LWORKL.  (OUTPUT/WORKSPACE)

c          Private (replicated) array on each PE or array allocated on

c          the front end.  See Data Distribution Note below.

c

c  LWORKL  Integer.  (INPUT)

c          LWORKL must be at least NCV**2 + 8*NCV .

c

c  INFO    Integer.  (INPUT/OUTPUT)

c          If INFO .EQ. 0, a randomly initial residual vector is used.

c          If INFO .NE. 0, RESID contains the initial residual vector,

c                          possibly from a previous run.

c          Error flag on output.

c          =  0: Normal exit.

c          =  1: Maximum number of iterations taken.

c                All possible eigenvalues of OP has been found. IPARAM(5)

c                returns the number of wanted converged Ritz values.

c          =  2: No longer an informational error. Deprecated starting

c                with release 2 of ARPACK.

c          =  3: No shifts could be applied during a cycle of the

c                Implicitly restarted Arnoldi iteration. One possibility

c                is to increase the size of NCV relative to NEV.

c                See remark 4 below.

c          = -1: N must be positive.

c          = -2: NEV must be positive.

c          = -3: NCV must be greater than NEV and less than or equal to N.

c          = -4: The maximum number of Arnoldi update iterations allowed

c                must be greater than zero.

c          = -5: WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'.

c          = -6: BMAT must be one of 'I' or 'G'.

c          = -7: Length of private work array WORKL is not sufficient.

c          = -8: Error return from trid. eigenvalue calculation;

c                Informatinal error from LAPACK routine dsteqr .

c          = -9: Starting vector is zero.

c          = -10: IPARAM(7) must be 1,2,3,4,5.

c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.

c          = -12: IPARAM(1) must be equal to 0 or 1.

c          = -13: NEV and WHICH = 'BE' are incompatable.

c          = -14: IPARAM(4) is nonzero and the user returned the number of

c             converged Ritz values less than zero or greater and ncv.

c          = -9999: Could not build an Arnoldi factorization.

c                   IPARAM(5) returns the size of the current Arnoldi

c                   factorization. The user is advised to check that

c                   enough workspace and array storage has been allocated.

c

c VERBOSE : Flag to print out the status of the computation per iteration

c

c\Remarks

c  1. The converged Ritz values are always returned in ascending

c     algebraic order.  The computed Ritz values are approximate

c     eigenvalues of OP.  The selection of WHICH should be made

c     with this in mind when Mode = 3,4,5.  After convergence,

c     approximate eigenvalues of the original problem may be obtained

c     with the ARPACK subroutine pdseupd .

c

c  2. If the Ritz vectors corresponding to the converged Ritz values

c     are needed, the user must call pdseupd  immediately following completion

c     of pdsaupd . This is new starting with version 2.1 of ARPACK.

c

c  3. If M can be factored into a Cholesky factorization M = LL`

c     then Mode = 2 should not be selected.  Instead one should use

c     Mode = 1 with  OP = inv(L)*A*inv(L`).  Appropriate triangular

c     linear systems should be solved with L and L` rather

c     than computing inverses.  After convergence, an approximate

c     eigenvector z of the original problem is recovered by solving

c     L`z = x  where x is a Ritz vector of OP.

c

c  4. At present there is no a-priori analysis to guide the selection

c     of NCV relative to NEV.  The only formal requrement is that NCV > NEV.

c     However, it is recommended that NCV .ge. 2*NEV.  If many problems of

c     the same type are to be solved, one should experiment with increasing

c     NCV while keeping NEV fixed for a given test problem.  This will

c     usually decrease the required number of OP*x operations but it

c     also increases the work and storage required to maintain the orthogonal

c     basis vectors.   The optimal "cross-over" with respect to CPU time

c     is problem dependent and must be determined empirically.

c

c  5. If IPARAM(7) = 2 then in the Reverse commuication interface the user

c     must do the following. When IDO = 1, Y = OP * X is to be computed.

c     When IPARAM(7) = 2 OP = inv(B)*A. After computing A*X the user

c     must overwrite X with A*X. Y is then the solution to the linear set

c     of equations B*Y = A*X.

c

c  6. When IPARAM(1) = 0, and IDO = 3, the user needs to provide the

c     NP = IPARAM(8) shifts in locations:

c     1   WORKL(IPNTR(11))

c     2   WORKL(IPNTR(11)+1)

c                        .

c                        .

c                        .

c     NP  WORKL(IPNTR(11)+NP-1).

c

c     The eigenvalues of the current tridiagonal matrix are located in

c     WORKL(IPNTR(6)) through WORKL(IPNTR(6)+NCV-1). They are in the

c     order defined by WHICH. The associated Ritz estimates are located in

c     WORKL(IPNTR(8)), WORKL(IPNTR(8)+1), ... , WORKL(IPNTR(8)+NCV-1).

c

c  7. When IPARAM(4)=USRCHK is not equal to zero and IDO=4, then the

c     user will set IPARAM(5)=NCONV with the number of Ritz values that

c     satisfy the users' convergence criterion before pdsaupd.f is

c     re-entered. Note that IPNTR(5) contains the pointer into WORKL that

c     contains the NCV by 2 tridiagonal matrix T.

c

c-----------------------------------------------------------------------

c

c\Data Distribution Note:

c

c  Fortran-D syntax:

c  ================

c  REAL       RESID(N), V(LDV,NCV), WORKD(3*N), WORKL(LWORKL)

c  DECOMPOSE  D1(N), D2(N,NCV)

c  ALIGN      RESID(I) with D1(I)

c  ALIGN      V(I,J)   with D2(I,J)

c  ALIGN      WORKD(I) with D1(I)     range (1:N)

c  ALIGN      WORKD(I) with D1(I-N)   range (N+1:2*N)

c  ALIGN      WORKD(I) with D1(I-2*N) range (2*N+1:3*N)

c  DISTRIBUTE D1(BLOCK), D2(BLOCK,:)

c  REPLICATED WORKL(LWORKL)

c

c  Cray MPP syntax:

c  ===============

c  REAL       RESID(N), V(LDV,NCV), WORKD(N,3), WORKL(LWORKL)

c  SHARED     RESID(BLOCK), V(BLOCK,:), WORKD(BLOCK,:)

c  REPLICATED WORKL(LWORKL)

c

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.N. Parlett, "The Symmetric Eigenvalue Problem". Prentice-Hall,

c     1980.

c  4. B.N. Parlett, B. Nour-Omid, "Towards a Black Box Lanczos Program",

c     Computer Physics Communications, 53 (1989), pp 169-179.

c  5. B. Nour-Omid, B.N. Parlett, T. Ericson, P.S. Jensen, "How to

c     Implement the Spectral Transformation", Math. Comp., 48 (1987),

c     pp 663-673.

c  6. R.G. Grimes, J.G. Lewis and H.D. Simon, "A Shifted Block Lanczos

c     Algorithm for Solving Sparse Symmetric Generalized Eigenproblems",

c     SIAM J. Matr. Anal. Apps.,  January (1993).

c  7. L. Reichel, W.B. Gragg, "Algorithm 686: FORTRAN Subroutines

c     for Updating the QR decomposition", ACM TOMS, December 1990,

c     Volume 16 Number 4, pp 369-377.

c  8. R.B. Lehoucq, D.C. Sorensen, "Implementation of Some Spectral

c     Transformations in a k-Step Arnoldi Method". In Preparation.

c

c\Routines called:

c     pdsaup2   Parallel ARPACK routine that implements the Implicitly Restarted

c              Arnoldi Iteration.

c     dstats    ARPACK routine that initializes timing and other statistics

c              variables.

c     pivout   Parallel ARPACK utility routine that prints integers.

c     second   ARPACK utility routine for timing.

c     pdvout    Parallel ARPACK utility routine that prints vectors.

c     pdlamch   ScaLAPACK routine that determines machine constants.

c

c\Authors

c     Kristi Maschhoff ( Parallel Code )

c     Danny Sorensen               Phuong Vu

c     Richard Lehoucq              CRPC / Rice University

c     Dept. of Computational &     Houston, Texas

c     Applied Mathematics

c     Rice University

c     Houston, Texas

c

c\Parallel Modifications

c     Kristi Maschhoff

c

c\Revision history:

c     Starting Point: Serial Code FILE: saupd.F   SID: 2.4

c

c\SCCS Information:

c FILE: saupd.F   SID: 1.7   DATE OF SID: 04/10/01

c

c\Remarks

c     1. None

c

c\EndLib

c

c-----------------------------------------------------------------------

c

      subroutine mypdsaupd

     &   ( comm, ido, bmat, n, which, nev, tol, resid, ncv, v, ldv,

     &     iparam, ipntr, workd, workl, lworkl, info, verbose )

c

      include  'mpif.h'

c

c     %------------------%

c     | MPI Variables    |

c     %------------------%

c

      integer    comm, myid

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*1, which*2

      integer    ido, info, ldv, lworkl, n, ncv, nev

      integer    verbose

      Double precision

     &           tol

c

c     %-----------------%

c     | Array Arguments |

c     %-----------------%

c

      integer    iparam(11), ipntr(11)

      Double precision

     &           resid(n), v(ldv,ncv), workd(3*n), workl(lworkl)

c

c     %------------%

c     | Parameters |

c     %------------%

c

      Double precision

     &           one, zero

      parameter(one = 1.0 , zero = 0.0 )

c

c     %---------------%

c     | Local Scalars |

c     %---------------%

c

      integer    bounds, ierr, ih, iq, ishift, iupd, iw,

     &           ldh, ldq, msglvl, mxiter, mode, nb,

     &           nev0, next, np, ritz, j, usrchk, unconv

      save       bounds, ierr, ih, iq, ishift, iupd, iw,

     &           ldh, ldq, msglvl, mxiter, mode, nb,

     &           nev0, next, np, ritz, usrchk, unconv

c

c     %----------------------%

c     | External Subroutines |

c     %----------------------%

c

      external   mypdsaup2 , pdvout , pivout, second, dstats

c

c     %--------------------%

c     | External Functions |

c     %--------------------%

c

      Double precision

     &           pdlamch

      external   pdlamch

c

c     %-----------------------%

c     | Executable Statements |

c     %-----------------------%

c

      if (ido .eq. 0) then

c

c        %-------------------------------%

c        | Initialize timing statistics  |

c        | & message level for debugging |

c        %-------------------------------%

c

         call dstats

         call second (t0)

         msglvl = msaupd

c

         ierr   = 0

         ishift = iparam(1)

         mxiter = iparam(3)

c         nb     = iparam(4)

         nb     = 1

         usrchk = iparam(4)

         unconv = 0

c

c        %--------------------------------------------%

c        | Revision 2 performs only implicit restart. |

c        %--------------------------------------------%

c

         iupd   = 1

         mode   = iparam(7)

c

c        %----------------%

c        | Error checking |

c        %----------------%

c

         if (n .le. 0) then

            ierr = -1

         else if (nev .le. 0) then

            ierr = -2

         else if (ncv .le. nev) then

            ierr = -3

         end if

c

c        %----------------------------------------------%

c        | NP is the number of additional steps to      |

c        | extend the length NEV Lanczos factorization. |

c        %----------------------------------------------%

c

         np     = ncv - nev

c

         if (mxiter .le. 0)                     ierr = -4

         if (which .ne. 'LM' .and.

     &       which .ne. 'SM' .and.

     &       which .ne. 'LA' .and.

     &       which .ne. 'SA' .and.

     &       which .ne. 'BE')                   ierr = -5

         if (bmat .ne. 'I' .and. bmat .ne. 'G') ierr = -6

c

         if (lworkl .lt. ncv**2 + 8*ncv)        ierr = -7

         if (mode .lt. 1 .or. mode .gt. 5) then

                                                ierr = -10

         else if (mode .eq. 1 .and. bmat .eq. 'G') then

                                                ierr = -11

         else if (ishift .lt. 0 .or. ishift .gt. 1) then

                                                ierr = -12

         else if (nev .eq. 1 .and. which .eq. 'BE') then

                                                ierr = -13

         end if

c

c        %------------%

c        | Error Exit |

c        %------------%

c

         if (ierr .ne. 0) then

            info = ierr

            ido  = 99

            go to 9000

         end if

c

c        %------------------------%

c        | Set default parameters |

c        %------------------------%

c

         if (nb .le. 0) nb = 1

         if (tol .le. zero .and. usrchk .eq. 0)

     &  tol = pdlamch(comm, 'EpsMach')

c

c        %----------------------------------------------%

c        | NP is the number of additional steps to      |

c        | extend the length NEV Lanczos factorization. |

c        | NEV0 is the local variable designating the   |

c        | size of the invariant subspace desired.      |

c        %----------------------------------------------%

c

         np     = ncv - nev

         nev0   = nev

c

c        %-----------------------------%

c        | Zero out internal workspace |

c        %-----------------------------%

c

         do 10 j = 1, ncv**2 + 8*ncv

            workl(j) = zero

 10      continue

c

c        %-------------------------------------------------------%

c        | Pointer into WORKL for address of H, RITZ, BOUNDS, 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:2*ncv) := generated tridiagonal matrix        |

c        | workl(2*ncv+1:2*ncv+ncv) := ritz values               |

c        | workl(3*ncv+1:3*ncv+ncv) := computed error bounds     |

c        | workl(4*ncv+1:4*ncv+ncv*ncv) := rotation matrix Q     |

c        | workl(4*ncv+ncv*ncv+1:7*ncv+ncv*ncv) := workspace     |

c        %-------------------------------------------------------%

c

         ldh    = ncv

         ldq    = ncv

         ih     = 1

         ritz   = ih     + 2*ldh

         bounds = ritz   + ncv

         iq     = bounds + ncv

         iw     = iq     + ncv**2

         next   = iw     + 3*ncv

c

         ipntr(4) = next

         ipntr(5) = ih

         ipntr(6) = ritz

         ipntr(7) = bounds

         ipntr(11) = iw

      end if

c

  if (ido .eq. 4) then

c

c     %------------------------------------------------------------%

c     | Returning from reverse communication; iparam(5) contains   |

c     | the number of Ritz values that satisfy the users criterion |

c     %------------------------------------------------------------%

c

     unconv = iparam(5)

     if (unconv .lt. 0 .or. unconv .gt. ncv ) then

        ierr = -14

            info = ierr

            ido  = 99

            go to 9000

     end if

  end if

c

c     %-------------------------------------------------------%

c     | Carry out the Implicitly restarted Lanczos Iteration. |

c     %-------------------------------------------------------%

c

      call mypdsaup2

     &   ( comm, ido, bmat, n, which, nev0, np, tol, resid, mode, iupd,

     &     ishift, mxiter, v, ldv, workl(ih), ldh, workl(ritz),

     &     workl(bounds), workl(iq), ldq, workl(iw), ipntr, workd,

     &     info, verbose, usrchk, unconv )

c

c     %--------------------------------------------------%

c     | ido .ne. 99 implies use of reverse communication |

c     | to compute operations involving OP or shifts.    |

c     %--------------------------------------------------%

c

      if (ido .eq. 3) iparam(8) = np

      if (ido .ne. 99) go to 9000

c

      iparam(3) = mxiter

      iparam(5) = np

      iparam(9) = nopx

      iparam(10) = nbx

      iparam(11) = nrorth

c

c     %------------------------------------%

c     | Exit if there was an informational |

c     | error within pdsaup2 .              |

c     %------------------------------------%

c

      if (info .lt. 0) go to 9000

      if (info .eq. 2) info = 3

c

      if (msglvl .gt. 0) then

         call pivout (comm, logfil, 1, mxiter, ndigit,

     &               '_saupd: number of update iterations taken')

         call pivout (comm, logfil, 1, np, ndigit,

     &               '_saupd: number of "converged" Ritz values')

         call pdvout  (comm, logfil, np, workl(ritz), ndigit,

     &               '_saupd: final Ritz values')

         call pdvout  (comm, logfil, np, workl(bounds), ndigit,

     &               '_saupd: corresponding error bounds')

      end if

c

      call second (t1)

      tsaupd = t1 - t0

c

      if (msglvl .gt. 0) then

         call mpi_comm_rank( comm, myid, ierr )

         if ( myid .eq. 0 ) then

c

c        %--------------------------------------------------------%

c        | Version Number & Version Date are defined in version.h |

c        %--------------------------------------------------------%

c

         write (6,1000)

         write (6,1100) mxiter, nopx, nbx, nrorth, nitref, nrstrt,

     &                  tmvopx, tmvbx, tsaupd, tsaup2, tsaitr, titref,

     &                  tgetv0, tseigt, tsgets, tsapps, tsconv

 1000    format (//,

     &      5x, '==========================================',/

     &      5x, '= Symmetric implicit Arnoldi update code =',/

     &      5x, '= Version Number:', ' 2.1' , 19x, ' =',/

     &      5x, '= Version Date:  ', ' 3/19/97' , 14x, ' =',/

     &      5x, '==========================================',/

     &      5x, '= Summary of timing statistics           =',/

     &      5x, '==========================================',//)

 1100    format (

     &      5x, 'Total number update iterations             = ', i5,/

     &      5x, 'Total number of OP*x operations            = ', i5,/

     &      5x, 'Total number of B*x operations             = ', i5,/

     &      5x, 'Total number of reorthogonalization steps  = ', i5,/

     &      5x, 'Total number of iterative refinement steps = ', i5,/

     &      5x, 'Total number of restart steps              = ', i5,/

     &      5x, 'Total time in user OP*x operation          = ', f12.6,/

     &      5x, 'Total time in user B*x operation           = ', f12.6,/

     &      5x, 'Total time in Arnoldi update routine       = ', f12.6,/

     &      5x, 'Total time in p_saup2 routine              = ', f12.6,/

     &      5x, 'Total time in basic Arnoldi iteration loop = ', f12.6,/

     &      5x, 'Total time in reorthogonalization phase    = ', f12.6,/

     &      5x, 'Total time in (re)start vector generation  = ', f12.6,/

     &      5x, 'Total time in trid eigenvalue subproblem   = ', f12.6,/

     &      5x, 'Total time in getting the shifts           = ', f12.6,/

     &      5x, 'Total time in applying the shifts          = ', f12.6,/

     &      5x, 'Total time in convergence testing          = ', f12.6)

         end if

      end if

c

 9000 continue

c

      return

c

c     %----------------%

c     | End of pdsaupd  |

c     %----------------%

c

      end