53
54
55
57 USE intbufdef_mod
59
60
61
62#include "implicit_f.inc"
63
64
65
66#if defined(MUMPS5)
67#include "dmumps_struc.h"
68#endif
69#include "com04_c.inc"
70#include "units_c.inc"
71
72
73
74
75 INTEGER N ,NNZ ,IADK(*) ,JDIK(*),,ISTOP,
76 . NI ,ITOK(*) ,IADI(*),JDII(*),N_MAX,
77 . NNZM ,IADM(*),JDIM(*),IPREC,ITASK,IPRINT,
78 . W_DDL(*),IBFV(*),IRBE3(*) ,LRBE3(*),NDDLI_G,
79 . IRBE2(*),LRBE2(*)
80 INTEGER NDOF(*),NE_IMP(*),NSREM ,NSL,INTP_C,
81 . IPARI(*) ,NUM_IMP(*),NS_IMP(*),IND_IMP(*)
82 INTEGER NMONV,IMONV(*),MONVOL(*),FR_MV(*)
83
85 . diag_k(*), lt_k(*) ,diag_m(*),lt_m(*) ,x(*), rtol,
86 . v(*) ,w(*) ,r(*) ,y(*),shift,kcond,lt_i(*),flm,f_x
88 . a(3,*),ar(3,*),ve(3,*),d(3,*),dr(3,*),xe(3,*),
89 . ms(*),volmon(*),skew(*) ,xframe(*),xi_c(*),r0
90
91 TYPE(INTBUF_STRUCT_) INTBUF_TAB(*)
92 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
93#ifdef MUMPS5
94
95
96
97 INTEGER I,J,ITN,IP,NLIM,IBID,IUPD
99 . r2(n),b(n),cs1(2),rr,
100 . r02,anorm,rnorm, ynorm,rbid
103 . alfa, b1, beta, beta1, bstep, cs,
104 . cgnorm, dbar, delta, denom, diag,
105 . eps, epsa, epsln, epsr, epsx,
106 . gamma, gbar, gmax, gmin, gpert,
107 . lqnorm, oldb, qrnorm, rhs1, rhs2,
108 . s, sn, snprod, t, tnorm,
109 . x1cg, x1lq, ynorm2, zbar, z
110 CHARACTER*16 EXIT
111 CHARACTER*11 WARR
112 CHARACTER*52 MSG(-1:8)
113
114 DATA EXIT /'TERMINATION WITH'/
115 DATA warr /'**WARRING**'/
116
117 DATA msg
118 . / 'BETA2 = 0. IF M = I, F AND X ARE EIGENVECTORS OF K',
119 . 'BETA1 = 0. THE EXACT SOLUTION IS X = 0',
120 . 'REQUESTED ACCURACY ACHIEVED, AS DETERMINED BY RTOL',
121 . 'REASONABLE ACCURACY ACHIEVED, GIVEN EPS',
122 . 'X HAS CONVERGED TO AN EIGENVECTOR',
123 . 'ACOND HAS EXCEEDED 0.1/EPS',
124 . 'THE ITERATION LIMIT WAS REACHED',
125 . 'APROD DOES NOT DEFINE A SYMMETRIC MATRIX',
126 . 'MSOLVE DOES NOT DEFINE A SYMMETRIC MATRIX',
127 . 'MSOLVE DOES NOT DEFINE A POS-DEF PRECONDITIONER' /
128 TYPE(PRGRAPH) :: GBID
129 TYPE(DMUMPS_STRUC) MBID
130
131
132
133
134
135
136
137 iupd = 0
138 nlim=n_max
139
140
141 IF(iprint/=0)THEN
142 WRITE(iout,*)' *** BEGIN PRECONDITION LANCZOS ITERATION ***'
143 WRITE(iout,*)
144 ENDIF
145 IF(iprint<0)THEN
146 WRITE(istdo,*)' *** BEGIN PRECONDITION LANCZOS ITERATION ***'
147 WRITE(istdo,*)
148 ENDIF
149 ip = iabs(iprint)
150 istop = 0
151 itn = 0
152 anorm = zero
153 kcond = zero
154 rnorm = zero
155 ynorm = zero
156
157
158
159
160
161
163 1 n ,ni ,iadk ,jdik ,diag_k,
164 2 lt_k ,iadi ,jdii ,itok ,lt_i ,
165 3 x ,y ,a ,ar ,
166 5 ve ,ms ,xe ,d ,dr ,
167 6 ndof ,ipari ,intbuf_tab ,num_imp,
168 7 ns_imp,ne_imp,nsrem ,nsl ,ibfv ,
169 8 skew ,xframe,monvol,volmon,igrsurf ,
170 9 fr_mv,nmonv ,imonv ,ind_imp,
171 a xi_c ,iupd ,irbe3 ,lrbe3 ,irbe2 ,lrbe2 )
172 DO i = 1, n
173 b(i) = r(i)
174 r(i) = r(i) - y(i)
175 ENDDO
177 1 gbid ,ibid ,ibid , diag_k,lt_k ,
178 2 iadk ,jdik ,ibid ,ibid ,ibid ,
179 3 ibid ,rbid ,ibid ,ibid ,mbid ,
180 4 ibid ,ibid ,ibid ,ibid ,ibid ,
181 5 ibid ,ibid ,
182 1 n ,nnzm,iadm ,jdim ,diag_m ,
183 2 lt_m ,r ,y )
184
185
186
187 b1 = zero
188
189 CALL produt_w( n ,r ,y ,w_ddl,beta1)
190
191
192
193 IF (beta1 < zero) THEN
194
195
196 beta1 = -beta1
197 WRITE(iout, 3000) warr, msg(8)
198 WRITE(istdo, 3000) warr, msg(8)
199 ENDIF
200
201
202
203 IF (beta1 == zero) THEN
204 istop=-1
205 GOTO 900
206 ENDIF
207
208
209
210 IF (itol==2) r02 = beta1
211 beta1 = sqrt( beta1 )
212 s = one / beta1
213
214 IF (itol==1) THEN
216 rr = sqrt(r02)
217 ELSE
218 rr = beta1
219 ENDIF
220 IF(iprint/=0)THEN
221 IF(iprint<0) WRITE(istdo,1000)itn,rr
222 WRITE(iout,1000)itn,rr
223 ENDIF
224 DO i = 1, n
225 v(i) = s * y(i)
226 ENDDO
227
229 1 n ,ni ,iadk ,jdik ,diag_k,
230 2 lt_k ,iadi ,jdii ,itok ,lt_i ,
231 3 v ,y ,a ,ar ,
232 5 ve ,ms ,xe ,d ,dr ,
233 6 ndof ,ipari ,intbuf_tab ,num_imp,
234 7 ns_imp,ne_imp,nsrem ,nsl ,ibfv ,
235 8 skew ,xframe,monvol,volmon,igrsurf ,
236 9 fr_mv,nmonv ,imonv ,ind_imp,
237 a xi_c ,iupd ,irbe3 ,lrbe3 ,irbe2 ,lrbe2 )
238
239
240
241
242
243 IF (shift/=zero) THEN
244 DO i = 1, n
245 y(i) = y(i)-shift*v(i)
246 ENDDO
247 ENDIF
249 t=alfa / beta1
250 DO i = 1, n
251 y(i) = y(i)-t*r(i)
252 ENDDO
253
254
255
258 t = z / s
259 DO i = 1, n
260 y(i) = y(i)-t*v(i)
261 ENDDO
262
263 DO i = 1, n
264 r2(i) = y(i)
265 ENDDO
267 1 gbid ,ibid ,ibid , diag_k,lt_k ,
268 2 iadk ,jdik ,ibid ,ibid ,ibid ,
269 3 ibid ,rbid ,ibid ,ibid ,mbid ,
270 4 ibid ,ibid ,ibid ,ibid ,ibid ,
271 5 ibid ,ibid ,
272 1 n ,nnzm,iadm ,jdim ,diag_m ,
273 2 lt_m ,r2 ,y )
274 oldb = beta1
275 CALL produt_w( n ,r2 ,y ,w_ddl,beta)
276
277 IF (beta < zero) THEN
278
279
280 beta=-beta
281 WRITE(iout, 3000) warr, msg(8)
282 WRITE(istdo, 3000) warr, msg(8)
283 ENDIF
284
285
286
287 beta = sqrt( beta )
288 IF (beta <= flm) istop = -1
289
290
291
292
293
294
295
296
297
298
299
300 cgnorm = beta1
301 gbar = alfa
302 dbar = beta
303 rhs1 = beta1
304 rhs2 = zero
305 bstep = zero
306 snprod = one
307 tnorm = alfa*alfa + beta*beta
308 ynorm2 = zero
309 gmax = abs( alfa ) + flm
310 gmin = gmax
311
312
313
314
315
316
317 DO i = 1, n
318 w(i) = v(i)
319 ENDDO
320
321
322
323
324
325
326
327
328 300 itn = itn + 1
329 anorm = sqrt( tnorm )
330 ynorm = sqrt( ynorm2 )
331 epsa = anorm * flm
332 epsx = anorm * ynorm * flm
333 epsr = anorm * ynorm * rtol
334 diag = gbar
335 IF (diag == zero) diag = epsa
336
337 lqnorm = sqrt( rhs1*rhs1 + rhs2*rhs2 )
338 qrnorm = snprod * beta1
339 cgnorm = qrnorm * beta / abs( diag )
340
341
342
343
344
345
346
347 IF (lqnorm <= cgnorm) THEN
348 kcond = gmax / gmin
349 ELSE
350 denom =
min( gmin, abs( diag ) )
351 kcond = gmax / denom
352 ENDIF
353
354
355
356
357 IF (istop == 0) THEN
358 IF (itn >= nlim ) istop = 5
359 IF (kcond >=em01/flm) istop = 4
360 IF (epsx >= beta1 ) istop = 3
361 IF (cgnorm <= epsx ) istop = 2
362 IF (cgnorm <= epsr ) istop = 1
363 ENDIF
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381 IF(iprint/=0)THEN
382 IF (mod(itn,ip)==0)THEN
383 WRITE(iout,1001)itn,cgnorm
384 IF(iprint<0) WRITE(istdo,1001)itn,cgnorm
385 ENDIF
386 ENDIF
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402 IF (istop /= 0) GO TO 800
403
404 s = one / beta
405
406 DO i = 1, n
407 v(i) = s * y(i)
408 ENDDO
409
411 1 n ,ni ,iadk ,jdik ,diag_k,
412 2 lt_k ,iadi ,jdii ,itok ,lt_i ,
413 3 v ,y ,a ,ar ,
414 5 ve ,ms ,xe ,d ,dr ,
415 6 ndof ,ipari ,intbuf_tab ,num_imp,
416 7 ns_imp,ne_imp,nsrem ,nsl ,ibfv ,
417 8 skew ,xframe,monvol,volmon,igrsurf ,
418 9 fr_mv,nmonv ,imonv ,ind_imp,
419 a xi_c ,iupd ,irbe3 ,lrbe3 ,irbe2 ,lrbe2 )
420 IF (shift/=zero) THEN
421 DO i = 1, n
422 y(i) = y(i)-shift*v(i)
423 ENDDO
424 ENDIF
425 t=beta / oldb
426 DO i = 1, n
427 y(i) = y(i)-t*r(i)
428 ENDDO
430 t = alfa / beta
431 DO i = 1, n
432 y(i) = y(i)-t*r2(i)
433 ENDDO
434 DO i = 1, n
435 r(i) = r2(i)
436 ENDDO
437 DO i = 1, n
438 r2(i) = y(i)
439 ENDDO
441 1 gbid ,ibid ,ibid , diag_k,lt_k ,
442 2 iadk ,jdik ,ibid ,ibid ,ibid ,
443 3 ibid ,rbid ,ibid ,ibid ,mbid ,
444 4 ibid ,ibid ,ibid ,ibid ,ibid ,
445 5 ibid ,ibid ,
446 1 n ,nnzm,iadm ,jdim ,diag_m ,
447 2 lt_m ,r2 ,y )
448 oldb = beta
449 CALL produt_w( n ,r2 ,y ,w_ddl,beta)
450 IF (beta < zero) THEN
451
452
453 beta=-beta
454 WRITE(iout, 3000) warr, msg(8)
455 WRITE(istdo, 3000) warr, msg(8)
456 ENDIF
457 beta = sqrt( beta )
458 tnorm = tnorm + alfa**2 + oldb**2 + beta**2
459
460
461
462 gamma = sqrt( gbar*gbar + oldb*oldb )
463 cs = gbar / gamma
464 sn = oldb / gamma
465 delta = cs * dbar + sn * alfa
466 gbar = sn * dbar - cs * alfa
467 epsln = sn * beta
468 dbar = - cs * beta
469
470
471
472 z = rhs1 / gamma
473 s = z * cs
474 t = z * sn
475
476 DO i = 1, n
477 x(i) = (w(i) * s + v(i) * t) + x(i)
478 w(i) = w(i) * sn - v(i) * cs
479 ENDDO
480
481
482
483
484 bstep = snprod * cs * z + bstep
485 snprod = snprod * sn
486 gmax =
max( gmax, gamma )
487 gmin =
min( gmin, gamma )
488 ynorm2 = z*z + ynorm2
489 rhs1 = rhs2 - delta * z
490 rhs2 = - epsln * z
491
492
493
494
495 GO TO 300
496
497
498
499
500
501
502
503
504
505
506 800 IF (cgnorm <= lqnorm) THEN
507 zbar = rhs1 / diag
508 bstep = snprod * zbar + bstep
509 ynorm = sqrt( ynorm2 + zbar*zbar )
510 rnorm = cgnorm
511 DO i = 1, n
512 x(i) = x(i)+zbar*w(i)
513 ENDDO
514 ELSE
515 rnorm = lqnorm
516 ENDIF
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537 900 CONTINUE
538 IF (iprint/=0) THEN
539 WRITE(iout,2000) EXIT, istop, itn,
540 . EXIT, anorm, kcond,
541 . EXIT, rnorm, ynorm
542 WRITE(iout, 3000) EXIT, msg(istop)
543 IF (iprint<0) THEN
544 WRITE(istdo,2000) EXIT, istop, itn,
545 . EXIT, anorm, kcond,
546 . EXIT, rnorm, ynorm
547 WRITE(istdo, 3000) EXIT, msg(istop)
548 ENDIF
549 ENDIF
550
552 1 n ,ni ,iadk ,jdik ,diag_k,
553 2 lt_k ,iadi ,jdii ,itok ,lt_i ,
554 3 x ,r ,a ,ar ,
555 5 ve ,ms ,xe ,d ,dr ,
556 6 ndof ,ipari ,intbuf_tab ,num_imp,
557 7 ns_imp,ne_imp,nsrem ,nsl ,ibfv ,
558 8 skew ,xframe,monvol,volmon,igrsurf ,
559 9 fr_mv,nmonv ,imonv ,ind_imp,
560 a xi_c ,iupd ,irbe3 ,lrbe3 ,irbe2 ,lrbe2 )
561 DO i = 1, n
562 r(i) = r(i) - b(i)
563 ENDDO
564 IF (itol==1) THEN
566 rnorm = sqrt(r02)/rr
567 ENDIF
568 IF (iprint/=0) THEN
569 WRITE(iout,1002)itn,rnorm
570 IF(iprint<0) WRITE(istdo,1002)itn,rnorm
571 ENDIF
572 IF (istop>0.AND.istop<4) istop=0
573
574 RETURN
575
576
577 1100 FORMAT(/ 1p, ' BETA1 =', e10.2, 3x, 'ALPHA1 =', e10.2
578 $ / ' (V1,V2) BEFORE AND AFTER ', e14.2
579 $ / ' LOCAL REORTHOGONALIZATION', e14.2)
580 1200 FORMAT(// 5x, 'ITN', 7x, 'X1(CG)', 10x,
581 $ 'NORM(R)', 5x, 'BSTEP', 7x, 'NORM(A)', 3x, 'COND(A)')
582 1300 FORMAT(1p, i8, e19.10, e11.2, e14.5, 2e10.2)
583 1500 FORMAT(1x)
584 2000 FORMAT(/ 1p, a, 6x, 'ISTOP =', i3, 15x, 'ITN =', i8
585 $ / a, 6x, 'ANORM =', e12.4, 6x, 'KCOND =', e12.4
586 $ / a, 6x, 'RNORM =', e12.4, 6x, 'YNORM =', e12.4)
587 3000 FORMAT( a, 6x, a )
588 1000 FORMAT(5x,'ITERATION=',i8,5x,' INITIAL RESIDUAL NORM=',e11.4)
589 1001 FORMAT(5x,'ITERATION=',i8,5x,' C.G. RESIDUAL NORM=',e11.4)
590 1002 FORMAT(3x,'TOTAL LANCZOS ITERATION=',i8,5x,
591 . ' RELATIVE RESIDUAL NORM=',e11.4)
592 2001 FORMAT(/ 5x, 'WITH', 2x, 'ANORM =', e12.4, 2x, 'YNORM =',
593 . e12.4,2x,'KCOND =', e12.4)
594 2002 FORMAT(/ 5x, 'WITH', 2x, 'ALFA =', e12.4, 2x, 'BETA =',
595 . e12.4,2x,'OLDB =', e12.4)
596
597
598#endif
double precision function ddot(n, dx, incx, dy, incy)
DDOT
real(wp) function dnrm2(n, x, incx)
DNRM2
subroutine prec_solvp(iprec, itask, graphe, iad_elem, fr_elem, diag_k, lt_k, iadk, jdik, itab, iprint, insolv, it, fac_k, ipiv_k, nk, mumps_par, cddlp, isolv, idsc, iddl, ikc, inloc, ndof, nddl, nnz, iadm, jdim, diag_m, lt_m, v, z)
subroutine mav_ltp(nddl, nddli, iadl, jdil, diag_k, lt_k, iadi, jdii, itok, lt_i, v, w, a, ar, ve, ms, x, d, dr, ndof, ipari, intbuf_tab, num_imp, ns_imp, ne_imp, nsrem, nsl, ibfv, skew, xframe, monvol, volmon, igrsurf, fr_mv, nmonv, imonv, index2, xi_c, iupd, irbe3, lrbe3, irbe2, lrbe2)
subroutine produt_w(nddl, x, y, w, r)