Esope source files

   1 : C GAUSSK    SOURCE    PV        05/12/01    21:15:16     5252           
   2 :       subroutine gaussk(a,n,np,b,m,mp)
   3 : c solution d un systeme lineaire par elimination de Gauss-Jordan
   4 : c A matrice carree N*N, concretement NP*NP. 
   5 : c B matrice N*M contenant M termes de second membre, concretement NP*MP
   6 : c en sortie A est inversee et B est la collection des vecteurs solutions
   7 :       parameter (nmax = 50)
   8 :       IMPLICIT INTEGER(I-N)
   9 :       IMPLICIT REAL*8(A-H,O-Z)
  10 :       real*8 a(np,np),b(np,mp)
  11 :       integer ipiv(nmax),indxr(nmax),indxc(nmax)
  12 : c suppose nmax > n
  13 : 
  14 :       do 11 j = 1,N
  15 :         ipiv(j) = 0
  16 :  11   continue
  17 : 
  18 : * boucle principale sur les colonnes
  19 :       do 22 i = 1,n
  20 :         big = 0.
  21 :         do 13 j = 1,n
  22 :            if (ipiv(j).ne.1) then
  23 :              do 12 k=1,n
  24 :                 if (ipiv(k).eq.0) then
  25 :                   if (abs(a(j,k)).ge.big) then
  26 :                     big = abs(a(j,k))
  27 :                     irow = j
  28 :                     icol = k
  29 :                   endif
  30 :                 else if (ipiv(k).gt.1) then
  31 :                   write(6,*) ' matrice singuliere ' 
  32 :                 endif
  33 :  12          continue
  34 :            endif
  35 :  13     continue
  36 :         ipiv(icol) = ipiv(icol) + 1
  37 : 
  38 : * on fait passer le pivot sur la diagonale et on reindexe les colonnes
  39 :         if (irow.ne.icol) then
  40 :           do 14 L=1,N
  41 :             DUM = A(IROW,L)
  42 :             a(irow,l) = a(icol,l)
  43 :             a(icol,l) = dum
  44 :  14       continue
  45 :           do 15 l = 1,m
  46 :             dum = b(irow,l)
  47 :             b(irow,l) = b(icol,l)
  48 :             b(icol,l) = dum
  49 :  15       continue
  50 :         endif
  51 : * divise la rangee par le pivot
  52 :         indxr(i) = irow
  53 :         indxc(i) = icol
  54 :         if (a(icol,icol).eq.0) then
  55 :            write(6,*) 'matrice singuliere '
  56 :            return
  57 :         endif 
  58 :         pivinv = 1./a(icol,icol)
  59 :         a(icol,icol) = 1.
  60 :         do 16 l = 1,n
  61 :           a(icol,l) = a(icol,l)*pivinv
  62 :  16     continue
  63 :         do 17 l = 1,m
  64 :           b(icol,l) = b(icol,l)*pivinv
  65 :  17     continue
  66 : * reduit ensuite les rangees
  67 :         do 21 ll=1,n
  68 : * a l exception des pivots
  69 :            if (ll.ne.icol) then
  70 :              dum = a(ll,icol)
  71 :              a(ll,icol) =0.
  72 :              do 18 l=1,n
  73 :                a(ll,l) = a(ll,l) - (a(icol,l)*dum)
  74 :  18          continue
  75 :              do 19 l= 1,m
  76 :                b(ll,l) = b(ll,l) - (b(icol,l)*dum)
  77 :  19          continue
  78 :            endif
  79 :  21     continue
  80 :  22   continue
  81 : * fin de la boucle principale
  82 : 
  83 : * remettre les colonnes dans l ordre par permutation inverse
  84 :       do 24 l=N,1,-1
  85 :         if(indxr(l).ne.indxc(l)) then
  86 :           do 23 k =1,n
  87 :             dum = a(k,indxr(l))
  88 :             a(k,indxr(l)) = a(k,indxc(l))
  89 :             a(k,indxc(l)) = dum
  90 :  23       continue
  91 :         endif
  92 :  24   continue
  93 : 
  94 :       return
  95 :       end
  96 :  
  97 :  
  98 :