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

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/02/15 17:48:49  mclareni
   6 * Kernlib
   7 *
   8 *
   9 #include "kernnum/pilot.h"
  10       SUBROUTINE DINV(N,A,IDIM,R,IFAIL)
  11       REAL R(N),T1,T2,T3
  12       DOUBLE PRECISION A(IDIM,N),DET,TEMP,S,
  13      $                 C11,C12,C13,C21,C22,C23,C31,C32,C33
  14       CHARACTER*6 NAME
  15       DATA NAME/'DINV'/,KPRNT/0/
  16 C
  17 C     ******************************************************************
  18 C
  19 C     REPLACES A BY ITS INVERSE.
  20 C
  21 C     (PARAMETERS AS FOR DEQINV.)
  22 C
  23 C     CALLS ... DFACT, DFINV, F010PR, ABEND.
  24 C
  25 C     ******************************************************************
  26 C
  27 C  TEST FOR PARAMETER ERRORS.
  28 C
  29       IF((N.LT.1).OR.(N.GT.IDIM)) GO TO 7
  30 C
  31 C  TEST FOR N.LE.3.
  32 C
  33       IF(N.GT.3) GO TO 6
  34       IFAIL=0
  35       IF(N.LT.3) GO TO 4
  36 C
  37 C  N=3 CASE.
  38 C
  39 C     COMPUTE COFACTORS.
  40       C11=A(2,2)*A(3,3)-A(2,3)*A(3,2)
  41       C12=A(2,3)*A(3,1)-A(2,1)*A(3,3)
  42       C13=A(2,1)*A(3,2)-A(2,2)*A(3,1)
  43       C21=A(3,2)*A(1,3)-A(3,3)*A(1,2)
  44       C22=A(3,3)*A(1,1)-A(3,1)*A(1,3)
  45       C23=A(3,1)*A(1,2)-A(3,2)*A(1,1)
  46       C31=A(1,2)*A(2,3)-A(1,3)*A(2,2)
  47       C32=A(1,3)*A(2,1)-A(1,1)*A(2,3)
  48       C33=A(1,1)*A(2,2)-A(1,2)*A(2,1)
  49       T1=ABS(SNGL(A(1,1)))
  50       T2=ABS(SNGL(A(2,1)))
  51       T3=ABS(SNGL(A(3,1)))
  52 C
  53 C     (SET TEMP=PIVOT AND DET=PIVOT*DET.)
  54       IF(T1.GE.T2) GO TO 1
  55          IF(T3.GE.T2) GO TO 2
  56 C        (PIVOT IS A21)
  57             TEMP=A(2,1)
  58             DET=C13*C32-C12*C33
  59             GO TO 3
  60     1 IF(T3.GE.T1) GO TO 2
  61 C     (PIVOT IS A11)
  62          TEMP=A(1,1)
  63          DET=C22*C33-C23*C32
  64          GO TO 3
  65 C     (PIVOT IS A31)
  66     2    TEMP=A(3,1)
  67          DET=C23*C12-C22*C13
  68 C
  69 C     SET ELEMENTS OF INVERSE IN A.
  70     3 IF(DET.EQ.0D0) GO TO 8
  71       S=TEMP/DET
  72       A(1,1)=S*C11
  73       A(1,2)=S*C21
  74       A(1,3)=S*C31
  75       A(2,1)=S*C12
  76       A(2,2)=S*C22
  77       A(2,3)=S*C32
  78       A(3,1)=S*C13
  79       A(3,2)=S*C23
  80       A(3,3)=S*C33
  81       RETURN
  82 C
  83     4 IF(N.LT.2) GO TO 5
  84 C
  85 C  N=2 CASE BY CRAMERS RULE.
  86 C
  87       DET=A(1,1)*A(2,2)-A(1,2)*A(2,1)
  88       IF(DET.EQ.0D0) GO TO 8
  89       S=1D0/DET
  90       C11   =S*A(2,2)
  91       A(1,2)=-S*A(1,2)
  92       A(2,1)=-S*A(2,1)
  93       A(2,2)=S*A(1,1)
  94       A(1,1)=C11
  95       RETURN
  96 C
  97 C  N=1 CASE.
  98 C
  99     5 IF(A(1,1).EQ.0D0) GO TO 8
 100       A(1,1)=1D0/A(1,1)
 101       RETURN
 102 C
 103 C  N.GT.3 CASES.  FACTORIZE MATRIX AND INVERT.
 104 C
 105     6 CALL DFACT(N,A,IDIM,R,IFAIL,DET,JFAIL)
 106       IF(IFAIL.NE.0) RETURN
 107       CALL DFINV(N,A,IDIM,R)
 108       RETURN
 109 C
 110 C  ERROR EXITS.
 111 C
 112     7 IFAIL=+1
 113       CALL F010PR(NAME,N,IDIM,K,KPRNT)
 114       RETURN
 115 C
 116     8 IFAIL=-1
 117       RETURN
 118 C
 119       END