F011: Repeated Solution of Linear Equations, Matrix Inversion, Determinant

Author(s): G.A. Erskine, H. Lipps	Library: KERNLIB
Submitter:	Submitted: 18.12.1979
Language: Fortran or Assembler or COMPASS	Revised: 27.11.1984

These subroutines provide a two-step procedure for solving sets of linear equations

which is faster than the library programs RINV (F010) when (*) must be solved repeatedly for the same matrix A with different sets of right-hand sides. The inverse matrix and the determinant det(A) may also be calculated.

Structure:

SUBROUTINE subprograms
User Entry Names: RFACT, RFEQN, RFINV, DFACT, DFEQN, DFINV, CFACT, CFEQN, CFINV
Internal Entry Names: TMPRNT
Files Referenced: Printer
External References: KERMTR, ABEND

Usage:

For (type REAL), (type DOUBLE PRECISION), (type COMPLEX):

    CALL tFACT(N,A,IDIM,IR,IFAIL,DET,JFAIL)
    CALL tFEQN(N,A,IDIM,IR,K,B)
    CALL tFINV(N,A,IDIM,IR)

N

(INTEGER) Order of the square matrix A.

A

(Type according to t) Two-dimensional array whose first dimension has the value IDIM.

IDIM

(INTEGER) First dimension of array A (and of array B if ).

IR

(INTEGER) Array of at least N elements, required as working space.

IFAIL

(INTEGER) On exit, IFAIL will be set to if A is found to be singular, and to 0 otherwise. (Singularity will often go undetected because of rounding errors during factorization even if the elements of A have integral values.)

DET

(Type according to t) On exit, DET will be set to the value det(A) unless JFAIL returns a non-zero value.

JFAIL

(INTEGER) On exit, JFAIL will be set to zero if det(A) can be safely evaluated. Otherwise JFAIL is set as follows:
if det(A) is probably too small,
if det(A) is probably too large.

K

(INTEGER) Number of columns of the matrices B and X.

B

(Type according to t) In general, a two-dimensional array whose first dimension has the value IDIM. B may be one-dimensional if .

Subroutine tFACT must be called with matrix A in array A prior to any calls to tFEQN and tFINV. On return the situation is as follows:

1. Provided A is non-singular, IFAIL will be set to 0, and A and R will be set in preparation for calls to tFEQN and tFINV.

   If A is singular, IFAIL will be set to , in which case any subsequent call to tFEQN or tFINV will give unpredictable results.

2. Provided det(A) can be safely evaluated within the range of the computer, JFAIL will be set to 0 and and DET will be set to det(A). In particular, if A is singular, both JFAIL and DET will be set to zero.

   If the evaluation of det(A) would probably cause underflow, JFAIL will be set to and DET will be set to zero.

   If the evaluation of det(A) would probably cause overflow, JFAIL will be set to and DET will be incorrect.

   Execution continues, and subsequent calls to tFEQN and tFINV will give correct results.

Subroutine tFEQN may be called only after tFACT has been called, with the contents of A and R unchanged, and with matrix B in array B. On return, B will contain the solution X, with A and R unchanged. Therefore a single call to tFACT may be followed by several calls to tFEQN with differing B.

Subroutine tFINV may be called only after tFACT has been called, with the contents of A and R unchanged. On return, A will contain the inverse of A. Therefore, once tFINV has been called, it is no longer meaningful to call tFEQN with A as parameter.

Method:

Triangular factorization with row interchanges. The inverse matrix is the product, in reverse order, of the in-place inverses of the triangular factors. The array R holds information specifying the row interchanges.

Accuracy:

On computers with IBM 370 architecture, inner products are accumulated using double-precision arithmetic internally for arrays of type REAL and COMPLEX.

Error handling:

If or or , a message is printed and program execution is terminated by calling ABEND (Z035).

Examples:

Assume that the matrix A, the matrix B, and the 10-element vector z are stored according to the Fortran convention in arrays A, B and Z respectively of a program containing the declarations

    DIMENSION IR(25)
    COMPLEX A(25,30),B(25,10),Z(25),DET

Then, unless A is singular (which is to cause a jump to statement 100), the following statements will set , replace B by , replace z by , and replace A by :

    CALL CFACT (10,A,25,IR,IFAIL,DET,JFAIL)
    IF(IFAIL .NE. 0) GO TO 100
    CALL CFEQN(10,A,25,IR,3,B)
    CALL CFEQN(10,A,25,IR,1,Z)
    CALL CFINV(10,A,25,IR)