Télécharger gaussk.eso

Retour à la liste

Numérotation des lignes :

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

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