12.1 Overview

12.2 Replacing calls to F01AAF

12.3 Replacing calls to F01ABF

12.4 Replacing calls to F04AAF

12.5 Replacing calls to F04AEF

12.6 Replacing calls to F04ANF and F01AXF

12.7 Replacing calls to F04ASF and F04ATF

12.8 Replacing calls to F04QAF

12.2 Replacing calls to F01AAF

12.3 Replacing calls to F01ABF

12.4 Replacing calls to F04AAF

12.5 Replacing calls to F04AEF

12.6 Replacing calls to F04ANF and F01AXF

12.7 Replacing calls to F04ASF and F04ATF

12.8 Replacing calls to F04QAF

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

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.

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

INTEGER IA, N, IX, IFAIL

DOUBLE PRECISION A(IA,N), X(IX,N), WORK(N)

IFAIL = 1

CALL F01AAF( A, IA, N, X, IX, WORK, IFAIL )

IF ( IFAIL .NE. 0 ) THEN

An error has occurred

END IF

DOUBLE PRECISION A(IA,N), X(IX,N), WORK(N)

IFAIL = 1

CALL F01AAF( A, IA, N, X, IX, WORK, IFAIL )

IF ( IFAIL .NE. 0 ) THEN

An error has occurred

END IF

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.

INTEGER I, J, IA, N, IX, IFAIL

INTEGER IPVT(N)

DOUBLE PRECISION A(IA,N), X(IX,N), WORK(N), DUMMY(2)

DO 2 J = 1, N

DO 1 I = 1, N

X(I,J) = A(I,J)

1 CONTINUE

2 CONTINUE

CALL PDA_DGEFA( X, IX, N, IPVT, IFAIL )

IF ( IFAIL .NE. 0 ) THEN

An error has occurred

ELSE

CALL PDA_DGEDI( X, IX, N, IPVT, DUMMY, WORK, 1 )

END IF

INTEGER IPVT(N)

DOUBLE PRECISION A(IA,N), X(IX,N), WORK(N), DUMMY(2)

DO 2 J = 1, N

DO 1 I = 1, N

X(I,J) = A(I,J)

1 CONTINUE

2 CONTINUE

CALL PDA_DGEFA( X, IX, N, IPVT, IFAIL )

IF ( IFAIL .NE. 0 ) THEN

An error has occurred

ELSE

CALL PDA_DGEDI( X, IX, N, IPVT, DUMMY, WORK, 1 )

END IF

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

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

INTEGER IA, IB, N, M, IX, IFAIL

DOUBLE PRECISION A(IA,N), B(IB,M), X(IX,M), WORK(N)

DO 2 J = 1, N

DO 1 I = 1, N

B(I,J) = 0D0

1 CONTINUE

2 CONTINUE

DO 3 I = 1, N

B(I,I) = 1D0

3 CONTINUE

IFAIL = 1

CALL F04AAF( A, IA, B, IB, N, M, X, IX, WORK, IFAIL )

IF ( IFAIL .NE. 0 ) THEN

An error has occurred

END IF

DOUBLE PRECISION A(IA,N), B(IB,M), X(IX,M), WORK(N)

DO 2 J = 1, N

DO 1 I = 1, N

B(I,J) = 0D0

1 CONTINUE

2 CONTINUE

DO 3 I = 1, N

B(I,I) = 1D0

3 CONTINUE

IFAIL = 1

CALL F04AAF( A, IA, B, IB, N, M, X, IX, WORK, IFAIL )

IF ( IFAIL .NE. 0 ) THEN

An error has occurred

END IF

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.

INTEGER I, J, IA, N, IX, IFAIL

INTEGER IPVT(N)

DOUBLE PRECISION A(IA,N), X(IX,N), WORK(N), DUMMY(2)

DO 2 J = 1, N

DO 1 I = 1, N

X(I,J) = A(I,J)

1 CONTINUE

2 CONTINUE

CALL PDA_DGEFA( X, IX, N, IPVT, IFAIL )

IF ( IFAIL .NE. 0 ) THEN

An error has occurred

ELSE

CALL PDA_DGEDI( X, IX, N, IPVT, DUMMY, WORK, 1 )

END IF

INTEGER IPVT(N)

DOUBLE PRECISION A(IA,N), X(IX,N), WORK(N), DUMMY(2)

DO 2 J = 1, N

DO 1 I = 1, N

X(I,J) = A(I,J)

1 CONTINUE

2 CONTINUE

CALL PDA_DGEFA( X, IX, N, IPVT, IFAIL )

IF ( IFAIL .NE. 0 ) THEN

An error has occurred

ELSE

CALL PDA_DGEDI( X, IX, N, IPVT, DUMMY, WORK, 1 )

END IF

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.

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.

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

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.