JDQR

From this page you can obtain a Matlab® implementation of the JDQR algorithm.

The JDQR algorithm can be used for computing a few selected eigenvalues, with the associated eigenvectors, of a matrix A. The matrix may be real or complex, Hermitian or non-Hermitian, respectively. The algorithm is effective especially when A is sparse and of large size.
The Jacobi-Davidson method is used to compute a partial Schur decomposition of A. This decomposition leads to the wanted eigenpairs.

The Jacobi-Davidson has been introduced in

G.L.G. Sleijpen and H.A. van der Vorst,
A Jacobi-Davidson iteration method for linear eigenvalue problems,
SIAM J. Matrix Anal. Appl. (SIMAX), 17 (1996), pp. 401-425 .

For the JDQR algorithme, see

D.R. Fokkema, G.L.G. Sleijpen, and H.A. van der Vorst,
Jacobi-Davidson style QR and QZ algorithms for the reduction of matrix pencils,
SIAM J. Sc. Comput., 20:1 (1998), pp. 94-125

The Matlab® code contains a number of function M-files (as a GMRES, Bi-CGSTAB, BiCGstab(ell), ... implementations for solving linear systems) that can be used t0 improve the performance of the JDQR algorithm. These function M-files are bundled in one file jdqr.m. For these functions, jdqr.m requires Matlab version 5.1 or higher.

We also provide a few "test" files that can help to see how the jdqr.m can be used. jdqr.m.tar.gz contains the M-files jdqr.m plus the test files in tarred and zipped form (apply gunzip jdqr.m.tar.gz and tar xvf jdqr.m.tar to unpack). jdqr.m and the test files can also be obtained separately: see below, where you can also find a description of their use.


File: jdqr.m
Requires: Matlab Version 5.1.
Function: [X,Lambda,Q,S,HISTORY] = JDQR(matrix/filename,K,SIGMA,OPTIONS);

Description: Lambda = JDQR(A)
returns the K absolute largest eigenvalues in the K vector Lambda. Here K=MIN(N,5) (unless K has been specified) and N=size(A,1).
JDQR(A) (without output argument) displays the eigenvalues on the screen.

[X,Lambda] = JDQR(A)

returns the eigenvectors in the N by K matrix X and the eigenvalues in the K by K diagonal matrix Lambda. Lambda contains the Jordan structure if there are multiple eigenvalues.

[X,Lambda,HISTORY] = JDQR(A)

returns also the convergence history (norms of the subsequential residuals).

[X,Lambda,Q,S] = JDQR(A)

returns also a partial Schur decomposition: S is a K by K upper triangular matrix and Q is an N by K orthonormal matrix such that A*Q = Q*S. The diagonal elements of S are eigenvalues of A.

[X,Lambda,Q,S,HISTORY] = JDQR(A)

returns also the convergence history.

... = JDQR('Afun')

... = JDQR('Afun',N)
The first input argument is either a square matrix A (which can be full or sparse, symmetric or nonsymmetric, real or complex), or a string containing the name of an M-file ('Afun') which applies a linear operator to a given column vector. In the latter case, the M-file must return the the order N of the problem with N = Afun([],'dimension') or N must be an input argument. For example, EIGS('fft',...) is much faster than EIGS(F,...) where F is the explicit FFT matrix.

The remaining input arguments are optional and can be given in practically any order:

... = JDQR(A,K,SIGMA,OPTIONS)
... = JDQR('Afun',K,SIGMA,OPTIONS)
where
K An integer, the number of eigenvalues desired.
SIGMA A scalar shift or a two letter string.
OPTIONS A structure containing additional parameters.

If K is not specified, then K = MIN(N,5) eigenvalues are computed.

If SIGMA is not specified, then the Kth eigenvalues largest in magnitude are computed. If SIGMA is zero, then the Kth eigenvalues smallest in magnitude are computed. If SIGMA is a real or complex scalar, then the Kth eigenvalues nearest SIGMA are computed. If SIGMA is one of the following strings, then it specifies the desired eigenvalues.

SIGMA Location wanted eigenvalues
'LM' Largest Magnitude (the default)
'SM'Smallest Magnitude (same as SIGMA = 0)
'LR'Largest Real part
'SR'Smallest Real part
'BE'Both Ends. Computes K/2 eigenvalues from each end of the spectrum (one more from the high end if K is odd.)

The OPTIONS structure specifies certain parameters in the algorithm.

Field name Parameter Default
OPTIONS.Tol Convergence tolerance: norm(A*Q-Q*S,1) <= tol * norm(A,1) 1e-8
OPTIONS.jmin minimum dimension search subspace k+5
OPTIONS.jmax maximum dimension search subspace jmin+5
OPTIONS.MaxIt Maximum number of iterations. 100
OPTIONS.v0 Starting space ones(N,1)+0.1*rand(N,1)
OPTIONS.Schur Gives schur decomposition also in case of 2 or 3 output arguments (X=Q, Lambda=R). 'no'
OPTIONS.TestSpace For using harmonic Ritz values If 'TestSpace'='Harmonic' then SIGMA=0 is the default value for SIGMA 'Standard'
OPTIONS.Disp Shows size of intermediate residuals and displays the approximate eigenvalues. 0
OPTIONS.LSolver Linear solver 'GMRES'
OPTIONS.LS_Tol Residual reduction linear solver 1,0.7,0.7^2,..
OPTIONS.LS_MaxIt Maximum number iteration steps of the linear solver 5
OPTIONS.LS_ell ell for BiCGstab(ell) 4
OPTIONS.Precond Preconditioner LU=[[],[]]

For instance

OPTIONS= STRUCT('Tol',1.0e-10,'LSolver','BiCGstab','LS_ell',2,'Precond',M);
changes the convergence tolerance to 1.0e-10, takes BiCGstab(2) as linear solver and the preconditioner defined in ILU.m if M is the string 'ILU', or M = L*U if M is an N by 2*N matrix: M = [L,U].

The preconditoner can be specified in the OPTIONS structure, but also in the argument list:

... = JDQR(...,M,...)
... = JDQR(...,L,U,..)
... = JDQR(...,'M',...)
... = JDQR(...,'L','U',...)
as an N by N matrix M (then M is the preconditioner), or an N by 2*N matrix M (then L*U is the preconditioner, where M = [L,U]), or as N by N matrices L and U (then L*U is the preconditioner), or as one or two strings containing the name of M-files ('M', or 'L' and 'U') which apply a linear operator to a given column vector.

JDQR (without input arguments) lists the options and the defaults.


File: testB.m
Requires: jdqr.m and Matlab Version 5.2.
Function:testB

Description: Contains some examples for using jdqr.m


File: testA.m
Requires: jdqr.m, ILU.m, Example1.m, Example2.m and Matlab Version 5.2.
Function:testA

Description: Contains some examples for using preconditioning in jdqr.m


File: Example1.m
Requires: Matlab Version 5.2.
Function:out=Example1(v)

Description: Example of a function file for jdqr of a linear operator. out is the result of the operator applied to the vector v.


File: Example2.m
Requires: Matlab Version 5.2.
Function:out=Example2(v,flag)

Description: Example of a function file for jdqr of a linear operator. out is the result of the operator applied to the vector v when applied with 1 input argument. If flag='dimension' then out is the dimension N. If flag='preconditioner' then out is the result of "a preconditoner".


File: ILU.m
Requires: Matlab Version 5.2.
Function:out=ILU(v)

Description: Example of a function file for jdqr of a preconditioner. out is the result of the preconditoner applied to the vector v.


Henk A. van der Vorst© vorst@math.uu.nl

Last update: April 19 1999