c\BeginDoc c c\Name: znaupd c c\Description: c Reverse communication interface for the Implicitly Restarted Arnoldi c iteration. This is intended to be used to find a few eigenpairs of a c complex linear operator OP with respect to a semi-inner product defined c by a hermitian positive semi-definite real matrix B. B may be the identity c matrix. NOTE: if both OP and B are real, then dsaupd or dnaupd should c be used. c c c The computed approximate eigenvalues are called Ritz values and c the corresponding approximate eigenvectors are called Ritz vectors. c c znaupd is usually called iteratively to solve one of the c following problems: c c Mode 1: A*x = lambda*x. c ===> OP = A and B = I. c c Mode 2: A*x = lambda*M*x, 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: A*x = lambda*M*x, M symmetric semi-definite c ===> OP = inv[A - sigma*M]*M and B = M. c ===> shift-and-invert mode c If OP*x = amu*x, then lambda = sigma + 1/amu. c 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 znaupd c ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, c IPNTR, WORKD, WORKL, LWORKL, RWORK, INFO ) c c\Arguments c IDO Integer. (INPUT/OUTPUT) c Reverse communication flag. IDO must be zero on the first c call to znaupd. 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 znaupd with the result. The operand is given in c WORKD(IPNTR(1)), the result must be put in WORKD(IPNTR(2)). 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 * Z and Z = B * X where c IPNTR(1) is the pointer into WORKD for X, c IPNTR(2) is the pointer into WORKD for Y, c IPNTR(3) is the pointer into WORKD for Z. c IDO = 2: compute Y = M * 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 and return the shifts in the first c NP locations of WORKL. c IDO = 4: compute Z = OP * X c IDO = 99: done c ------------------------------------------------------------- c After the initialization phase, when the routine is used in c the "shift-and-invert" mode, the vector M * X is already c available and does not need to be recomputed in forming OP*X. 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 BMAT = 'I' -> standard eigenvalue problem A*x = lambda*x c BMAT = 'G' -> generalized eigenvalue problem A*x = lambda*M*x c c N Integer. (INPUT) c Dimension of the eigenproblem. c c WHICH Character*2. (INPUT) c 'LM' -> want the NEV eigenvalues of largest magnitude. c 'SM' -> want the NEV eigenvalues of smallest magnitude. c 'LR' -> want the NEV eigenvalues of largest real part. c 'SR' -> want the NEV eigenvalues of smallest real part. c 'LI' -> want the NEV eigenvalues of largest imaginary part. c 'SI' -> want the NEV eigenvalues of smallest imaginary part. c c NEV Integer. (INPUT) c Number of eigenvalues of OP to be computed. 0 < NEV < N-1. c c TOL Double precision scalar. (INPUT) c Stopping criteria: the relative accuracy of the Ritz value c is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I)) c where ABS(RITZ(I)) is the magnitude when RITZ(I) is complex. c DEFAULT = dlamch('EPS') (machine precision as computed c by the LAPACK auxiliary subroutine dlamch). c c RESID Complex*16 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. NCV must satisfy the two c inequalities 2 <= NCV-NEV and NCV <= N. c This will indicate how many Arnoldi vectors are generated c at each iteration. After the startup phase in which NEV c Arnoldi vectors are generated, the algorithm generates c approximately NCV-NEV Arnoldi vectors at each subsequent update c iteration. Most of the cost in generating each Arnoldi vector is c in the matrix-vector operation OP*x. c NOTE: 2 <= NCV-NEV in order that complex conjugate pairs of Ritz c values are kept together. (See remark 4 below) c c V Complex*16 array N by NCV. (OUTPUT) c Contains the final set of Arnoldi basis vectors. c c LDV Integer. (INPUT) c Leading dimension of V exactly as declared in the calling 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 filter out c the components of the unwanted eigenvector. c ------------------------------------------------------------- c ISHIFT = 0: the shifts are to be provided by the user via c reverse communication. The NCV eigenvalues of c the Hessenberg matrix H are returned in the part c of WORKL array corresponding to RITZ. c ISHIFT = 1: exact shifts with respect to the current c Hessenberg matrix H. This is equivalent to c restarting the iteration from the beginning c after updating the starting vector with a linear c combination of Ritz vectors associated with the c "wanted" eigenvalues. c ISHIFT = 2: other choice of internal shift to be defined. c ------------------------------------------------------------- c c IPARAM(2) = No longer referenced 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) = NB: blocksize to be used in the recurrence. c The code currently works only for NB = 1. c c IPARAM(5) = NCONV: number of "converged" Ritz values. c This represents the number of Ritz values that satisfy c the 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; See under \Description of znaupd for the c four modes available. c c IPARAM(8) = NP c When ido = 3 and the user provides shifts through reverse c communication (IPARAM(1)=0), _naupd returns NP, the number c of shifts the user is to provide. 0 < NP < NCV-NEV. 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 14. (OUTPUT) c Pointer to mark the starting locations in the WORKD and WORKL c arrays for matrices/vectors used by the Arnoldi 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 NCV upper Hessenberg c matrix H in WORKL. c IPNTR(6): pointer to the ritz value array RITZ c IPNTR(7): pointer to the (projected) ritz vector array Q c IPNTR(8): pointer to the error BOUNDS array in WORKL. c Note: IPNTR(9:13) is only referenced by zneupd. See Remark 2 below. 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 bounds. c IPNTR(14): pointer to the NP shifts in WORKL. See Remark 5 below. c ------------------------------------------------------------- c c WORKD Complex*16 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 !!!!!!!!!! c See Data Distribution Note below. c c WORKL Complex*16 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 3*NCV**2 + 5*NCV. c c RWORK Double precision work array of length NCV (WORKSPACE) c Private (replicated) array on each PE or array allocated on c the front end. c 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-NEV >= 2 and less than or equal to N. c = -4: The maximum number of Arnoldi update iteration c must be greater than zero. 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 array is not sufficient. c = -8: Error return from LAPACK eigenvalue calculation; c = -9: Starting vector is zero. c = -10: IPARAM(7) must be 1,2,3. c = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable. c = -12: IPARAM(1) must be equal to 0 or 1. c = -9999: Could not build an Arnoldi factorization. c User input error highly likely. Please c check actual array dimensions and layout. c IPARAM(5) returns the size of the current Arnoldi c factorization. c c\Remarks c 1. The computed Ritz values are approximate eigenvalues of OP. The c selection of WHICH should be made with this in mind when using c Mode = 3. When operating in Mode = 3 setting WHICH = 'LM' will c compute the NEV eigenvalues of the original problem that are c closest to the shift SIGMA . After convergence, approximate eigenvalues c of the original problem may be obtained with the ARPACK subroutine zneupd. c c 2. If a basis for the invariant subspace corresponding to the converged Ritz c values is needed, the user must call zneupd immediately following c completion of znaupd. This is new starting with release 2 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 of NCV c relative to NEV. The only formal requirement is that NCV > NEV + 2. c However, it is recommended that NCV .ge. 2*NEV+1. 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 See Chapter 8 of Reference 2 for further information. c c 5. When IPARAM(1) = 0, and IDO = 3, the user needs to provide the c NP = IPARAM(8) complex shifts in locations c WORKL(IPNTR(14)), WORKL(IPNTR(14)+1), ... , WORKL(IPNTR(14)+NP). c Eigenvalues of the current upper Hessenberg matrix are located in c WORKL(IPNTR(6)) through WORKL(IPNTR(6)+NCV-1). They are ordered c according to the order defined by WHICH. The associated Ritz estimates c are located in WORKL(IPNTR(8)), WORKL(IPNTR(8)+1), ... , c WORKL(IPNTR(8)+NCV-1). c c-----------------------------------------------------------------------