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