12 Matrices

 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

The routines for matrix operations in this library are:

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:

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

        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

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

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

        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

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

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.