OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
imp_lanz.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| imp_lanzp ../engine/source/implicit/imp_lanz.F
25!||--- called by ------------------------------------------------------
26!|| lin_solvh0 ../engine/source/implicit/lin_solv.F
27!|| lin_solvh1 ../engine/source/implicit/lin_solv.F
28!|| lin_solvih2 ../engine/source/implicit/lin_solv.F
29!||--- calls -----------------------------------------------------
30!|| mav_ltp ../engine/source/implicit/produt_v.F
31!|| prec_solvp ../engine/source/implicit/prec_solv.F
32!|| produt_w ../engine/source/implicit/produt_v.F
33!||--- uses -----------------------------------------------------
34!|| dsgraph_mod ../engine/share/modules/dsgraph_mod.F
35!|| groupdef_mod ../common_source/modules/groupdef_mod.f
36!|| intbufdef_mod ../common_source/modules/interfaces/intbufdef_mod.f90
37!||====================================================================
38 SUBROUTINE imp_lanzp(IPREC ,
39 1 N ,NNZ ,IADK ,JDIK ,DIAG_K ,
40 2 LT_K ,NI ,ITOK ,IADI ,JDII ,
41 3 LT_I ,NNZM ,IADM ,JDIM ,DIAG_M ,
42 3 LT_M ,X ,R ,ITOL ,RTOL ,
43 4 V ,W ,Y ,ITASK ,IPRINT ,
44 5 SHIFT ,KCOND ,N_MAX ,FLM ,F_X ,
45 6 ISTOP ,W_DDL,A ,AR ,
46 9 VE ,MS ,XE ,D ,DR ,
47 A NDOF ,IPARI ,INTBUF_TAB ,NUM_IMP,
48 B NS_IMP,NE_IMP,NSREM ,NSL ,NMONV ,
49 C IMONV ,MONVOL ,IGRSURF ,VOLMON,
50 D FR_MV ,IBFV ,SKEW ,XFRAME,IND_IMP,
51 H XI_C ,R0 ,NDDLI_G,INTP_C,IRBE3 ,
52 E LRBE3 ,IRBE2 ,LRBE2 )
53C-----------------------------------------------
54C M o d u l e s
55C-----------------------------------------------
56 USE dsgraph_mod
57 USE intbufdef_mod
58 USE groupdef_mod
59C-----------------------------------------------
60C I m p l i c i t T y p e s
61C-----------------------------------------------
62#include "implicit_f.inc"
63C-----------------------------------------------
64C C o m m o n B l o c k s
65C-----------------------------------------------
66#if defined(MUMPS5)
67#include "dmumps_struc.h"
68#endif
69#include "com04_c.inc"
70#include "units_c.inc"
71C-----------------------------------------------
72C D u m m y A r g u m e n t s
73C-----------------------------------------------
74C----------resol [K]{X}={F}---------
75 INTEGER N ,NNZ ,IADK(*) ,JDIK(*),ITOL,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(*)
83C REAL
84 my_real
85 . DIAG_K(*), LT_K(*) ,DIAG_M(*),LT_M(*) ,X(*), RTOL,
86 . V(*) ,W(*) ,R(*) ,Y(*),SHIFT,KCOND,LT_I(*),FLM,F_X
87 my_real
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
94C-----------------------------------------------
95C L o c a l V a r i a b l e s
96C-----------------------------------------------
97 INTEGER I,J,ITN,IP,NLIM,IBID,IUPD
98 my_real
99 . r2(n),b(n),cs1(2),rr,
100 . r02,anorm,rnorm, ynorm,rbid
101 my_real
102 . ddot, dnrm2,
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
130C*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
131C* based on File symmlq.f
132C*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
133
134C machine precision minimum dependant du type REAL ou DOUBLE
135c write(*,*)'flm=',flm
136C-------add if (RTOL<FLM) RTOL=FLM
137 iupd = 0
138 nlim=n_max
139C* Print heading and initialize.
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
157C* Set up y for the first Lanczos vector v1.
158C* y is really beta1 * P * v1 where P = C**(-1).
159C* y and beta1 will be zero if b = 0.
160
161C --R=B-A.X--- ---
162 CALL mav_ltp(
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
176 CALL prec_solvp(iprec,itask ,
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 )
184C IF ( GOODB ) THEN
185C B1 = Y(1)
186C ELSE
187 b1 = zero
188C ENDIF
189 CALL produt_w( n ,r ,y ,w_ddl,beta1)
190
191C* TEST FOR AN INDEFINITE PRECONDITIONER.
192
193 IF (beta1 < zero) THEN
194c ISTOP = 8
195c GO TO 900
196 beta1 = -beta1
197 WRITE(iout, 3000) warr, msg(8)
198 WRITE(istdo, 3000) warr, msg(8)
199 ENDIF
200
201C* IF B = 0 EXACTLY, STOP WITH X = 0.
202
203 IF (beta1 == zero) THEN
204 istop=-1
205 GOTO 900
206 ENDIF
207
208C* HERE AND LATER, V IS REALLY P * (THE LANCZOS V).
209
210 IF (itol==2) r02 = beta1
211 beta1 = sqrt( beta1 )
212 s = one / beta1
213C
214 IF (itol==1) THEN
215 CALL produt_w(n,r,r,w_ddl,r02)
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
228 CALL mav_ltp(
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
239C* SET UP Y FOR THE SECOND LANCZOS VECTOR.
240C* AGAIN, Y IS BETA * P * V2 WHERE P = C**(-1).
241C* Y AND BETA WILL BE ZERO OR VERY SMALL IF B IS AN EIGENVECTOR.
242
243 IF (shift/=zero) THEN
244 DO i = 1, n
245 y(i) = y(i)-shift*v(i)
246 ENDDO
247 ENDIF
248 CALL produt_w( n ,v ,y ,w_ddl,alfa)
249 t=alfa / beta1
250 DO i = 1, n
251 y(i) = y(i)-t*r(i)
252 ENDDO
253
254C* MAKE SURE R2 WILL BE ORTHOGONAL TO THE FIRST V.
255
256 CALL produt_w( n ,v ,y ,w_ddl,z)
257 CALL produt_w( n ,v ,v ,w_ddl,s)
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
266 CALL prec_solvp(iprec,itask ,
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
278c ISTOP = 8
279c GO TO 900
280 beta=-beta
281 WRITE(iout, 3000) warr, msg(8)
282 WRITE(istdo, 3000) warr, msg(8)
283 ENDIF
284
285C* CAUSE TERMINATION (LATER) IF BETA IS ESSENTIALLY ZERO.
286
287 beta = sqrt( beta )
288 IF (beta <= flm) istop = -1
289
290C* SEE IF THE LOCAL REORTHOGONALIZATION ACHIEVED ANYTHING.
291
292C CALL PRODUT_V( N ,R2 ,R2 ,T)
293C DENOM = ONE/(SQRT( S*T ) + FLM)
294C S = Z * DENOM
295C CALL PRODUT_V( N ,V ,R2 ,T)
296C T = T * DENOM
297
298C* INITIALIZE OTHER QUANTITIES.
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
312C IF ( GOODB ) THEN
313C DO I = 1, N
314C W(I) = ZERO
315C ENDDO
316C ELSE
317 DO i = 1, n
318 w(i) = v(i)
319 ENDDO
320C ENDIF
321
322C* ------------------------------------------------------------------
323C* MAIN ITERATION LOOP.
324C* ------------------------------------------------------------------
325
326C* ESTIMATE VARIOUS NORMS AND TEST FOR CONVERGENCE.
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
341C* ESTIMATE COND(A).
342C* IN THIS VERSION WE LOOK AT THE DIAGONALS OF L IN THE
343C* FACTORIZATION OF THE TRIDIAGONAL MATRIX, T = L*Q.
344C* SOMETIMES, T(K) CAN BE MISLEADINGLY ILL-CONDITIONED WHEN
345C* T(K+1) IS NOT, SO WE MUST BE CAREFUL NOT TO OVERESTIMATE KCOND.
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
354C* SEE IF ANY OF THE STOPPING CRITERIA ARE SATISFIED.
355C* IN RARE CASES, ISTOP IS ALREADY -1 FROM ABOVE (ABAR = CONST * I).
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
364C* ==================================================================
365
366C* SEE IF IT IS TIME TO PRINT SOMETHING.
367
368C* IF (NOUT <= 0) GO TO 600
369C* IF (N <= 40) GO TO 400
370C* IF (ITN <= 10) GO TO 400
371C* IF (ITN >= ITNLIM - 10) GO TO 400
372C* IF (MOD(ITN,10) == 0) GO TO 400
373C* IF (CGNORM <= TEN*EPSX) GO TO 400
374C* IF (CGNORM <= TEN*EPSR) GO TO 400
375C* IF (KCOND >= EM2/FLM ) GO TO 400
376C* IF (ISTOP /= 0) GO TO 400
377C* GOTO 600
378C*
379C* PRINT A LINE FOR THIS ITERATION.
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
388C* 400 ZBAR = RHS1 / DIAG
389C Z = (SNPROD * ZBAR + BSTEP) / BETA1
390C X1LQ = X(1) + B1 * BSTEP / BETA1
391C X1CG = X(1) + W(1) * ZBAR + B1 * Z
392C
393C IF ( ITN == 0) WRITE(NOUT, 1200)
394C WRITE(NOUT, 1300) ITN, X1CG, CGNORM, BSTEP/BETA1, ANORM, KCOND
395C IF (MOD(ITN,10) == 0) WRITE(NOUT, 1500)
396C* ==================================================================
397
398
399C* OBTAIN THE CURRENT LANCZOS VECTOR V = (1 / BETA)*Y
400C* AND SET UP Y FOR THE NEXT ITERATION.
401
402 IF (istop /= 0) GO TO 800
403C------place CRIT_STOP ICI-----
404 s = one / beta
405
406 DO i = 1, n
407 v(i) = s * y(i)
408 ENDDO
409
410 CALL mav_ltp(
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
429 CALL produt_w( n ,v ,y ,w_ddl,alfa)
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
440 CALL prec_solvp(iprec,itask ,
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
451c ISTOP = 6
452c GOTO 800
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
460C* COMPUTE THE NEXT PLANE ROTATION FOR Q.
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
470C* UPDATE X.
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
481C* ACCUMULATE THE STEP ALONG THE DIRECTION B,
482C* AND GO ROUND AGAIN.
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
491c WRITE(IOUT,2001)
492c . ANORM, SQRT(YNORM2),GMAX/GMIN
493c WRITE(IOUT,2002)
494c . ALFA, BETA,OLDB
495 GO TO 300
496
497C* ------------------------------------------------------------------
498C* END OF MAIN ITERATION LOOP.
499C* ------------------------------------------------------------------
500
501C* MOVE TO THE CG POINT IF IT SEEMS BETTER.
502C* IN THIS VERSION OF SYMMLQ, THE CONVERGENCE TESTS INVOLVE
503C* ONLY CGNORM, SO WE'RE UNLIKELY TO STOP AT AN LQ POINT,
504C* EXCEPT IF THE ITERATION LIMIT INTERFERES.
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
518C* IF ( GOODB ) THEN
519
520C* ADD THE STEP ALONG B.
521
522C* BSTEP = BSTEP / BETA1
523C* DO I = 1, N
524C* Y(I) = B(I)
525C* ENDDO
526C* CALL PREC_SOLV(IPREC,
527C* 1 N ,NNZM,IADM ,JDIM ,DIAG_M ,
528C* 2 LT_M ,B ,Y )
529C* DO I = 1, N
530C* X(I) = X(I)+BSTEP*Y(I)
531C* ENDDO
532C* END IF
533
534C* ==================================================================
535C* DISPLAY FINAL STATUS.
536C* ==================================================================
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
550C --R=B-A.X--- ---
551 CALL mav_ltp(
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
565 CALL produt_w( n ,r ,r ,w_ddl,r02)
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
573C--------POUR F*X--------
574 RETURN
575
576C* ------------------------------------------------------------------
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)
596C* ------------------------------------------------------------------
597C* END OF SYMMLQ
598#endif
599 END
ddot
double precision function ddot(n, dx, incx, dy, incy)
DDOT
Definition ddot.f:82
dnrm2
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
imp_lanzp
subroutine imp_lanzp(iprec, n, nnz, iadk, jdik, diag_k, lt_k, ni, itok, iadi, jdii, lt_i, nnzm, iadm, jdim, diag_m, lt_m, x, r, itol, rtol, v, w, y, itask, iprint, shift, kcond, n_max, flm, f_x, istop, w_ddl, a, ar, ve, ms, xe, d, dr, ndof, ipari, intbuf_tab, num_imp, ns_imp, ne_imp, nsrem, nsl, nmonv, imonv, monvol, igrsurf, volmon, fr_mv, ibfv, skew, xframe, ind_imp, xi_c, r0, nddli_g, intp_c, irbe3, lrbe3, irbe2, lrbe2)
Definition imp_lanz.F:53
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
dsgraph_mod
Definition dsgraph_mod.F:31
groupdef_mod
Definition groupdef_mod.F:662
prec_solvp
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)
Definition prec_solv.F:147
mav_ltp
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)
Definition produt_v.F:566
produt_w
subroutine produt_w(nddl, x, y, w, r)
Definition produt_v.F:106
dsgraph_mod::prgraph
Definition dsgraph_mod.F:45