Zgsev Zgsgv Znaupd Zneupd

Znaupd() Sub Znaupd ( Ido As  Long, Bmat As  String, N As  Long, Which As  String, Nev As  Long, Resid() As  Complex, Ncv As  Long, V() As  Complex, IParam() As  Long, Ipntr() As  Long, Workd() As  Complex, Workl() As  Complex, RWork() As  Double, Info As  Long, Optional Tol As  Double = 0  ) Arnoldi factorization of a complex sparse matrix (Arpack) Purpose The eigenvalues and eigenvectors of a complex sparse matrix are computed by Implicitly restarted Arnoldi method (IRAM). This routine computes the Arnoldi factorization, which is the first half of the algorithm. The approximate eigenvalues computed from the factorization are called Ritz values and the corresponding approximate eigenvectors are called Ritz vectors. The routine Zneupd computes the Ritz values and the Ritz vectors, which is the second half of the algorithm, from the output of this routine. This routine computes approximations to a few eigenpairs of a complex linear operator "OP" with respect to a semi-inner product defined by a Hermitian positive semi-definite matrix B. This routine is iteratively called until the convergence to solve one of the following problems. Regular mode (mode = 1): A*x = λ*x. (OP = A, B = I) Regular inverse mode (mode = 2): A*x = λ*M*x, M is Hermitian positive definite. (OP = M^(-1)*A, B = M) Shift-and-invert mode (mode = 3): A*x = λ*M*x, M is Hermitian semi-definite. (OP = (A - σ*M)^(-1)*M, B = M) Note - mode is soecified by Iparam(6). Parameters [in,out] Ido Reverse communication flag. [in] Ido must be zero on the first call. [out] Ido will be set to indicate the type of operation to be performed. Control is then given back to the calling routine which has the responsibility to carry out the requested operation and call again with the result. = -1: Compute y = OP * x. (for the initialization phase) = 1: Compute y = OP * x. = 2: Compute y = B * x. = 3: Compute np shifts. (Note 3) = 99: done. Note 1 - x is Workd(Ipntr(0) - 1), Workd(Ipntr(0)), ..., Workd(Ipntr(0) + n - 2). Note 2 - y is Workd(Ipntr(1) - 1), Workd(Ipntr(1)), ..., Workd(Ipntr(1) + n - 2). Note 3 - When Iparam(0) = 0, and Ido = 3, the user needs to provide the np (= Iparam(7)) shifts in locations: Workl(Ipntr(10) - 1), Workl(Ipntr(10)), ..., Workl(Ipntr(10) + np - 2). [in] Bmat = "I": Standard eigenvalue problem A*x = λ*x = "G": Generalized eigenvalue problem A*x = λ*B*x [in] N Dimension of the matrix. (N >= 0) (If N = 0, returns without computation) [in] Which = "LM": Compute the Nev eigenvalues of largest magnitude. = "SM": Compute the Nev eigenvalues of smallest magnitude. = "LR": Compute the Nev eigenvalues of largest real part. = "SR": Compute the Nev eigenvalues of smallest real part. = "LI": Compute the Nev eigenvalues of largest imaginary part. = "SI": Compute the Nev eigenvalues of smallest imaginary part. [in] Nev Number of eigenvalues to be computed. (0 < Nev < N) [in] Tol Stopping criterion: the acceptable relative accuracy of the Ritz value. If Tol <= 0, machine epsilon is assumed. [in,out] Resid() Array Resid(LResid - 1) (LResid >= N) [in] If Info <> 0, the initial residual vector possibly from a previous run. [out] The final residual vector. [in] Ncv Number of columns of the matrix V. (Nev < Ncv <= N) This will indicate how many Arnoldi vectors are generated at each iteration. [out] V() Array V(Ldv - 1, LV - 1) (Ldv >= N, LV >= Ncv) Contains the final set of Arnoldi basis vectors. [in,out] Iparam() Array Iparam(LIparam - 1) (LIparam >= 11) Iparam(0) (ishift): [in] Method for selecting the implicit shifts. = 0: The shifts are provided by the user via reverse communication. = 1: Exact shifts with respect to the current Hessenberg matrix H. Iparam(1) Not used. Iparam(2) (mxiter): [in] Maximum number of Arnoldi update iterations allowed. [out] Actual number of iterations taken. Iparam(3) (nb): [in] Blocksize to be used in the recurrence. (currently nb must be 1) Iparam(4) (nconv): [out] The number of Ritz values that satisfy the convergence criterion. Iparam(5) Not used. Iparam(6) (mode): [in] Type of eigenvalue problem (mode = 1, ..., 5) Iparam(7) (np): [out] When Ido = 3 and Iparam(0) = 0, returns the number of shifts the user is to provide. (0 < np <= Ncv - Nev) Iparam(8) (numop): [out] Total number of OP*x operations. Iparam(9) (numopb): [out] Total number of B*x operations if Bmat = "G". Iparam(10) (numreo): [out] Total number of steps of re-orthogonalization. [out] Ipntr() Array Ipntr(LIpntr - 1) (LIpntr >= 14) Pointer to mark the starting locations in the work array for matrices/vectors used by the Arnoldi iteration (1-based index). Ipntr(0): Pointer to the current operand vector x in Workd. Ipntr(1): Pointer to the current result vector y in Workd. Ipntr(2): Pointer to the vector B*x in Workd when used in the shift-and-invert mode. Ipntr(3): Pointer to the next available location in Workl that is untouched by the program. Ipntr(4): Pointer to the Ncv by Ncv upper Hessenberg matrix H in Workl. Ipntr(5): Pointer to the real part of the Ritz value array in Workl. Ipntr(6): Pointer to the imaginary part of the Ritz value array in Workl. Ipntr(7): Pointer to the Ritz estimates in array Workl associated with the Ritz values in Workl. Ipntr(8), ..., Ipntr(12): Referenced only by Dneupd. Ipntr(13): Pointer to the np shifts in Workl. [out] Workd() Array Workd(LWorkd - 1) (LWorkd >= 3*N) Distributed work array used for reverse communication. [out] Workl() Array Workl(LWorkl - 1) (LWorkl >= 3*Ncv^2 + 5*Ncv) Local work array. [out] RWork() Array RWork[LRWork] (LRWork >= Ncv) Private (replicated) work array. [in,out] Info [in] = 0: A random initial residual vector is used. <> 0: Resid() contains the initial residual vector, possibly from a previous run. [out] Return code. = 0: Successful exit. < 0: The (-Info)-th argument is invalid. = 1: Maximum number of iterations taken. Iparam(4) returns the number of converged Ritz values. = 3: No shifts could be applied during a cycle of the implicitly restarted Arnoldi iteration. A possible remedy is to increase Ncv relative to Nev (Ncv >= 2*Nev is recommended). = 11: Initial residual vector is zero. = 12: Failed to build a Arnoldi factorization. Iparam(4) returns the size of the current factorization. = 13: Error return from LAPACK eigenvalue calculation. Reference ARPACK-NG (https://github.com/opencollab/arpack-ng) Example Program Computes the largest eigenvalue and its eigenvector of the matrix A where ( 0.20-0.11i -0.93-0.32i 0.81+0.37i ) A = ( -0.80-0.92i -0.29+0.86i 0.64+0.51i ) ( 0.71+0.59i -0.15+0.19i 0.20+0.94i ) Sub Ex_Znaupd() Const N = 3, Nnz = 9, Ncv = 3, Nev = 2 Const Rvec = 1, Bmat = "I", Which = "LM" Dim A(Nnz - 1) As Complex, Ia(N) As Long, Ja(Nnz - 1) As Long Dim D(Nev - 1) As Complex, Z(N - 1, Nev - 1) As Complex Dim V(N - 1, Ncv - 1) As Complex, Resid(N - 1) As Complex Dim Workd(3 * N - 1) As Complex, Workl(3 * Ncv * Ncv + 5 * Ncv - 1) As Complex Dim Workev(3 * Ncv - 1) As Complex, XX(N - 1) As Complex, YY(N - 1) As Complex Dim Selct(Ncv - 1) As Long, IParam(10) As Long, Ipntr(13) As Long Dim Sigma As Complex, RWork(Ncv - 1) As Double Dim Ido As Long, Info As Long A(0) = Cmplx(0.2, -0.11): A(1) = Cmplx(-0.93, -0.32): A(2) = Cmplx(0.81, 0.37): A(3) = Cmplx(-0.8, -0.92): A(4) = Cmplx(-0.29, 0.86): A(5) = Cmplx(0.64, 0.51): A(6) = Cmplx(0.71, 0.59): A(7) = Cmplx(-0.15, 0.19): A(8) = Cmplx(0.2, 0.94) Ia(0) = 0: Ia(1) = 3: Ia(2) = 6: Ia(3) = 9 Ja(0) = 0: Ja(1) = 1: Ja(2) = 2: Ja(3) = 0: Ja(4) = 1: Ja(5) = 2: Ja(6) = 0: Ja(7) = 1: Ja(8) = 2 IParam(0) = 1: IParam(2) = 300: IParam(3) = 1: IParam(6) = 1 Info = 0 Ido = 0 Do Call Znaupd(Ido, Bmat, N, Which, Nev, Resid(), Ncv, V(), IParam(), Ipntr(), Workd(), Workl(), RWork(), Info) If Ido = -1 Or Ido = 1 Then Call ZCopy(N, Workd(Ipntr(0) - 1), XX(0)) Call CsrZusmv("N", N, N, Cmplx(1), A(), Ia(), Ja(), XX(), Cmplx(0), YY()) Call ZCopy(N, YY(0), Workd(Ipntr(1) - 1)) End If Loop While (Ido <> 99) If Info <> 0 Then Debug.Print "Znaupd Info =" + Str(Info) Call Zneupd(Rvec, "A", Selct(), D(), Z(), Sigma, Workev(), Bmat, N, Which, Nev, Resid(), Ncv, V(), IParam(), Ipntr(), Workd(), Workl(), RWork(), Info) If Info <> 0 Then Debug.Print "Zneupd Info =" + Str(Info) Debug.Print "D =" Debug.Print Creal(D(0)), Cimag(D(0)), Creal(D(1)), Cimag(D(1)) Debug.Print "Z =" Debug.Print Creal(Z(0, 0)), Cimag(Z(0, 0)), Creal(Z(0, 1)), Cimag(Z(0, 1)) Debug.Print Creal(Z(1, 0)), Cimag(Z(1, 0)), Creal(Z(1, 1)), Cimag(Z(1, 1)) Debug.Print Creal(Z(2, 0)), Cimag(Z(2, 0)), Creal(Z(2, 1)), Cimag(Z(2, 1)) Debug.Print "Nconv =" + Str(IParam(4)) + ", Info =" + Str(Info) End Sub Example Results D = 1.05593587167591 0.900255855387814 0.213005352563271 1.29637306909393 Z = -0.626960983365491 0.172271418863805 7.9679684322005E-03 0.150485895565997 0.268918491094688 -5.24016295711404E-02 0.3557286693982 0.561279089020254 -0.662645215037481 -0.251158748471167 0.241195085578284 0.691041486022676 Nconv = 2, Info = 0