Esope source files

   1 : C GAUSSJ    SOURCE    PV        08/09/23    21:15:04     6164           
   2 :       subroutine gaussj(a,n,np,b,m,mp)
   3 : c 
   4 : c  resolution d'un systeme lineaire par le pivot de gauss
   5 : c  voir numerical recipee en fortran
   6 : c  en fortran pur
   7 : c  a est la matrice
   8 : c  b contient les econd membre en entree 
   9 : c
  10 :       IMPLICIT INTEGER(I-N)
  11 :       IMPLICIT REAL*8(A-H,O-Z)
  12 :       parameter (nmax=50)
  13 :       dimension a(np,*),b(np,*),ipiv(nmax),indxr(nmax),indxc(nmax)
  14 : *
  15 : *      write (*,*) 'Entrée dans GAUSSJ.'
  16 :       do 11 j=1,n
  17 :         ipiv(j)=0
  18 : 11    continue
  19 :       do 22 i=1,n
  20 :         big=0.d0
  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,*)  'singular matrix in gaussj.eso'
  32 :               endif
  33 : 12          continue
  34 :           endif
  35 : 13      continue
  36 :         ipiv(icol)=ipiv(icol)+1
  37 :         if (irow.ne.icol) then
  38 :           do 14 l=1,n
  39 :             dum=a(irow,l)
  40 :             a(irow,l)=a(icol,l)
  41 :             a(icol,l)=dum
  42 : 14        continue
  43 :           do 15 l=1,m
  44 :             dum=b(irow,l)
  45 :             b(irow,l)=b(icol,l)
  46 :             b(icol,l)=dum
  47 : 15        continue
  48 :         endif
  49 :         indxr(i)=irow
  50 :         indxc(i)=icol
  51 :         if (a(icol,icol).eq.0.) a(icol,icol)=1.d50
  52 :         pivinv=1.d0/a(icol,icol)
  53 :         a(icol,icol)=1.d0
  54 :         do 16 l=1,n
  55 :           a(icol,l)=a(icol,l)*pivinv
  56 : 16      continue
  57 :         do 17 l=1,m
  58 :           b(icol,l)=b(icol,l)*pivinv
  59 : 17      continue
  60 :         do 21 ll=1,n
  61 :           if(ll.ne.icol)then
  62 :             dum=a(ll,icol)
  63 :             a(ll,icol)=0.d0
  64 :             do 18 l=1,n
  65 :               a(ll,l)=a(ll,l)-a(icol,l)*dum
  66 : 18          continue
  67 :             do 19 l=1,m
  68 :               b(ll,l)=b(ll,l)-b(icol,l)*dum
  69 : 19          continue
  70 :           endif
  71 : 21      continue
  72 : 22    continue
  73 :       do 24 l=n,1,-1
  74 :         if(indxr(l).ne.indxc(l))then
  75 :           do 23 k=1,n
  76 :             dum=a(k,indxr(l))
  77 :             a(k,indxr(l))=a(k,indxc(l))
  78 :             a(k,indxc(l))=dum
  79 : 23        continue
  80 :         endif
  81 : 24    continue    
  82 : 
  83 : *      write (*,*) 'Sortie de GAUSSJ.'
  84 :       return
  85 :       end
  86 : 
  87 : 
  88 : 
  89 :  
  90 :  
  91 :