The TELEMAC-MASCARET system  trunk
gmres.f
Go to the documentation of this file.
1 ! ****************
2  SUBROUTINE gmres
3 ! ****************
4 !
5  & (x,a,b,mesh,r0,v,av,cfg,infogr,aux)
6 !
7 !***********************************************************************
8 ! BIEF V7P1
9 !***********************************************************************
10 !
11 !brief SOLVES A LINEAR SYSTEM A X = B
12 !+ USING THE GMRES (GENERALISED MINIMUM RESIDUAL) METHOD.
13 !code
14 !+ ALGORITHM:
15 !+
16 !+ |
17 !+ | INITIALISATION
18 !+ | ---------------
19 !+ |
20 !+ | R0 = B - A X0
21 !+ |
22 !+ | V1 = R0 / IIR0II
23 !+ |
24 !+ | K = DIMENSION OF THE KRYLOV SPACE
25 !+ |
26 !+ |
27 !+ |
28 !+ |
29 !+ | ITERATIONS
30 !+ | ----------
31 !+ |
32 !+ | 1) BUILDS AN ORTHONORMAL BASE (VI) OF DIMENSION K+1
33 !+ |
34 !+ | I = 1,...,K
35 !+ |
36 !+ | VI+1 = A * VI
37 !+ |
38 !+ | J = 1,...,I
39 !+ |
40 !+ | HI+1,J = ( VI+1 , VJ ) HESSENBERG MATRIX (K+1,K)
41 !+ |
42 !+ | VI+1 <--- VI+1 - HI+1,J * VJ
43 !+ |
44 !+ | VI+1 = VI+1 / IIVI+1II
45 !+ |
46 !+ |
47 !+ | 2) MULTIPLIES H BY Q = RK * RK-1 * ... * R1
48 !+ |
49 !+ | GIVENS' ROTATION MATRIX RI :
50 !+ | - -
51 !+ | I ID(J-1) I
52 !+ | I I
53 !+ | I CJ SJ I
54 !+ | RI = I -SJ CJ I
55 !+ | I I
56 !+ | I ID(K-J) I
57 !+ | - -
58 !+ |
59 !+ | 2 2
60 !+ | WITH CJ + SJ = 1
61 !+ |
62 !+ | CJ AND SJ SUCH THAT H BECOMES TRIANGULAR
63 !+ |
64 !+ | ID(J-1) : IDENTITY MATRIX J-1*J-1
65 !+ |
66 !+ |
67 !+ | 3) SOLVES THE SYSTEM H Y = Q E
68 !+ | - -
69 !+ | E VECTOR (NR0,0,0,0,0,0,.....)
70 !+ |
71 !+ | NR0 NORM OF THE RESIDUAL
72 !+ |
73 !+ |
74 !+ |
75 !+ | 4) COMPUTES X(M+1) = X(M) + V * Y
76 !+ |
77 !+ | V : MATRIX (3*NPOIN,K) WHICH COLUMNS ARE THE VJ
78 !+ |
79 !+ | 5) CHECKS THE RESIDUAL...
80 !+ |
81 !
82 !warning CROUT AND GSEBE PRECONDITIONING ARE TREATED HERE LIKE A
83 !+ LU PRECONDITIONING. FOR CROUT IT MEANS THAT THE DIAGONAL
84 !+ OF THE PRECONDITIONING MATRIX IS NOT TAKEN INTO ACCOUNT
85 !code
86 !+-----------------------------------------------------------------------
87 !+ PRECONDITIONING
88 !+-----------------------------------------------------------------------
89 !+ PRECON VALUE I MEANING
90 !+-----------------------------------------------------------------------
91 !+ 0 OR 1 I NO PRECONDITIONING
92 !+ I
93 !+ 2 I DIAGONAL PRECONDITIONING USING THE MATRIX
94 !+ I DIAGONAL
95 !+ I
96 !+ 3 I BLOCK-DIAGONAL PRECONDITIONING
97 !+ I
98 !+ 5 I DIAGONAL PRECONDITIONING USING THE ABSOLUTE
99 !+ I VALUE OF THE MATRIX DIAGONAL
100 !+ I
101 !+ 7 I CROUT EBE PRECONDITIONING
102 !+ I
103 !+ 11 I GAUSS-SEIDEL EBE PRECONDITIONING
104 !+ I
105 !+-----------------------------------------------------------------------
106 !
107 !history J-M HERVOUET & C. MOULIN (LNH)
108 !+ 26/08/1993
109 !+ FIRST IMPLEMENTATION
110 !+ First version.
111 !
112 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
113 !+ 13/07/2010
114 !+ V6P0
115 !+ Translation of French comments within the FORTRAN sources into
117 !
118 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
119 !+ 21/08/2010
120 !+ V6P0
121 !+ Creation of DOXYGEN tags for automated documentation and
122 !+ cross-referencing of the FORTRAN sources
123 !
124 !history J-M HERVOUET (EDF LAB, LNHE)
125 !+ 12/06/2014
126 !+ V7P0
127 !+ The case A=diagonal caused a crash, because in that case the Krylov
128 !+ space cannot be a basis. This case is now treated since the
129 !+ solution is easily found : X=A%D/B
130 !
131 !history J-M HERVOUET (EDF LAB, LNHE)
132 !+ 23/07/2015
133 !+ V7P1
134 !+ Now when the algorithm fails the dimension of the Krylov space
135 !+ is reduced. This is more general than the previous solution
136 !+ that assumed the matrix diagonal. Various cases with diagonal
137 !+ matrix also tested.
138 !
139 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
140 !| A |-->| MATRIX OF THE SYSTEM
141 !| AUX |-->| MATRIX FOR PRECONDITIONING.
142 !| AV |<->| WORK ARRAY
143 !| B |-->| RIGHT-HAND SIDE OF THE SYSTEM
144 !| CFG |-->| STRUCTURE OF SOLVER CONFIGURATION
145 !| INFOGR |-->| IF YES, PRINT A LOG.
146 !| MESH |-->| MESH STRUCTURE.
147 !| R |<->| RESIDUAL
148 !| V |<->| WORK ARRAY
149 !| X |<->| INITIAL VALUE, THEN SOLUTION
150 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
151 !
152  USE bief, ex_gmres => gmres
153 !
155  IMPLICIT NONE
156 !
157 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
158 !
159  TYPE(slvcfg), INTENT(INOUT) :: CFG
160  TYPE(bief_obj), INTENT(INOUT) :: B
161  TYPE(bief_obj), INTENT(INOUT) :: X,V,AV,R0
162  TYPE(bief_mesh), INTENT(INOUT) :: MESH
163  TYPE(bief_obj), INTENT(IN) :: A
164  TYPE(bief_obj), INTENT(INOUT) :: AUX
165  LOGICAL, INTENT(IN) :: INFOGR
166 !
167 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
168 !
169 ! HARDCODED MAXIMUM CFG%KRYLOV=20
170 !
171  DOUBLE PRECISION H(21,20),C(20),S(20),E(21),BID
172 !
173  DOUBLE PRECISION R,ZZ,NB,PREC,NR0
174 !
175  INTEGER I,J,K,L,M
176 !
177  LOGICAL RELAT,CROUT,GSEB,PRECO
178 !
179  INTRINSIC sqrt,abs
180 !
181 !-----------------------------------------------------------------------
182 !
183  crout=.false.
184  IF(mod(cfg%PRECON,7).EQ.0) crout=.true.
185  gseb=.false.
186  IF(mod(cfg%PRECON,11).EQ.0) gseb=.true.
187  preco=.false.
188  IF(crout.OR.gseb.OR.mod(cfg%PRECON,13).EQ.0) preco=.true.
189 !
190  IF(preco) THEN
191 ! -1
192 ! COMPUTES L B
193  CALL godown(b, aux,b,'D',mesh,.false.)
194  ENDIF
195 !
196 ! INITIALISES
197 !
198  nb = p_dots(b,b,mesh)
199  nb = sqrt(nb)
200 !
201  relat = .true.
202  IF(nb.LT.1.d0) THEN
203  nb = 1.d0
204  relat = .false.
205  ENDIF
206 !
207  m = 0
208 !
209 ! INITIALISES THE RESIDUAL R : A X0 - B
210 !
211  CALL matrbl( 'X=AY ',r0,a,x,bid,mesh)
212 !
213  IF(preco) THEN
214  CALL godown(r0, aux,r0,'D',mesh,.false.)
215  CALL puog(x,aux,x, 'D',mesh,.false.)
216  ENDIF
217 !
218  CALL os('X=X-Y ', x=r0, y=b)
219 !
220 ! CHECKS THAT THE ACCURACY IS NOT ALREADY REACHED
221 !
222  nr0=p_dots(r0,r0,mesh)
223  nr0=sqrt(nr0)
224  prec = nr0/nb
225 !
226  IF(prec.LE.cfg%EPS) GO TO 3000
227 !
228 !-----------------------------------------------------------------------
229 ! ITERATIONS LOOP
230 !-----------------------------------------------------------------------
231 !
232 20 CONTINUE
233 !
234  m = m+1
235 !
236 ! K CAN BE REDUCED DURING THE ITERATION IF THE ALGORITHM FAILS
237 !
238  k = cfg%KRYLOV
239 !
240 ! COMPUTES THE VECTOR V1 = - R0 / NORM(R0)
241 ! (- SIGN BECAUSE R = AX - B INSTEAD OF B - AX)
242 !
243  CALL os('X=CY ', x=v%ADR(1)%P, y=r0, c=-1.d0/nr0)
244 !
245 !-----------------------------------------------------------------------
246 ! COMPUTES THE ORTHONORMAL BASE AND MATRIX H
247 !-----------------------------------------------------------------------
248 !
249 ! K-1 1ST COLUMNS
250 !
251  DO j=1,k-1
252 !
253  IF(preco) THEN
254  CALL goup(b , aux , v%ADR(j)%P ,'D',mesh,.true.)
257  ELSE
259  ENDIF
260 !
261 ! NEXT VECTOR IN THE BASIS : V%ADR(J+1)%P
262 !
264 !
265 ! ORTHOGONALISATION PROCESS
266 !
267  DO i = 1,j
270  ENDDO
271 !
273 !
274  IF(h(j+1,j).LT.1.d-20) THEN
275  IF(infogr) THEN
276  WRITE(lu,*) 'GMRES: NORM OF VECTOR ',j+1
277  WRITE(lu,*) 'IN THE BASIS EQUALS ',h(j+1,j)
278  WRITE(lu,*) 'SIZE OF KRYLOV SPACE REDUCED TO ',j
279  ENDIF
280 ! REDUCTION OF KRYLOV SPACE
281  k=j
282  GO TO 2000
283  ENDIF
284  CALL os('X=CX ', x=v%ADR(j+1)%P, c=1.d0/h(j+1,j))
285 !
286  ENDDO ! J
287 !
288 ! K-TH COLUMN (VECTOR V(K+1) IS NOT COMPLETELY BUILT)
289 !
290 2000 CONTINUE
291 !
292  IF(preco) THEN
293  CALL goup(b , aux , v%ADR(k)%P , 'D' ,mesh,.true.)
296  ELSE
298  ENDIF
299 !
301 !
302  DO i = 1,k
304  h(k+1,k) = h(k+1,k) - h(i,k)**2
305  ENDDO
306 !
307 ! IN THEORY H(K+1,K) IS POSITIVE
308 ! TO MACHINE ACCURACY
309  h(k+1,k) = sqrt( max(h(k+1,k),0.d0) )
310 !
311 !-----------------------------------------------------------------------
312 ! BUILDS GIVENS' ROTATIONS AND APPLIES TO H AND E
313 !-----------------------------------------------------------------------
314 !
315 ! OTHER COMPONENTS FOR E ARE 0
316  e(1) = nr0
317 !
318 ! ROTATIONS IN RANGE [1 - K]
319 !
320  DO i = 1 , k
321 !
322 ! MODIFIES COLUMN I OF H BY THE PREVIOUS ROTATIONS
323  IF(i.GE.2) THEN
324  DO j = 1,i-1
325  zz = c(j) * h(j,i) + s(j) * h(j+1,i)
326  h(j+1,i) = -s(j) * h(j,i) + c(j) * h(j+1,i)
327  h(j,i) = zz
328  ENDDO ! J
329  ENDIF
330 ! MODIFIES COLUMN I OF H BY ROTATION I
331  r = sqrt( h(i,i)**2 + h(i+1,i)**2 )
332  IF(abs(r).LT.1.d-6) THEN
333  IF(infogr) THEN
334  WRITE(lu,92) r
335  ENDIF
336  GO TO 3000
337  ENDIF
338  c(i) = h(i,i) / r
339  s(i) = h(i+1,i) / r
340  h(i,i) = r
341 ! H(I+1,I) = 0.D0 (WILL NOT BE USED AGAIN)
342 ! MODIFIES VECTOR E
343  e(i+1) = -s(i) * e(i)
344  e(i ) = c(i) * e(i)
345 !
346  ENDDO ! I
347 !
348 !-----------------------------------------------------------------------
349 ! SOLVES SYSTEM H*Y = E (H UPPER TRIANGULAR OF DIMENSION K)
350 ! Y SAME AS E
351 !
352 ! H(I,I) <> 0 HAS ALREADY BEEN CHECKED ON R
353 !-----------------------------------------------------------------------
354 !
355  e(k) = e(k) / h(k,k)
356  DO j = k-1,1,-1
357  DO l = j+1,k
358  e(j) = e(j) - h(j,l) * e(l)
359  ENDDO
360  e(j) = e(j) / h(j,j)
361  ENDDO
362 !
363 !-----------------------------------------------------------------------
364 ! BUILDS THE SOLUTION FOR STEP M : X(M+1) = X(M) + VK * Y(K)
365 !-----------------------------------------------------------------------
366 !
367  DO l = 1,k
368 !
369 ! COMPUTES THE NEW SOLUTION X
370 !
371  CALL os('X=X+CY ', x=x , y=v%ADR(l)%P , c=e(l))
372 !
373 ! COMPUTES THE RESIDUAL : RM+1 = RM + A * VK * ZK
374 !
375  CALL os('X=X+CY ', x=r0, y=av%ADR(l)%P, c= e(l))
376 !
377  ENDDO
378 !
379 ! CHECKS THAT THE ACCURACY IS NOT ALREADY REACHED
380 ! THE RESIDUAL NORM IS GIVEN BY ABS(E(K+1))
381 !
382 ! NR0 =NORM OF R0
383  nr0 = abs(e(k+1))
384  prec = nr0/nb
385 !
386  IF (prec.GE.cfg%EPS.AND.m.LT.cfg%NITMAX) GOTO 20
387 !
388 3000 CONTINUE
389 !
390 !-----------------------------------------------------------------------
391 !
392  IF(preco) THEN
393  CALL goup( x,aux,x,'D',mesh,.false.)
394  ENDIF
395 !
396 !-----------------------------------------------------------------------
397 !
398 ! IF(INFOGR) THEN
399  IF(m.LT.cfg%NITMAX) THEN
400  IF(infogr) THEN
401  IF(relat) THEN
402  WRITE(lu,102) m,prec
403  ELSE
404  WRITE(lu,202) m,prec
405  ENDIF
406  ENDIF
407  ELSE
408  IF(relat) THEN
409  WRITE(lu,104) m,prec
410  ELSE
411  WRITE(lu,204) m,prec
412  ENDIF
413  ENDIF
414 ! ENDIF
415 !
416  RETURN
417 !
418 !-----------------------------------------------------------------------
419 !
420  92 FORMAT(1x,'GMRES (BIEF) : ALGORITHM FAILED R=',g16.7,/,
421  & 1x,' IF THE MATRIX IS DIAGONAL, CHOOSE',/,
422  & 1x,' THE CONJUGATE GRADIENT SOLVER.')
423 102 FORMAT(1x,'GMRES (BIEF) : ',
424  & 1i8,' ITERATIONS, RELATIVE PRECISION:',g16.7)
425 202 FORMAT(1x,'GMRES (BIEF) : ',
426  & 1i8,' ITERATIONS, ABSOLUTE PRECISION:',g16.7)
427 104 FORMAT(1x,'GMRES (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
428  & 1i8,' RELATIVE PRECISION:',g16.7)
429 204 FORMAT(1x,'GMRES (BIEF) : EXCEEDING MAXIMUM ITERATIONS:',
430  & 1i8,' ABSOLUTE PRECISION:',g16.7)
431 !
432 !-----------------------------------------------------------------------
433 !
434  END
435
subroutine godown(X, A, B, DITR, MESH, COPY)
Definition: godown.f:7
subroutine goup(X, A, B, DITR, MESH, COPY)
Definition: goup.f:7
subroutine gmres(X, A, B, MESH, R0, V, AV, CFG, INFOGR, AUX)
Definition: gmres.f:7
subroutine puog(X, A, B, DITR, MESH, COPY)
Definition: puog.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