|
|
◆ dsaupd()
| void dsaupd |
( |
int * |
ido, |
|
|
char |
bmat, |
|
|
int |
n, |
|
|
const char * |
which, |
|
|
int |
nev, |
|
|
double |
tol, |
|
|
double |
resid[], |
|
|
int |
ncv, |
|
|
int |
ldv, |
|
|
double |
v[], |
|
|
int |
iparam[], |
|
|
int |
ipntr[], |
|
|
double |
workd[], |
|
|
double |
workl[], |
|
|
int |
lworkl, |
|
|
int * |
info |
|
) |
| |
Lanczos factorization of a symmetric sparse matrix
- 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) (Note 4)
= 1: Compute y = OP * x. (Note 4)
= 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].
Note 4 - When mode = 2, after computing OP * x the user must overwrite x with OP * x. |
| [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] (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. |
| [in] | ldv | Leading dimension of the array v. (ldv >= n) |
| [out] | v[] | Array v[ldv * lv] (lv >= ncv)
Matrix V containing ncv Lanczos basis vectors. |
| [in,out] | iparam[] | Array iparam[liparam] (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] (lipntr >= 11)
Pointer to mark the starting locations in the work array for matrices/vectors used by the Lanczos 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 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] (lworkd >= 3*n)
Distributed work array used for reverse communication. |
| [out] | workl[] | Array workl[lworkl]
Local work array. |
| [in] | lworkl | Size of array workl[]. (lworkl >= ncv^2 + 8*ncv) |
| [in,out] | info | [in]
= 0: A random initial residual vector is used.
!= 0: resid[] contains the initial residual vector.
[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 LAPACK eigenvalue calculation. |
- Reference
- ARPACK-NG (https://github.com/opencollab/arpack-ng)
|