V3FIT
dPackgmres.f
1 !*
2 !* Copyright (C) CERFACS 1998
3 !*
4 !* SOFTWARE LICENSE AGREEMENT NOTICE - THIS SOFTWARE IS BEING PROVIDED TO
5 !* YOU BY CERFACS UNDER THE FOLLOWING LICENSE. BY DOWN-LOADING, INSTALLING
6 !* AND/OR USING THE SOFTWARE YOU AGREE THAT YOU HAVE READ, UNDERSTOOD AND
7 !* WILL COMPLY WITH THESE FOLLOWING TERMS AND CONDITIONS.
8 !*
9 !* 1 - This software program provided in source code format ("the " Source
10 !* Code ") and any associated documentation (the " Documentation ") are
11 !* licensed, not sold, to you.
12 !*
13 !* 2 - CERFACS grants you a personal, non-exclusive, non-transferable and
14 !* royalty-free right to use, copy or modify the Source Code and
15 !* Documentation, provided that you agree to comply with the terms and
16 !* restrictions of this agreement. You may modify the Source Code and
17 !* Documentation to make source code derivative works, object code
18 !* derivative works and/or documentation derivative Works (called "
19 !* Derivative Works "). The Source Code, Documentation and Derivative
20 !* Works (called " Licensed Software ") may be used by you for personal
21 !* and non-commercial use only. " non-commercial use " means uses that are
22 !* not or will not result in the sale, lease or rental of the Licensed
23 !* Software and/or the use of the Licensed Software in any commercial
24 !* product or service. CERFACS reserves all rights not expressly granted
25 !* to you. No other licenses are granted or implied.
26 !*
27 !* 3 - The Source Code and Documentation are and will remain the sole
28 !* property of CERFACS. The Source Code and Documentation are copyrighted
29 !* works. You agree to treat any modification or derivative work of the
30 !* Licensed Software as if it were part of the Licensed Software itself.
31 !* In return for this license, you grant CERFACS a non-exclusive perpetual
32 !* paid-up royalty-free license to make, sell, have made, copy, distribute
33 !* and make derivative works of any modification or derivative work you
34 !* make of the Licensed Software.
35 !*
36 !* 4- The licensee shall acknowledge the contribution of the Source Code
37 !* (using the reference [1]) in any publication of material dependent upon
38 !* upon the use of the Source Code. The licensee shall use reasonable
39 !* endeavours to notify the authors of the package of this publication.
40 !*
41 !* [1] V. Frayssé, L. Giraud, S. Gratton, and J. Langou, A set of GMRES
42 !* routines for real and complex arithmetics on high performance
43 !* computers, CERFACS Technical Report TR/PA/03/3, public domain software
44 !* available on www.cerfacs/algor/Softs, 2003
45 !*
46 !* 5- CERFACS has no obligation to support the Licensed Software it is
47 !* providing under this license.
48 !*
49 !* THE LICENSED SOFTWARE IS PROVIDED " AS IS " AND CERFACS MAKE NO
50 !* REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED. BY WAY OF EXAMPLE,
51 !* BUT NOT LIMITATION, CERFACS MAKE NO REPRESENTATIONS OR WARRANTIES OF
52 !* MERCHANTIBILY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE OF
53 !* THE LICENSED SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD
54 !* PARTY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS. CERFACS WILL NOT
55 !* BE LIABLE FOR ANY CONSEQUENTIAL, INCIDENTAL, OR SPECIAL DAMAGES, OR ANY
56 !* OTHER RELIEF, OR FOR ANY CLAIM BY ANY THIRD PARTY, ARISING FROM YOUR
57 !* USE OF THE LICENSED SOFTWARE.
58 !*
59 !* 6- For information regarding a commercial license for the Source Code
60 !* and Documentation, please contact Mrs Campassens (campasse@cerfacs.fr)
61 !*
62 !* 7- This license is effective until terminated. You may terminate this
63 !* license at any time by destroying the Licensed Software.
64 !*
65 !* I agree all the terms and conditions of the above license agreement
66 !*
67  MODULE gmres_dim
68  INTEGER, PARAMETER :: dp = selected_real_kind(12,100)
69  INTEGER :: icheck=0
70 ! DATA icheck /0/
71  INTEGER :: iterOut, jH,
72  & retlbl=0, j, northo, comprsd
73  REAL(dp) :: dnormres, beta, bn, sA, sb, sPA, sPb,
74  & dloo, dnormx, truenormres, bea, be
75  END MODULE gmres_dim
76 
77  SUBROUTINE drive_dgmres(n,nloc,m,lwork,work,
78  & irc,icntl,cntl,info,rinfo)
79  USE gmres_dim, ONLY: dp, icheck
80  IMPLICIT NONE
81 
82 !* Purpose
83 !* =======
84 !* drive_dgmres is the driver routine for solving the linear system
85 !* Ax = b using the !* Generalized Minimal Residual iterative method
86 !* with preconditioning.
87 !* This solver is implemented with a reverse communication scheme: control
88 !* is returned to the user for computing the (preconditioned)
89 !* matrix-vector product.
90 !* See the User's Guide for an example of use.
91 !*
92 !*
93 !*Written : June 1996
94 !*Authors : Luc Giraud, Serge Gratton, V. Fraysse
95 !* Parallel Algorithms - CERFACS
96 !*
97 !*Updated : April 1997
98 !*Authors : Valerie Fraysse, Luc Giraud, Serge Gratton
99 !* Parallel Algorithms - CERFACS
100 !*
101 !*Updated : June 1998
102 !*Authors : Valerie Fraysse, Luc Giraud, Serge Gratton
103 !* Parallel Algorithms - CERFACS
104 !*Purpose : Make clear that the warning and error messages come from the
105 !* dgmres modules.
106 !*
107 !*Updated : December 2002 - L. Giraud, J.Langou
108 !*Purpose : Add the capability to avoid explicit residual calculation at restart
109 !*
110 !*Updated : March 2005 - L. Giraud
111 !*Purpose : small bug in the computation of the restart if too large vs the
112 !* size of the workspace
113 !*
114 !* Arguments
115 !* =========
116 !*
117 !* n (input) INTEGER.
118 !* On entry, the dimension of the problem.
119 !* Unchanged on exit.
120 !*
121 !* nloc (input) INTEGER.
122 !* On entry, the dimension of the local problem.
123 !* In a parallel distributed envirionment, this corresponds
124 !* to the size of the subset of entries of the right hand side
125 !* and solution allocated to the calling process.
126 !* Unchanged on exit.
127 !*
128 !*
129 !* m (input/output) INTEGER
130 !* Restart parameter, <= N. This parameter controls the amount
131 !* of memory required for matrix H (see WORK and H).
132 !* Unchanged on exit.
133 !*
134 !* lwork (input) INTEGER
135 !* size of the workspace
136 !* if (icntl(5) = 0 or 1 )
137 !* lwork >= m*m + m*(n+5) + 5*n+2, if icntl(8) = 1
138 !* lwork >= m*m + m*(n+5) + 6*n+2, if icntl(8) = 0
139 !* if (icntl(5) = 2 or 3 )
140 !* lwork >= m*m + m*(n+5) + 5*n+m+1, if icntl(8) = 1
141 !* lwork >= m*m + m*(n+5) + 6*n+m+1, if icntl(8) = 0
142 !*
143 !* work (workspace) REAL(dp)/REAL(dp) array, length lwork
144 !* work contains the required vector and matrices stored in the
145 !* following order :
146 !* x (n,1) : computed solution.
147 !* b (n,1) : right hand side.
148 !* r0 (n,1) : vector workspace.
149 !* w (n,1) : vector workspace.
150 !* V (n,m) : Krylov basis.
151 !* H (m+1,m+1) : Hessenberg matrix (full storage).
152 !* yCurrent (m,1) : solution of the current LS
153 !* xCurrent (n,1) : current iterate
154 !* rotSin (m,1) : Sine of the Givens rotation
155 !* rotCos (m,1) : Cosine of the Givens rotation
156 !*
157 !* irc (input/output) INTEGER array. length 5
158 !* irc(1) : REVCOM used for reverse communication
159 !* (type of external operation)
160 !* irc(2) : COLX used for reverse communication
161 !* irc(3) : COLY used for reverse communication
162 !* irc(4) : COLZ used for reverse communication
163 !* irc(5) : NBSCAL used for reverse communication
164 !*
165 !* icntl (input/output) INTEGER array. length 9
166 !* icntl(1) : stdout for error messages
167 !* icntl(2) : stdout for warnings
168 !* icntl(3) : stdout for convergence history
169 !* icntl(4) : 0 - no preconditioning
170 !* 1 - left preconditioning
171 !* 2 - right preconditioning
172 !* 3 - double side preconditioning
173 !* 4 - error, default set in Init
174 !* icntl(5) : 0 - modified Gram-Schmidt
175 !* 1 - iterative modified Gram-Schmidt
176 !* 2 - classical Gram-Schmidt
177 !* 3 - iterative classical Gram-Schmidt
178 !* icntl(6) : 0 - default initial guess x_0 = 0 (to be set)
179 !* 1 - user supplied initial guess
180 !* icntl(7) : maximum number of iterations
181 !* icntl(8) : 0 - use recurrence formula at restart
182 !* 1 - default compute the true residual at each restart
183 !* icntl(9) : frequency of updated solution return to user
184 !*
185 !* cntl (input) REAL(dp) array, length 5
186 !* cntl(1) : tolerance for convergence
187 !* cntl(2) : scaling factor for normwise perturbation on A
188 !* cntl(3) : scaling factor for normwise perturbation on b
189 !* cntl(4) : scaling factor for normwise perturbation on the
190 !* preconditioned matrix
191 !* cntl(5) : scaling factor for normwise perturbation on
192 !* preconditioned right hand side
193 !*
194 !* info (output) INTEGER array, length 3
195 !* info(1) : 0 - normal exit
196 !* -1 - n < 1
197 !* -2 - m < 1
198 !* -3 - lwork too small
199 !* -4 - convergence not achieved after icntl(7) iterations
200 !* -5 - precondition type not set by user
201 !* info(2) : if info(1)=0 - number of iteration to converge
202 !* if info(1)=-3 - minimum workspace size necessary
203 !* info(3) : optimal size for the workspace
204 !*
205 !* rinfo (output) REAL(dp) array, length 2
206 !* if info(1)=0
207 !* rinfo(1) : backward error for the preconditioned system
208 !* rinfo(2) : backward error for the unpreconditioned system
209 !*
210 !*Input variables
211 !*---------------
212  INTEGER, INTENT(IN) :: n, nloc, lwork
213  REAL(dp),INTENT(IN) :: cntl(5)
214  REAL(dp) sA, sb, sPA, sPb
215 !*Output variables
216 !*----------------
217  integer, INTENT(OUT) :: info(3)
218  REAL(dp),INTENT(OUT) :: rinfo(2)
219 !*Input/Output variables
220 !*----------------------
221  integer, INTENT(INOUT) :: icntl(9), irc(5), m
222  REAL(dp),INTENT(INOUT) :: work(lwork)
223 !*Local variables
224 !*---------------
225  integer xptr, bptr, wptr, r0ptr, Vptr, Hptr, dotptr
226  integer yCurrent,rotSin, rotCos, xCurrent
227  integer sizeWrk, newRestart
228  integer iwarn, ierr, ihist, compRsd
229 ! real rn,rx,rc
230  REAL(dp) :: rn,rx,rc
231  REAL(dp), PARAMETER :: DZRO=0
232 ! parameter (DZRO = 0.0d0)
233 !*
234 !*
235 !SPH - use INT, REAL intrinsic ifix, float
236 !*
237 !* Executable statements :
238 !*
239  ierr = icntl(1)
240  iwarn = icntl(2)
241  ihist = icntl(3)
242  comprsd = icntl(8)
243 
244  if (ierr.lt.0) ierr = 6
245 
246  if (comprsd.eq.1) then
247  sizewrk = m*m + m*(nloc+5) + 5*nloc+ 1
248  else
249  sizewrk = m*m + m*(nloc+5) + 6*nloc+ 1
250  endif
251 
252  if (icheck.eq.0) then
253 !*Check the value of the arguments
254  if ((n.lt.1).or.(nloc.lt.1)) then
255  write(ierr,*)
256  write(ierr,*)' ERROR GMRES : '
257  write(ierr,*)' N < 1 '
258  write(ierr,*)
259  info(1) = -1
260  irc(1) = 0
261  return
262  endif
263  if (m.lt.1) then
264  write(ierr,*)
265  write(ierr,*)' ERROR GMRES :'
266  write(ierr,*)' M < 1 '
267  write(ierr,*)
268  info(1) = -2
269  irc(1) = 0
270  return
271  endif
272  if ((icntl(4).ne.0).and.(icntl(4).ne.1).and.
273  & (icntl(4).ne.2).and.(icntl(4).ne.3)) then
274  write(ierr,*)
275  write(ierr,*)' ERROR GMRES : '
276  write(ierr,*)' Undefined preconditioner '
277  write(ierr,*)
278  info(1) = -5
279  irc(1) = 0
280  return
281  endif
282 
283  if ((icntl(5).lt.0).or.(icntl(5).gt.3)) then
284  icntl(5) = 0
285  if (iwarn.ne.0) then
286  write(iwarn,*)
287  write(iwarn,*) ' WARNING GMRES : '
288  write(iwarn,*) ' Undefined orthogonalisation '
289  write(iwarn,*) ' Default MGS '
290  write(iwarn,*)
291  endif
292  endif
293 
294  if ((icntl(5).eq.2).or.(icntl(5).eq.3)) then
295 !*the workspace should be large enough to store the m dot-products
296  sizewrk = sizewrk + m
297  else
298  sizewrk = sizewrk + 1
299  endif
300 
301  if (iwarn.ne.0) then
302  write(iwarn,*)
303  write(iwarn,*) ' WARNING GMRES : '
304  write(iwarn,*) ' For M = ',m,' optimal value '
305  write(iwarn,*) ' for LWORK = ', sizewrk
306  write(iwarn,*)
307  endif
308 
309  if ((icntl(6).ne.0).and.(icntl(6).ne.1)) then
310  icntl(6) = 0
311  if (iwarn.ne.0) then
312  write(iwarn,*)
313  write(iwarn,*) ' WARNING GMRES : '
314  write(iwarn,*) ' Undefined intial guess '
315  write(iwarn,*) ' Default x0 = 0 '
316  write(iwarn,*)
317  endif
318  endif
319  if (icntl(7).le.0) then
320  icntl(7) = n
321  if (iwarn.ne.0) then
322  write(iwarn,*)
323  write(iwarn,*) ' WARNING GMRES :'
324  write(iwarn,*) ' Negative max number of iterations'
325  write(iwarn,*) ' Default N '
326  write(iwarn,*)
327  endif
328  endif
329  if ((icntl(8).ne.0).and.(icntl(8).ne.1)) then
330  icntl(8) = 1
331  write(iwarn,*)
332  write(iwarn,*) ' WARNING GMRES :'
333  write(iwarn,*) ' Undefined strategy for the residual'
334  write(iwarn,*) ' at restart'
335  write(iwarn,*) ' Default 1 '
336  write(iwarn,*)
337  endif
338 !*Check if the restart parameter is correct and if the size of the
339 !* workspace is big enough for the restart.
340 !*If not try to fix correctly the parameters
341 
342  if ((m .gt. n).or.(lwork.lt.sizewrk)) then
343  if (m .gt. n) then
344  m = n
345  if (iwarn.ne.0) then
346  write(iwarn,*)
347  write(iwarn,*) ' WARNING GMRES : '
348  write(iwarn,*) ' Parameter M bigger than N'
349  write(iwarn,*) ' New value for M ',m
350  write(iwarn,*)
351  endif
352  if (comprsd.eq.1) then
353  sizewrk = m*m + m*(nloc+5) + 5*nloc+1
354  else
355  sizewrk = m*m + m*(nloc+5) + 6*nloc+1
356  endif
357  if ((icntl(5).eq.2).or.(icntl(5).eq.3)) then
358 !*the workspace should be large enough to store the m dot-products
359  sizewrk = sizewrk + m
360  else
361  sizewrk = sizewrk + 1
362  endif
363  endif
364  if ((lwork.lt.sizewrk).and.(n.eq.nloc)) then
365 !*Compute the maximum size of the m according to the memory space
366  rn = n !real(n)
367  rx = rn + 5 !rn + 5.0
368  rc = 5*rn+1-lwork !5.0*rn + 1 - real(lwork)
369 
370 !*Update the linear part of the second order equation to be solved
371  if ((icntl(5).eq.2).or.(icntl(5).eq.3)) then
372  rx = rx + 1
373  endif
374 !*Update the constant part of the second order equation to be solved
375 !*
376  if (icntl(8).eq.0) then
377  rc = rc + rn
378  endif
379 !*Solution of the the second order equation
380  newrestart = int((-rx+sqrt(rx*rx-4*rc))/2._dp)
381 ! newRestart = int((-rx+sqrt(rx**2-4.0*rc))/2.0)
382  if (newrestart.gt.0) then
383  m = newrestart
384  if (iwarn.ne.0) then
385  write(iwarn,*)
386  write(iwarn,*)' WARNING GMRES : '
387  write(iwarn,*)' Workspace too small for M'
388  write(iwarn,*)' New value for M ',m
389  write(iwarn,*)
390  endif
391  else
392  write(ierr,*)
393  write(ierr,*)' ERROR GMRES : '
394  write(ierr,*)' Not enough space for the problem'
395  write(ierr,*)' the space does not permit any m'
396  write(ierr,*)
397  info(1) = -3
398  irc(1) = 0
399  return
400  endif
401  endif
402  if ((lwork.lt.sizewrk).and.(n.ne.nloc)) then
403  write(ierr,*)
404  write(ierr,*)' ERROR GMRES : '
405  write(ierr,*)' Not enough space for the problem'
406  write(ierr,*)
407  info(1) = -3
408  irc(1) = 0
409  return
410  endif
411  endif
412 
413  info(3) = sizewrk
414  icheck = 1
415 
416 !*save the parameters in the history file
417 
418  if (ihist.ne.0) then
419  write(ihist,'(10x,A39)') 'CONVERGENCE HISTORY FOR GMRES'
420  write(ihist,*)
421  write(ihist,'(A30,I2)') 'Errors are displayed in unit: ',ierr
422  if (iwarn.eq.0) then
423  write(ihist,'(A27)') 'Warnings are not displayed:'
424  else
425  write(ihist,'(A32,I2)') 'Warnings are displayed in unit: ',
426  & iwarn
427  endif
428  write(ihist,'(A13,I7)') 'Matrix size: ',n
429  write(ihist,'(A19,I7)') 'Local matrix size: ',nloc
430  write(ihist,'(A9,I7)') 'Restart: ',m
431  if (icntl(4).eq.0) then
432  write(ihist,'(A18)') 'No preconditioning'
433  elseif (icntl(4).eq.1) then
434  write(ihist,'(A20)') 'Left preconditioning'
435  elseif (icntl(4).eq.2) then
436  write(ihist,'(A21)') 'Right preconditioning'
437  elseif (icntl(4).eq.3) then
438  write(ihist,'(A30)') 'Left and right preconditioning'
439  endif
440  if (icntl(5).eq.0) then
441  write(ihist,'(A21)') 'Modified Gram-Schmidt'
442  elseif (icntl(5).eq.1) then
443  write(ihist,'(A31)') 'Iterative modified Gram-Schmidt'
444  elseif (icntl(5).eq.2) then
445  write(ihist,'(A22)') 'Classical Gram-Schmidt'
446  else
447  write(ihist,'(A32)') 'Iterative classical Gram-Schmidt'
448  endif
449  if (icntl(6).eq.0) then
450  write(ihist,'(A29)') 'Default initial guess x_0 = 0'
451  else
452  write(ihist,'(A27)') 'User supplied initial guess'
453  endif
454  if (icntl(8).eq.1) then
455  write(ihist,'(A33)') 'True residual computed at restart'
456  else
457  write(ihist,'(A30)') 'Recurrence residual at restart'
458  endif
459  write(ihist,'(A30,I5)') 'Maximum number of iterations: ',
460  & icntl(7)
461  write(ihist,'(A27,E8.2)') 'Tolerance for convergence: ',
462  & cntl(1)
463 
464  write(ihist,'(A53)')
465  & 'Backward error on the unpreconditioned system Ax = b:'
466  sa = cntl(2)
467  sb = cntl(3)
468  if ((sa.eq.dzro).and.(sb.eq.dzro)) then
469  write(ihist,'(A39)')
470  & ' the residual is normalised by ||b||'
471  else
472  write(ihist,1) sa,sb
473 1 format(' the residual is normalised by ',e8.2,
474  & ' !*||x|| + ',e8.2)
475  endif
476  spa = cntl(4)
477  spb = cntl(5)
478  write(ihist,2)
479 2 format('Backward error on the preconditioned system',
480  & ' (P1)A(P2)y = (P1)b:')
481  if ((spa.eq.dzro).and.(spb.eq.dzro)) then
482  write(ihist,3)
483 3 format(' the preconditioned residual is normalised ',
484  & 'by ||(P1)b||')
485  else
486  write(ihist,4) spa, spb
487 4 format(' the preconditioned residual is normalised by ', e8.2,
488  & ' !*||(P2)y|| + ',e8.2)
489  endif
490 
491  write(ihist,5) info(3)
492 5 format('Optimal size for the local workspace:',i7)
493  write(ihist,*)
494  write(ihist,6)
495 6 format('Convergence history: b.e. on the preconditioned system')
496  write(ihist,7)
497 7 format(' Iteration Arnoldi b.e. True b.e.')
498  endif
499 
500  endif
501 !*setup some pointers on the workspace
502  xptr = 1
503  bptr = xptr + nloc
504  r0ptr = bptr + nloc
505  wptr = r0ptr + nloc
506  vptr = wptr + nloc
507  if (comprsd.eq.1) then
508  hptr = vptr + m*nloc
509  else
510  hptr = vptr + (m+1)*nloc
511  endif
512  dotptr = hptr + (m+1)*(m+1)
513  if ((icntl(5).eq.2).or.(icntl(5).eq.3)) then
514  ycurrent = dotptr + m
515  else
516  ycurrent = dotptr + 1
517  endif
518  xcurrent = ycurrent + m
519  rotsin = xcurrent + nloc
520  rotcos = rotsin + m
521 
522  call dgmres(nloc,m,work(bptr),work(xptr),
523  & work(hptr),work(wptr),work(r0ptr),work(vptr),
524  & work(dotptr),work(ycurrent),work(xcurrent),
525  & work(rotsin),work(rotcos),irc,icntl,cntl,info,rinfo)
526 
527  if (irc(1).eq.0) then
528  icheck = 0
529  endif
530 
531  return
532  end
533 
534  subroutine dgmres(n,m,b,x,H,w,r0,V,dot,yCurrent,xCurrent,rotSin,
535  & rotCos,irc,icntl,cntl,info,rinfo)
536  USE gmres_dim
537  IMPLICIT NONE
538 !*
539 !*
540 !* Purpose
541 !* =======
542 !* dgmres solves the linear system Ax = b using the
543 !* Generalized Minimal Residual iterative method
544 !*
545 !*When preconditioning is used we solve :
546 !* M_1^{-1} A M_2^{-1} y = M_1^{-1} b
547 !* x = M_2^{-1} y
548 !*
549 !* Convergence test based on the normwise backward error for
550 !* the preconditioned system
551 !*
552 !*Written : June 1996
553 !*Authors : Luc Giraud, Serge Gratton, V. Fraysse
554 !* Parallel Algorithms - CERFACS
555 !*
556 !*Updated : April 1997
557 !*Authors : Valerie Fraysse, Luc Giraud, Serge Gratton
558 !* Parallel Algorithms - CERFACS
559 !*
560 !*Updated : March 1998
561 !*Purpose : Pb with F90 on DEC ws
562 !* cure : remove "ZDSCAL" when used to initialize vectors to zero
563 !*
564 !*Updated : May 1998
565 !*Purpose : r0(1) <-- r0'r0 : pb when used with DGEMV for the dot product
566 !* cure : w(1) <-- r0'r0
567 !*
568 !*Updated : June 1998
569 !*Purpose : Make clear that the warning and error messages come from the
570 !* dgmres modules.
571 !*
572 !*Updated : February 2001 - L. Giraud
573 !*Purpose : In complex version, initializations to zero performed in complex
574 !* arithmetic to avoid implicit conversion by the compiler.
575 !*
576 !*Updated : July 2001 - L. Giraud, J. Langou
577 !*Purpose : Avoid to compute the approximate solution at each step of
578 !* the Krylov space construction when spA is zero.
579 !*
580 !*Updated : November 2002 - S. Gratton
581 !*Purpose : Use Givens rotations conform to the classical definition.
582 !* No impact one the convergence history.
583 !*
584 !*Updated : November 2002 - L. Giraud
585 !*Purpose : Properly handle the situation when the convergence is obtained
586 !* exactly at the "IterMax" iteration
587 !*
588 !*Updated : December 2002 - L. Giraud, J.Langou
589 !*Purpose : Add the capability to avoid explicit residual calculation at restart
590 !*
591 !*Updated : January 2003 - L. Giraud, S. Gratton
592 !*Purpose : Use Givens rotations from BLAS.
593 !*
594 !*Updated : March 2003 - L. Giraud
595 !*Purpose : Set back retlbl to zero, if initial guess is solution
596 !* or right-hand side is zero
597 !*
598 !*Updated : September 2003 - L. Giraud
599 !*Purpose : Include room in the workspace to store the results of the dot products
600 !* Fix the bugs that appeared when M > Nloc
601 !*
602 !* Arguments
603 !* =========
604 !*
605 !* n (input) INTEGER.
606 !* On entry, the dimension of the problem.
607 !* Unchanged on exit.
608 !*
609 !* m (input) INTEGER
610 !* Restart parameter, <= N. This parameter controls the amount
611 !* of memory required for matrix H (see WORK and H).
612 !* Unchanged on exit.
613 !*
614 !* b (input) REAL(dp)/REAL(dp)
615 !* Right hand side of the linear system.
616 !*
617 !* x (output) REAL(dp)/REAL(dp)
618 !* Computed solution of the linear system.
619 !*
620 !* H (workspace) REAL(dp)/REAL(dp)
621 !* Hessenberg matrix built within dgmres
622 !*
623 !* w (workspace) REAL(dp)/REAL(dp)
624 !* Vector used as temporary storage
625 !*
626 !* r0 (workspace) REAL(dp)/REAL(dp)
627 !* Vector used as temporary storage
628 !*
629 !* V (workspace) REAL(dp)/REAL(dp)
630 !* Basis computed by the Arnoldi's procedure.
631 !*
632 !* dot (workspace) REAL(dp)/REAL(dp)
633 !* Store the results of the dot product calculation
634 !*
635 !* yCurrent (workspace) REAL(dp)/REAL(dp)
636 !* solution of the current LS
637 !*
638 !* xCurrent (workspace) REAL(dp)/REAL(dp)
639 !* current iterate
640 !*
641 !* rotSin (workspace) REAL(dp)/REAL(dp)
642 !* Sine of the Givens rotation
643 !*
644 !* rotCos (workspace) REAL(dp)
645 !* Cosine of the Givens rotation
646 !*
647 !* irc (input/output) INTEGER array. length 5
648 !* irc(1) : REVCOM used for reverse communication
649 !* (type of external operation)
650 !* irc(2) : COLX used for reverse communication
651 !* irc(3) : COLY used for reverse communication
652 !* irc(4) : COLZ used for reverse communication
653 !* irc(5) : NBSCAL used for reverse communication
654 !*
655 !* icntl (input) INTEGER array. length 9
656 !* icntl(1) : stdout for error messages
657 !* icntl(2) : stdout for warnings
658 !* icntl(3) : stdout for convergence history
659 !* icntl(4) : 0 - no preconditioning
660 !* 1 - left preconditioning
661 !* 2 - right preconditioning
662 !* 3 - double side preconditioning
663 !* 4 - error, default set in Init
664 !* icntl(5) : 0 - modified Gram-Schmidt
665 !* 1 - iterative modified Gram-Schmidt
666 !* 2 - classical Gram-Schmidt
667 !* 3 - iterative classical Gram-Schmidt
668 !* icntl(6) : 0 - default initial guess x_0 = 0 (to be set)
669 !* 1 - user supplied initial guess
670 !* icntl(7) : maximum number of iterations
671 !* icntl(8) : 1 - default compute the true residual at each restart
672 !* 0 - use recurence formula at restart
673 !* icntl(9) : frequency of updated solution return to user
674 !*
675 !* cntl (input) REAL(dp) array, length 5
676 !* cntl(1) : tolerance for convergence
677 !* cntl(2) : scaling factor for normwise perturbation on A
678 !* cntl(3) : scaling factor for normwise perturbation on b
679 !* cntl(4) : scaling factor for normwise perturbation on the
680 !* preconditioned matrix
681 !* cntl(5) : scaling factor for normwise perturbation on
682 !* preconditioned right hand side
683 !*
684 !* info (output) INTEGER array, length 3
685 !* info(1) : 0 - normal exit
686 !* -1 - n < 1
687 !* -2 - m < 1
688 !* -3 - lwork too small
689 !* -4 - convergence not achieved after icntl(7) iterations
690 !* -5 - precondition type not set by user
691 !* info(2) : if info(1)=0 - number of iteration to converge
692 !* if info(1)=-3 - minimum workspace size necessary
693 !* info(3) : optimal size for the workspace
694 !*
695 !*rinfo (output) REAL(dp) array, length 2
696 !* if info(1)=0
697 !* rinfo(1) : backward error for the preconditioned system
698 !* rinfo(2) : backward error for the unpreconditioned system
699 !*
700 !*Input variables
701 !*---------------
702  integer, INTENT(IN) :: n, m, icntl(9)
703  REAL(dp),INTENT(IN) :: b(n), cntl(5)
704 
705 !*Output variables
706 !*----------------
707  integer, INTENT(OUT) :: info(3)
708  REAL(dp),INTENT(OUT) :: rinfo(2)
709 
710 !*Input/Output variables
711 !*----------------------
712  integer, INTENT(INOUT) :: irc(5)
713  REAL(dp),INTENT(INOUT) :: x(n), H(m+1,m+1), w(n), r0(n)
714  REAL(dp) :: V(n,m), dot(*), yCurrent(m)
715  REAL(dp) :: xCurrent(n), rotSin(m), rotCos(m)
716 
717 !*Local variables
718 !*---------------
719  INTEGER iterMax, initGuess, iOrthog
720  integer xptr, bptr, wptr, r0ptr, Vptr, Hptr, yptr, xcuptr
721  integer dotptr,itMod
722  integer typePrec, iwarn, ihist
723  REAL(dp) temp, dnormw
724  REAL(dp) dVi, aux
725  REAL(dp) auxHjj, auxHjp1j
726 
727  integer, parameter :: noPrec = 0, leftprec = 1, rightprec = 2,
728  & dbleprec = 3, peek = 5
729 
730  REAL(dp), PARAMETER :: ZERO=0, one=1
731 ! parameter (ZERO = 0.0d0, ONE = 1.0d0)
732  REAL(dp), PARAMETER :: DZRO=0, done=1
733 ! parameter (DZRO = 0.0d0, DONE = 1.0d0)
734 
735 
736 !*External functions
737 !*------------------
738  REAL(dp) :: dnrm2
739  external dnrm2
740 
741 !*Saved variables (moved to module GMRES_DIM)
742 !*---------------
743 ! save iterOut, jH, beta, bn, dnormres, retlbl, j
744 ! save sA, sb, sPA, sPb, dnormx, trueNormRes, bea, be
745 ! save dloo, nOrtho, compRsd
746 
747 
748 !*Reverse communication variables
749 !*-------------------------------
750  integer, parameter :: matvec=1, precondleft=2, precondright=3,
751  & prosca=4
752 ! parameter(matvec=1, precondLeft=2, precondRight=3, prosca=4)
753 ! integer retlbl
754 ! DATA retlbl /0/
755 
756 !* Executable statements
757 
758 !*setup some pointers on the workspace
759 
760 !!!SPH ADD
761  info = 0
762 
763  xptr = 1
764  bptr = xptr + n
765  r0ptr = bptr + n
766  wptr = r0ptr + n
767  vptr = wptr + n
768  if (icntl(8).eq.1) then
769  hptr = vptr + m*n
770  else
771  hptr = vptr + (m+1)*n
772  endif
773  dotptr = hptr + (m+1)*(m+1)
774  if ((icntl(5).eq.2).or.(icntl(5).eq.3)) then
775  yptr = dotptr + m
776  else
777  yptr = dotptr + 1
778  endif
779  xcuptr = yptr + m
780 
781  iwarn = icntl(2)
782  ihist = icntl(3)
783  typeprec = icntl(4)
784  iorthog = icntl(5)
785  initguess = icntl(6)
786  itermax = icntl(7)
787 
788  if (retlbl.eq.0) then
789  comprsd = icntl(8)
790  endif
791 !*ORNL
792  itmod = icntl(9)
793  if (itmod .eq. 0) itmod = huge(itmod)-1
794 
795  if (retlbl.ne.0) then
796  if (retlbl.eq.5) then
797  goto 5
798  else if (retlbl.eq.6) then
799  goto 6
800  else if (retlbl.eq.8) then
801  goto 8
802  else if (retlbl.eq.11) then
803  goto 11
804  else if (retlbl.eq.16) then
805  goto 16
806  else if (retlbl.eq.18) then
807  goto 18
808  else if (retlbl.eq.21) then
809  goto 21
810  else if (retlbl.eq.26) then
811  goto 26
812  else if (retlbl.eq.31) then
813  goto 31
814  else if (retlbl.eq.32) then
815  goto 32
816  else if (retlbl.eq.33) then
817  goto 33
818  else if (retlbl.eq.34) then
819  goto 34
820  else if (retlbl.eq.36) then
821  goto 36
822  else if (retlbl.eq.37) then
823  goto 37
824  else if (retlbl.eq.38) then
825  goto 38
826  else if (retlbl.eq.41) then
827  goto 41
828  else if (retlbl.eq.43) then
829  goto 43
830  else if (retlbl.eq.46) then
831  goto 46
832  else if (retlbl.eq.48) then
833  goto 48
834  else if (retlbl.eq.51) then
835  goto 51
836  else if (retlbl.eq.52) then
837  goto 52
838  else if (retlbl.eq.61) then
839  goto 61
840  else if (retlbl.eq.66) then
841  goto 66
842  else if (retlbl.eq.68) then
843  goto 68
844 !*ORNL
845  else if (retlbl.eq.69) then
846  goto 69
847  endif
848  endif
849 
850 
851 
852 !*intialization of various variables
853 
854  iterout = 0
855  beta = dzro
856 
857  if (initguess.eq.0) then
858  do j=1,n
859  x(j) = zero
860  enddo
861  endif
862 
863 !* bn = dnrm2(n,b,1)
864 
865  irc(1) = prosca
866  irc(2) = bptr
867  irc(3) = bptr
868  irc(4) = dotptr
869  irc(5) = 1
870  retlbl = 5
871  return
872  5 continue
873  bn = sqrt((dot(1)))
874 
875  if (bn.eq.dzro) then
876  do j=1,n
877  x(j) = zero
878  enddo
879  if (iwarn.ne.0) then
880  write(iwarn,*)
881  write(iwarn,*) ' WARNING GMRES : '
882  write(iwarn,*) ' Null right hand side'
883  write(iwarn,*) ' solution set to zero'
884  write(iwarn,*)
885  endif
886  jh = 0
887  bea = dzro
888  be = dzro
889  write(ihist,'(I5,11x,E8.2,$)') jh,bea
890  write(ihist,'(7x,E8.2)') be
891  info(1) = 0
892  info(2) = 0
893  rinfo(1) = dzro
894  rinfo(2) = dzro
895  irc(1) = 0
896  retlbl = 0
897  return
898  endif
899 
900 !*Compute the scaling factor for the backward error on the
901 !* unpreconditioned sytem
902 
903  sa = cntl(2)
904  sb = cntl(3)
905  if ((sa.eq.dzro).and.(sb.eq.dzro)) then
906  sb = bn
907  endif
908 !*Compute the scaling factor for the backward error on the
909 !* preconditioned sytem
910 
911  spa = cntl(4)
912  spb = cntl(5)
913  if ((spa.eq.dzro).and.(spb.eq.dzro)) then
914  if ((typeprec.eq.noprec).or.(typeprec.eq.rightprec)) then
915  spb = bn
916  else
917  irc(1) = precondleft
918  irc(2) = bptr
919  irc(4) = r0ptr
920  retlbl = 6
921  return
922  endif
923  endif
924  6 continue
925  if ((spa.eq.dzro).and.(spb.eq.dzro)) then
926  if ((typeprec.eq.dbleprec).or.(typeprec.eq.leftprec)) then
927 
928 !* sPb = dnrm2(n,r0,1)
929 
930  irc(1) = prosca
931  irc(2) = r0ptr
932  irc(3) = r0ptr
933  irc(4) = dotptr
934  irc(5) = 1
935  retlbl = 8
936  return
937  endif
938  endif
939  8 continue
940  if ((spa.eq.dzro).and.(spb.eq.dzro)) then
941  if ((typeprec.eq.dbleprec).or.(typeprec.eq.leftprec)) then
942  spb = sqrt((dot(1)))
943 !*
944  endif
945  endif
946 
947 
948 !*Compute the first residual
949 !* Y = AX : r0 <-- A x
950 
951 !*The residual is computed only if the initial guess is not zero
952 
953  if (initguess.ne.0) then
954  irc(1) = matvec
955  irc(2) = xptr
956  irc(4) = r0ptr
957  retlbl = 11
958  return
959  endif
960  11 continue
961  if (initguess.ne.0) then
962  do j=1,n
963  r0(j) = b(j)-r0(j)
964  enddo
965  else
966  call dcopy(n,b,1,r0,1)
967  endif
968 
969 !*Compute the preconditioned residual if necessary
970 !* M_1Y = X : w <-- M_1^{-1} r0
971 
972  if ((typeprec.eq.noprec).or.(typeprec.eq.rightprec)) then
973  call dcopy(n,r0,1,w,1)
974  else
975  irc(1) = precondleft
976  irc(2) = r0ptr
977  irc(4) = wptr
978  retlbl = 16
979  return
980  endif
981  16 continue
982 
983 
984 !* beta = dnrm2(n,w,1)
985 
986 
987  irc(1) = prosca
988  irc(2) = wptr
989  irc(3) = wptr
990  irc(4) = dotptr
991  irc(5) = 1
992  retlbl = 18
993  return
994  18 continue
995  beta = sqrt((dot(1)))
996 
997  if (beta .eq. dzro) then
998 !* The residual is exactly zero : x is the exact solution
999  info(1) = 0
1000  info(2) = 0
1001  rinfo(1) = dzro
1002  rinfo(2) = dzro
1003  irc(1) = 0
1004  retlbl = 0
1005  jh = 0
1006  bea = dzro
1007  be = dzro
1008  write(ihist,'(I5,11x,E8.2,$)') jh,bea
1009  write(ihist,'(7x,E8.2)') be
1010  if (iwarn.ne.0) then
1011  write(iwarn,*)
1012  write(iwarn,*) ' WARNING GMRES : '
1013  write(iwarn,*) ' Intial residual is zero'
1014  write(iwarn,*) ' initial guess is solution'
1015  write(iwarn,*)
1016  endif
1017  return
1018  endif
1019 
1020  aux = one/beta
1021  do j=1,n
1022  v(j,1) = zero
1023  enddo
1024  call daxpy(n,aux,w,1,v(1,1),1)
1025 
1026 !* Most outer loop : dgmres iteration
1027 
1028 !* REPEAT
1029  7 continue
1030 
1031 
1032  h(1,m+1)=beta
1033  do j=1,m
1034  h(j+1,m+1) = zero
1035  enddo
1036 
1037 !* Construction of the hessenberg matrix WORK and of the orthogonal
1038 !* basis V such that AV=VH
1039 
1040  jh = 1
1041  10 continue
1042 !*Remark : this do loop has been written with a while do
1043 !* because the
1044 !* " do jH=1,restart "
1045 !* fails with the reverse communication.
1046 !* do jH=1,restart
1047 
1048 
1049 !*Compute the preconditioned residual if necessary
1050 
1051  if ((typeprec.eq.rightprec).or.(typeprec.eq.dbleprec)) then
1052 
1053 !* Y = M_2^{-1}X : w <-- M_2^{-1} V(1,jH)
1054 
1055  irc(1) = precondright
1056  irc(2) = vptr + (jh-1)*n
1057  irc(4) = wptr
1058  retlbl = 21
1059  return
1060  else
1061  call dcopy(n,v(1,jh),1,w,1)
1062  endif
1063  21 continue
1064 
1065 !* Y = AX : r0 <-- A w
1066 
1067  irc(1) = matvec
1068  irc(2) = wptr
1069  irc(4) = r0ptr
1070  retlbl = 26
1071  return
1072  26 continue
1073 
1074 !* MY = X : w <-- M_1^{-1} r0
1075 
1076  if ((typeprec.eq.noprec).or.(typeprec.eq.rightprec)) then
1077  call dcopy(n,r0,1,w,1)
1078  else
1079  irc(1) = precondleft
1080  irc(2) = r0ptr
1081  irc(4) = wptr
1082  retlbl = 31
1083  return
1084  endif
1085  31 continue
1086 
1087 !*Orthogonalization using either MGS or IMGS
1088 !*
1089 !*initialize the Hessenberg matrix to zero in order to be able to use
1090 !* IMGS as orthogonalization procedure.
1091  do j=1,jh
1092  h(j,jh) = zero
1093  enddo
1094  northo = 0
1095  19 continue
1096  northo = northo +1
1097  dloo = dzro
1098 
1099  if ((iorthog.eq.0).or.(iorthog.eq.1)) then
1100 !*MGS
1101 
1102 !* do j=1,jH
1103 
1104  j = 1
1105 !* REPEAT
1106  endif
1107  23 continue
1108  if ((iorthog.eq.0).or.(iorthog.eq.1)) then
1109 
1110 !* dVi = ddot(n,V(1,j),1,w,1)
1111 
1112  irc(1) = prosca
1113  irc(2) = vptr + (j-1)*n
1114  irc(3) = wptr
1115  irc(4) = dotptr
1116  irc(5) = 1
1117  retlbl = 32
1118  return
1119  endif
1120  32 continue
1121  if ((iorthog.eq.0).or.(iorthog.eq.1)) then
1122  dvi = dot(1)
1123  h(j,jh) = h(j,jh) + dvi
1124  dloo = dloo + dvi**2
1125  aux = -dvi
1126  call daxpy(n,aux,v(1,j),1,w,1)
1127  j = j + 1
1128  if (j.le.jh) goto 23
1129 !* enddo_j
1130  else
1131 !*CGS
1132 !*Gathered dot product calculation
1133 
1134 !* call dgemv('C',n,jH,ONE,V(1,1),n,w,1,ZERO,r0,1)
1135 
1136  irc(1) = prosca
1137  irc(2) = vptr
1138  irc(3) = wptr
1139  irc(4) = dotptr
1140  irc(5) = jh
1141  retlbl = 34
1142  return
1143  endif
1144  34 continue
1145  if ((iorthog.eq.2).or.(iorthog.eq.3)) then
1146 
1147  call daxpy(jh,one,dot,1,h(1,jh),1)
1148  call dgemv('N',n,jh,-one,v(1,1),n,dot,1,one,w,1)
1149  dloo = dnrm2(jh,dot,1)**2
1150  endif
1151 
1152 !* dnormw = dnrm2(n,w,1)
1153 
1154  irc(1) = prosca
1155  irc(2) = wptr
1156  irc(3) = wptr
1157  irc(4) = dotptr
1158  irc(5) = 1
1159  retlbl = 33
1160  return
1161  33 continue
1162  dnormw = sqrt((dot(1)))
1163 
1164  if ((iorthog.eq.1).or.(iorthog.eq.3)) then
1165 !*IMGS / CGS orthogonalisation
1166  dloo = sqrt(dloo)
1167 !*check the orthogonalization quality
1168  if (((2*dnormw).le.dloo).and.(northo.lt.3)) then
1169  goto 19
1170  endif
1171  endif
1172 
1173  h(jh+1,jh) = dnormw
1174  if ((jh.lt.m).or.(icntl(8).eq.0)) then
1175  aux = one/dnormw
1176  do j=1,n
1177  v(j,jh+1) = zero
1178  enddo
1179  call daxpy(n,aux,w,1,v(1,jh+1),1)
1180  endif
1181 !*Apply previous Givens rotations to the new column of H
1182  do j=1,jh-1
1183  call drot(1, h(j,jh), 1, h(j+1,jh), 1,(rotcos(j)),
1184  & rotsin(j))
1185  enddo
1186  auxhjj = h(jh,jh)
1187  auxhjp1j= h(jh+1,jh)
1188  call drotg(auxhjj, auxhjp1j,temp,rotsin(jh))
1189  rotcos(jh) = temp
1190 !*Apply current rotation to the rhs of the least squares problem
1191  call drot(1, h(jh,m+1), 1, h(jh+1,m+1), 1, (rotcos(jh)),
1192  & rotsin(jh))
1193 
1194 !*zabs(H(jH+1,m+1)) is the residual computed using the least squares
1195 !* solver
1196 !*Complete the QR factorisation of the Hessenberg matrix by apply the current
1197 !*rotation to the last entry of the collumn
1198  call drot(1, h(jh,jh), 1, h(jh+1,jh), 1, (rotcos(jh)),
1199  & rotsin(jh))
1200  h(jh+1,jh) = zero
1201 
1202 !*Get the Least square residual
1203 
1204  dnormres = abs(h(jh+1,m+1))
1205 !*ORNL
1206  if ((spa.ne.dzro).or.(mod(iterout*m+jh,itmod).eq.0)) then
1207 
1208 !*Compute the solution of the current linear least squares problem
1209 
1210  call dcopy(jh,h(1,m+1),1,ycurrent,1)
1211  call dtrsv('U','N','N',jh,h,m+1,ycurrent,1)
1212 
1213 !*Compute the value of the new iterate
1214 
1215  call dgemv('N',n,jh,one,v,n,
1216  & ycurrent,1,zero,xcurrent,1)
1217 
1218  if ((typeprec.eq.rightprec).or.(typeprec.eq.dbleprec)) then
1219 
1220 !* Y = M_2^{-1}X : r0 <-- M_2^{-1} xCurrent
1221 
1222  irc(1) = precondright
1223  irc(2) = xcuptr
1224  irc(4) = r0ptr
1225  retlbl = 36
1226  return
1227  else
1228  call dcopy(n,xcurrent,1,r0,1)
1229  endif
1230  endif
1231  36 continue
1232 
1233 
1234  if (spa.ne.dzro) then
1235 !*Update the current solution
1236  call dcopy(n,x,1,xcurrent,1)
1237  call daxpy(n,one,r0,1,xcurrent,1)
1238 
1239 !* dnormx = dnrm2(n,xCurrent,1)
1240 
1241  irc(1) = prosca
1242  irc(2) = xcuptr
1243  irc(3) = xcuptr
1244  irc(4) = dotptr
1245  irc(5) = 1
1246  retlbl = 38
1247  return
1248  else
1249  dnormx = done
1250  endif
1251  38 continue
1252  if (spa.ne.dzro) then
1253  dnormx = sqrt((dot(1)))
1254  endif
1255 
1256  bea = dnormres/(spa*dnormx+spb)
1257 !*ORNL
1258  if (mod(iterout*m+jh,itmod).eq.0) then
1259  call dcopy(n,x,1,xcurrent,1)
1260  call daxpy(n,one,r0,1,xcurrent,1)
1261  irc(1) = peek
1262  irc(2) = xcuptr
1263  info(1) = iterout*m+jh
1264  rinfo(1) = bea
1265  retlbl = 69
1266  return
1267  endif
1268  69 continue
1269 
1270 !*Check the convergence based on the Arnoldi Backward error for the
1271 !*preconditioned system
1272  if ((bea.le.cntl(1)).or.(iterout*m+jh.ge.itermax)) then
1273 !*
1274 !*The Arnoldi Backward error indicates that dgmres might have converge
1275 !*enforce the calculation of the true residual at next restart
1276  comprsd = 1
1277 
1278 !* If the update of X has not yet been performed
1279  if (spa.eq.dzro) then
1280 
1281 !*Compute the solution of the current linear least squares problem
1282 
1283  call dcopy(jh,h(1,m+1),1,ycurrent,1)
1284  call dtrsv('U','N','N',jh,h,m+1,ycurrent,1)
1285 
1286 !*Compute the value of the new iterate
1287 
1288  call dgemv('N',n,jh,one,v,n,
1289  & ycurrent,1,zero,xcurrent,1)
1290 
1291  if ((typeprec.eq.rightprec).or.(typeprec.eq.dbleprec)) then
1292 !*
1293 !* Y = M_2^{-1}X : r0 <-- M_2^{-1} xCurrent
1294 
1295  irc(1) = precondright
1296  irc(2) = xcuptr
1297  irc(4) = r0ptr
1298  retlbl = 37
1299  return
1300  else
1301  call dcopy(n,xcurrent,1,r0,1)
1302  endif
1303  endif
1304  endif
1305  37 continue
1306  if ((bea.le.cntl(1)).or.(iterout*m+jh.ge.itermax)) then
1307  if (spa.eq.dzro) then
1308 !*Update the current solution
1309  call dcopy(n,x,1,xcurrent,1)
1310  call daxpy(n,one,r0,1,xcurrent,1)
1311  endif
1312 
1313  call dcopy(n,xcurrent,1,r0,1)
1314 !*Compute the true residual, the Arnoldi one may be unaccurate
1315 
1316 !* Y = AX : w <-- A r0
1317 
1318  irc(1) = matvec
1319  irc(2) = r0ptr
1320  irc(4) = wptr
1321  retlbl = 41
1322  return
1323  endif
1324  41 continue
1325  if ((bea.le.cntl(1)).or.(iterout*m+jh.ge.itermax)) then
1326 
1327  do j=1,n
1328  w(j) = b(j) - w(j)
1329  enddo
1330 !*Compute the norm of the unpreconditioned residual
1331 
1332 !* trueNormRes = dnrm2(n,w,1)
1333 
1334  irc(1) = prosca
1335  irc(2) = wptr
1336  irc(3) = wptr
1337  irc(4) = dotptr
1338  irc(5) = 1
1339  retlbl = 43
1340  return
1341  endif
1342  43 continue
1343  if ((bea.le.cntl(1)).or.(iterout*m+jh.ge.itermax)) then
1344  truenormres = sqrt((dot(1)))
1345 
1346  if ((typeprec.eq.leftprec).or.(typeprec.eq.dbleprec)) then
1347 
1348 !* MY = X : r0 <-- M_1^{-1} w
1349 
1350  irc(1) = precondleft
1351  irc(2) = wptr
1352  irc(4) = r0ptr
1353  retlbl = 46
1354  return
1355  else
1356  call dcopy(n,w,1,r0,1)
1357  endif
1358  endif
1359  46 continue
1360  if ((bea.le.cntl(1)).or.(iterout*m+jh.ge.itermax)) then
1361 
1362 !* dnormres = dnrm2(n,r0,1)
1363 
1364  irc(1) = prosca
1365  irc(2) = r0ptr
1366  irc(3) = r0ptr
1367  irc(4) = dotptr
1368  irc(5) = 1
1369  retlbl = 48
1370  return
1371  endif
1372  48 continue
1373  if ((bea.le.cntl(1)).or.(iterout*m+jh.ge.itermax)) then
1374  dnormres = sqrt((dot(1)))
1375 
1376  be = dnormres/(spa*dnormx+spb)
1377 !*Save the backward error on a file if convergence history requested
1378  if (ihist.ne.0) then
1379  write(ihist,1000)iterout*m+jh,bea,be
1380 1000 format(i5,11x,e8.2,7x,e8.2)
1381  endif
1382 
1383  endif
1384 
1385 
1386 !*Check again the convergence
1387  if ((bea.le.cntl(1)).or.(iterout*m+jh.ge.itermax)) then
1388  if ((be.le.cntl(1)).or.(iterout*m+jh.ge.itermax)) then
1389 !*The convergence has been achieved, we restore the solution in x
1390 !*and compute the two backward errors.
1391  call dcopy(n,xcurrent,1,x,1)
1392 
1393  if (sa.ne.dzro) then
1394 
1395 !* dnormx = dnrm2(n,x,1)
1396 
1397  irc(1) = prosca
1398  irc(2) = xptr
1399  irc(3) = xptr
1400  irc(4) = dotptr
1401  irc(5) = 1
1402  retlbl = 51
1403  return
1404  endif
1405  endif
1406  endif
1407  51 continue
1408  if ((bea.le.cntl(1)).or.(iterout*m+jh.ge.itermax)) then
1409  if ((be.le.cntl(1)).or.(iterout*m+jh.ge.itermax)) then
1410  if (sa.ne.dzro) then
1411  dnormx = sqrt((dot(1)))
1412 
1413  else
1414  dnormx = done
1415  endif
1416 !*Return the backward errors
1417  rinfo(1) = be
1418  rinfo(2) = truenormres/(sa*dnormx+sb)
1419  if (be.le.cntl(1)) then
1420  info(1) = 0
1421  if (ihist.ne.0) then
1422  write(ihist,*)
1423  write(ihist,'(A20)') 'Convergence achieved'
1424  endif
1425  else if (be.gt.cntl(1)) then
1426  if (iwarn.ne.0) then
1427  write(iwarn,*)
1428  write(iwarn,*) ' WARNING GMRES : '
1429  write(iwarn,*) ' No convergence after '
1430  write(iwarn,*) iterout*m+jh,' iterations '
1431  write(iwarn,*)
1432  endif
1433  if (ihist.ne.0) then
1434  write(ihist,*)
1435  write(ihist,*) ' WARNING GMRES :'
1436  write(ihist,*) ' No convergence after '
1437  write(ihist,*) iterout*m+jh,' iterations '
1438  write(ihist,*)
1439  endif
1440  info(1) = -4
1441  endif
1442  if (ihist.ne.0) then
1443  write(ihist,1010) rinfo(1)
1444  write(ihist,1011) rinfo(2)
1445 1010 format('B.E. on the preconditioned system: ',e8.2)
1446 1011 format('B.E. on the unpreconditioned system: ',e8.2)
1447  endif
1448  info(2) = iterout*m+jh
1449  if (ihist.ne.0) then
1450  write(ihist,'(A10,I2)') 'info(1) = ',info(1)
1451  write(ihist,'(A32,I5)')
1452  & 'Number of iterations (info(2)): ',info(2)
1453  endif
1454  irc(1) = 0
1455  retlbl = 0
1456  return
1457  endif
1458  else
1459 !*Save the backward error on a file if convergence history requested
1460  if (ihist.ne.0) then
1461  write(ihist,1001)iterout*m+jh,bea
1462 1001 format(i5,11x,e8.2,7x,'--')
1463  endif
1464 
1465  endif
1466 
1467  jh = jh + 1
1468  if (jh.le.m) then
1469  goto 10
1470  endif
1471 
1472  iterout = iterout + 1
1473 
1474 !*we have completed the Krylov space construction, we restart if
1475 !*we have not yet exceeded the maximum number of iterations allowed.
1476 
1477  if ((spa.eq.dzro).and.(bea.gt.cntl(1))) then
1478 
1479 !*Compute the solution of the current linear least squares problem
1480 
1481  jh = jh - 1
1482  call dcopy(jh,h(1,m+1),1,ycurrent,1)
1483  call dtrsv('U','N','N',jh,h,m+1,ycurrent,1)
1484 
1485 !*Compute the value of the new iterate
1486 
1487  call dgemv('N',n,jh,one,v,n,
1488  & ycurrent,1,zero,xcurrent,1)
1489 
1490  if ((typeprec.eq.rightprec).or.(typeprec.eq.dbleprec)) then
1491 
1492 !* Y = M_2^{-1}X : r0 <-- M_2^{-1} xCurrent
1493 
1494  irc(1) = precondright
1495  irc(2) = xcuptr
1496  irc(4) = r0ptr
1497  retlbl = 52
1498  return
1499  else
1500  call dcopy(n,xcurrent,1,r0,1)
1501  endif
1502  endif
1503  52 continue
1504  if ((spa.eq.dzro).and.(bea.gt.cntl(1))) then
1505 !*Update the current solution
1506  call dcopy(n,x,1,xcurrent,1)
1507  call daxpy(n,one,r0,1,xcurrent,1)
1508  endif
1509 
1510  call dcopy(n,xcurrent,1,x,1)
1511 
1512  if (comprsd.eq.1) then
1513 
1514 !*Compute the true residual
1515 
1516  call dcopy(n,x,1,w,1)
1517  irc(1) = matvec
1518  irc(2) = wptr
1519  irc(4) = r0ptr
1520  retlbl = 61
1521  return
1522  endif
1523  61 continue
1524  if (comprsd.eq.1) then
1525  do j=1,n
1526  r0(j) = b(j) - r0(j)
1527  enddo
1528 
1529 !*Precondition the new residual if necessary
1530 
1531  if ((typeprec.eq.leftprec).or.(typeprec.eq.dbleprec)) then
1532 
1533 !* MY = X : w <-- M_1^{-1} r0
1534 
1535  irc(1) = precondleft
1536  irc(2) = r0ptr
1537  irc(4) = wptr
1538  retlbl = 66
1539  return
1540  else
1541  call dcopy(n,r0,1,w,1)
1542  endif
1543  endif
1544  66 continue
1545 
1546 !* beta = dnrm2(n,w,1)
1547 
1548  if (comprsd.eq.1) then
1549  irc(1) = prosca
1550  irc(2) = wptr
1551  irc(3) = wptr
1552  irc(4) = dotptr
1553  irc(5) = 1
1554  retlbl = 68
1555  return
1556  endif
1557  68 continue
1558  if (comprsd.eq.1) then
1559  beta = sqrt((dot(1)))
1560 
1561  else
1562 !*Use recurrence to approximate the residual at restart
1563  beta = abs(h(m+1,m+1))
1564 !*Apply the Givens rotation is the reverse order
1565  do j=m,1,-1
1566  h(j,m+1) = zero
1567  call drot(1, h(j,m+1), 1, h(j+1,m+1), 1,
1568  & (rotcos(j)), -rotsin(j))
1569  enddo
1570 
1571 !*On applique les vecteurs V
1572 
1573  call dgemv('N',n,m+1,one,v,n,h(1,m+1),1,zero,w,1)
1574 
1575  endif
1576  do j=1,n
1577  v(j,1) = zero
1578  enddo
1579  aux = one/beta
1580  call daxpy(n,aux,w,1,v(1,1),1)
1581 
1582  goto 7
1583 
1584  end
1585 
1586 
1587  subroutine init_dgmres(icntl,cntl)
1588  USE gmres_dim
1589  IMPLICIT NONE
1590 !*
1591 !* Purpose
1592 !* =======
1593 !* Set default values for the parameters defining the characteristics
1594 !*of the Gmres algorithm.
1595 !* See the User's Guide for an example of use.
1596 !*
1597 !*
1598 !*Written : April 1997
1599 !*Authors : Valerie Fraysse, Luc Giraud, Serge Gratton
1600 !* Parallel Algorithms - CERFACS
1601 !*
1602 !*
1603 !* Arguments
1604 !* =========
1605 !*
1606 !*icntl (output) INTEGER array. length 9
1607 !* icntl(1) : stdout for error messages
1608 !* icntl(2) : stdout for warnings
1609 !* icntl(3) : stdout for convergence history
1610 !* icntl(4) : 0 - no preconditioning
1611 !* 1 - left preconditioning
1612 !* 2 - right preconditioning
1613 !* 3 - double side preconditioning
1614 !* 4 - error, default set in Init
1615 !* icntl(5) : 0 - modified Gram-Schmidt
1616 !* 1 - iterative modified Gram-Schmidt
1617 !* 2 - classical Gram-Schmidt
1618 !* 3 - iterative classical Gram-Schmidt
1619 !* icntl(6) : 0 - default initial guess x_0 = 0 (to be set)
1620 !* 1 - user supplied initial guess
1621 !* icntl(7) : maximum number of iterations
1622 !* icntl(8) : 1 - default compute the true residual at each restart
1623 !* 0 - use recurence formaula at restart
1624 !* icntl(9) : frequency of updated solution return to user
1625 !*
1626 !*cntl (output) REAL(dp) array, length 5
1627 !* cntl(1) : tolerance for convergence
1628 !* cntl(2) : scaling factor for normwise perturbation on A
1629 !* cntl(3) : scaling factor for normwise perturbation on b
1630 !* cntl(4) : scaling factor for normwise perturbation on the
1631 !* preconditioned matrix
1632 !* cntl(5) : scaling factor for normwise perturbation on
1633 !* preconditioned right hand side
1634 !*
1635 !*Output variables
1636 !*----------------
1637  integer, INTENT(OUT) :: icntl(9)
1638  REAL(dp),INTENT(OUT) :: cntl(5)
1639 
1640  icntl(1) = 6
1641  icntl(2) = 6
1642  icntl(3) = 0
1643  icntl(4) = 4
1644  icntl(5) = 0
1645  icntl(6) = 0
1646  icntl(7) = -1
1647  icntl(8) = 1
1648  icntl(9) = huge(icntl(9))-1
1649 
1650  cntl(1) = 1.e-5_dp !1.0 d -5
1651  cntl(2) = 0 !0.0 d 0
1652  cntl(3) = 0 !0.0 d 0
1653  cntl(4) = 0 !0.0 d 0
1654  cntl(5) = 0 !0.0 d 0
1655 
1656  return
1657  end
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17