MINICERN/packlib/kernlib/kernnum/f010fort/deqinv.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/02/15 17:48:48  mclareni
   6 * Kernlib
   7 *
   8 *
   9 #include "kernnum/pilot.h"
  10       SUBROUTINE DEQINV(N,A,IDIM,R,IFAIL,K,B)
  11       REAL R(N),T1,T2,T3
  12       DOUBLE PRECISION A(IDIM,N),B(IDIM,K),DET,TEMP,S,
  13      $                 B1,B2,C11,C12,C13,C21,C22,C23,C31,C32,C33
  14       CHARACTER*6 NAME
  15       DATA NAME/'DEQINV'/,KPRNT/1/
  16 C
  17 C     ******************************************************************
  18 C
  19 C     REPLACES B BY THE SOLUTION X OF A*X=B, AND REPLACES A BY ITS IN-
  20 C     VERSE.
  21 C
  22 C     N            ORDER OF THE SQUARE MATRIX IN ARRAY A.
  23 C
  24 C     A            (DOUBLE PRECISION) TWO-DIMENSIONAL ARRAY CONTAINING
  25 C                  AN N BY N MATRIX.
  26 C
  27 C     IDIM         FIRST DIMENSION PARAMETER OF ARRAYS A AND B.
  28 C
  29 C     R            (REAL) WORKING VECTOR OF LENGTH NOT LESS THAN N.
  30 C
  31 C     IFAIL        OUTPUT PARAMETER.   IFAIL= 0 ... NORMAL EXIT.
  32 C                                      IFAIL=-1 ... SINGULAR MATRIX.
  33 C
  34 C     K            NUMBER OF COLUMNS OF THE MATRIX IN ARRAY B.
  35 C
  36 C     B            (DOUBLE PRECISION) TWO-DIMENSIONAL ARRAY CONTAINING
  37 C                  AN N BY K MATRIX.
  38 C
  39 C     CALLS ... DFACT, DFINV, F010PR, ABEND.
  40 C
  41 C     ******************************************************************
  42 C
  43 C  TEST FOR PARAMETER ERRORS.
  44 C
  45       IF((N.LT.1).OR.(N.GT.IDIM).OR.(K.LT.1)) GO TO 10
  46 C
  47 C  TEST FOR N.LE.3.
  48 C
  49       IF(N.GT.3) GO TO 9
  50       IFAIL=0
  51       IF(N.LT.3) GO TO 5
  52 C
  53 C  N=3 CASE.
  54 C
  55 C     COMPUTE COFACTORS.
  56       C11=A(2,2)*A(3,3)-A(2,3)*A(3,2)
  57       C12=A(2,3)*A(3,1)-A(2,1)*A(3,3)
  58       C13=A(2,1)*A(3,2)-A(2,2)*A(3,1)
  59       C21=A(3,2)*A(1,3)-A(3,3)*A(1,2)
  60       C22=A(3,3)*A(1,1)-A(3,1)*A(1,3)
  61       C23=A(3,1)*A(1,2)-A(3,2)*A(1,1)
  62       C31=A(1,2)*A(2,3)-A(1,3)*A(2,2)
  63       C32=A(1,3)*A(2,1)-A(1,1)*A(2,3)
  64       C33=A(1,1)*A(2,2)-A(1,2)*A(2,1)
  65       T1=ABS(SNGL(A(1,1)))
  66       T2=ABS(SNGL(A(2,1)))
  67       T3=ABS(SNGL(A(3,1)))
  68 C
  69 C     (SET TEMP=PIVOT AND DET=PIVOT*DET.)
  70       IF(T1.GE.T2) GO TO 1
  71          IF(T3.GE.T2) GO TO 2
  72 C        (PIVOT IS A21)
  73             TEMP=A(2,1)
  74             DET=C13*C32-C12*C33
  75             GO TO 3
  76     1 IF(T3.GE.T1) GO TO 2
  77 C     (PIVOT IS A11)
  78          TEMP=A(1,1)
  79          DET=C22*C33-C23*C32
  80          GO TO 3
  81 C     (PIVOT IS A31)
  82     2    TEMP=A(3,1)
  83          DET=C23*C12-C22*C13
  84 C
  85 C     SET ELEMENTS OF INVERSE IN A.
  86     3 IF(DET.EQ.0D0) GO TO 11
  87       S=TEMP/DET
  88       A(1,1)=S*C11
  89       A(1,2)=S*C21
  90       A(1,3)=S*C31
  91       A(2,1)=S*C12
  92       A(2,2)=S*C22
  93       A(2,3)=S*C32
  94       A(3,1)=S*C13
  95       A(3,2)=S*C23
  96       A(3,3)=S*C33
  97 C
  98 C     REPLACE B BY AINV*B.
  99       DO 4 J=1,K
 100          B1=B(1,J)
 101          B2=B(2,J)
 102          B(1,J)=A(1,1)*B1+A(1,2)*B2+A(1,3)*B(3,J)
 103          B(2,J)=A(2,1)*B1+A(2,2)*B2+A(2,3)*B(3,J)
 104          B(3,J)=A(3,1)*B1+A(3,2)*B2+A(3,3)*B(3,J)
 105     4 CONTINUE
 106       RETURN
 107 C
 108     5 IF(N.LT.2) GO TO 7
 109 C
 110 C  N=2 CASE BY CRAMERS RULE.
 111 C
 112       DET=A(1,1)*A(2,2)-A(1,2)*A(2,1)
 113       IF(DET.EQ.0D0) GO TO 11
 114       S=1D0/DET
 115       C11   =S*A(2,2)
 116       A(1,2)=-S*A(1,2)
 117       A(2,1)=-S*A(2,1)
 118       A(2,2)=S*A(1,1)
 119       A(1,1)=C11
 120       DO 6 J=1,K
 121          B1=B(1,J)
 122          B(1,J)=C11*B1+A(1,2)*B(2,J)
 123          B(2,J)=A(2,1)*B1+A(2,2)*B(2,J)
 124     6 CONTINUE
 125       RETURN
 126 C
 127 C  N=1 CASE.
 128 C
 129     7 IF(A(1,1).EQ.0D0) GO TO 11
 130       A(1,1)=1D0/A(1,1)
 131       DO 8 J=1,K
 132          B(1,J)=A(1,1)*B(1,J)
 133     8 CONTINUE
 134       RETURN
 135 C
 136 C  N.GT.3 CASES.  FACTORIZE MATRIX, INVERT AND SOLVE SYSTEM.
 137 C
 138     9 CALL DFACT(N,A,IDIM,R,IFAIL,DET,JFAIL)
 139       IF(IFAIL.NE.0) RETURN
 140       CALL DFEQN(N,A,IDIM,R,K,B)
 141       CALL DFINV(N,A,IDIM,R)
 142       RETURN
 143 C
 144 C  ERROR EXITS.
 145 C
 146    10 IFAIL=+1
 147       CALL F010PR(NAME,N,IDIM,K,KPRNT)
 148       RETURN
 149 C
 150    11 IFAIL=-1
 151       RETURN
 152 C
 153       END