OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
imp_pcg.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/.
23C-------------------------
24!||====================================================================
25!|| crit_stop ../engine/source/implicit/imp_pcg.F
26!||--- called by ------------------------------------------------------
27!|| imp_pcg1 ../engine/source/implicit/imp_fsa_inv.F
28!|| imp_pcgh ../engine/source/implicit/imp_pcg.F
29!|| imp_ppcgh ../engine/source/implicit/imp_pcg.F
30!||====================================================================
31 INTEGER FUNCTION crit_stop(IT,R2,IT_MAX,TOL)
32C-----------------------------------------------
33C I m p l i c i t T y p e s
34C-----------------------------------------------
35#include "implicit_f.inc"
36C-----------------------------------------------
37C D u m m y A r g u m e n t s
38C-----------------------------------------------
39 INTEGER it,it_max
40 my_real
41 . r2,tol
42C-----------------------------------------------
43C L o c a l V a r i a b l e s
44C-----------------------------------------------
45C-----------------------------
46 IF (it>=it_max) THEN
47 crit_stop = 0
48 RETURN
49 ENDIF
50 IF (r2<=tol) THEN
51 crit_stop = 0
52 ELSE
53 crit_stop = 1
54 ENDIF
55C--------------------------------------------
56 RETURN
57 END
58!||====================================================================
59!|| imp_isave ../engine/source/implicit/imp_pcg.F
60!||--- called by ------------------------------------------------------
61!|| imp_pcgh ../engine/source/implicit/imp_pcg.F
62!||====================================================================
63 SUBROUTINE imp_isave(
64 1 NDDL,X,R,DIAG_K,DIAG_M,
65 2 NNZ,LT_K,LT_K0,LT_M,LT_M0,
66 3 IADK,JDIK,IADK0,JDIK0,IADM,
67 4 JDIM,IADM0,JDIM0,NNZM,TOLS,
68 5 NLIM,itol,EPS_M,IPREC)
69C-----------------------------------------------
70C I m p l i c i t T y p e s
71C-----------------------------------------------
72#include "implicit_f.inc"
73C-----------------------------------------------
74C D u m m y A r g u m e n t s
75C-----------------------------------------------
76 my_real
77 . x(*), r(*), tols, eps_m,
78 . diag_k(*), diag_m(*),
79 . lt_k(*), lt_k0(*), lt_m(*), lt_m0(*)
80 INTEGER NDDL, NNZ, NNZM, NLIM, ITOL, IPREC,
81 . iadk(*), jdik(*), iadk0(*), jdik0(*),
82 . iadm(*), jdim(*), iadm0(*), jdim0(*)
83C-----------------------------------------------
84C L o c a l V a r i a b l e s
85C-----------------------------------------------
86 OPEN(66,file="input_mat",form="unformatted")
87 WRITE(66)nddl,nnz,nnzm,nlim,itol,iprec
88 WRITE(66)x(1:nddl)
89 WRITE(66)r(1:nddl)
90 WRITE(66)diag_k(1:nddl)
91 WRITE(66)diag_m(1:nddl)
92 WRITE(66)lt_k(1:nnz)
93 WRITE(66)lt_k0(1:nnz)
94 WRITE(66)lt_m(1:nnzm)
95 WRITE(66)lt_m0(1:nnzm)
96 WRITE(66)iadk(1:nddl+1)
97 WRITE(66)jdik(1:iadk(nddl+1)-1)
98 WRITE(66)iadk0(1:nddl+1)
99 WRITE(66)jdik0(1:iadk0(nddl+1)-1)
100 WRITE(66)iadm(1:nddl+1)
101 WRITE(66)jdim(1:iadm(nddl+1)-1)
102 WRITE(66)iadm0(1:nddl+1)
103 WRITE(66)jdim0(1:iadm0(nddl+1)-1)
104 WRITE(66)eps_m,tols
105 RETURN
106 END
107
108!||====================================================================
109!|| imp_isave2 ../engine/source/implicit/imp_pcg.F
110!||--- called by ------------------------------------------------------
111!|| imp_pcgh ../engine/source/implicit/imp_pcg.F
112!||====================================================================
113 SUBROUTINE imp_isave2(
114 1 NDDL,X,DIAG_K,NNZ,LT_K,
115 3 IADK,JDIK,LT_K0,IADK0,JDIK0)
116C-----------------------------------------------
117C I m p l i c i t T y p e s
118C-----------------------------------------------
119#include "implicit_f.inc"
120C-----------------------------------------------
121C D u m m y A r g u m e n t s
122C-----------------------------------------------
123 my_real
124 . x(*),
125 . diag_k(*),
126 . lt_k(*),lt_k0(*)
127 INTEGER NDDL, NNZ,
128 . IADK(*), JDIK(*),IADK0(*), JDIK0(*)
129C-----------------------------------------------
130C L o c a l V a r i a b l e s
131C-----------------------------------------------
132 INTEGER I, J, K
133
134 OPEN(66,file="input_matA")
135 WRITE(66,*)nddl,nddl,nnz+nddl
136 DO i = 1, nddl
137 DO k = iadk(i), iadk(i+1)-1
138 j=jdik(k)
139 IF(j < i) THEN
140 WRITE(66,*)i,j,lt_k(k)
141 else
142 print*,'error1'
143 END IF
144 END DO
145 WRITE(66,*)i,i,diag_k(i)
146 DO k = iadk0(i), iadk0(i+1)-1
147 j=jdik0(k)
148 IF(j > i) THEN
149 WRITE(66,*)i,j,lt_k0(k)
150 else
151 print*,'error2'
152 END IF
153 END DO
154 END DO
155 OPEN(67,file="input_vecX")
156 WRITE(67,*)nddl
157 DO i = 1, nddl
158 WRITE(67,*)x(i)
159 END DO
160 RETURN
161 END
162
163!||====================================================================
164!|| imp_rsave ../engine/source/implicit/imp_pcg.F
165!||--- called by ------------------------------------------------------
166!|| imp_pcgh ../engine/source/implicit/imp_pcg.F
167!||====================================================================
168 SUBROUTINE imp_rsave(
169 1 NDDL,X,R)
170C-----------------------------------------------
171C I m p l i c i t T y p e s
172C-----------------------------------------------
173#include "implicit_f.inc"
174C-----------------------------------------------
175C D u m m y A r g u m e n t s
176C-----------------------------------------------
177 my_real
178 . x(*), r(*)
179 INTEGER NDDL
180C-----------------------------------------------
181C L o c a l V a r i a b l e s
182C-----------------------------------------------
183
184 WRITE(66)X(1:NDDL)
185 WRITE(66)r(1:nddl)
186 CLOSE(66)
187 RETURN
188 END
189
190
191
192!||====================================================================
193!|| imp_pcgh ../engine/source/implicit/imp_pcg.F
194!||--- called by ------------------------------------------------------
195!|| lin_solvh0 ../engine/source/implicit/lin_solv.F
196!|| lin_solvh1 ../engine/source/implicit/lin_solv.F
197!|| lin_solvhm ../engine/source/implicit/lin_solv.F
198!|| lin_solvih2 ../engine/source/implicit/lin_solv.F
199!||--- calls -----------------------------------------------------
200!|| ancmsg ../engine/source/output/message/message.F
201!|| arret ../engine/source/system/arret.F
202!|| crit_stop ../engine/source/implicit/imp_pcg.F
203!|| imp_isave ../engine/source/implicit/imp_pcg.F
204!|| imp_isave2 ../engine/source/implicit/imp_pcg.F
205!|| imp_rsave ../engine/source/implicit/imp_pcg.F
206!|| mav_lth ../engine/source/implicit/produt_v.f
207!|| mixedsol ../engine/source/implicit/imp_pcg.F
208!|| my_barrier ../engine/source/system/machine.F
209!|| nlim_mix ../engine/source/implicit/imp_pcg.F
210!|| omp_get_num_threads ../engine/source/engine/openmp_stub.F90
211!|| prec_solvh ../engine/source/implicit/prec_solv.F
212!|| produt_h ../engine/source/implicit/produt_v.F
213!|| pupd ../engine/source/implicit/imp_pcg.F
214!|| spmd_sum_s ../engine/source/mpi/implicit/imp_spmd.f
215!|| spmd_sumfc_v ../engine/source/mpi/implicit/imp_spmd.F
216!|| startime ../engine/source/system/timer_mod.f90
217!|| stoptime ../engine/source/system/timer_mod.F90
218!||--- uses -----------------------------------------------------
219!|| dsgraph_mod ../engine/share/modules/dsgraph_mod.F
220!|| groupdef_mod ../common_source/modules/groupdef_mod.F
221!|| imp_frk ../engine/share/modules/impbufdef_mod.F
222!|| imp_intm ../engine/share/modules/imp_intm.F
223!|| imp_workh ../engine/share/modules/impbufdef_mod.F
224!|| intbufdef_mod ../common_source/modules/interfaces/intbufdef_mod.F90
225!|| message_mod ../engine/share/message_module/message_mod.F
226!||====================================================================
227 SUBROUTINE imp_pcgh( IPREC ,
228 1 NDDL ,NNZ ,IADK ,JDIK ,DIAG_K ,
229 2 LT_K ,NDDLI ,ITOK ,IADI ,JDII ,
230 3 LT_I ,NNZM ,IADM ,JDIM ,DIAG_M ,
231 4 LT_M ,X ,R ,ITOL ,TOL ,
232 5 P ,Z ,Y ,ITASK ,IPRINT ,
233 6 N_MAX ,EPS_M ,F_X ,ISTOP ,W_DDL ,
234 8 A ,AR ,VE ,MS ,XE ,
235 9 D ,DR ,NDOF ,IPARI ,INTBUF_TAB,
236 A NUM_IMP,NS_IMP,NE_IMP,NSREM ,
237 B NSL ,NMONV ,IMONV ,MONVOL,IGRSURF ,
238 C VOLMON,FR_MV ,IBFV ,SKEW ,
239 D XFRAME ,GRAPHE,IAD_ELEM,FR_ELEM,ITAB ,
240 E INSOLV ,ITN ,FAC_K ,IPIV_K,NK ,
241 F MUMPS_PAR,CDDLP,ISOLV,IDSC ,IDDL ,
242 G IKC ,INLOC ,IND_IMP,XI_C ,R0 ,
243 H NDDLI_G,INTP_C,IRBE3 ,LRBE3 ,IRBE2 ,
244 I LRBE2 )
245C-----------------------------------------------
246C M o d u l e s
247C-----------------------------------------------
248 USE dsgraph_mod
249 USE imp_workh
250 USE imp_frk
251 USE imp_intm
252 USE message_mod
253 USE intbufdef_mod
254 USE groupdef_mod
255C-----------------------------------------------
256C I m p l i c i t T y p e s
257C-----------------------------------------------
258#include "implicit_f.inc"
259C-----------------------------------------------
260C C o m m o n B l o c k s
261C-----------------------------------------------
262#include "comlock.inc"
263#include "com04_c.inc"
264#include "units_c.inc"
265#include "task_c.inc"
266#if defined(MUMPS5)
267#include "dmumps_struc.h"
268#endif
269C-----------------------------------------------
270C D u m m y A r g u m e n t s
271C-----------------------------------------------
272C----------resol [K]{X}={F}---------
273 INTEGER NDDL ,NNZ ,ITOL,
274 . nddli ,itok(*) ,iadi(*),jdii(*),n_max,
275 . nnzm ,iprec,itask,iprint,
276 . istop,w_ddl(*),ibfv(*),intp_c,irbe3(*) ,lrbe3(*),
277 . irbe2(*) ,lrbe2(*)
278 INTEGER NDOF(*),NE_IMP(*),NSREM ,NSL,
279 . ipari(*) ,num_imp(*),ns_imp(*) ,ind_imp(*)
280 INTEGER NMONV,IMONV(*),MONVOL(*),FR_MV(*)
281 INTEGER IAD_ELEM(2,*), FR_ELEM(*), ITAB(*),
282 . INSOLV, ITN, IPIV_K(*), NK, CDDLP(*), ISOLV, IDSC,
283 . IDDL(*), IKC(*), INLOC(*),NDDLI_G, MICID
284 my_real
285 . LT_I(*), TOL, EPS_M, F_X, XI_C(*)
286#if defined knf
287 INTEGER , target, intent(inout) ::
288 . IADK(NDDL+1) ,JDIK(NNZ), IADM(NDDL+1), JDIM(NNZM)
289 my_real , target, intent(inout) ::
290 . DIAG_K(NDDL), LT_K(NNZ) ,DIAG_M(NDDL),LT_M(NNZM) ,X(NDDL),
291 . P(NDDL) ,Z(NDDL) ,R(NDDL) ,Y(NDDL)
292#elif 1
293 INTEGER
294 . IADK(*) ,JDIK(*), IADM(*),JDIM(*)
295 my_real
296 . DIAG_K(*), LT_K(*) ,DIAG_M(*),LT_M(*) ,X(*),
297 . P(*) ,Z(*) ,R(*) ,Y(*)
298#endif
299 my_real
300 . a(3,*),ar(3,*),ve(3,*),d(3,*),dr(3,*),xe(3,*),
301 . ms(*),volmon(*),skew(*) ,xframe(*),fac_k(*),
302 . r0(*)
303 TYPE(prgraph) :: GRAPHE(*)
304
305#ifdef MUMPS5
306 TYPE(dmumps_struc) MUMPS_PAR
307#else
308 ! Fake declaration as DMUMPS_STRUC is shipped with MUMPS
309 INTEGER MUMPS_PAR
310#endif
311
312 TYPE(intbuf_struct_) INTBUF_TAB(*)
313 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
314C-----------------------------------------------
315C External function
316C-----------------------------------------------
317 LOGICAL MIXEDSOL
318 EXTERNAL mixedsol
319#ifdef MUMPS5
320C-----------------------------------------------
321C L o c a l V a r i a b l e s
322C-----------------------------------------------
323 INTEGER I,J,IT,IP,NLIM,ND,IERR,IPRI,IERROR,NNZI,LOC_PROC,
324 . ick,isave,n, tb
325 parameter(nd=4)
326 CHARACTER*4 EXIT
327 CHARACTER*11 WARR
328 my_real
329 . s , r2, r02,alpha,beta,g0,g1,rr,tols,toln,tols2
330 INTEGER CRIT_STOP,IUPD,IUPD0,F_DDL,L_DDL,IFLG,GPUR4R8,
331 . itp,lcom,nddl1,lcomi,lcomx,k,
332 . ntag(nddl),igather(nddl)
333 EXTERNAL crit_stop
334
335 my_real
336 . anorm2,xnorm2,l_a,l_b2,l_b,a_old,b_old,tmp,r2_old,rmax
337 my_real
338 . cs,dbar, delta, denom, kcond,snprod,qrnorm,
339 . gamma, gbar, gmax, gmin, epsln,lqnorm,diag,cgnorm,
340 . oldb, rhs1, rhs2,sn, zbar, zl ,oldb2,tnorm2,eps(4)
341 double precision
342 . tt, tf
343C--------------GPU SPECIFIC--------------------------
344 real*4, DIMENSION(:),ALLOCATABLE ::
345 . lt_k_sp, lt_k0_sp, lt_m_sp, lt_m0_sp, diag_k_sp, diag_m_sp,
346 . x_sp, r_sp,w_sp, lt_i0_sp
347#if defined knf
348 my_real, DIMENSION(:),ALLOCATABLE :: wr
349 INTEGER, DIMENSION(:),ALLOCATABLE :: IFRTMP, TABLE, INDTMP
350#elif 1
351 my_real, DIMENSION(:),ALLOCATABLE :: w, wr, vgat, vsca, xir, yir
352 INTEGER, DIMENSION(:),ALLOCATABLE :: IFRTMP, INDEX, TABLE, INDTMP
353#endif
354
355c#if defined knf
356c!dec$ attributes offload: mic ::
357c & M_IADK, M_IADM, M_JDIM, M_JDIK
358c INTEGER, pointer, save, dimension(:) ::
359c . M_IADK, M_IADM, M_JDIM, M_JDIK
360c!dec$ attributes offload: mic ::
361c & M_LT_K, M_LT_M, M_DIAG_K, M_DIAG_M, M_P,
362c & M_X, M_R, M_Y, M_Z
363c DOUBLE PRECISION, pointer, save, dimension(:) ::
364c . M_LT_K, M_LT_M, M_DIAG_K, M_DIAG_M, M_P,
365c . M_X, M_R, M_Y, M_Z
366c#endif
367C--------------END GPU SPECIFIC----------------------
368 DATA EXIT /'WITH'/
369 DATA warr /'**WARRING**'/
370C--------------INITIALISATION--------------------------
371 isave=0
372C ISAVE=1 ! sauvegarde input/res
373C ISAVE=2 ! sauvegarde input simple
374#if defined knf
375 IF(nddli > 0)THEN
376 nnzi = iadi0(nddli+1)-iadi0(1)
377 ELSE
378 nnzi = 0
379 END IF
380C
381 lcom=0
382 lcomi=0
383 lcomx=0
384 IF(nspmd > 1)THEN
385 DO i = 1, nspmd
386 lcom = lcom + nd_fr(i)
387 END DO
388 IF(nddli_g > 0)THEN
389 ntag(1:nddl)=0
390 DO i = 1, nspmd
391 DO j = iad_sl(i),iad_sl(i+1)-1
392 k=iddl_sl(j)
393 IF(ntag(k)==0)THEN
394 lcomi=lcomi+1
395 ntag(k)=lcomi
396 igather(lcomi)=k-1
397 ENDIF
398 ENDDO
399 END DO
400 DO i = 1, nspmd
401 lcomx=lcomx+f_ddl+iad_srem(i+1)-iad_srem(i)+1
402 END DO
403 DO i=1,nddl_si
404 DO j =iad_si(i),iad_si(i+1)-1
405 k =jdi_si(j)
406 IF(ntag(k)==0)THEN
407 lcomi=lcomi+1
408 ntag(k)=lcomi
409 igather(lcomi)=k-1
410 ENDIF
411 ENDDO
412 ENDDO
413C
414 DO i=1,nddl_sl
415 k = iddl_sl(i)
416 IF(ntag(k)==0)THEN
417 lcomi=lcomi+1
418 ntag(k)=lcomi
419 igather(lcomi)=k-1
420 ENDIF
421C
422 DO j =iad_ss(i),iad_ss(i+1)-1
423 k =jdi_sl(j)
424 IF(ntag(k)==0)THEN
425 lcomi=lcomi+1
426 ntag(k)=lcomi
427 igather(lcomi)=k-1
428 ENDIF
429 ENDDO
430 ENDDO
431 END IF
432
433C adding stop test in case of SPMD and null domain => bad because need to communicate V
434C IF (NDDL==0 .AND. NNZ==0 .AND. LCOM==0 .AND. LCOMI==0
435C . .AND. LCOMX==0) RETURN
436 END IF
437#endif
438 loc_proc=ispmd+1
439C-----Istop=1 : NIT>=N_LIM+ RR>TOL
440C------X(I)=ZERO----******N_MAX,EPS_M a calculer avant----
441 f_ddl=1+itask*nddl/nthread
442 l_ddl=(itask+1)*nddl/nthread
443C
444 istop=0
445 nlim=n_max
446 CALL nlim_mix(n_max,nddli_g,nlim)
447 nlim=max(nlim,2)
448 tols=max(tol,eps_m)
449 ipri = iprint
450 IF (ispmd/=0) ipri = 0
451C--------
452C routine speciale pour sauvegarder les input
453 IF(isave==1)CALL imp_isave(
454 1 nddl,x,r,diag_k,diag_m,
455 2 nnz,lt_k,lt_k0,lt_m,lt_m0,
456 3 iadk,jdik,iadk0,jdik0,iadm,
457 4 jdim,iadm0,jdim0,nnzm,tols,
458 5 nlim,itol,eps_m,iprec)
459 IF(isave==2)CALL imp_isave2(
460 1 nddl,r,diag_k,nnz,lt_k,
461 3 iadk,jdik,lt_k0,iadk0,jdik0)
462C--------
463 IF (itask == 0) THEN
464 IF(ipri/=0)THEN
465 WRITE(iout,*)' *** begin conjugate gradient iteration ***'
466 WRITE(IOUT,*)
467 ENDIF
468 IF(IPRI<0)THEN
469 WRITE(ISTDO,*)' *** begin conjugate gradient iteration ***'
470 WRITE(ISTDO,*)
471 ENDIF
472 END IF !(ITASK == 0) THEN
473 IP = IABS(IPRI)
474 IT=0
475C
476 IUPD = 0
477 IUPD0 = 0
478C----------------------
479 CALL MY_BARRIER
480C---------------------
481C 1 : no gpu
482C -1 : nvidia dp
483C -3 : ocl dp
484cc NOGPU=-1
485cc IF(NOGPU==1)THEN
486C variable to switch between double (GPUR4R8=2)or single (GPUR4R8=1)precision GPU
487 GPUR4R8=2
488
489#if !defined knf
490
491 CALL PREC_SOLVH(IPREC, ITASK ,
492 1 GRAPHE,IAD_ELEM,FR_ELEM,DIAG_K,LT_K ,
493 2 IADK ,JDIK ,ITAB ,IPRI,INSOLV ,
494 3 ITN ,FAC_K , IPIV_K, NK ,MUMPS_PAR,
495 4 CDDLP ,ISOLV , IDSC , IDDL ,IKC ,
496 5 INLOC ,NDOF , NDDL ,NNZM ,IADM ,
497 6 JDIM ,DIAG_M , LT_M ,R ,Z ,
498 7 F_DDL ,L_DDL )
499C----------------------
500 CALL MY_BARRIER
501C---------------------
502 DO I = F_DDL,L_DDL
503 P(I) = Z(I)
504 ENDDO
505 CALL PRODUT_H(F_DDL ,L_DDL ,R ,Z ,W_DDL , G0 ,ITASK )
506C---------------------
507 CALL MAV_LTH(
508 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K,
509 2 LT_K ,IADI ,JDII ,ITOK ,LT_I ,
510 3 P ,Y ,A ,AR ,
511 5 VE ,MS ,XE ,D ,DR ,
512 6 NDOF ,IPARI ,INTBUF_TAB ,NUM_IMP,
513 7 NS_IMP,NE_IMP,NSREM ,NSL ,IBFV ,
514 8 SKEW ,XFRAME,MONVOL,VOLMON,IGRSURF ,
515 9 FR_MV,NMONV ,IMONV ,IND_IMP,
516 A XI_C ,IUPD ,IRBE3 ,LRBE3 ,IRBE2 ,
517 B LRBE2 ,F_DDL ,L_DDL ,ITASK )
518C----------------------
519 CALL MY_BARRIER
520C---------------------
521 CALL PRODUT_H(F_DDL ,L_DDL ,P ,Y ,W_DDL , S ,ITASK )
522C---------------------
523 IF (G0==ZERO) THEN
524 IF (ITOL>1) THEN
525 ISTOP=-1
526 GOTO 200
527 ELSE
528 ALPHA = ZERO
529 ENDIF
530 ELSE
531 ALPHA = G0/S
532 ENDIF
533 TOLS2=TOLS*TOLS
534 IF (ITOL==1) THEN
535 CALL PRODUT_H( F_DDL ,L_DDL ,R ,R ,W_DDL , R02 ,ITASK )
536C---------------------
537 R2 =R02
538 ELSEIF (ITOL==3) THEN
539C------ R2<TOL*TOL*ANORM2*XNORM2------
540 R02=ABS(G0)
541 R2 =R02
542 L_A=ONE/ALPHA
543C L_B2=R02
544 TNORM2=L_A*L_A
545 ANORM2=ZERO
546 XNORM2=ZERO
547 A_OLD=L_A
548 B_OLD=ZERO
549 OLDB = SQRT(R02)
550 ELSEIF (ITOL==4) THEN
551 R02=ALPHA*ALPHA*ABS(G0)
552 EPS(1)=R02
553 R2=R02
554 ELSE
555 R02=ABS(G0)
556 R2 =R02
557 ENDIF
558 IF (R02==ZERO) THEN
559 ISTOP=-1
560 GOTO 200
561 ENDIF
562 TOLN=R02*TOLS2
563.AND. IF(IPRI/=0ITASK==0)THEN
564 RR = SQRT(R2)
565 IF(IPRI<0) WRITE(ISTDO,1000)IT,RR
566 WRITE(IOUT,1000)IT,RR
567 ENDIF
568C-------pour etre coherent avec lanzos for linear
569C----------------------
570 CALL MY_BARRIER
571C---------------------
572 IT=1
573 DO I=F_DDL,L_DDL
574 X(I) = X(I) + ALPHA*P(I)
575 R(I) = R(I) - ALPHA*Y(I)
576 ENDDO
577C----------------------
578 CALL MY_BARRIER
579C---------------------
580 CALL PREC_SOLVH(IPREC, ITASK ,
581 1 GRAPHE,IAD_ELEM,FR_ELEM,DIAG_K,LT_K ,
582 2 IADK ,JDIK ,ITAB ,IPRI,INSOLV ,
583 3 ITN ,FAC_K , IPIV_K, NK ,MUMPS_PAR,
584 4 CDDLP ,ISOLV , IDSC , IDDL ,IKC ,
585 5 INLOC ,NDOF , NDDL ,NNZM ,IADM ,
586 6 JDIM ,DIAG_M , LT_M ,R ,Z ,
587 7 F_DDL ,L_DDL )
588C----------------------
589 CALL MY_BARRIER
590C---------------------
591 CALL PRODUT_H( F_DDL ,L_DDL ,R ,Z ,W_DDL , G1 ,ITASK )
592C---------------------
593 BETA=G1/G0
594 IF (ITOL==1) THEN
595 CALL PRODUT_H( F_DDL ,L_DDL ,R ,R ,W_DDL , R2 ,ITASK )
596 ELSEIF (ITOL==3) THEN
597C------ R2<TOL*TOL*ANORM2*XNORM2------
598 R2=ABS(G1)
599 L_B2=ABS(BETA)*A_OLD*A_OLD
600 L_B=SQRT(L_B2)
601 TNORM2=TNORM2+L_B2
602 B_OLD=BETA
603C* INITIALIZE OTHER QUANTITIES.
604 GBAR = L_A
605 DBAR = L_B
606 RHS1 = OLDB
607 RHS2 = ZERO
608 GMAX = ABS( L_A ) + EPS_M
609 GMIN = GMAX
610 OLDB2 = L_B2
611 OLDB = L_B
612 ELSEIF (ITOL==4) THEN
613 R2=R02
614 ELSE
615 R2=ABS(G1)
616 ENDIF
617 G0 = G1
618 IF (ITOL==3) TOLN=TOLN*TNORM2
619C----------------------
620 CALL MY_BARRIER
621C---------------------
622#include "lockon.inc"
623 ISTOP=CRIT_STOP(IT,R2,NLIM,TOLN)
624#include "lockoff.inc"
625.AND. IF(NDDLI_G>0INTP_C==-1)IUPD0 = -INTP_C
626 DO WHILE (ISTOP==1)
627 DO I=F_DDL,L_DDL
628 P(I) = Z(I) + BETA*P(I)
629 ENDDO
630.AND. IF(IUPD0>0IT>NDDLI_G/2) IUPD = IUPD0
631C----------------------
632 CALL MY_BARRIER
633C---------------------
634 CALL MAV_LTH(
635 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K,
636 2 LT_K ,IADI ,JDII ,ITOK ,LT_I ,
637 3 P ,Y ,A ,AR ,
638 5 VE ,MS ,XE ,D ,DR ,
639 6 NDOF ,IPARI ,INTBUF_TAB ,NUM_IMP,
640 7 NS_IMP,NE_IMP,NSREM ,NSL ,IBFV ,
641 8 SKEW ,XFRAME,MONVOL,VOLMON,IGRSURF ,
642 9 FR_MV,NMONV ,IMONV ,IND_IMP,
643 A XI_C ,IUPD ,IRBE3 ,LRBE3 ,IRBE2 ,
644 B LRBE2 ,F_DDL ,L_DDL ,ITASK )
645C----------------------
646 CALL MY_BARRIER
647C---------------------
648 CALL PRODUT_H(F_DDL ,L_DDL ,P ,Y ,W_DDL , S ,ITASK )
649C---------------------
650 ALPHA=G0/S
651 DO I=F_DDL,L_DDL
652 X(I) = X(I) + ALPHA*P(I)
653 R(I) = R(I) - ALPHA*Y(I)
654 ENDDO
655C----------------------
656 CALL MY_BARRIER
657C---------------------
658 CALL PREC_SOLVH(IPREC, ITASK ,
659 1 GRAPHE,IAD_ELEM,FR_ELEM,DIAG_K,LT_K ,
660 2 IADK ,JDIK ,ITAB ,IPRI,INSOLV ,
661 3 ITN ,FAC_K , IPIV_K, NK ,MUMPS_PAR,
662 4 CDDLP ,ISOLV , IDSC , IDDL ,IKC ,
663 5 INLOC ,NDOF , NDDL ,NNZM ,IADM ,
664 6 JDIM ,DIAG_M , LT_M ,R ,Z ,
665 7 F_DDL ,L_DDL )
666C----------------------
667 CALL MY_BARRIER
668C---------------------
669 CALL PRODUT_H(F_DDL ,L_DDL ,R ,Z ,W_DDL , G1 ,ITASK )
670C---------------------
671 BETA=G1/G0
672 G0 = G1
673 IF (ITOL==1) THEN
674 CALL PRODUT_H( F_DDL ,L_DDL ,R ,R ,W_DDL , R2 ,ITASK )
675 ELSEIF (ITOL==3) THEN
676 R2 =ABS(G1)
677 S=ONE/ALPHA
678 L_A=S+B_OLD*A_OLD
679C----------L_B2 : beta(j-1)^2-A_OLD :1/alpha(j-1)-------
680 L_B2=ABS(BETA)*S*S
681 L_B=SQRT(L_B2)
682 A_OLD=S
683 B_OLD=BETA
684 ANORM2=TNORM2
685 TNORM2=TNORM2+L_A*L_A+OLDB2+L_B2
686 GAMMA = SQRT( GBAR*GBAR + OLDB2 )
687 CS = GBAR / GAMMA
688 SN = OLDB / GAMMA
689 DELTA = CS * DBAR + SN * L_A
690 GBAR = SN * DBAR - CS * L_A
691 EPSLN = SN * L_B
692 DBAR = - CS * L_B
693 ZL = RHS1 / GAMMA
694 XNORM2 = XNORM2+ZL*ZL
695 GMAX = MAX( GMAX, GAMMA )
696 GMIN = MIN( GMIN, GAMMA )
697 RHS1 = RHS2 - DELTA * ZL
698 RHS2 = - EPSLN * ZL
699 TOLN=TOLS2*ANORM2*XNORM2
700 OLDB2 = L_B2
701 OLDB = L_B
702 ELSEIF (ITOL==4) THEN
703 TMP=ALPHA*ALPHA*ABS(G1)
704 IF (IT>=ND) THEN
705 DO J=1,ND-1
706 EPS(J)=EPS(J+1)
707 ENDDO
708 EPS(ND)=TMP
709 R2=ZERO
710 DO J=1,ND
711 R2=R2+EPS(J)
712 ENDDO
713 ELSE
714 EPS(IT+1)=TMP
715 R2=R2+TMP
716 ENDIF
717 ELSE
718 R2 =ABS(G1)
719 ENDIF
720C
721.AND. IF(IPRI/=0ITASK==0)THEN
722 IF(MOD(IT,IP)==0)THEN
723 RR = SQRT(R2/R02)
724 WRITE(IOUT,1001)IT,RR
725 IF(IPRI<0) WRITE(ISTDO,1001)IT,RR
726 ENDIF
727 ENDIF
728C----------------------
729 CALL MY_BARRIER
730C---------------------
731#include "lockon.inc"
732 ISTOP=CRIT_STOP(IT,R2,NLIM,TOLN)
733#include "lockoff.inc"
734
735.AND..OR. IF((IUPD>0IT==NLIM)
736.AND..AND. . (IUPD==0ISTOP/=1IUPD0>0))THEN
737C------- re-iter with linear Kint------
738 IUPD0 = 0
739 IF(IUPD>0) IUPD = 0
740#include "lockon.inc"
741 ISTOP = 1
742#include "lockoff.inc"
743
744 NLIM = NLIM + NLIM
745 DO I=F_DDL,L_DDL
746 X(I) = ZERO
747 ENDDO
748C----------------------
749 CALL MY_BARRIER
750C---------------------
751 CALL MAV_LTH(
752 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K,
753 2 LT_K ,IADI ,JDII ,ITOK ,LT_I ,
754 3 X ,Z ,A ,AR ,
755 5 VE ,MS ,XE ,D ,DR ,
756 6 NDOF ,IPARI ,INTBUF_TAB ,NUM_IMP,
757 7 NS_IMP,NE_IMP,NSREM ,NSL ,IBFV ,
758 8 SKEW ,XFRAME,MONVOL,VOLMON,IGRSURF ,
759 9 FR_MV,NMONV ,IMONV ,IND_IMP,
760 A XI_C ,IUPD ,IRBE3 ,LRBE3 ,IRBE2 ,
761 B LRBE2 ,F_DDL ,L_DDL ,ITASK )
762C----------------------
763 CALL MY_BARRIER
764C---------------------
765 DO I=F_DDL,L_DDL
766 R(I) = R0(I)-Z(I)
767 ENDDO
768C----------------------
769 CALL MY_BARRIER
770C---------------------
771 CALL PREC_SOLVH(IPREC, ITASK ,
772 1 GRAPHE,IAD_ELEM,FR_ELEM,DIAG_K,LT_K ,
773 2 IADK ,JDIK ,ITAB ,IPRI,INSOLV ,
774 3 ITN ,FAC_K , IPIV_K, NK ,MUMPS_PAR,
775 4 CDDLP ,ISOLV , IDSC , IDDL ,IKC ,
776 5 INLOC ,NDOF , NDDL ,NNZM ,IADM ,
777 6 JDIM ,DIAG_M , LT_M ,R ,Z ,
778 7 F_DDL ,L_DDL )
779C----------------------
780 CALL MY_BARRIER
781C---------------------
782 CALL PRODUT_H( F_DDL ,L_DDL ,R ,Z ,W_DDL , G0 ,ITASK )
783 BETA = ZERO
784C------ R2<TOL*TOL*ANORM2*XNORM2------
785 IF (ITOL==3) THEN
786 TNORM2=L_A*L_A
787 ANORM2=ZERO
788 XNORM2=ZERO
789 A_OLD=ZERO
790 B_OLD=ZERO
791 L_B=ZERO
792C* INITIALIZE OTHER QUANTITIES.
793 GBAR = L_A
794 DBAR = ZERO
795 RHS1 = SQRT(R02)
796 RHS2 = ZERO
797 GMAX = ABS( L_A ) + EPS_M
798 GMIN = GMAX
799 OLDB2 = ZERO
800 OLDB = L_B
801 ELSEIF (ITOL==4) THEN
802 DO J=1,ND
803 EPS(J)=R02
804 ENDDO
805 ENDIF
806 ENDIF
807C----------------------
808 CALL MY_BARRIER
809C---------------------
810 IT = IT +1
811 ENDDO
812 200 CONTINUE
813
814#elif defined knf
815 IF(GPUR4R8==1)THEN
816 ELSEIF(GPUR4R8==2)THEN
817 IF(ITASK==0)THEN
818C Code MIC simple precision
819
820C debut code specifique MIC double precision
821CAD
822CAD
823CAD DIRECTIVES POUR MIC
824CAD
825C init MIC
826 CALL MIC_INIT(
827 . NSPMD,LOC_PROC,L_SPMD(LOC_PROC),IDEVICE,MICID,IERR)
828 IF(IERR < 0) THEN
829 IF(IERR==-1)THEN
830 CALL ANCMSG(MSGID=223,ANMODE=ANINFO_BLIND)
831 ELSEIF(IERR==-2)THEN
832 CALL ANCMSG(MSGID=224,ANMODE=ANINFO_BLIND)
833 ELSEIF(IERR==-3)THEN
834 CALL ANCMSG(MSGID=225,ANMODE=ANINFO_BLIND)
835 END IF
836 CALL ARRET(2)
837 END IF
838 NDDL1=NDDL+1
839 print *,' init mic card number ',MICID
840C dummy variables to be replaced by pointer later
841c M_LT_K => LT_K(1:NNZ)
842c M_JDIK => JDIK(1:NNZ)
843c M_LT_M => LT_M(1:NNZM)
844c M_JDIM => JDIM(1:NNZM)
845c M_DIAG_K => DIAG_K(1:NDDL)
846c M_DIAG_M => DIAG_M(1:NDDL)
847c M_IADK => IADK(1:NDDL1)
848c M_IADM => IADM(1:NDDL1)
849c M_P => P(1:NDDL)
850c M_X => X(1:NDDL)
851c M_R => R(1:NDDL)
852c M_Y => Y(1:NDDL)
853c M_Z => Z(1:NDDL)
854cc call omp_set_num_threads_target (TARGET_MIC, MICID,116)
855 i= omp_get_num_threads()
856 print *,'number of threads(host):',i
857cc call kmp_set_blocktime_target (TARGET_MIC, MICID,INFINITE)
858 i= kmp_get_blocktime_target (TARGET_MIC, MICID)
859!dec$ attributes offload: mic :: ONE,ZERO
860!dec$ attributes offload: mic :: MIC_DSCAL,MIC_GATHER,MIC_SCATTER
861!dec$ attributes offload: mic :: MIC_DDOT,MIC_MV,MIC_SPMV,MIC_DAXPY
862!dec$ attributes offload: mic :: MIC_DCOPY
863 TF=OMP_GET_WTIME()
864 TB=0
865!dir$ omp offload target(mic:MICID)
866 & in (LT_K:length(nnz), free_if(.false.), align(512))
867!$OMP PARALLEL
868!$OMP END PARALLEL
869 TB = TB+NNZ
870 print *,' variables1= ',nnz
871!dir$ omp offload target(mic:MICID)
872 & in (JDIK:length(nnz), free_if(.false.), align(512))
873!$OMP PARALLEL
874!$OMP END PARALLEL
875 TB = TB+NNZ
876 print *,' variables2= ',nnz
877!dir$ omp offload target(mic:MICID)
878 & in (LT_K0:length(nnz), free_if(.false.), align(512))
879!$OMP PARALLEL
880!$OMP END PARALLEL
881 TB = TB+NNZ
882 print *,' variables3= ',nnz
883!dir$ omp offload target(mic:MICID)
884 & in (JDIK0:length(nnz), free_if(.false.), align(512))
885!$OMP PARALLEL
886!$OMP END PARALLEL
887 TB = TB+NNZ
888 print *,' variables4= ',nnz
889 IF(IPREC == 5) THEN
890!dir$ omp offload target(mic:MICID)
891& in (lt_m:length(nnzm), free_if(.false.), align(512))
892!$OMP PARALLEL
893!$OMP END PARALLEL
894 TB = TB+NNZM
895 print *,' variables5= ',nnzm
896!dir$ omp offload target(mic:MICID)
897 & in (JDIM:length(nnzm), free_if(.false.), align(512))
898!$OMP PARALLEL
899!$OMP END PARALLEL
900 TB = TB+NNZM
901 print *,' variables6= ',nnzm
902!dir$ omp offload target(mic:MICID)
903 & in (lt_m0:length(nnzm), free_if(.false.), align(512))
904!$OMP PARALLEL
905!$OMP END PARALLEL
906 TB = TB+NNZM
907 print *,' variables7= ',nnzm
908!dir$ omp offload target(mic:MICID)
909 & in (JDIM0:length(nnzm), free_if(.false.), align(512))
910!$OMP PARALLEL
911!$OMP END PARALLEL
912 TB = TB+NNZ
913 print *,' variables8= ',nnzm
914 END IF
915!dir$ omp offload target(mic:MICID)
916 & in (DIAG_K,DIAG_M,x,y,z,r,p:length(nddl),
917 & free_if(.false.))
918!$OMP PARALLEL
919!$OMP END PARALLEL
920 TB = TB+NDDL*7
921 print *,' variables9= ',nddl*7
922!dir$ omp offload target(mic:MICID)
923 & in (IADK,IADK0,IADM,IADM0:length(NDDL1),
924 & free_if(.false.), align(512))
925!$OMP PARALLEL
926!$OMP END PARALLEL
927 TB = TB+nddl1*4
928 print *,' variables10= ',nddl1*4
929 TF=OMP_GET_WTIME()-TF
930 print *,'transfer time cpu => mic(s) =',TF
931 print *,'transfer rate cpu => mic(mb/s)=',
932 . ((TB/(1024*1024))*8)/TF
933 IF(NSPMD > 1)THEN
934 ALLOCATE(W(NDDL),STAT=IERROR)
935 IERROR=IERROR+IERR
936 ALLOCATE(VGAT(LCOM),STAT=IERR)
937 IERROR=IERROR+IERR
938 ALLOCATE(VSCA(LCOM),STAT=IERR)
939 IERROR=IERROR+IERR
940 ALLOCATE(INDEX(LCOM),STAT=IERR)
941 IERROR=IERROR+IERR
942 ALLOCATE(TABLE(NDDL),STAT=IERR)
943 IERROR=IERROR+IERR
944 IF(IERROR /= 0)THEN
945C code erreur
946 END IF
947
948 DO I = 1, NDDL
949C transfo entier vers my_real
950 W(I) = W_DDL(I)
951 END DO
952
953 TABLE(1:NDDL)=0
954 DO I = 1, LCOM
955 N = IFR2K(I)
956 J = TABLE(N)
957 IF(J==0)THEN
958 INDEX(I)=0
959 TABLE(N)=I
960 ELSE
961 INDEX(I)=J
962 END IF
963 END DO
964 DEALLOCATE(TABLE)
965!dir$ omp offload target(mic:MICID)
966 & in (W:length(NDDL), free_if(.false.), align(512))
967!$OMP PARALLEL
968!$OMP END PARALLEL
969 print *,' variables11= ',i
970!dir$ omp offload target(mic:MICID)
971 & in (IFR2K,INDEX:length(LCOM), free_if(.false.))
972c & free_if(.false.), align(512))
973!$OMP PARALLEL
974!$OMP END PARALLEL
975 print *,' variables12= ',i
976 END IF
977
978 IF(NSPMD>1)THEN
979!dir$ omp offload target(mic:MICID)
980 & nocopy(IADK,IADK0,IADM,IADM0,LT_K,LT_K0,
981 & lt_m,lt_m0,JDIK,JDIK0,JDIM,JDIM0,
982 & DIAG_K,DIAG_M,x,y,z,r,p,w)
983 & inout(tt)
984!$OMP PARALLEL default(shared)
985
986
987!$OMP single
988 tt=OMP_GET_WTIME()
989 i= omp_get_num_threads()
990 print *,'number of threads(mic):',i
991!$OMP END single
992
993C-------------IT=0--------
994C----------{Z}=[K]{X}--------
995C----------{Z}=[Dk]{X}+[Lk]{X}+[Lk]^t{X}--------
996
997C D_Z = D_X * DIAG_K
998 CALL MIC_MV(NDDL,X,DIAG_K,Z)
999
1000C D_Z = D_Z + LT_K * D_X
1001 CALL MIC_SPMV(NDDL ,Z ,X ,LT_K ,IADK ,JDIK )
1002C D_Z = D_Z + LT_K0 * D_X
1003 CALL MIC_SPMV(NDDL ,Z ,X ,LT_K0,IADK0,JDIK0)
1004 IF(NDDLI > 0) THEN
1005C D_Z = D_Z + LT_I0 * D_X
1006ctempo CALL MIC_SPMV(NDDL ,Z ,X ,LT_I0,IADI0,JDII0)
1007 END IF
1008C D_R = D_R - D_Z
1009 CALL MIC_DAXPY(NDDL, -ONE, Z, R)
1010
1011C---------------------
1012C----------{Z}=[M]{R}-([M]=[Lm][Dm][Lm]^t)-------
1013C----------->{V}={R}+[Lm]^t{R}-(the first term presents the [I]{R} with is not included in Lm)-------
1014C----------->{W}=[Dm]{V}-------(the first term ->[I]{W}-same as above)
1015C----------->{Z}={W}+[Lm]{W}--------
1016C CALL PREC_SOLVGH
1017 IF(IPREC == 1)THEN
1018C D_Z = D_R
1019 CALL MIC_DCOPY(NDDL, R, Z)
1020 ELSEIF(IPREC == 2)THEN
1021C D_Z = D_R * DIAG_M
1022 CALL MIC_MV(NDDL,R,DIAG_M,Z)
1023 ELSEIF(IPREC == 5)THEN
1024C D_Y = D_R
1025 CALL MIC_DCOPY(NDDL, R, Y)
1026C D_Y = D_Y + LT_M * D_R
1027 CALL MIC_SPMV(NDDL ,Y ,R ,LT_M ,IADM ,JDIM )
1028C D_Z = D_Y * DIAG_M
1029 CALL MIC_MV(NDDL,Y,DIAG_M,Z)
1030C D_Y = D_Z
1031 CALL MIC_DCOPY(NDDL, Z, Y)
1032C D_Z = D_Z + LT_M0 * D_Y
1033 CALL MIC_SPMV(NDDL ,Z ,Y ,LT_M0,IADM0,JDIM0)
1034
1035 ELSE
1036 END IF
1037!$OMP END PARALLEL
1038C----------------------
1039.AND. IF(NSPMD > 1 IPREC > 1)THEN
1040 IF(IMONM > 0) CALL STARTIME(TIMERS,68)
1041!dir$ omp offload target(mic:MICID)
1042 & nocopy(z,ifr2k)
1043 & out(vgat:length(lcom))
1044!$OMP PARALLEL default(shared)
1045 CALL MIC_GATHER(NDDL, LCOM, Z, VGAT, IFR2K)
1046!$OMP END PARALLEL
1047 IF(IMONM > 0) CALL STOPTIME(TIMERS,68)
1048 IF(IMONM > 0) CALL STARTIME(TIMERS,66)
1049 CALL SPMD_SUMFC_V(VGAT,VSCA,INDEX,LCOM)
1050 IF(IMONM > 0) CALL STOPTIME(TIMERS,66)
1051 IF(IMONM > 0) CALL STARTIME(TIMERS,69)
1052!dir$ omp offload target(mic:MICID)
1053 & nocopy(z,ifr2k,index)
1054 & in(vsca:length(lcom))
1055!$OMP PARALLEL default(shared)
1056 CALL MIC_SCATTER(NDDL, LCOM, Z, VSCA, IFR2K, INDEX)
1057!$OMP END PARALLEL
1058 IF(IMONM > 0) CALL STOPTIME(TIMERS,69)
1059 END IF
1060C---------------------
1061!dir$ omp offload target(mic:MICID)
1062 & nocopy(z,p)
1063!$OMP PARALLEL default(shared)
1064C D_P = D_Z
1065 CALL MIC_DCOPY(NDDL, Z, P)
1066!$OMP END PARALLEL
1067C
1068 IF(NSPMD == 1) THEN
1069!dir$ omp offload target(mic:MICID)
1070 & nocopy(r,z)
1071 & inout(G0)
1072!$OMP PARALLEL default(shared) shared(g0)
1073 CALL MIC_DDOT(NDDL, R, Z, G0)
1074!$OMP END PARALLEL
1075 ELSE
1076C si spmd penser a faire avant D_Y = D_Z*W pour les poids
1077C D_Y = D_Z * D_W
1078!dir$ omp offload target(mic:MICID)
1079 & nocopy(z,y,w,r)
1080 & inout(G0)
1081!$OMP PARALLEL default(shared) shared(g0)
1082 CALL MIC_MV(NDDL,Z,W,Y)
1083 CALL MIC_DDOT(NDDL, R, Y, G0)
1084!$OMP END PARALLEL
1085 IF(IMONM > 0) CALL STARTIME(TIMERS,67)
1086 CALL SPMD_SUM_S(G0)
1087 IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
1088 ENDIF
1089
1090!dir$ omp offload target(mic:MICID)
1091 & nocopy(iadk,iadk0,iadm,iadm0,lt_k,lt_k0,lt_m)
1092 & nocopy(lt_m0,jdik,jdik0,jdim,jdim0,diag_k,diag_m,y,p,x,r,z)
1093!$OMP PARALLEL default(shared)
1094C---------------------
1095C CALL MAV_LTGH(
1096C 1 NDDL ,IADK ,JDIK ,DIAG_K,LT_K ,
1097C 2 P ,Y ,F_DDL ,L_DDL ,ITASK )
1098
1099C D_Y = D_P * DIAG_K
1100 CALL MIC_MV(NDDL,P,DIAG_K,Y)
1101C D_Y = D_Y + LT_K * D_P
1102 CALL MIC_SPMV(NDDL ,Y ,P ,LT_K ,IADK ,JDIK )
1103C D_Y = D_Y + LT_K0 * D_P
1104 CALL MIC_SPMV(NDDL ,Y ,P ,LT_K0,IADK0,JDIK0)
1105 IF(NDDLI > 0) THEN
1106C D_Y = D_Y + LT_I0 * D_P
1107ctempo CALL MIC_SPMV(NDDL ,Y ,P ,LT_I0,IADI0,JDII0)
1108 END IF
1109!$OMP END PARALLEL
1110C----------------------
1111 IF(NSPMD > 1)THEN
1112 IF(IMONM > 0) CALL STARTIME(TIMERS,68)
1113!dir$ omp offload target(mic:MICID)
1114 & nocopy(y,ifr2k)
1115 & out(vgat:length(lcom))
1116!$OMP PARALLEL default(shared)
1117 CALL MIC_GATHER(NDDL, LCOM, Y, VGAT, IFR2K)
1118!$OMP END PARALLEL
1119 IF(IMONM > 0) CALL STOPTIME(TIMERS,68)
1120 IF(IMONM > 0) CALL STARTIME(TIMERS,66)
1121 CALL SPMD_SUMFC_V(VGAT,VSCA,INDEX,LCOM)
1122 IF(IMONM > 0) CALL STOPTIME(TIMERS,66)
1123 IF(IMONM > 0) CALL STARTIME(TIMERS,69)
1124!dir$ omp offload target(mic:MICID)
1125 & nocopy(y,ifr2k,index)
1126 & in(vsca:length(lcom))
1127!$OMP PARALLEL default(shared)
1128 CALL MIC_SCATTER(NDDL, LCOM, Y, VSCA, IFR2K, INDEX)
1129!$OMP END PARALLEL
1130 IF(IMONM > 0) CALL STOPTIME(TIMERS,69)
1131 END IF
1132C---------------------
1133 IF(NSPMD == 1) THEN
1134!dir$ omp offload target(mic:MICID)
1135 & nocopy(p,y)
1136 & inout(S)
1137!$OMP PARALLEL default(shared) shared(s)
1138 CALL MIC_DDOT(NDDL, P, Y, S)
1139!$OMP END PARALLEL
1140 ELSE
1141C si spmd penser a faire avant D_Z = D_Y*D_W pour les poids
1142!dir$ omp offload target(mic:MICID)
1143 & nocopy(y,w,p,z)
1144 & inout(S)
1145!$OMP PARALLEL default(shared) shared(s)
1146 CALL MIC_MV(NDDL,Y,W,Z)
1147 CALL MIC_DDOT(NDDL, P, Z, S)
1148!$OMP END PARALLEL
1149 IF(IMONM > 0) CALL STARTIME(TIMERS,67)
1150 CALL SPMD_SUM_S(S)
1151 IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
1152 END IF
1153C---------------------
1154C---------------------
1155 IF (G0==ZERO) THEN
1156 IF (ITOL>1) THEN
1157 ISTOP=-1
1158 GOTO 206
1159 ELSE
1160 ALPHA = ZERO
1161 ENDIF
1162 ELSE
1163 ALPHA = G0/S
1164 ENDIF
1165 TOLS2=TOLS*TOLS
1166 IF (ITOL==1) THEN
1167 IF(NSPMD == 1) THEN
1168!dir$ omp offload target(mic:MICID)
1169 & nocopy(r)
1170 & inout(r02)
1171!$OMP PARALLEL default(shared) shared(r02)
1172 CALL MIC_DDOT(NDDL, R, R, R02)
1173!$OMP END PARALLEL
1174 ELSE
1175C si spmd penser a faire avant D_Z = D_R*W pour les poids
1176!dir$ omp offload target(mic:MICID)
1177 & nocopy(r,w,z)
1178 & inout(r02)
1179!$OMP PARALLEL default(shared) shared(r02)
1180 CALL MIC_MV(NDDL,R,W,Z)
1181 CALL MIC_DDOT(NDDL, Z, Z, R02)
1182!$OMP END PARALLEL
1183 IF(IMONM > 0) CALL STARTIME(TIMERS,67)
1184 CALL SPMD_SUM_S(R02)
1185 IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
1186 END IF
1187C---------------------
1188C---------------------
1189 R2 =R02
1190 ELSEIF (ITOL==3) THEN
1191C------ R2<TOL*TOL*ANORM2*XNORM2------
1192 R02=ABS(G0)
1193 R2 =R02
1194 L_A=ONE/ALPHA
1195C L_B2=R02
1196 TNORM2=L_A*L_A
1197 ANORM2=ZERO
1198 XNORM2=ZERO
1199 A_OLD=L_A
1200 B_OLD=ZERO
1201 OLDB = SQRT(R02)
1202 ELSEIF (ITOL==4) THEN
1203 R02=ALPHA*ALPHA*ABS(G0)
1204 EPS(1)=R02
1205 R2=R02
1206 ELSE
1207 R02=ABS(G0)
1208 R2 =R02
1209 ENDIF
1210 IF (R02==ZERO) THEN
1211 ISTOP=-1
1212 GOTO 206
1213 ENDIF
1214 TOLN=R02*TOLS2
1215.AND. IF(IPRI/=0ITASK==0)THEN
1216 RR = SQRT(R2)
1217 IF(IPRI<0) WRITE(ISTDO,1000)IT,RR
1218 WRITE(IOUT,1000)IT,RR
1219 ENDIF
1220C-------pour etre coherent avec lanzos for linear
1221C----------------------
1222c CALL MY_BARRIER
1223C---------------------
1224 IT=1
1225!dir$ omp offload target(mic:MICID)
1226 & nocopy(iadk,iadk0,iadm,iadm0,lt_k,lt_k0,lt_m)
1227 & nocopy(lt_m0,jdik,jdik0,jdim,jdim0,diag_k,diag_m,y,p,x,r,z)
1228!$OMP PARALLEL default(shared)
1229 CALL MIC_DAXPY(NDDL, ALPHA, P, X)
1230 CALL MIC_DAXPY(NDDL,-ALPHA, Y, R)
1231
1232C----------------------
1233c CALL MY_BARRIER
1234C---------------------
1235C CALL PREC_SOLVGH(IPREC, ITASK ,NDDL ,IADM , JDIM ,
1236C 6 DIAG_M , LT_M ,R ,Z ,
1237C 7 F_DDL ,L_DDL )
1238 IF(IPREC == 1)THEN
1239C D_Z = D_R
1240 CALL MIC_DCOPY(NDDL, R, Z)
1241 ELSEIF(IPREC == 2)THEN
1242C D_Z = D_R * DIAG_M
1243 CALL MIC_MV(NDDL,R,DIAG_M,Z)
1244 ELSEIF(IPREC == 5)THEN
1245C D_Y = D_R
1246 CALL MIC_DCOPY(NDDL, R, Y)
1247C D_Y = D_Y + LT_M * D_R
1248 CALL MIC_SPMV(NDDL ,Y ,R ,LT_M ,IADM ,JDIM )
1249C D_Z = D_Y * DIAG_M
1250 CALL MIC_MV(NDDL,Y,DIAG_M,Z)
1251 CALL MIC_DCOPY(NDDL, Z, Y)
1252C D_Z = D_Z + LT_M0 * D_Y
1253 CALL MIC_SPMV(NDDL ,Z ,Y ,LT_M0,IADM0,JDIM0)
1254 ELSE
1255 END IF
1256!$OMP END PARALLEL
1257C----------------------
1258.AND. IF(NSPMD > 1 IPREC > 1)THEN
1259 IF(IMONM > 0) CALL STARTIME(TIMERS,68)
1260!dir$ omp offload target(mic:MICID)
1261 & nocopy(z,ifr2k)
1262 & out(vgat:length(lcom))
1263!$OMP PARALLEL default(shared)
1264 CALL MIC_GATHER(NDDL, LCOM, Z, VGAT, IFR2K)
1265!$OMP END PARALLEL
1266 IF(IMONM > 0) CALL STOPTIME(TIMERS,68)
1267 IF(IMONM > 0) CALL STARTIME(TIMERS,66)
1268 CALL SPMD_SUMFC_V(VGAT,VSCA,INDEX,LCOM)
1269 IF(IMONM > 0) CALL STOPTIME(TIMERS,66)
1270 IF(IMONM > 0) CALL STARTIME(TIMERS,69)
1271!dir$ omp offload target(mic:MICID)
1272 & nocopy(z,ifr2k,index)
1273 & in(vsca:length(lcom))
1274!$OMP PARALLEL default(shared)
1275 CALL MIC_SCATTER(NDDL, LCOM, Z, VSCA, IFR2K, INDEX)
1276!$OMP END PARALLEL
1277 IF(IMONM > 0) CALL STOPTIME(TIMERS,69)
1278 END IF
1279C---------------------
1280 IF(NSPMD == 1) THEN
1281!dir$ omp offload target(mic:MICID)
1282 & nocopy(r,z)
1283 & inout(G1)
1284!$OMP PARALLEL default(shared) shared(G1)
1285 CALL MIC_DDOT(NDDL, R, Z, G1)
1286!$OMP END PARALLEL
1287 ELSE
1288C si spmd penser a faire avant D_Y = D_Z*D_W pour les poids
1289!dir$ omp offload target(mic:MICID)
1290 & nocopy(z,w,r,y)
1291 & inout(G1)
1292!$OMP PARALLEL default(shared) shared(G1)
1293 CALL MIC_MV(NDDL,Z,W,Y)
1294 CALL MIC_DDOT(NDDL, R, Y, G1)
1295!$OMP END PARALLEL
1296 IF(IMONM > 0) CALL STARTIME(TIMERS,67)
1297 CALL SPMD_SUM_S(G1)
1298 IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
1299 END IF
1300C---------------------
1301 BETA=G1/G0
1302 IF (ITOL==1) THEN
1303 IF(NSPMD == 1) THEN
1304!dir$ omp offload target(mic:MICID)
1305 & nocopy(r)
1306 & inout(R2)
1307!$OMP PARALLEL default(shared) shared(R2)
1308 CALL MIC_DDOT(NDDL, R, R, R2)
1309!$OMP END PARALLEL
1310 ELSE
1311C si spmd penser a faire avant D_Y = D_R*W pour les poids
1312!dir$ omp offload target(mic:MICID)
1313 & nocopy(r,w,y)
1314 & inout(R2)
1315!$OMP PARALLEL default(shared) shared(R2)
1316 CALL MIC_MV(NDDL,R,W,Y)
1317 CALL MIC_DDOT(NDDL, Y, Y, R2)
1318!$OMP END PARALLEL
1319 IF(IMONM > 0) CALL STARTIME(TIMERS,67)
1320 CALL SPMD_SUM_S(R2)
1321 IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
1322 END IF
1323 ELSEIF (ITOL==3) THEN
1324C------ R2<TOL*TOL*ANORM2*XNORM2------
1325 R2=ABS(G1)
1326 L_B2=ABS(BETA)*A_OLD*A_OLD
1327 L_B=SQRT(L_B2)
1328 TNORM2=TNORM2+L_B2
1329 B_OLD=BETA
1330C INITIALIZE OTHER QUANTITIES.
1331 GBAR = L_A
1332 DBAR = L_B
1333 RHS1 = OLDB
1334 RHS2 = ZERO
1335 GMAX = ABS( L_A ) + EPS_M
1336 GMIN = GMAX
1337 OLDB2 = L_B2
1338 OLDB = L_B
1339 ELSEIF (ITOL==4) THEN
1340 R2=R02
1341 ELSE
1342 R2=ABS(G1)
1343 ENDIF
1344 G0 = G1
1345 IF (ITOL==3) TOLN=TOLN*TNORM2
1346C----------------------
1347c CALL MY_BARRIER
1348C---------------------
1349 ISTOP=CRIT_STOP(IT,R2,NLIM,TOLN)
1350
1351C----------------------
1352c LOOP OVER ITERATIONS
1353C---------------------
1354
1355 DO WHILE (ISTOP==1)
1356
1357!dir$ omp offload target(mic:MICID)
1358 & nocopy(iadk,iadk0,iadm,iadm0,lt_k,lt_k0,lt_m)
1359 & nocopy(lt_m0,jdik,jdik0,jdim,jdim0,diag_k,diag_m,y,p,x,r,z)
1360!$OMP PARALLEL default(shared)
1361
1362C D_P = D_Z + BETA*D_P in 2 step : D_P = BETA*D_P ; D_P = D_P + D_Z
1363 CALL MIC_DSCAL(NDDL, BETA, P)
1364 CALL MIC_DAXPY(NDDL, ONE, Z, P)
1365C----------------------
1366c CALL MY_BARRIER
1367C---------------------
1368C CALL MAV_LTGH(
1369C 1 NDDL ,IADK ,JDIK ,DIAG_K,LT_K ,
1370C 2 P ,Y ,F_DDL ,L_DDL ,ITASK )
1371C D_Y = D_P * DIAG_K
1372 CALL MIC_MV(NDDL,P,DIAG_K,Y)
1373C D_Y = D_Y + LT_K * D_P
1374 CALL MIC_SPMV(NDDL ,Y ,P ,LT_K ,IADK ,JDIK )
1375C D_Y = D_Y + LT_K0 * D_P
1376 CALL MIC_SPMV(NDDL ,Y ,P ,LT_K0,IADK0,JDIK0)
1377 IF(NDDLI > 0) THEN
1378C D_Y = D_Y + LT_I0 * D_P
1379ctempo CALL MIC_SPMV(NDDL ,Y ,P ,LT_I0,IADI0,JDII0)
1380 END IF
1381!$OMP END PARALLEL
1382 IF(NSPMD > 1)THEN
1383 IF(IMONM > 0) CALL STARTIME(TIMERS,68)
1384!dir$ omp offload target(mic:MICID)
1385 & nocopy(y,ifr2k)
1386 & out(vgat:length(lcom))
1387!$OMP PARALLEL default(shared)
1388 CALL MIC_GATHER(NDDL, LCOM, Y, VGAT, IFR2K)
1389!$OMP END PARALLEL
1390 IF(IMONM > 0) CALL STOPTIME(TIMERS,68)
1391 IF(IMONM > 0) CALL STARTIME(TIMERS,66)
1392 CALL SPMD_SUMFC_V(VGAT,VSCA,INDEX,LCOM)
1393 IF(IMONM > 0) CALL STOPTIME(TIMERS,66)
1394 IF(IMONM > 0) CALL STARTIME(TIMERS,69)
1395!dir$ omp offload target(mic:MICID)
1396 & nocopy(y,ifr2k,index)
1397 & in(vsca:length(lcom))
1398!$OMP PARALLEL default(shared)
1399 CALL MIC_SCATTER(NDDL, LCOM, Y, VSCA, IFR2K, INDEX)
1400!$OMP END PARALLEL
1401 IF(IMONM > 0) CALL STOPTIME(TIMERS,69)
1402 END IF
1403C---------------------
1404 IF(NSPMD == 1) THEN
1405!dir$ omp offload target(mic:MICID)
1406 & nocopy(p,y)
1407 & inout(S)
1408!$OMP PARALLEL default(shared) shared(S)
1409 CALL MIC_DDOT(NDDL, P, Y, S)
1410!$OMP END PARALLEL
1411 ELSE
1412C si spmd penser a faire avant D_Y = D_Y*W pour les poids
1413!dir$ omp offload target(mic:MICID)
1414 & nocopy(y,w,p,z)
1415 & inout(S)
1416!$OMP PARALLEL default(shared) shared(S)
1417 CALL MIC_MV(NDDL,Y,W,Z)
1418 CALL MIC_DDOT(NDDL, P, Z, S)
1419!$OMP END PARALLEL
1420 IF(IMONM > 0) CALL STARTIME(TIMERS,67)
1421 CALL SPMD_SUM_S(S)
1422 IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
1423 END IF
1424C---------------------
1425 ALPHA=G0/S
1426!dir$ omp offload target(mic:MICID)
1427 & nocopy(iadk,iadk0,iadm,iadm0,lt_k,lt_k0,lt_m)
1428 & nocopy(lt_m0,jdik,jdik0,jdim,jdim0,diag_k,diag_m,y,p,x,r,z)
1429!$OMP PARALLEL default(shared)
1430 CALL MIC_DAXPY(NDDL, ALPHA, P, X)
1431 CALL MIC_DAXPY(NDDL,-ALPHA, Y, R)
1432C----------------------
1433C CALL PREC_SOLVGH(IPREC , ITASK ,NDDL ,IADM , JDIM ,
1434C 6 DIAG_M, LT_M ,R ,Z ,
1435C 7 F_DDL ,L_DDL )
1436 IF(IPREC == 1)THEN
1437C D_Z = D_R
1438 CALL MIC_DCOPY(NDDL, R, Z)
1439 ELSEIF(IPREC == 2)THEN
1440C D_Z = D_R * DIAG_M
1441 CALL MIC_MV(NDDL,R,DIAG_M,Z)
1442 ELSEIF(IPREC == 5)THEN
1443C D_Y = D_R
1444 CALL MIC_DCOPY(NDDL, R, Y)
1445C D_Y = D_Y + LT_M * D_R
1446 CALL MIC_SPMV(NDDL ,Y ,R ,LT_M ,IADM ,JDIM )
1447C D_Z = D_Y * DIAG_M
1448 CALL MIC_MV(NDDL,Y,DIAG_M,Z)
1449 CALL MIC_DCOPY(NDDL, Z, Y)
1450C D_Z = D_Z + LT_M0 * D_Y
1451 CALL MIC_SPMV(NDDL ,Z ,Y ,LT_M0,IADM0,JDIM0)
1452 ELSE
1453 END IF
1454!$OMP END PARALLEL
1455C----------------------
1456.AND. IF(NSPMD > 1 IPREC > 1)THEN
1457 IF(IMONM > 0) CALL STARTIME(TIMERS,68)
1458!dir$ omp offload target(mic:MICID)
1459 & nocopy(z,ifr2k)
1460 & out(vgat:length(lcom))
1461!$OMP PARALLEL default(shared)
1462 CALL MIC_GATHER(NDDL, LCOM, Z, VGAT, IFR2K)
1463!$OMP END PARALLEL
1464 IF(IMONM > 0) CALL STOPTIME(TIMERS,68)
1465 IF(IMONM > 0) CALL STARTIME(TIMERS,66)
1466 CALL SPMD_SUMFC_V(VGAT,VSCA,INDEX,LCOM)
1467 IF(IMONM > 0) CALL STOPTIME(TIMERS,66)
1468 IF(IMONM > 0) CALL STARTIME(TIMERS,69)
1469!dir$ omp offload target(mic:MICID)
1470 & nocopy(z,ifr2k,index)
1471 & in(vsca:length(lcom))
1472!$OMP PARALLEL default(shared)
1473 CALL MIC_SCATTER(NDDL, LCOM, Z, VSCA, IFR2K, INDEX)
1474!$OMP END PARALLEL
1475 IF(IMONM > 0) CALL STOPTIME(TIMERS,69)
1476 END IF
1477C---------------------
1478 IF(NSPMD == 1) THEN
1479!dir$ omp offload target(mic:MICID)
1480 & nocopy(r,z)
1481 & inout(G1)
1482!$OMP PARALLEL default(shared) shared(G1)
1483 CALL MIC_DDOT(NDDL, R, Z, G1)
1484!$OMP END PARALLEL
1485 ELSE
1486C si spmd penser a faire avant D_Y = D_Z*W pour les poids
1487!dir$ omp offload target(mic:MICID)
1488 & nocopy(z,w,r,y)
1489 & inout(G1)
1490!$OMP PARALLEL default(shared) shared(G1)
1491 CALL MIC_MV(NDDL,Z,W,Y)
1492 CALL MIC_DDOT(NDDL, R, Y, G1)
1493!$OMP END PARALLEL
1494 IF(IMONM > 0) CALL STARTIME(TIMERS,67)
1495 CALL SPMD_SUM_S(G1)
1496 IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
1497 END IF
1498C---------------------
1499cc!$OMP single
1500 BETA=G1/G0
1501 G0 = G1
1502cc!$OMP end single
1503 IF (ITOL==1) THEN
1504 IF(NSPMD == 1) THEN
1505!dir$ omp offload target(mic:MICID)
1506 & nocopy(r)
1507 & inout(R2)
1508!$OMP PARALLEL default(shared) shared(R2)
1509 CALL MIC_DDOT(NDDL, R, R, R2)
1510!$OMP END PARALLEL
1511 ELSE
1512C si spmd penser a faire avant D_Y = D_R*W pour les poids
1513!dir$ omp offload target(mic:MICID)
1514 & nocopy(r,w,y)
1515 & inout(R2)
1516!$OMP PARALLEL default(shared) shared(R2)
1517 CALL MIC_MV(NDDL,R,W,Y)
1518 CALL MIC_DDOT(NDDL, Y, Y, R2)
1519!$OMP END PARALLEL
1520 IF(IMONM > 0) CALL STARTIME(TIMERS,67)
1521 CALL SPMD_SUM_S(R2)
1522 IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
1523 END IF
1524 ELSEIF (ITOL==3) THEN
1525 R2 =ABS(G1)
1526 S=ONE/ALPHA
1527 L_A=S+B_OLD*A_OLD
1528C----------L_B2 : beta(j-1)^2-A_OLD :1/alpha(j-1)-------
1529 L_B2=ABS(BETA)*S*S
1530 L_B=SQRT(L_B2)
1531 A_OLD=S
1532 B_OLD=BETA
1533 ANORM2=TNORM2
1534 TNORM2=TNORM2+L_A*L_A+OLDB2+L_B2
1535 GAMMA = SQRT( GBAR*GBAR + OLDB2 )
1536 CS = GBAR / GAMMA
1537 SN = OLDB / GAMMA
1538 DELTA = CS * DBAR + SN * L_A
1539 GBAR = SN * DBAR - CS * L_A
1540 EPSLN = SN * L_B
1541 DBAR = - CS * L_B
1542 ZL = RHS1 / GAMMA
1543 XNORM2 = XNORM2+ZL*ZL
1544 GMAX = MAX( GMAX, GAMMA )
1545 GMIN = MIN( GMIN, GAMMA )
1546 RHS1 = RHS2 - DELTA * ZL
1547 RHS2 = - EPSLN * ZL
1548 TOLN=TOLS2*ANORM2*XNORM2
1549 OLDB2 = L_B2
1550 OLDB = L_B
1551 ELSEIF (ITOL==4) THEN
1552 TMP=ALPHA*ALPHA*ABS(G1)
1553 IF (IT>=ND) THEN
1554 DO J=1,ND-1
1555 EPS(J)=EPS(J+1)
1556 ENDDO
1557 EPS(ND)=TMP
1558 R2=ZERO
1559 DO J=1,ND
1560 R2=R2+EPS(J)
1561 ENDDO
1562 ELSE
1563 EPS(IT+1)=TMP
1564 R2=R2+TMP
1565 ENDIF
1566 ELSE
1567 R2 =ABS(G1)
1568 ENDIF
1569.AND. IF(IPRI/=0ITASK==0)THEN
1570 IF(MOD(IT,IP)==0)THEN
1571 RR = SQRT(R2/R02)
1572 WRITE(IOUT,1001)IT,RR
1573 IF(IPRI<0) WRITE(ISTDO,1001)IT,RR
1574 ENDIF
1575 ENDIF
1576C----------------------
1577c CALL MY_BARRIER
1578C---------------------
1579 ISTOP=CRIT_STOP(IT,R2,NLIM,TOLN)
1580C----------------------
1581c CALL MY_BARRIER
1582C---------------------
1583 IT = IT +1
1584 ENDDO
1585 206 CONTINUE
1586
1587 WRITE(IOUT,1001)IT-1,RR
1588
1589C!dec$ omp offload target(mic:MICID) out(X,R:alloc_if(.false.))
1590!$OMP PARALLEL
1591!$OMP single
1592 tt=OMP_GET_WTIME()-tt
1593 i= omp_get_num_threads()
1594 print *,' '
1595 print *,' '
1596 print *,' execution time on mic with',i,' threads'
1597 print *,' ==> ',tt,' seconds'
1598 print *,' '
1599!$OMP END single
1600!$OMP END PARALLEL
1601 DEALLOCATE(W)
1602 DEALLOCATE(VGAT)
1603 DEALLOCATE(VSCA)
1604 DEALLOCATE(INDEX)
1605c+1 tempo
1606 ELSE ! CODE NSPMD = 1
1607C debut code provisoire
1608!dec$ attributes offload: mic :: ONE,ZERO,ISTDO,IOUT
1609
1610!dir$ omp offload target(mic:MICID)
1611 & nocopy(IADK,IADK0,IADM,IADM0,LT_K,LT_K0,
1612 & lt_m,lt_m0,JDIK,JDIK0,JDIM,JDIM0,
1613 & DIAG_K,DIAG_M,x,y,z,r,p,w)
1614 & inout(tt)
1615
1616!$OMP PARALLEL default(shared)
1617
1618!$OMP single
1619 tt=OMP_GET_WTIME()
1620 i= omp_get_num_threads()
1621 print *,'number of threads(mic):',i
1622!$OMP END single
1623
1624C-------------IT=0--------
1625C----------{Z}=[K]{X}--------
1626C----------{Z}=[Dk]{X}+[Lk]{X}+[Lk]^t{X}--------
1627
1628C D_Z = D_X * DIAG_K
1629 CALL MIC_MV(NDDL,X,DIAG_K,Z)
1630
1631C D_Z = D_Z + LT_K * D_X
1632 CALL MIC_SPMV(NDDL ,Z ,X ,LT_K ,IADK ,JDIK )
1633C D_Z = D_Z + LT_K0 * D_X
1634 CALL MIC_SPMV(NDDL ,Z ,X ,LT_K0,IADK0,JDIK0)
1635 IF(NDDLI > 0) THEN
1636C D_Z = D_Z + LT_I0 * D_X
1637ctempo CALL MIC_SPMV(NDDL ,Z ,X ,LT_I0,IADI0,JDII0)
1638 END IF
1639C D_R = D_R - D_Z
1640 CALL MIC_DAXPY(NDDL, -ONE, Z, R)
1641C---------------------
1642C----------{Z}=[M]{R}-([M]=[Lm][Dm][Lm]^t)-------
1643C----------->{V}={R}+[Lm]^t{R}-(the first term presents the [I]{R} with is not included in Lm)-------
1644C----------->{W}=[Dm]{V}-------(the first term ->[I]{W}-same as above)
1645C----------->{Z}={W}+[Lm]{W}--------
1646C CALL PREC_SOLVGH
1647 IF(IPREC == 1)THEN
1648C D_Z = D_R
1649 CALL MIC_DCOPY(NDDL, R, Z)
1650 ELSEIF(IPREC == 2)THEN
1651C D_Z = D_R * DIAG_M
1652 CALL MIC_MV(NDDL,R,DIAG_M,Z)
1653 ELSEIF(IPREC == 5)THEN
1654C D_Y = D_R
1655 CALL MIC_DCOPY(NDDL, R, Y)
1656C D_Y = D_Y + LT_M * D_R
1657 CALL MIC_SPMV(NDDL ,Y ,R ,LT_M ,IADM ,JDIM )
1658C D_Z = D_Y * DIAG_M
1659 CALL MIC_MV(NDDL,Y,DIAG_M,Z)
1660C D_Y = D_Z
1661 CALL MIC_DCOPY(NDDL, Z, Y)
1662C D_Z = D_Z + LT_M0 * D_Y
1663 CALL MIC_SPMV(NDDL ,Z ,Y ,LT_M0,IADM0,JDIM0)
1664
1665 ELSE
1666 END IF
1667c!$OMP END PARALLEL
1668C----------------------
1669c IF(NSPMD > 1 .AND. IPREC > 1)THEN
1670c IF(IMONM > 0) CALL STARTIME(TIMERS,68)
1671c!dir$ omp offload target(mic:MICID)
1672c & nocopy(z,ifr2k)
1673c & out(vgat:length(lcom))
1674c!$OMP PARALLEL default(shared)
1675c CALL MIC_GATHER(NDDL, LCOM, Z, VGAT, IFR2K)
1676c!$OMP END PARALLEL
1677c IF(IMONM > 0) CALL STOPTIME(TIMERS,68)
1678c IF(IMONM > 0) CALL STARTIME(TIMERS,66)
1679c CALL SPMD_SUMFC_V(VGAT,VSCA,INDEX,LCOM)
1680c IF(IMONM > 0) CALL STOPTIME(TIMERS,66)
1681c IF(IMONM > 0) CALL STARTIME(TIMERS,69)
1682c!dir$ omp offload target(mic:MICID)
1683c & nocopy(z,ifr2k,index)
1684c & in(vsca:length(lcom))
1685c!$OMP PARALLEL default(shared)
1686c CALL MIC_SCATTER(NDDL, LCOM, Z, VSCA, IFR2K, INDEX)
1687c!$OMP END PARALLEL
1688c IF(IMONM > 0) CALL STOPTIME(TIMERS,69)
1689c END IF
1690C---------------------
1691c!dir$ omp offload target(mic:MICID)
1692c & nocopy(z,p)
1693c!$OMP PARALLEL default(shared)
1694C D_P = D_Z
1695 CALL MIC_DCOPY(NDDL, Z, P)
1696c!$OMP END PARALLEL
1697C
1698c IF(NSPMD == 1) THEN
1699c!dir$ omp offload target(mic:MICID)
1700c & nocopy(r,z)
1701c & inout(G0)
1702c!$OMP PARALLEL default(shared) shared(g0)
1703 CALL MIC_DDOT(NDDL, R, Z, G0)
1704c!$OMP END PARALLEL
1705c ELSE
1706C si spmd penser a faire avant D_Y = D_Z*W pour les poids
1707C D_Y = D_Z * D_W
1708c!dir$ omp offload target(mic:MICID)
1709c & nocopy(z,y,w,r)
1710c & inout(G0)
1711c!$OMP PARALLEL default(shared) shared(g0)
1712c CALL MIC_MV(NDDL,Z,W,Y)
1713c CALL MIC_DDOT(NDDL, R, Y, G0)
1714c!$OMP END PARALLEL
1715c IF(IMONM > 0) CALL STARTIME(TIMERS,67)
1716c CALL SPMD_SUM_S(G0)
1717c IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
1718c ENDIF
1719
1720c!dir$ omp offload target(mic:MICID)
1721c & nocopy(iadk,iadk0,iadm,iadm0,lt_k,lt_k0,lt_m)
1722c & nocopy(lt_m0,jdik,jdik0,jdim,jdim0,diag_k,diag_m,y,p,x,r,z)
1723c!$OMP PARALLEL default(shared)
1724C---------------------
1725C CALL MAV_LTGH(
1726C 1 NDDL ,IADK ,JDIK ,DIAG_K,LT_K ,
1727C 2 P ,Y ,F_DDL ,L_DDL ,ITASK )
1728
1729C D_Y = D_P * DIAG_K
1730 CALL MIC_MV(NDDL,P,DIAG_K,Y)
1731C D_Y = D_Y + LT_K * D_P
1732 CALL MIC_SPMV(NDDL ,Y ,P ,LT_K ,IADK ,JDIK )
1733C D_Y = D_Y + LT_K0 * D_P
1734 CALL MIC_SPMV(NDDL ,Y ,P ,LT_K0,IADK0,JDIK0)
1735 IF(NDDLI > 0) THEN
1736C D_Y = D_Y + LT_I0 * D_P
1737ctempo CALL MIC_SPMV(NDDL ,Y ,P ,LT_I0,IADI0,JDII0)
1738 END IF
1739c!$OMP END PARALLEL
1740C----------------------
1741c IF(NSPMD > 1)THEN
1742c IF(IMONM > 0) CALL STARTIME(TIMERS,68)
1743c!dir$ omp offload target(mic:MICID)
1744c & nocopy(y,ifr2k)
1745c & out(vgat:length(lcom))
1746c!$OMP PARALLEL default(shared)
1747c CALL MIC_GATHER(NDDL, LCOM, Y, VGAT, IFR2K)
1748c!$OMP END PARALLEL
1749c IF(IMONM > 0) CALL STOPTIME(TIMERS,68)
1750c IF(IMONM > 0) CALL STARTIME(TIMERS,66)
1751c CALL SPMD_SUMFC_V(VGAT,VSCA,INDEX,LCOM)
1752c IF(IMONM > 0) CALL STOPTIME(TIMERS,66)
1753c IF(IMONM > 0) CALL STARTIME(TIMERS,69)
1754c!dir$ omp offload target(mic:MICID)
1755c & nocopy(y,ifr2k,index)
1756c & in(vsca:length(lcom))
1757c!$OMP PARALLEL default(shared)
1758c CALL MIC_SCATTER(NDDL, LCOM, Y, VSCA, IFR2K, INDEX)
1759c!$OMP END PARALLEL
1760c IF(IMONM > 0) CALL STOPTIME(TIMERS,69)
1761c END IF
1762C---------------------
1763c IF(NSPMD == 1) THEN
1764c!dir$ omp offload target(mic:MICID)
1765c & nocopy(p,y)
1766c & inout(S)
1767c!$OMP PARALLEL default(shared) shared(s)
1768 CALL MIC_DDOT(NDDL, P, Y, S)
1769c!$OMP END PARALLEL
1770c ELSE
1771C si spmd penser a faire avant D_Z = D_Y*D_W pour les poids
1772c!dir$ omp offload target(mic:MICID)
1773c & nocopy(y,w,p,z)
1774c & inout(S)
1775c!$OMP PARALLEL default(shared) shared(s)
1776c CALL MIC_MV(NDDL,Y,W,Z)
1777c CALL MIC_DDOT(NDDL, P, Z, S)
1778c!$OMP END PARALLEL
1779c IF(IMONM > 0) CALL STARTIME(TIMERS,67)
1780c CALL SPMD_SUM_S(S)
1781c IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
1782c END IF
1783C---------------------
1784C---------------------
1785 IF (G0==ZERO) THEN
1786 IF (ITOL>1) THEN
1787 ISTOP=-1
1788 GOTO 2060
1789 ELSE
1790 ALPHA = ZERO
1791 ENDIF
1792 ELSE
1793 ALPHA = G0/S
1794 ENDIF
1795 TOLS2=TOLS*TOLS
1796 IF (ITOL==1) THEN
1797c IF(NSPMD == 1) THEN
1798c!dir$ omp offload target(mic:MICID)
1799c & nocopy(r)
1800c & inout(r02)
1801c!$OMP PARALLEL default(shared) shared(r02)
1802 CALL MIC_DDOT(NDDL, R, R, R02)
1803c!$OMP END PARALLEL
1804c ELSE
1805C si spmd penser a faire avant D_Z = D_R*W pour les poids
1806c!dir$ omp offload target(mic:MICID)
1807c & nocopy(r,w,z)
1808c & inout(r02)
1809c!$OMP PARALLEL default(shared) shared(r02)
1810c CALL MIC_MV(NDDL,R,W,Z)
1811c CALL MIC_DDOT(NDDL, Z, Z, R02)
1812c!$OMP END PARALLEL
1813c IF(IMONM > 0) CALL STARTIME(TIMERS,67)
1814c CALL SPMD_SUM_S(R02)
1815c IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
1816c END IF
1817C---------------------
1818C---------------------
1819 R2 =R02
1820 ELSEIF (ITOL==3) THEN
1821C------ R2<TOL*TOL*ANORM2*XNORM2------
1822c+1
1823!$OMP SINGLE
1824 R02=ABS(G0)
1825 R2 =R02
1826 L_A=ONE/ALPHA
1827C L_B2=R02
1828 TNORM2=L_A*L_A
1829 ANORM2=ZERO
1830 XNORM2=ZERO
1831 A_OLD=L_A
1832 B_OLD=ZERO
1833 OLDB = SQRT(R02)
1834c+1
1835!$OMP END SINGLE
1836 ELSEIF (ITOL==4) THEN
1837 R02=ALPHA*ALPHA*ABS(G0)
1838 EPS(1)=R02
1839 R2=R02
1840 ELSE
1841 R02=ABS(G0)
1842 R2 =R02
1843 ENDIF
1844 IF (R02==ZERO) THEN
1845 ISTOP=-1
1846 GOTO 2060
1847 ENDIF
1848c+1
1849!$OMP SINGLE
1850 TOLN=R02*TOLS2
1851.AND. IF(IPRI/=0ITASK==0)THEN
1852 RR = SQRT(R2)
1853 IF(IPRI<0) WRITE(ISTDO,1000)IT,RR
1854 WRITE(IOUT,1000)IT,RR
1855 ENDIF
1856c+1
1857!$OMP END SINGLE
1858C-------pour etre coherent avec lanzos for linear
1859C----------------------
1860c CALL MY_BARRIER
1861C---------------------
1862 IT=1
1863c!dir$ omp offload target(mic:MICID)
1864c & nocopy(iadk,iadk0,iadm,iadm0,lt_k,lt_k0,lt_m)
1865c & nocopy(lt_m0,jdik,jdik0,jdim,jdim0,diag_k,diag_m,y,p,x,r,z)
1866c!$OMP PARALLEL default(shared)
1867 CALL MIC_DAXPY(NDDL, ALPHA, P, X)
1868 CALL MIC_DAXPY(NDDL,-ALPHA, Y, R)
1869
1870C----------------------
1871c CALL MY_BARRIER
1872C---------------------
1873C CALL PREC_SOLVGH(IPREC, ITASK ,NDDL ,IADM , JDIM ,
1874C 6 DIAG_M , LT_M ,R ,Z ,
1875C 7 F_DDL ,L_DDL )
1876 IF(IPREC == 1)THEN
1877C D_Z = D_R
1878 CALL MIC_DCOPY(NDDL, R, Z)
1879 ELSEIF(IPREC == 2)THEN
1880C D_Z = D_R * DIAG_M
1881 CALL MIC_MV(NDDL,R,DIAG_M,Z)
1882 ELSEIF(IPREC == 5)THEN
1883C D_Y = D_R
1884 CALL MIC_DCOPY(NDDL, R, Y)
1885C D_Y = D_Y + LT_M * D_R
1886 CALL MIC_SPMV(NDDL ,Y ,R ,LT_M ,IADM ,JDIM )
1887C D_Z = D_Y * DIAG_M
1888 CALL MIC_MV(NDDL,Y,DIAG_M,Z)
1889 CALL MIC_DCOPY(NDDL, Z, Y)
1890C D_Z = D_Z + LT_M0 * D_Y
1891 CALL MIC_SPMV(NDDL ,Z ,Y ,LT_M0,IADM0,JDIM0)
1892 ELSE
1893 END IF
1894c!$OMP END PARALLEL
1895C----------------------
1896c IF(NSPMD > 1 .AND. IPREC > 1)THEN
1897c IF(IMONM > 0) CALL STARTIME(TIMERS,68)
1898c!dir$ omp offload target(mic:MICID)
1899c & nocopy(z,ifr2k)
1900c & out(vgat:length(lcom))
1901c!$OMP PARALLEL default(shared)
1902c CALL MIC_GATHER(NDDL, LCOM, Z, VGAT, IFR2K)
1903c!$OMP END PARALLEL
1904c IF(IMONM > 0) CALL STOPTIME(TIMERS,68)
1905c IF(IMONM > 0) CALL STARTIME(TIMERS,66)
1906c CALL SPMD_SUMFC_V(VGAT,VSCA,INDEX,LCOM)
1907c IF(IMONM > 0) CALL STOPTIME(TIMERS,66)
1908c IF(IMONM > 0) CALL STARTIME(TIMERS,69)
1909c!dir$ omp offload target(mic:MICID)
1910c & nocopy(z,ifr2k,index)
1911c & in(vsca:length(lcom))
1912c!$OMP PARALLEL default(shared)
1913c CALL MIC_SCATTER(NDDL, LCOM, Z, VSCA, IFR2K, INDEX)
1914c!$OMP END PARALLEL
1915c IF(IMONM > 0) CALL STOPTIME(TIMERS,69)
1916c END IF
1917C---------------------
1918c IF(NSPMD == 1) THEN
1919c!dir$ omp offload target(mic:MICID)
1920c & nocopy(r,z)
1921c & inout(G1)
1922c!$OMP PARALLEL default(shared) shared(G1)
1923 CALL MIC_DDOT(NDDL, R, Z, G1)
1924c!$OMP END PARALLEL
1925c ELSE
1926C si spmd penser a faire avant D_Y = D_Z*D_W pour les poids
1927c!dir$ omp offload target(mic:MICID)
1928c & nocopy(z,w,r,y)
1929c & inout(G1)
1930c!$OMP PARALLEL default(shared) shared(G1)
1931c CALL MIC_MV(NDDL,Z,W,Y)
1932c CALL MIC_DDOT(NDDL, R, Y, G1)
1933c!$OMP END PARALLEL
1934c IF(IMONM > 0) CALL STARTIME(TIMERS,67)
1935c CALL SPMD_SUM_S(G1)
1936c IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
1937c END IF
1938C---------------------
1939 BETA=G1/G0
1940 IF (ITOL==1) THEN
1941c IF(NSPMD == 1) THEN
1942c!dir$ omp offload target(mic:MICID)
1943c & nocopy(r)
1944c & inout(R2)
1945c!$OMP PARALLEL default(shared) shared(R2)
1946 CALL MIC_DDOT(NDDL, R, R, R2)
1947c!$OMP END PARALLEL
1948c ELSE
1949C si spmd penser a faire avant D_Y = D_R*W pour les poids
1950c!dir$ omp offload target(mic:MICID)
1951c & nocopy(r,w,y)
1952c & inout(R2)
1953c!$OMP PARALLEL default(shared) shared(R2)
1954c CALL MIC_MV(NDDL,R,W,Y)
1955c CALL MIC_DDOT(NDDL, Y, Y, R2)
1956c!$OMP END PARALLEL
1957c IF(IMONM > 0) CALL STARTIME(TIMERS,67)
1958c CALL SPMD_SUM_S(R2)
1959c IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
1960c END IF
1961 ELSEIF (ITOL==3) THEN
1962C------ R2<TOL*TOL*ANORM2*XNORM2------
1963c+1
1964!$OMP SINGLE
1965 R2=ABS(G1)
1966 L_B2=ABS(BETA)*A_OLD*A_OLD
1967 L_B=SQRT(L_B2)
1968 TNORM2=TNORM2+L_B2
1969 B_OLD=BETA
1970C INITIALIZE OTHER QUANTITIES.
1971 GBAR = L_A
1972 DBAR = L_B
1973 RHS1 = OLDB
1974 RHS2 = ZERO
1975 GMAX = ABS( L_A ) + EPS_M
1976 GMIN = GMAX
1977 OLDB2 = L_B2
1978 OLDB = L_B
1979c+1
1980!$OMP END SINGLE
1981 ELSEIF (ITOL==4) THEN
1982 R2=R02
1983 ELSE
1984 R2=ABS(G1)
1985 ENDIF
1986 G0 = G1
1987c+1
1988!$OMP SINGLE
1989 IF (ITOL==3) TOLN=TOLN*TNORM2
1990
1991
1992C----------------------
1993c CALL MY_BARRIER
1994C---------------------
1995C ISTOP=CRIT_STOP(IT,R2,NLIM,TOLN)
1996 IF (IT>=NLIM) THEN
1997 ISTOP = 0
1998 ELSE
1999 IF (R2<=TOLN) THEN
2000 ISTOP = 0
2001 ELSE
2002 ISTOP = 1
2003 ENDIF
2004 ENDIF
2005c+1
2006!$OMP END SINGLE
2007C----------------------
2008c LOOP OVER ITERATIONS
2009C---------------------
2010
2011 DO WHILE (ISTOP==1)
2012
2013c!dir$ omp offload target(mic:MICID)
2014c & nocopy(iadk,iadk0,iadm,iadm0,lt_k,lt_k0,lt_m)
2015c & nocopy(lt_m0,jdik,jdik0,jdim,jdim0,diag_k,diag_m,y,p,x,r,z)
2016c!$OMP PARALLEL default(shared)
2017
2018C D_P = D_Z + BETA*D_P in 2 step : D_P = BETA*D_P ; D_P = D_P + D_Z
2019 CALL MIC_DSCAL(NDDL, BETA, P)
2020 CALL MIC_DAXPY(NDDL, ONE, Z, P)
2021C----------------------
2022c CALL MY_BARRIER
2023C---------------------
2024C CALL MAV_LTGH(
2025C 1 NDDL ,IADK ,JDIK ,DIAG_K,LT_K ,
2026C 2 P ,Y ,F_DDL ,L_DDL ,ITASK )
2027C D_Y = D_P * DIAG_K
2028 CALL MIC_MV(NDDL,P,DIAG_K,Y)
2029C D_Y = D_Y + LT_K * D_P
2030 CALL MIC_SPMV(NDDL ,Y ,P ,LT_K ,IADK ,JDIK )
2031C D_Y = D_Y + LT_K0 * D_P
2032 CALL MIC_SPMV(NDDL ,Y ,P ,LT_K0,IADK0,JDIK0)
2033 IF(NDDLI > 0) THEN
2034C D_Y = D_Y + LT_I0 * D_P
2035ctempo CALL MIC_SPMV(NDDL ,Y ,P ,LT_I0,IADI0,JDII0)
2036 END IF
2037c!$OMP END PARALLEL
2038c IF(NSPMD > 1)THEN
2039c IF(IMONM > 0) CALL STARTIME(TIMERS,68)
2040c!dir$ omp offload target(mic:MICID)
2041c & nocopy(y,ifr2k)
2042c & out(vgat:length(lcom))
2043c!$OMP PARALLEL default(shared)
2044c CALL MIC_GATHER(NDDL, LCOM, Y, VGAT, IFR2K)
2045c!$OMP END PARALLEL
2046c IF(IMONM > 0) CALL STOPTIME(TIMERS,68)
2047c IF(IMONM > 0) CALL STARTIME(TIMERS,66)
2048c CALL SPMD_SUMFC_V(VGAT,VSCA,INDEX,LCOM)
2049c IF(IMONM > 0) CALL STOPTIME(TIMERS,66)
2050c IF(IMONM > 0) CALL STARTIME(TIMERS,69)
2051c!dir$ omp offload target(mic:MICID)
2052c & nocopy(y,ifr2k,index)
2053c & in(vsca:length(lcom))
2054c!$OMP PARALLEL default(shared)
2055c CALL MIC_SCATTER(NDDL, LCOM, Y, VSCA, IFR2K, INDEX)
2056c!$OMP END PARALLEL
2057c IF(IMONM > 0) CALL STOPTIME(TIMERS,69)
2058c END IF
2059C---------------------
2060c IF(NSPMD == 1) THEN
2061c!dir$ omp offload target(mic:MICID)
2062c & nocopy(p,y)
2063c & inout(S)
2064c!$OMP PARALLEL default(shared) shared(S)
2065 CALL MIC_DDOT(NDDL, P, Y, S)
2066c!$OMP END PARALLEL
2067c ELSE
2068C si spmd penser a faire avant D_Y = D_Y*W pour les poids
2069c!dir$ omp offload target(mic:MICID)
2070c & nocopy(y,w,p,z)
2071c & inout(S)
2072c!$OMP PARALLEL default(shared) shared(S)
2073c CALL MIC_MV(NDDL,Y,W,Z)
2074c CALL MIC_DDOT(NDDL, P, Z, S)
2075c!$OMP END PARALLEL
2076c IF(IMONM > 0) CALL STARTIME(TIMERS,67)
2077c CALL SPMD_SUM_S(S)
2078c IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
2079c END IF
2080C---------------------
2081 ALPHA=G0/S
2082c!dir$ omp offload target(mic:MICID)
2083c & nocopy(iadk,iadk0,iadm,iadm0,lt_k,lt_k0,lt_m)
2084c & nocopy(lt_m0,jdik,jdik0,jdim,jdim0,diag_k,diag_m,y,p,x,r,z)
2085c!$OMP PARALLEL default(shared)
2086 CALL MIC_DAXPY(NDDL, ALPHA, P, X)
2087 CALL MIC_DAXPY(NDDL,-ALPHA, Y, R)
2088C----------------------
2089C CALL PREC_SOLVGH(IPREC , ITASK ,NDDL ,IADM , JDIM ,
2090C 6 DIAG_M, LT_M ,R ,Z ,
2091C 7 F_DDL ,L_DDL )
2092 IF(IPREC == 1)THEN
2093C D_Z = D_R
2094 CALL MIC_DCOPY(NDDL, R, Z)
2095 ELSEIF(IPREC == 2)THEN
2096C D_Z = D_R * DIAG_M
2097 CALL MIC_MV(NDDL,R,DIAG_M,Z)
2098 ELSEIF(IPREC == 5)THEN
2099C D_Y = D_R
2100 CALL MIC_DCOPY(NDDL, R, Y)
2101C D_Y = D_Y + LT_M * D_R
2102 CALL MIC_SPMV(NDDL ,Y ,R ,LT_M ,IADM ,JDIM )
2103C D_Z = D_Y * DIAG_M
2104 CALL MIC_MV(NDDL,Y,DIAG_M,Z)
2105 CALL MIC_DCOPY(NDDL, Z, Y)
2106C D_Z = D_Z + LT_M0 * D_Y
2107 CALL MIC_SPMV(NDDL ,Z ,Y ,LT_M0,IADM0,JDIM0)
2108 ELSE
2109 END IF
2110c!$OMP END PARALLEL
2111C----------------------
2112c IF(NSPMD > 1 .AND. IPREC > 1)THEN
2113c IF(IMONM > 0) CALL STARTIME(TIMERS,68)
2114c!dir$ omp offload target(mic:MICID)
2115c & nocopy(z,ifr2k)
2116c & out(vgat:length(lcom))
2117c!$OMP PARALLEL default(shared)
2118c CALL MIC_GATHER(NDDL, LCOM, Z, VGAT, IFR2K)
2119c!$OMP END PARALLEL
2120c IF(IMONM > 0) CALL STOPTIME(TIMERS,68)
2121c IF(IMONM > 0) CALL STARTIME(TIMERS,66)
2122c CALL SPMD_SUMFC_V(VGAT,VSCA,INDEX,LCOM)
2123c IF(IMONM > 0) CALL STOPTIME(TIMERS,66)
2124c IF(IMONM > 0) CALL STARTIME(TIMERS,69)
2125c!dir$ omp offload target(mic:MICID)
2126c & nocopy(z,ifr2k,index)
2127c & in(vsca:length(lcom))
2128c!$OMP PARALLEL default(shared)
2129c CALL MIC_SCATTER(NDDL, LCOM, Z, VSCA, IFR2K, INDEX)
2130c!$OMP END PARALLEL
2131c IF(IMONM > 0) CALL STOPTIME(TIMERS,69)
2132c END IF
2133C---------------------
2134c IF(NSPMD == 1) THEN
2135c!dir$ omp offload target(mic:MICID)
2136c & nocopy(r,z)
2137c & inout(G1)
2138c!$OMP PARALLEL default(shared) shared(G1)
2139 CALL MIC_DDOT(NDDL, R, Z, G1)
2140c!$OMP END PARALLEL
2141c ELSE
2142C si spmd penser a faire avant D_Y = D_Z*W pour les poids
2143c!dir$ omp offload target(mic:MICID)
2144c & nocopy(z,w,r,y)
2145c & inout(G1)
2146c!$OMP PARALLEL default(shared) shared(G1)
2147c CALL MIC_MV(NDDL,Z,W,Y)
2148c CALL MIC_DDOT(NDDL, R, Y, G1)
2149c!$OMP END PARALLEL
2150c IF(IMONM > 0) CALL STARTIME(TIMERS,67)
2151c CALL SPMD_SUM_S(G1)
2152c IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
2153c END IF
2154C---------------------
2155c+1
2156!$OMP single
2157 BETA=G1/G0
2158 G0 = G1
2159c+1
2160!$OMP end single
2161 IF (ITOL==1) THEN
2162c IF(NSPMD == 1) THEN
2163c!dir$ omp offload target(mic:MICID)
2164c & nocopy(r)
2165c & inout(R2)
2166c!$OMP PARALLEL default(shared) shared(R2)
2167 CALL MIC_DDOT(NDDL, R, R, R2)
2168c!$OMP END PARALLEL
2169c ELSE
2170C si spmd penser a faire avant D_Y = D_R*W pour les poids
2171c!dir$ omp offload target(mic:MICID)
2172c & nocopy(r,w,y)
2173c & inout(R2)
2174c!$OMP PARALLEL default(shared) shared(R2)
2175c CALL MIC_MV(NDDL,R,W,Y)
2176c CALL MIC_DDOT(NDDL, Y, Y, R2)
2177c!$OMP END PARALLEL
2178c IF(IMONM > 0) CALL STARTIME(TIMERS,67)
2179c CALL SPMD_SUM_S(R2)
2180c IF(IMONM > 0) CALL STOPTIME(TIMERS,67)
2181c END IF
2182 ELSEIF (ITOL==3) THEN
2183c+1
2184!$OMP single
2185 R2 =ABS(G1)
2186 S=ONE/ALPHA
2187 L_A=S+B_OLD*A_OLD
2188C----------L_B2 : beta(j-1)^2-A_OLD :1/alpha(j-1)-------
2189 L_B2=ABS(BETA)*S*S
2190 L_B=SQRT(L_B2)
2191 A_OLD=S
2192 B_OLD=BETA
2193 ANORM2=TNORM2
2194 TNORM2=TNORM2+L_A*L_A+OLDB2+L_B2
2195 GAMMA = SQRT( GBAR*GBAR + OLDB2 )
2196 CS = GBAR / GAMMA
2197 SN = OLDB / GAMMA
2198 DELTA = CS * DBAR + SN * L_A
2199 GBAR = SN * DBAR - CS * L_A
2200 EPSLN = SN * L_B
2201 DBAR = - CS * L_B
2202 ZL = RHS1 / GAMMA
2203 XNORM2 = XNORM2+ZL*ZL
2204 GMAX = MAX( GMAX, GAMMA )
2205 GMIN = MIN( GMIN, GAMMA )
2206 RHS1 = RHS2 - DELTA * ZL
2207 RHS2 = - EPSLN * ZL
2208 TOLN=TOLS2*ANORM2*XNORM2
2209 OLDB2 = L_B2
2210 OLDB = L_B
2211c+1
2212!$OMP end single
2213 ELSEIF (ITOL==4) THEN
2214 TMP=ALPHA*ALPHA*ABS(G1)
2215 IF (IT>=ND) THEN
2216 DO J=1,ND-1
2217 EPS(J)=EPS(J+1)
2218 ENDDO
2219 EPS(ND)=TMP
2220 R2=ZERO
2221 DO J=1,ND
2222 R2=R2+EPS(J)
2223 ENDDO
2224 ELSE
2225 EPS(IT+1)=TMP
2226 R2=R2+TMP
2227 ENDIF
2228 ELSE
2229 R2 =ABS(G1)
2230 ENDIF
2231c+1
2232!$OMP single
2233.AND. IF(IPRI/=0ITASK==0)THEN
2234 IF(MOD(IT,IP)==0)THEN
2235 RR = SQRT(R2/R02)
2236 WRITE(IOUT,1001)IT,RR
2237 IF(IPRI<0) WRITE(ISTDO,1001)IT,RR
2238 ENDIF
2239 ENDIF
2240C----------------------
2241c CALL MY_BARRIER
2242C---------------------
2243C ISTOP=CRIT_STOP(IT,R2,NLIM,TOLN)
2244 IF (IT>=NLIM) THEN
2245 ISTOP = 0
2246 ELSE
2247 IF (R2<=TOLN) THEN
2248 ISTOP = 0
2249 ELSE
2250 ISTOP = 1
2251 ENDIF
2252 ENDIF
2253C----------------------
2254c CALL MY_BARRIER
2255C---------------------
2256 IT = IT +1
2257c+1
2258!$OMP end single
2259 ENDDO
2260 2060 CONTINUE
2261!$OMP single
2262 tt=OMP_GET_WTIME()-tt
2263 i= omp_get_num_threads()
2264 print *,' '
2265 print *,' '
2266 print *,' execution time on mic with',i,' threads'
2267 print *,' ==> ',tt,' seconds'
2268 print *,' '
2269!$OMP END single
2270c+1
2271!$OMP END PARALLEL
2272 WRITE(IOUT,1001)IT-1,RR
2273 END IF ! fin NSPMD ==1
2274 TF=OMP_GET_WTIME()
2275!dec$ omp offload target(mic:MICID) out(X,R:alloc_if(.false.))
2276!$OMP PARALLEL
2277c!$OMP single
2278c tt=OMP_GET_WTIME()-tt
2279c i= omp_get_num_threads()
2280c print *,' '
2281c print *,' '
2282c print *,' Execution time on MIC with',i,' THREADS'
2283c print *,' ==> ',tt,' SECONDS'
2284c print *,' '
2285c!$OMP END single
2286!$OMP END PARALLEL
2287 TF=OMP_GET_WTIME()-TF
2288 print *,'transfer time mic=>cpu(s) =',TF
2289 print *,'transfer rate mic=>cpu(mb/s)=',
2290 . ((2*NDDL/(1024*1024))*8)/TF
2291c END IF ! fin NSPMD ==1
2292C fin code provisoire
2293
2294c IF(NSPMD > 1) THEN
2295c DEALLOCATE(W)
2296c DEALLOCATE(VGAT)
2297c DEALLOCATE(VSCA)
2298c DEALLOCATE(INDEX)
2299c END IF
2300
2301 END IF ! fin GPUR4R8==2 (MIC double precision)
2302
2303 END IF !only task0
2304C fin code specifique MIC double precision
2305#endif
2306C fin code specifique GPU
2307C routine speciale pour sauvegarder les resultats
2308 IF(ISAVE==1)CALL IMP_RSAVE(NDDL,X,R)
2309
2310C
2311C------due to ISTOP local variable
2312C IF(IT>=NLIM.AND.ISOLV==7) ISTOP = 3
2313 IF(ITASK == 0) THEN
2314#include "lockon.inc"
2315 ISTOP_H = ISTOP
2316#include "lockoff.inc"
2317 CALL PUPD(N_MAX,NDDLI_G,IT)
2318C
2319 IF(IT>=NLIM)THEN
2320 IF (ISOLV==7) THEN
2321#include "lockon.inc"
2322 ISTOP_H = 3
2323#include "lockoff.inc"
2324 ELSE
2325#include "lockon.inc"
2326 ISTOP_H = 1
2327#include "lockoff.inc"
2328
2329.NOT. IF (MIXEDSOL()) THEN
2330 WRITE(IOUT,*)
2331 WRITE(IOUT,1003)NLIM
2332 WRITE(IOUT,*)
2333 WRITE(ISTDO,*)
2334 WRITE(ISTDO,1003)NLIM
2335 WRITE(ISTDO,*)
2336 IF(ITOL==3)THEN
2337 WRITE(IOUT,2000)
2338 . EXIT, SQRT(ANORM2), SQRT(XNORM2),GMAX/GMIN
2339 WRITE(ISTDO,2000)
2340 . EXIT, SQRT(ANORM2), SQRT(XNORM2),GMAX/GMIN
2341 ENDIF
2342C----add condition to remove wrong message with mix solver
2343 IF (R2>HUNDRED*TOLN) THEN
2344 ISTOP_H = 2
2345 RR = SQRT(R2/R02)
2346 WRITE(IOUT,*)
2347 WRITE(IOUT,1004)RR
2348 WRITE(IOUT,*)
2349 WRITE(ISTDO,*)
2350 WRITE(ISTDO,1004)RR
2351 WRITE(ISTDO,*)
2352 ENDIF
2353.NOT. END IF !((MIXEDSOL)) THEN
2354
2355 END IF !(ISOLV==7) THEN
2356 ENDIF
2357 IF(IPRI/=0)THEN
2358 RR = SQRT(R2/R02)
2359 WRITE(IOUT,*)
2360 WRITE(IOUT,1002)IT,RR
2361 IF(ITOL==3)WRITE(IOUT,2000)
2362 . EXIT, SQRT(ANORM2), SQRT(XNORM2),GMAX/GMIN
2363 WRITE(IOUT,*)
2364 IF(IPRI<0) THEN
2365 WRITE(ISTDO,*)
2366 WRITE(ISTDO,1002)IT,RR
2367 IF(ITOL==3)WRITE(ISTDO,2000)
2368 . EXIT, SQRT(ANORM2), SQRT(XNORM2),GMAX/GMIN
2369 WRITE(ISTDO,*)
2370 ENDIF
2371 ENDIF
2372 IF(ITOL==3) THEN
2373 K_LAMDA0= GMIN
2374 K_LAMDA1= GMAX
2375 ELSE
2376 K_LAMDA0= ZERO
2377 K_LAMDA1= ZERO
2378 ENDIF
2379 END IF !(ITASK == 0) THEN
2380C----------------------
2381 CALL MY_BARRIER
2382C----------------------
2383#include "lockon.inc"
2384 ISTOP = ISTOP_H
2385#include "lockoff.inc"
2386
2387 CALL MY_BARRIER
2388C--------------------------------------------
2389 1000 FORMAT(5X,'iteration=',I8,5X,' initial residual norm=',E11.4)
2390 1001 FORMAT(5X,'iteration=',I8,5X,' relative residual norm=',E11.4)
2391 1002 FORMAT(3X,'total c.g. iteration=',I8,5X,
2392 . ' relative residual norm=',E11.4)
2393 1003 FORMAT(5X,
2394 . '---warning : the iteration limit number was reached',I8)
2395 1004 FORMAT(5X,'warning:c.g stop with relative residual norm=',E11.4)
2396 2000 FORMAT(/ 5X, A, 2X, 'anorm =', E12.4, 2X, 'xnorm =', E12.4,2X,
2397 . 'kcond =', E12.4)
2398 2002 FORMAT(/ 5X, 'with', 2X, 'alfa =', E12.4, 2X, 'beta =',
2399 . E12.4,2X,'oldb =', E12.4)
2400 RETURN
2401#endif
2402 END
2403
2404!||====================================================================
2405!|| pupd ../engine/source/implicit/imp_pcg.F
2406!||--- called by ------------------------------------------------------
2407!|| imp_pcgh ../engine/source/implicit/imp_pcg.F
2408!|| imp_ppcgh ../engine/source/implicit/imp_pcg.F
2409!||====================================================================
2410 SUBROUTINE PUPD(NDDL,NDDLI,ITN)
2411C-----------------------------------------------
2412C I m p l i c i t T y p e s
2413C-----------------------------------------------
2414#include "implicit_f.inc"
2415C-----------------------------------------------
2416C C o m m o n B l o c k s
2417C-----------------------------------------------
2418#include "impl1_c.inc"
2419C-----------------------------------------------
2420C D u m m y A r g u m e n t s
2421C-----------------------------------------------
2422 INTEGER NDDL,NDDLI,ITN
2423C-----------------------------------------------
2424C L o c a l V a r i a b l e s
2425C-----------------------------------------------
2426 INTEGER I,NLIM
2427C-----------------------------
2428C
2429 IT_PCG = IT_PCG + ITN
2430.OR. IF (ISOLV==5ISOLV==6) THEN
2431 IF (IPUPD>0) THEN
2432 NLIM=IPUPD+NDDLI/10
2433 ELSE
2434 NLIM=MAX(5,NDDL/10000)+NDDLI/10
2435 ENDIF
2436 IF (ITN>NLIM) THEN
2437c WRITE(IOUT,*) '****RE-FRACTO. DUE TO ITN=',ITN
2438 IT_PCG = -IT_PCG
2439 ENDIF
2440 END IF
2441C--------------------------------------------
2442 RETURN
2443 END
2444!||====================================================================
2445!|| nlim_mix ../engine/source/implicit/imp_pcg.F
2446!||--- called by ------------------------------------------------------
2447!|| imp_pcgh ../engine/source/implicit/imp_pcg.F
2448!|| imp_ppcgh ../engine/source/implicit/imp_pcg.F
2449!||====================================================================
2450 SUBROUTINE NLIM_MIX(NDDL,NDDLI,NLIM)
2451C-----------------------------------------------
2452C I m p l i c i t T y p e s
2453C-----------------------------------------------
2454#include "implicit_f.inc"
2455C-----------------------------------------------
2456C C o m m o n B l o c k s
2457C-----------------------------------------------
2458#include "impl1_c.inc"
2459C-----------------------------------------------
2460C D u m m y A r g u m e n t s
2461C-----------------------------------------------
2462 INTEGER NDDL,NDDLI,NLIM
2463C-----------------------------------------------
2464C L o c a l V a r i a b l e s
2465C-----------------------------------------------
2466 INTEGER I
2467C-----------------------------
2468.AND. IF (ISOLV/=5ISOLV/=6) RETURN
2469 IF (IPUPD>0) THEN
2470 NLIM=10*IPUPD+NDDLI
2471 ELSE
2472 NLIM=10+NDDL/100+NDDLI
2473 ENDIF
2474 NLIM=MIN(NLIM,NDDL)
2475C--------------------------------------------
2476 RETURN
2477 END
2478!||====================================================================
2479!|| mixedsol ../engine/source/implicit/imp_pcg.F
2480!||--- called by ------------------------------------------------------
2481!|| imp_pcgh ../engine/source/implicit/imp_pcg.F
2482!||====================================================================
2483 LOGICAL FUNCTION MIXEDSOL()
2484C-----------------------------------------------
2485C I m p l i c i t T y p e s
2486C-----------------------------------------------
2487#include "implicit_f.inc"
2488C-----------------------------------------------
2489C C o m m o n B l o c k s
2490C-----------------------------------------------
2491#include "impl1_c.inc"
2492C-----------------------------------------------
2493C L o c a l V a r i a b l e s
2494C-----------------------------------------------
2495 INTEGER I,J
2496C-----------------------------------------------
2497 MIXEDSOL=.FALSE.
2498.OR..AND. IF ((ISOLV==5ISOLV==6)IT_PCG<0) MIXEDSOL=.TRUE.
2499
2500 RETURN
2501 END
2502
2503!||====================================================================
2504!|| imp_inix ../engine/source/implicit/imp_pcg.F
2505!||--- called by ------------------------------------------------------
2506!|| imp_ppcgh ../engine/source/implicit/imp_pcg.F
2507!||--- calls -----------------------------------------------------
2508!|| mav_mn ../engine/source/implicit/produt_v.F
2509!|| mav_nm ../engine/source/implicit/produt_v.F
2510!|| my_barrier ../engine/source/system/machine.F
2511!||--- uses -----------------------------------------------------
2512!|| imp_pcg_proj ../engine/share/modules/impbufdef_mod.F
2513!||====================================================================
2514 SUBROUTINE IMP_INIX(F_DDL,L_DDL,NDDL,X ,R ,W_DDL,ITASK )
2515C-----------------------------------------------
2516C M o d u l e s
2517C-----------------------------------------------
2518 USE IMP_PCG_PROJ
2519C-----------------------------------------------
2520C I m p l i c i t T y p e s
2521C-----------------------------------------------
2522#include "implicit_f.inc"
2523C-----------------------------------------------
2524C C o m m o n B l o c k s
2525C-----------------------------------------------
2526#include "impl1_c.inc"
2527C-----------------------------------------------
2528C D u m m y A r g u m e n t s
2529C-----------------------------------------------
2530 INTEGER F_DDL,L_DDL,NDDL,ITASK,W_DDL(*)
2531 my_real
2532 . X(*), R(*)
2533C-----------------------------------------------
2534c FUNCTION: initialization of X0=[S][LAMDA]^-1[S]^t{R}
2535c
2536c Note:
2537c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
2538c
2539c TYPE NAME FUNCTION
2540c I NDDL - equation dimension
2541c I NPV - projection vector number
2542c I W_DDL(*) - itag for each id of subdomains
2543c I F_ND,L_ND,ITASK - id in each ITASK:thread id (//)
2544c I R(NDDL) - right-hand vector
2545c O X(NDDL) - X0
2546C-----------------------------------------------
2547C L o c a l V a r i a b l e s
2548C-----------------------------------------------
2549 INTEGER I,J,M,NPV
2550C----------------------
2551 IF (NPCGPV == 0) RETURN
2552 NPV = NPCGPV
2553C----------------------
2554C {W}=[S]^t{R}
2555C---------------------
2556 CALL MAV_NM(F_DDL,L_DDL,NDDL ,NPV ,PROJ_S ,
2557 . R ,PROJ_W,W_DDL,ITASK )
2558C----------------------
2559 CALL MY_BARRIER
2560C---------------------
2561C [LAMDA]^-1{W}
2562C---------------------
2563 IF (ITASK ==0) THEN
2564 DO I=1,NPV
2565 PROJ_W(I)=PROJ_W(I)*PROJ_LA_1(I)
2566 ENDDO
2567C---------------------
2568C {X0}=[S]{W}
2569C---------------------
2570 CALL MAV_MN(NDDL ,NPV ,PROJ_S ,PROJ_W ,X ,ITASK )
2571C----------------------
2572 END IF !(ITASK ==0) THEN
2573C----------------------
2574 RETURN
2575 END
2576!||====================================================================
2577!|| imp_pro_p ../engine/source/implicit/imp_pcg.F
2578!||--- called by ------------------------------------------------------
2579!|| imp_ppcgh ../engine/source/implicit/imp_pcg.F
2580!||--- calls -----------------------------------------------------
2581!|| mav_mn ../engine/source/implicit/produt_v.F
2582!|| mav_nm ../engine/source/implicit/produt_v.F
2583!|| my_barrier ../engine/source/system/machine.F
2584!||--- uses -----------------------------------------------------
2585!|| imp_pcg_proj ../engine/share/modules/impbufdef_mod.F
2586!||====================================================================
2587 SUBROUTINE IMP_PRO_P(F_DDL,L_DDL,NDDL ,P ,W_DDL,ITASK )
2588C-----------------------------------------------
2589C M o d u l e s
2590C-----------------------------------------------
2591 USE IMP_PCG_PROJ
2592C-----------------------------------------------
2593C I m p l i c i t T y p e s
2594C-----------------------------------------------
2595#include "implicit_f.inc"
2596C-----------------------------------------------
2597C C o m m o n B l o c k s
2598C-----------------------------------------------
2599#include "impl1_c.inc"
2600C-----------------------------------------------
2601C D u m m y A r g u m e n t s
2602C-----------------------------------------------
2603 INTEGER F_DDL ,L_DDL ,NDDL,ITASK,W_DDL(*)
2604 my_real
2605 . P(*)
2606C-----------------------------------------------
2607c FUNCTION: Projection of {p}={p}-[S][LAMDA]^-1[T]^t{p}
2608c
2609c Note:
2610c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
2611c
2612c TYPE NAME FUNCTION
2613c I NDDL - equation dimension
2614c I NPV - projection vector number
2615c I W_DDL(*) - itag for each id of subdomains
2616c I F_ND,L_ND,ITASK - id in each ITASK:thread id (//)
2617c IO P(NDDL) - PCG p vector
2618C-----------------------------------------------
2619C L o c a l V a r i a b l e s
2620C-----------------------------------------------
2621 INTEGER I,J,NPV
2622 my_real
2623 . P1(NDDL)
2624C----------------------
2625 IF (NPCGPV == 0) RETURN
2626 NPV = NPCGPV
2627C----------------------
2628C {W}=[T]^t{p}
2629C---------------------
2630 CALL MAV_NM(F_DDL,L_DDL,NDDL ,NPV ,PROJ_T ,
2631 . P ,PROJ_W ,W_DDL,ITASK )
2632C----------------------
2633 CALL MY_BARRIER
2634C---------------------
2635C [LAMDA]^-1{W}
2636C---------------------
2637 IF (ITASK==0) THEN
2638 DO I=1,NPV
2639 PROJ_W(I)=PROJ_W(I)*PROJ_LA_1(I)
2640 ENDDO
2641C----------------------
2642C {P1}=[S]{W}
2643C---------------------
2644 CALL MAV_MN(NDDL ,NPV ,PROJ_S ,PROJ_W ,P1 ,ITASK )
2645 DO I=1,NDDL
2646 P(I)=P(I)-P1(I)
2647 ENDDO
2648 END IF
2649C
2650 RETURN
2651 END
2652!||====================================================================
2653!|| imp_inis ../engine/source/implicit/imp_pcg.F
2654!||--- called by ------------------------------------------------------
2655!|| imp_inisi ../engine/source/implicit/imp_pcg.F
2656!||====================================================================
2657 SUBROUTINE IMP_INIS(ND ,F_ND ,L_ND, NPF, NPL, S )
2658C-----------------------------------------------
2659C I m p l i c i t T y p e s
2660C-----------------------------------------------
2661#include "implicit_f.inc"
2662C-----------------------------------------------
2663C D u m m y A r g u m e n t s
2664C-----------------------------------------------
2665 INTEGER ND,F_ND ,L_ND,NPF, NPL
2666 my_real
2667 . S(ND,*),ALEAT
2668 EXTERNAL ALEAT
2669C-----------------------------------------------
2670c FUNCTION: initialization of [S]
2671c
2672c Note:
2673c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
2674c
2675c TYPE NAME FUNCTION
2676c I F_ND ,L_ND - equation dimension divised by Nthead
2677c I NPF,NPL - projection vector number (first to last)
2678c I ND - equation dimension
2679c I S(ND,NPV) - S-Matrix
2680C-----------------------------------------------
2681C L o c a l V a r i a b l e s
2682C-----------------------------------------------
2683 INTEGER I,J,M,INORM
2684C---------------------
2685 DO J=NPF,NPL
2686 DO I=F_ND ,L_ND
2687 S(I,J)=ALEAT()
2688 ENDDO
2689 ENDDO
2690C
2691 RETURN
2692 END
2693!||====================================================================
2694!|| imp_updv2 ../engine/source/implicit/imp_pcg.F
2695!||--- called by ------------------------------------------------------
2696!|| imp_ppcgh ../engine/source/implicit/imp_pcg.F
2697!||--- calls -----------------------------------------------------
2698!|| arret ../engine/source/system/arret.F
2699!|| mmav_lth ../engine/source/implicit/produt_v.F
2700!|| produt_h ../engine/source/implicit/produt_v.F
2701!||--- uses -----------------------------------------------------
2702!|| groupdef_mod ../common_source/modules/groupdef_mod.F
2703!|| imp_pcg_proj ../engine/share/modules/impbufdef_mod.F
2704!|| intbufdef_mod ../common_source/modules/interfaces/intbufdef_mod.F90
2705!||====================================================================
2706 SUBROUTINE IMP_UPDV2(
2707 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K ,
2708 2 LT_K ,ITOK ,IADI ,JDII ,LT_I ,
2709 3 ITASK ,W_DDL ,A ,AR ,VE ,
2710 4 D ,DR ,NDOF ,IPARI ,INTBUF_TAB,
2711 5 NUM_IMP,NS_IMP,NE_IMP,NSREM ,
2712 6 NSL ,NMONV ,IMONV ,MONVOL,IGRSURF ,
2713 7 VOLMON,FR_MV ,IBFV ,SKEW ,
2714 8 XFRAME ,IND_IMP,XI_C ,MS ,XE ,
2715 9 IRBE3 ,LRBE3 ,IRBE2 ,LRBE2 ,U_TMP ,
2716 A P ,Y ,F_DDL ,L_DDL ,IT ,
2717 B ICHECK,Z ,IADM ,JDIM ,DIAG_M,LT_M)
2718C-----------------------------------------------
2719C M o d u l e s
2720C-----------------------------------------------
2721 USE IMP_PCG_PROJ
2722 USE INTBUFDEF_MOD
2723 USE GROUPDEF_MOD
2724C-----------------------------------------------
2725C I m p l i c i t T y p e s
2726C-----------------------------------------------
2727#include "implicit_f.inc"
2728C-----------------------------------------------
2729C C o m m o n B l o c k s
2730C-----------------------------------------------
2731#include "com04_c.inc"
2732#include "impl1_c.inc"
2733#include "units_c.inc"
2734C-----------------------------------------------
2735C D u m m y A r g u m e n t s
2736C-----------------------------------------------
2737 INTEGER F_DDL,L_DDL,NDDL ,IADK(*) ,JDIK(*),
2738 . NDDLI ,ITOK(*) ,IADI(*),JDII(*),
2739 . ITASK,W_DDL(*),IBFV(*),IRBE3(*) ,LRBE3(*),
2740 . IRBE2(*) ,LRBE2(*),IT,ICHECK,IADM(*),JDIM(*)
2741 INTEGER NDOF(*),NE_IMP(*),NSREM ,NSL,
2742 . IPARI(*) ,NUM_IMP(*),NS_IMP(*) ,IND_IMP(*)
2743 INTEGER NMONV,IMONV(*),MONVOL(*),FR_MV(*)
2744 my_real
2745 . DIAG_K(*), LT_K(*),LT_I(*) ,XI_C(*),U_TMP(*),P(*),Y(*)
2746 my_real
2747 . A(3,*),AR(3,*),VE(3,*),D(3,*),DR(3,*),XE(3,*),
2748 . MS(*),VOLMON(*),SKEW(*) ,XFRAME(*),
2749 . Z(*),DIAG_M(*),LT_M(*)
2750
2751 TYPE(INTBUF_STRUCT_) INTBUF_TAB(*)
2752 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
2753C-----------------------------------------------
2754c FUNCTION: update S,T of Projection
2755c
2756c Note:
2757c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
2758c
2759c TYPE NAME FUNCTION
2760c I NDDL,NNZ - dimension of [K] and number of non zero (strict triangular)
2761c I IADK,JDIK - indice arrays for compressed row(col.) format of [K]
2762c I DIAG_K(NDDL) - diagonal terms of [K]
2763c I LT_K(NNZ) - strict triangular of [K]
2764c O Proj_S(NDDL,M_VS) - [S] reduced small Eigenvectors
2765c O Proj_T(NDDL,M_VS) - [T] =[K][S]
2766C-----------------------------------------------
2767C L o c a l V a r i a b l e s
2768C-----------------------------------------------
2769#if defined(MYREAL8) && !defined(WITHOUT_LINALG)
2770 CHARACTER JOBZ, UPLO
2771 INTEGER I,J,IP,NLIM,ND,IUPD,IPRI,IERROR,NNZI,M,
2772 . F_DDLI,L_DDLI,INFO,LW,M_VS1,INORM,NP
2773 my_real
2774 . WORK(3*M_VS+3),W(M_VS+1),P2,PY,S,VP,VY,P_K2(2,2),
2775 . WORK_T(NDDL)
2776C--- COMPUTE another lowest eigenvecor V(U_TMP)
2777 IF (IPRO_S0/=4) RETURN
2778C
2779 CALL PRODUT_H( F_DDL ,L_DDL ,P ,P ,W_DDL , P2 ,ITASK )
2780 CALL PRODUT_H( F_DDL ,L_DDL ,P ,Y ,W_DDL , PY ,ITASK )
2781 IF (IT==0 ) THEN
2782 M_TMP(1,1)=ONE
2783 S = ONE/SQRT(P2)
2784 K_TMP(1,1)=PY/P2
2785 DO I=F_DDL,L_DDL
2786 U_TMP(I)=S*P(I)
2787 ENDDO
2788 ELSE
2789 CALL PRODUT_H( F_DDL ,L_DDL ,P ,U_TMP,W_DDL , VP ,ITASK )
2790 CALL PRODUT_H( F_DDL ,L_DDL ,Y ,U_TMP,W_DDL , VY ,ITASK )
2791 M_TMP(1,2)=VP
2792 M_TMP(2,2)=P2
2793 K_TMP(1,2)=VY
2794 K_TMP(2,2)=PY
2795 ND =2
2796 JOBZ='v'
2797 UPLO='u'
2798 LW=3*ND
2799 CALL DSYGV( 1, JOBZ, UPLO, ND, K_TMP, ND, M_TMP, ND, W, WORK,
2800 $ LW, INFO )
2801 DO I=F_DDL,L_DDL
2802 U_TMP(I)=K_TMP(1,1)*U_TMP(I)+K_TMP(2,1)*P(I)
2803 ENDDO
2804 M_TMP(1,1)=ONE
2805 K_TMP(1,1)=W(1)
2806 ENDIF
2807C-------checking V(U_TMP)
2808.OR. if (icheck==1IMPDEB>0) then
2809 CALL PRODUT_H( F_DDL ,L_DDL ,U_TMP,U_TMP,W_DDL , S ,ITASK )
2810 CALL MMAV_LTH(
2811 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K,
2812 2 LT_K ,IADI ,JDII ,ITOK ,LT_I ,
2813 3 U_TMP ,WORK_T,A ,AR ,
2814 5 VE ,MS ,XE ,D ,DR ,
2815 6 NDOF ,IPARI ,INTBUF_TAB ,NUM_IMP,
2816 7 NS_IMP,NE_IMP,NSREM ,NSL ,IBFV ,
2817 8 SKEW ,XFRAME,MONVOL,VOLMON,IGRSURF ,
2818 9 FR_MV,NMONV ,IMONV ,IND_IMP,
2819 A XI_C ,IUPD ,IRBE3 ,LRBE3 ,IRBE2 ,
2820 B LRBE2 ,IADM ,JDIM ,DIAG_M,LT_M ,
2821 C F_DDL ,L_DDL ,ITASK ,Z )
2822 CALL PRODUT_H( F_DDL ,L_DDL ,U_TMP,WORK_T,W_DDL ,P2,ITASK)
2823 write(iout,*)'inf,|v|,vav,lamda1=',info,s,p2,K_TMP(1,1)
2824 endif
2825
2826#else
2827 CALL ARRET(5)
2828#endif
2829C
2830 RETURN
2831 END
2832!||====================================================================
2833!|| imp_updv1 ../engine/source/implicit/imp_pcg.F
2834!||--- called by ------------------------------------------------------
2835!|| imp_ppcgh ../engine/source/implicit/imp_pcg.F
2836!||====================================================================
2837 SUBROUTINE IMP_UPDV1(
2838 1 R2 ,R2_OLD ,RMAX ,PROJ_U,P ,
2839 2 IT ,F_DDL,L_DDL )
2840C-----------------------------------------------
2841C I m p l i c i t T y p e s
2842C-----------------------------------------------
2843#include "implicit_f.inc"
2844C-----------------------------------------------
2845C C o m m o n B l o c k s
2846C-----------------------------------------------
2847#include "impl1_c.inc"
2848C-----------------------------------------------
2849C D u m m y A r g u m e n t s
2850C-----------------------------------------------
2851 INTEGER IT,F_DDL,L_DDL
2852 my_real
2853 . R2, R2_OLD ,RMAX ,PROJ_U(*),P(*)
2854C-----------------------------------------------
2855C L o c a l V a r i a b l e s
2856C-----------------------------------------------
2857 INTEGER I
2858 my_real
2859 . RR
2860C-----------------------------------
2861 IF (IPRO_S0==4) RETURN
2862C
2863 IF (IT==1) THEN
2864 R2_OLD=R2
2865 ELSE
2866 RR=R2/R2_OLD
2867 R2_OLD=R2
2868 IF (RMAX<RR) THEN
2869 DO I=F_DDL,L_DDL
2870 PROJ_U(I) = P(I)
2871 ENDDO
2872 RMAX=RR
2873 ENDIF
2874 END IF
2875C
2876 RETURN
2877 END
2878!||====================================================================
2879!|| imp_ppcgh ../engine/source/implicit/imp_pcg.F
2880!||--- called by ------------------------------------------------------
2881!|| lin_solvh0 ../engine/source/implicit/lin_solv.F
2882!|| lin_solvh1 ../engine/source/implicit/lin_solv.F
2883!|| lin_solvih2 ../engine/source/implicit/lin_solv.F
2884!||--- calls -----------------------------------------------------
2885!|| crit_stop ../engine/source/implicit/imp_pcg.F
2886!|| imp_inisi ../engine/source/implicit/imp_pcg.F
2887!|| imp_inist ../engine/source/implicit/imp_pcg.F
2888!|| imp_inix ../engine/source/implicit/imp_pcg.F
2889!|| imp_pro_p ../engine/source/implicit/imp_pcg.F
2890!|| imp_updst ../engine/source/implicit/imp_pcg.F
2891!|| imp_updv1 ../engine/source/implicit/imp_pcg.F
2892!|| imp_updv2 ../engine/source/implicit/imp_pcg.F
2893!|| mmav_lth ../engine/source/implicit/produt_v.F
2894!|| mmv_lh ../engine/source/implicit/produt_v.F
2895!|| mmv_lth ../engine/source/implicit/produt_v.F
2896!|| my_barrier ../engine/source/system/machine.F
2897!|| nlim_mix ../engine/source/implicit/imp_pcg.F
2898!|| produt_h ../engine/source/implicit/produt_v.F
2899!|| pupd ../engine/source/implicit/imp_pcg.F
2900!||--- uses -----------------------------------------------------
2901!|| dsgraph_mod ../engine/share/modules/dsgraph_mod.F
2902!|| groupdef_mod ../common_source/modules/groupdef_mod.F
2903!|| imp_pcg_proj ../engine/share/modules/impbufdef_mod.F
2904!|| imp_workh ../engine/share/modules/impbufdef_mod.F
2905!|| intbufdef_mod ../common_source/modules/interfaces/intbufdef_mod.F90
2906!||====================================================================
2907 SUBROUTINE IMP_PPCGH( IPREC ,
2908 1 NDDL ,NNZ ,IADK ,JDIK ,DIAG_K ,
2909 2 LT_K ,NDDLI ,ITOK ,IADI ,JDII ,
2910 3 LT_I ,NNZM ,IADM ,JDIM ,DIAG_M ,
2911 4 LT_M ,X ,R ,ITOL ,TOL ,
2912 5 P ,Z ,Y ,ITASK ,IPRINT ,
2913 6 N_MAX ,EPS_M ,F_X ,ISTOP ,W_DDL ,
2914 8 A ,AR ,VE ,MS ,XE ,
2915 9 D ,DR ,NDOF ,IPARI ,INTBUF_TAB,
2916 A NUM_IMP,NS_IMP,NE_IMP,NSREM ,
2917 B NSL ,NMONV ,IMONV ,MONVOL,IGRSURF ,
2918 C VOLMON,FR_MV ,IBFV ,SKEW ,
2919 D XFRAME ,GRAPHE,IAD_ELEM,FR_ELEM,ITAB ,
2920 E INSOLV ,ITN ,FAC_K ,IPIV_K,NK ,
2921 F MUMPS_PAR,CDDLP,ISOLV,IDSC ,IDDL ,
2922 G IKC ,INLOC ,IND_IMP,XI_C ,R0 ,
2923 H NDDLI_G,INTP_C,IRBE3 ,LRBE3 ,IRBE2 ,
2924 I LRBE2 )
2925C-----------------------------------------------
2926C M o d u l e s
2927C-----------------------------------------------
2928 USE DSGRAPH_MOD
2929 USE IMP_WORKH
2930 USE IMP_PCG_PROJ
2931 USE INTBUFDEF_MOD
2932 USE GROUPDEF_MOD
2933C-----------------------------------------------
2934C I m p l i c i t T y p e s
2935C-----------------------------------------------
2936#include "implicit_f.inc"
2937C-----------------------------------------------
2938C C o m m o n B l o c k s
2939C-----------------------------------------------
2940#if defined(MUMPS5)
2941#include "dmumps_struc.h"
2942#endif
2943#include "com04_c.inc"
2944#include "units_c.inc"
2945#include "task_c.inc"
2946C-----------------------------------------------
2947C D u m m y A r g u m e n t s
2948C-----------------------------------------------
2949C----------resol [K]{X}={F}---------
2950 INTEGER NDDL ,NNZ ,IADK(*) ,JDIK(*),ITOL,
2951 . NDDLI ,ITOK(*) ,IADI(*),JDII(*),N_MAX,
2952 . NNZM ,IADM(*),JDIM(*),IPREC,ITASK,IPRINT,
2953 . ISTOP,W_DDL(*),IBFV(*),INTP_C,IRBE3(*) ,LRBE3(*),
2954 . IRBE2(*) ,LRBE2(*)
2955 INTEGER NDOF(*),NE_IMP(*),NSREM ,NSL,
2956 . IPARI(*) ,NUM_IMP(*),NS_IMP(*) ,IND_IMP(*)
2957 INTEGER NMONV,IMONV(*),MONVOL(*),FR_MV(*)
2958 INTEGER IAD_ELEM(2,*), FR_ELEM(*), ITAB(*),
2959 . INSOLV, ITN, IPIV_K(*), NK, CDDLP(*), ISOLV, IDSC,
2960 . IDDL(*), IKC(*), INLOC(*),NDDLI_G
2961 my_real DIAG_K(*), LT_K(*) ,DIAG_M(*),LT_M(*) ,X(*), TOL, P(*) ,Z(*) ,R(*) ,Y(*),LT_I(*) ,EPS_M,F_X,XI_C(*)
2962 my_real A(3,*),AR(3,*),VE(3,*),D(3,*),DR(3,*),XE(3,*), MS(*),VOLMON(*),SKEW(*) ,XFRAME(*),FAC_K(*), R0(*)
2963 TYPE(PRGRAPH) :: GRAPHE(*)
2964#ifdef MUMPS5
2965 TYPE(DMUMPS_STRUC) MUMPS_PAR
2966#else
2967 ! Fake declaration as DMUMPS_STRUC is shipped with MUMPS
2968 INTEGER MUMPS_PAR
2969#endif
2970 TYPE(INTBUF_STRUCT_) INTBUF_TAB(*)
2971 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
2972#ifdef MUMPS5
2973C-----------------------------------------------
2974C L o c a l V a r i a b l e s
2975C-----------------------------------------------
2976 INTEGER I,J,IT,IP,NLIM,ND,IERR,IPRI,IERROR,NNZI,LOC_PROC,ICK
2977 PARAMETER (ND=4)
2978 CHARACTER*4 EXIT
2979 CHARACTER*11 WARR
2980 my_real S , R2, R02,ALPHA,BETA,G0,G1,RR,TOLS,TOLN,TOLS2
2981 INTEGER CRIT_STOP,IUPD,IUPD0,F_DDL,L_DDL,IFLG,GPUR4R8,ITP
2982 EXTERNAL CRIT_STOP
2983 my_real ANORM2,XNORM2,L_A,L_B2,L_B,A_OLD,B_OLD,TMP,G00,R2_OLD,RMAX
2984 my_real CS,DBAR, DELTA, DENOM, KCOND,SNPROD,QRNORM,GAMMA, GBAR, GMAX, GMIN, EPSLN,LQNORM,DIAG,CGNORM,
2985 . OLDB, RHS1, RHS2,SN, ZBAR, ZL ,OLDB2,TNORM2,EPS(4)
2986C--------------P-PCG Usage--------------------------
2987 my_real, DIMENSION(:),ALLOCATABLE :: SQ_DIAG_M
2988 SAVE SQ_DIAG_M
2989C--------------END GPU SPECIFIC----------------------
2990 DATA EXIT /'with'/
2991 DATA WARR /'**warring**'/
2992C-----------------------------------------------
2993 LOC_PROC=ISPMD+1
2994C--------------INITIALISATION--------------------------
2995C-----Istop=1 : NIT>=N_LIM+ RR>TOL
2996C------X(I)=ZERO----******N_MAX,EPS_M a calculer avant----
2997 F_DDL=1+ITASK*NDDL/NTHREAD
2998 L_DDL=(ITASK+1)*NDDL/NTHREAD
2999C
3000 ISTOP=0
3001 NLIM=N_MAX
3002 CALL NLIM_MIX(N_MAX,NDDLI_G,NLIM)
3003 NLIM=MAX(NLIM,2)
3004 TOLS=MAX(TOL,EPS_M)
3005 IPRI = IPRINT
3006 IF (ISPMD/=0) IPRI = 0
3007C--------
3008 IF (ITASK == 0) THEN
3009 IF(IPRI/=0)THEN
3010 WRITE(IOUT,*)'*** begin projected conjugate gradient iteration ***'
3011 WRITE(IOUT,*)
3012 ENDIF
3013 IF(IPRI<0)THEN
3014 WRITE(ISTDO,*)'*** begin projected conjugate gradient iteration ***'
3015 WRITE(ISTDO,*)
3016 ENDIF
3017 END IF !(ITASK == 0) THEN
3018 IP = IABS(IPRI)
3019 IT=0
3020C-------------Vector Z is no more necessary-> work array for A'*v (MMAV_LTH)
3021 IUPD = 0
3022 IUPD0 = 0
3023 IF (IPREC>1) THEN
3024 IF (ITASK == 0) ALLOCATE (SQ_DIAG_M(NDDL))
3025C----------------------
3026 CALL MY_BARRIER
3027C--------For S ini only-------------
3028 DO I = F_DDL,L_DDL
3029 SQ_DIAG_M(I) = SQRT(DIAG_M(I))
3030 ENDDO
3031 END IF !(IPREC>1) THEN
3032C----------------------
3033 CALL MY_BARRIER
3034C--------For S ini only----X is in place of R to not change R---------
3035 RMAX=ZERO
3036 CALL IMP_INISI(
3037 1 NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
3038 2 NDDLI ,ITOK ,IADI ,JDII ,LT_I ,
3039 3 NNZM ,IADM ,JDIM ,DIAG_M ,LT_M ,
3040 4 X ,P ,Y ,
3041 8 A ,AR ,VE ,MS ,XE ,
3042 9 D ,DR ,NDOF ,IPARI ,INTBUF_TAB,
3043 A NUM_IMP,NS_IMP,NE_IMP,NSREM ,
3044 B NSL ,NMONV ,IMONV ,MONVOL,IGRSURF ,
3045 C VOLMON,FR_MV ,IBFV ,SKEW ,
3046 D XFRAME ,Z ,IDDL ,IKC ,INLOC ,
3047 G IND_IMP,XI_C ,IRBE3 ,LRBE3 ,IRBE2 ,
3048 H LRBE2 ,ITASK ,W_DDL ,F_DDL ,L_DDL ,
3049 I SQ_DIAG_M)
3050C /---------------/
3051 CALL MY_BARRIER
3052C /---------------/
3053 CALL IMP_INIST(
3054 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K ,
3055 2 LT_K ,ITOK ,IADI ,JDII ,LT_I ,
3056 3 ITASK ,W_DDL ,A ,AR ,VE ,
3057 4 D ,DR ,NDOF ,IPARI ,INTBUF_TAB,
3058 5 NUM_IMP,NS_IMP,NE_IMP,NSREM ,
3059 6 NSL ,NMONV ,IMONV ,MONVOL,IGRSURF ,
3060 7 VOLMON,FR_MV ,IBFV ,SKEW ,
3061 8 XFRAME ,IND_IMP,XI_C ,MS ,XE ,
3062 9 IRBE3 ,LRBE3 ,IRBE2 ,LRBE2 ,IADM ,
3063 A JDIM ,SQ_DIAG_M, LT_M ,F_DDL ,L_DDL,
3064 B Z )
3065C /---------------/
3066 CALL MY_BARRIER
3067C /---------------/
3068C--------b'=Lm^t b--------------
3069 CALL MMV_LTH(NDDL ,IADM ,JDIM ,SQ_DIAG_M,LT_M ,
3070 2 R ,Z ,F_DDL ,L_DDL ,ITASK )
3071C----------------------
3072 CALL MY_BARRIER
3073C---------------------
3074 DO I = F_DDL,L_DDL
3075 R(I) = Z(I)
3076 ENDDO
3077 CALL PRODUT_H(F_DDL ,L_DDL ,R ,R ,W_DDL , G0 ,ITASK )
3078C------x' ini
3079 CALL IMP_INIX(F_DDL,L_DDL,NDDL,X ,R ,W_DDL,ITASK )
3080C----------------------
3081 CALL MY_BARRIER
3082C---------------------
3083 CALL MMAV_LTH(
3084 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K,
3085 2 LT_K ,IADI ,JDII ,ITOK ,LT_I ,
3086 3 X ,Y ,A ,AR ,
3087 5 VE ,MS ,XE ,D ,DR ,
3088 6 NDOF ,IPARI ,INTBUF_TAB ,NUM_IMP,
3089 7 NS_IMP,NE_IMP,NSREM ,NSL ,IBFV ,
3090 8 SKEW ,XFRAME,MONVOL,VOLMON,IGRSURF ,
3091 9 FR_MV,NMONV ,IMONV ,IND_IMP,
3092 A XI_C ,IUPD ,IRBE3 ,LRBE3 ,IRBE2 ,
3093 B LRBE2 ,IADM ,JDIM ,SQ_DIAG_M,LT_M ,
3094 C F_DDL ,L_DDL ,ITASK ,Z )
3095C----------------------
3096 CALL MY_BARRIER
3097C---------------------
3098 IF (ITOL==1) THEN
3099 R02 = G0
3100 ELSEIF (ITOL==3) THEN
3101C----------initialization of XNORM2
3102 CALL PRODUT_H( F_DDL ,L_DDL ,X ,X ,W_DDL , XNORM2 ,ITASK )
3103 END IF
3104C
3105 DO I=F_DDL,L_DDL
3106 R(I) = R(I)-Y(I)
3107 ENDDO
3108C----------------------
3109 CALL MY_BARRIER
3110C---------------------
3111C CALL PREC_SOLVH(IPREC, ITASK ,
3112c 1 GRAPHE,IAD_ELEM,FR_ELEM,DIAG_K,LT_K ,
3113c 2 IADK ,JDIK ,ITAB ,IPRI,INSOLV ,
3114c 3 ITN ,FAC_K , IPIV_K, NK ,MUMPS_PAR,
3115c 4 CDDLP ,ISOLV , IDSC , IDDL ,IKC ,
3116c 5 INLOC ,NDOF , NDDL ,NNZM ,IADM ,
3117c 6 JDIM ,DIAG_M , LT_M ,R ,Z ,
3118c 7 F_DDL ,L_DDL )
3119 DO I = F_DDL,L_DDL
3120 P(I) = R(I)
3121 ENDDO
3122C----------------------
3123 CALL MY_BARRIER
3124C---------------------
3125 CALL IMP_PRO_P(F_DDL,L_DDL,NDDL ,P ,W_DDL,ITASK )
3126C----------------------
3127 CALL MY_BARRIER
3128C---------------------
3129 CALL PRODUT_H(F_DDL ,L_DDL ,R ,R ,W_DDL , G0 ,ITASK )
3130C---------------------
3131 CALL MMAV_LTH(
3132 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K,
3133 2 LT_K ,IADI ,JDII ,ITOK ,LT_I ,
3134 3 P ,Y ,A ,AR ,
3135 5 VE ,MS ,XE ,D ,DR ,
3136 6 NDOF ,IPARI ,INTBUF_TAB ,NUM_IMP,
3137 7 NS_IMP,NE_IMP,NSREM ,NSL ,IBFV ,
3138 8 SKEW ,XFRAME,MONVOL,VOLMON,IGRSURF ,
3139 9 FR_MV,NMONV ,IMONV ,IND_IMP,
3140 A XI_C ,IUPD ,IRBE3 ,LRBE3 ,IRBE2 ,
3141 B LRBE2 ,IADM ,JDIM ,SQ_DIAG_M,LT_M ,
3142 C F_DDL ,L_DDL ,ITASK ,Z )
3143C----------------------
3144 CALL MY_BARRIER
3145C---------------------
3146 CALL IMP_UPDV2(
3147 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K ,
3148 2 LT_K ,ITOK ,IADI ,JDII ,LT_I ,
3149 3 ITASK ,W_DDL ,A ,AR ,VE ,
3150 4 D ,DR ,NDOF ,IPARI ,INTBUF_TAB,
3151 5 NUM_IMP,NS_IMP,NE_IMP,NSREM ,
3152 6 NSL ,NMONV ,IMONV ,MONVOL,IGRSURF ,
3153 7 VOLMON,FR_MV ,IBFV ,SKEW ,
3154 8 XFRAME ,IND_IMP,XI_C ,MS ,XE ,
3155 9 IRBE3 ,LRBE3 ,IRBE2 ,LRBE2 ,PROJ_V,
3156 A P ,Y ,F_DDL ,L_DDL ,0 ,
3157 B 0 ,Z ,IADM ,JDIM ,SQ_DIAG_M,
3158 C LT_M )
3159 CALL PRODUT_H(F_DDL ,L_DDL ,P ,Y ,W_DDL , S ,ITASK )
3160C---------------------
3161 IF (G0==ZERO) THEN
3162 IF (ITOL>1) THEN
3163 ISTOP=-1
3164 GOTO 200
3165 ELSE
3166 ALPHA = ZERO
3167 ENDIF
3168 ELSE
3169 ALPHA = G0/S
3170 ENDIF
3171 TOLS2=TOLS*TOLS
3172 IF (ITOL==1) THEN
3173C---------------------
3174 R2 =G0
3175 ELSEIF (ITOL==3) THEN
3176C------ R2<TOL*TOL*ANORM2*XNORM2------
3177 R02=ABS(G0)
3178 R2 =R02
3179 L_A=ONE/ALPHA
3180C L_B2=R02
3181 TNORM2=L_A*L_A
3182 ANORM2=ZERO
3183 IF (M_VS == 0) XNORM2=ZERO
3184 A_OLD=L_A
3185 B_OLD=ZERO
3186 OLDB = SQRT(R02)
3187 ELSEIF (ITOL==4) THEN
3188 R02=ALPHA*ALPHA*ABS(G0)
3189 EPS(1)=R02
3190 R2=R02
3191 ELSE
3192 R02=ABS(G0)
3193 R2 =R02
3194 ENDIF
3195 IF (R02==ZERO) THEN
3196 ISTOP=-1
3197 GOTO 200
3198 ENDIF
3199 TOLN=R02*TOLS2
3200.AND. IF(IPRI/=0ITASK==0)THEN
3201 RR = SQRT(R2)
3202 IF(IPRI<0) WRITE(ISTDO,1000)IT,RR
3203 WRITE(IOUT,1000)IT,RR
3204 ENDIF
3205C-------pour etre coherent avec lanzos for linear
3206C----------------------
3207 CALL MY_BARRIER
3208C---------------------
3209 IT=1
3210 DO I=F_DDL,L_DDL
3211 X(I) = X(I) + ALPHA*P(I)
3212 R(I) = R(I) - ALPHA*Y(I)
3213 ENDDO
3214C----------------------
3215 CALL MY_BARRIER
3216C---------------------
3217c CALL PREC_SOLVH(IPREC, ITASK ,
3218c 1 GRAPHE,IAD_ELEM,FR_ELEM,DIAG_K,LT_K ,
3219c 2 IADK ,JDIK ,ITAB ,IPRI,INSOLV ,
3220c 3 ITN ,FAC_K , IPIV_K, NK ,MUMPS_PAR,
3221c 4 CDDLP ,ISOLV , IDSC , IDDL ,IKC ,
3222c 5 INLOC ,NDOF , NDDL ,NNZM ,IADM ,
3223c 6 JDIM ,DIAG_M , LT_M ,R ,Z ,
3224c 7 F_DDL ,L_DDL )
3225C---------------------
3226 CALL PRODUT_H( F_DDL ,L_DDL ,R ,R ,W_DDL , G1 ,ITASK )
3227C---------------------
3228 BETA=G1/G0
3229 IF (ITOL==1) THEN
3230 R2 = G1
3231 ELSEIF (ITOL==3) THEN
3232C------ R2<TOL*TOL*ANORM2*XNORM2------
3233 R2=ABS(G1)
3234 L_B2=ABS(BETA)*A_OLD*A_OLD
3235 L_B=SQRT(L_B2)
3236 TNORM2=TNORM2+L_B2
3237 B_OLD=BETA
3238C* INITIALIZE OTHER QUANTITIES.
3239 GBAR = L_A
3240 DBAR = L_B
3241 RHS1 = OLDB
3242 RHS2 = ZERO
3243 GMAX = ABS( L_A ) + EPS_M
3244 GMIN = GMAX
3245 OLDB2 = L_B2
3246 OLDB = L_B
3247 ELSEIF (ITOL==4) THEN
3248 R2=R02
3249 ELSE
3250 R2=ABS(G1)
3251 ENDIF
3252 G0 = G1
3253 IF (ITOL==3) TOLN=TOLN*TNORM2
3254C----------------------
3255 CALL MY_BARRIER
3256C---------------------
3257 ISTOP=CRIT_STOP(IT,R2,NLIM,TOLN)
3258.AND. IF(NDDLI_G>0INTP_C==-1)IUPD0 = -INTP_C
3259 DO WHILE (ISTOP==1)
3260 DO I=F_DDL,L_DDL
3261 P(I) = R(I) + BETA*P(I)
3262 ENDDO
3263.AND. IF(IUPD0>0IT>NDDLI_G/2) IUPD = IUPD0
3264C
3265C------PCG(PROJECTION)----
3266 CALL IMP_PRO_P(F_DDL,L_DDL,NDDL ,P ,W_DDL,ITASK )
3267C
3268C----------------------
3269 CALL MY_BARRIER
3270C---------------------
3271 CALL MMAV_LTH(
3272 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K,
3273 2 LT_K ,IADI ,JDII ,ITOK ,LT_I ,
3274 3 P ,Y ,A ,AR ,
3275 5 VE ,MS ,XE ,D ,DR ,
3276 6 NDOF ,IPARI ,INTBUF_TAB ,NUM_IMP,
3277 7 NS_IMP,NE_IMP,NSREM ,NSL ,IBFV ,
3278 8 SKEW ,XFRAME,MONVOL,VOLMON,IGRSURF ,
3279 9 FR_MV,NMONV ,IMONV ,IND_IMP,
3280 A XI_C ,IUPD ,IRBE3 ,LRBE3 ,IRBE2 ,
3281 B LRBE2 ,IADM ,JDIM ,SQ_DIAG_M,LT_M ,
3282 C F_DDL ,L_DDL ,ITASK ,Z )
3283C----------------------
3284 CALL MY_BARRIER
3285C---------------------
3286 ICK=0
3287 CALL IMP_UPDV2(
3288 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K ,
3289 2 LT_K ,ITOK ,IADI ,JDII ,LT_I ,
3290 3 ITASK ,W_DDL ,A ,AR ,VE ,
3291 4 D ,DR ,NDOF ,IPARI ,INTBUF_TAB,
3292 5 NUM_IMP,NS_IMP,NE_IMP,NSREM ,
3293 6 NSL ,NMONV ,IMONV ,MONVOL,IGRSURF ,
3294 7 VOLMON,FR_MV ,IBFV ,SKEW ,
3295 8 XFRAME ,IND_IMP,XI_C ,MS ,XE ,
3296 9 IRBE3 ,LRBE3 ,IRBE2 ,LRBE2 ,PROJ_V,
3297 A P ,Y ,F_DDL ,L_DDL ,IT ,
3298 B ICK ,Z ,IADM ,JDIM ,SQ_DIAG_M,
3299 C LT_M )
3300C----------------------
3301 CALL MY_BARRIER
3302C---------------------
3303 CALL PRODUT_H(F_DDL ,L_DDL ,P ,Y ,W_DDL , S ,ITASK )
3304C---------------------
3305 ALPHA=G0/S
3306 IF(MOD(IT,10)<0)THEN
3307 DO I=F_DDL,L_DDL
3308 X(I) = X(I) + ALPHA*P(I)
3309 ENDDO
3310 CALL MMAV_LTH(
3311 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K,
3312 2 LT_K ,IADI ,JDII ,ITOK ,LT_I ,
3313 3 X ,Y ,A ,AR ,
3314 5 VE ,MS ,XE ,D ,DR ,
3315 6 NDOF ,IPARI ,INTBUF_TAB ,NUM_IMP,
3316 7 NS_IMP,NE_IMP,NSREM ,NSL ,IBFV ,
3317 8 SKEW ,XFRAME,MONVOL,VOLMON,IGRSURF ,
3318 9 FR_MV,NMONV ,IMONV ,IND_IMP,
3319 A XI_C ,IUPD ,IRBE3 ,LRBE3 ,IRBE2 ,
3320 B LRBE2 ,IADM ,JDIM ,SQ_DIAG_M,LT_M ,
3321 C F_DDL ,L_DDL ,ITASK ,Z )
3322 DO I=F_DDL,L_DDL
3323 R(I) = R0(I) - Y(I)
3324 ENDDO
3325C----------------------
3326 CALL MY_BARRIER
3327C---------------------
3328 ELSE
3329 DO I=F_DDL,L_DDL
3330 X(I) = X(I) + ALPHA*P(I)
3331 R(I) = R(I) - ALPHA*Y(I)
3332 ENDDO
3333 END IF ! IF(MOD(IT,IP)==0)THEN
3334C----------------------
3335 CALL MY_BARRIER
3336C---------------------
3337c CALL PREC_SOLVH(IPREC, ITASK ,
3338c 1 GRAPHE,IAD_ELEM,FR_ELEM,DIAG_K,LT_K ,
3339c 2 IADK ,JDIK ,ITAB ,IPRI,INSOLV ,
3340c 3 ITN ,FAC_K , IPIV_K, NK ,MUMPS_PAR,
3341c 4 CDDLP ,ISOLV , IDSC , IDDL ,IKC ,
3342c 5 INLOC ,NDOF , NDDL ,NNZM ,IADM ,
3343c 6 JDIM ,DIAG_M , LT_M ,R ,Z ,
3344c 7 F_DDL ,L_DDL )
3345C---------------------
3346 CALL PRODUT_H(F_DDL ,L_DDL ,R ,R ,W_DDL , G1 ,ITASK )
3347C---------------------
3348 BETA=G1/G0
3349 G0 = G1
3350 IF (ITOL==1) THEN
3351 R2 = G1
3352 ELSEIF (ITOL==3) THEN
3353 R2 =ABS(G1)
3354 S=ONE/ALPHA
3355 L_A=S+B_OLD*A_OLD
3356C----------L_B2 : beta(j-1)^2-A_OLD :1/alpha(j-1)-------
3357 L_B2=ABS(BETA)*S*S
3358 L_B=SQRT(L_B2)
3359 A_OLD=S
3360 B_OLD=BETA
3361 ANORM2=TNORM2
3362 TNORM2=TNORM2+L_A*L_A+OLDB2+L_B2
3363 GAMMA = SQRT( GBAR*GBAR + OLDB2 )
3364 CS = GBAR / GAMMA
3365 SN = OLDB / GAMMA
3366 DELTA = CS * DBAR + SN * L_A
3367 GBAR = SN * DBAR - CS * L_A
3368 EPSLN = SN * L_B
3369 DBAR = - CS * L_B
3370 ZL = RHS1 / GAMMA
3371 XNORM2 = XNORM2+ZL*ZL
3372 GMAX = MAX( GMAX, GAMMA )
3373 GMIN = MIN( GMIN, GAMMA )
3374 RHS1 = RHS2 - DELTA * ZL
3375 RHS2 = - EPSLN * ZL
3376 TOLN=TOLS2*ANORM2*XNORM2
3377 OLDB2 = L_B2
3378 OLDB = L_B
3379 ELSEIF (ITOL==4) THEN
3380 TMP=ALPHA*ALPHA*ABS(G1)
3381 IF (IT>=ND) THEN
3382 DO J=1,ND-1
3383 EPS(J)=EPS(J+1)
3384 ENDDO
3385 EPS(ND)=TMP
3386 R2=ZERO
3387 DO J=1,ND
3388 R2=R2+EPS(J)
3389 ENDDO
3390 ELSE
3391 EPS(IT+1)=TMP
3392 R2=R2+TMP
3393 ENDIF
3394 ELSE
3395 R2 =ABS(G1)
3396 ENDIF
3397C
3398 CALL IMP_UPDV1(
3399 1 R2 ,R2_OLD ,RMAX ,PROJ_V,P ,
3400 2 IT ,F_DDL,L_DDL )
3401C
3402.AND. IF(IPRI/=0ITASK==0)THEN
3403 IF(MOD(IT,IP)==0)THEN
3404 RR = SQRT(R2/R02)
3405 WRITE(IOUT,1001)IT,RR
3406 IF(IPRI<0) WRITE(ISTDO,1001)IT,RR
3407 ENDIF
3408 ENDIF
3409C----------------------
3410 CALL MY_BARRIER
3411C---------------------
3412 ISTOP=CRIT_STOP(IT,R2,NLIM,TOLN)
3413.AND..OR. IF((IUPD>0IT==NLIM)
3414.AND..AND. . (IUPD==0ISTOP/=1IUPD0>0))THEN
3415C------- re-iter with linear Kint------
3416 IUPD0 = 0
3417 IF(IUPD>0) IUPD = 0
3418 ISTOP = 1
3419 NLIM = NLIM + NLIM
3420 DO I=F_DDL,L_DDL
3421 X(I) = ZERO
3422 ENDDO
3423C----------------------
3424 CALL MY_BARRIER
3425C---------------------
3426 CALL MMAV_LTH(
3427 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K,
3428 2 LT_K ,IADI ,JDII ,ITOK ,LT_I ,
3429 3 X ,Y ,A ,AR ,
3430 5 VE ,MS ,XE ,D ,DR ,
3431 6 NDOF ,IPARI ,INTBUF_TAB ,NUM_IMP,
3432 7 NS_IMP,NE_IMP,NSREM ,NSL ,IBFV ,
3433 8 SKEW ,XFRAME,MONVOL,VOLMON,IGRSURF ,
3434 9 FR_MV,NMONV ,IMONV ,IND_IMP,
3435 A XI_C ,IUPD ,IRBE3 ,LRBE3 ,IRBE2 ,
3436 B LRBE2 ,IADM ,JDIM ,SQ_DIAG_M,LT_M ,
3437 C F_DDL ,L_DDL ,ITASK ,Z )
3438C----------------------
3439 CALL MY_BARRIER
3440C---------------------
3441 DO I=F_DDL,L_DDL
3442 R(I) = R0(I)-Y(I)
3443 ENDDO
3444C----------------------
3445 CALL MY_BARRIER
3446c CALL PREC_SOLVH(IPREC, ITASK ,
3447c 1 GRAPHE,IAD_ELEM,FR_ELEM,DIAG_K,LT_K ,
3448c 2 IADK ,JDIK ,ITAB ,IPRI,INSOLV ,
3449c 3 ITN ,FAC_K , IPIV_K, NK ,MUMPS_PAR,
3450c 4 CDDLP ,ISOLV , IDSC , IDDL ,IKC ,
3451c 5 INLOC ,NDOF , NDDL ,NNZM ,IADM ,
3452c 6 JDIM ,DIAG_M , LT_M ,R ,Z ,
3453c 7 F_DDL ,L_DDL )
3454C---------------------
3455 CALL PRODUT_H( F_DDL ,L_DDL ,R ,R ,W_DDL , G0 ,ITASK )
3456 BETA = ZERO
3457C------ R2<TOL*TOL*ANORM2*XNORM2------
3458 IF (ITOL==3) THEN
3459 TNORM2=L_A*L_A
3460 ANORM2=ZERO
3461 XNORM2=ZERO
3462 A_OLD=ZERO
3463 B_OLD=ZERO
3464 L_B=ZERO
3465C* INITIALIZE OTHER QUANTITIES.
3466 GBAR = L_A
3467 DBAR = ZERO
3468 RHS1 = SQRT(R02)
3469 RHS2 = ZERO
3470 GMAX = ABS( L_A ) + EPS_M
3471 GMIN = GMAX
3472 OLDB2 = ZERO
3473 OLDB = L_B
3474 ELSEIF (ITOL==4) THEN
3475 DO J=1,ND
3476 EPS(J)=R02
3477 ENDDO
3478 ENDIF
3479 ENDIF
3480C----------------------
3481 CALL MY_BARRIER
3482C---------------------
3483 IT = IT +1
3484 ENDDO
3485 200 CONTINUE
3486C--------PCG (PROJECTION)
3487 IF (ISTOP==0) THEN
3488 CALL IMP_UPDST(
3489 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K ,
3490 2 LT_K ,ITOK ,IADI ,JDII ,LT_I ,
3491 3 ITASK ,W_DDL ,A ,AR ,VE ,
3492 4 D ,DR ,NDOF ,IPARI ,INTBUF_TAB,
3493 5 NUM_IMP,NS_IMP,NE_IMP,NSREM ,
3494 6 NSL ,NMONV ,IMONV ,MONVOL,IGRSURF ,
3495 7 VOLMON,FR_MV ,IBFV ,SKEW ,
3496 8 XFRAME ,IND_IMP,XI_C ,MS ,XE ,
3497 9 IRBE3 ,LRBE3 ,IRBE2 ,LRBE2 ,X ,
3498 A P ,Y ,IADM ,JDIM ,SQ_DIAG_M ,
3499 B LT_M ,F_DDL ,L_DDL )
3500C----------------------
3501 CALL MY_BARRIER
3502C----------------------
3503C--------x=Lm x'--------------
3504 CALL MMV_LH(NDDL ,IADM0 ,JDIM0 ,SQ_DIAG_M,LT_M0 ,
3505 2 X ,Z ,F_DDL ,L_DDL ,ITASK )
3506 DO I=F_DDL,L_DDL
3507 X(I) = Z(I)
3508 ENDDO
3509 IF (ITASK == 0) THEN
3510 NCG_RUN = NCG_RUN + 1
3511 IF (IPREC>1) DEALLOCATE(SQ_DIAG_M)
3512 END IF
3513C----------------------
3514 CALL MY_BARRIER
3515C---------------------
3516 END IF
3517C------due to ISTOP local variable
3518C IF(IT>=NLIM.AND.ISOLV==7) ISTOP = 3
3519c IF(ITOL==3)then
3520c DO I=F_DDL,L_DDL
3521c Z(I)=X(I)/DIAG_M(I)
3522c ENDDO
3523C----------------------
3524c CALL MY_BARRIER
3525C---------------------
3526c CALL PRODUT_H( F_DDL ,L_DDL ,X ,Z ,W_DDL , XNORM2 ,ITASK )
3527c END IF
3528 IF(ITASK == 0) THEN
3529 ISTOP_H = ISTOP
3530 CALL PUPD(N_MAX,NDDLI_G,IT)
3531C
3532 IF(IT>=NLIM)THEN
3533 IF (ISOLV==7) THEN
3534 ISTOP_H = 3
3535 ELSE
3536
3537 WRITE(IOUT,*)
3538 WRITE(IOUT,1003)NLIM
3539 WRITE(IOUT,*)
3540 WRITE(ISTDO,*)
3541 WRITE(ISTDO,1003)NLIM
3542 WRITE(ISTDO,*)
3543 IF(ITOL==3)THEN
3544 WRITE(IOUT,2000) EXIT, SQRT(ANORM2), SQRT(XNORM2),GMAX/GMIN
3545 WRITE(ISTDO,2000) EXIT, SQRT(ANORM2), SQRT(XNORM2),GMAX/GMIN
3546 ENDIF
3547 ISTOP_H = 1
3548 IF (R2>HUNDRED*TOLN) THEN
3549 ISTOP_H = 2
3550 RR = SQRT(R2/R02)
3551 WRITE(IOUT,*)
3552 WRITE(IOUT,1004)RR
3553 WRITE(IOUT,*)
3554 WRITE(ISTDO,*)
3555 WRITE(ISTDO,1004)RR
3556 WRITE(ISTDO,*)
3557 ENDIF
3558
3559 END IF !(ISOLV==7) THEN
3560 ENDIF
3561 IF(IPRI/=0)THEN
3562 RR = SQRT(R2/R02)
3563 WRITE(IOUT,*)
3564 WRITE(IOUT,1002)IT,RR
3565 IF(ITOL==3)WRITE(IOUT,2000) EXIT, SQRT(ANORM2), SQRT(XNORM2),GMAX/GMIN
3566C WRITE(IOUT,*) 'TOLN,TOLS2=',TOLN,TOLS2
3567 WRITE(IOUT,*)
3568 IF(IPRI<0) THEN
3569 WRITE(ISTDO,*)
3570 WRITE(ISTDO,1002)IT,RR
3571 IF(ITOL==3)WRITE(ISTDO,2000)EXIT, SQRT(ANORM2), SQRT(XNORM2),GMAX/GMIN
3572 WRITE(ISTDO,*)
3573 ENDIF
3574 ENDIF
3575 IF(ITOL==3) THEN
3576 K_LAMDA0= GMIN
3577 K_LAMDA1= GMAX
3578 ELSE
3579 K_LAMDA0= ZERO
3580 K_LAMDA1= ZERO
3581 ENDIF
3582 END IF !(ITASK == 0) THEN
3583C----------------------
3584 CALL MY_BARRIER
3585C----------------------
3586 ISTOP = ISTOP_H
3587C--------------------------------------------
3588 1000 FORMAT(5X,'iteration=',I8,5X,' initial residual norm=',E11.4)
3589 1001 FORMAT(5X,'iteration=',I8,5X,' relative residual norm=',E11.4)
3590 1002 FORMAT(3X,'total c.g. iteration=',I8,5X,
3591 . ' relative residual norm=',E11.4)
3592 1003 FORMAT(5X,
3593 . '---warning : the iteration limit number was reached',I8)
3594 1004 FORMAT(5X,'warning:c.g stop with relative residual norm=',E11.4)
3595 2000 FORMAT(/ 5X, A, 2X, 'anorm =', E12.4, 2X, 'xnorm =', E12.4,2X,
3596 . 'kcond =', E12.4)
3597 2002 FORMAT(/ 5X, 'with', 2X, 'alfa =', E12.4, 2X, 'beta =',
3598 . E12.4,2X,'oldb =', E12.4)
3599 RETURN
3600#endif
3601 END
3602!||====================================================================
3603!|| imp_inist ../engine/source/implicit/imp_pcg.F
3604!||--- called by ------------------------------------------------------
3605!|| imp_ppcgh ../engine/source/implicit/imp_pcg.F
3606!||--- calls -----------------------------------------------------
3607!|| arret ../engine/source/system/arret.F
3608!|| mam_nm ../engine/source/implicit/produt_v.F
3609!|| mav_mm ../engine/source/implicit/produt_v.F
3610!|| mmav_lth ../engine/source/implicit/produt_v.F
3611!|| my_barrier ../engine/source/system/machine.F
3612!||--- uses -----------------------------------------------------
3613!|| groupdef_mod ../common_source/modules/groupdef_mod.F
3614!|| imp_pcg_proj ../engine/share/modules/impbufdef_mod.F
3615!|| intbufdef_mod ../common_source/modules/interfaces/intbufdef_mod.F90
3616!||====================================================================
3617 SUBROUTINE IMP_INIST(
3618 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K ,
3619 2 LT_K ,ITOK ,IADI ,JDII ,LT_I ,
3620 3 ITASK ,W_DDL ,A ,AR ,VE ,
3621 4 D ,DR ,NDOF ,IPARI ,INTBUF_TAB ,
3622 5 NUM_IMP,NS_IMP,NE_IMP,NSREM ,
3623 6 NSL ,NMONV ,IMONV ,MONVOL,IGRSURF ,
3624 7 VOLMON,FR_MV ,IBFV ,SKEW ,
3625 8 XFRAME ,IND_IMP,XI_C ,MS ,XE ,
3626 9 IRBE3 ,LRBE3 ,IRBE2 ,LRBE2 ,IADM ,
3627 A JDIM ,DIAG_M, LT_M ,F_DDL ,L_DDL ,
3628 B V_W )
3629C-----------------------------------------------
3630C M o d u l e s
3631C-----------------------------------------------
3632 USE IMP_PCG_PROJ
3633 USE INTBUFDEF_MOD
3634 USE GROUPDEF_MOD
3635C-----------------------------------------------
3636C I m p l i c i t T y p e s
3637C-----------------------------------------------
3638#include "implicit_f.inc"
3639C-----------------------------------------------
3640C C o m m o n B l o c k s
3641C-----------------------------------------------
3642#include "com04_c.inc"
3643#include "impl1_c.inc"
3644#include "units_c.inc"
3645C-----------------------------------------------
3646C D u m m y A r g u m e n t s
3647C-----------------------------------------------
3648 INTEGER NDDL ,IADK(*) ,JDIK(*), NDDLI ,ITOK(*) ,IADI(*),JDII(*),
3649 . ITASK,W_DDL(*),IBFV(*),IRBE3(*) ,LRBE3(*), IRBE2(*) ,LRBE2(*),F_DDL,L_DDL,IADM(*) ,JDIM(*)
3650 INTEGER NDOF(*),NE_IMP(*),NSREM ,NSL,IPARI(*) ,NUM_IMP(*),NS_IMP(*) ,IND_IMP(*)
3651 INTEGER NMONV,IMONV(*),MONVOL(*),FR_MV(*)
3652 my_real DIAG_K(*), LT_K(*),LT_I(*) ,XI_C(*),DIAG_M(*), LT_M(*)
3653 my_real A(3,*),AR(3,*),VE(3,*),D(3,*),DR(3,*),XE(3,*),MS(*),VOLMON(*),SKEW(*) ,XFRAME(*),V_W(*)
3654
3655 TYPE(INTBUF_STRUCT_) INTBUF_TAB(*)
3656 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
3657C-----------------------------------------------
3658c FUNCTION: initialization of S,T of Projection taking into account to [M] preconditioner
3659c
3660c Note:
3661c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
3662c
3663c TYPE NAME FUNCTION
3664c I NDDL,NNZ - dimension of [K] and number of non zero (strict triangular)
3665c I IADK,JDIK - indice arrays for compressed row(col.) format of [K]
3666c I DIAG_K(NDDL) - diagonal terms of [K]
3667c I LT_K(NNZ) - strict triangular of [K]
3668c O Proj_S(NDDL,M_VS) - [S] reduced small Eigenvectors
3669c O Proj_T(NDDL,M_VS) - [T] =[K][S]
3670c I W1,W2 - work arrays
3671C-----------------------------------------------
3672C L o c a l V a r i a b l e s
3673C-----------------------------------------------
3674#if defined(MYREAL8) && !defined(WITHOUT_LINALG)
3675 CHARACTER JOBZ, UPLO
3676 INTEGER I,J,IT,IP,NLIM,ND,IUPD,IPRI,IERROR,NNZI,M,
3677 . F_DDLI,L_DDLI,INFO,LW,INORM,M_VS1
3678 my_real
3679 . WORK(3*M_VS+3),W(M_VS+1)
3680CC---------------------
3681C IF (NPCGPV >=0) RETURN
3682C----------------------
3683C CALL MY_BARRIER
3684C---------------------
3685 IF (NCG_RUN==0) THEN
3686 M_VS1 =M_VS
3687 ELSE
3688 M_VS1 =M_VS+1
3689 END IF
3690 IF (ITASK==0) NPCGPV = M_VS
3691C----------------------
3692 CALL MY_BARRIER
3693C---------------------
3694C---------------------
3695C T=[K][S]
3696C---------------------
3697 IUPD = 0
3698 DO M = 1,M_VS1
3699 CALL MMAV_LTH(
3700 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K,
3701 2 LT_K ,IADI ,JDII ,ITOK ,LT_I ,
3702 3 PROJ_S(1,M),PROJ_T(1,M),A ,AR ,
3703 5 VE ,MS ,XE ,D ,DR ,
3704 6 NDOF ,IPARI ,INTBUF_TAB ,NUM_IMP,
3705 7 NS_IMP,NE_IMP,NSREM ,NSL ,IBFV ,
3706 8 SKEW ,XFRAME,MONVOL,VOLMON,IGRSURF ,
3707 9 FR_MV,NMONV ,IMONV ,IND_IMP,
3708 A XI_C ,IUPD ,IRBE3 ,LRBE3 ,IRBE2 ,
3709 B LRBE2 ,IADM ,JDIM ,DIAG_M,LT_M ,
3710 C F_DDL ,L_DDL ,ITASK ,V_W )
3711C----------------------
3712 CALL MY_BARRIER
3713C---------------------
3714 END DO !M = 1,M_VS1
3715C----------------------
3716C [m0]=[S]^t[S]
3717C---------------------
3718C CALL MAM_NM(F_DDL ,L_DDL ,NDDL ,NPCGPV ,PROJ_S ,PROJ_S ,P_M0 ,WDDL , ITASK )
3719C----------------------
3720C [k0]=[S]^t[T]
3721C---------------------
3722 CALL MAM_NM(F_DDL ,L_DDL ,NDDL ,M_VS1 ,PROJ_S ,
3723 . PROJ_T ,PROJ_K ,W_DDL ,ITASK )
3724C----------------------
3725 CALL MY_BARRIER
3726C---------------------
3727C ([k0]-lamda(i)[m0])*phi(i)=0; Rather ([m0]-lamda(i)[I])*phi(i)=0
3728C---------------------
3729 IF (ITASK==0) THEN
3730 JOBZ='v'
3731 UPLO='u'
3732 LW=3*M_VS1
3733 CALL DSYEV(JOBZ, UPLO, M_VS1,PROJ_K, M_VS1, W, WORK, LW, INFO)
3734C [S]<-[S][phi]
3735C---------------------
3736 CALL MAV_MM(NDDL ,M_VS1 ,PROJ_S ,PROJ_K ,ITASK )
3737C----------------------
3738C [T]<-[T][phi]
3739C---------------------
3740 CALL MAV_MM(NDDL ,M_VS1 ,PROJ_T ,PROJ_K ,ITASK )
3741C----------------------
3742C [LAMDA]^-1<- 1.0/lamda(i)
3743C---------------------
3744 DO I=1,M_VS1
3745 PROJ_LA_1(I)= ONE/SIGN(MAX(EM20,ABS(W(I))),W(I))
3746 ENDDO
3747.AND. if (iline==1IMPDEB>0) then
3748 IF (NCG_RUN==0) THEN
3749 write(iout,*)'info,ini eigenvalue=',info,(ONE/PROJ_LA_1(I),I=1,M_VS1)
3750 ELSEIF (iline==1) THEN
3751 write(iout,*)'info,upd eigenvalue=',info,(ONE/PROJ_LA_1(I),I=1,M_VS1)
3752 END IF
3753 endif
3754 END IF !(ITASK==0) THEN
3755#else
3756 CALL ARRET(5)
3757#endif
3758 RETURN
3759 END
3760!||====================================================================
3761!|| imp_updst ../engine/source/implicit/imp_pcg.F
3762!||--- called by ------------------------------------------------------
3763!|| imp_ppcgh ../engine/source/implicit/imp_pcg.F
3764!||--- calls -----------------------------------------------------
3765!|| mortho_gs ../engine/source/implicit/produt_v.F
3766!|| my_barrier ../engine/source/system/machine.F
3767!||--- uses -----------------------------------------------------
3768!|| groupdef_mod ../common_source/modules/groupdef_mod.F
3769!|| imp_pcg_proj ../engine/share/modules/impbufdef_mod.F
3770!|| intbufdef_mod ../common_source/modules/interfaces/intbufdef_mod.F90
3771!||====================================================================
3772 SUBROUTINE IMP_UPDST(
3773 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K ,
3774 2 LT_K ,ITOK ,IADI ,JDII ,LT_I ,
3775 3 ITASK ,W_DDL ,A ,AR ,VE ,
3776 4 D ,DR ,NDOF ,IPARI ,INTBUF_TAB,
3777 5 NUM_IMP,NS_IMP,NE_IMP,NSREM ,
3778 6 NSL ,NMONV ,IMONV ,MONVOL,IGRSURF ,
3779 7 VOLMON,FR_MV ,IBFV ,SKEW ,
3780 8 XFRAME ,IND_IMP,XI_C ,MS ,XE ,
3781 9 IRBE3 ,LRBE3 ,IRBE2 ,LRBE2 ,U ,
3782 A P ,Y ,IADM ,JDIM ,DIAG_M ,
3783 B LT_M ,F_DDL ,L_DDL )
3784C-----------------------------------------------
3785C M o d u l e s
3786C-----------------------------------------------
3787 USE IMP_PCG_PROJ
3788 USE INTBUFDEF_MOD
3789 USE GROUPDEF_MOD
3790C-----------------------------------------------
3791C I m p l i c i t T y p e s
3792C-----------------------------------------------
3793#include "implicit_f.inc"
3794C-----------------------------------------------
3795C C o m m o n B l o c k s
3796C-----------------------------------------------
3797#include "com04_c.inc"
3798#include "impl1_c.inc"
3799C-----------------------------------------------
3800C D u m m y A r g u m e n t s
3801C-----------------------------------------------
3802 INTEGER F_DDL,L_DDL,NDDL ,IADK(*) ,JDIK(*),
3803 . NDDLI ,ITOK(*) ,IADI(*),JDII(*),
3804 . ITASK,W_DDL(*),IBFV(*),IRBE3(*) ,LRBE3(*),
3805 . IRBE2(*) ,LRBE2(*),IADM(*) ,JDIM(*)
3806 INTEGER NDOF(*),NE_IMP(*),NSREM ,NSL,
3807 . IPARI(*) ,NUM_IMP(*),NS_IMP(*) ,IND_IMP(*)
3808 INTEGER NMONV,IMONV(*),MONVOL(*),FR_MV(*)
3809 my_real
3810 . DIAG_K(*), LT_K(*),LT_I(*) ,XI_C(*),U(*),P(*),Y(*)
3811 my_real
3812 . A(3,*),AR(3,*),VE(3,*),D(3,*),DR(3,*),XE(3,*),
3813 . MS(*),VOLMON(*),SKEW(*) ,XFRAME(*),
3814 . DIAG_M(*), LT_M(*)
3815
3816 TYPE(INTBUF_STRUCT_) INTBUF_TAB(*)
3817 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
3818C-----------------------------------------------
3819c FUNCTION: update S,T of Projection taking into account to preconditioner [M]
3820c
3821c Note:
3822c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
3823c
3824c TYPE NAME FUNCTION
3825c I NDDL,NNZ - dimension of [K] and number of non zero (strict triangular)
3826c I IADK,JDIK - indice arrays for compressed row(col.) format of [K]
3827c I DIAG_K(NDDL) - diagonal terms of [K]
3828c I LT_K(NNZ) - strict triangular of [K]
3829c O Proj_S(NDDL,M_VS) - [S] reduced small Eigenvectors
3830c O Proj_T(NDDL,M_VS) - [T] =[K][S]
3831c I W1,W2 - work arrays
3832C-----------------------------------------------
3833C L o c a l V a r i a b l e s
3834C-----------------------------------------------
3835 CHARACTER JOBZ, UPLO
3836 INTEGER I,J,IT,IP,NLIM,ND,IUPD,IPRI,IERROR,NNZI,M,
3837 . INFO,LW,M_VS1,INORM,NP
3838 my_real
3839 . WORK(3*M_VS+3),W(M_VS+1)
3840C------M_VS input one -NPCGPV: activated ---------------
3841 IF (NPCGPV == 0) RETURN
3842C----------------------
3843 CALL MY_BARRIER
3844C---------------------
3845 M_VS1 = NPCGPV+ 1
3846C------add previous solution U ; default : IPRO_S0=2
3847C---IPRO_S0=1 : alertory updated w/ 1 vector x
3848C---IPRO_S0=2 : ini S w/ first p vectors updated w/ 2 vectors p,x
3849C---IPRO_S0=3 : ini S w/ first p vectors updated w/ 2 vectors p,x (p<-RR_max)
3850C---IPRO_S0=4 : ini S w/ first p vectors updated w/ 2 vectors v,x (v update each ite)
3851c IF (IPRO_S0 > 2.AND.IPRO_S0 < 4.AND.NCG_RUN==0) THEN
3852c CALL IMP_INIS(NDDL ,F_DDL ,L_DDL , 1,NPCGPV-1 ,PROJ_S )
3853C----------------------
3854c CALL MY_BARRIER
3855C---------------------
3856c CALL MORTHO_GS(F_DDL ,L_DDL ,NDDL ,1 ,NPCGPV-1 ,
3857c . PROJ_S ,W_DDL ,ITASK )
3858C----------------------
3859c CALL MY_BARRIER
3860C---------------------
3861c END IF !(IPRO_S0 > 1.AND.NCG_RUN==0) THEN
3862C
3863 IF (IPRO_S0==2) THEN
3864 DO I=F_DDL,L_DDL
3865 PROJ_S(I,NPCGPV)=U(I)
3866 PROJ_S(I,M_VS1)=P(I)
3867 ENDDO
3868.OR. ELSEIF (IPRO_S0==3IPRO_S0==4) THEN
3869 DO I=F_DDL,L_DDL
3870 PROJ_S(I,NPCGPV)=U(I)
3871 PROJ_S(I,M_VS1)=PROJ_V(I)
3872 ENDDO
3873 ELSE
3874 DO I=F_DDL,L_DDL
3875 PROJ_S(I,M_VS1)=U(I)
3876 ENDDO
3877 END IF !(MOD(IPRO_S0,2)==0) THEN
3878C----------------------
3879 CALL MY_BARRIER
3880C---------------------
3881 CALL MORTHO_GS(F_DDL ,L_DDL ,NDDL ,NPCGPV ,M_VS1 ,
3882 . PROJ_S ,W_DDL ,ITASK )
3883C-------will do (T=[K]S ...) in IMP_INIST
3884C
3885 RETURN
3886 END
3887!||====================================================================
3888!|| imp_inisi ../engine/source/implicit/imp_pcg.F
3889!||--- called by ------------------------------------------------------
3890!|| imp_ppcgh ../engine/source/implicit/imp_pcg.F
3891!||--- calls -----------------------------------------------------
3892!|| imp_inis ../engine/source/implicit/imp_pcg.F
3893!|| mmav_lth ../engine/source/implicit/produt_v.F
3894!|| mortho_gs ../engine/source/implicit/produt_v.F
3895!|| my_barrier ../engine/source/system/machine.F
3896!|| produt_h ../engine/source/implicit/produt_v.F
3897!||--- uses -----------------------------------------------------
3898!|| groupdef_mod ../common_source/modules/groupdef_mod.F
3899!|| imp_pcg_proj ../engine/share/modules/impbufdef_mod.F
3900!|| intbufdef_mod ../common_source/modules/interfaces/intbufdef_mod.F90
3901!||====================================================================
3902 SUBROUTINE IMP_INISI(
3903 1 NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
3904 2 NDDLI ,ITOK ,IADI ,JDII ,LT_I ,
3905 3 NNZM ,IADM ,JDIM ,DIAG_M ,LT_M ,
3906 4 R ,P ,Y ,
3907 8 A ,AR ,VE ,MS ,XE ,
3908 9 D ,DR ,NDOF ,IPARI ,INTBUF_TAB ,
3909 A NUM_IMP,NS_IMP,NE_IMP,NSREM ,
3910 B NSL ,NMONV ,IMONV ,MONVOL,IGRSURF ,
3911 C VOLMON,FR_MV ,IBFV ,SKEW ,
3912 D XFRAME ,Z ,IDDL ,IKC ,INLOC ,
3913 G IND_IMP,XI_C ,IRBE3 ,LRBE3 ,IRBE2 ,
3914 H LRBE2 ,ITASK ,W_DDL ,F_DDL ,L_DDL ,
3915 I SQ_DIAG_M)
3916C-----------------------------------------------
3917C M o d u l e s
3918C-----------------------------------------------
3919 USE IMP_PCG_PROJ
3920 USE INTBUFDEF_MOD
3921 USE GROUPDEF_MOD
3922C-----------------------------------------------
3923C I m p l i c i t T y p e s
3924C-----------------------------------------------
3925#include "implicit_f.inc"
3926C-----------------------------------------------
3927C C o m m o n B l o c k s
3928C-----------------------------------------------
3929#include "com04_c.inc"
3930#include "impl1_c.inc"
3931C-----------------------------------------------
3932C D u m m y A r g u m e n t s
3933C-----------------------------------------------
3934 INTEGER NDDL ,IADK(*) ,JDIK(*),
3935 . NDDLI ,ITOK(*) ,IADI(*),JDII(*),
3936 . NNZM ,IADM(*),JDIM(*),ITASK,
3937 . W_DDL(*),IBFV(*),IRBE3(*) ,LRBE3(*),
3938 . IRBE2(*) ,LRBE2(*),F_DDL ,L_DDL
3939 INTEGER NDOF(*),NE_IMP(*),NSREM ,NSL,
3940 . IPARI(*) ,NUM_IMP(*),NS_IMP(*) ,IND_IMP(*)
3941 INTEGER NMONV,IMONV(*),MONVOL(*),FR_MV(*)
3942 INTEGER IDDL(*), IKC(*), INLOC(*)
3943 my_real
3944 . DIAG_K(*), LT_K(*) ,DIAG_M(*),LT_M(*) ,
3945 . P(*) ,R(*) ,Y(*),LT_I(*) ,XI_C(*),Z(*),SQ_DIAG_M(*)
3946 my_real
3947 . A(3,*),AR(3,*),VE(3,*),D(3,*),DR(3,*),XE(3,*),
3948 . MS(*),VOLMON(*),SKEW(*) ,XFRAME(*)
3949
3950 TYPE(INTBUF_STRUCT_) INTBUF_TAB(*)
3951 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
3952C-----------------------------------------------
3953c FUNCTION: initialization of S,by M_VS PCG iterations
3954c
3955c Note:
3956c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
3957c
3958c TYPE NAME FUNCTION
3959c I NDDL,NNZ - dimension of [K] and number of non zero (strict triangular)
3960c I IADK,JDIK - indice arrays for compressed row(col.) format of [K]
3961c I DIAG_K(NDDL) - diagonal terms of [K]
3962c I LT_K(NNZ) - strict triangular of [K]
3963c O Proj_S(NDDL,M_VS) - [S] reduced small Eigenvectors
3964c IO P,Y,X,R(NDDL) -work vectors
3965C-----------------------------------------------
3966C L o c a l V a r i a b l e s
3967C-----------------------------------------------
3968 INTEGER I,J,IT,IP,NLIM,ND,IUPD,IPRI,IERROR,NNZI,M,
3969 . F_DDLI,L_DDLI,INFO,LW,INORM,NPV,ITP
3970 my_real
3971 . ALPHA,BETA, G0,G1,S
3972 my_real
3973 . ALEAT
3974 EXTERNAL ALEAT
3975C---------------------
3976.OR. IF (NPCGPV >=0NCG_RUN>0) RETURN
3977C
3978 NPV=MIN(NDDL-1,M_VS)
3979C---IPRO_S0=1 : alertory updated w/ 1 vector x
3980C---IPRO_S0=2 : ini S w/ first p vectors updated w/ 2 vectors p,x
3981C---IPRO_S0=3 : alertory ini S vectors updated w/ 2 vectors p,x (p<-RR_max)
3982C---------------------
3983 IUPD = 0
3984 IF (IPRO_S0 > 1) THEN
3985C------seems more stable-----
3986 DO I=F_DDL,L_DDL
3987 R(I) = ALEAT()
3988 ENDDO
3989C----------------------
3990 CALL MY_BARRIER
3991C---------------------
3992 DO I = F_DDL,L_DDL
3993 P(I) = R(I)
3994 ENDDO
3995 CALL PRODUT_H(F_DDL ,L_DDL ,R ,R ,W_DDL , G0 ,ITASK )
3996c
3997 DO ITP=1,NPV
3998 CALL MMAV_LTH(
3999 1 NDDL ,NDDLI ,IADK ,JDIK ,DIAG_K,
4000 2 LT_K ,IADI ,JDII ,ITOK ,LT_I ,
4001 3 P ,Y ,A ,AR ,
4002 5 VE ,MS ,XE ,D ,DR ,
4003 6 NDOF ,IPARI ,INTBUF_TAB ,NUM_IMP,
4004 7 NS_IMP,NE_IMP,NSREM ,NSL ,IBFV ,
4005 8 SKEW ,XFRAME,MONVOL,VOLMON,IGRSURF ,
4006 9 FR_MV,NMONV ,IMONV ,IND_IMP,
4007 A XI_C ,IUPD ,IRBE3 ,LRBE3 ,IRBE2 ,
4008 B LRBE2 ,IADM ,JDIM ,SQ_DIAG_M,LT_M ,
4009 C F_DDL ,L_DDL ,ITASK ,Z )
4010C----------------------
4011 CALL MY_BARRIER
4012C---------------------
4013C---------------------
4014 CALL PRODUT_H(F_DDL ,L_DDL ,P ,Y ,W_DDL , S ,ITASK )
4015C---------------------
4016 ALPHA = G0/S
4017C----------------------
4018 CALL MY_BARRIER
4019C---------------------
4020 DO I=F_DDL,L_DDL
4021 R(I) = R(I) - ALPHA*Y(I)
4022 ENDDO
4023C----------------------
4024 CALL MY_BARRIER
4025C---------------------
4026 CALL PRODUT_H( F_DDL ,L_DDL ,R ,R ,W_DDL , G1 ,ITASK )
4027C---------------------
4028 BETA=G1/G0
4029 G0 = G1
4030C----------------------
4031 CALL MY_BARRIER
4032C---------------------
4033 DO I=F_DDL,L_DDL
4034 P(I) = R(I) + BETA*P(I)
4035 ENDDO
4036C----------------------
4037 CALL MY_BARRIER
4038C---------------------
4039 DO I=F_DDL,L_DDL
4040 PROJ_S(I,ITP) = P(I)
4041 ENDDO
4042C----------------------
4043 CALL MY_BARRIER
4044C---------------------
4045 END DO !ITP=1,NPC
4046 ELSE
4047 CALL IMP_INIS(NDDL ,F_DDL ,L_DDL ,1 ,NPV ,PROJ_S )
4048C----------------------
4049 CALL MY_BARRIER
4050C---------------------
4051 END IF !(IPRO_S0 > 1) THEN
4052C--------matrix GS orthonalization and normalized
4053 CALL MORTHO_GS(F_DDL ,L_DDL ,NDDL ,1 ,NPV ,
4054 . PROJ_S ,W_DDL ,ITASK )
4055C
4056 RETURN
4057 END
my_real
#define my_real
Definition cppsort.cpp:32
the
end diagonal values have been computed in the(sparse) matrix id.SOL
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
alpha
#define alpha
Definition eval.h:35
imp_isave
subroutine imp_isave(nddl, x, r, diag_k, diag_m, nnz, lt_k, lt_k0, lt_m, lt_m0, iadk, jdik, iadk0, jdik0, iadm, jdim, iadm0, jdim0, nnzm, tols, nlim, itol, eps_m, iprec)
Definition imp_pcg.F:69
imp_rsave
subroutine imp_rsave(nddl, x, r)
Definition imp_pcg.F:170
nlim_mix
subroutine nlim_mix(nddl, nddli, nlim)
Definition imp_pcg.F:2451
imp_pcgh
subroutine imp_pcgh(iprec, nddl, nnz, iadk, jdik, diag_k, lt_k, nddli, itok, iadi, jdii, lt_i, nnzm, iadm, jdim, diag_m, lt_m, x, r, itol, tol, p, z, y, itask, iprint, n_max, eps_m, 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, graphe, iad_elem, fr_elem, itab, insolv, itn, fac_k, ipiv_k, nk, mumps_par, cddlp, isolv, idsc, iddl, ikc, inloc, ind_imp, xi_c, r0, nddli_g, intp_c, irbe3, lrbe3, irbe2, lrbe2)
Definition imp_pcg.F:245
crit_stop
integer function crit_stop(it, r2, it_max, tol)
Definition imp_pcg.F:32
imp_isave2
subroutine imp_isave2(nddl, x, diag_k, nnz, lt_k, iadk, jdik, lt_k0, iadk0, jdik0)
Definition imp_pcg.F:116
spmd_sum_s
subroutine spmd_sum_s(s)
Definition imp_spmd.F:1037
max
#define max(a, b)
Definition macros.h:21
dsgraph_mod
Definition dsgraph_mod.F:31
groupdef_mod
Definition groupdef_mod.F:662
imp_frk
Definition impbufdef_mod.F:740
imp_frk::nd_fr
integer, dimension(:), allocatable nd_fr
Definition impbufdef_mod.F:754
imp_intm
Definition imp_intm.F:118
imp_intm::jdi_si
integer, dimension(:), allocatable jdi_si
Definition imp_intm.F:174
imp_intm::iddl_sl
integer, dimension(:), allocatable iddl_sl
Definition imp_intm.F:178
imp_intm::iad_ss
integer, dimension(:), allocatable iad_ss
Definition imp_intm.F:175
imp_intm::iad_si
integer, dimension(:), allocatable iad_si
Definition imp_intm.F:174
imp_intm::jdi_sl
integer, dimension(:), allocatable jdi_sl
Definition imp_intm.F:175
imp_intm::iad_sl
integer, dimension(:), allocatable iad_sl
Definition imp_intm.F:145
imp_intm::nddl_si
integer nddl_si
Definition imp_intm.F:173
imp_intm::nddl_sl
integer nddl_sl
Definition imp_intm.F:173
imp_intm::iad_srem
integer, dimension(:), allocatable iad_srem
Definition imp_intm.F:145
imp_workh
Definition impbufdef_mod.F:571
imp_workh::jdik0
integer, dimension(:), allocatable jdik0
Definition impbufdef_mod.F:584
imp_workh::iadi0
integer, dimension(:), allocatable iadi0
Definition impbufdef_mod.F:588
imp_workh::jdim0
integer, dimension(:), allocatable jdim0
Definition impbufdef_mod.F:584
imp_workh::iadk0
integer, dimension(:), allocatable iadk0
Definition impbufdef_mod.F:584
imp_workh::iadm0
integer, dimension(:), allocatable iadm0
Definition impbufdef_mod.F:584
message_mod
Definition message_mod.F:1249
mav_lth
subroutine mav_lth(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, f_ddl, l_ddl, itask)
Definition produt_v.F:1380
produt_v
subroutine produt_v(nddl, x, y, r)
Definition produt_v.F:33
startime
subroutine startime(event, itask)
Definition timer.F:93
dsgraph_mod::prgraph
Definition dsgraph_mod.F:45