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