NAG Library Routine Document
f04qaf (real_gen_sparse_lsqsol)
1
Purpose
f04qaf solves sparse nonsymmetric equations, sparse linear least squares problems and sparse damped linear least squares problems, using a Lanczos algorithm.
2
Specification
Fortran Interface
Subroutine f04qaf ( 
m, n, b, x, se, aprod, damp, atol, btol, conlim, itnlim, msglvl, itn, anorm, acond, rnorm, arnorm, xnorm, work, ruser, lruser, iuser, liuser, inform, ifail) 
Integer, Intent (In)  ::  m, n, msglvl, lruser, liuser  Integer, Intent (Inout)  ::  itnlim, iuser(liuser), ifail  Integer, Intent (Out)  ::  itn, inform  Real (Kind=nag_wp), Intent (In)  ::  damp, atol, btol, conlim  Real (Kind=nag_wp), Intent (Inout)  ::  b(m), ruser(lruser)  Real (Kind=nag_wp), Intent (Out)  ::  x(n), se(n), anorm, acond, rnorm, arnorm, xnorm, work(n,2)  External  ::  aprod 

C Header Interface
#include nagmk26.h
void 
f04qaf_ (const Integer *m, const Integer *n, double b[], double x[], double se[], void (NAG_CALL *aprod)(Integer *mode, const Integer *m, const Integer *n, double x[], double y[], double ruser[], const Integer *lruser, Integer iuser[], const Integer *liuser), const double *damp, const double *atol, const double *btol, const double *conlim, Integer *itnlim, const Integer *msglvl, Integer *itn, double *anorm, double *acond, double *rnorm, double *arnorm, double *xnorm, double work[], double ruser[], const Integer *lruser, Integer iuser[], const Integer *liuser, Integer *inform, Integer *ifail) 

3
Description
f04qaf can be used to solve a system of linear equations
where
$A$ is an
$n$ by
$n$ sparse nonsymmetric matrix, or can be used to solve linear least squares problems, so that
f04qaf minimizes the value
$\rho $ given by
where
$A$ is an
$m$ by
$n$ sparse matrix and
$\Vert r\Vert $ denotes the Euclidean length of
$r$ so that
${\Vert r\Vert}^{2}={r}^{\mathrm{T}}r$. A damping argument,
$\lambda $, may be included in the least squares problem in which case
f04qaf minimizes the value
$\rho $ given by
$\lambda $ is supplied as the argument
damp and should of course be zero if the solution to problems
(1) or
(2) is required. Minimizing
$\rho $ in
(3) is often called ridge regression.
f04qaf is based upon algorithm LSQR (see
Paige and Saunders (1982a) and
Paige and Saunders (1982b)) and solves the problems by an algorithm based upon the Lanczos process. The routine does not require
$A$ explicitly, but
$A$ is specified via
aprod which must perform the operations
$\left(y+Ax\right)$ and
$\left(x+{A}^{\mathrm{T}}y\right)$ for a given
$n$element vector
$x$ and
$m$ element vector
$y$. A argument to
aprod specifies which of the two operations is required on a given entry.
The routine also returns estimates of the standard errors of the sample regression coefficients (
${x}_{\mathit{i}}$, for
$\mathit{i}=1,2,\dots ,n$) given by the diagonal elements of the estimated variancecovariance matrix
$V$. When problem
(2) is being solved and
$A$ is of full rank, then
$V$ is given by
and when problem
(3) is being solved then
$V$ is given by
Let
$\stackrel{}{A}$ denote the matrix
let
$\stackrel{}{r}$ denote the residual vector
corresponding to an iterate
$x$, so that
$\rho =\Vert \stackrel{}{r}\Vert $ is the function being minimized, and let
$\Vert A\Vert $ denote the Frobenius (Euclidean) norm of
$A$. Then the routine accepts
$x$ as a solution if it is estimated that one of the following two conditions is satisfied:
where
${\mathit{tol}}_{1}$ and
${\mathit{tol}}_{2}$ are usersupplied tolerances which estimate the relative errors in
$A$ and
$b$ respectively. Condition
(6) is appropriate for compatible problems where, in theory, we expect the residual to be zero and will be satisfied by an acceptable solution
$x$ to a compatible problem. Condition
(7) is appropriate for incompatible systems where we do not expect the residual to be zero and is based on the observation that, in theory,
when
$x$ is a solution to the least squares problem, and so
(7) will be satisfied by an acceptable solution
$x$ to a linear least squares problem.
The routine also includes a test to prevent convergence to solutions,
$x$, with unacceptably large elements. This can happen if
$A$ is nearly singular or is nearly rank deficient. If we let the singular values of
$\stackrel{}{A}$ be
then the condition number of
$\stackrel{}{A}$ is defined as
where
${\sigma}_{k}$ is the smallest nonzero singular value of
$\stackrel{}{A}$ and hence
$k$ is the rank of
$\stackrel{}{A}$. When
$k<n$, then
$\stackrel{}{A}$ is rank deficient, the least squares solution is not unique and
f04qaf will normally converge to the minimal length solution. In practice
$\stackrel{}{A}$ will not have exactly zero singular values, but may instead have small singular values that we wish to regard as zero.
The routine provides for this possibility by terminating if
where
${c}_{\mathrm{lim}}$ is a usersupplied limit on the condition number of
$\stackrel{}{A}$. For problem
(1) termination with this condition indicates that
$A$ is nearly singular and for problem
(2) indicates that
$A$ is nearly rank deficient and so has near linear dependencies in its columns. In this case inspection of
$\Vert r\Vert $,
$\Vert {A}^{\mathrm{T}}r\Vert $ and
$\Vert x\Vert $, which are all returned by the routine, will indicate whether or not an acceptable solution has been found. Condition
(8), perhaps in conjunction with
$\lambda \ne 0$, can be used to try and ‘regularize’ least squares solutions. A full discussion of the stopping criteria is given in Section 6 of
Paige and Saunders (1982a).
Introduction of a nonzero damping argument
$\lambda $ tends to reduce the size of the computed solution and to make its components less sensitive to changes in the data, and
f04qaf is applicable when a value of
$\lambda $ is known
a priori. To have an effect,
$\lambda $ should normally be at least
$\sqrt{\epsilon}\Vert A\Vert $ where
$\epsilon $ is the
machine precision. For further discussion see
Paige and Saunders (1982b) and the references given there.
Whenever possible the matrix $A$ should be scaled so that the relative errors in the elements of $A$ are all of comparable size. Such a scaling helps to prevent the least squares problem from being unnecessarily sensitive to data errors and will normally reduce the number of iterations required. At the very least, in the absence of better information, the columns of $A$ should be scaled to have roughly equal column length.
4
References
Paige C C and Saunders M A (1982a) LSQR: An algorithm for sparse linear equations and sparse least squares ACM Trans. Math. Software 8 43–71
Paige C C and Saunders M A (1982b) Algorithm 583 LSQR: Sparse linear equations and least squares problems ACM Trans. Math. Software 8 195–209
5
Arguments
 1: $\mathbf{m}$ – IntegerInput

On entry: $m$, the number of rows of the matrix $A$.
Constraint:
${\mathbf{m}}\ge 1$.
 2: $\mathbf{n}$ – IntegerInput

On entry: $n$, the number of columns of the matrix $A$.
Constraint:
${\mathbf{n}}\ge 1$.
 3: $\mathbf{b}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: the righthand side vector $b$.
On exit:
b is overwritten.
 4: $\mathbf{x}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the solution vector $x$.
 5: $\mathbf{se}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the estimates of the standard errors of the components of
$x$. Thus
${\mathbf{se}}\left(i\right)$ contains an estimate of
$\sqrt{{\nu}_{ii}}$, where
${\nu}_{ii}$ is the
$i$th diagonal element of the estimated variancecovariance matrix
$V$. The estimates returned in
se will be the lower bounds on the actual estimated standard errors, but will usually be correct to at least one significant figure.
 6: $\mathbf{aprod}$ – Subroutine, supplied by the user.External Procedure

aprod must perform the operations
$y:=y+Ax$ and
$x:=x+{A}^{\mathrm{T}}y$ for given vectors
$x$ and
$y$.
The specification of
aprod is:
Fortran Interface
Integer, Intent (In)  ::  m, n, lruser, liuser  Integer, Intent (Inout)  ::  mode, iuser(liuser)  Real (Kind=nag_wp), Intent (Inout)  ::  x(n), y(m), ruser(lruser) 

C Header Interface
#include nagmk26.h
void 
aprod (Integer *mode, const Integer *m, const Integer *n, double x[], double y[], double ruser[], const Integer *lruser, Integer iuser[], const Integer *liuser) 

 1: $\mathbf{mode}$ – IntegerInput/Output

On entry: specifies which operation is to be performed.
 ${\mathbf{mode}}=1$
 aprod must compute $y+Ax$.
 ${\mathbf{mode}}=2$
 aprod must compute $x+{A}^{\mathrm{T}}y$.
On exit: may be used as a flag to indicate a failure in the computation of
$y+Ax$ or
$x+{A}^{\mathrm{T}}y$. If
mode is negative on exit from
aprod,
f04qaf will exit immediately with
ifail set to
mode.
 2: $\mathbf{m}$ – IntegerInput

On entry: $m$, the number of rows of $A$.
 3: $\mathbf{n}$ – IntegerInput

On entry: $n$, the number of columns of $A$.
 4: $\mathbf{x}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: the vector $x$.
On exit: if
${\mathbf{mode}}=1$,
x must be unchanged.
If
${\mathbf{mode}}=2$,
x must contain
$x+{A}^{\mathrm{T}}y$.
 5: $\mathbf{y}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: the vector $y$.
On exit: if
${\mathbf{mode}}=1$,
y must contain
$y+Ax$.
If
${\mathbf{mode}}=2$,
y must be unchanged.
 6: $\mathbf{ruser}\left({\mathbf{lruser}}\right)$ – Real (Kind=nag_wp) arrayUser Workspace
 7: $\mathbf{lruser}$ – IntegerInput
 8: $\mathbf{iuser}\left({\mathbf{liuser}}\right)$ – Integer arrayUser Workspace
 9: $\mathbf{liuser}$ – IntegerInput

aprod is called with the arguments
ruser,
lruser,
iuser and
liuser as supplied to
f04qaf. You should use the arrays
ruser,
lruser,
iuser and
liuser to supply information to
aprod.
aprod must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
f04qaf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: aprod should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
f04qaf. If your code inadvertently
does return any NaNs or infinities,
f04qaf is likely to produce unexpected results.
 7: $\mathbf{damp}$ – Real (Kind=nag_wp)Input

On entry: the value
$\lambda $. If either problem
(1) or problem
(2) is to be solved,
damp must be supplied as zero.
 8: $\mathbf{atol}$ – Real (Kind=nag_wp)Input

On entry: the tolerance,
${\mathit{tol}}_{1}$, of the convergence criteria
(6) and
(7); it should be an estimate of the largest relative error in the elements of
$A$. For example, if the elements of
$A$ are correct to about
$4$ significant figures, then
atol should be set to about
$5\times {10}^{4}$. If
atol is supplied as less than
$\epsilon $, where
$\epsilon $ is the
machine precision, the value
$\epsilon $ is used instead of
atol.
 9: $\mathbf{btol}$ – Real (Kind=nag_wp)Input

On entry: the tolerance,
${\mathit{tol}}_{2}$, of the convergence criterion
(6); it should be an estimate of the largest relative error in the elements of
$B$. For example, if the elements of
$B$ are correct to about
$4$ significant figures, then
btol should be set to about
$5\times {10}^{4}$. If
btol is supplied as less than
$\epsilon $, the value
$\epsilon $ is used instead of
btol.
 10: $\mathbf{conlim}$ – Real (Kind=nag_wp)Input

On entry: the value
${c}_{\mathrm{lim}}$ of equation
(8); it should be an upper limit on the condition number of
$\stackrel{}{A}$.
conlim should not normally be chosen much larger than
$1.0/{\mathbf{atol}}$. If
conlim is supplied as zero, the value
$1.0/\epsilon $ is used instead of
conlim.
 11: $\mathbf{itnlim}$ – IntegerInput/Output

On entry: an upper limit on the number of iterations. If
${\mathbf{itnlim}}\le 0$, the value
n is used in place of
itnlim, but for illconditioned problems a higher value of
itnlim is likely to be necessary.
On exit: unchanged unless
${\mathbf{itnlim}}\le 0$ on entry, in which case it is set to
n.
 12: $\mathbf{msglvl}$ – IntegerInput

On entry: the level of printing from
f04qaf. If
${\mathbf{msglvl}}\le 0$, then no printing occurs, but otherwise messages will be output on the advisory message channel (see
x04abf). A description of the printed output is given in
Section 9.1. The level of printing is determined as follows:
 ${\mathbf{msglvl}}\le 0$
 No printing.
 ${\mathbf{msglvl}}=1$
 A brief summary is printed just prior to return from f04qaf.
 ${\mathbf{msglvl}}\ge 2$
 A summary line is printed periodically to monitor the progress of f04qaf, together with a brief summary just prior to return from f04qaf.
 13: $\mathbf{itn}$ – IntegerOutput

On exit: the number of iterations performed.
 14: $\mathbf{anorm}$ – Real (Kind=nag_wp)Output

On exit: an estimate of
$\Vert \stackrel{}{A}\Vert $ for the matrix
$\stackrel{}{A}$ of
(4).
 15: $\mathbf{acond}$ – Real (Kind=nag_wp)Output

On exit: an estimate of $\mathrm{cond}\left(\stackrel{}{A}\right)$ which is a lower bound.
 16: $\mathbf{rnorm}$ – Real (Kind=nag_wp)Output

On exit: an estimate of
$\Vert \stackrel{}{r}\Vert $ for the residual,
$\stackrel{}{r}$, of
(5) corresponding to the solution
$x$ returned in
x. Note that
$\Vert \stackrel{}{r}\Vert $ is the function being minimized.
 17: $\mathbf{arnorm}$ – Real (Kind=nag_wp)Output

On exit: an estimate of the
$\Vert {\stackrel{}{A}}^{\mathrm{T}}\stackrel{}{r}\Vert $ corresponding to the solution
$x$ returned in
x.
 18: $\mathbf{xnorm}$ – Real (Kind=nag_wp)Output

On exit: an estimate of
$\Vert x\Vert $ for the solution
$x$ returned in
x.
 19: $\mathbf{work}\left({\mathbf{n}},2\right)$ – Real (Kind=nag_wp) arrayWorkspace

 20: $\mathbf{ruser}\left({\mathbf{lruser}}\right)$ – Real (Kind=nag_wp) arrayUser Workspace

ruser is not used by
f04qaf, but is passed directly to
aprod and may be used to pass information to this routine.
 21: $\mathbf{lruser}$ – IntegerInput

On entry: the dimension of the array
ruser as declared in the (sub)program from which
f04qaf is called.
Constraint:
${\mathbf{lruser}}\ge 1$.
 22: $\mathbf{iuser}\left({\mathbf{liuser}}\right)$ – Integer arrayUser Workspace

iuser is not used by
f04qaf, but is passed directly to
aprod and may be used to pass information to this routine.
 23: $\mathbf{liuser}$ – IntegerInput

On entry: the dimension of the array
iuser as declared in the (sub)program from which
f04qaf is called.
Constraint:
${\mathbf{liuser}}\ge 1$.
 24: $\mathbf{inform}$ – IntegerOutput

On exit: the reason for termination of
f04qaf.
 ${\mathbf{inform}}=0$
 The exact solution is $x=0$. No iterations are performed in this case.
 ${\mathbf{inform}}=1$
 The termination criterion of (6) has been satisfied with ${\mathit{tol}}_{1}$ and ${\mathit{tol}}_{2}$ as the values supplied in atol and btol respectively.
 ${\mathbf{inform}}=2$
 The termination criterion of (7) has been satisfied with ${\mathit{tol}}_{1}$ as the value supplied in atol.
 ${\mathbf{inform}}=3$
 The termination criterion of (6) has been satisfied with ${\mathit{tol}}_{1}$ and/or ${\mathit{tol}}_{2}$ as the value $\epsilon $, where $\epsilon $ is the machine precision. One or both of the values supplied in atol and btol must have been less than $\epsilon $ and was too small for this machine.
 ${\mathbf{inform}}=4$
 The termination criterion of (7) has been satisfied with ${\mathit{tol}}_{1}$ as the value $\epsilon $. The value supplied in atol must have been less than $\epsilon $ and was too small for this machine.
The values
${\mathbf{inform}}=5$,
$6$ and
$7$ correspond to failure with
${\mathbf{ifail}}={\mathbf{2}}$,
${\mathbf{3}}$ and
${\mathbf{4}}$ respectively (see
Section 6) and when
ifail is negative
inform will be set to the same negative value.
 25: $\mathbf{ifail}$ – IntegerInput/Output

On entry:
ifail must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is
$0$.
When the value $\mathbf{1}\text{ or}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}<0$

A negative value of
ifail indicates an exit from
f04qaf because you have set
mode negative in
aprod. The value of
ifail will be the same as your setting of
mode.
 ${\mathbf{ifail}}=1$

On entry,  ${\mathbf{m}}<1$, 
or  ${\mathbf{n}}<1$, 
or  ${\mathbf{lruser}}<1$, 
or  ${\mathbf{liuser}}<1$. 
 ${\mathbf{ifail}}=2$

The condition of
(8) has been satisfied for the value of
${c}_{\mathrm{lim}}$ supplied in
conlim. If this failure is unexpected you should check that
aprod is working correctly. Although conditions
(6) or
(7) have not been satisfied, the values returned in
rnorm,
arnorm and
xnorm may nevertheless indicate that an acceptable solution has been reached.
 ${\mathbf{ifail}}=3$

The condition of
(8) has been satisfied for the value
${c}_{\mathrm{lim}}=1.0/\epsilon $, where
$\epsilon $ is the
machine precision. The matrix
$\stackrel{}{A}$ is nearly singular or rank deficient and the problem is too illconditioned for this machine. If this failure is unexpected, you should check that
aprod is working correctly.
 ${\mathbf{ifail}}=4$

The limit on the number of iterations has been reached. The number of iterations required by f04qaf and the condition of the matrix $\stackrel{}{A}$ can depend strongly on the scaling of the problem. Poor scaling of the rows and columns of $A$ should be avoided whenever possible.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
When the problem is compatible, the computed solution
$x$ will satisfy the equation
where an estimate of
$\Vert r\Vert $ is returned in the argument
rnorm. When the problem is incompatible, the computed solution
$x$ will satisfy the equation
where an estimate of
$\Vert e\Vert $ is returned in the argument
arnorm. See also Section 6.2 of
Paige and Saunders (1982b).
8
Parallelism and Performance
f04qaf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
The time taken by
f04qaf is likely to be principally determined by the time taken in
aprod, which is called twice on each iteration, once with
${\mathbf{mode}}=1$ and once with
${\mathbf{mode}}=2$. The time taken per iteration by the remaining operations in
f04qaf is approximately proportional to
$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(m,n\right)$.
The Lanczos process will usually converge more quickly if
$A$ is preconditioned by a nonsingular matrix
$M$ that approximates
$A$ in some sense and is also chosen so that equations of the form
$My=c$ can efficiently be solved for
$y$. For a discussion of preconditioning, see the
F11 Chapter Introduction. In the context of
f04qaf, problem
(1) is equivalent to
and problem
(2) is equivalent to minimizing
Note that the normal matrix
${\left(A{M}^{1}\right)}^{\mathrm{T}}\left(A{M}^{1}\right)={M}^{\mathrm{T}}\left({A}^{\mathrm{T}}A\right){M}^{1}$ so that the preconditioning
$A{M}^{1}$ is equivalent to the preconditioning
${M}^{\mathrm{T}}\left({A}^{\mathrm{T}}A\right){M}^{1}$ of the normal matrix
${A}^{\mathrm{T}}A$.
Preconditioning can be incorporated into
f04qaf simply by coding
aprod to compute
$y+A{M}^{1}x$ and
$x+{M}^{\mathrm{T}}{A}^{\mathrm{T}}y$ in place of
$y+Ax$ and
$x+{A}^{\mathrm{T}}y$ respectively, and then solving the equations
$Mx=y$ for
$x$ on return from
f04qaf. The quantity
$y+A{M}^{1}x$ should be computed by solving
$Mz=x$ for
$z$ and then computing
$y+Az$, and
$x+{M}^{\mathrm{T}}{A}^{\mathrm{T}}y$ should be computed by solving
${M}^{\mathrm{T}}z={A}^{\mathrm{T}}y$ for
$z$ and then forming
$x+z$.
9.1
Description of the Printed Output
When
${\mathbf{msglvl}}>0$, then
f04qaf will produce output (except in the case where the routine fails with
${\mathbf{ifail}}={\mathbf{1}}$) on the advisory message channel (see
x04abf).
When
${\mathbf{msglvl}}\ge 2$ then a summary line is printed periodically giving the following information:
Output 
Meaning 
ITN 
Iteration number, $k$. 
X(1) 
The first element of the current iterate ${x}_{k}$. 
FUNCTION 
The current value of the function, $\rho $, being minimized. 
COMPAT 
An estimate of $\Vert {\stackrel{}{r}}_{k}\Vert /\Vert b\Vert $, where ${\stackrel{}{r}}_{k}$ is the residual corresponding to ${x}_{k}$. This value should converge to zero (in theory) if and only if the problem is compatible. COMPAT decreases monotonically. 
INCOMPAT 
An estimate of $\Vert {\stackrel{}{A}}^{\mathrm{T}}{\stackrel{}{r}}_{k}\Vert /\left(\Vert \stackrel{}{A}\Vert \Vert {\stackrel{}{r}}_{k}\Vert \right)$ which should converge to zero if and only if at the solution $\rho $ is nonzero. INCOMPAT is not usually monotonic. 
NRM(ABAR) 
A monotonically increasing estimate of $\Vert \stackrel{}{A}\Vert $. 
COND(ABAR) 
A monotonically increasing estimate of the condition number $\mathrm{cond}\left(\stackrel{}{A}\right)$. 
10
Example
This example solves the linear least squares problem
where
$A$ is the
$13$ by
$12$ matrix and
$b$ is the
$13$ element vector given by
with
$h=0.1$.
Such a problem can arise by considering the Neumann problem on a rectangle
where
$C$ is the boundary of the rectangle, and discretizing as illustrated below with the square mesh
Figure 1
The $12$ by $12$ symmetric part of $A$ represents the difference equations and the final row comes from the normalizing condition. The example program has $g\left(x,y\right)=1$ at all the internal mesh points, but apart from this is written in a general manner so that the number of rows (NROWS) and columns (NCOLS) in the grid can readily be altered.
10.1
Program Text
Program Text (f04qafe.f90)
10.2
Program Data
Program Data (f04qafe.d)
10.3
Program Results
Program Results (f04qafe.r)