Sources Esope | Cast3M

gaussk
C GAUSSK    SOURCE    PV        05/12/01    21:15:16     5252
      subroutine gaussk(a,n,np,b,m,mp)
c solution d un systeme lineaire par elimination de Gauss-Jordan
c A matrice carree N*N, concretement NP*NP.
c B matrice N*M contenant M termes de second membre, concretement NP*MP
c en sortie A est inversee et B est la collection des vecteurs solutions
      parameter (nmax = 50)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      real*8 a(np,np),b(np,mp)
      integer ipiv(nmax),indxr(nmax),indxc(nmax)
c suppose nmax > n
 
      do 11 j = 1,N
        ipiv(j) = 0
 11   continue
 
* boucle principale sur les colonnes
      do 22 i = 1,n
        big = 0.
        do 13 j = 1,n
           if (ipiv(j).ne.1) then
             do 12 k=1,n
                if (ipiv(k).eq.0) then
                  if (abs(a(j,k)).ge.big) then
                    big = abs(a(j,k))
                    irow = j
                    icol = k
                  endif
                else if (ipiv(k).gt.1) then
                  write(6,*) ' matrice singuliere '
                endif
 12          continue
           endif
 13     continue
        ipiv(icol) = ipiv(icol) + 1
 
* on fait passer le pivot sur la diagonale et on reindexe les colonnes
        if (irow.ne.icol) then
          do 14 L=1,N
            DUM = A(IROW,L)
            a(irow,l) = a(icol,l)
            a(icol,l) = dum
 14       continue
          do 15 l = 1,m
            dum = b(irow,l)
            b(irow,l) = b(icol,l)
            b(icol,l) = dum
 15       continue
        endif
* divise la rangee par le pivot
        indxr(i) = irow
        indxc(i) = icol
        if (a(icol,icol).eq.0) then
           write(6,*) 'matrice singuliere '
           return
        endif
        pivinv = 1./a(icol,icol)
        a(icol,icol) = 1.
        do 16 l = 1,n
          a(icol,l) = a(icol,l)*pivinv
 16     continue
        do 17 l = 1,m
          b(icol,l) = b(icol,l)*pivinv
 17     continue
* reduit ensuite les rangees
        do 21 ll=1,n
* a l exception des pivots
           if (ll.ne.icol) then
             dum = a(ll,icol)
             a(ll,icol) =0.
             do 18 l=1,n
               a(ll,l) = a(ll,l) - (a(icol,l)*dum)
 18          continue
             do 19 l= 1,m
               b(ll,l) = b(ll,l) - (b(icol,l)*dum)
 19          continue
           endif
 21     continue
 22   continue
* fin de la boucle principale
 
* remettre les colonnes dans l ordre par permutation inverse
      do 24 l=N,1,-1
        if(indxr(l).ne.indxc(l)) then
          do 23 k =1,n
            dum = a(k,indxr(l))
            a(k,indxr(l)) = a(k,indxc(l))
            a(k,indxc(l)) = dum
 23       continue
        endif
 24   continue
 
      return
      end