MINICERN/packlib/kernlib/kernnum/f010fort/ceqn.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 CEQN(N,A,IDIM,R,IFAIL,K,B)
  11       REAL R(N),T1,T2,T3
  12       COMPLEX A(IDIM,N),B(IDIM,K),ONE,DET,S,TEMP,
  13      $        B1,Y1,Y2,L11,L21,L22,L31,L32,L33,U12,U13,U23
  14       CHARACTER*6 NAME
  15       DATA NAME/'CEQN'/,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, AFTER WHICH A IS UNDEFINED.
  21 C
  22 C     (PARAMETERS AS FOR CEQINV.)
  23 C
  24 C     CALLS ... CFACT, CFEQN, F010PR, ABEND.
  25 C
  26 C     ******************************************************************
  27 C
  28 C  TEST FOR PARAMETER ERRORS.
  29 C
  30       IF((N.LT.1).OR.(N.GT.IDIM).OR.(K.LT.1)) GO TO 11
  31 C
  32 C  TEST FOR N.LE.3.
  33 C
  34       IF(N.GT.3) GO TO 10
  35       IFAIL=0
  36       IF(N.LT.3) GO TO 6
  37 C
  38 C  N=3 CASE.
  39 C
  40 C     FACTORIZE MATRIX A=L*U.
  41 C     (FIRST PIVOT SEARCH)
  42       T1=ABS(REAL(A(1,1)))+ABS(AIMAG(A(1,1)))
  43       T2=ABS(REAL(A(2,1)))+ABS(AIMAG(A(2,1)))
  44       T3=ABS(REAL(A(3,1)))+ABS(AIMAG(A(3,1)))
  45       IF(T1.GE.T2) GO TO 1
  46          IF(T3.GE.T2) GO TO 2
  47 C        (PIVOT IS A21)
  48             M1=2
  49             M2=1
  50             M3=3
  51             GO TO 3
  52     1 IF(T3.GE.T1) GO TO 2
  53 C     (PIVOT IS A11)
  54          M1=1
  55          M2=2
  56          M3=3
  57          GO TO 3
  58 C     (PIVOT IS A31)
  59     2    M1=3
  60          M2=2
  61          M3=1
  62     3 TEMP=A(M1,1)
  63       IF( REAL(TEMP).EQ.0. .AND. AIMAG(TEMP).EQ.0. ) GO TO 10
  64       L11=ONE/TEMP
  65       U12=L11*A(M1,2)
  66       U13=L11*A(M1,3)
  67       L22=A(M2,2)-A(M2,1)*U12
  68       L32=A(M3,2)-A(M3,1)*U12
  69 C     (SECOND PIVOT SEARCH)
  70       T2=ABS(REAL(L22))+ABS(AIMAG(L22))
  71       T3=ABS(REAL(L32))+ABS(AIMAG(L32))
  72       IF(T2.GE.T3) GO TO 4
  73          I=M2
  74          M2=M3
  75          M3=I
  76          TEMP=L22
  77          L22=L32
  78          L32=TEMP
  79     4 L21=A(M2,1)
  80       L31=A(M3,1)
  81       IF( REAL(L22).EQ.0. .AND. AIMAG(L22).EQ.0. ) GO TO 10
  82       L22=ONE/L22
  83       U23=L22*(A(M2,3)-L21*U13)
  84       TEMP=A(M3,3)-L31*U13-L32*U23
  85       IF( REAL(TEMP).EQ.0. .AND. AIMAG(TEMP).EQ.0. ) GO TO 10
  86       L33=ONE/TEMP
  87 C
  88 C     SOLVE L*Y=B AND U*X=Y.
  89       DO 5 J=1,K
  90          Y1=L11*B(M1,J)
  91          Y2=L22*(B(M2,J)-L21*Y1)
  92          B(3,J)=L33*(B(M3,J)-L31*Y1-L32*Y2)
  93          B(2,J)=Y2-U23*B(3,J)
  94          B(1,J)=Y1-U12*B(2,J)-U13*B(3,J)
  95     5 CONTINUE
  96       RETURN
  97 C
  98     6 IF(N.LT.2) GO TO 8
  99 C
 100 C  N=2 CASE BY CRAMERS RULE.
 101 C
 102       DET=A(1,1)*A(2,2)-A(1,2)*A(2,1)
 103       IF( REAL(DET).EQ.0. .AND. AIMAG(DET).EQ.0. ) GO TO 12
 104       S=ONE/DET
 105       DO 7 J=1,K
 106          B1=B(1,J)
 107          B(1,J)=S*(A(2,2)*B1-A(1,2)*B(2,J))
 108          B(2,J)=S*(-A(2,1)*B1+A(1,1)*B(2,J))
 109     7 CONTINUE
 110       RETURN
 111 C
 112 C  N=1 CASE.
 113 C
 114     8 IF( REAL(A(1,1)).EQ.0. .AND. AIMAG(A(1,1)).EQ.0. ) GO TO 12
 115       S=ONE/A(1,1)
 116       DO 9 J=1,K
 117          B(1,J)=S*B(1,J)
 118     9 CONTINUE
 119       RETURN
 120 C
 121 C  N.GT.3 CASES.  FACTORIZE MATRIX AND SOLVE SYSTEM.
 122 C
 123    10 CALL CFACT(N,A,IDIM,R,IFAIL,DET,JFAIL)
 124       IF(IFAIL.NE.0) RETURN
 125       CALL CFEQN(N,A,IDIM,R,K,B)
 126       RETURN
 127 C
 128 C  ERROR EXITS.
 129 C
 130    11 IFAIL=+1
 131       CALL F010PR(NAME,N,IDIM,K,KPRNT)
 132       RETURN
 133 C
 134    12 IFAIL=-1
 135       RETURN
 136 C
 137       END