Sources Esope | Cast3M

gauss
C GAUSS     SOURCE    CHAT      05/01/13    00:16:08     5004
c     Sous-Programme du module MISTRAL0
C     --------------------------------------------------------------------------
      SUBROUTINE GAUSS (A,B,N,NSY,X,ISOL)
C     --------------------------------------------------------------------------
C     Resolution de NSY systemes de N equations lineaires a N inconnues :
C                                  A.X = B
C     de meme matrice carree A et de matrices colonnes B differentes.
C     ..........................................................................
C     Pour tout ISY compris entre 1 et NSY inclus,
C          tout  J  compris entre 1 et  N  inclus :
C            N
C          SOMME  A(I,J) * X(I,ISY) = B(J,ISY)
C           I=1
C     ..........................................................................
C     Methode de Gauss :
C     (version 2.2 FORTRAN du 17/03/1998)
C     transformation en matrice triangulaire inferieure (virtuellement :
C     permutation interne des indices uniquement).
C
C     Pour tout  M  compris entre 1 et  N  inclus,
C     il existe (au debut de l'etape M) :
C     - une permutation IM de {1,N},
C     - une matrice carree AM (stockee dans A),
C     - des matrices colonnes BM (stockees dans B),
C     telles que, pour tout ISY compris entre 1 et NSY inclus :
C     - pour tout J compris entre 1 et J0=N-M+1 inclus :
C        J0
C       SOMME AM(IM(K),J) * X(IM(K),ISY) = BM(J,ISY)
C        K=1
C     - pour tout J compris entre J0+1 et N inclus :
C                                   J-1
C       X(IM(J),ISY) = BM(J,ISY) - SOMME AM(IM(K),J) * X(IM(K),ISY)
C                                   K=1
C
C     Les matrices A et B sont donc modifiees par GAUSS.
C     --------------------------------------------------------------------------
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H, O-Z)
c     std :
      PARAMETER ( NMAX = 7 )
      DIMENSION A(1:N,1:*),B(1:N,1:*),IM(1:NMAX),X(1:N,1:*),ISOL(1:*)
c     std.
C
C     Initialisation de la permutation
C
      DO K = 1,N
       IM(K) = K
      END DO
C     ------------------------------------------------------------------
C     Debut de la boucle sur M de reduction du nombre d'equations
C     -----------------------------------------------------------
      DO M = 1,N-1
       J0M1 = N-M
       J0 = J0M1+1
C     ----------------------------------
C     Recherche du pivot sur la ligne J0
C
       K0 = J0
       AI0J0 = A(IM(K0),J0)
       DO K = J0M1,1,-1
        AIJ0 = A(IM(K),J0)
        IF (ABS(AIJ0).GT.ABS(AI0J0)) THEN
         K0 = K
         AI0J0 = AIJ0
        END IF
       END DO
       I0 = IM(K0)
C
C     Systemes singuliers
C
       IF (AI0J0.EQ.0.) THEN
        GOTO 1
       END IF
C     ------------------------------------
C     Permutation eventuelle des inconnues
C
       IF (K0.LT.J0) THEN
        DO K = K0,J0M1
         IM(K) = IM(K+1)
        END DO
        IM(J0) = I0
       END IF
C     -----------------------------
C     Coefficients de X(IM(J0),ISY)
C
       DO K = 1,J0M1
        I = IM(K)
        A(I,J0) = A(I,J0)/AI0J0
       END DO
       DO ISY = 1,NSY
        B(J0,ISY) = B(J0,ISY)/AI0J0
       END DO
C     -----------------------------------------
C     Etablissement des N-M equations restantes
C
       DO J = 1,J0M1
        AI0J = A(I0,J)
        DO K = 1,J0M1
         I = IM(K)
         A(I,J) = A(I,J)-AI0J*A(I,J0)
        END DO
        DO ISY = 1,NSY
         B(J,ISY) = B(J,ISY)-AI0J*B(J0,ISY)
        END DO
       END DO
C     ------
      END DO
C     ---------------------------------------------------------
C     Fin de la boucle sur M de reduction du nombre d'equations
C     ------------------------------------------------------------------
      IM1 = IM(1)
      AIM11 = A(IM1,1)
C
C     Systemes singuliers
C
      IF (AIM11.EQ.0.) THEN
       J0 = 1
       GOTO 1
      END IF
C     ------------------------------------------------------------------
C     Calculs des solutions X(IM(J),ISY)
C
      DO ISY = 1,NSY
       X(IM1,ISY) = B(1,ISY)/AIM11
       DO J = 2,N
        XINJ = B(J,ISY)
        DO K = 1,J-1
         I = IM(K)
         XINJ = XINJ-A(I,J)*X(I,ISY)
        END DO
        X(IM(J),ISY) = XINJ
       END DO
       ISOL(ISY) = 1
      END DO
      RETURN
C     ------------------------------------------------------------------
C     Systemes singuliers
C
1     DO ISY = 1,NSY
       IF (B(J0,ISY).EQ.0.) THEN
C       Infinite de solutions
        ISOL(ISY) = 2
       ELSE
C       Aucune solution
        ISOL(ISY) = 0
       END IF
      END DO
        RETURN
      END