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