Dgsev Dgsgv Dnaupd Dneupd Dsaupd Dseupd Dssev Dssgv

Dsaupd() Sub Dsaupd ( 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  ) Lanczos factorization of a symmetric sparse matrix (Arpack) Purpose The eigenvalues and eigenvectors of a real symmetric sparse matrix are computed by Implicitly restarted Lanczos method (IRLM). This routine computes the Lanczos factorization, which is the first half of the algorithm. The approximate eigenvalues computed from the Lanczos factorization are called Ritz values and the corresponding approximate eigenvectors are called Ritz vectors. The routine Dseupd 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 that is real and symmetric with respect to a real positive semi-definite symmetric matrix B, i.e. B*OP = (OP^T)*B or <x, OPy> = <OPx, y> where <z, w> = z^T*B*w. This routine is iteratively called until the convergence to solve one of the following problems. Regular mode (mode = 1): A*x = λ*x, A is symmetric. (OP = A, B = I) Regular inverse mode (mode = 2): A*x = λ*M*x, A is symmetric, M is symmetric positive definite. (OP = M^(-1)*A, B = M) Shift-and-invert mode (mode = 3): K*x = λ*M*x, K is symmetric, M is symmetric positive semi-definite. (OP = (K - σ*M)^(-1)*M, B = M) Buckling mode (mode = 4): K*x = λ*KG*x, K is symmetric positive semi-definite, KG is symmetric indefinite. (OP = (K - σ*KG)^(-1)*K, B = K) Cayley transformed mode (mode = 5): A*x = λ*M*x, A is symmetric, M is symmetric positive semi-definite. (OP = (A - σ*M)^(-1)*(A + σ*M), B = M) Note 1 - mode is soecified by Iparam(6). Note 2 - If M can be factored into a Cholesky factorization M = L*L^T then mode = 2 should not be selected. Instead one should use mode = 1 with OP = L^(-1)*A*L^(-T). After convergence, an approximate eigenvector z of the original problem is recovered by solving L^T*z = x where x is a Ritz vector of OP. 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 = "LA": Compute the Nev largest (algebraic) eigenvalues. = "SA": Compute the Nev smallest (algebraic) eigenvalues. = "LM": Compute the Nev largest (in magnitude) eigenvalues. = "SM": Compute the Nev smallest (in magnitude) eigenvalues. = "BE": Compute Nev eigenvalues, half from each end of the spectrum. When Nev is odd, compute one more from the high end than from the low end. [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 precision 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 Lanczos vectors are generated at each iteration. [out] V() Array V(Ldv - 1, LV - 1) (Ldv >= N, LV >= Ncv) Matrix V containing Ncv Lanczos 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 reduced tridiagonal matrix T. Iparam(1) Not used. Iparam(2) (mxiter): [in] Maximum number of Lanczos 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 >= 11) Pointer to mark the starting locations in the work array for matrices/vectors used by the Lanczos iteration (1-based). 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 2 tridiagonal matrix T in Workl. Ipntr(5): Pointer to the Ncv Ritz values array in Workl. Ipntr(6): Pointer to the Ritz estimates in array Workl associated with the Ritz values located in Ritz in Workl. Ipntr(7), Ipntr(8), Ipntr(9): Referenced only by Dseupd. Ipntr(10): 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 >= Ncv^2 + 8*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 Lanczos 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 Lanczos factorization. Iparam(4) returns the size of the current factorization. = 13: Error return from tridiagonal eigenvalue calculation (LAPACK routine DSTEQR). Reference ARPACK-NG (https://github.com/opencollab/arpack-ng) Example Program Computes largest two eigenvalues and their eigenvectors of the symmetric matrix A where ( 2.2 -0.11 -0.32 ) A = ( -0.11 2.93 0.81 ) ( -0.32 0.81 2.37 ) Sub Ex_Dsaupd() Const N = 3, Nnz = 6, Ncv = N, Nev = Ncv - 1 Const Rvec = 1, Bmat = "I", Which = "LM" Dim A(Nnz - 1) As Double, Ia(N) As Long, Ja(Nnz - 1) As Long Dim D(Nev - 1) As Double, Z(N - 1, Nev - 1) As Double Dim V(N - 1, Ncv - 1) As Double, Resid(N - 1) As Double Dim Workd(3 * N - 1) As Double, Workl(Ncv * (Ncv + 8) - 1) As Double Dim XX(N - 1) As Double, YY(N - 1) As Double Dim Selct(Ncv - 1) As Long, IParam(10) As Long, Ipntr(10) As Long, Sigma As Double Dim Ido As Long, Nconv As Long, Info As Long A(0) = 2.2: A(1) = -0.11: A(2) = 2.93: A(3) = -0.32: A(4) = 0.81: A(5) = 2.37 Ia(0) = 0: Ia(1) = 1: Ia(2) = 3: Ia(3) = 6 Ja(0) = 0: Ja(1) = 0: Ja(2) = 1: Ja(3) = 0: Ja(4) = 1: Ja(5) = 2 IParam(0) = 1: IParam(2) = 300: IParam(3) = 1: IParam(6) = 1 Info = 0 Ido = 0 Do Call Dsaupd(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 SsrDusmv("L", 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 "Dsaupd Info =" + Str(Info) Call Dseupd(Rvec, "A", Selct(), D(), Z(), Sigma, Bmat, N, Which, Nev, Resid(), Ncv, V(), IParam(), Ipntr(), Workd(), Workl(), Info) If Info <> 0 Then Debug.Print "Dseupd Info =" + Str(Info) Debug.Print "D =", D(0), D(1) 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 = 2.22943643244226 3.56350401844728 Z = -0.894521385341403 -0.200931256445799 -0.390994588677271 0.78468897753835 0.216690384632057 0.586421212707156 Nconv = 2