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