Télécharger la plaquette Cast3M


fr en

Télécharger gauscl.eso

Retour à la liste

Numérotation des lignes : 

gauscl

  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