V3FIT
gmres_lib.f
1  MODULE gmres_lib
2  USE stel_kinds, ONLY: dp
3  USE stel_constants, ONLY: one, zero
4  USE mpi_inc
5  IMPLICIT NONE
6 
7  TYPE gmres_info
8  INTEGER :: m, mblk_size, icntl(9), info(3)
9  INTEGER :: endglobrow, startglobrow, iam, nprocs,
10  1 ngmres_type=2
11 #if defined(mpi_opt)
12  INTEGER :: my_comm=mpi_comm_world,
13  1 my_comm_world=mpi_comm_world
14 #else
15  INTEGER :: my_comm=0,
16  1 my_comm_world=0
17 #endif
18  INTEGER, POINTER :: rcounts(:), disp(:)
19  LOGICAL :: lactive = .true.
20  REAL(dp) :: cntl(5), ftol
21  LOGICAL :: lverbose = .true.
22  LOGICAL :: l_nonlinear = .true.
23  END TYPE gmres_info
24 
25  INTEGER, PARAMETER :: finish=0, matveci=1, precondleft=2,
26  1 precondright=3, dotprod=4, peek=5
27 
28  CONTAINS
29 
30  SUBROUTINE gmres_ser (n, gi, yAx, apply_precond,
31  & getnlforce, x0, b)
32 !-----------------------------------------------
33 ! D u m m y A r g u m e n t s
34 !-----------------------------------------------
35  INTEGER :: n
36  TYPE(GMRES_INFO) :: gi
37  REAL(dp), INTENT(IN) :: b(n)
38  REAL(dp), INTENT(INOUT) :: x0(n)
39 ! REAL(dp) :: cntl(5)
40 !-----------------------------------------------
41 ! L o c a l V a r i a b l e s
42 !-----------------------------------------------
43  INTEGER :: revcom, colx, coly, colz, nbscal, m, iam
44  INTEGER :: irc(5)
45  INTEGER :: nout, lwork, icount, jcount, kout, kprint
46  REAL(dp) :: rinfo(2), fsq_nl=-1, fsq_lin, fsq_min,
47  1 delfsq, fsq_last, bnorm, xmod, xmax,
48  2 fsq_min_lin, gsum
49  REAL(dp), TARGET, ALLOCATABLE :: work(:)
50  REAL(dp), ALLOCATABLE :: xmin(:)
51  REAL(dp), POINTER :: sx(:), sy(:), sz(:)
52  LOGICAL :: lprint, lstore, l_nonlinear
53 !-----------------------------------------------
54  EXTERNAL yax, apply_precond, getnlforce
55 !-----------------------------------------------
56 !
57 ! EASY-TO-USE WRAPPER FOR GMRES DRIVER CALL
58 !
59 ! X0: on input, initial guess if icntl(6) == 1
60 ! on output, solution of Ax = b
61 ! NOTE: it is not overwritten until the END of this routine
62 !
63  bnorm = sqrt(sum(b*b))
64  IF (bnorm .EQ. 0) RETURN
65 
66 ! PRINT *,'SERIAL BNORM: ', bnorm
67 
68  fsq_min = 1; fsq_min_lin = 1
69  delfsq = 1
70  jcount = 0; icount = 0
71  kout = 0; kprint = 0
72  l_nonlinear = gi%l_nonlinear
73 
74  m = gi%m
75  iam = gi%iam
76  lprint = (iam.EQ.0 .AND. gi%lverbose)
77 
78 ! lwork = m**2 + m*(n+6) + 5*n + 1
79  lwork = m**2 + m*(n+6) + 6*n + 1 !Additional space for peek revcom (5 -> 6)
80 
81  ALLOCATE (work(lwork), stat=nout)
82  IF (nout .NE. 0) stop 'Allocation error in gmres!'
83  work = 0
84  IF (gi%icntl(6) .EQ. 1) work(1:n) = x0/bnorm
85  work(n+1:2*n) = b(1:n)/bnorm
86 
87  IF (lprint) print 900
88  900 FORMAT(/,1x,'GMRES (SERIAL) CONVERGENCE SUMMARY',/,1x,15('-'))
89 
90  !****************************************
91  !* Reverse communication implementation
92  !****************************************
93 
94  10 CONTINUE
95  CALL drive_dgmres(n, n, m, lwork, work, irc,
96  & gi%icntl, gi%cntl, gi%info, rinfo)
97  revcom = irc(1)
98  colx = irc(2)
99  coly = irc(3)
100  colz = irc(4)
101  nbscal = irc(5)
102  sx => work(colx:); sy => work(coly:); sz => work(colz:)
103 
104  IF (revcom .EQ. matveci) THEN
105 ! perform the matrix vector product work(colz) <-- A * work(colx)
106  CALL yax (sx, sz, n)
107 !Debug
108 ! gsum = SUM(work(colz:colz+n-1)**2)
109 ! IF (lprint) PRINT *,'CALL MATVECI, |Ax| = ', SQRT(gsum)
110 !End Debug
111  GOTO 10
112 
113  ELSE IF (revcom.eq.precondleft) THEN
114 ! perform the left preconditioning
115 ! IF (lprint) PRINT *,'CALL PRECONDL'
116 ! work(colz) <-- M^{-1} * work(colx)
117 ! WRITE(10000,*) "left_dcopy"; CALL FLUSH(10000)
118  CALL dcopy(n,sx,1,sz,1)
119 ! WRITE(10000,*) "precondLeft"; CALL FLUSH(10000)
120  CALL apply_precond(sz)
121  GOTO 10
122 
123  ELSE IF (revcom .EQ. precondright) THEN
124 ! perform the right preconditioning
125 ! IF (lprint) PRINT *,'CALL PRECONDR'
126 ! WRITE(10000,*) "right_dcopy"; CALL FLUSH(10000)
127  CALL dcopy(n,sx,1,sz,1)
128 ! WRITE(10000,*) "precondRight"; CALL FLUSH(10000)
129  CALL apply_precond(sz)
130  GOTO 10
131 
132  ELSE IF (revcom .EQ. dotprod) THEN
133 ! perform the scalar product
134 ! work(colz) <-- work(colx) work(coly)
135 !
136 ! CALL dgemv('C',n,nbscal,ONE,sx,n,sy,1,ZERO,sz,1)
137 ! WRITE(10000,*) "dotProd/truncate"; CALL FLUSH(10000)
138  DO nout = 0, nbscal-1
139  work(colz+nout) = sum(work(colx:colx+n-1)*work(coly:coly+n-1))
140  CALL truncate(work(colz+nout), 15)
141  colx = colx+n
142  END DO
143 
144 ! Debug
145 ! IF (lprint) THEN
146 ! gsum = SUM(work(colz:colz+nbscal-1)**2)
147 ! PRINT 200,' CALL DOTPROD: nbscal = ', nbscal,
148 ! 1 ' |WORK|: ', SQRT(gsum)
149 ! 200 FORMAT(a,i4,a,1p,e14.6)
150 ! END IF
151 ! End Debug
152  GOTO 10
153 
154  ELSE IF (revcom .EQ. peek) THEN
155  fsq_last = fsq_nl !need for delfsq criteria
156  IF (l_nonlinear) CALL getnlforce(work(colx), fsq_nl, bnorm) !get nonlinear force
157 
158  fsq_lin = (rinfo(1)*bnorm)**2
159  icount = icount+1
160  lstore = (fsq_nl.LT.fsq_min .OR. icount.EQ.1)
161  fsq_min = min(fsq_min, fsq_nl)
162  fsq_min_lin = min(fsq_min_lin, fsq_lin)
163  IF (fsq_lin .LT. 1.e-30_dp) gi%icntl(7)=gi%info(1)
164 
165  IF (lprint .AND. icount.EQ.1) THEN
166  IF (l_nonlinear) THEN
167  print 905
168  ELSE
169  print 906
170  END IF
171  END IF
172 
173  kout = gi%info(1)
174  IF (icount.EQ.1 .OR. mod(icount,10).EQ.0) THEN
175  kprint = kout
176  xmod = bnorm*sqrt(sum(work(colx:colx+n-1)**2))
177  xmax = bnorm*maxval(abs(work(colx:colx+n-1)))
178  IF (icount .EQ. 1) THEN
179  fsq_min = fsq_nl; fsq_min_lin = fsq_lin
180  END IF
181  IF (lprint) THEN
182  IF (l_nonlinear) THEN
183  print 910, kout, fsq_nl, xmod, xmax, fsq_lin
184  ELSE
185  print 911, kout, xmod, xmax, fsq_lin
186  END IF
187  END IF
188  END IF
189 
190  IF (gi%ngmres_type .EQ. 1) GOTO 10 !OLDSTYLE:IGNORE LOGIC BELOW
191 
192  delfsq = (fsq_last-fsq_nl)/fsq_min
193  IF (delfsq .LT. 0.05_dp) THEN
194  jcount = jcount+1
195  ELSE
196  jcount = 0
197  END IF
198 
199 ! STOPPING CRITERIA (reset max iterations to current iteration)
200  IF (fsq_nl.GT.(3*fsq_min) .OR. jcount.GT.3
201  & .OR. fsq_nl.LE.gi%ftol) gi%icntl(7) = gi%info(1)
202  IF (lstore) THEN
203  IF (.NOT.ALLOCATED(xmin)) ALLOCATE (xmin(n))
204  xmin = work(colx:colx+n-1)
205  END IF
206  GOTO 10
207 
208  ENDIF
209  905 FORMAT(1x,'ITER',7x,'FSQ_NL',10x,'||X||',9x,'MAX|X|',9x,'FSQ_ARN')
210  906 FORMAT(1x,'ITER',7x,'||X||',11x,'MAX|X|',9x,'FSQ_ARN')
211  910 FORMAT(i5, 4(3x,1pe12.3))
212  911 FORMAT(i5, 3(3x,1pe12.3))
213 
214  100 CONTINUE
215 
216  IF (ALLOCATED (xmin)) THEN
217  x0(1:n) = xmin(1:n)
218  DEALLOCATE (xmin)
219  ELSE
220  x0(1:n) = work(1:n)
221  fsq_min = fsq_nl
222  IF (.NOT.l_nonlinear) fsq_min = fsq_lin
223  END IF
224 
225  gi%ftol = fsq_min
226  x0 = bnorm*x0
227  IF (lprint .AND. kprint.NE.kout) THEN
228  xmod = sqrt(sum(x0*x0)); xmax = maxval(abs(x0))
229  IF (l_nonlinear) THEN
230  print 910, kout, fsq_min, xmod, xmax, fsq_min_lin
231  ELSE
232  print 911, kout, xmod, xmax, fsq_lin
233  END IF
234  END IF
235 
236  DEALLOCATE (work)
237 
238 !Sanity check
239 ! CALL GetNLForce(x0, fsq_nl, one)
240 ! IF (iam .EQ. 0) PRINT *,'FSQ_MIN: ', fsq_min,' FSQ_NL: ', fsq_nl
241 
242  END SUBROUTINE gmres_ser
243 
244  SUBROUTINE gmres_par (n, gi, yAx, apply_precond, GetNLForce,
245  & x0, b)
246 !-----------------------------------------------
247 ! D u m m y A r g u m e n t s
248 !-----------------------------------------------
249  INTEGER, INTENT(IN) :: n
250  TYPE(GMRES_INFO) :: gi
251  REAL(dp), INTENT(IN) :: b(n)
252  REAL(dp), INTENT(INOUT) :: x0(n)
253  EXTERNAL apply_precond, yax, getnlforce
254 #if defined(MPI_OPT)
255 !-----------------------------------------------
256 ! L o c a l V a r i a b l e s
257 !-----------------------------------------------
258  INTEGER :: revcom, colx, coly, colz, nbscal, m, iam
259  INTEGER :: irc(5), jcount, icount
260  INTEGER :: nout, lwork, kout, kprint, MPI_ERR,
261  & my_comm_world, my_comm
262  INTEGER, ALLOCATABLE :: itest(:)
263  REAL(dp) :: rinfo(2), fsq_nl, fsq_min, fsq_last, delfsq, fsq_lin,
264  & fsq_min_lin
265  REAL(dp), ALLOCATABLE :: work(:), xmin(:)
266  INTEGER :: nloc, myrowstart, myrowend
267  REAL(dp), ALLOCATABLE, DIMENSION(:) :: aux, tmpbuf
268  REAL(dp) :: skston, skstoff, xmax, xmod, bnorm, lsum, gsum
269  LOGICAL, PARAMETER :: lfast=.false.
270  LOGICAL :: lactive = .true., lprint, lstore, l_nonlinear
271 !-----------------------------------------------
272 !
273 ! EASY-TO-USE WRAPPER FOR GMRES PARALLEL DRIVER CALL
274 !
275 ! X0: on input, initial guess if icntl(6) == 1
276 ! on output, solution of Ax = b
277 ! NOTE: it is not overwritten UNTIL the end of this routine
278 
279  lactive = gi%lactive
280  my_comm_world = gi%MY_COMM_WORLD
281  my_comm = gi%MY_COMM
282  l_nonlinear = gi%l_nonlinear
283 
284  fsq_min = 1; fsq_min_lin = 1
285  delfsq = 1
286  jcount = 0; icount = 0; kout = 0; kprint = 0
287 
288  nloc=(gi%endglobrow-gi%startglobrow+1)*gi%mblk_size
289  myrowstart=(gi%startglobrow-1)*gi%mblk_size+1
290  myrowend=myrowstart+nloc-1
291 
292  IF (.NOT.lactive) THEN
293  nloc = 1
294  myrowstart = 1; myrowend = myrowstart
295  END IF
296 
297  ALLOCATE(tmpbuf(n), aux(n), itest(gi%nprocs), stat=nout)
298  IF (nout .NE. 0) stop 'Allocation error in gmres_fun!'
299 
300  m = gi%m; iam = gi%iam;
301  lprint = (iam.EQ.0 .AND. gi%lverbose .AND. lactive)
302  lwork = m**2 + m*(nloc+6) + 6*nloc + 1 !Additional space for peek revcom (5 -> 6)
303  ALLOCATE (work(lwork), stat=nout)
304  IF (nout .NE. 0) stop 'Allocation error in gmres_fun!'
305  work = 0
306 
307 ! PRINT *,'In gmres_par, iam: ',iam,' lactive: ',lactive
308 
309  lactive0: IF (lactive) THEN
310  lsum=sum(b(myrowstart:myrowend)**2)
311  CALL mpi_allreduce(lsum,gsum,1,mpi_real8,mpi_sum,my_comm,mpi_err)
312  bnorm=sqrt(gsum)
313  IF (bnorm .EQ. 0) RETURN
314 
315  IF (gi%icntl(6) .EQ. 1)
316  & work(1:nloc) = x0(myrowstart:myrowend)/bnorm
317 
318  work(nloc+1:2*nloc) = b(myrowstart:myrowend)/bnorm
319 
320  IF (lprint) print 900
321 900 FORMAT(/,1x,'GMRES (PARALLEL) CONVERGENCE SUMMARY',/,1x,15('-'))
322 
323  ENDIF lactive0
324 
325  CALL mpi_bcast(bnorm, 1, mpi_real8, 0, my_comm_world,
326  & mpi_err)
327 
328 !****************************************
329 !* Reverse communication implementation
330 !****************************************
331 
332  10 CONTINUE
333 
334  IF (lactive) THEN
335  CALL drive_dgmres(n, nloc, m, lwork, work, irc,
336  & gi%icntl, gi%cntl, gi%info, rinfo)
337  ELSE
338  irc = 1
339  ENDIF
340 
341 !BCast to all processors in the local world group
342  CALL mpi_bcast(irc(1), 1, mpi_integer, 0, my_comm_world,
343  & mpi_err)
344  CALL mpi_bcast(rinfo(1), 1, mpi_real8, 0, my_comm_world,
345  & mpi_err)
346 
347  revcom = irc(1)
348  colx = irc(2)
349  coly = irc(3)
350  colz = irc(4)
351  nbscal = irc(5)
352 
353  revcom_loop: IF (revcom .EQ. matveci) THEN
354 ! perform the matrix vector product work(colz) <-- A * work(colx)
355 ! WRITE(10000+iam,*) "matvec"; CALL FLUSH(10000+iam)
356  CALL yax (work(colx), work(colz), nloc)
357  GOTO 10
358 
359  ELSE IF (revcom .EQ. precondleft) THEN
360 ! IF (lprint) PRINT *,'CALL PRECONL'
361 ! perform the left preconditioning work(colz) <-- M^{-1} * work(colx)
362 ! WRITE(10000+iam,*) "left_dcopy"; CALL FLUSH(10000+iam)
363  IF (lactive) CALL dcopy(nloc,work(colx),1,work(colz),1)
364  GOTO 10
365 
366  ELSE IF (revcom .EQ. precondright) THEN
367  CALL dcopy(nloc,work(colx),1,work(colz),1)
368 
369  IF (lactive) THEN
370  aux(myrowstart:myrowend) = work(colz:colz+nloc-1)
371 ! WRITE(10000+iam,*) "precondRight"; CALL FLUSH(10000+iam)
372  CALL apply_precond(aux)
373  work(colz:colz+nloc-1)=aux(myrowstart:myrowend)
374  END IF
375  GOTO 10
376 
377  ELSE IF (revcom .EQ. dotprod) THEN
378 ! perform the scalar product (uses nbscal columns of A starting at work(colx))
379  ! work(colz) <-- work(colx) work(coly)
380 
381  lactive_dp: IF (lactive) THEN
382  !MAKE SURE nbscal is the same on all processors -
383  CALL mpi_allgather(nbscal, 1, mpi_integer, itest, 1,
384  & mpi_integer, my_comm, mpi_err)
385  IF (lprint .AND. any(itest(1:gi%nprocs) .NE. nbscal)) THEN
386  print *,'itest: ',itest(1:gi%nprocs)
387  stop 'nbscal not same on all procs!'
388  END IF
389 
390 ! THIS IS FASTER THAN ALLGATHERV LOOP, BUT MAY LEAD TO SLIGHTLY DIFFERENT CONVERGENCE
391 ! SEQUENCES FOR DIFFERENT # PROCESSORS. THE DGEMV CALL IS SLIGHTLY SLOWER THAN THE LOOP
392  IF (lfast) THEN
393 
394 !! CALL dgemv('C',nloc,nbscal,one,work(colx),nloc,work(coly),1,zero,aux,1)
395 
396  DO nout=1,nbscal
397  aux(nout)= sum(work(colx:colx+nloc-1)*work(coly:coly+nloc-1))
398  colx = colx+nloc
399  END DO
400 
401  CALL mpi_allreduce(aux,work(colz),nbscal,mpi_real8,mpi_sum,
402  & my_comm,mpi_err)
403 
404  ELSE
405 
406 ! THE TIMING FOR THIS DOES NOT SCALE AS WELL WITH # PROCESSORS, BUT
407 ! LEADS TO A CONVERGENCE SEQUENCE THAT IS INDEPENDENT OF # PROCESSORS
408  CALL mpi_gatherv(work(coly), nloc, mpi_real8, aux, gi%rcounts,
409  & gi%disp, mpi_real8, 0, my_comm, mpi_err)
410 
411  DO nout=0,nbscal-1
412  CALL mpi_gatherv(work(colx),nloc,mpi_real8,tmpbuf,gi%rcounts,
413  & gi%disp,mpi_real8,0,my_comm,mpi_err)
414  IF (iam .EQ. 0) work(colz+nout) = sum(tmpbuf*aux) !DO FOR ALL PROCS IF ALLGATHER FORM USED
415  colx = colx+nloc
416  END DO
417 
418  CALL mpi_bcast(work(colz), nbscal, mpi_real8, 0, my_comm,
419  & mpi_err)
420  END IF
421 
422  DO nout = 0,nbscal-1
423  CALL truncate(work(colz+nout), 15)
424  END DO
425  END IF lactive_dp
426 
427  GOTO 10
428 
429  ELSE IF (revcom .EQ. peek) THEN
430  IF (lactive) THEN
431  fsq_last = fsq_nl !need for delfsq criteria
432  aux(myrowstart:myrowend) = work(colx:colx+nloc-1)
433  END IF
434 !get nonlinear force
435  IF (l_nonlinear) CALL getnlforce(aux, fsq_nl, bnorm)
436 
437  fsq_lin = (rinfo(1)*bnorm)**2
438  icount = icount+1
439  lstore = (fsq_nl.LT.fsq_min .OR. icount.EQ.1)
440  fsq_min = min(fsq_min, fsq_nl)
441  fsq_min_lin = min(fsq_min_lin, fsq_lin)
442 
443  IF (fsq_lin .LT. 1.e-30_dp) gi%icntl(7)=gi%info(1)
444 
445  IF (lprint .AND. icount.EQ.1) THEN
446  IF (l_nonlinear) THEN
447  print 905
448  ELSE
449  print 906
450  END IF
451  END IF
452  kout = gi%info(1)
453 
454  lactive1: IF (lactive) THEN
455 
456  IF (icount.EQ.1 .OR. mod(icount,10).EQ.0) THEN
457  lsum=sum(aux(myrowstart:myrowend)**2)
458  CALL mpi_reduce(lsum, gsum, 1, mpi_real8, mpi_sum, 0,
459  1 my_comm, mpi_err)
460  xmod=bnorm*sqrt(gsum)
461  lsum=maxval(abs(aux(myrowstart:myrowend)))
462  CALL mpi_reduce(lsum, xmax, 1, mpi_real8, mpi_max, 0,
463  1 my_comm, mpi_err)
464  xmax=bnorm*xmax
465  kprint = kout
466  IF (lprint) THEN
467  IF (l_nonlinear) THEN
468  print 910, kout, fsq_nl, xmod, xmax, fsq_lin
469  ELSE
470  print 911, kout, xmod, xmax, fsq_lin
471  END IF
472  END IF
473  END IF
474 
475  END IF lactive1
476 
477  IF (gi%ngmres_type .EQ. 1) GOTO 10 !OLDSTYLE:IGNORE LOGIC BELOW
478 
479 ! LACTIVE1: IF (lactive) THEN
480  delfsq = (fsq_last-fsq_nl)/fsq_min
481  IF (delfsq .LT. 0.05_dp) THEN !PROGRESS COUNTER
482  jcount = jcount+1
483  ELSE
484  jcount = 0
485  END IF
486 
487 ! STOPPING CRITERIA (reset max iterations to current iteration)
488  IF (fsq_nl.GT.(3*fsq_min) .OR. jcount.GT.3
489  & .OR. fsq_nl.LE.gi%ftol) gi%icntl(7) = gi%info(1)
490  IF (lstore) THEN
491  IF (.NOT.ALLOCATED(xmin)) ALLOCATE (xmin(nloc))
492  xmin = work(colx:colx+nloc-1)
493  END IF
494  GOTO 10
495 
496  ENDIF revcom_loop
497 
498  905 FORMAT(1x,'ITER',7x,'FSQ_NL',10x,'||X||',9x,'MAX|X|',9x,'FSQ_ARN')
499  906 FORMAT(1x,'ITER',7x,'||X||',11x,'MAX|X|',9x,'FSQ_ARN')
500  910 FORMAT(i5, 4(3x,1pe12.3))
501  911 FORMAT(i5, 3(3x,1pe12.3))
502 
503 !*******************************
504 ! end reverse loop: dump the solution to a file for debugging
505 !******************************
506  GOTO 100
507 
508  lactive2: IF (lprint .AND. lactive) THEN
509 
510  nout = 11
511  OPEN(nout,file='sol_dTestgmres',status='unknown')
512  IF (gi%icntl(5).EQ.0) THEN
513  WRITE(nout,*) 'Orthogonalisation : MGS'
514  ELSEIF (gi%icntl(5).eq.1) THEN
515  WRITE(nout,*) 'Orthogonalisation : IMGS'
516  ELSEIF (gi%icntl(5).eq.2) THEN
517  WRITE(nout,*) 'Orthogonalisation : CGS'
518  ELSEIF (gi%icntl(5).eq.3) THEN
519  WRITE(nout,*) 'Orthogonalisation : ICGS'
520  ENDIF
521  WRITE(nout,*) 'Restart : ', m
522  WRITE(nout,*) 'info(1) = ',gi%info(1),' info(2) = ',gi%info(2)
523  WRITE(nout,*) 'rinfo(1) = ',rinfo(1),' rinfo(2) = ',rinfo(2)
524  WRITE(nout,*) 'Optimal workspace = ', gi%info(3)
525  WRITE(nout,*) 'Solution : '
526  DO jcount=1,n
527  WRITE(nout,*) work(jcount)
528  ENDDO
529  WRITE(nout,*) ' '
530 
531  END IF lactive2
532 
533  100 CONTINUE
534 
535  lactive3: IF (lactive) THEN
536  IF (ALLOCATED (xmin)) THEN ! otherwise, work(1:nloc) holds final state
537  work(1:nloc) = xmin(1:nloc)
538  DEALLOCATE (xmin)
539  ELSE
540  fsq_min=fsq_nl
541  END IF
542 
543  CALL mpi_allgatherv(work, nloc, mpi_real8, tmpbuf, gi%rcounts,
544  & gi%disp, mpi_real8, my_comm, mpi_err)
545  x0 = tmpbuf*bnorm
546 
547  IF (lprint .AND. kout.NE.kprint) THEN
548  xmod = sqrt(sum(x0*x0)); xmax = maxval(abs(x0))
549  IF (l_nonlinear) THEN
550  print 910, kout, fsq_min, xmod, xmax, fsq_min_lin
551  ELSE
552  print 911, kout, xmod, xmax, fsq_lin
553  END IF
554  END IF
555 
556  gi%ftol = fsq_min
557 
558 !Sanity check
559 ! CALL GetNLForce(tmpbuf, fsq_nl, bnorm)
560 ! IF (iam .EQ. 0) PRINT *,'FSQ_MIN: ',fsq_min,' FSQ_NL: ', fsq_nl
561 
562  END IF lactive3
563 
564 
565  DEALLOCATE (work, tmpbuf, aux, itest)
566 #endif
567  END SUBROUTINE gmres_par
568 
569  SUBROUTINE truncate(num, iprec0)
570  USE stel_kinds, ONLY: dp
571 !-----------------------------------------------
572 ! D u m m y A r g u m e n t s
573 !-----------------------------------------------
574  INTEGER, INTENT(IN) :: iprec0
575  REAL(dp), INTENT(INOUT) :: num
576 !-----------------------------------------------
577 ! L o c a l V a r i a b l e s
578 !-----------------------------------------------
579  CHARACTER*24 :: chnum, tchnum
580 !-----------------------------------------------
581 !
582 ! TRUNCATES double-precision to precision iprec0 digits, keeping exponent range of double
583 ! WRITE TO INTERNAL FILE TO DO TRUNCATION
584 !
585 ! RETURN
586 
587  WRITE (chnum, '(a,i2,a,i2,a)') '(1p,e',iprec0+7,'.',iprec0,')'
588  WRITE (tchnum, chnum) num
589 
590  READ (tchnum, chnum) num
591 
592  END SUBROUTINE truncate
593 
594 
595  END MODULE gmres_lib
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
gmres_lib::gmres_info
Definition: gmres_lib.f:7
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17