C GAUSCL    SOURCE    CHAT      05/01/13    00:16:05     5004
c ============================================================
c
c       Resolution par la methode de GAUSS pivot total
c                     A*X=B
c       A(dimini,dimini)  Matrice
c       X(dimini)       Vecteur inconnu
c       B(dimini)       Second membre
c       Dima             dimension du probleme
c       iprint        =  0 pas de trace
c                     >= 4 Test de qualite de la solution
c                     >= 5 Affiche la matrice et le second membre
c ============================================================

      SUBROUTINE GAUSCL(A,X,B,DIMINI,DIMA,IPRINT)
      IMPLICIT INTEGER(I-N)
      integer dimmat,dima,dimini
      parameter(dimmat=10)
      REAL*8 a(1:dimini,1:dimini),x(1:dimini),b(1:dimini)
      REAL*8 mat(1:dimmat,1:dimmat+1),sol(1:dimmat)
      integer ind(1:dimmat)
      integer lig,col,i,j,k,l,iprint
      REAL*8 xmax,coeff
      character*80 fmt1

      if ((dima.gt.dimmat).or.(dima.gt.dimini)) then
        write(*,*)'Pb dans les dimensions'
        stop
      endif

      write(fmt1,1000) dima
1000  format('(2X,',I2,'(E13.5),6X,E13.5)')

      do i=1,dima
        do j=1,dima
          mat(i,j)=a(i,j)
        enddo
        mat(i,dima+1)=b(i)
        ind(i)=i
      enddo

      if (iprint.ge.5) then
        write(*,*)'Matrice Origine.'
        do i=1,dima
          write(*,fmt1) (mat(i,j),j=1,dima+1)
        enddo
        write(*,*)
      endif

c ... Triangularisation de Mat
      do l=1,dima-1
        lig=0
        col=0
        xmax=0.0
c ..... recherche du plus grand pivot
        do i=l,dima
          do j=l,dima
            if (abs(mat(i,j)).gt.xmax) then
              xmax=abs(mat(i,j))
              lig=i
              col=j
            endif
          enddo
        enddo
        if (lig.eq.0) then
          write(*,*) 'Detection d''un pivot nul. DET = 0'
          stop 'Dans gauss'
        endif
c ..... inversion des lignes
        if (lig.ne.l) then
          do j=l,dima+1
            xmax=mat(l,j)
            mat(l,j)=mat(lig,j)
            mat(lig,j)=xmax
          enddo
        endif
c ..... inversion des colonnes
        if (col.ne.l) then
          do i=1,dima
            xmax=mat(i,l)
            mat(i,l)=mat(i,col)
            mat(i,col)=xmax
          enddo
          k=ind(l)
          ind(l)=ind(col)
          ind(col)=k
        endif

        if (iprint.ge.6) then
          write(*,*)'Apres inversion des lig & col. Matrice No : ',l
          do i=1,dima
            write(*,fmt1) (mat(i,j),j=1,dima+1)
          enddo
          write(*,*)
        endif

c ..... diagonalisation
        do i=l+1,dima
          coeff=mat(i,l)/mat(l,l)
          mat(i,l)=0.0
          do j=l+1,dima+1
            mat(i,j)=mat(i,j)-mat(l,j)*coeff
          enddo
        enddo
        if (iprint.ge.5) then
          write(*,*)'Apres diagonalisation. Matrice No : ',l
          do i=1,dima
            write(*,fmt1) (mat(i,j),j=1,dima+1)
          enddo
          write(*,*)
        endif
      enddo

c ... resolution du systeme
      sol(dima)=mat(dima,dima+1)/mat(dima,dima)
      do l=dima-1,1,-1
        xmax=0.0
        do j=l+1,dima
          xmax=xmax+mat(l,j)*sol(j)
        enddo
        sol(l)=(mat(l,dima+1)-xmax)/mat(l,l)
      enddo
c ... reclassement des colonnes
      do j=1,dima
        x(ind(j))=sol(j)
      enddo

      if (iprint.ge.5) then
        write(*,*)'Solution de la matrice origine.'
        do i=1,dima
          write(*,*) x(i)
        enddo
        write(*,*)
      endif

c ... Test de qualite
      if (iprint.ge.4) then
        do i=1,dima
          xmax=0.0
          do j=1,dima
            xmax=xmax+a(i,j)*x(j)
          enddo
          write(*,*)'    ligne No ',i,' . Ecart de : ',xmax-b(i)
        enddo
      endif

      return
      end