next up previous index
Next: F012 Symmetric Positive-Definite Up: CERNLIB Previous: F010 Linear Equations

F011 Repeated Solution of Linear Equations, Matrix Inversion, Determinant

Routine ID: F011
Author(s): G.A. Erskine, H. LippsLibrary: KERNLIB
Submitter: Submitted: 18.12.1979
Language: Fortran or Assembler or COMPASSRevised: 27.11.1984

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

AX=B&sp;(*)

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 A-1 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 (N001), ABEND (Z035)

Usage:

For t=R (type REAL), t=D (type DOUBLE PRECISION), t=C (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 K > 1 ).
IR
( INTEGER) Array of at least N elements, required as working space.
IFAIL
( INTEGER) On exit, IFAIL will be set to -1 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:
= -1 if det( A) is probably too small,
= +1 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 K = 1 .

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 -1 , 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 -1 and DET will be set to zero.

    If the evaluation of det( A) would probably cause overflow, JFAIL will be set to +1 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 A-1 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 A-1 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 N < 1 or IDIM < N or K < 1 , a message is printed and program execution is terminated by calling ABEND (Z035).

Examples:

Assume that the 10 x10 matrix A, the 10 x3 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 DET=det(A) , replace B by A-1B , replace z by A-1z , and replace A by A-1 :
    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)

F012



next up previous index
Next: F012 Symmetric Positive-Definite Up: CERNLIB Previous: F010 Linear Equations


Janne Saarela
Mon Apr 3 15:06:23 METDST 1995