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 :
© Cast3M 2003 - All rights reserved.
Disclaimer