Télécharger la plaquette Cast3M


fr en

Télécharger gaussj.eso

Retour à la liste

Numérotation des lignes : 

gaussj

  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.  
  92.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales