NAG Toolbox: nag_eigen_real_symm_sparse_eigsys (f02fj)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{m}$ – int64int32nag_int scalar

$m$, the number of eigenvalues required.

Constraint: $\mathbf{m}\ge 1$.

2:     $\mathrm{noits}$ – int64int32nag_int scalar

The maximum number of major iterations (eigenvalue sub-problems) to be performed. If $\mathbf{noits}\le 0$, the value $100$ is used in place of noits.

3:     $\mathrm{tol}$ – double scalar

A relative tolerance to be used in accepting eigenvalues and eigenvectors. If the eigenvalues are required to about $t$ significant figures, tol should be set to about ${10}^{-t}$. $d_i$ is accepted as an eigenvalue as soon as two successive approximations to $d_i$ differ by less than $\left(\left|{\tilde{d}}_i\right|\times \mathbf{tol}\right)/10$, where ${\tilde{d}}_i$ is the latest approximation to $d_i$. Once an eigenvalue has been accepted, an eigenvector is accepted as soon as $\left(d_i{f}_i\right)/\left(d_i-d_k\right)<\mathbf{tol}$, where $f_i$ is the normalized residual of the current approximation to the eigenvector (see Further Comments for further information). The values of the $f_i$ and $d_i$ can be printed from monit. If tol is supplied outside the range ($\epsilon ,1.0$), where $\epsilon$ is the machine precision, the value $\epsilon$ is used in place of tol.

4:     $\mathrm{dot}$ – function handle or string containing name of m-file

dot must return the value $w^{\mathrm{T}}Bz$ for given vectors $w$ and $z$. For the standard eigenvalue problem, where $B=I$, dot must return the dot product $w^{\mathrm{T}}z$.

[result, iflag, user] = dot(iflag, n, z, w, user)

Input Parameters

1:     $\mathrm{iflag}$ – int64int32nag_int scalar

Is always non-negative.

2:     $\mathrm{n}$ – int64int32nag_int scalar

The number of elements in the vectors $z$ and $w$ and the order of the matrix $B$.

3:     $\mathrm{z}\left(\mathbf{n}\right)$ – double array

The vector $z$ for which $w^{\mathrm{T}}Bz$ is required.

4:     $\mathrm{w}\left(\mathbf{n}\right)$ – double array

The vector $w$ for which $w^{\mathrm{T}}Bz$ is required.

5:     $\mathrm{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:     $\mathrm{result}$ – double scalar

result returns the value $w^{\mathrm{T}}Bz$ for given vectors $w$ and $z$.

2:     $\mathrm{iflag}$ – int64int32nag_int scalar

May be used as a flag to indicate a failure in the computation of $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:     $\mathrm{user}$ – Any MATLAB object

5:     $\mathrm{image}$ – function handle or string containing name of m-file

image must return the vector $w=Cz$ for a given vector $z$.

[iflag, w, user] = image(iflag, n, z, user)

Input Parameters

1:     $\mathrm{iflag}$ – int64int32nag_int scalar

Is always non-negative.

2:     $\mathrm{n}$ – int64int32nag_int scalar

$n$, the number of elements in the vectors $w$ and $z$, and the order of the matrix $C$.

3:     $\mathrm{z}\left(\mathbf{n}\right)$ – double array

The vector $z$ for which $Cz$ is required.

4:     $\mathrm{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:     $\mathrm{iflag}$ – int64int32nag_int scalar

May be used as a flag to indicate a failure in the computation of $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:     $\mathrm{w}\left(\mathbf{n}\right)$ – double array

The vector $w=Cz$.

3:     $\mathrm{user}$ – Any MATLAB object

6:     $\mathrm{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 arguments istate and nextit allow selective printing by monit.

monit(istate, nextit, nevals, nevecs, k, f, d)

Input Parameters

1:     $\mathrm{istate}$ – int64int32nag_int scalar

Specifies the state of nag_eigen_real_symm_sparse_eigsys (f02fj).

$\mathbf{istate}=0$

No eigenvalue or eigenvector has just been accepted.

$\mathbf{istate}=1$

One or more eigenvalues have been accepted since the last call to monit.

$\mathbf{istate}=2$

One or more eigenvectors have been accepted since the last call to monit.

$\mathbf{istate}=3$

One or more eigenvalues and eigenvectors have been accepted since the last call to monit.

$\mathbf{istate}=4$

Return from nag_eigen_real_symm_sparse_eigsys (f02fj) is about to occur.

2:     $\mathrm{nextit}$ – int64int32nag_int scalar

The number of the next iteration.

3:     $\mathrm{nevals}$ – int64int32nag_int scalar

The number of eigenvalues accepted so far.

4:     $\mathrm{nevecs}$ – int64int32nag_int scalar

The number of eigenvectors accepted so far.

5:     $\mathrm{k}$ – int64int32nag_int scalar

$k$, the number of simultaneous iteration vectors.

6:     $\mathrm{f}\left(\mathbf{k}\right)$ – double array

A vector of error quantities measuring the state of convergence of the simultaneous iteration vectors. See tol and Further Comments for further details. Each element of f is initially set to the value $4.0$ and an element remains at $4.0$ until the corresponding vector is tested.

7:     $\mathrm{d}\left(\mathbf{k}\right)$ – double array

$\mathbf{d}\left(i\right)$ contains the latest approximation to the absolute value of the $i$th eigenvalue of $C$.

7:     $\mathrm{novecs}$ – int64int32nag_int scalar

The number of approximate vectors that are being supplied in x. If novecs is outside the range $\left(0,\mathbf{k}\right)$, the value $0$ is used in place of novecs.

8:     $\mathrm{x}\left(\mathit{ldx},\mathbf{k}\right)$ – double array

ldx, the first dimension of the array, must satisfy the constraint $\mathit{ldx}\ge \mathbf{n}$.

If $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$. 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:     $\mathrm{n}$ – int64int32nag_int scalar

Default: the first dimension of the array x.

$n$, the order of the matrix $C$.

Constraint: $\mathbf{n}\ge 1$.

2:     $\mathrm{k}$ – int64int32nag_int scalar

Suggested value: $\mathbf{k}=\mathbf{m}+4$ will often be a reasonable choice in the absence of better information.

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.

Constraint: $\mathbf{m}<\mathbf{k}\le \mathbf{n}$.

3:     $\mathrm{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.

Output Parameters

1:     $\mathrm{m}$ – int64int32nag_int scalar

$m^{\prime }$, the number of eigenvalues actually found. It is equal to $m$ if $\mathbf{ifail}=\mathbf{0}$ on exit, and is less than $m$ if $\mathbf{ifail}=\mathbf{2}$, $\mathbf{3}$ or $\mathbf{4}$. See Error Indicators and Warnings and Further Comments for further information.

2:     $\mathrm{noits}$ – int64int32nag_int scalar

The number of iterations actually performed.

3:     $\mathrm{x}\left(\mathit{ldx},\mathbf{k}\right)$ – double array

If $\mathbf{ifail}=\mathbf{0}$, $\mathbf{2}$, $\mathbf{3}$ or $\mathbf{4}$, the first $m^{\prime }$ columns contain the eigenvectors corresponding to the eigenvalues returned in the first $m^{\prime }$ elements of d; and the next $k-m^{\prime }-1$ columns contain approximations to the eigenvectors corresponding to the approximate eigenvalues returned in the next $k-m^{\prime }-1$ elements of d. Here $m^{\prime }$ is the value returned in m, the number of eigenvalues actually found. The $k$th column is used as workspace.

4:     $\mathrm{d}\left(\mathbf{k}\right)$ – double array

If $\mathbf{ifail}=\mathbf{0}$, $\mathbf{2}$, $\mathbf{3}$ or $\mathbf{4}$, the first $m^{\prime }$ elements contain the first $m^{\prime }$ eigenvalues in decreasing order of magnitude; and the next $k-m^{\prime }-1$ elements contain approximations to the next $k-m^{\prime }-1$ eigenvalues. Here $m^{\prime }$ is the value returned in m, the number of eigenvalues actually found. $\mathbf{d}\left(k\right)$ contains the value $e$ where $\left(-e,e\right)$ is the latest interval over which Chebyshev acceleration is performed.

5:     $\mathrm{user}$ – Any MATLAB object

6:     $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

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

W  $\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.

   $\mathbf{ifail}=1$

On entry,	$\mathbf{n}<1$,
or	$\mathbf{m}<1$,
or	$\mathbf{m}\ge \mathbf{k}$,
or	$\mathbf{k}>\mathbf{n}$,
or	$\mathit{ldx}<\mathbf{n}$,
or	$\mathit{lwork}<3\times \mathbf{k}+\max \left(\mathbf{k}\times \mathbf{k},2\times \mathbf{n}\right)$,
or	$\mathit{lruser}<1$,
or	$\mathit{liuser}<1$.

W  $\mathbf{ifail}=2$

Not all the requested eigenvalues and vectors have been obtained. Approximations to the $r$th eigenvalue are oscillating rapidly indicating that severe cancellation is occurring in the $r$th eigenvector and so m is returned as $\left(r-1\right)$. A restart with a larger value of k may permit convergence.

W  $\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  $\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.

   $\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.

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Example

Open in the MATLAB editor: f02fj_example

function f02fj_example


fprintf('f02fj example results\n\n');

n = 16;
m = int64(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);

[m, noits, x, d, user, ifail] = ...
  f02fj(...
        m, noits, tol, @dot, @image, @monit, novecs, x, 'user', {a, b});

fprintf('\nFinal results\n\n');
disp('Eigenvalues');
disp(1./d(1:m)');
disp('Eigenvectors');
% Normalize eigenvectors before printing
disp(x(:,1:m)/diag(x(1,1:m)));



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

f02fj example results


  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

Final results

Eigenvalues
    0.5488    0.5900    0.5994    0.6850

Eigenvectors
    1.0000    1.0000    1.0000    1.0000
   -1.1586   -0.8089    1.1274   -1.2368
    1.1682   -0.7555   -1.0699    1.9252
   -1.1298    0.7444   -1.3512   -1.3185
    1.6919    1.4943    1.8266    0.8027
   -1.8797   -1.2826    1.7928   -0.4766
    1.8854   -1.2505   -1.7589    1.4809
   -1.7604    1.3538   -2.0148   -0.6525
    1.7604    1.3538    2.0148   -0.6525
   -1.8854   -1.2505    1.7589    1.4809
    1.8797   -1.2826   -1.7928   -0.4766
   -1.6919    1.4943   -1.8266    0.8027
    1.1298    0.7444    1.3512   -1.3185
   -1.1682   -0.7555    1.0699    1.9252
    1.1586   -0.8089   -1.1274   -1.2368
   -1.0000    1.0000   -1.0000    1.0000