Function File: [`blockVectorX`, `lambda`] = **lobpcg** (`blockVectorX, operatorA`)

Function File: [`blockVectorX`, `lambda`, `failureFlag`] = **lobpcg** (`blockVectorX, operatorA`)

Function File: [`blockVectorX`, `lambda`, `failureFlag`, `lambdaHistory`, `residualNormsHistory`] = **lobpcg** (`blockVectorX, operatorA, operatorB, operatorT, blockVectorY, residualTolerance, maxIterations, verbosityLevel`)

Solves Hermitian partial eigenproblems using preconditioning.

The first form outputs the array of algebraic smallest eigenvalues

lambdaand corresponding matrix of orthonormalized eigenvectorsblockVectorXof the Hermitian (full or sparse) operatoroperatorAusing input matrixblockVectorXas an initial guess, without preconditioning, somewhat similar to:# for real symmetric operator operatorA opts.issym = 1; opts.isreal = 1; K = size (blockVectorX, 2); [blockVectorX, lambda] = eigs (operatorA, K, 'SR', opts); # for Hermitian operator operatorA K = size (blockVectorX, 2); [blockVectorX, lambda] = eigs (operatorA, K, 'SR');The second form returns a convergence flag. If

failureFlagis 0 then all the eigenvalues converged; otherwise not all converged.The third form computes smallest eigenvalues

lambdaand corresponding eigenvectorsblockVectorXof the generalized eigenproblem Ax=lambda Bx, where Hermitian operatorsoperatorAandoperatorBare given as functions, as well as a preconditioner,operatorT. The operatorsoperatorBandoperatorTmust be in additionpositive definite. To compute the largest eigenpairs ofoperatorA, simply apply the code tooperatorAmultiplied by -1. The code does not involveanymatrix factorizations ofoperatorAandoperatorB, thus, e.g., it preserves the sparsity and the structure ofoperatorAandoperatorB.

residualToleranceandmaxIterationscontrol tolerance and max number of steps, andverbosityLevel= 0, 1, or 2 controls the amount of printed info.lambdaHistoryis a matrix with all iterative lambdas, andresidualNormsHistoryare matrices of the history of 2-norms of residualsRequired input:

blockVectorX(class numeric) - initial approximation to eigenvectors, full or sparse matrix n-by-blockSize.blockVectorXmust be full rank.operatorA(class numeric, char, or function_handle) - the main operator of the eigenproblem, can be a matrix, a function name, or handleOptional function input:

operatorB(class numeric, char, or function_handle) - the second operator, if solving a generalized eigenproblem, can be a matrix, a function name, or handle; by default if empty,`operatorB = I`

.operatorT(class char or function_handle) - the preconditioner, by default`operatorT(blockVectorX) = blockVectorX`

.Optional constraints input:

blockVectorY(class numeric) - a full or sparse n-by-sizeY matrix of constraints, where sizeY < n.blockVectorYmust be full rank. The iterations will be performed in the (operatorB-) orthogonal complement of the column-space ofblockVectorY.Optional scalar input parameters:

residualTolerance(class numeric) - tolerance, by default,`residualTolerance = n * sqrt (eps)`

maxIterations- max number of iterations, by default,`maxIterations = min (n, 20)`

verbosityLevel- either 0 (no info), 1, or 2 (with pictures); by default,`verbosityLevel = 0`

.Required output:

blockVectorXandlambda(class numeric) both are computed blockSize eigenpairs, where`blockSize = size (blockVectorX, 2)`

for the initial guessblockVectorXif it is full rank.Optional output:

failureFlag(class integer) as described above.lambdaHistory(class numeric) as described above.residualNormsHistory(class numeric) as described above.Functions

`operatorA(blockVectorX)`

,`operatorB(blockVectorX)`

and`operatorT(blockVectorX)`

must supportblockVectorXbeing a matrix, not just a column vector.Every iteration involves one application of

operatorAandoperatorB, and one ofoperatorT.Main memory requirements: 6 (9 if

`isempty(operatorB)=0`

) matrices of the same size asblockVectorX, 2 matrices of the same size asblockVectorY(if present), and two square matrices of the size 3*blockSize.In all examples below, we use the Laplacian operator in a 20x20 square with the mesh size 1 which can be generated in MATLAB by running:

A = delsq (numgrid ('S', 21)); n = size (A, 1);or in MATLAB and Octave by:

[~,~,A] = laplacian ([19, 19]); n = size (A, 1);Note that

`laplacian`

is a function of the specfun octave-forge package.The following Example:

[blockVectorX, lambda, failureFlag] = lobpcg (randn (n, 8), A, 1e-5, 50, 2);attempts to compute 8 first eigenpairs without preconditioning, but not all eigenpairs converge after 50 steps, so failureFlag=1.

The next Example:

blockVectorY = []; lambda_all = []; for j = 1:4 [blockVectorX, lambda] = lobpcg (randn (n, 2), A, blockVectorY, 1e-5, 200, 2); blockVectorY = [blockVectorY, blockVectorX]; lambda_all = [lambda_all' lambda']'; pause; endattemps to compute the same 8 eigenpairs by calling the code 4 times with blockSize=2 using orthogonalization to the previously founded eigenvectors.

The following Example:

R = ichol (A, struct('michol', 'on')); precfun = @(x)R\(R'\x); [blockVectorX, lambda, failureFlag] = lobpcg (randn (n, 8), A, [], @(x)precfun(x), 1e-5, 60, 2);computes the same eigenpairs in less then 25 steps, so that failureFlag=0 using the preconditioner function

`precfun`

, defined inline. If`precfun`

is defined as an octave function in a file, the function handle`@(x)precfun(x)`

can be equivalently replaced by the function name`precfun`

. Running:[blockVectorX, lambda, failureFlag] = lobpcg (randn (n, 8), A, speye (n), @(x)precfun(x), 1e-5, 50, 2);produces similar answers, but is somewhat slower and needs more memory as technically a generalized eigenproblem with B=I is solved here.

The following example for a mostly diagonally dominant sparse matrix A demonstrates different types of preconditioning, compared to the standard use of the main diagonal of A:

clear all; close all; n = 1000; M = spdiags ([1:n]', 0, n, n); precfun = @(x)M\x; A = M + sprandsym (n, .1); Xini = randn (n, 5); maxiter = 15; tol = 1e-5; [~,~,~,~,rnp] = lobpcg (Xini, A, tol, maxiter, 1); [~,~,~,~,r] = lobpcg (Xini, A, [], @(x)precfun(x), tol, maxiter, 1); subplot (2,2,1), semilogy (r'); hold on; semilogy (rnp', ':>'); title ('No preconditioning (top)'); axis tight; M(1,2) = 2; precfun = @(x)M\x; % M is no longer symmetric [~,~,~,~,rns] = lobpcg (Xini, A, [], @(x)precfun(x), tol, maxiter, 1); subplot (2,2,2), semilogy (r'); hold on; semilogy (rns', '--s'); title ('Nonsymmetric preconditioning (square)'); axis tight; M(1,2) = 0; precfun = @(x)M\(x+10*sin(x)); % nonlinear preconditioning [~,~,~,~,rnl] = lobpcg (Xini, A, [], @(x)precfun(x), tol, maxiter, 1); subplot (2,2,3), semilogy (r'); hold on; semilogy (rnl', '-.*'); title ('Nonlinear preconditioning (star)'); axis tight; M = abs (M - 3.5 * speye (n, n)); precfun = @(x)M\x; [~,~,~,~,rs] = lobpcg (Xini, A, [], @(x)precfun(x), tol, maxiter, 1); subplot (2,2,4), semilogy (r'); hold on; semilogy (rs', '-d'); title ('Selective preconditioning (diamond)'); axis tight;## References

This main function

`lobpcg`

is a version of the preconditioned conjugate gradient method (Algorithm 5.1) described in A. V. Knyazev, Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient Method, SIAM Journal on Scientific Computing 23 (2001), no. 2, pp. 517-541. http://dx.doi.org/10.1137/S1064827500366124## Known bugs/features

- an excessively small requested tolerance may result in often restarts and instability. The code is not written to produce an eps-level accuracy! Use common sense.
- the code may be very sensitive to the number of eigenpairs computed, if there is a cluster of eigenvalues not completely included, cf.
operatorA = diag ([1 1.99 2:99]); [blockVectorX, lambda] = lobpcg (randn (100, 1),operatorA, 1e-10, 80, 2); [blockVectorX, lambda] = lobpcg (randn (100, 2),operatorA, 1e-10, 80, 2); [blockVectorX, lambda] = lobpcg (randn (100, 3),operatorA, 1e-10, 80, 2);## Distribution

The main distribution site: http://math.ucdenver.edu/~aknyazev/

A C-version of this code is a part of the http://code.google.com/p/blopex/ package and is directly available, e.g., in PETSc and HYPRE.

Package: linear-algebra