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 special routine to save inputs
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 IF(ipri/=0.AND.itask==0)THEN
564 rr = sqrt(r2)
565 IF(ipri<0) WRITE(istdo,1000)it,rr
566 WRITE(iout,1000)it,rr
567 ENDIF
568C-------to be consistent with Lanczos 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 IF(nddli_g>0.AND.intp_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 IF(iupd0>0.AND.it>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 IF(ipri/=0.AND.itask==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 IF((iupd>0.AND.it==nlim).OR.
736 . (iupd==0.AND.istop/=1.AND.iupd0>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 FOR 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 Error code
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 If SPMD think to do before d_y = d_z*w for weights
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 If SPMD think to do before d_z = d_y*d_w for weights
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 If SPMD think to do before d_z = d_r*w for weights
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-------to be consistent with Lanczos 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 If SPMD think to do before d_y = d_z*d_w for weights
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 If SPMD think to do before d_y = d_r*w for weights
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 If SPMD think about doing before d_y = d_y*w for weights
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 If SPMD think to do before d_y = d_z*w for weights
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 If SPMD think to do before d_y = d_r*w for weights
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 provisional code
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 If SPMD think to do before d_y = d_z*w for weights
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 If SPMD think to do before d_z = d_y*d_w for weights
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 If SPMD think to do before d_z = d_r*w for weights
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-------to be consistent with Lanczos 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 If SPMD think to do before d_y = d_z*d_w for weights
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 If SPMD think to do before d_y = d_r*w for weights
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 If SPMD think about doing before d_y = d_y*w for weights
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 If SPMD think to do before d_y = d_z*w for weights
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 If SPMD think to do before d_y = d_r*w for weights
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 End Provisional Code
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 End SPECIFIC GPU code
2307C special routine to save results
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-------to be consistent with Lanczos 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
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
dsgraph_mod
Definition dsgraph_mod.F:31
groupdef_mod
Definition groupdef_mod.F:668
imp_frk
Definition impbufdef_mod.F:740
imp_frk::ifr2k
integer, dimension(:), allocatable ifr2k
Definition impbufdef_mod.F:758
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:1257
prec_solvh
subroutine prec_solvh(iprec, itask, graphe, iad_elem, fr_elem, diag_k, lt_k, iadk, jdik, itab, iprint, insolv, it, fac_k, ipiv_k, nk, mumps_par, cddlp, isolv, idsc, iddl, ikc, inloc, ndof, nddl, nnz, iadm, jdim, diag_m, lt_m, v, z, f_ddl, l_ddl)
Definition prec_solv.F:590
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_h
subroutine produt_h(f_ddl, l_ddl, x, y, w, r, itask)
Definition produt_v.F:1533
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:895
arret
subroutine arret(nn)
Definition arret.F:86
my_barrier
subroutine my_barrier
Definition machine.F:31
dsgraph_mod::prgraph
Definition dsgraph_mod.F:45