OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
imp_fsa_inv.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| sp_stat0 ../engine/source/implicit/imp_fsa_inv.F
25!||--- called by ------------------------------------------------------
26!|| imp_fsa_inv2 ../engine/source/implicit/imp_fsa_inv.F
27!|| imp_fsa_invh ../engine/source/implicit/imp_fsa_inv.F
28!|| imp_fsa_invh2 ../engine/source/implicit/imp_fsa_inv.F
29!|| imp_fsa_invp ../engine/source/implicit/imp_fsa_inv.F
30!|| imp_fsa_invp2 ../engine/source/implicit/imp_fsa_inv.F
31!|| m_lnz ../engine/source/implicit/imp_solv.F
32!|| sms_fsa_invh ../engine/source/ams/sms_fsa_inv.F
33!||====================================================================
34 SUBROUTINE sp_stat0(IL ,IADK ,JDIK ,NC ,JM )
35C-----------------------------------------------
36C I m p l i c i t T y p e s
37C-----------------------------------------------
38#include "implicit_f.inc"
39C-----------------------------------------------
40C D u m m y A r g u m e n t s
41C-----------------------------------------------
42 INTEGER IL, IADK(*) ,JDIK(*),NC,JM(*)
43C REAL
44C-----------------------------------------------
45C L o c a l V a r i a b l e s
46C-----------------------------------------------
47 INTEGER I,J,K
48C-----------------------------------------------
49 nc=0
50 DO k =iadk(il),iadk(il+1)-1
51 nc=nc+1
52 jm(nc)=jdik(k)
53 ENDDO
54 nc=nc+1
55 jm(nc)=il
56C--------------------------------------------
57 RETURN
58 END
59!||====================================================================
60!|| sp_static ../engine/source/implicit/imp_fsa_inv.F
61!||--- called by ------------------------------------------------------
62!|| imp_pc_inv ../engine/source/implicit/imp_pc_inv.F
63!||--- calls -----------------------------------------------------
64!|| arret ../engine/source/system/arret.F
65!|| sp_a2 ../engine/source/implicit/imp_fsa_inv.F
66!||====================================================================
67 SUBROUTINE sp_static(NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
68 . IADM ,JDIM ,NNZM ,NC ,JM ,
69 . MAXC ,PSI ,IP )
70C-----------------------------------------------
71C I m p l i c i t T y p e s
72C-----------------------------------------------
73#include "implicit_f.inc"
74C-----------------------------------------------
75C D u m m y A r g u m e n t s
76C-----------------------------------------------
77C--IP: 0 standard;Apaire: with pre-filter;2.3 Level 2;> 4 FSAI (IP-4: 6.7: Level 2)
78 INTEGER NDDL ,MAXC,IADK(*) ,JDIK(*)
79 INTEGER NNZM,IADM(*) ,JDIM(*),NC(*),JM(MAXC,*),IP
80C REAL
81 my_real
82 . lt_k(*),diag_k(*),psi
83C-----------------------------------------------
84C L o c a l V a r i a b l e s
85C-----------------------------------------------
86 INTEGER I,J,N,K,I1,IFSAI,IOPT
87 my_real
88 . psr
89C-----------------------------------------------
90 IF (ip>=4) THEN
91 ifsai=1
92 iopt=ip-4
93 ELSE
94 ifsai=0
95 iopt=ip
96 ENDIF
97C------define sparse static pattern:----
98 nnzm = 0
99 iadm(1)=1
100C-------symmetry firstly---
101 IF (mod(iopt,2)>0) THEN
102C-------pre-filtrage----
103 DO i=1,nddl
104 DO k =iadk(i),iadk(i+1)-1
105 j=jdik(k)
106 psr = psi*sqrt(diag_k(i)*diag_k(j))
107 IF (abs(lt_k(k))>=psr) THEN
108 nnzm = nnzm+1
109 jdim(nnzm)=j
110 ENDIF
111 ENDDO
112 iadm(i+1)=nnzm+1
113 ENDDO
114 ELSE
115 DO i=2,nddl+1
116 iadm(i)=iadk(i)
117 ENDDO
118 nnzm= iadk(nddl+1)-1
119 DO i=1,nnzm
120 jdim(i)=jdik(i)
121 ENDDO
122 ENDIF
123C
124 DO i=1,nddl
125 nc(i)=0
126 ENDDO
127C
128 IF (iopt>1) THEN
129 DO i=1,nddl
130C------------with diag----
131 nc(i) = nc(i)+1
132 jm(nc(i),i)=i
133 DO k =iadm(i),iadm(i+1)-1
134 j=jdim(k)
135 nc(j) = nc(j)+1
136 jm(nc(j),j)=i
137 nc(i) = nc(i)+1
138 jm(nc(i),i)=j
139 ENDDO
140 ENDDO
141 CALL sp_a2(nddl,nc,jm,maxc,ifsai)
142 ELSE
143 IF (ifsai==1) THEN
144 DO i=1,nddl
145 nc(i) = nc(i)+1
146 jm(nc(i),i)=i
147 DO k =iadm(i),iadm(i+1)-1
148 j=jdim(k)
149 nc(j) = nc(j)+1
150 jm(nc(j),j)=i
151 ENDDO
152 IF (nc(i)>maxc) THEN
153 WRITE(*,*)'N>MAXB',nc(i),maxc,i
154 CALL arret(2)
155 ENDIF
156 ENDDO
157 ELSE
158 DO i=1,nddl
159 nc(i) = nc(i)+1
160 jm(nc(i),i)=i
161 DO k =iadm(i),iadm(i+1)-1
162 j=jdim(k)
163 nc(j) = nc(j)+1
164 jm(nc(j),j)=i
165 nc(i) = nc(i)+1
166 jm(nc(i),i)=j
167 ENDDO
168 IF (nc(i)>maxc) THEN
169 WRITE(*,*)'N>MAXB',nc(i),maxc,i
170 CALL arret(2)
171 ENDIF
172 ENDDO
173 ENDIF
174 ENDIF
175C--------------------------------------------
176 RETURN
177 END
178!||====================================================================
179!|| sp_a2 ../engine/source/implicit/imp_fsa_inv.f
180!||--- called by ------------------------------------------------------
181!|| sp_static ../engine/source/implicit/imp_fsa_inv.F
182!||--- calls -----------------------------------------------------
183!|| arret ../engine/source/system/arret.F
184!|| intab2 ../engine/source/implicit/imp_fsa_inv.F
185!||====================================================================
186 SUBROUTINE sp_a2(NDDL,NC,JM,MAXC,IFSAI)
187C-----------------------------------------------
188C I m p l i c i t T y p e s
189C-----------------------------------------------
190#include "implicit_f.inc"
191C-----------------------------------------------
192C D u m m y A r g u m e n t s
193C-----------------------------------------------
194 INTEGER NDDL,NC(*),MAXC,IFSAI
195 INTEGER JM(MAXC,*)
196C-----------------------------------------------
197C External function
198C-----------------------------------------------
199 INTEGER INTAB2
200 EXTERNAL INTAB2
201C REAL
202C-----------------------------------------------
203C L o c a l V a r i a b l e s
204C-----------------------------------------------
205 INTEGER I,J,NN(NDDL),JN(MAXC,NDDL)
206C-----------------------------------------------
207 DO i=1,nddl
208 nn(i)=nc(i)
209 nc(i)=0
210 DO j=1,nn(i)
211 jn(j,i)=jm(j,i)
212 ENDDO
213 ENDDO
214 IF (ifsai==1) THEN
215 DO i=1,nddl
216 nc(i) = nc(i)+1
217 jm(nc(i),i)=i
218 DO j=i+1,nddl
219 IF(intab2(nn(i),jn(1,i),nn(j),jn(1,j))>0) THEN
220 nc(j) = nc(j)+1
221 jm(nc(j),j)=i
222 IF (nc(j)>maxc) THEN
223 WRITE(*,*)'N>MAXB',nc(j),maxc,j
224 CALL arret(2)
225 ENDIF
226 ENDIF
227 ENDDO
228 ENDDO
229 ELSE
230 DO i=1,nddl
231 nc(i) = nc(i)+1
232 jm(nc(i),i)=i
233 DO j=i+1,nddl
234 IF(intab2(nn(i),jn(1,i),nn(j),jn(1,j))>0) THEN
235 nc(i) = nc(i)+1
236 jm(nc(i),i)=j
237 nc(j) = nc(j)+1
238 jm(nc(j),j)=i
239 IF (nc(i)>maxc) THEN
240 WRITE(*,*)'N>MAXB',nc(i),maxc,i
241 CALL arret(2)
242 ENDIF
243 ENDIF
244 ENDDO
245 ENDDO
246 ENDIF
247C--------------------------------------------
248 RETURN
249 END
250C-------------resol A(N,N).MJ=B ---B=MJ(input)-
251!||====================================================================
252!|| imp_fsai ../engine/source/implicit/imp_fsa_inv.F
253!||--- called by ------------------------------------------------------
254!|| fsa_solv ../engine/source/implicit/imp_fsa_inv.F
255!|| imp_fsa_inv2 ../engine/source/implicit/imp_fsa_inv.F
256!|| imp_fsa_invh ../engine/source/implicit/imp_fsa_inv.F
257!|| imp_fsa_invh2 ../engine/source/implicit/imp_fsa_inv.F
258!|| imp_fsa_invp ../engine/source/implicit/imp_fsa_inv.F
259!|| imp_fsa_invp2 ../engine/source/implicit/imp_fsa_inv.F
260!|| sms_fsa_invh ../engine/source/ams/sms_fsa_inv.f
261!||--- calls -----------------------------------------------------
262!|| imp_fac_icj ../engine/source/implicit/imp_fac_ic.f
263!|| prec0_solv ../engine/source/implicit/prec_solv.F
264!||====================================================================
265 SUBROUTINE imp_fsai(N ,IADA ,JDIA ,DIAG_A ,LT_A,
266 . MAXA ,MJ )
267C-----------------------------------------------
268C I m p l i c i t T y p e s
269C-----------------------------------------------
270#include "implicit_f.inc"
271C-----------------------------------------------
272C D u m m y A r g u m e n t s
273C-----------------------------------------------
274 INTEGER N ,IADA(*) ,JDIA(*),MAXA
275C REAL
276 my_real
277 . diag_a(*),lt_a(*),mj(*)
278C-----------------------------------------------
279C L o c a l V a r i a b l e s
280C-----------------------------------------------
281 INTEGER I,IADL(N+1),JDIL(MAXA),NNZL,NNE,IWA1(N)
282 my_real
283 . LT_L(MAXA),WA1(N),DIAG_L(N)
284C-----------------------------------------------
285 CALL imp_fac_icj(
286 1 n ,maxa ,iada ,jdia ,diag_a ,
287 2 lt_a ,iadl ,jdil ,diag_l,lt_l ,
288 3 zero ,nnzl ,maxa ,iwa1 ,wa1 ,
289 4 nne )
290C
291#ifdef MUMPS5
292 CALL prec0_solv(n ,nnzl ,iadl ,jdil ,diag_l ,
293 1 lt_l ,mj ,wa1 )
294#endif
295 DO i=1,n
296 mj(i)=wa1(i)
297 ENDDO
298C--------------------------------------------
299 RETURN
300 END
301C-------------set submatrix A(N,N) Format m.c.c.s. for FSAI ----
302!||====================================================================
303!|| get_subs0 ../engine/source/implicit/imp_fsa_inv.F
304!||--- called by ------------------------------------------------------
305!|| imp_fsa_inv2 ../engine/source/implicit/imp_fsa_inv.F
306!||--- calls -----------------------------------------------------
307!|| ind_lt2ln ../engine/source/implicit/imp_fsa_inv.F
308!|| intab0 ../engine/source/implicit/imp_fsa_inv.F
309!||====================================================================
310 SUBROUTINE get_subs0(NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
311 . NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
312 . JM ,MAXA )
313C-----------------------------------------------
314C I m p l i c i t T y p e s
315C-----------------------------------------------
316#include "implicit_f.inc"
317C-----------------------------------------------
318C D u m m y A r g u m e n t s
319C-----------------------------------------------
320 INTEGER NDDL ,IADK(*) ,JDIK(*),IADA(*) ,JDIA(*)
321 INTEGER NC ,JM(*),MAXA
322C REAL
323 my_real
324 . LT_K(*),DIAG_K(*),LT_A(*),DIAG_A(*)
325C-----------------------------------------------
326C External function
327C-----------------------------------------------
328 INTEGER INTAB0
329 EXTERNAL INTAB0
330C-----------------------------------------------
331C L o c a l V a r i a b l e s
332C-----------------------------------------------
333 INTEGER I,J,K,JJ,NNZA,N
334C--------------------------------------------
335 nnza=0
336 iada(1)=1
337 DO i=1,nc
338 j=jm(i)
339 diag_a(i)=diag_k(j)
340 DO k=iadk(j),iadk(j+1)-1
341 jj=jdik(k)
342 n=intab0(nc,jm,jj)
343 IF (n>0) THEN
344 nnza=nnza+1
345 jdia(nnza)=n
346 lt_a(nnza)=lt_k(k)
347 ENDIF
348 ENDDO
349 iada(i+1)=nnza+1
350 ENDDO
351 CALL ind_lt2ln(nc,iada ,jdia ,lt_a, nnza )
352C
353 RETURN
354 END
355C----------version spmd---set submatrix A(N,N) Format m.c.c.s. for FSAI ----
356!||====================================================================
357!|| get_subsp ../engine/source/implicit/imp_fsa_inv.F
358!||--- called by ------------------------------------------------------
359!|| fsa_solv ../engine/source/implicit/imp_fsa_inv.F
360!|| imp_fsa_invh ../engine/source/implicit/imp_fsa_inv.F
361!|| imp_fsa_invh2 ../engine/source/implicit/imp_fsa_inv.F
362!|| imp_fsa_invp ../engine/source/implicit/imp_fsa_inv.F
363!|| imp_fsa_invp2 ../engine/source/implicit/imp_fsa_inv.F
364!||--- calls -----------------------------------------------------
365!|| ind_lt2ln ../engine/source/implicit/imp_fsa_inv.F
366!|| intab0 ../engine/source/implicit/imp_fsa_inv.F
367!||====================================================================
368 SUBROUTINE get_subsp(NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
369 . NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
370 . JM ,MAXA ,IDLFT0,IDLFT1 ,DIAG_C,
371 . LT_C ,DIAG_M ,LT_M )
372C-----------------------------------------------
373C I m p l i c i t T y p e s
374C-----------------------------------------------
375#include "implicit_f.inc"
376C-----------------------------------------------
377C D u m m y A r g u m e n t s
378C-----------------------------------------------
379 INTEGER NDDL ,IADK(*) ,JDIK(*),IADA(*) ,JDIA(*)
380 INTEGER NC ,JM(*),MAXA,IDLFT0,IDLFT1
381C REAL
382 my_real
383 . LT_K(*),DIAG_K(*),LT_A(*),DIAG_A(*),LT_C(*),DIAG_C(*),
384 . DIAG_M(*) ,LT_M(*)
385C-----------------------------------------------
386C External function
387C-----------------------------------------------
388 INTEGER INTAB0
389 EXTERNAL INTAB0
390C-----------------------------------------------
391C L o c a l V a r i a b l e s
392C-----------------------------------------------
393 INTEGER I,J,K,JJ,NNZA,N,K0
394C--------------------------------------------
395 NNZA=0
396 iada(1)=1
397 k0=iadk(idlft1+1)-1
398#include "vectorize.inc"
399 DO i=1,nc
400 j=jm(i)
401 IF (j<=idlft0) THEN
402 diag_a(i)=diag_m(j)
403 DO k=iadk(j),iadk(j+1)-1
404 jj=jdik(k)
405 n=intab0(nc,jm,jj)
406 IF (n>0) THEN
407 nnza=nnza+1
408 jdia(nnza)=n
409 lt_a(nnza)=lt_m(k)
410 ENDIF
411 ENDDO
412 ELSEIF (j>idlft1) THEN
413 diag_a(i)=diag_c(j-idlft1)
414 DO k=iadk(j),iadk(j+1)-1
415 jj=jdik(k)
416C IF (JJ>IDLFT1) THEN
417 n=intab0(nc,jm,jj)
418 IF (n>0) THEN
419 nnza=nnza+1
420 jdia(nnza)=n
421 lt_a(nnza)=lt_c(k-k0)
422 ENDIF
423C ENDIF
424 ENDDO
425 ELSE
426 diag_a(i)=diag_k(j)
427 DO k=iadk(j),iadk(j+1)-1
428 jj=jdik(k)
429 n=intab0(nc,jm,jj)
430 IF (n>0) THEN
431 nnza=nnza+1
432 jdia(nnza)=n
433 lt_a(nnza)=lt_k(k)
434 ENDIF
435 ENDDO
436 ENDIF
437 iada(i+1)=nnza+1
438 ENDDO
439C CALL IND_LT2L(NC,IADA ,JDIA ,LT_A, NNZA )
440 CALL ind_lt2ln(nc,iada ,jdia ,lt_a, nnza )
441C
442 RETURN
443 END
444C-------------set submatrix A(N,N) Format m.c.r.s. for FSAI ----
445!||====================================================================
446!|| get_subsa ../engine/source/implicit/imp_fsa_inv.F
447!||--- calls -----------------------------------------------------
448!|| intab0 ../engine/source/implicit/imp_fsa_inv.F
449!||====================================================================
450 SUBROUTINE get_subsa(NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
451 . NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
452 . JM )
453C-----------------------------------------------
454C I m p l i c i t T y p e s
455C-----------------------------------------------
456#include "implicit_f.inc"
457C-----------------------------------------------
458C D u m m y A r g u m e n t s
459C-----------------------------------------------
460 INTEGER NDDL ,IADK(*) ,JDIK(*),IADA(*) ,JDIA(*)
461 INTEGER NC ,JM(*)
462C REAL
463 my_real
464 . LT_K(*),DIAG_K(*),LT_A(*),DIAG_A(*)
465C-----------------------------------------------
466C External function
467C-----------------------------------------------
468 INTEGER INTAB0
469 EXTERNAL INTAB0
470C-----------------------------------------------
471C L o c a l V a r i a b l e s
472C-----------------------------------------------
473 INTEGER I,J,K,JJ,NNZA,N
474C--------------------------------------------
475 NNZA=0
476 iada(1)=1
477 DO i=1,nc
478 j=jm(i)
479 diag_a(i)=diag_k(j)
480 DO k=iadk(j),iadk(j+1)-1
481 jj=jdik(k)
482 n=intab0(nc,jm,jj)
483 IF (n>0) THEN
484 nnza=nnza+1
485 jdia(nnza)=n
486 lt_a(nnza)=lt_k(k)
487 ENDIF
488 ENDDO
489 iada(i+1)=nnza+1
490 ENDDO
491C
492 RETURN
493 END
494!||====================================================================
495!|| intab0 ../engine/source/implicit/imp_fsa_inv.F
496!||--- called by ------------------------------------------------------
497!|| dim_subnz ../engine/source/implicit/imp_solv.F
498!|| get_kijs ../engine/source/implicit/imp_pc_inv.F
499!|| get_subs0 ../engine/source/implicit/imp_fsa_inv.F
500!|| get_subsa ../engine/source/implicit/imp_fsa_inv.F
501!|| get_subsn ../engine/source/implicit/imp_fsa_inv.F
502!|| get_subsp ../engine/source/implicit/imp_fsa_inv.F
503!|| get_subsp_sms ../engine/source/ams/sms_fsa_inv.F
504!|| intab2 ../engine/source/implicit/imp_fsa_inv.F
505!|| spc_fr_k ../engine/source/mpi/implicit/imp_fri.F
506!|| upd_fr_k ../engine/source/mpi/implicit/imp_fri.F
507!||====================================================================
508 INTEGER FUNCTION intab0(NIC,IC,N)
509C----6---------------------------------------------------------------7---------8
510C I m p l i c i t T y p e s
511C-----------------------------------------------
512#include "implicit_f.inc"
513C-----------------------------------------------------------------
514C D u m m y A r g u m e n t s
515C-----------------------------------------------
516 INTEGER n ,nic,ic(*)
517C-----------------------------------------------
518C L o c a l V a r i a b l e s
519C-----------------------------------------------
520 INTEGER i,j
521C---- 6 ----- IC is already in an increasing order ----------------------------------- 7 -------- 8
522 intab0=0
523 if (n<ic(1).OR.n>ic(nic)) RETURN
524 IF (n<nic/2) THEN
525 DO i =1,nic
526 IF (n==ic(i)) THEN
527 intab0=i
528 RETURN
529 ENDIF
530 ENDDO
531 ELSE
532 DO i =nic,1,-1
533 IF (n==ic(i)) THEN
534 intab0=i
535 RETURN
536 ENDIF
537 ENDDO
538 ENDIF
539C
540 RETURN
541 END
542!||====================================================================
543!|| intab2 ../engine/source/implicit/imp_fsa_inv.F
544!||--- called by ------------------------------------------------------
545!|| get_suba ../engine/source/implicit/imp_pc_inv.F
546!|| sp_a2 ../engine/source/implicit/imp_fsa_inv.F
547!||--- calls -----------------------------------------------------
548!|| intab0 ../engine/source/implicit/imp_fsa_inv.F
549!||====================================================================
550 INTEGER FUNCTION intab2(NIC1,IC1,NIC2,IC2)
551C----6---------------------------------------------------------------7---------8
552C I m p l i c i t T y p e s
553C-----------------------------------------------
554#include "implicit_f.inc"
555C-----------------------------------------------------------------
556C D u m m y A r g u m e n t s
557C-----------------------------------------------
558 INTEGER nic1,ic1(*),nic2,ic2(*)
559C-----------------------------------------------
560C External function
561C-----------------------------------------------
562 INTEGER intab0
563 EXTERNAL intab0
564C-----------------------------------------------
565C L o c a l V a r i a b l e s
566C-----------------------------------------------
567 INTEGER i,j
568C----6-----ICi est deja en ordre croissante--return indice in IC1--7---------8
569 intab2=0
570 IF (ic1(nic1)<ic2(1).OR.ic2(nic2)<ic1(1)) RETURN
571C IF (NIC1>NIC2) THEN
572 DO i =1,nic2
573 intab2=intab0(nic1,ic1,ic2(i))
574 IF (intab2>0) RETURN
575 ENDDO
576 RETURN
577 END
578C ------factorized sparse approximate inverse version spmd-------
579!||====================================================================
580!|| imp_fsa_invp ../engine/source/implicit/imp_fsa_inv.F
581!||--- calls -----------------------------------------------------
582!|| ancmsg ../engine/source/output/message/message.F
583!|| arret ../engine/source/system/arret.F
584!|| get_subsp ../engine/source/implicit/imp_fsa_inv.F
585!|| imp_fsai ../engine/source/implicit/imp_fsa_inv.F
586!|| imp_pcg1 ../engine/source/implicit/imp_fsa_inv.F
587!|| sp_stat0 ../engine/source/implicit/imp_fsa_inv.F
588!||--- uses -----------------------------------------------------
589!|| message_mod ../engine/share/message_module/message_mod.F
590!||====================================================================
591 SUBROUTINE imp_fsa_invp(
592 1 NDDL ,NNZ ,IADK ,JDIK ,DIAG_K ,
593 2 LT_K ,DIAG_M ,LT_M ,MAXC ,MAX_A ,
594 3 NNE ,IDLFT0 ,IDLFT1,MAX_D )
595C-----------------------------------------------
596C M o d u l e s
597C-----------------------------------------------
598 USE message_mod
599C-----------------------------------------------
600C I m p l i c i t T y p e s
601C-----------------------------------------------
602#include "implicit_f.inc"
603C-----------------------------------------------
604C D u m m y A r g u m e n t s
605C-----------------------------------------------
606 INTEGER NDDL ,NNZ ,IADK(*),JDIK(*),MAXC ,MAX_A ,NNE,
607 . IDLFT0 ,IDLFT1,MAX_D
608C REAL
609 my_real
610 . DIAG_K(*), LT_K(*),DIAG_M(*), LT_M(*)
611C-----------------------------------------------
612C L o c a l V a r i a b l e s
613C-----------------------------------------------
614C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
615 INTEGER I,J,K,M,N,NC,IADA(MAXC+1),JDIA(MAX_A),JM(MAXC+1)
616 INTEGER MAX_L,IERR,I_CHK
617 my_real
618 . DIAG_A(MAXC),MJ(MAXC),
619 . diag_c(nddl-idlft1+1),lt_c(max_d+1)
620 my_real, DIMENSION(:),ALLOCATABLE :: lt_a
621C-----------------------------
622 ALLOCATE(lt_a(max_a),stat=ierr)
623 IF (ierr/=0) THEN
624 CALL ancmsg(msgid=19,anmode=aninfo,
625 . c1='FOR IMPLICIT PRECONDITION')
626 CALL arret(2)
627 ENDIF
628C--------------copy the useful part------
629 i_chk = nne
630 nne = 0
631 k=iadk(idlft1+1)-1
632 DO i=idlft1+1,nddl
633 diag_c(i-idlft1) = diag_m(i)
634 DO j=iadk(i),iadk(i+1)-1
635 lt_c(j-k)=lt_m(j)
636 ENDDO
637 ENDDO
638 DO i=idlft0+1,nddl
639 CALL sp_stat0(i ,iadk ,jdik ,nc ,jm )
640 CALL get_subsp(nddl ,iadk ,jdik ,diag_k ,lt_k ,
641 . nc ,iada ,jdia ,diag_a ,lt_a ,
642 . jm ,max_a ,idlft0,idlft1 ,diag_c,
643 . lt_c ,diag_m ,lt_m )
644 DO j=1,nc-1
645 mj(j)=zero
646 ENDDO
647 mj(nc)=one
648 IF (nc>10000) THEN
649 max_l=max_a
650 CALL imp_pcg1(
651 1 nc ,max_l ,iada ,jdia ,diag_a ,
652 2 lt_a ,mj ,ierr )
653 IF (i_chk>0.AND.ierr<0) nne = i
654 ELSE
655 max_l=1+(nc*(nc-1))/2
656 CALL imp_fsai(nc ,iada ,jdia ,diag_a ,lt_a ,
657 . max_l ,mj )
658 ENDIF
659C------------filtering----Diagonal is in LT_M (last one)--
660 diag_m(i)=mj(nc)
661 IF (diag_m(i)<em20) THEN
662 IF (nne==0.AND.i_chk==0) nne = i
663 diag_m(i)=abs(diag_m(i))
664 diag_m(i)=max(em20,diag_m(i))
665 ENDIF
666 DO k =1,nc-1
667 m=iadk(i)+k-1
668 lt_m(m)=mj(k)/diag_m(i)
669 ENDDO
670 IF (i_chk>0.AND.mj(nc)<em20) diag_m(i)=mj(nc)
671 ENDDO
672 DEALLOCATE(lt_a)
673C
674 RETURN
675 END
676!||====================================================================
677!|| imp_pcg1 ../engine/source/implicit/imp_fsa_inv.F
678!||--- called by ------------------------------------------------------
679!|| fsa_solv ../engine/source/implicit/imp_fsa_inv.F
680!|| imp_fsa_inv2 ../engine/source/implicit/imp_fsa_inv.F
681!|| imp_fsa_invh ../engine/source/implicit/imp_fsa_inv.F
682!|| imp_fsa_invh2 ../engine/source/implicit/imp_fsa_inv.F
683!|| imp_fsa_invp ../engine/source/implicit/imp_fsa_inv.F
684!|| imp_fsa_invp2 ../engine/source/implicit/imp_fsa_inv.F
685!||--- calls -----------------------------------------------------
686!|| crit_stop ../engine/source/implicit/imp_pcg.F
687!|| mav_lt ../engine/source/implicit/produt_v.F
688!|| produt_v0 ../engine/source/implicit/produt_v.F
689!||====================================================================
690 SUBROUTINE imp_pcg1(
691 1 NDDL ,NNZ ,IADK ,JDIK ,DIAG_K ,
692 2 LT_K ,R ,ISP )
693C-----------------------------------------------
694C I m p l i c i t T y p e s
695C-----------------------------------------------
696#include "implicit_f.inc"
697C-----------------------------------------------
698C C o m m o n B l o c k s
699C-----------------------------------------------
700#include "impl2_c.inc"
701C-----------------------------------------------
702C D u m m y A r g u m e n t s
703C-----------------------------------------------
704C----------resol [K]{X}={F}---------
705 INTEGER NDDL ,NNZ ,IADK(*) ,JDIK(*)
706C REAL
707 my_real
708 . DIAG_K(*), LT_K(*) ,R(*)
709C-----------------------------------------------
710C L o c a l V a r i a b l e s
711C-----------------------------------------------
712 INTEGER I,J,IT,IP,NLIM,ND,IPRE,NNZM,ISTOP,ITOL,ISP
713 my_real
714 . S , R2, R02,ALPHA,BETA,G0,G1,RR,TOLS,TOLN,TOLS2
715 my_real
716 . X(NDDL) ,P(NDDL) ,Z(NDDL) ,Y(NDDL),DIAG_M(NDDL)
717 INTEGER CRIT_STOP
718 EXTERNAL CRIT_STOP
719 my_real
720 . anorm2,xnorm2,l_a,l_b2,l_b,a_old,b_old,tmp,eps_m
721 my_real
722 . cs,dbar, delta, denom, kcond,snprod,qrnorm,
723 . gamma, gbar, gmax, gmin, epsln,lqnorm,diag,cgnorm,
724 . oldb, rhs1, rhs2,sn, zbar, zl ,oldb2,tnorm2,eps(4)
725C--------------INITIALISATION--------------------------
726 itol = 1
727 ipre = 2
728 nnzm = 0
729 nlim=max(nddl,2)
730 tols=sqrt(p_mach)
731 eps_m = p_mach
732 it=0
733 nd = nddl
734C-------------IT=0--------
735C------X(I)=ZERO--------
736 DO i=1,nddl
737 x(i) = zero
738 diag_m(i)=one/max(em20,diag_k(i))
739 ENDDO
740 CALL mav_lt(
741 1 nddl ,nnz ,iadk ,jdik ,diag_k,
742 2 lt_k ,x ,z )
743 DO i=1,nddl
744 r(i) = r(i)-z(i)
745 ENDDO
746 DO i=1,nddl
747 z(i) = r(i) *diag_m(i)
748 ENDDO
749 DO i=1,nddl
750 p(i) = z(i)
751 ENDDO
752 CALL produt_v0(nddl,r,z,g0)
753 CALL mav_lt(
754 1 nddl ,nnz ,iadk ,jdik ,diag_k,
755 2 lt_k ,p ,y )
756 CALL produt_v0(nddl,p,y,s)
757 alpha = g0/s
758 tols2=tols*tols
759 IF (itol==1) THEN
760 CALL produt_v0(nddl,r,r,r02)
761 r2 =r02
762 ELSEIF (itol==3) THEN
763C------ R2<TOL*TOL*ANORM2*XNORM2------
764 r02=abs(g0)
765 r2 =r02
766 l_a=one/alpha
767C L_B2=R02
768 tnorm2=l_a*l_a
769 anorm2=zero
770 xnorm2=zero
771 a_old=l_a
772 b_old=zero
773 oldb = sqrt(r02)
774 ELSEIF (itol==4) THEN
775 r02=alpha*alpha*abs(g0)
776 eps(1)=r02
777 r2=r02
778 ELSE
779 r02=abs(g0)
780 r2 =r02
781 ENDIF
782 IF (r02==zero) GOTO 200
783 toln=r02*tols2
784C-------to be consistent with lanzos for linear
785 it=1
786 DO i=1,nddl
787 x(i) = x(i) + alpha*p(i)
788 r(i) = r(i) - alpha*y(i)
789 ENDDO
790 DO i=1,nddl
791 z(i) = r(i) *diag_m(i)
792 ENDDO
793 CALL produt_v0(nddl,r,z,g1)
794 beta=g1/g0
795 IF (itol==1) THEN
796 CALL produt_v0(nddl,r,r,r2)
797 ELSEIF (itol==3) THEN
798C------ R2<TOL*TOL*ANORM2*XNORM2------
799 r2=abs(g1)
800 l_b2=abs(beta)*a_old*a_old
801 l_b=sqrt(l_b2)
802 tnorm2=tnorm2+l_b2
803 b_old=beta
804C* INITIALIZE OTHER QUANTITIES.
805 gbar = l_a
806 dbar = l_b
807 rhs1 = oldb
808 rhs2 = zero
809 gmax = abs( l_a ) + eps_m
810 gmin = gmax
811 oldb2 = l_b2
812 oldb = l_b
813 ELSEIF (itol==4) THEN
814 r2=r02
815 ELSE
816 r2=abs(g1)
817 ENDIF
818 g0 = g1
819 IF (itol==3) toln=toln*anorm2
820 istop=crit_stop(it,r2,nlim,toln)
821 DO WHILE (istop==1)
822 DO i=1,nddl
823 p(i) = z(i) + beta*p(i)
824 ENDDO
825 CALL mav_lt(
826 1 nddl ,nnz ,iadk ,jdik ,diag_k,
827 2 lt_k ,p ,y )
828 CALL produt_v0(nddl,p,y,s)
829 alpha=g0/s
830 DO i=1,nddl
831 x(i) = x(i) + alpha*p(i)
832 r(i) = r(i) - alpha*y(i)
833 ENDDO
834 DO i=1,nddl
835 z(i) = r(i) *diag_m(i)
836 ENDDO
837 CALL produt_v0(nddl,r,z,g1)
838 beta=g1/g0
839 g0 = g1
840 IF (itol==1) THEN
841 CALL produt_v0(nddl,r,r,r2)
842 ELSEIF (itol==3) THEN
843 r2 =abs(g1)
844 s=one/alpha
845 l_a=s+b_old*a_old
846C----------L_B2 : beta(j-1)^2-A_OLD :1/alpha(j-1)-------
847 l_b2=abs(beta)*s*s
848 l_b=sqrt(l_b2)
849 a_old=s
850 b_old=beta
851 anorm2=tnorm2
852 tnorm2=tnorm2+l_a*l_a+oldb2+l_b2
853 gamma = sqrt( gbar*gbar + oldb2 )
854 cs = gbar / gamma
855 sn = oldb / gamma
856 delta = cs * dbar + sn * l_a
857 gbar = sn * dbar - cs * l_a
858 epsln = sn * l_b
859 dbar = - cs * l_b
860 zl = rhs1 / gamma
861 xnorm2 = xnorm2+zl*zl
862 gmax = max( gmax, gamma )
863 gmin = min( gmin, gamma )
864 rhs1 = rhs2 - delta * zl
865 rhs2 = - epsln * zl
866 toln=tols2*anorm2*xnorm2
867 oldb2 = l_b2
868 oldb = l_b
869 ELSEIF (itol==4) THEN
870 tmp=alpha*alpha*abs(g1)
871 IF (it>=nd) THEN
872 DO j=1,nd-1
873 eps(j)=eps(j+1)
874 ENDDO
875 eps(nd)=tmp
876 r2=zero
877 DO j=1,nd
878 r2=r2+eps(j)
879 ENDDO
880 ELSE
881 eps(it+1)=tmp
882 r2=r2+tmp
883 ENDIF
884 ELSE
885 r2 =abs(g1)
886 ENDIF
887 istop=crit_stop(it,r2,nlim,toln)
888 it = it +1
889 ENDDO
890 200 CONTINUE
891 IF(it>=nlim)THEN
892 isp =-1
893 ELSE
894 isp = 0
895 ENDIF
896C RR = SQRT(R2/R02)
897C WRITE(*,1002)IT,RR
898C--------X->R--------
899 DO i=1,nddl
900 r(i) = x(i)
901 ENDDO
902C--------------------------------------------
903 RETURN
904 END
905!||====================================================================
906!|| imp_fsa_inv2 ../engine/source/implicit/imp_fsa_inv.F
907!||--- calls -----------------------------------------------------
908!|| ancmsg ../engine/source/output/message/message.F
909!|| arret ../engine/source/system/arret.F
910!|| get_subs0 ../engine/source/implicit/imp_fsa_inv.F
911!|| imp_fsai ../engine/source/implicit/imp_fsa_inv.F
912!|| imp_kfiltr ../engine/source/implicit/imp_fsa_inv.F
913!|| imp_pcg1 ../engine/source/implicit/imp_fsa_inv.F
914!|| sp_stat0 ../engine/source/implicit/imp_fsa_inv.F
915!||--- uses -----------------------------------------------------
916!|| message_mod ../engine/share/message_module/message_mod.F
917!||====================================================================
918 SUBROUTINE imp_fsa_inv2(
919 1 NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
920 2 IADM ,JDIM ,DIAG_M ,LT_M ,MAXC ,
921 3 MAX_A ,NNE ,D_TOL ,P_MACH)
922C-----------------------------------------------
923C M o d u l e s
924C-----------------------------------------------
925 USE message_mod
926C-----------------------------------------------
927C I m p l i c i t T y p e s
928C-----------------------------------------------
929#include "implicit_f.inc"
930C-----------------------------------------------
931C D u m m y A r g u m e n t s
932C-----------------------------------------------
933 INTEGER NDDL ,IADK(*),JDIK(*),MAXC ,MAX_A ,NNE,
934 . IADM(*),JDIM(*)
935C REAL
936 my_real
937 . diag_k(*), lt_k(*),diag_m(*), lt_m(*) ,d_tol ,p_mach
938C-----------------------------------------------
939C L o c a l V a r i a b l e s
940C-----------------------------------------------
941C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
942 INTEGER I,J,K,M,N,NC,IERR
943 INTEGER MAX_L,I_CHK
944 INTEGER, DIMENSION(:),ALLOCATABLE :: IADA,JDIA,JM
945 my_real, DIMENSION(:),ALLOCATABLE :: DIAG_A,LT_A,MJ
946C-----------------------------
947 I_CHK = nne
948 nne = 0
949 ALLOCATE(iada(maxc+1))
950 ALLOCATE(jdia(max_a))
951 ALLOCATE(jm(maxc+1))
952 ALLOCATE(diag_a(maxc))
953 ALLOCATE(mj(maxc))
954 ALLOCATE(lt_a(max_a),stat=ierr)
955 IF (ierr/=0) THEN
956 CALL ancmsg(msgid=19,anmode=aninfo,
957 . c1='FOR IMPLICIT PRECONDITION')
958 CALL arret(2)
959 ENDIF
960 DO i=1,nddl
961 CALL sp_stat0(i ,iadm ,jdim ,nc ,jm )
962 CALL get_subs0(nddl ,iadk ,jdik ,diag_k ,lt_k ,
963 . nc ,iada ,jdia ,diag_a ,lt_a ,
964 . jm ,max_a )
965 DO j=1,nc-1
966 mj(j)=zero
967 ENDDO
968 mj(nc)=one
969 IF (nc>10000) THEN
970 max_l=max_a
971 CALL imp_pcg1(
972 1 nc ,max_l ,iada ,jdia ,diag_a ,
973 2 lt_a ,mj ,ierr )
974 IF (i_chk>0.AND.ierr<0) nne = i
975 ELSE
976 max_l=1+(nc*(nc-1))/2
977 CALL imp_fsai(nc ,iada ,jdia ,diag_a ,lt_a ,
978 . max_l ,mj )
979 ENDIF
980C-----------------
981 diag_m(i)=mj(nc)
982 IF (diag_m(i)<em20) THEN
983 IF (nne==0.AND.i_chk==0) nne = i
984 diag_m(i)=abs(diag_m(i))
985 diag_m(i)=max(em20,diag_m(i))
986 ENDIF
987 DO k =1,nc-1
988 m=iadm(i)+k-1
989 lt_m(m)=mj(k)/diag_m(i)
990 ENDDO
991 IF (i_chk>0.AND.mj(nc)<em20) diag_m(i)=mj(nc)
992 ENDDO
993C
994 DEALLOCATE(iada,jdia)
995 DEALLOCATE(jm)
996 DEALLOCATE(diag_a,lt_a)
997 DEALLOCATE(mj)
998 k = 1
999 IF (d_tol>zero)
1000 . CALL imp_kfiltr(k ,nddl ,iadm ,jdim ,diag_m ,
1001 . lt_m ,d_tol ,p_mach,diag_k)
1002C
1003 RETURN
1004 END
1005!||====================================================================
1006!|| imp_fsa_invp2 ../engine/source/implicit/imp_fsa_inv.F
1007!||--- calls -----------------------------------------------------
1008!|| ancmsg ../engine/source/output/message/message.F
1009!|| arret ../engine/source/system/arret.F
1010!|| get_subsp ../engine/source/implicit/imp_fsa_inv.F
1011!|| imp_fsai ../engine/source/implicit/imp_fsa_inv.F
1012!|| imp_kfiltr ../engine/source/implicit/imp_fsa_inv.F
1013!|| imp_pcg1 ../engine/source/implicit/imp_fsa_inv.F
1014!|| sp_stat0 ../engine/source/implicit/imp_fsa_inv.F
1015!||--- uses -----------------------------------------------------
1016!|| message_mod ../engine/share/message_module/message_mod.F
1017!||====================================================================
1018 SUBROUTINE imp_fsa_invp2(
1019 1 NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
1020 2 IADM ,JDIM ,DIAG_M ,LT_M ,MAXC ,
1021 3 MAX_A ,NNE ,IDLFT0 ,IDLFT1 ,MAX_D ,
1022 4 D_TOL ,P_MACH)
1023C-----------------------------------------------
1024C M o d u l e s
1025C-----------------------------------------------
1026 USE message_mod
1027C-----------------------------------------------
1028C I m p l i c i t T y p e s
1029C-----------------------------------------------
1030#include "implicit_f.inc"
1031C-----------------------------------------------
1032C D u m m y A r g u m e n t s
1033C-----------------------------------------------
1034 INTEGER NDDL ,NNZ ,IADK(*),JDIK(*),MAXC ,MAX_A ,NNE,
1035 . IDLFT0 ,IDLFT1,MAX_D,IADM(*),JDIM(*)
1036C REAL
1037 my_real
1038 . DIAG_K(*), LT_K(*),DIAG_M(*), LT_M(*),D_TOL ,P_MACH
1039C-----------------------------------------------
1040C L o c a l V a r i a b l e s
1041C-----------------------------------------------
1042C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
1043 INTEGER I,J,K,M,N,NC,IADA(MAXC+1),JDIA(MAX_A),JM(MAXC+1)
1044 INTEGER MAX_L,IERR,I_CHK
1045 my_real
1046 . DIAG_A(MAXC),MJ(MAXC),
1047 . DIAG_C(NDDL-IDLFT1+1),LT_C(MAX_D+1)
1048 my_real, DIMENSION(:),ALLOCATABLE :: lt_a
1049C-----------------------------
1050 ALLOCATE(lt_a(max_a),stat=ierr)
1051 IF (ierr/=0) THEN
1052 CALL ancmsg(msgid=19,anmode=aninfo,
1053 . c1='FOR IMPLICIT PRECONDITION')
1054 CALL arret(2)
1055 ENDIF
1056C--------------copy the useful part------
1057 i_chk = nne
1058 nne = 0
1059 k=iadk(idlft1+1)-1
1060 DO i=idlft1+1,nddl
1061 diag_c(i-idlft1) = diag_m(i)
1062 DO j=iadk(i),iadk(i+1)-1
1063 lt_c(j-k)=lt_m(j)
1064 ENDDO
1065 ENDDO
1066 DO i=idlft0+1,nddl
1067 CALL sp_stat0(i ,iadm ,jdim ,nc ,jm )
1068 CALL get_subsp(nddl ,iadk ,jdik ,diag_k ,lt_k ,
1069 . nc ,iada ,jdia ,diag_a ,lt_a ,
1070 . jm ,max_a ,idlft0,idlft1 ,diag_c,
1071 . lt_c ,diag_m,lt_m )
1072 DO j=1,nc-1
1073 mj(j)=zero
1074 ENDDO
1075 mj(nc)=one
1076 IF (nc>10000) THEN
1077 max_l=max_a
1078 CALL imp_pcg1(
1079 1 nc ,max_l ,iada ,jdia ,diag_a ,
1080 2 lt_a ,mj ,ierr )
1081 IF (i_chk>0.AND.ierr<0) nne = i
1082 ELSE
1083 max_l=1+(nc*(nc-1))/2
1084 CALL imp_fsai(nc ,iada ,jdia ,diag_a ,lt_a ,
1085 . max_l ,mj )
1086 ENDIF
1087C------------Diagonal is in LT_M (last one)--
1088 diag_m(i)=mj(nc)
1089 IF (diag_m(i)<em20) THEN
1090 IF (nne==0.AND.i_chk==0) nne = i
1091 diag_m(i)=abs(diag_m(i))
1092 diag_m(i)=max(em20,diag_m(i))
1093 ENDIF
1094 DO k =1,nc-1
1095 m=iadm(i)+k-1
1096 lt_m(m)=mj(k)/diag_m(i)
1097 ENDDO
1098 IF (i_chk>0.AND.mj(nc)<em20) diag_m(i)=mj(nc)
1099 ENDDO
1100C
1101 DEALLOCATE(lt_a)
1102 k = idlft0+1
1103 IF (d_tol>zero)
1104 . CALL imp_kfiltr(k ,nddl ,iadm ,jdim ,diag_m ,
1105 . lt_m ,d_tol ,p_mach,diag_k)
1106C
1107 RETURN
1108 END
1109!||====================================================================
1110!|| imp_kfiltr ../engine/source/implicit/imp_fsa_inv.F
1111!||--- called by ------------------------------------------------------
1112!|| imp_fsa_inv2 ../engine/source/implicit/imp_fsa_inv.F
1113!|| imp_fsa_inv2hp ../engine/source/implicit/imp_fsa_inv.F
1114!|| imp_fsa_invh2 ../engine/source/implicit/imp_fsa_inv.F
1115!|| imp_fsa_invp2 ../engine/source/implicit/imp_fsa_inv.F
1116!||--- calls -----------------------------------------------------
1117!|| ancmsg ../engine/source/output/message/message.F
1118!|| arret ../engine/source/system/arret.F
1119!|| cp_int ../engine/source/implicit/produt_v.F
1120!|| cp_real ../engine/source/implicit/produt_v.F
1121!||--- uses -----------------------------------------------------
1122!|| message_mod ../engine/share/message_module/message_mod.F
1123!||====================================================================
1124 SUBROUTINE imp_kfiltr(NDF ,ND ,IADA ,JDIA ,DIAG_A ,
1125 . LT_A ,TOL ,E_PS ,DIAG_K )
1126C-----------------------------------------------
1127C M o d u l e s
1128C-----------------------------------------------
1129 USE message_mod
1130C-----------------------------------------------
1131C I m p l i c i t T y p e s
1132C-----------------------------------------------
1133#include "implicit_f.inc"
1134C-----------------------------------------------
1135C D u m m y A r g u m e n t s
1136C-----------------------------------------------
1137 INTEGER NDF,ND ,IADA(*) ,JDIA(*)
1138C REAL
1139 my_real
1140 . DIAG_A(*),LT_A(*),TOL,E_PS,DIAG_K(*)
1141C-----------------------------------------------
1142C L o c a l V a r i a b l e s
1143C-----------------------------------------------
1144 INTEGER I,J,K,NZ,IERR,MNZ,INORM
1145 INTEGER, DIMENSION(:),ALLOCATABLE :: IADL,JDIL
1146 my_real
1147 . min_d,max_d,mtol,dd,taux
1148 my_real, DIMENSION(:),ALLOCATABLE :: lt_l
1149C-----------------------------
1150 print *,'D_tol,p_mach=',tol,e_ps
1151 nz = iada(nd+1)-iada(1)
1152 DO i = 1, ndf-1
1153 diag_a(i) = zero
1154 DO j = iada(i), iada(i+1)-1
1155 lt_a(j) = zero
1156 ENDDO
1157 iada(i+1) = 1
1158 ENDDO
1159 IF (nz>0.AND.tol>zero) THEN
1160 ALLOCATE(iadl(nd+1))
1161 ALLOCATE(jdil(nz),lt_l(nz),stat=ierr)
1162 IF (ierr/=0) THEN
1163 CALL ancmsg(msgid=19,anmode=aninfo,
1164 . c1='FOR IMPLICIT PRECONDITION')
1165 CALL arret(2)
1166 ENDIF
1167C-----------------------------------------------
1168 CALL cp_int(nd+1,iada,iadl)
1169 CALL cp_int(nz,jdia,jdil)
1170 CALL cp_real(nz,lt_a,lt_l)
1171 max_d = zero
1172 min_d = ep20
1173 DO i = ndf, nd
1174 max_d = max(max_d,diag_a(i))
1175 min_d = min(min_d,diag_a(i))
1176 ENDDO
1177 print *,'max_d,min_d=',max_d,min_d
1178c MTOL = TOL*MIN_D
1179c MTOL = MAX(E_PS,MTOL)
1180C------post-filtration------
1181 mnz = 0
1182 DO i = ndf, nd
1183 DO j = iadl(i), iadl(i+1)-1
1184 mtol = tol*min(diag_a(jdil(j)),diag_a(i))
1185 mtol = max(e_ps,mtol)
1186 IF (abs(lt_l(j))>mtol) THEN
1187 mnz = mnz + 1
1188 jdia(mnz) = jdil(j)
1189 lt_a(mnz) = lt_l(j)
1190 ENDIF
1191 ENDDO
1192 iada(i+1) = mnz + 1
1193 ENDDO
1194 DEALLOCATE(iadl,jdil)
1195 DEALLOCATE(lt_l)
1196 taux = one*mnz/nz
1197c print *,'mnz,nz=',mnz,nz,mtol
1198 print *,'filtrage factor=',taux
1199C------retrun from M normalized------
1200 ENDIF
1201 inorm = 0
1202 IF (inorm>zero) THEN
1203 DO i = ndf, nd
1204 diag_a(i) = diag_a(i)/diag_k(i)
1205 DO j = iada(i), iada(i+1)-1
1206 dd = sqrt(diag_k(i)/diag_k(jdia(j)))
1207 lt_a(j) = dd*lt_a(j)
1208 ENDDO
1209 ENDDO
1210 ENDIF
1211C--------------------------------------------
1212 RETURN
1213 END
1214!||====================================================================
1215!|| get_subsn ../engine/source/implicit/imp_fsa_inv.f
1216!||--- calls -----------------------------------------------------
1217!|| ind_lt2ln ../engine/source/implicit/imp_fsa_inv.F
1218!|| intab0 ../engine/source/implicit/imp_fsa_inv.F
1219!||====================================================================
1220 SUBROUTINE get_subsn(NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
1221 . NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
1222 . JM ,MAXA )
1223C-----------------------------------------------
1224C I m p l i c i t T y p e s
1225C-----------------------------------------------
1226#include "implicit_f.inc"
1227C-----------------------------------------------
1228C D u m m y A r g u m e n t s
1229C-----------------------------------------------
1230 INTEGER NDDL ,IADK(*) ,JDIK(*),IADA(*) ,JDIA(*)
1231 INTEGER NC ,JM(*),MAXA
1232C REAL
1233 my_real
1234 . lt_k(*),diag_k(*),lt_a(*),diag_a(*)
1235C-----------------------------------------------
1236C External function
1237C-----------------------------------------------
1238 INTEGER INTAB0
1239 EXTERNAL INTAB0
1240C-----------------------------------------------
1241C L o c a l V a r i a b l e s
1242C-----------------------------------------------
1243 INTEGER I,J,K,JJ,NNZA,N
1244 my_real
1245 . dd
1246C--------------------------------------------
1247 nnza=0
1248 iada(1)=1
1249 DO i=1,nc
1250 j=jm(i)
1251 diag_a(i)=one
1252 DO k=iadk(j),iadk(j+1)-1
1253 jj=jdik(k)
1254 n=intab0(nc,jm,jj)
1255 IF (n>0) THEN
1256 dd = sqrt(diag_k(j)*diag_k(jj))
1257 nnza=nnza+1
1258 jdia(nnza)=n
1259 lt_a(nnza)=lt_k(k)/dd
1260 ENDIF
1261 ENDDO
1262 iada(i+1)=nnza+1
1263 ENDDO
1264 CALL ind_lt2ln(nc,iada ,jdia ,lt_a, nnza )
1265C
1266 RETURN
1267 END
1268!||====================================================================
1269!|| ind_lt2ln ../engine/source/implicit/imp_fsa_inv.F
1270!||--- called by ------------------------------------------------------
1271!|| get_subs0 ../engine/source/implicit/imp_fsa_inv.F
1272!|| get_subsn ../engine/source/implicit/imp_fsa_inv.F
1273!|| get_subsp ../engine/source/implicit/imp_fsa_inv.F
1274!|| get_subsp_sms ../engine/source/ams/sms_fsa_inv.F
1275!||--- calls -----------------------------------------------------
1276!|| cp_int ../engine/source/implicit/produt_v.F
1277!|| cp_real ../engine/source/implicit/produt_v.F
1278!||====================================================================
1279 SUBROUTINE ind_lt2ln(NDDL,IADK ,JDIK ,LT_K, MAXL )
1280C-------------
1281C-----------------------------------------------
1282C I m p l i c i t T y p e s
1283C-----------------------------------------------
1284#include "implicit_f.inc"
1285C-----------------------------------------------
1286C D u m m y A r g u m e n t s
1287C-----------------------------------------------
1288 INTEGER NDDL, IADK(*),JDIK(*),MAXL
1289 my_real
1290 . lt_k(*)
1291C REAL
1292C-----------------------------------------------
1293C L o c a l V a r i a b l e s
1294C-----------------------------------------------
1295 INTEGER IADM(NDDL+1),JDIM(MAXL),ICOL(NDDL)
1296 INTEGER I,JD,J,K,N,NM
1297 my_real
1298 . LT_M(MAXL)
1299C-----------------------------------------------
1300 CALL CP_INT(NDDL+1,IADK,IADM)
1301 CALL CP_INT(MAXL,JDIK,JDIM)
1302 CALL CP_REAL(MAXL,LT_K,LT_M)
1303C
1304 DO i = 1, nddl
1305 icol(i) = 0
1306 DO j = iadm(i),iadm(i+1)-1
1307 jd = jdim(j)
1308 icol(jd) = icol(jd) + 1
1309 ENDDO
1310 ENDDO
1311C
1312 iadk(1) = 1
1313 DO i = 1,nddl
1314 iadk(i+1) = iadk(i)+icol(i)
1315 ENDDO
1316C
1317 DO i = 1,nddl
1318 icol(i) = 0
1319 DO j=iadm(i),iadm(i+1)-1
1320 jd = jdim(j)
1321 k = iadk(jd) + icol(jd)
1322 jdik(k) = i
1323 lt_k(k) = lt_m(j)
1324 icol(jd) = icol(jd) + 1
1325 ENDDO
1326 ENDDO
1327C--------------------------------------------
1328 RETURN
1329 END
1330C ------factorized sparse approximate inverse version hybrid-------
1331!||====================================================================
1332!|| imp_fsa_invh ../engine/source/implicit/imp_fsa_inv.F
1333!||--- calls -----------------------------------------------------
1334!|| ancmsg ../engine/source/output/message/message.F
1335!|| arret ../engine/source/system/arret.F
1336!|| get_subsp ../engine/source/implicit/imp_fsa_inv.F
1337!|| imp_fsai ../engine/source/implicit/imp_fsa_inv.F
1338!|| imp_pcg1 ../engine/source/implicit/imp_fsa_inv.F
1339!|| my_barrier ../engine/source/system/machine.F
1340!|| sp_stat0 ../engine/source/implicit/imp_fsa_inv.F
1341!||--- uses -----------------------------------------------------
1342!|| message_mod ../engine/share/message_module/message_mod.F
1343!||====================================================================
1344 SUBROUTINE imp_fsa_invh(
1345 1 NDDL ,NNZ ,IADK ,JDIK ,DIAG_K ,
1346 2 LT_K ,DIAG_M ,LT_M ,MAXC ,MAX_A ,
1347 3 NNE ,IDLFT0 ,IDLFT1,MAX_D ,ITASK )
1348C-----------------------------------------------
1349C M o d u l e s
1350C-----------------------------------------------
1351 USE message_mod
1352C-----------------------------------------------
1353C I m p l i c i t T y p e s
1354C-----------------------------------------------
1355#include "implicit_f.inc"
1356C-----------------------------------------------
1357C D u m m y A r g u m e n t s
1358C-----------------------------------------------
1359 INTEGER NDDL ,NNZ ,IADK(*),JDIK(*),MAXC ,MAX_A ,NNE,
1360 . IDLFT0 ,IDLFT1,MAX_D,ITASK
1361C REAL
1362 my_real
1363 . DIAG_K(*), LT_K(*),DIAG_M(*), LT_M(*)
1364C-----------------------------------------------
1365C L o c a l V a r i a b l e s
1366C-----------------------------------------------
1367C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
1368 INTEGER I,J,K,M,N,NC,MAX_L,IERR,I_CHK,IER1,
1369 . JM(MAXC+1)
1370 my_real
1371 . DIAG_C(NDDL-IDLFT1+1),LT_C(MAX_D+1)
1372 INTEGER, DIMENSION(:),ALLOCATABLE :: IADA,JDIA
1373 my_real, DIMENSION(:),ALLOCATABLE :: DIAG_A,LT_A,MJ
1374C-----------------------------
1375 I_CHK = nne
1376 nne = 0
1377 IF ((idlft0+1)>nddl) RETURN
1378C
1379 ALLOCATE(iada(maxc+1),diag_a(maxc),mj(maxc),stat=ier1)
1380 ALLOCATE(lt_a(max_a),jdia(max_a),stat=ierr)
1381
1382 IF ((ierr+ier1)/=0) THEN
1383 IF (itask == 0 ) THEN
1384 CALL ancmsg(msgid=19,anmode=aninfo,
1385 . c1='FOR IMPLICIT PRECONDITION')
1386 CALL arret(2)
1387 END IF !(ITASK == 0 ) THEN
1388 ENDIF
1389
1390C--------------copy the useful part------
1391 k=iadk(idlft1+1)-1
1392 DO i=idlft1+1,nddl
1393 diag_c(i-idlft1) = diag_m(i)
1394 DO j=iadk(i),iadk(i+1)-1
1395 lt_c(j-k)=lt_m(j)
1396 ENDDO
1397 ENDDO
1398C----------------------
1399 CALL my_barrier
1400C---------------------
1401C
1402C SMP dynamic parallel loop
1403C
1404!$OMP DO SCHEDULE(DYNAMIC,1)
1405 DO i=idlft0+1,nddl
1406 CALL sp_stat0(i ,iadk ,jdik ,nc ,jm )
1407 CALL get_subsp(nddl ,iadk ,jdik ,diag_k ,lt_k ,
1408 . nc ,iada ,jdia ,diag_a ,lt_a ,
1409 . jm ,max_a ,idlft0,idlft1 ,diag_c,
1410 . lt_c ,diag_m ,lt_m )
1411 DO j=1,nc-1
1412 mj(j)=zero
1413 ENDDO
1414 mj(nc)=one
1415C
1416 IF (nc>10000) THEN
1417 max_l=max_a
1418 CALL imp_pcg1(
1419 1 nc ,max_l ,iada ,jdia ,diag_a ,
1420 2 lt_a ,mj ,ierr )
1421C
1422 IF (i_chk>0.AND.ierr<0) nne = i
1423 ELSE
1424C
1425 max_l=1+(nc*(nc-1))/2
1426 CALL imp_fsai(nc ,iada ,jdia ,diag_a ,lt_a ,
1427 . max_l ,mj )
1428C
1429 ENDIF
1430C------------filtering----Diagonal is in LT_M (last one)--
1431 diag_m(i)=mj(nc)
1432 IF (diag_m(i)<em20) THEN
1433 IF (nne==0.AND.i_chk==0) nne = i
1434 diag_m(i)=abs(diag_m(i))
1435 diag_m(i)=max(em20,diag_m(i))
1436 ENDIF
1437 DO k =1,nc-1
1438 m=iadk(i)+k-1
1439 lt_m(m)=mj(k)/diag_m(i)
1440 ENDDO
1441C
1442 IF (i_chk>0.AND.mj(nc)<em20) diag_m(i)=mj(nc)
1443 ENDDO
1444
1445!$OMP END DO
1446C
1447 DEALLOCATE(iada,diag_a,mj)
1448 DEALLOCATE(lt_a,jdia)
1449C
1450 RETURN
1451 END
1452!||====================================================================
1453!|| imp_fsa_invh2 ../engine/source/implicit/imp_fsa_inv.F
1454!||--- calls -----------------------------------------------------
1455!|| ancmsg ../engine/source/output/message/message.F
1456!|| arret ../engine/source/system/arret.F
1457!|| get_subsp ../engine/source/implicit/imp_fsa_inv.F
1458!|| imp_fsai ../engine/source/implicit/imp_fsa_inv.F
1459!|| imp_kfiltr ../engine/source/implicit/imp_fsa_inv.F
1460!|| imp_pcg1 ../engine/source/implicit/imp_fsa_inv.F
1461!|| my_barrier ../engine/source/system/machine.F
1462!|| sp_stat0 ../engine/source/implicit/imp_fsa_inv.F
1463!||--- uses -----------------------------------------------------
1464!|| message_mod ../engine/share/message_module/message_mod.F
1465!||====================================================================
1466 SUBROUTINE imp_fsa_invh2(
1467 1 NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
1468 2 IADM ,JDIM ,DIAG_M ,LT_M ,MAXC ,
1469 3 MAX_A ,NNE ,IDLFT0 ,IDLFT1 ,MAX_D ,
1470 4 D_TOL ,P_MACH ,ITASK )
1471C-----------------------------------------------
1472C M o d u l e s
1473C-----------------------------------------------
1474 USE message_mod
1475C-----------------------------------------------
1476C I m p l i c i t T y p e s
1477C-----------------------------------------------
1478#include "implicit_f.inc"
1479C-----------------------------------------------
1480C D u m m y A r g u m e n t s
1481C-----------------------------------------------
1482 INTEGER NDDL ,NNZ ,IADK(*),JDIK(*),MAXC ,MAX_A ,NNE,
1483 . idlft0 ,idlft1,max_d,iadm(*),jdim(*),itask
1484C REAL
1485 my_real
1486 . diag_k(*), lt_k(*),diag_m(*), lt_m(*),d_tol ,p_mach
1487C-----------------------------------------------
1488C L o c a l V a r i a b l e s
1489C-----------------------------------------------
1490C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
1491 INTEGER I,J,K,M,N,NC,MAX_L,IERR,I_CHK,IER1,
1492 . JM(MAXC+1)
1493 my_real
1494 . DIAG_C(NDDL-IDLFT1+1),LT_C(MAX_D+1)
1495 INTEGER, DIMENSION(:),ALLOCATABLE :: IADA,JDIA
1496 my_real, DIMENSION(:),ALLOCATABLE :: DIAG_A,LT_A,MJ
1497C-----------------------------
1498 i_chk = nne
1499 nne = 0
1500 IF ((idlft0+1)>nddl) RETURN
1501C
1502 ALLOCATE(iada(maxc+1),diag_a(maxc),mj(maxc),stat=ier1)
1503 ALLOCATE(lt_a(max_a),jdia(max_a),stat=ierr)
1504C--------
1505 IF ((ierr+ier1)/=0) THEN
1506 IF (itask == 0 ) THEN
1507 CALL ancmsg(msgid=19,anmode=aninfo,
1508 . c1='FOR IMPLICIT PRECONDITION')
1509 CALL arret(2)
1510 END IF !(ITASK == 0 ) THEN
1511 ENDIF
1512C--------------copy the useful part------
1513 k=iadk(idlft1+1)-1
1514 DO i=idlft1+1,nddl
1515 diag_c(i-idlft1) = diag_m(i)
1516 DO j=iadk(i),iadk(i+1)-1
1517 lt_c(j-k)=lt_m(j)
1518 ENDDO
1519 ENDDO
1520C----------------------
1521 CALL my_barrier
1522C---------------------
1523C SMP dynamic parallel loop
1524C
1525!$OMP DO SCHEDULE(DYNAMIC,1)
1526 DO i=idlft0+1,nddl
1527 CALL sp_stat0(i ,iadm ,jdim ,nc ,jm )
1528 CALL get_subsp(nddl ,iadk ,jdik ,diag_k ,lt_k ,
1529 . nc ,iada ,jdia ,diag_a ,lt_a ,
1530 . jm ,max_a ,idlft0,idlft1 ,diag_c,
1531 . lt_c ,diag_m,lt_m )
1532 DO j=1,nc-1
1533 mj(j)=zero
1534 ENDDO
1535 mj(nc)=one
1536 IF (nc>10000) THEN
1537 max_l=max_a
1538 CALL imp_pcg1(
1539 1 nc ,max_l ,iada ,jdia ,diag_a ,
1540 2 lt_a ,mj ,ierr )
1541 IF (i_chk>0.AND.ierr<0) nne = i
1542 ELSE
1543 max_l=1+(nc*(nc-1))/2
1544 CALL imp_fsai(nc ,iada ,jdia ,diag_a ,lt_a ,
1545 . max_l ,mj )
1546 ENDIF
1547C------------Diagonal is in LT_M (last one)--
1548 diag_m(i)=mj(nc)
1549 IF (diag_m(i)<em20) THEN
1550 IF (nne==0.AND.i_chk==0) nne = i
1551 diag_m(i)=abs(diag_m(i))
1552 diag_m(i)=max(em20,diag_m(i))
1553 ENDIF
1554 DO k =1,nc-1
1555 m=iadm(i)+k-1
1556 lt_m(m)=mj(k)/diag_m(i)
1557 ENDDO
1558C
1559 IF (i_chk>0.AND.mj(nc)<em20) diag_m(i)=mj(nc)
1560 ENDDO
1561C
1562!$OMP END DO
1563C
1564 DEALLOCATE(iada,diag_a,mj)
1565 DEALLOCATE(lt_a,jdia)
1566C
1567 IF (itask == 0 ) THEN
1568 k = idlft0+1
1569 IF (d_tol>zero)
1570 . CALL imp_kfiltr(k ,nddl ,iadm ,jdim ,diag_m ,
1571 . lt_m ,d_tol ,p_mach,diag_k)
1572 END IF
1573C
1574 RETURN
1575 END
1576!||====================================================================
1577!|| sp_dim ../engine/source/implicit/imp_fsa_inv.F
1578!||--- called by ------------------------------------------------------
1579!|| imp_fsa_inv2hp ../engine/source/implicit/imp_fsa_inv.F
1580!|| imp_fsa_invhp ../engine/source/implicit/imp_fsa_inv.F
1581!||====================================================================
1582 SUBROUTINE sp_dim(IL ,IADK ,JDIK ,NC ,MAX_A ,MAX_L )
1583C-----------------------------------------------
1584C I m p l i c i t T y p e s
1585C-----------------------------------------------
1586#include "implicit_f.inc"
1587C-----------------------------------------------
1588C D u m m y A r g u m e n t s
1589C-----------------------------------------------
1590 INTEGER IL, IADK(*) ,JDIK(*),NC,MAX_A ,MAX_L
1591C REAL
1592C-----------------------------------------------
1593C L o c a l V a r i a b l e s
1594C-----------------------------------------------
1595 INTEGER I,J,K
1596C-----------------------------------------------
1597 nc=0
1598 DO k =iadk(il),iadk(il+1)-1
1599 nc=nc+1
1600 ENDDO
1601 nc=nc+1
1602C------MAX_L<- MAX_A
1603 IF (nc <= 10000) THEN
1604 max_l = 1+(nc*(nc-1))/2
1605 ELSE
1606 max_l = max_a
1607 END IF
1608C--------------------------------------------
1609 RETURN
1610 END
1611!||====================================================================
1612!|| fsa_solv ../engine/source/implicit/imp_fsa_inv.F
1613!||--- called by ------------------------------------------------------
1614!|| imp_fsa_inv2hp ../engine/source/implicit/imp_fsa_inv.f
1615!|| imp_fsa_invhp ../engine/source/implicit/imp_fsa_inv.F
1616!||--- calls -----------------------------------------------------
1617!|| get_subsp ../engine/source/implicit/imp_fsa_inv.F
1618!|| imp_fsai ../engine/source/implicit/imp_fsa_inv.F
1619!|| imp_pcg1 ../engine/source/implicit/imp_fsa_inv.f
1620!||====================================================================
1621 SUBROUTINE fsa_solv(
1622 1 NDDL ,NC ,IADK ,JDIK ,DIAG_K ,
1623 2 LT_K ,DIAG_M,LT_M ,DIAG_C,LT_C ,
1624 3 MAX_A ,IDLFT0,IDLFT1,NNE ,I_CHK ,
1625 4 IADM ,JDIM ,I )
1626C-----------------------------------------------
1627C I m p l i c i t T y p e s
1628C-----------------------------------------------
1629#include "implicit_f.inc"
1630C-----------------------------------------------
1631C D u m m y A r g u m e n t s
1632C-----------------------------------------------
1633 INTEGER I,NDDL ,NC ,IADK(*),JDIK(*),MAX_A ,NNE,
1634 . idlft0,idlft1 ,i_chk,iadm(*),jdim(*)
1635C REAL
1636 my_real
1637 . diag_k(*), lt_k(*),diag_m(*), lt_m(*) ,diag_c(*),lt_c(*)
1638C-----------------------------------------------
1639C L o c a l V a r i a b l e s
1640C-----------------------------------------------
1641C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
1642 INTEGER J,K,M,N,MAX_L,IERR,IER1,IADA(NC+1),JM(NC)
1643c . IADA(NC+1),JDIA(MAX_A),JM(NC)
1644c my_real
1645c . DIAG_A(NC),LT_A(MAX_A),MJ(NC)
1646 INTEGER, DIMENSION(:),ALLOCATABLE :: JDIA
1647 my_real, DIMENSION(:),ALLOCATABLE :: DIAG_A,LT_A,MJ
1648C-----------------------------
1649C
1650 ALLOCATE(DIAG_A(NC),MJ(NC),STAT=ier1)
1651 ALLOCATE(lt_a(max_a),jdia(max_a),stat=ierr)
1652C-----------------------------
1653 j=0
1654 DO k =iadm(i),iadm(i+1)-1
1655 j=j+1
1656 jm(j)=jdim(k)
1657 ENDDO
1658 j=j+1
1659 jm(j)=i
1660C
1661 CALL get_subsp(nddl ,iadk ,jdik ,diag_k ,lt_k ,
1662 . nc ,iada ,jdia ,diag_a ,lt_a ,
1663 . jm ,max_a ,idlft0,idlft1 ,diag_c,
1664 . lt_c ,diag_m ,lt_m )
1665 DO j=1,nc-1
1666 mj(j)=zero
1667 ENDDO
1668 mj(nc)=one
1669C
1670 IF (nc > 10000) THEN
1671 max_l=max_a
1672 CALL imp_pcg1(
1673 1 nc ,max_l ,iada ,jdia ,diag_a ,
1674 2 lt_a ,mj ,ierr )
1675C
1676 IF (i_chk>0.AND.ierr<0) nne = i
1677 ELSE
1678C
1679 max_l=1+(nc*(nc-1))/2
1680 CALL imp_fsai(nc ,iada ,jdia ,diag_a ,lt_a ,
1681 . max_l ,mj )
1682C
1683 ENDIF
1684C------------filtering----Diagonal is in LT_M (last one)--
1685 diag_m(i)=mj(nc)
1686 IF (diag_m(i)<em20) THEN
1687 IF (nne==0.AND.i_chk==0) nne = i
1688 diag_m(i)=abs(diag_m(i))
1689 diag_m(i)=max(em20,diag_m(i))
1690 ENDIF
1691 DO k =1,nc-1
1692 m=iadm(i)+k-1
1693 lt_m(m)=mj(k)/diag_m(i)
1694 ENDDO
1695 IF (i_chk>0.AND.mj(nc)<em20) diag_m(i)=mj(nc)
1696C
1697 DEALLOCATE(diag_a,mj)
1698 DEALLOCATE(lt_a,jdia)
1699 RETURN
1700 END
1701C ------factorized sparse approximate inverse version hybrid SMP inside
1702!||====================================================================
1703!|| imp_fsa_invhp ../engine/source/implicit/imp_fsa_inv.F
1704!||--- called by ------------------------------------------------------
1705!|| lin_solvh1 ../engine/source/implicit/lin_solv.F
1706!||--- calls -----------------------------------------------------
1707!|| fsa_solv ../engine/source/implicit/imp_fsa_inv.F
1708!|| my_barrier ../engine/source/system/machine.F
1709!|| omp_get_thread_num ../engine/source/engine/openmp_stub.F90
1710!|| sp_dim ../engine/source/implicit/imp_fsa_inv.f
1711!||--- uses -----------------------------------------------------
1712!|| message_mod ../engine/share/message_module/message_mod.F
1713!||====================================================================
1714 SUBROUTINE imp_fsa_invhp(
1715 1 NDDL ,NNZ ,IADK ,JDIK ,DIAG_K ,
1716 2 LT_K ,DIAG_M ,LT_M ,MAXC ,MAX_A ,
1717 3 NNE ,IDLFT0 ,IDLFT1,MAX_D )
1718C-----------------------------------------------
1719C M o d u l e s
1720C-----------------------------------------------
1721 USE message_mod
1722C-----------------------------------------------
1723C I m p l i c i t T y p e s
1724C-----------------------------------------------
1725#include "implicit_f.inc"
1726C-----------------------------------------------
1727C C o m m o n B l o c k s
1728C-----------------------------------------------
1729#include "task_c.inc"
1730C-----------------------------------------------
1731C D u m m y A r g u m e n t s
1732C-----------------------------------------------
1733 INTEGER NDDL ,NNZ ,IADK(*),JDIK(*),MAXC ,MAX_A ,NNE,
1734 . idlft0 ,idlft1,max_d
1735C REAL
1736 my_real
1737 . diag_k(*), lt_k(*),diag_m(*), lt_m(*)
1738C-----------------------------------------------
1739C L o c a l V a r i a b l e s
1740C-----------------------------------------------
1741C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
1742 INTEGER I,J,K,M,N,NC,MAX_L,IERR,I_CHK,IER1,
1743 . ITSK,F_DDL,L_DDL,N1
1744 my_real
1745 . DIAG_C(NDDL-IDLFT1+1),LT_C(MAX_D+1)
1746 INTEGER OMP_GET_THREAD_NUM
1747 EXTERNAL OMP_GET_THREAD_NUM
1748C-----------------------------
1749 I_CHK = nne
1750 nne = 0
1751 IF ((idlft0+1)>nddl) RETURN
1752C
1753C--------------copy the useful part------
1754 k=iadk(idlft1+1)-1
1755!$OMP PARALLEL PRIVATE(ITSK,F_DDL,L_DDL,NC,MAX_L,I,J,N1)
1756 itsk = omp_get_thread_num()
1757 n1 = nddl-idlft1
1758 f_ddl = idlft1+1+itsk*n1/ nthread
1759 l_ddl = idlft1+(itsk+1)*n1/ nthread
1760C
1761C----------------------
1762 CALL my_barrier
1763C---------------------
1764 DO i=f_ddl,l_ddl
1765 diag_c(i-idlft1) = diag_m(i)
1766 DO j=iadk(i),iadk(i+1)-1
1767 lt_c(j-k)=lt_m(j)
1768 ENDDO
1769 ENDDO
1770C----------------------
1771 CALL my_barrier
1772C---------------------
1773C
1774C SMP dynamic parallel loop
1775C
1776!$OMP DO SCHEDULE(DYNAMIC,1)
1777 DO i=idlft0+1,nddl
1778 CALL sp_dim(i ,iadk ,jdik ,nc ,max_a ,max_l )
1779 CALL fsa_solv(
1780 1 nddl ,nc ,iadk ,jdik ,diag_k ,
1781 2 lt_k ,diag_m,lt_m ,diag_c,lt_c ,
1782 3 max_l ,idlft0,idlft1,nne ,i_chk ,
1783 4 iadk ,jdik ,i )
1784 ENDDO
1785
1786!$OMP END DO
1787!$OMP END PARALLEL
1788C
1789 RETURN
1790 END
1791!||====================================================================
1792!|| imp_fsa_inv2hp ../engine/source/implicit/imp_fsa_inv.F
1793!||--- called by ------------------------------------------------------
1794!|| lin_solvih2 ../engine/source/implicit/lin_solv.F
1795!||--- calls -----------------------------------------------------
1796!|| fsa_solv ../engine/source/implicit/imp_fsa_inv.F
1797!|| imp_kfiltr ../engine/source/implicit/imp_fsa_inv.F
1798!|| my_barrier ../engine/source/system/machine.F
1799!|| omp_get_thread_num ../engine/source/engine/openmp_stub.F90
1800!|| sp_dim ../engine/source/implicit/imp_fsa_inv.F
1801!||--- uses -----------------------------------------------------
1802!|| message_mod ../engine/share/message_module/message_mod.F
1803!||====================================================================
1804 SUBROUTINE imp_fsa_inv2hp(
1805 1 NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
1806 2 IADM ,JDIM ,DIAG_M ,LT_M ,MAXC ,
1807 3 MAX_A ,NNE ,IDLFT0 ,IDLFT1 ,MAX_D ,
1808 4 D_TOL ,P_MACH )
1809C-----------------------------------------------
1810C M o d u l e s
1811C-----------------------------------------------
1812 USE message_mod
1813C-----------------------------------------------
1814C I m p l i c i t T y p e s
1815C-----------------------------------------------
1816#include "implicit_f.inc"
1817C-----------------------------------------------
1818C C o m m o n B l o c k s
1819C-----------------------------------------------
1820#include "task_c.inc"
1821C-----------------------------------------------
1822C D u m m y A r g u m e n t s
1823C-----------------------------------------------
1824 INTEGER NDDL ,NNZ ,IADK(*),JDIK(*),MAXC ,MAX_A ,NNE,
1825 . IDLFT0 ,IDLFT1,MAX_D,IADM(*),JDIM(*)
1826C REAL
1827 my_real
1828 . diag_k(*), lt_k(*),diag_m(*), lt_m(*),d_tol ,p_mach
1829C-----------------------------------------------
1830C L o c a l V a r i a b l e s
1831C-----------------------------------------------
1832C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
1833 INTEGER I,J,K,M,N,NC,MAX_L,IERR,I_CHK,IER1,
1834 . ITSK,F_DDL,L_DDL,N1
1835 my_real
1836 . DIAG_C(NDDL-IDLFT1+1),LT_C(MAX_D+1)
1837 INTEGER OMP_GET_THREAD_NUM
1838 EXTERNAL OMP_GET_THREAD_NUM
1839C-----------------------------
1840 I_CHK = nne
1841 nne = 0
1842 IF ((idlft0+1)>nddl) RETURN
1843C
1844C--------------copy the useful part------
1845 k=iadk(idlft1+1)-1
1846!$OMP PARALLEL PRIVATE(ITSK,F_DDL,L_DDL,NC,MAX_L,N1,J,I)
1847 itsk = omp_get_thread_num()
1848 n1 = nddl-idlft1
1849 f_ddl = idlft1+1+itsk*n1/ nthread
1850 l_ddl = idlft1+(itsk+1)*n1/ nthread
1851 DO i=idlft1+1,nddl
1852 diag_c(i-idlft1) = diag_m(i)
1853 DO j=iadk(i),iadk(i+1)-1
1854 lt_c(j-k)=lt_m(j)
1855 ENDDO
1856 ENDDO
1857C----------------------
1858 CALL my_barrier
1859C---------------------
1860C SMP dynamic parallel loop
1861C
1862!$OMP DO SCHEDULE(DYNAMIC,1)
1863 DO i=idlft0+1,nddl
1864 CALL sp_dim(i ,iadm ,jdim ,nc ,max_a ,max_l )
1865 CALL fsa_solv(
1866 1 nddl ,nc ,iadk ,jdik ,diag_k ,
1867 2 lt_k ,diag_m,lt_m ,diag_c,lt_c ,
1868 3 max_l ,idlft0,idlft1,nne ,i_chk ,
1869 4 iadm ,jdim ,i )
1870 ENDDO
1871C
1872!$OMP END DO
1873!$OMP END PARALLEL
1874C
1875Citask0 IF (ITASK == 0 ) THEN
1876 k = idlft0+1
1877 IF (d_tol>zero)
1878 . CALL imp_kfiltr(k ,nddl ,iadm ,jdim ,diag_m ,
1879 . lt_m ,d_tol ,p_mach,diag_k)
1880Citask0 END IF
1881C
1882 RETURN
1883 END
1884
1885
my_real
#define my_real
Definition cppsort.cpp:32
if
if(complex_arithmetic) id
imp_fac_icj
subroutine imp_fac_icj(nddl, nnz, iadk, jdik, diag_k, lt_k, iadm, jdim, diag_m, lt_m, psi, nnzm, max_l, isky, li, nne)
Definition imp_fac_ic.F:36
imp_fsa_invp
subroutine imp_fsa_invp(nddl, nnz, iadk, jdik, diag_k, lt_k, diag_m, lt_m, maxc, max_a, nne, idlft0, idlft1, max_d)
Definition imp_fsa_inv.F:595
sp_stat0
subroutine sp_stat0(il, iadk, jdik, nc, jm)
Definition imp_fsa_inv.F:35
imp_fsa_invh
subroutine imp_fsa_invh(nddl, nnz, iadk, jdik, diag_k, lt_k, diag_m, lt_m, maxc, max_a, nne, idlft0, idlft1, max_d, itask)
Definition imp_fsa_inv.F:1348
imp_fsa_invp2
subroutine imp_fsa_invp2(nddl, iadk, jdik, diag_k, lt_k, iadm, jdim, diag_m, lt_m, maxc, max_a, nne, idlft0, idlft1, max_d, d_tol, p_mach)
Definition imp_fsa_inv.F:1023
intab0
integer function intab0(nic, ic, n)
Definition imp_fsa_inv.F:509
get_subsp
subroutine get_subsp(nddl, iadk, jdik, diag_k, lt_k, nc, iada, jdia, diag_a, lt_a, jm, maxa, idlft0, idlft1, diag_c, lt_c, diag_m, lt_m)
Definition imp_fsa_inv.F:372
get_subsa
subroutine get_subsa(nddl, iadk, jdik, diag_k, lt_k, nc, iada, jdia, diag_a, lt_a, jm)
Definition imp_fsa_inv.F:453
imp_pcg1
subroutine imp_pcg1(nddl, nnz, iadk, jdik, diag_k, lt_k, r, isp)
Definition imp_fsa_inv.F:693
imp_fsai
subroutine imp_fsai(n, iada, jdia, diag_a, lt_a, maxa, mj)
Definition imp_fsa_inv.F:267
get_subs0
subroutine get_subs0(nddl, iadk, jdik, diag_k, lt_k, nc, iada, jdia, diag_a, lt_a, jm, maxa)
Definition imp_fsa_inv.F:313
ind_lt2ln
subroutine ind_lt2ln(nddl, iadk, jdik, lt_k, maxl)
Definition imp_fsa_inv.F:1280
sp_dim
subroutine sp_dim(il, iadk, jdik, nc, max_a, max_l)
Definition imp_fsa_inv.F:1583
sp_static
subroutine sp_static(nddl, iadk, jdik, diag_k, lt_k, iadm, jdim, nnzm, nc, jm, maxc, psi, ip)
Definition imp_fsa_inv.F:70
intab2
integer function intab2(nic1, ic1, nic2, ic2)
Definition imp_fsa_inv.F:551
imp_fsa_inv2hp
subroutine imp_fsa_inv2hp(nddl, iadk, jdik, diag_k, lt_k, iadm, jdim, diag_m, lt_m, maxc, max_a, nne, idlft0, idlft1, max_d, d_tol, p_mach)
Definition imp_fsa_inv.F:1809
imp_fsa_invh2
subroutine imp_fsa_invh2(nddl, iadk, jdik, diag_k, lt_k, iadm, jdim, diag_m, lt_m, maxc, max_a, nne, idlft0, idlft1, max_d, d_tol, p_mach, itask)
Definition imp_fsa_inv.F:1471
imp_fsa_invhp
subroutine imp_fsa_invhp(nddl, nnz, iadk, jdik, diag_k, lt_k, diag_m, lt_m, maxc, max_a, nne, idlft0, idlft1, max_d)
Definition imp_fsa_inv.F:1718
fsa_solv
subroutine fsa_solv(nddl, nc, iadk, jdik, diag_k, lt_k, diag_m, lt_m, diag_c, lt_c, max_a, idlft0, idlft1, nne, i_chk, iadm, jdim, i)
Definition imp_fsa_inv.F:1626
imp_kfiltr
subroutine imp_kfiltr(ndf, nd, iada, jdia, diag_a, lt_a, tol, e_ps, diag_k)
Definition imp_fsa_inv.F:1126
sp_a2
subroutine sp_a2(nddl, nc, jm, maxc, ifsai)
Definition imp_fsa_inv.F:187
get_subsn
subroutine get_subsn(nddl, iadk, jdik, diag_k, lt_k, nc, iada, jdia, diag_a, lt_a, jm, maxa)
Definition imp_fsa_inv.F:1223
imp_fsa_inv2
subroutine imp_fsa_inv2(nddl, iadk, jdik, diag_k, lt_k, iadm, jdim, diag_m, lt_m, maxc, max_a, nne, d_tol, p_mach)
Definition imp_fsa_inv.F:922
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
message_mod
Definition message_mod.F:1257
prec0_solv
subroutine prec0_solv(nddl, nnz, iadm, jdim, diag_m, lt_m, v, z)
Definition prec_solv.F:239
produt_v0
subroutine produt_v0(nddl, x, y, r)
Definition produt_v.F:944
cp_real
subroutine cp_real(n, x, xc)
Definition produt_v.F:871
mav_lt
subroutine mav_lt(nddl, nnz, iadl, jdil, diag_k, lt_k, v, w)
Definition produt_v.F:340
cp_int
subroutine cp_int(n, x, xc)
Definition produt_v.F:916
sms_fsa_invh
subroutine sms_fsa_invh(nnzm, iadm, jdim, diag_m, lt_m, nndft0, nndft1, itask, diag_inv)
Definition sms_fsa_inv.F:344
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