JDQR

(This document in ps,pdf)

jdqr computes eigenpairs of a square matrix or operator.

Lambda = jdqr(A) returns the absolute largest eigenvalues of A in a k vector Lambda. Here k = min(5,n) and n = size(A,1). jdqr(A) (without output argument) displays the k eigenvalues.

[X,Jordan] = jdqr(A,B) returns the eigenvectors X and the Jordan structure Jordan: A$\ast$X = X$\ast$Jordan. The diagonal of Jordan contains the eigenvalues: Lambda = diag(Jordan). Jordan is an k by k matrix with the eigenvalues on the diagonal and possibly non-zeros on the first upper diagonal elements. The other entries are zero. The columns of X have norm 1.

[X,Jordan,Q,S] = jdqr(A)
If four or more output arguments are required then Q is n by k orthonormal, and S is k by k upper triangular such that they form a partial generalized Schur decomposition: A$\ast$Q = Q$\ast$S. Then Lambda = diag(S) and X = Q$\ast$Y with Y the eigenvectors of the pair S: S$\ast$Y = Y$\ast$Jordan (see also Options.Schur).

... = jdqr(A)
... = jdqr('Afun')
... = jdqr('Afun',n)
The first input argument is either a square matrix (which can be full or sparse, symmetric or nonsymmetric, real or complex), or a string containing the name of an M-file which applies a linear operator to the columns of a given matrix. In the latter case, the M-file, say Afun.m, must return the dimensionn of the problem with n = Afun([ ],'dimension') orn must be specified in the list of input arguments. For example, jdqr('fft',...) is much faster than jdqr(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,M)
... = jdqr('Afun',k,Sigma,Options,M),

where

	k	an integer, the number of desired eigenvalues.
	Sigma	a scalar shift or a two letter string.
	Options	a structure containing additional parameters.
	M	a string or a matrix that specifies the preconditioner.

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 a real or complex scalar, then the kth eigenvalues nearest Sigma are computed. If Sigma is row vector of size (1,m), then the jth eigenvalue nearest to Sigma(1,min(m,j)) is computed for j = 1:k. Sigma is the ``target'' for the desired eigenvalues. If Sigma is one of the following strings, then it specifies the desired eigenvalues.

	Sigma	Specified eigenvalues
	'LM'	Largest Magnitude
	'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.)

If 'TestSpace' is 'Harmonic' (see Options), then Sigma = 0 is the default, otherwise Sigma = 'LM' is the default.

The Options structure specifies certain parameters in the algorithm.

	Field name	Parameter	Default
	Options.Tol	Convergence tolerance: norm(r) $<=$ tol/sqrt(k)	1e-8
	Options.jmin	Minimum dimension search subspace V	k+5
	Options.jmax	Maximum dimension search subspace V	jmin+5
	Options.MaxIt	Maximum number of iterations.	100
	Options.v0	Starting space	ones+0.1$\ast$rand
	Options.Schur	Gives schur decomposition If 'yes', then X and Jordan are not computed and [Q,S,history] is the list of output arguments.	'no'
	Options.TestSpace	Defines the test subspace W 'Standard': W = V 'Harmonic': W = A$\ast$V-sigma$\ast$V	'standard'
	Options.Disp	If Disp=1 or Disp='yes', then, the input parameters, the residual size at each step, and the eigenvalues at detection are displayed, and the convergence history is plotted. If Disp=2 then, in addition, approximate eigenvalues are plotted at each step.	'no'
	Options.LSolver	Linear solver	'GMRES'
	Options.LS$\underline{~}$Tol	Residual reduction linear solver	1,0.7,0.7 $\,\widehat{~}\,$2,...
	Options.LS$\underline{~}$MaxIt	Maximum number it. linear solver	5
	Options.LS$\underline{~}$ell	ell for BiCGstab(ell)	4
	Options.Precond	Preconditioner.	identity.

For instance,

Options = struct('Tol',1.0e-8,'LSolver','BiCGstab','LS$\underline{~}$ell',4,'Precond',M);

changes the convergence tolerance to 1.0e-8, takes BiCGstab as linear solver, and takes M as preconditioner (for ways of defining M, see below).

There are a few other Options that can be specified. They are listed here.

[X,Jordan,history] = jdqr(A,...)
[X,Jordan,Q,S,history] = jdqr(A,...)
returns also the convergence history.
history is an array of 3 columns: history(i,1) is the residual norm at step j = history(i,2), history(i,3) is the cumulative number of multiplications by A at step j. If a search for a new eigenvalue is started at step J then j = history(i,2) = history(i+1,2), history(i,1) is the norm of the "old" residual, and history(i+1,1) is the norm of the "new" one. history is empty if the required number of eigenvalues are detected in the initialization phase.

jdqr without input arguments returns the options and its defaults.

Subsections

• Preconditioning in jdqr
• Way of using a preconditioner
• Additional Options
• Accuracy