The TELEMAC-MASCARET system  trunk
cgsqua.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE cgsqua
3 ! *****************
4 !
5  &(x,a,b,mesh, g,g0,p,k,h,ahpk,cfg,infogr)
6 !
7 !***********************************************************************
8 ! BIEF V6P1 21/08/2010
9 !***********************************************************************
10 !
11 !brief SOLVES THE LINEAR SYSTEM A X = B
12 !+ USING THE SQUARED CONJUGATE GRADIENT METHOD.
13 !code
14 !+ ALGORITHM:
15 !+
16 !+ |
17 !+ | INITIALISATION
18 !+ | ---------------
19 !+ |
20 !+ | 0 N
21 !+ | X VECTOR IN R , APPROXIMATION OF THE SOLUTION
22 !+ |
23 !+ | 0 0
24 !+ | G = A X - B
25 !+ |
26 !+ | K0 = P0 = G0
27 !+ |
28 !+ |
29 !+ | ITERATIONS
30 !+ | ----------
31 !+ |
32 !+ | M 0
33 !+ | M ( K , G )
34 !+ | RO = ------------
35 !+ | M 0
36 !+ | (A P , G )
37 !+ |
38 !+ | M M M M
39 !+ | H = K - RO * A P
40 !+ |
41 !+ | M+1 M M M M
42 !+ | G = G - RO * A ( H + K )
43 !+ |
44 !+ | M+1 M M M M
45 !+ | X = X - RO * ( H + K )
46 !+ |
47 !+ | 0 M+1
48 !+ | ( G , G )
49 !+ | BETA = ------------
50 !+ | 0 M
51 !+ | ( G , G )
52 !+ |
53 !+ | M+1 M+1 M M M M
54 !+ | P = G + 2*BETA * H + BETA**2 * P
55 !+ |
56 !+ | M+1 M+1 M M
57 !+ | K = G + BETA * H
58 !+ |
59 !
60 !code
61 !+-----------------------------------------------------------------------
62 !+ PRECONDITIONING
63 !+ (ALSO SEE SOLV01)
64 !+-----------------------------------------------------------------------
65 !+ PRECON VALUE I MEANING
66 !+-----------------------------------------------------------------------
67 !+ 0 I NO PRECONDITIONING
68 !+ 2 I DIAGONAL PRECONDITIONING WITH THE MATRIX
69 !+ I DIAGONAL
70 !+ 3 I DIAGONAL PRECONDITIONING WITH THE CONDENSED
71 !+ I MATRIX
72 !+ 5 I OTHER
73 !+-----------------------------------------------------------------------
74 !code
75 !+ MEANING OF IELM :
76 !+
77 !+ TYPE OF ELEMENT NUMBER OF POINTS CODED IN THIS SUBROUTINE
78 !+
79 !+ 11 : P1 TRIANGLE 3 YES
80 !+ 12 : P2 TRIANGLE 6 YES
81 !+ 13 : P1-ISO P1 TRIANGLE 6 YES
82 !+ 14 : P2 TRIANGLE 7
83 !+ 21 : Q1 QUADRILATERAL 4 YES
84 !+ 22 : Q2 QUADRILATERAL 8
85 !+ 24 : Q2 QUADRILATERAL 9
86 !+ 31 : P1 TETRAHEDRON 4 YES
87 !+ 32 : P2 TETRAHEDRON 10
88 !+ 41 : MITHRIDATE PRISMS 6 YES
89 !
90 !history J-M HERVOUET (LNH) ; F LEPEINTRE (LNH)
91 !+ 24/04/97
92 !+ V5P1
93 !+
94 !
95 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
96 !+ 13/07/2010
97 !+ V6P0
98 !+ Translation of French comments within the FORTRAN sources into
100 !
101 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
102 !+ 21/08/2010
103 !+ V6P0
104 !+ Creation of DOXYGEN tags for automated documentation and
105 !+ cross-referencing of the FORTRAN sources
106 !
107 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
108 !| A |-->| MATRIX OF THE SYSTEM
109 !| AHPK |<->| WORK STRUCTURE
110 !| B |-->| RIGHT-HAND SIDE OF SYSTEM
111 !| CFG |-->| CFG(1): STORAGE OF MATRIX
112 !| | | CFG(2): MATRIX VECTOR PRODUCT
114 !| G0 |<->| INITIAL GRADIENT
115 !| H |<->| WORK STRUCTURE
116 !| INFOGR |-->| IF YES, INFORMATION PRINTED
117 !| K |<->| WORK STRUCTURE
118 !| MESH |-->| MESH STRUCTURE
119 !| P |<->| WORK STRUCTURE
120 !| X |<--| INITIAL VALUE, THEN SOLUTION
121 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
122 !
123  USE bief, ex_cgsqua => cgsqua
124 !
126  IMPLICIT NONE
127 !
128 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
129 !
130  TYPE(bief_obj) , INTENT(INOUT) :: X,G,G0,P,K,H,AHPK
131  TYPE(bief_obj) , INTENT(IN ) :: B
132  TYPE(bief_obj) , INTENT(IN) :: A
133  TYPE(slvcfg) , INTENT(INOUT) :: CFG
134  LOGICAL , INTENT(IN) :: INFOGR
135  TYPE(bief_mesh) , INTENT(INOUT) :: MESH
136 !
137 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
138 !
139  DOUBLE PRECISION XL,RO,TESTL,RL,GMP1G0,BETA,GMG0,C
140 !
141  INTEGER M
142 !
143  LOGICAL RELAT
144 !
145  INTRINSIC sqrt
146 !
147 !-----------------------------------------------------------------------
148 !
149 ! COMPUTES THE NORM OF THE SECOND MEMBER
150 !
151  xl = p_dots(b,b,mesh)
152 !
153  IF(xl.LT.1.d0) THEN
154  xl = 1.d0
155  relat = .false.
156  ELSE
157  relat = .true.
158  ENDIF
159 !
160  m = 0
161 !
162 ! INITIALISES G : A X0 - B
163 !
164  CALL matrbl( 'X=AY ',g,a,x,c, mesh)
165 !
166  CALL os('X=X-Y ', x=g, y=b)
167 !
168 ! CHECKS THAT THE ACCURACY HAS NOT ALREADY BEEN REACHED
169 !
170  rl = p_dots(g,g,mesh)
171 !
172  IF(rl.LT.cfg%EPS**2*xl) THEN
173  testl = sqrt(rl/xl)
174  IF (infogr) THEN
175  IF(relat) THEN
176  WRITE(lu,101) m,testl
177  ELSE
178  WRITE(lu,201) m,testl
179  ENDIF
180  ENDIF
181  GOTO 1000
182  ENDIF
183 !
184 ! INITIALISES G0 , P , AND K
185 !
186  CALL os('X=Y ', x=g0, y=g)
187  CALL os('X=Y ', x=p , y=g)
188  CALL os('X=Y ', x=k , y=g)
189 !
190  m = 1
191 20 CONTINUE
192 !
193 ! COMPUTES AP (IN H, RECOMPUTED LATER)
194 !
195  CALL matrbl( 'X=AY ',h,a,p,c, mesh)
196 !
197 ! COMPUTES RO
198 !
199  ro = p_dots(k,g0,mesh) / p_dots(h,g0,mesh)
200 !
201 ! COMPUTES H+K (IN H, AP ALREADY BEING IN H)
202 !
203  CALL os('X=CX ', x=h, c=-ro)
204  CALL os('X=X+CY ', x=h, y=k, c=2.d0)
205 !
206 ! M+1 M M M
207 ! COMPUTES X = X - RO * (H+K) (H+K IN H)
208 !
209  CALL os('X=X+CY ', x=x, y=h, c=-ro)
210 !
211 ! COMPUTES A(H+K) (HERE H+K IN H, A(H+K) IN AHPK)
212 !
213  CALL matrbl( 'X=AY ',ahpk,a,h,c, mesh)
214 !
215 ! M 0
216 ! COMPUTES ( G , G )
217 !
218  gmg0 = p_dots(g,g0,mesh)
219 !
220 ! COMPUTES GM
221 !
222  CALL os('X=X+CY ', x=g, y=ahpk, c=-ro)
223 !
224  rl = p_dots(g,g,mesh)
225  IF (rl.GT.cfg%EPS**2*xl) THEN
226  IF (m.GE.cfg%NITMAX) THEN
227  testl=sqrt(rl/xl)
228  IF(relat) THEN
229  WRITE(lu,103) m,testl
230  ELSE
231  WRITE(lu,203) m,testl
232  ENDIF
233  GOTO 1000
234  ELSE
235  m = m + 1
236  ENDIF
237  ELSE
238  IF(infogr) THEN
239  testl=sqrt(rl/xl)
240  IF(relat) THEN
241  WRITE(lu,101) m,testl
242  ELSE
243  WRITE(lu,201) m,testl
244  ENDIF
245  ENDIF
246  GOTO 1000
247  ENDIF
248 !
249 ! M+1 0
250 ! COMPUTES ( G , G )
251 !
252  gmp1g0 = p_dots(g,g0,mesh)
253 !
254 ! COMPUTES BETA
255 !
256  beta = gmp1g0 / gmg0
257 !
258 ! M
259 ! COMPUTES H (H+K IN H)
260 !
261  CALL os('X=X-Y ', x=h, y=k)
262 !
263 ! M+1
264 ! COMPUTES P
265 !
266  CALL os('X=CX ', x=p, c=beta**2 )
267  CALL os('X=X+Y ', x=p, y=g)
268  CALL os('X=X+CY ', x=p, y=h, c=2*beta)
269 !
270 ! M+1
271 ! COMPUTES K
272 !
273  CALL os('X=Y ', x=k, y=g)
274  CALL os('X=X+CY ', x=k, y=h, c=beta)
275 !
276  GOTO 20
277 !
278 1000 RETURN
279 !
280 !-----------------------------------------------------------------------
281 !
282 ! FORMATS
283 !
284 101 FORMAT(1x,'CGSQUA (BIEF) : ',1i8,' ITERATIONS',
285  & ' RELATIVE PRECISION: ',g16.7)
286 201 FORMAT(1x,'CGSQUA (BIEF) : ',1i8,' ITERATIONS',
287  & ' ABSOLUTE PRECISION: ',g16.7)
288 103 FORMAT(1x,'CGSQUA (BIEF) : EXCEEDING MAXIMUM ITERATIONS ',1i8,
289  & ' RELATIVE PRECISION: ',g16.7)
290 203 FORMAT(1x,'CGSQUA (BIEF) : EXCEEDING MAXIMUM ITERATIONS ',1i8,
291  & ' ABSOLUTE PRECISION:',g16.7)
292 !
293 !-----------------------------------------------------------------------
294 !
295  END
subroutine cgsqua(X, A, B, MESH, G, G0, P, K, H, AHPK, CFG, INFOGR)
Definition: cgsqua.f:7
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
subroutine matrbl(OP, X, A, Y, C, MESH)
Definition: matrbl.f:7
double precision function p_dots(X, Y, MESH)
Definition: p_dots.f:7
Definition: bief.f:3