12 Matrices

The routines for matrix operations in this library are:

• PDA_DGEFA (SLATEC/CAMSUN) 
  Factor a matrix using Gaussian elimination. This is needed before the determinant and inverse can be calculated by PDA_DGEDI.
• PDA_DGEDI (SLATEC/CAMSUN) 
  Compute the determinant and inverse of a matrix using the factors computed by PDA_DGECO or PDA_DGEFA.
• PDA_DGEFS (SLATEC/CAMSUN) 
  Solve a general system of linear equations. This solves the problem A * x = b. A is a square matrix, x and b are vectors. The problem A * X = B where all are matrices can be solved as several systems A * x = b. This routine supports such an undertaking since it is able to re-use a previous factorisation of A.
• PDA_DBOLS (SLATEC/CAMSUN) 
  Solve the problem E * x = f (in the least squares sense) with bounds on selected x values. E is a matrix, x and f are vectors.
• PDA_LSQR
  Solves sparse unsymmetric, linear least squares and damped least squares problems

12.1 Overview

There is some excess baggage in this field in NAG, since it distinguishes approximate from accurate routines. The approximation is not in the analytical, but in the numeric sense.

Leaving F04QAF aside, the need is for:

• A matrix inverter (F04AAF is used only as an inverter).
• A solver for A * x = b where A is square and x and b are vectors. The problem A * X = B where A is n by n, X and B are n by m, can be split into m problems A * x = b all with the same A.
• The third need is a least-squares solver for A * x = b where A is not square and b has more dimensions than x.

Looking at SLATEC, four routines are needed. PDA_DGEDI computes the determinant and/or inverse of a matrix, but must be preceded by a factoriser, namely PDA_DGEFA. PDA_DGEFS solves A * x = b where A is square. It can re-use a factorisation of A from a previous run, and in that sense supports the solution of A * X = B. Finally, PDA_DBOLS solves the over-determined problem A * x = b in the least-squares sense.

12.2 Replacing calls to F01AAF

This routine is an approximate matrix inverter and would be replaced by PDA_DGEFA followed by PDA_DGEDI. The NAG code might look like

The SLATEC routines invert the matrix in situ, so A must be copied to X first. It is important that PDA_DGEDI be called only if PDA_DGEFA signals correct processing as IFAIL = 0. Otherwise PDA_DGEDI will encounter a division by zero. The last argument to PDA_DGEDI determines the returned information, 1 chooses inverted matrix but no determinant.

12.3 Replacing calls to F01ABF

This routine is a matrix inverter and would be replaced by PDA_DGEFA followed by PDA_DGEDI. Migration hints are not yet available.

12.4 Replacing calls to F04AAF

Although this routine solves A * X = B, where all three are matrices, it is actually used only to find the inverse of a matrix. In that use all three matrices are square, A is given, B is unity, and X is the inverse of A. The NAG code might look like

The SLATEC routines invert the matrix in situ, so A must be copied to X first. It is important that PDA_DGEDI be called only if PDA_DGEFA signals correct processing as IFAIL = 0. Otherwise PDA_DGEDI will encounter a division by zero. The last argument to PDA_DGEDI determines the returned information, 1 chooses inverted matrix but no determinant.

12.5 Replacing calls to F04AEF

This routine solves A * X = B. The problem would be split into a number of problems A * x_i = b_i, where x_i and b_i are corresponding columns of X and B. Each problem is solved by PDA_DGEFS, which can re-use a factorisation of A that it worked out in the first call. Migration hints are not yet available.

12.6 Replacing calls to F04ANF and F01AXF

F04ANF solves an over-determined problem A * x = b and would be replaced by PDA_DBOLS. F01AXF is an auxiliary routine. Migration hints are not yet available.

12.7 Replacing calls to F04ASF and F04ATF

These routines solve A * x = b where A is square and would be replaced by PDA_DGEFS. Migration hints are not yet available.

12.8 Replacing calls to F04QAF

The routine PDA_LSQR solves exactly the same problems as F04QAF, the only differences are that the workspace requirements and argument ordering are slightly different.