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