Sources Esope | Cast3M

invma2
C INVMA2    SOURCE    FD218221  23/03/09    21:15:04     11622          
      SUBROUTINE INVMA2(A,NZ,ICRIT)
C
C  ====================================================================
C        INVERSION DE MATRICE A PAR LA METHODE DES DETERMINANTS
C            matrice limitee a la dimension 4 georgin insa
C    A  MATRICE (NZ*NZ) A INVERSER EN ENTREE, MATRICE INVERSEE EN SORTIE
C    ICRIT=1 SI MATRICE SINGULIERE, 0 SINON
C
C  ====================================================================
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(4,4),O(4,4)
      CALL ZERO(O,4,4)
      ICRIT=0
C
C ************************** TRANSPOSITION ****************************
C
      A11=A(1,1)
      A12=A(1,2)
      A13=A(1,3)
      A14=A(1,4)
      A21=A(2,1)
      A22=A(2,2)
      A23=A(2,3)
      A24=A(2,4)
      A31=A(3,1)
      A32=A(3,2)
      A33=A(3,3)
      A34=A(3,4)
      A41=A(4,1)
      A42=A(4,2)
      A43=A(4,3)
      A44=A(4,4)
      CALL ZERO(A,4,4)
C
C ************************** TEST SUR LA DIMENSION ********************
C
      GOTO (1,200,300,400),NZ
    1 CONTINUE
      DELTA=A11
      O(1,1)=1.D0
      GOTO 500
  200 CONTINUE
      DELTA=A11*A22-A21*A12
      O(1,1)=A22
      O(1,2)=-A21
      O(2,1)=-A12
      O(2,2)=A11
      GOTO 500
  300 CONTINUE
      DELTA1=A11*A22*A33+A12*A23*A31+A21*A32*A13
      DELTA2=A31*A22*A13+A32*A23*A11+A21*A12*A33
      DELTA=DELTA1-DELTA2
      O(1,1)=A22*A33-A32*A23
      O(1,2)=A31*A23-A21*A33
      O(1,3)=A21*A32-A31*A22
      O(2,1)=A32*A13-A12*A33
      O(2,2)=A11*A33-A31*A13
      O(2,3)=A31*A12-A11*A32
      O(3,1)=A12*A23-A22*A13
      O(3,2)=A21*A13-A11*A23
      O(3,3)=A11*A22-A21*A12
      GOTO 500
  400 CONTINUE
       O(1,1)=(A22*A33*A44+A32*A43*A24+A42*A23*A34)
      O(1,1)=O(1,1)-(A42*A33*A24+A32*A23*A44+A43*A34*A22)
       O(1,2)=-(A21*A33*A44+A31*A43*A24+A41*A23*A34)
      O(1,2)=O(1,2)+(A41*A33*A24+A31*A23*A44+A43*A34*A21)
       O(1,3)=(A21*A32*A44+A31*A42*A24+A41*A22*A34)
      O(1,3)=O(1,3)-(A41*A32*A24+A42*A34*A21+A31*A22*A44)
       O(1,4)=-(A21*A32*A43+A22*A33*A41+A31*A42*A23)
      O(1,4)=O(1,4)+(A41*A32*A23+A42*A33*A21+A31*A22*A43)
       O(2,1)=-(A12*A33*A44+A32*A43*A14+A42*A13*A34)
      O(2,1)=O(2,1)+(A42*A33*A14+A32*A13*A44+A12*A43*A34)
       O(2,2)=(A11*A33*A44+A31*A43*A14+A41*A13*A34)
      O(2,2)=O(2,2)-(A41*A33*A14+A31*A13*A44+A11*A43*A34)
       O(2,3)=-(A11*A32*A44+A31*A42*A14+A41*A12*A34)
      O(2,3)=O(2,3)+(A41*A32*A14+A42*A34*A11+A31*A12*A44)
       O(2,4)=(A11*A32*A43+A31*A42*A13+A41*A12*A33)
      O(2,4)=O(2,4)-(A41*A32*A13+A31*A12*A43+A42*A33*A11)
       O(3,1)=(A12*A23*A44+A22*A43*A14+A13*A24*A42)
      O(3,1)=O(3,1)-(A42*A23*A14+A22*A13*A44+A43*A24*A12)
       O(3,2)=-(A11*A23*A44+A21*A43*A14+A41*A13*A24)
      O(3,2)=O(3,2)+(A41*A23*A14+A43*A24*A11+A21*A13*A44)
       O(3,3)=(A11*A22*A44+A21*A42*A14+A41*A12*A24)
      O(3,3)=O(3,3)-(A41*A22*A14+A21*A12*A44+A11*A42*A24)
       O(3,4)=-(A11*A22*A43+A21*A42*A13+A41*A12*A23)
      O(3,4)=O(3,4)+(A41*A22*A13+A21*A12*A43+A11*A42*A23)
       O(4,1)=-(A12*A23*A34+A22*A33*A14+A13*A24*A32)
      O(4,1)=O(4,1)+(A32*A23*A14+A22*A13*A34+A12*A33*A24)
       O(4,2)=(A11*A23*A34+A21*A33*A14+A31*A13*A24)
      O(4,2)=O(4,2)-(A31*A23*A14+A21*A13*A34+A33*A24*A11)
       O(4,3)=-(A11*A22*A34+A12*A24*A31+A21*A32*A14)
      O(4,3)=O(4,3)+(A31*A22*A14+A32*A24*A11+A21*A12*A34)
       O(4,4)=(A11*A22*A33+A21*A32*A13+A31*A12*A23)
      O(4,4)=O(4,4)-(A31*A22*A13+A32*A23*A11+A21*A12*A33)
      DELTA=A11*O(1,1)+A12*O(1,2)+A13*O(1,3)+A14*O(1,4)
      DELTA2=A21*O(2,1)+A22*O(2,2)+A23*O(2,3)+A14*O(2,4)
      GOTO 500
  500 CONTINUE
      IF (ABS(DELTA).LT.1.D-50) THEN
       ICRIT=1
       RETURN
      ENDIF
      DO I=1,NZ
       DO J=1,NZ
        A(I,J)=O(I,J)/DELTA
       ENDDO
      ENDDO
      RETURN
      END