Télécharger gauscl.eso

Retour à la liste

Numérotation des lignes :

  1. C GAUSCL SOURCE CHAT 05/01/13 00:16:05 5004
  2. c ============================================================
  3. c
  4. c Resolution par la methode de GAUSS pivot total
  5. c A*X=B
  6. c A(dimini,dimini) Matrice
  7. c X(dimini) Vecteur inconnu
  8. c B(dimini) Second membre
  9. c Dima dimension du probleme
  10. c iprint = 0 pas de trace
  11. c >= 4 Test de qualite de la solution
  12. c >= 5 Affiche la matrice et le second membre
  13. c ============================================================
  14.  
  15. SUBROUTINE GAUSCL(A,X,B,DIMINI,DIMA,IPRINT)
  16. IMPLICIT INTEGER(I-N)
  17. integer dimmat,dima,dimini
  18. parameter(dimmat=10)
  19. REAL*8 a(1:dimini,1:dimini),x(1:dimini),b(1:dimini)
  20. REAL*8 mat(1:dimmat,1:dimmat+1),sol(1:dimmat)
  21. integer ind(1:dimmat)
  22. integer lig,col,i,j,k,l,iprint
  23. REAL*8 xmax,coeff
  24. character*80 fmt1
  25.  
  26. if ((dima.gt.dimmat).or.(dima.gt.dimini)) then
  27. write(*,*)'Pb dans les dimensions'
  28. stop
  29. endif
  30.  
  31. write(fmt1,1000) dima
  32. 1000 format('(2X,',I2,'(E13.5),6X,E13.5)')
  33.  
  34. do i=1,dima
  35. do j=1,dima
  36. mat(i,j)=a(i,j)
  37. enddo
  38. mat(i,dima+1)=b(i)
  39. ind(i)=i
  40. enddo
  41.  
  42. if (iprint.ge.5) then
  43. write(*,*)'Matrice Origine.'
  44. do i=1,dima
  45. write(*,fmt1) (mat(i,j),j=1,dima+1)
  46. enddo
  47. write(*,*)
  48. endif
  49.  
  50. c ... Triangularisation de Mat
  51. do l=1,dima-1
  52. lig=0
  53. col=0
  54. xmax=0.0
  55. c ..... recherche du plus grand pivot
  56. do i=l,dima
  57. do j=l,dima
  58. if (abs(mat(i,j)).gt.xmax) then
  59. xmax=abs(mat(i,j))
  60. lig=i
  61. col=j
  62. endif
  63. enddo
  64. enddo
  65. if (lig.eq.0) then
  66. write(*,*) 'Detection d''un pivot nul. DET = 0'
  67. stop 'Dans gauss'
  68. endif
  69. c ..... inversion des lignes
  70. if (lig.ne.l) then
  71. do j=l,dima+1
  72. xmax=mat(l,j)
  73. mat(l,j)=mat(lig,j)
  74. mat(lig,j)=xmax
  75. enddo
  76. endif
  77. c ..... inversion des colonnes
  78. if (col.ne.l) then
  79. do i=1,dima
  80. xmax=mat(i,l)
  81. mat(i,l)=mat(i,col)
  82. mat(i,col)=xmax
  83. enddo
  84. k=ind(l)
  85. ind(l)=ind(col)
  86. ind(col)=k
  87. endif
  88.  
  89. if (iprint.ge.6) then
  90. write(*,*)'Apres inversion des lig & col. Matrice No : ',l
  91. do i=1,dima
  92. write(*,fmt1) (mat(i,j),j=1,dima+1)
  93. enddo
  94. write(*,*)
  95. endif
  96.  
  97. c ..... diagonalisation
  98. do i=l+1,dima
  99. coeff=mat(i,l)/mat(l,l)
  100. mat(i,l)=0.0
  101. do j=l+1,dima+1
  102. mat(i,j)=mat(i,j)-mat(l,j)*coeff
  103. enddo
  104. enddo
  105. if (iprint.ge.5) then
  106. write(*,*)'Apres diagonalisation. Matrice No : ',l
  107. do i=1,dima
  108. write(*,fmt1) (mat(i,j),j=1,dima+1)
  109. enddo
  110. write(*,*)
  111. endif
  112. enddo
  113.  
  114. c ... resolution du systeme
  115. sol(dima)=mat(dima,dima+1)/mat(dima,dima)
  116. do l=dima-1,1,-1
  117. xmax=0.0
  118. do j=l+1,dima
  119. xmax=xmax+mat(l,j)*sol(j)
  120. enddo
  121. sol(l)=(mat(l,dima+1)-xmax)/mat(l,l)
  122. enddo
  123. c ... reclassement des colonnes
  124. do j=1,dima
  125. x(ind(j))=sol(j)
  126. enddo
  127.  
  128. if (iprint.ge.5) then
  129. write(*,*)'Solution de la matrice origine.'
  130. do i=1,dima
  131. write(*,*) x(i)
  132. enddo
  133. write(*,*)
  134. endif
  135.  
  136. c ... Test de qualite
  137. if (iprint.ge.4) then
  138. do i=1,dima
  139. xmax=0.0
  140. do j=1,dima
  141. xmax=xmax+a(i,j)*x(j)
  142. enddo
  143. write(*,*)' ligne No ',i,' . Ecart de : ',xmax-b(i)
  144. enddo
  145. endif
  146.  
  147. return
  148. end
  149.  
  150.  

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