Dgsev Dgsgv Dnaupd Dneupd Dsaupd Dseupd Dssev Dssgv

Dnaupd() Sub Dnaupd ( Ido As  Long, Bmat As  String, N As  Long, Which As  String, Nev As  Long, Resid() As  Double, Ncv As  Long, V() As  Double, IParam() As  Long, Ipntr() As  Long, Workd() As  Double, Workl() As  Double, Info As  Long, Optional Tol As  Double = 0  ) Arnoldi factorization of a general sparse matrix (Arpack) Purpose The eigenvalues and eigenvectors of a general 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 Dneupd 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 linear operator "OP" with respect to a semi-inner product defined by a symmetric positive semi-definite real matrix B. Note - If the linear operator "OP" is real and symmetric with respect to the real positive semi-definite symmetric matrix B, i.e. B*OP = (OP^T)*B, then subroutine Dsaupd should be used instead. 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 symmetric positive definite. (OP = M^(-1)*A, B = M) Shift-and-invert mode (in real arithmetic) (mode = 3): A*x = λ*M*x, M is symmetric semi-definite. (OP = real part of (A - σ*M)^(-1)*M, B = M) Shift-and-invert mode (in real arithmetic) (mode = 4): A*x = λ*M*x, M is symmetric semi-definite. (OP = imaginary part of (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 - 1) [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 + 2 < 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 + 6*Ncv) Local 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.2 -0.11 -0.93 ) A = ( -0.32 0.81 0.37 ) ( -0.8 -0.92 -0.29 ) Sub Ex_Dnaupd() Const N = 3, Nnz = 9, Ncv = N, Nev = Ncv - 2 Const Rvec = 1, Bmat = "I", Which = "LM" Dim A(Nnz - 1) As Double, Ia(N) As Long, Ja(Nnz - 1) As Long Dim Dr(Nev) As Double, Di(Nev) As Double, Z(N - 1, Nev) As Double Dim V(N - 1, Ncv - 1) As Double, Resid(N - 1) As Double Dim Workd(3 * N - 1) As Double, Workl(3 * Ncv * (Ncv + 2) - 1) As Double Dim Workev(3 * Ncv - 1) As Double, XX(N - 1) As Double, YY(N - 1) As Double Dim Selct(Ncv - 1) As Long, IParam(10) As Long, Ipntr(13) As Long Dim Sigmar As Double, Sigmai As Double Dim Ido As Long, Nconv As Long, Info As Long A(0) = 0.2: A(1) = -0.11: A(2) = -0.93: A(3) = -0.32: A(4) = 0.81: A(5) = 0.37: A(6) = -0.8: A(7) = -0.92: A(8) = -0.29 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 Dnaupd(Ido, Bmat, N, Which, Nev, Resid(), Ncv, V(), IParam(), Ipntr(), Workd(), Workl(), Info) If Ido = -1 Or Ido = 1 Then Call DCopy(N, Workd(Ipntr(0) - 1), XX(0)) Call CsrDusmv("N", N, N, 1, A(), Ia(), Ja(), XX(), 0, YY()) Call DCopy(N, YY(0), Workd(Ipntr(1) - 1)) End If Loop While (Ido <> 99) Nconv = IParam(4) If Info <> 0 Then Debug.Print "Dnaupd Info =" + Str(Info) Call Dneupd(Rvec, "A", Selct(), Dr(), Di(), Z(), Sigmar, Sigmai, Workev(), Bmat, N, Which, Nev, Resid(), Ncv, V(), IParam(), Ipntr(), Workd(), Workl(), Info) If Info <> 0 Then Debug.Print "Dneupd Info =" + Str(Info) Debug.Print "D =", Dr(0), Di(0) Debug.Print "Z =" Debug.Print Z(0, 0), Z(0, 1) Debug.Print Z(1, 0), Z(1, 1) Debug.Print Z(2, 0), Z(2, 1) Debug.Print "Nconv =" + Str(Nconv) End Sub Example Results D = 0.812065011925673 0.48915757543818 Z = -0.365053973674333 -0.416255621362224 0.64300369896322 -0.194693188992647 -5.47400351005566E-02 0.488989966998924 Nconv = 2