JDQR

From this page you can get a Matlab® implementation of the JDQR algorithm.
The JDQR algorithm can be used for computing a few selected eigenvalues with some desirable property together with the associated eigenvectors of a matrix A. The matrix can be real or complex, Hermitian or non-Hermitian, .... The algorithm is effective especially in case A is sparse and of large size.
The Jacobi-Davidson method is used to compute a partial Schur decomposition of A. The 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 algorithm, 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® implementation here is based on the algorithms as discussed in
Zhaojun Bai, James Demmel, Jack Dongarra, Axel Ruhe, and Henk Van der Vorst (eds.),
Templates for the Solution of Algebraic Eigenvalue Problems: a Practical Guide (project in Progress),
Chapter 3.6 (for symmetric matrices A) and
Chapter 5.11 (for general type of matrices A).
The Matlab® code here contains a number of function M-files (as a GMRES, Bi-CGSTAB, BiCGstab(ell), ... implementation for solving linear systems) that can be used the improve the performance of the JDQZ algorithm. These functions M-files are bundled in one file jdqr.m. In this form, 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 non symmetric, 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 preconditioner 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.1.
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.1.
Function: testA
Description: Contains some examples for using preconditioning in jdqr.m

File: Example1.m
Requires: Matlab Version 5.1.
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.1.
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 preconditioner". 

File: ILU.m
Requires: Matlab Version 5.1.
Function: out=ILU(v)
Description: Example of a function file for jdqr of a preconditioner. out is the result of the preconditioner applied to the vector v.

Gerard L. G. Sleijpen© G.L.G.Sleijpen@uu.nl

Last update: April 8 1999