Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_eigen_real_symm_sparse_eigsys (f02fj)

Purpose

nag_eigen_real_symm_sparse_eigsys (f02fj) finds eigenvalues and eigenvectors of a real sparse symmetric or generalized symmetric eigenvalue problem.

Syntax

[m, noits, x, d, user, ifail] = f02fj(m, noits, tol, dot, image, monit, novecs, x, 'n', n, 'k', k, 'user', user)
[m, noits, x, d, user, ifail] = nag_eigen_real_symm_sparse_eigsys(m, noits, tol, dot, image, monit, novecs, x, 'n', n, 'k', k, 'user', user)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 22: n has been made optional
.

Description

nag_eigen_real_symm_sparse_eigsys (f02fj) finds the m$m$ eigenvalues of largest absolute value and the corresponding eigenvectors for the real eigenvalue problem
 Cx = λ x $Cx=λ x$ (1)
where C$C$ is an n$n$ by n$n$ matrix such that
 BC = CTB $BC=CTB$ (2)
for a given positive definite matrix B$B$. C$C$ is said to be B$B$-symmetric. Different specifications of C$C$ allow for the solution of a variety of eigenvalue problems. For example, when
 C = A  and  B = I  where  A = AT $C=A and B=I where A=AT$
the function finds the m$m$ eigenvalues of largest absolute magnitude for the standard symmetric eigenvalue problem
 Ax = λx. $Ax=λx.$ (3)
The function is intended for the case where A$A$ is sparse.
As a second example, when
 C = B − 1A $C=B-1A$
where
 A = AT $A=AT$
the function finds the m$m$ eigenvalues of largest absolute magnitude for the generalized symmetric eigenvalue problem
 Ax = λBx. $Ax=λBx.$ (4)
The function is intended for the case where A$A$ and B$B$ are sparse.
The function does not require C$C$ explicitly, but C$C$ is specified via image which, given an n$n$-element vector z$z$, computes the image w$w$ given by
 w = Cz. $w=Cz.$
For instance, in the above example, where C = B1A$C={B}^{-1}A$, image will need to solve the positive definite system of equations Bw = Az$Bw=Az$ for w$w$.
To find the m$m$ eigenvalues of smallest absolute magnitude of (3) we can choose C = A1$C={A}^{-1}$ and hence find the reciprocals of the required eigenvalues, so that image will need to solve Aw = z$Aw=z$ for w$w$, and correspondingly for (4) we can choose C = A1B$C={A}^{-1}B$ and solve Aw = Bz$Aw=Bz$ for w$w$.
A table of examples of choice of image is given in Table 1. It should be remembered that the function also returns the corresponding eigenvectors and that B$B$ is positive definite. Throughout A$A$ is assumed to be symmetric and, where necessary, nonsingularity is also assumed.
 EigenvaluesRequired Problem Ax = λx (B = I)$Ax=\lambda x\left(B=I\right)$ Ax = λBx$Ax=\lambda Bx$ ABx = λx$ABx=\lambda x$ Largest Compute w = Az$w=Az$ Solve Bw = Az$Bw=Az$ Compute w = ABz$w=ABz$ Smallest (Find 1 / λ$1/\lambda$) Solve Aw = z$Aw=z$ Solve Aw = Bz$Aw=Bz$ Solve Av = z$Av=z$, Bw = v$Bw=v$ Furthest from σ$\sigma$  (Find λ − σ$\lambda -\sigma$) Compute w = (A − σI)z$w=\left(A-\sigma I\right)z$ Solve Bw = (A − σB)z$Bw=\left(A-\sigma B\right)z$ Compute w = (AB − σI)z$w=\left(AB-\sigma I\right)z$ Closest to σ$\sigma$  (Find 1 / (λ − σ)$1/\left(\lambda -\sigma \right)$) Solve (A − σI)w = z$\left(A-\sigma I\right)w=z$ Solve (A − σB)w = Bz$\left(A-\sigma B\right)w=Bz$ Solve (AB − σI)w = z$\left(AB-\sigma I\right)w=z$
Table 1
The Requirement of image for Various Problems.
The matrix B$B$ also need not be supplied explicitly, but is specified via dot which, given n$n$-element vectors z$z$ and w$w$, computes the generalized dot product wTBz${w}^{\mathrm{T}}Bz$.
nag_eigen_real_symm_sparse_eigsys (f02fj) is based upon routine SIMITZ (see Nikolai (1979)), which is itself a derivative of the Algol procedure ritzit (see Rutishauser (1970)), and uses the method of simultaneous (subspace) iteration. (See Parlett (1998) for a description, analysis and advice on the use of the method.)
The function performs simultaneous iteration on k > m$k>m$ vectors. Initial estimates to pk$p\le k$ eigenvectors, corresponding to the p$p$ eigenvalues of C$C$ of largest absolute value, may be supplied to nag_eigen_real_symm_sparse_eigsys (f02fj). When possible k$k$ should be chosen so that the k$k$th eigenvalue is not too close to the m$m$ required eigenvalues, but if k$k$ is initially chosen too small then nag_eigen_real_symm_sparse_eigsys (f02fj) may be re-entered, supplying approximations to the k$k$ eigenvectors found so far and with k$k$ then increased.
At each major iteration nag_eigen_real_symm_sparse_eigsys (f02fj) solves an r$r$ by r$r$ (rk$r\le k$) eigenvalue sub-problem in order to obtain an approximation to the eigenvalues for which convergence has not yet occurred. This approximation is refined by Chebyshev acceleration.

References

Nikolai P J (1979) Algorithm 538: Eigenvectors and eigenvalues of real generalized symmetric matrices by simultaneous iteration ACM Trans. Math. Software 5 118–125
Parlett B N (1998) The Symmetric Eigenvalue Problem SIAM, Philadelphia
Rutishauser H (1969) Computational aspects of F L Bauer's simultaneous iteration method Numer. Math. 13 4–13
Rutishauser H (1970) Simultaneous iteration method for symmetric matrices Numer. Math. 16 205–223

Parameters

Compulsory Input Parameters

1:     m – int64int32nag_int scalar
m$m$, the number of eigenvalues required.
Constraint: m1${\mathbf{m}}\ge 1$.
2:     noits – int64int32nag_int scalar
The maximum number of major iterations (eigenvalue sub-problems) to be performed. If noits0${\mathbf{noits}}\le 0$, the value 100$100$ is used in place of noits.
3:     tol – double scalar
A relative tolerance to be used in accepting eigenvalues and eigenvectors. If the eigenvalues are required to about t$t$ significant figures, tol should be set to about 10t${10}^{-t}$. di${d}_{i}$ is accepted as an eigenvalue as soon as two successive approximations to di${d}_{i}$ differ by less than (|i| × tol) / 10$\left(|{\stackrel{~}{d}}_{i}|×{\mathbf{tol}}\right)/10$, where i${\stackrel{~}{d}}_{i}$ is the latest approximation to di${d}_{i}$. Once an eigenvalue has been accepted, an eigenvector is accepted as soon as (difi) / (didk) < tol$\left({d}_{i}{f}_{i}\right)/\left({d}_{i}-{d}_{k}\right)<{\mathbf{tol}}$, where fi${f}_{i}$ is the normalized residual of the current approximation to the eigenvector (see Section [Further Comments] for further information). The values of the fi${f}_{i}$ and di${d}_{i}$ can be printed from monit. If tol is supplied outside the range (ε,1.0$\epsilon ,1.0$), where ε$\epsilon$ is the machine precision, the value ε$\epsilon$ is used in place of tol.
4:     dot – function handle or string containing name of m-file
dot must return the value wTBz${w}^{\mathrm{T}}Bz$ for given vectors w$w$ and z$z$. For the standard eigenvalue problem, where B = I$B=I$, dot must return the dot product wTz${w}^{\mathrm{T}}z$.
[result, iflag, user] = dot(iflag, n, z, w, user)

Input Parameters

1:     iflag – int64int32nag_int scalar
Is always non-negative.
2:     n – int64int32nag_int scalar
The number of elements in the vectors z$z$ and w$w$ and the order of the matrix B$B$.
3:     z(n) – double array
The vector z$z$ for which wTBz${w}^{\mathrm{T}}Bz$ is required.
4:     w(n) – double array
The vector w$w$ for which wTBz${w}^{\mathrm{T}}Bz$ is required.
5:     user – Any MATLAB object
dot is called from nag_eigen_real_symm_sparse_eigsys (f02fj) with the object supplied to nag_eigen_real_symm_sparse_eigsys (f02fj).

Output Parameters

1:     result – double scalar
The result of the function.
2:     iflag – int64int32nag_int scalar
May be used as a flag to indicate a failure in the computation of wTBz${w}^{\mathrm{T}}Bz$. If iflag is negative on exit from dot, nag_eigen_real_symm_sparse_eigsys (f02fj) will exit immediately with ifail set to iflag. Note that in this case dot must still be assigned a value.
3:     user – Any MATLAB object
5:     image – function handle or string containing name of m-file
image must return the vector w = Cz$w=Cz$ for a given vector z$z$.
[iflag, w, user] = image(iflag, n, z, user)

Input Parameters

1:     iflag – int64int32nag_int scalar
Is always non-negative.
2:     n – int64int32nag_int scalar
n$n$, the number of elements in the vectors w$w$ and z$z$, and the order of the matrix C$C$.
3:     z(n) – double array
The vector z$z$ for which Cz$Cz$ is required.
4:     user – Any MATLAB object
image is called from nag_eigen_real_symm_sparse_eigsys (f02fj) with the object supplied to nag_eigen_real_symm_sparse_eigsys (f02fj).

Output Parameters

1:     iflag – int64int32nag_int scalar
May be used as a flag to indicate a failure in the computation of w$w$. If iflag is negative on exit from image, nag_eigen_real_symm_sparse_eigsys (f02fj) will exit immediately with ifail set to iflag.
2:     w(n) – double array
The vector w = Cz$w=Cz$.
3:     user – Any MATLAB object
6:     monit – function handle or string containing name of m-file
monit is used to monitor the progress of nag_eigen_real_symm_sparse_eigsys (f02fj). monit may be the dummy function nag_eigen_real_symm_sparse_eigsys_dummy_monit (f02fjz) if no monitoring is actually required. (nag_eigen_real_symm_sparse_eigsys_dummy_monit (f02fjz) is included in the NAG Toolbox.) monit is called after the solution of each eigenvalue sub-problem and also just prior to return from nag_eigen_real_symm_sparse_eigsys (f02fj). The parameters istate and nextit allow selective printing by monit.
monit(istate, nextit, nevals, nevecs, k, f, d)

Input Parameters

1:     istate – int64int32nag_int scalar
Specifies the state of nag_eigen_real_symm_sparse_eigsys (f02fj).
istate = 0${\mathbf{istate}}=0$
No eigenvalue or eigenvector has just been accepted.
istate = 1${\mathbf{istate}}=1$
One or more eigenvalues have been accepted since the last call to monit.
istate = 2${\mathbf{istate}}=2$
One or more eigenvectors have been accepted since the last call to monit.
istate = 3${\mathbf{istate}}=3$
One or more eigenvalues and eigenvectors have been accepted since the last call to monit.
istate = 4${\mathbf{istate}}=4$
Return from nag_eigen_real_symm_sparse_eigsys (f02fj) is about to occur.
2:     nextit – int64int32nag_int scalar
The number of the next iteration.
3:     nevals – int64int32nag_int scalar
The number of eigenvalues accepted so far.
4:     nevecs – int64int32nag_int scalar
The number of eigenvectors accepted so far.
5:     k – int64int32nag_int scalar
k$k$, the number of simultaneous iteration vectors.
6:     f(k) – double array
A vector of error quantities measuring the state of convergence of the simultaneous iteration vectors. See tol and Section [Further Comments] for further details. Each element of f is initially set to the value 4.0$4.0$ and an element remains at 4.0$4.0$ until the corresponding vector is tested.
7:     d(k) – double array
d(i)${\mathbf{d}}\left(i\right)$ contains the latest approximation to the absolute value of the i$i$th eigenvalue of C$C$.
7:     novecs – int64int32nag_int scalar
The number of approximate vectors that are being supplied in x. If novecs is outside the range (0,k)$\left(0,{\mathbf{k}}\right)$, the value 0$0$ is used in place of novecs.
8:     x(ldx,k) – double array
ldx, the first dimension of the array, must satisfy the constraint ldxn$\mathit{ldx}\ge {\mathbf{n}}$.
If 0 < novecsk$0<{\mathbf{novecs}}\le {\mathbf{k}}$, the first novecs columns of x must contain approximations to the eigenvectors corresponding to the novecs eigenvalues of largest absolute value of C$C$. Supplying approximate eigenvectors can be useful when reasonable approximations are known, or when nag_eigen_real_symm_sparse_eigsys (f02fj) is being restarted with a larger value of k. Otherwise it is not necessary to supply approximate vectors, as simultaneous iteration vectors will be generated randomly by nag_eigen_real_symm_sparse_eigsys (f02fj).

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the array x.
n$n$, the order of the matrix C$C$.
Constraint: n1${\mathbf{n}}\ge 1$.
2:     k – int64int32nag_int scalar
Default: The second dimension of the array x.
The number of simultaneous iteration vectors to be used. Too small a value of k may inhibit convergence, while a larger value of k incurs additional storage and additional work per iteration.
Default: m + 4${\mathbf{m}}+4$
Constraint: m < kn${\mathbf{m}}<{\mathbf{k}}\le {\mathbf{n}}$.
3:     user – Any MATLAB object
user is not used by nag_eigen_real_symm_sparse_eigsys (f02fj), but is passed to dot and image. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

Input Parameters Omitted from the MATLAB Interface

ldx work lwork ruser lruser iuser liuser

Output Parameters

1:     m – int64int32nag_int scalar
m${m}^{\prime }$, the number of eigenvalues actually found. It is equal to m$m$ if ${\mathbf{ifail}}={\mathbf{0}}$ on exit, and is less than m$m$ if ${\mathbf{ifail}}={\mathbf{2}}$, 3${\mathbf{3}}$ or 4${\mathbf{4}}$. See Sections [Error Indicators and Warnings] and [Further Comments] for further information.
2:     noits – int64int32nag_int scalar
The number of iterations actually performed.
3:     x(ldx,k) – double array
ldxn$\mathit{ldx}\ge {\mathbf{n}}$.
If ${\mathbf{ifail}}={\mathbf{0}}$, 2${\mathbf{2}}$, 3${\mathbf{3}}$ or 4${\mathbf{4}}$, the first m${m}^{\prime }$ columns contain the eigenvectors corresponding to the eigenvalues returned in the first m${m}^{\prime }$ elements of d; and the next km1$k-{m}^{\prime }-1$ columns contain approximations to the eigenvectors corresponding to the approximate eigenvalues returned in the next km1$k-{m}^{\prime }-1$ elements of d. Here m${m}^{\prime }$ is the value returned in m, the number of eigenvalues actually found. The k$k$th column is used as workspace.
4:     d(k) – double array
If ${\mathbf{ifail}}={\mathbf{0}}$, 2${\mathbf{2}}$, 3${\mathbf{3}}$ or 4${\mathbf{4}}$, the first m${m}^{\prime }$ elements contain the first m${m}^{\prime }$ eigenvalues in decreasing order of magnitude; and the next km1$k-{m}^{\prime }-1$ elements contain approximations to the next km1$k-{m}^{\prime }-1$ eigenvalues. Here m${m}^{\prime }$ is the value returned in m, the number of eigenvalues actually found. d(k)${\mathbf{d}}\left(k\right)$ contains the value e$e$ where (e,e)$\left(-e,e\right)$ is the latest interval over which Chebyshev acceleration is performed.
5:     user – Any MATLAB object
6:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

W ifail < 0${\mathbf{ifail}}<0$
A negative value of ifail indicates an exit from nag_eigen_real_symm_sparse_eigsys (f02fj) because you have set iflag negative in dot or image. The value of ifail will be the same as your setting of iflag.
ifail = 1${\mathbf{ifail}}=1$
 On entry, n < 1${\mathbf{n}}<1$, or m < 1${\mathbf{m}}<1$, or m ≥ k${\mathbf{m}}\ge {\mathbf{k}}$, or k > n${\mathbf{k}}>{\mathbf{n}}$, or ldx < n$\mathit{ldx}<{\mathbf{n}}$, or lwork < 3 × k + max (k × k,2 × n)$\mathit{lwork}<3×{\mathbf{k}}+\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{k}}×{\mathbf{k}},2×{\mathbf{n}}\right)$, or lruser < 1$\mathit{lruser}<1$, or liuser < 1$\mathit{liuser}<1$.
W ifail = 2${\mathbf{ifail}}=2$
Not all the requested eigenvalues and vectors have been obtained. Approximations to the r$r$th eigenvalue are oscillating rapidly indicating that severe cancellation is occurring in the r$r$th eigenvector and so m is returned as (r1)$\left(r-1\right)$. A restart with a larger value of k may permit convergence.
W ifail = 3${\mathbf{ifail}}=3$
Not all the requested eigenvalues and vectors have been obtained. The rate of convergence of the remaining eigenvectors suggests that more than noits iterations would be required and so the input value of m has been reduced. A restart with a larger value of k may permit convergence.
W ifail = 4${\mathbf{ifail}}=4$
Not all the requested eigenvalues and vectors have been obtained. noits iterations have been performed. A restart, possibly with a larger value of k, may permit convergence.
ifail = 5${\mathbf{ifail}}=5$
This error is very unlikely to occur, but indicates that convergence of the eigenvalue sub-problem has not taken place. Restarting with a different set of approximate vectors may allow convergence. If this error occurs you should check carefully that nag_eigen_real_symm_sparse_eigsys (f02fj) is being called correctly.

Accuracy

Eigenvalues and eigenvectors will normally be computed to the accuracy requested by the parameter tol, but eigenvectors corresponding to small or to close eigenvalues may not always be computed to the accuracy requested by the parameter tol. Use of the monit to monitor acceptance of eigenvalues and eigenvectors is recommended.

The time taken by nag_eigen_real_symm_sparse_eigsys (f02fj) will be principally determined by the time taken to solve the eigenvalue sub-problem and the time taken by dot and image. The time taken to solve an eigenvalue sub-problem is approximately proportional to nk2$n{k}^{2}$. It is important to be aware that several calls to dot and image may occur on each major iteration.
As can be seen from Table 1, many applications of nag_eigen_real_symm_sparse_eigsys (f02fj) will require the image to solve a system of linear equations. For example, to find the smallest eigenvalues of Ax = λ Bx$Ax=\lambda Bx$, image needs to solve equations of the form Aw = Bz$Aw=Bz$ for w$w$ and functions from Chapters F01 and F04 will frequently be useful in this context. In particular, if A$A$ is a positive definite variable band matrix, nag_linsys_real_posdef_vband_solve (f04mc) may be used after A$A$ has been factorized by nag_matop_real_vband_posdef_fac (f01mc). Thus factorization need be performed only once prior to calling nag_eigen_real_symm_sparse_eigsys (f02fj). An illustration of this type of use is given in the example program.
An approximation h${\stackrel{~}{d}}_{h}$, to the i$i$th eigenvalue, is accepted as soon as h${\stackrel{~}{d}}_{h}$ and the previous approximation differ by less than |h| × tol / 10$|{\stackrel{~}{d}}_{h}|×{\mathbf{tol}}/10$. Eigenvectors are accepted in groups corresponding to clusters of eigenvalues that are equal, or nearly equal, in absolute value and that have already been accepted. If dr${d}_{r}$ is the last eigenvalue in such a group and we define the residual rj${r}_{j}$ as
 rj = Cxj − yr $rj=Cxj-yr$
where yr${y}_{r}$ is the projection of Cxj$C{x}_{j}$, with respect to B$B$, onto the space spanned by x1,x2,, xr${x}_{1},{x}_{2},\dots ,{x}_{r}$, and xj${x}_{j}$ is the current approximation to the j$j$th eigenvector, then the value fi${f}_{i}$ returned in monit is given by
 fi = max‖rj‖B / ‖Cxj‖B   ‖x‖B2 = xTBx $fi = max‖rj‖B / ‖Cxj‖B ‖x‖B2 = xTBx$
and each vector in the group is accepted as an eigenvector if
 (|dr|fr) / (|dr| − e) < tol, $(|dr| fr)/(|dr|-e)
where e$e$ is the current approximation to |k|$|{\stackrel{~}{d}}_{k}|$. The values of the fi${f}_{i}$ are systematically increased if the convergence criteria appear to be too strict. See Rutishauser (1970) for further details.
The algorithm implemented by nag_eigen_real_symm_sparse_eigsys (f02fj) differs slightly from SIMITZ (see Nikolai (1979)) in that the eigenvalue sub-problem is solved using the singular value decomposition of the upper triangular matrix R$R$ of the Gram–Schmidt factorization of Cxr$C{x}_{r}$, rather than forming RTR${R}^{\mathrm{T}}R$.

Example

```function nag_eigen_real_symm_sparse_eigsys_example
n = 16;
m = 4;
noits = int64(1000);
tol = 0.0001;
novecs = int64(0);
x = zeros(n, m+2);
% a, b will be passed in user
a = diag(ones(n,1))+diag(-0.25*ones(n-1,1),1)+diag(-0.25*ones(n-1,1),-1)+...
diag(-0.25*ones(n-4,1),4)+diag(-0.25*ones(n-4,1),-4);
b = diag(ones(n,1))+diag(-0.5*ones(n-1,1),1)+diag(-0.5*ones(n-1,1),-1);
[mOut, noitsOut, xOut, d, user, ifail] = ...
nag_eigen_real_symm_sparse_eigsys(int64(m), noits, tol, ...
@dot, @image, @monit, novecs, x, 'user', {a, b})

function [result, iflag, user] = dot(iflag, n, z, w, user)
b = user{2};
result=transpose(w)*b*z;
function [iflag, w, user] = image(iflag, n, z, user)

a=user{1};
b=user{2};

w=inv(a)*b*z;
function monit(istate, nextit, nevals, nevecs, k, f, d)

if (istate ~= 0)
fprintf('\n  istate = %d nextit = %d\n', istate, nextit);
fprintf('  nevals = %d nevecs = %d\n', nevals, nevecs);
fprintf('       f           d \n');
for i=1:double(k)
fprintf('%11.3f %11.3f\n',f(i), d(i));
end
end
```
```

istate = 3 nextit = 17
nevals = 1 nevecs = 1
f           d
0.000       1.822
4.000       1.695
4.000       1.668
4.000       1.460
4.000       1.275
4.000       1.132

istate = 4 nextit = 30
nevals = 4 nevecs = 4
f           d
0.000       1.822
0.000       1.695
0.000       1.668
0.000       1.460
4.000       1.275
4.000       1.153

mOut =

4

noitsOut =

30

xOut =

-0.1189    0.2153    0.1648    0.1561   -0.4139    0.5277
0.1378   -0.1741    0.1858   -0.1931    0.2347   -0.2991
-0.1389   -0.1626   -0.1763    0.3005    0.1768   -0.2255
0.1343    0.1602   -0.2227   -0.2058   -0.1628    0.2075
-0.2012    0.3217    0.3010    0.1253   -0.2233    0.2848
0.2235   -0.2761    0.2954   -0.0744    0.0672   -0.0856
-0.2242   -0.2692   -0.2899    0.2312    0.1936   -0.2468
0.2093    0.2914   -0.3320   -0.1018   -0.1668    0.2126
-0.2093    0.2914    0.3320   -0.1018    0.1668   -0.2126
0.2242   -0.2692    0.2899    0.2312   -0.1936    0.2468
-0.2235   -0.2761   -0.2954   -0.0744   -0.0667    0.0850
0.2012    0.3217   -0.3010    0.1253    0.2241   -0.2856
-0.1343    0.1602    0.2227   -0.2058    0.1621   -0.2067
0.1389   -0.1626    0.1763    0.3005   -0.1770    0.2256
-0.1378   -0.1741   -0.1858   -0.1931   -0.2338    0.2981
0.1189    0.2153   -0.1648    0.1561    0.4142   -0.5280

d =

1.8223
1.6949
1.6684
1.4600
1.2748
1.1528

user =

[16x16 double]    [16x16 double]

ifail =

0

```
```function f02fj_example
n = 16;
m = 4;
noits = int64(1000);
tol = 0.0001;
novecs = int64(0);
x = zeros(n, m+2);
% a, b will be passed in user
a = diag(ones(n,1))+diag(-0.25*ones(n-1,1),1)+diag(-0.25*ones(n-1,1),-1)+...
diag(-0.25*ones(n-4,1),4)+diag(-0.25*ones(n-4,1),-4);
b = diag(ones(n,1))+diag(-0.5*ones(n-1,1),1)+diag(-0.5*ones(n-1,1),-1);
[mOut, noitsOut, xOut, d, user, ifail] = ...
f02fj(int64(m), noits, tol, @dot, @image, @monit, novecs, x, 'user', {a, b})

function [result, iflag, user] = dot(iflag, n, z, w, user)
b = user{2};
result=transpose(w)*b*z;
function [iflag, w, user] = image(iflag, n, z, user)

a=user{1};
b=user{2};

w=inv(a)*b*z;
function monit(istate, nextit, nevals, nevecs, k, f, d)

if (istate ~= 0)
fprintf('\n  istate = %d nextit = %d\n', istate, nextit);
fprintf('  nevals = %d nevecs = %d\n', nevals, nevecs);
fprintf('       f           d \n');
for i=1:double(k)
fprintf('%11.3f %11.3f\n',f(i), d(i));
end
end
```
```

istate = 3 nextit = 17
nevals = 1 nevecs = 1
f           d
0.000       1.822
4.000       1.695
4.000       1.668
4.000       1.460
4.000       1.275
4.000       1.132

istate = 4 nextit = 30
nevals = 4 nevecs = 4
f           d
0.000       1.822
0.000       1.695
0.000       1.668
0.000       1.460
4.000       1.275
4.000       1.153

mOut =

4

noitsOut =

30

xOut =

-0.1189    0.2153    0.1648    0.1561   -0.4139    0.5277
0.1378   -0.1741    0.1858   -0.1931    0.2347   -0.2991
-0.1389   -0.1626   -0.1763    0.3005    0.1768   -0.2255
0.1343    0.1602   -0.2227   -0.2058   -0.1628    0.2075
-0.2012    0.3217    0.3010    0.1253   -0.2233    0.2848
0.2235   -0.2761    0.2954   -0.0744    0.0672   -0.0856
-0.2242   -0.2692   -0.2899    0.2312    0.1936   -0.2468
0.2093    0.2914   -0.3320   -0.1018   -0.1668    0.2126
-0.2093    0.2914    0.3320   -0.1018    0.1668   -0.2126
0.2242   -0.2692    0.2899    0.2312   -0.1936    0.2468
-0.2235   -0.2761   -0.2954   -0.0744   -0.0667    0.0850
0.2012    0.3217   -0.3010    0.1253    0.2241   -0.2856
-0.1343    0.1602    0.2227   -0.2058    0.1621   -0.2067
0.1389   -0.1626    0.1763    0.3005   -0.1770    0.2256
-0.1378   -0.1741   -0.1858   -0.1931   -0.2338    0.2981
0.1189    0.2153   -0.1648    0.1561    0.4142   -0.5280

d =

1.8223
1.6949
1.6684
1.4600
1.2748
1.1528

user =

[16x16 double]    [16x16 double]

ifail =

0

```