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; unpaire:avec pre-filtre; 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------------avec 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 est deja en ordre croissante---------------------------7---------8
522 INTAB0=0
523.OR. IF (N<IC(1)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 dans IC1--7---------8
569 INTAB2=0
570.OR. IF (IC1(NIC1)<IC2(1)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 la partie utile------
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.AND. IF (I_CHK>0IERR<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------------filtrage----Diagonal est dans LT_M (last one)--
660 DIAG_M(I)=MJ(NC)
661 IF (DIAG_M(I)<EM20) THEN
662.AND. IF (NNE==0I_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.AND. IF (I_CHK>0MJ(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-------pour etre coherent avec 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 1002 FORMAT(3X,'total c.g. iteration=',I8,5X,
904 . ' relative residual norm=',E11.4)
905 1003 FORMAT(5X,
906 . '---warning : the iteration limit number was reached',
907 . 1X,'in preconditioner')
908 RETURN
909 END
910!||====================================================================
911!|| imp_fsa_inv2 ../engine/source/implicit/imp_fsa_inv.F
912!||--- calls -----------------------------------------------------
913!|| ancmsg ../engine/source/output/message/message.F
914!|| arret ../engine/source/system/arret.F
915!|| get_subs0 ../engine/source/implicit/imp_fsa_inv.F
916!|| imp_fsai ../engine/source/implicit/imp_fsa_inv.F
917!|| imp_kfiltr ../engine/source/implicit/imp_fsa_inv.F
918!|| imp_pcg1 ../engine/source/implicit/imp_fsa_inv.F
919!|| sp_stat0 ../engine/source/implicit/imp_fsa_inv.F
920!||--- uses -----------------------------------------------------
921!|| message_mod ../engine/share/message_module/message_mod.F
922!||====================================================================
923 SUBROUTINE IMP_FSA_INV2(
924 1 NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
925 2 IADM ,JDIM ,DIAG_M ,LT_M ,MAXC ,
926 3 MAX_A ,NNE ,D_TOL ,P_MACH)
927C-----------------------------------------------
928C M o d u l e s
929C-----------------------------------------------
930 USE MESSAGE_MOD
931C-----------------------------------------------
932C I m p l i c i t T y p e s
933C-----------------------------------------------
934#include "implicit_f.inc"
935C-----------------------------------------------
936C D u m m y A r g u m e n t s
937C-----------------------------------------------
938 INTEGER NDDL ,IADK(*),JDIK(*),MAXC ,MAX_A ,NNE,
939 . IADM(*),JDIM(*)
940C REAL
941 my_real
942 . DIAG_K(*), LT_K(*),DIAG_M(*), LT_M(*) ,D_TOL ,P_MACH
943C-----------------------------------------------
944C L o c a l V a r i a b l e s
945C-----------------------------------------------
946C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
947 INTEGER I,J,K,M,N,NC,IERR
948 INTEGER MAX_L,I_CHK
949 INTEGER, DIMENSION(:),ALLOCATABLE :: IADA,JDIA,JM
950 my_real, DIMENSION(:),ALLOCATABLE :: DIAG_A,LT_A,MJ
951C-----------------------------
952 I_CHK = NNE
953 NNE = 0
954 ALLOCATE(IADA(MAXC+1))
955 ALLOCATE(JDIA(MAX_A))
956 ALLOCATE(JM(MAXC+1))
957 ALLOCATE(DIAG_A(MAXC))
958 ALLOCATE(MJ(MAXC))
959 ALLOCATE(LT_A(MAX_A),STAT=IERR)
960 IF (IERR/=0) THEN
961 CALL ANCMSG(MSGID=19,ANMODE=ANINFO,
962 . C1='for IMPLICIT precondition')
963 CALL ARRET(2)
964 ENDIF
965 DO I=1,NDDL
966 CALL SP_STAT0(I ,IADM ,JDIM ,NC ,JM )
967 CALL GET_SUBS0(NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
968 . NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
969 . JM ,MAX_A )
970 DO J=1,NC-1
971 MJ(J)=ZERO
972 ENDDO
973 MJ(NC)=ONE
974 IF (NC>10000) THEN
975 MAX_L=MAX_A
976 CALL IMP_PCG1(
977 1 NC ,MAX_L ,IADA ,JDIA ,DIAG_A ,
978 2 LT_A ,MJ ,IERR )
979.AND. IF (I_CHK>0IERR<0) NNE = I
980 ELSE
981 MAX_L=1+(NC*(NC-1))/2
982 CALL IMP_FSAI(NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
983 . MAX_L ,MJ )
984 ENDIF
985C-----------------
986 DIAG_M(I)=MJ(NC)
987 IF (DIAG_M(I)<EM20) THEN
988.AND. IF (NNE==0I_CHK==0) NNE = I
989 DIAG_M(I)=ABS(DIAG_M(I))
990 DIAG_M(I)=MAX(EM20,DIAG_M(I))
991 ENDIF
992 DO K =1,NC-1
993 M=IADM(I)+K-1
994 LT_M(M)=MJ(K)/DIAG_M(I)
995 ENDDO
996.AND. IF (I_CHK>0MJ(NC)<EM20) DIAG_M(I)=MJ(NC)
997 ENDDO
998C
999 DEALLOCATE(IADA,JDIA)
1000 DEALLOCATE(JM)
1001 DEALLOCATE(DIAG_A,LT_A)
1002 DEALLOCATE(MJ)
1003 K = 1
1004 IF (D_TOL>ZERO)
1005 . CALL IMP_KFILTR(K ,NDDL ,IADM ,JDIM ,DIAG_M ,
1006 . LT_M ,D_TOL ,P_MACH,DIAG_K)
1007C
1008 RETURN
1009 END
1010!||====================================================================
1011!|| imp_fsa_invp2 ../engine/source/implicit/imp_fsa_inv.F
1012!||--- calls -----------------------------------------------------
1013!|| ancmsg ../engine/source/output/message/message.F
1014!|| arret ../engine/source/system/arret.F
1015!|| get_subsp ../engine/source/implicit/imp_fsa_inv.F
1016!|| imp_fsai ../engine/source/implicit/imp_fsa_inv.F
1017!|| imp_kfiltr ../engine/source/implicit/imp_fsa_inv.F
1018!|| imp_pcg1 ../engine/source/implicit/imp_fsa_inv.F
1019!|| sp_stat0 ../engine/source/implicit/imp_fsa_inv.F
1020!||--- uses -----------------------------------------------------
1021!|| message_mod ../engine/share/message_module/message_mod.F
1022!||====================================================================
1023 SUBROUTINE IMP_FSA_INVP2(
1024 1 NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
1025 2 IADM ,JDIM ,DIAG_M ,LT_M ,MAXC ,
1026 3 MAX_A ,NNE ,IDLFT0 ,IDLFT1 ,MAX_D ,
1027 4 D_TOL ,P_MACH)
1028C-----------------------------------------------
1029C M o d u l e s
1030C-----------------------------------------------
1031 USE MESSAGE_MOD
1032C-----------------------------------------------
1033C I m p l i c i t T y p e s
1034C-----------------------------------------------
1035#include "implicit_f.inc"
1036C-----------------------------------------------
1037C D u m m y A r g u m e n t s
1038C-----------------------------------------------
1039 INTEGER NDDL ,NNZ ,IADK(*),JDIK(*),MAXC ,MAX_A ,NNE,
1040 . IDLFT0 ,IDLFT1,MAX_D,IADM(*),JDIM(*)
1041C REAL
1042 my_real
1043 . DIAG_K(*), LT_K(*),DIAG_M(*), LT_M(*),D_TOL ,P_MACH
1044C-----------------------------------------------
1045C L o c a l V a r i a b l e s
1046C-----------------------------------------------
1047C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
1048 INTEGER I,J,K,M,N,NC,IADA(MAXC+1),JDIA(MAX_A),JM(MAXC+1)
1049 INTEGER MAX_L,IERR,I_CHK
1050 my_real
1051 . DIAG_A(MAXC),MJ(MAXC),
1052 . DIAG_C(NDDL-IDLFT1+1),LT_C(MAX_D+1)
1053 my_real, DIMENSION(:),ALLOCATABLE :: LT_A
1054C-----------------------------
1055 ALLOCATE(LT_A(MAX_A),STAT=IERR)
1056 IF (IERR/=0) THEN
1057 CALL ANCMSG(MSGID=19,ANMODE=ANINFO,
1058 . C1='for IMPLICIT precondition')
1059 CALL ARRET(2)
1060 ENDIF
1061C--------------copy la partie utile------
1062 I_CHK = NNE
1063 NNE = 0
1064 K=IADK(IDLFT1+1)-1
1065 DO I=IDLFT1+1,NDDL
1066 DIAG_C(I-IDLFT1) = DIAG_M(I)
1067 DO J=IADK(I),IADK(I+1)-1
1068 LT_C(J-K)=LT_M(J)
1069 ENDDO
1070 ENDDO
1071 DO I=IDLFT0+1,NDDL
1072 CALL SP_STAT0(I ,IADM ,JDIM ,NC ,JM )
1073 CALL GET_SUBSP(NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
1074 . NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
1075 . JM ,MAX_A ,IDLFT0,IDLFT1 ,DIAG_C,
1076 . LT_C ,DIAG_M,LT_M )
1077 DO J=1,NC-1
1078 MJ(J)=ZERO
1079 ENDDO
1080 MJ(NC)=ONE
1081 IF (NC>10000) THEN
1082 MAX_L=MAX_A
1083 CALL IMP_PCG1(
1084 1 NC ,MAX_L ,IADA ,JDIA ,DIAG_A ,
1085 2 LT_A ,MJ ,IERR )
1086.AND. IF (I_CHK>0IERR<0) NNE = I
1087 ELSE
1088 MAX_L=1+(NC*(NC-1))/2
1089 CALL IMP_FSAI(NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
1090 . MAX_L ,MJ )
1091 ENDIF
1092C------------Diagonal est dans LT_M (last one)--
1093 DIAG_M(I)=MJ(NC)
1094 IF (DIAG_M(I)<EM20) THEN
1095.AND. IF (NNE==0I_CHK==0) NNE = I
1096 DIAG_M(I)=ABS(DIAG_M(I))
1097 DIAG_M(I)=MAX(EM20,DIAG_M(I))
1098 ENDIF
1099 DO K =1,NC-1
1100 M=IADM(I)+K-1
1101 LT_M(M)=MJ(K)/DIAG_M(I)
1102 ENDDO
1103.AND. IF (I_CHK>0MJ(NC)<EM20) DIAG_M(I)=MJ(NC)
1104 ENDDO
1105C
1106 DEALLOCATE(LT_A)
1107 K = IDLFT0+1
1108 IF (D_TOL>ZERO)
1109 . CALL IMP_KFILTR(K ,NDDL ,IADM ,JDIM ,DIAG_M ,
1110 . LT_M ,D_TOL ,P_MACH,DIAG_K)
1111C
1112 RETURN
1113 END
1114!||====================================================================
1115!|| imp_kfiltr ../engine/source/implicit/imp_fsa_inv.F
1116!||--- called by ------------------------------------------------------
1117!|| imp_fsa_inv2 ../engine/source/implicit/imp_fsa_inv.F
1118!|| imp_fsa_inv2hp ../engine/source/implicit/imp_fsa_inv.F
1119!|| imp_fsa_invh2 ../engine/source/implicit/imp_fsa_inv.F
1120!|| imp_fsa_invp2 ../engine/source/implicit/imp_fsa_inv.F
1121!||--- calls -----------------------------------------------------
1122!|| ancmsg ../engine/source/output/message/message.F
1123!|| arret ../engine/source/system/arret.F
1124!|| cp_int ../engine/source/implicit/produt_v.F
1125!|| cp_real ../engine/source/implicit/produt_v.F
1126!||--- uses -----------------------------------------------------
1127!|| message_mod ../engine/share/message_module/message_mod.F
1128!||====================================================================
1129 SUBROUTINE IMP_KFILTR(NDF ,ND ,IADA ,JDIA ,DIAG_A ,
1130 . LT_A ,TOL ,E_PS ,DIAG_K )
1131C-----------------------------------------------
1132C M o d u l e s
1133C-----------------------------------------------
1134 USE MESSAGE_MOD
1135C-----------------------------------------------
1136C I m p l i c i t T y p e s
1137C-----------------------------------------------
1138#include "implicit_f.inc"
1139C-----------------------------------------------
1140C D u m m y A r g u m e n t s
1141C-----------------------------------------------
1142 INTEGER NDF,ND ,IADA(*) ,JDIA(*)
1143C REAL
1144 my_real
1145 . DIAG_A(*),LT_A(*),TOL,E_PS,DIAG_K(*)
1146C-----------------------------------------------
1147C L o c a l V a r i a b l e s
1148C-----------------------------------------------
1149 INTEGER I,J,K,NZ,IERR,MNZ,INORM
1150 INTEGER, DIMENSION(:),ALLOCATABLE :: IADL,JDIL
1151 my_real
1152 . MIN_D,MAX_D,MTOL,DD,TAUX
1153 my_real, DIMENSION(:),ALLOCATABLE :: LT_L
1154C-----------------------------
1155 print *,'d_tol,p_mach=',tol,e_ps
1156 NZ = IADA(ND+1)-IADA(1)
1157 DO I = 1, NDF-1
1158 DIAG_A(I) = ZERO
1159 DO J = IADA(I), IADA(I+1)-1
1160 LT_A(J) = ZERO
1161 ENDDO
1162 IADA(I+1) = 1
1163 ENDDO
1164.AND. IF (NZ>0TOL>ZERO) THEN
1165 ALLOCATE(IADL(ND+1))
1166 ALLOCATE(JDIL(NZ),LT_L(NZ),STAT=IERR)
1167 IF (IERR/=0) THEN
1168 CALL ANCMSG(MSGID=19,ANMODE=ANINFO,
1169 . C1='for IMPLICIT precondition')
1170 CALL ARRET(2)
1171 ENDIF
1172C-----------------------------------------------
1173 CALL CP_INT(ND+1,IADA,IADL)
1174 CALL CP_INT(NZ,JDIA,JDIL)
1175 CALL CP_REAL(NZ,LT_A,LT_L)
1176 MAX_D = ZERO
1177 MIN_D = EP20
1178 DO I = NDF, ND
1179 MAX_D = MAX(MAX_D,DIAG_A(I))
1180 MIN_D = MIN(MIN_D,DIAG_A(I))
1181 ENDDO
1182 print *,'max_d,min_d=',max_d,min_d
1183c MTOL = TOL*MIN_D
1184c MTOL = MAX(E_PS,MTOL)
1185C------post-filtration------
1186 MNZ = 0
1187 DO I = NDF, ND
1188 DO J = IADL(I), IADL(I+1)-1
1189 MTOL = TOL*MIN(DIAG_A(JDIL(J)),DIAG_A(I))
1190 MTOL = MAX(E_PS,MTOL)
1191 IF (ABS(LT_L(J))>MTOL) THEN
1192 MNZ = MNZ + 1
1193 JDIA(MNZ) = JDIL(J)
1194 LT_A(MNZ) = LT_L(J)
1195 ENDIF
1196 ENDDO
1197 IADA(I+1) = MNZ + 1
1198 ENDDO
1199 DEALLOCATE(IADL,JDIL)
1200 DEALLOCATE(LT_L)
1201 TAUX = ONE*MNZ/NZ
1202c print *,'mnz,nz=',mnz,nz,mtol
1203 print *,'filtrage factor=',TAUX
1204C------retrun from M normalized------
1205 ENDIF
1206 INORM = 0
1207 IF (INORM>ZERO) THEN
1208 DO I = NDF, ND
1209 DIAG_A(I) = DIAG_A(I)/DIAG_K(I)
1210 DO J = IADA(I), IADA(I+1)-1
1211 DD = SQRT(DIAG_K(I)/DIAG_K(JDIA(J)))
1212 LT_A(J) = DD*LT_A(J)
1213 ENDDO
1214 ENDDO
1215 ENDIF
1216C--------------------------------------------
1217 RETURN
1218 END
1219!||====================================================================
1220!|| get_subsn ../engine/source/implicit/imp_fsa_inv.F
1221!||--- calls -----------------------------------------------------
1222!|| ind_lt2ln ../engine/source/implicit/imp_fsa_inv.F
1223!|| intab0 ../engine/source/implicit/imp_fsa_inv.F
1224!||====================================================================
1225 SUBROUTINE GET_SUBSN(NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
1226 . NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
1227 . JM ,MAXA )
1228C-----------------------------------------------
1229C I m p l i c i t T y p e s
1230C-----------------------------------------------
1231#include "implicit_f.inc"
1232C-----------------------------------------------
1233C D u m m y A r g u m e n t s
1234C-----------------------------------------------
1235 INTEGER NDDL ,IADK(*) ,JDIK(*),IADA(*) ,JDIA(*)
1236 INTEGER NC ,JM(*),MAXA
1237C REAL
1238 my_real
1239 . LT_K(*),DIAG_K(*),LT_A(*),DIAG_A(*)
1240C-----------------------------------------------
1241C External function
1242C-----------------------------------------------
1243 INTEGER INTAB0
1244 EXTERNAL INTAB0
1245C-----------------------------------------------
1246C L o c a l V a r i a b l e s
1247C-----------------------------------------------
1248 INTEGER I,J,K,JJ,NNZA,N
1249 my_real
1250 . DD
1251C--------------------------------------------
1252 NNZA=0
1253 IADA(1)=1
1254 DO I=1,NC
1255 J=JM(I)
1256 DIAG_A(I)=ONE
1257 DO K=IADK(J),IADK(J+1)-1
1258 JJ=JDIK(K)
1259 N=INTAB0(NC,JM,JJ)
1260 IF (N>0) THEN
1261 DD = SQRT(DIAG_K(J)*DIAG_K(JJ))
1262 NNZA=NNZA+1
1263 JDIA(NNZA)=N
1264 LT_A(NNZA)=LT_K(K)/DD
1265 ENDIF
1266 ENDDO
1267 IADA(I+1)=NNZA+1
1268 ENDDO
1269 CALL IND_LT2LN(NC,IADA ,JDIA ,LT_A, NNZA )
1270C
1271 RETURN
1272 END
1273!||====================================================================
1274!|| ind_lt2ln ../engine/source/implicit/imp_fsa_inv.F
1275!||--- called by ------------------------------------------------------
1276!|| get_subs0 ../engine/source/implicit/imp_fsa_inv.F
1277!|| get_subsn ../engine/source/implicit/imp_fsa_inv.F
1278!|| get_subsp ../engine/source/implicit/imp_fsa_inv.F
1279!|| get_subsp_sms ../engine/source/ams/sms_fsa_inv.F
1280!||--- calls -----------------------------------------------------
1281!|| cp_int ../engine/source/implicit/produt_v.F
1282!|| cp_real ../engine/source/implicit/produt_v.F
1283!||====================================================================
1284 SUBROUTINE IND_LT2LN(NDDL,IADK ,JDIK ,LT_K, MAXL )
1285C-------------
1286C-----------------------------------------------
1287C I m p l i c i t T y p e s
1288C-----------------------------------------------
1289#include "implicit_f.inc"
1290C-----------------------------------------------
1291C D u m m y A r g u m e n t s
1292C-----------------------------------------------
1293 INTEGER NDDL, IADK(*),JDIK(*),MAXL
1294 my_real
1295 . LT_K(*)
1296C REAL
1297C-----------------------------------------------
1298C L o c a l V a r i a b l e s
1299C-----------------------------------------------
1300 INTEGER IADM(NDDL+1),JDIM(MAXL),ICOL(NDDL)
1301 INTEGER I,JD,J,K,N,NM
1302 my_real
1303 . LT_M(MAXL)
1304C-----------------------------------------------
1305 CALL CP_INT(NDDL+1,IADK,IADM)
1306 CALL CP_INT(MAXL,JDIK,JDIM)
1307 CALL CP_REAL(MAXL,LT_K,LT_M)
1308C
1309 DO I = 1, NDDL
1310 ICOL(I) = 0
1311 DO J = IADM(I),IADM(I+1)-1
1312 JD = JDIM(J)
1313 ICOL(JD) = ICOL(JD) + 1
1314 ENDDO
1315 ENDDO
1316C
1317 IADK(1) = 1
1318 DO I = 1,NDDL
1319 IADK(I+1) = IADK(I)+ICOL(I)
1320 ENDDO
1321C
1322 DO I = 1,NDDL
1323 ICOL(I) = 0
1324 DO J=IADM(I),IADM(I+1)-1
1325 JD = JDIM(J)
1326 K = IADK(JD) + ICOL(JD)
1327 JDIK(K) = I
1328 LT_K(K) = LT_M(J)
1329 ICOL(JD) = ICOL(JD) + 1
1330 ENDDO
1331 ENDDO
1332C--------------------------------------------
1333 RETURN
1334 END
1335C ------factorized sparse approximate inverse version hybrid-------
1336!||====================================================================
1337!|| imp_fsa_invh ../engine/source/implicit/imp_fsa_inv.F
1338!||--- calls -----------------------------------------------------
1339!|| ancmsg ../engine/source/output/message/message.F
1340!|| arret ../engine/source/system/arret.F
1341!|| get_subsp ../engine/source/implicit/imp_fsa_inv.F
1342!|| imp_fsai ../engine/source/implicit/imp_fsa_inv.F
1343!|| imp_pcg1 ../engine/source/implicit/imp_fsa_inv.F
1344!|| my_barrier ../engine/source/system/machine.F
1345!|| sp_stat0 ../engine/source/implicit/imp_fsa_inv.F
1346!||--- uses -----------------------------------------------------
1347!|| message_mod ../engine/share/message_module/message_mod.F
1348!||====================================================================
1349 SUBROUTINE IMP_FSA_INVH(
1350 1 NDDL ,NNZ ,IADK ,JDIK ,DIAG_K ,
1351 2 LT_K ,DIAG_M ,LT_M ,MAXC ,MAX_A ,
1352 3 NNE ,IDLFT0 ,IDLFT1,MAX_D ,ITASK )
1353C-----------------------------------------------
1354C M o d u l e s
1355C-----------------------------------------------
1356 USE MESSAGE_MOD
1357C-----------------------------------------------
1358C I m p l i c i t T y p e s
1359C-----------------------------------------------
1360#include "implicit_f.inc"
1361C-----------------------------------------------
1362C D u m m y A r g u m e n t s
1363C-----------------------------------------------
1364 INTEGER NDDL ,NNZ ,IADK(*),JDIK(*),MAXC ,MAX_A ,NNE,
1365 . IDLFT0 ,IDLFT1,MAX_D,ITASK
1366C REAL
1367 my_real
1368 . DIAG_K(*), LT_K(*),DIAG_M(*), LT_M(*)
1369C-----------------------------------------------
1370C L o c a l V a r i a b l e s
1371C-----------------------------------------------
1372C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
1373 INTEGER I,J,K,M,N,NC,MAX_L,IERR,I_CHK,IER1,
1374 . JM(MAXC+1)
1375 my_real
1376 . DIAG_C(NDDL-IDLFT1+1),LT_C(MAX_D+1)
1377 INTEGER, DIMENSION(:),ALLOCATABLE :: IADA,JDIA
1378 my_real, DIMENSION(:),ALLOCATABLE :: DIAG_A,LT_A,MJ
1379C-----------------------------
1380 I_CHK = NNE
1381 NNE = 0
1382 IF ((IDLFT0+1)>NDDL) RETURN
1383C
1384 ALLOCATE(IADA(MAXC+1),DIAG_A(MAXC),MJ(MAXC),STAT=IER1)
1385 ALLOCATE(LT_A(MAX_A),JDIA(MAX_A),STAT=IERR)
1386
1387 IF ((IERR+IER1)/=0) THEN
1388 IF (ITASK == 0 ) THEN
1389 CALL ANCMSG(MSGID=19,ANMODE=ANINFO,
1390 . C1='for IMPLICIT precondition')
1391 CALL ARRET(2)
1392 END IF !(ITASK == 0 ) THEN
1393 ENDIF
1394
1395C--------------copy la partie utile------
1396 K=IADK(IDLFT1+1)-1
1397 DO I=IDLFT1+1,NDDL
1398 DIAG_C(I-IDLFT1) = DIAG_M(I)
1399 DO J=IADK(I),IADK(I+1)-1
1400 LT_C(J-K)=LT_M(J)
1401 ENDDO
1402 ENDDO
1403C----------------------
1404 CALL MY_BARRIER
1405C---------------------
1406C
1407C Boucle parallele dynamique SMP
1408C
1409!$OMP DO SCHEDULE(DYNAMIC,1)
1410 DO I=IDLFT0+1,NDDL
1411 CALL SP_STAT0(I ,IADK ,JDIK ,NC ,JM )
1412 CALL GET_SUBSP(NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
1413 . NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
1414 . JM ,MAX_A ,IDLFT0,IDLFT1 ,DIAG_C,
1415 . LT_C ,DIAG_M ,LT_M )
1416 DO J=1,NC-1
1417 MJ(J)=ZERO
1418 ENDDO
1419 MJ(NC)=ONE
1420C
1421 IF (NC>10000) THEN
1422 MAX_L=MAX_A
1423 CALL IMP_PCG1(
1424 1 NC ,MAX_L ,IADA ,JDIA ,DIAG_A ,
1425 2 LT_A ,MJ ,IERR )
1426C
1427.AND. IF (I_CHK>0IERR<0) NNE = I
1428 ELSE
1429C
1430 MAX_L=1+(NC*(NC-1))/2
1431 CALL IMP_FSAI(NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
1432 . MAX_L ,MJ )
1433C
1434 ENDIF
1435C------------filtrage----Diagonal est dans LT_M (last one)--
1436 DIAG_M(I)=MJ(NC)
1437 IF (DIAG_M(I)<EM20) THEN
1438.AND. IF (NNE==0I_CHK==0) NNE = I
1439 DIAG_M(I)=ABS(DIAG_M(I))
1440 DIAG_M(I)=MAX(EM20,DIAG_M(I))
1441 ENDIF
1442 DO K =1,NC-1
1443 M=IADK(I)+K-1
1444 LT_M(M)=MJ(K)/DIAG_M(I)
1445 ENDDO
1446C
1447.AND. IF (I_CHK>0MJ(NC)<EM20) DIAG_M(I)=MJ(NC)
1448 ENDDO
1449
1450!$OMP END DO
1451C
1452 DEALLOCATE(IADA,DIAG_A,MJ)
1453 DEALLOCATE(LT_A,JDIA)
1454C
1455 RETURN
1456 END
1457!||====================================================================
1458!|| imp_fsa_invh2 ../engine/source/implicit/imp_fsa_inv.F
1459!||--- calls -----------------------------------------------------
1460!|| ancmsg ../engine/source/output/message/message.F
1461!|| arret ../engine/source/system/arret.F
1462!|| get_subsp ../engine/source/implicit/imp_fsa_inv.F
1463!|| imp_fsai ../engine/source/implicit/imp_fsa_inv.F
1464!|| imp_kfiltr ../engine/source/implicit/imp_fsa_inv.F
1465!|| imp_pcg1 ../engine/source/implicit/imp_fsa_inv.F
1466!|| my_barrier ../engine/source/system/machine.F
1467!|| sp_stat0 ../engine/source/implicit/imp_fsa_inv.F
1468!||--- uses -----------------------------------------------------
1469!|| message_mod ../engine/share/message_module/message_mod.F
1470!||====================================================================
1471 SUBROUTINE IMP_FSA_INVH2(
1472 1 NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
1473 2 IADM ,JDIM ,DIAG_M ,LT_M ,MAXC ,
1474 3 MAX_A ,NNE ,IDLFT0 ,IDLFT1 ,MAX_D ,
1475 4 D_TOL ,P_MACH ,ITASK )
1476C-----------------------------------------------
1477C M o d u l e s
1478C-----------------------------------------------
1479 USE MESSAGE_MOD
1480C-----------------------------------------------
1481C I m p l i c i t T y p e s
1482C-----------------------------------------------
1483#include "implicit_f.inc"
1484C-----------------------------------------------
1485C D u m m y A r g u m e n t s
1486C-----------------------------------------------
1487 INTEGER NDDL ,NNZ ,IADK(*),JDIK(*),MAXC ,MAX_A ,NNE,
1488 . IDLFT0 ,IDLFT1,MAX_D,IADM(*),JDIM(*),ITASK
1489C REAL
1490 my_real
1491 . DIAG_K(*), LT_K(*),DIAG_M(*), LT_M(*),D_TOL ,P_MACH
1492C-----------------------------------------------
1493C L o c a l V a r i a b l e s
1494C-----------------------------------------------
1495C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
1496 INTEGER I,J,K,M,N,NC,MAX_L,IERR,I_CHK,IER1,
1497 . JM(MAXC+1)
1498 my_real
1499 . DIAG_C(NDDL-IDLFT1+1),LT_C(MAX_D+1)
1500 INTEGER, DIMENSION(:),ALLOCATABLE :: IADA,JDIA
1501 my_real, DIMENSION(:),ALLOCATABLE :: DIAG_A,LT_A,MJ
1502C-----------------------------
1503 I_CHK = NNE
1504 NNE = 0
1505 IF ((IDLFT0+1)>NDDL) RETURN
1506C
1507 ALLOCATE(IADA(MAXC+1),DIAG_A(MAXC),MJ(MAXC),STAT=IER1)
1508 ALLOCATE(LT_A(MAX_A),JDIA(MAX_A),STAT=IERR)
1509C--------
1510 IF ((IERR+IER1)/=0) THEN
1511 IF (ITASK == 0 ) THEN
1512 CALL ANCMSG(MSGID=19,ANMODE=ANINFO,
1513 . C1='for IMPLICIT precondition')
1514 CALL ARRET(2)
1515 END IF !(ITASK == 0 ) THEN
1516 ENDIF
1517C--------------copy la partie utile------
1518 K=IADK(IDLFT1+1)-1
1519 DO I=IDLFT1+1,NDDL
1520 DIAG_C(I-IDLFT1) = DIAG_M(I)
1521 DO J=IADK(I),IADK(I+1)-1
1522 LT_C(J-K)=LT_M(J)
1523 ENDDO
1524 ENDDO
1525C----------------------
1526 CALL MY_BARRIER
1527C---------------------
1528C Boucle parallele dynamique SMP
1529C
1530!$OMP DO SCHEDULE(DYNAMIC,1)
1531 DO I=IDLFT0+1,NDDL
1532 CALL SP_STAT0(I ,IADM ,JDIM ,NC ,JM )
1533 CALL GET_SUBSP(NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
1534 . NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
1535 . JM ,MAX_A ,IDLFT0,IDLFT1 ,DIAG_C,
1536 . LT_C ,DIAG_M,LT_M )
1537 DO J=1,NC-1
1538 MJ(J)=ZERO
1539 ENDDO
1540 MJ(NC)=ONE
1541 IF (NC>10000) THEN
1542 MAX_L=MAX_A
1543 CALL IMP_PCG1(
1544 1 NC ,MAX_L ,IADA ,JDIA ,DIAG_A ,
1545 2 LT_A ,MJ ,IERR )
1546.AND. IF (I_CHK>0IERR<0) NNE = I
1547 ELSE
1548 MAX_L=1+(NC*(NC-1))/2
1549 CALL IMP_FSAI(NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
1550 . MAX_L ,MJ )
1551 ENDIF
1552C------------Diagonal est dans LT_M (last one)--
1553 DIAG_M(I)=MJ(NC)
1554 IF (DIAG_M(I)<EM20) THEN
1555.AND. IF (NNE==0I_CHK==0) NNE = I
1556 DIAG_M(I)=ABS(DIAG_M(I))
1557 DIAG_M(I)=MAX(EM20,DIAG_M(I))
1558 ENDIF
1559 DO K =1,NC-1
1560 M=IADM(I)+K-1
1561 LT_M(M)=MJ(K)/DIAG_M(I)
1562 ENDDO
1563C
1564.AND. IF (I_CHK>0MJ(NC)<EM20) DIAG_M(I)=MJ(NC)
1565 ENDDO
1566C
1567!$OMP END DO
1568C
1569 DEALLOCATE(IADA,DIAG_A,MJ)
1570 DEALLOCATE(LT_A,JDIA)
1571C
1572 IF (ITASK == 0 ) THEN
1573 K = IDLFT0+1
1574 IF (D_TOL>ZERO)
1575 . CALL IMP_KFILTR(K ,NDDL ,IADM ,JDIM ,DIAG_M ,
1576 . LT_M ,D_TOL ,P_MACH,DIAG_K)
1577 END IF
1578C
1579 RETURN
1580 END
1581!||====================================================================
1582!|| sp_dim ../engine/source/implicit/imp_fsa_inv.F
1583!||--- called by ------------------------------------------------------
1584!|| imp_fsa_inv2hp ../engine/source/implicit/imp_fsa_inv.F
1585!|| imp_fsa_invhp ../engine/source/implicit/imp_fsa_inv.F
1586!||====================================================================
1587 SUBROUTINE SP_DIM(IL ,IADK ,JDIK ,NC ,MAX_A ,MAX_L )
1588C-----------------------------------------------
1589C I m p l i c i t T y p e s
1590C-----------------------------------------------
1591#include "implicit_f.inc"
1592C-----------------------------------------------
1593C D u m m y A r g u m e n t s
1594C-----------------------------------------------
1595 INTEGER IL, IADK(*) ,JDIK(*),NC,MAX_A ,MAX_L
1596C REAL
1597C-----------------------------------------------
1598C L o c a l V a r i a b l e s
1599C-----------------------------------------------
1600 INTEGER I,J,K
1601C-----------------------------------------------
1602 NC=0
1603 DO K =IADK(IL),IADK(IL+1)-1
1604 NC=NC+1
1605 ENDDO
1606 NC=NC+1
1607C------MAX_L<- MAX_A
1608 IF (NC <= 10000) THEN
1609 MAX_L = 1+(NC*(NC-1))/2
1610 ELSE
1611 MAX_L = MAX_A
1612 END IF
1613C--------------------------------------------
1614 RETURN
1615 END
1616!||====================================================================
1617!|| fsa_solv ../engine/source/implicit/imp_fsa_inv.F
1618!||--- called by ------------------------------------------------------
1619!|| imp_fsa_inv2hp ../engine/source/implicit/imp_fsa_inv.F
1620!|| imp_fsa_invhp ../engine/source/implicit/imp_fsa_inv.F
1621!||--- calls -----------------------------------------------------
1622!|| get_subsp ../engine/source/implicit/imp_fsa_inv.F
1623!|| imp_fsai ../engine/source/implicit/imp_fsa_inv.F
1624!|| imp_pcg1 ../engine/source/implicit/imp_fsa_inv.F
1625!||====================================================================
1626 SUBROUTINE FSA_SOLV(
1627 1 NDDL ,NC ,IADK ,JDIK ,DIAG_K ,
1628 2 LT_K ,DIAG_M,LT_M ,DIAG_C,LT_C ,
1629 3 MAX_A ,IDLFT0,IDLFT1,NNE ,I_CHK ,
1630 4 IADM ,JDIM ,I )
1631C-----------------------------------------------
1632C I m p l i c i t T y p e s
1633C-----------------------------------------------
1634#include "implicit_f.inc"
1635C-----------------------------------------------
1636C D u m m y A r g u m e n t s
1637C-----------------------------------------------
1638 INTEGER I,NDDL ,NC ,IADK(*),JDIK(*),MAX_A ,NNE,
1639 . IDLFT0,IDLFT1 ,I_CHK,IADM(*),JDIM(*)
1640C REAL
1641 my_real
1642 . DIAG_K(*), LT_K(*),DIAG_M(*), LT_M(*) ,DIAG_C(*),LT_C(*)
1643C-----------------------------------------------
1644C L o c a l V a r i a b l e s
1645C-----------------------------------------------
1646C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
1647 INTEGER J,K,M,N,MAX_L,IERR,IER1,IADA(NC+1),JM(NC)
1648c . IADA(NC+1),JDIA(MAX_A),JM(NC)
1649c my_real
1650c . DIAG_A(NC),LT_A(MAX_A),MJ(NC)
1651 INTEGER, DIMENSION(:),ALLOCATABLE :: JDIA
1652 my_real, DIMENSION(:),ALLOCATABLE :: DIAG_A,LT_A,MJ
1653C-----------------------------
1654C
1655 ALLOCATE(DIAG_A(NC),MJ(NC),STAT=IER1)
1656 ALLOCATE(LT_A(MAX_A),JDIA(MAX_A),STAT=IERR)
1657C-----------------------------
1658 J=0
1659 DO K =IADM(I),IADM(I+1)-1
1660 J=J+1
1661 JM(J)=JDIM(K)
1662 ENDDO
1663 J=J+1
1664 JM(J)=I
1665C
1666 CALL GET_SUBSP(NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
1667 . NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
1668 . JM ,MAX_A ,IDLFT0,IDLFT1 ,DIAG_C,
1669 . LT_C ,DIAG_M ,LT_M )
1670 DO J=1,NC-1
1671 MJ(J)=ZERO
1672 ENDDO
1673 MJ(NC)=ONE
1674C
1675 IF (NC > 10000) THEN
1676 MAX_L=MAX_A
1677 CALL IMP_PCG1(
1678 1 NC ,MAX_L ,IADA ,JDIA ,DIAG_A ,
1679 2 LT_A ,MJ ,IERR )
1680C
1681.AND. IF (I_CHK>0IERR<0) NNE = I
1682 ELSE
1683C
1684 MAX_L=1+(NC*(NC-1))/2
1685 CALL IMP_FSAI(NC ,IADA ,JDIA ,DIAG_A ,LT_A ,
1686 . MAX_L ,MJ )
1687C
1688 ENDIF
1689C------------filtrage----Diagonal est dans LT_M (last one)--
1690 DIAG_M(I)=MJ(NC)
1691 IF (DIAG_M(I)<EM20) THEN
1692.AND. IF (NNE==0I_CHK==0) NNE = I
1693 DIAG_M(I)=ABS(DIAG_M(I))
1694 DIAG_M(I)=MAX(EM20,DIAG_M(I))
1695 ENDIF
1696 DO K =1,NC-1
1697 M=IADM(I)+K-1
1698 LT_M(M)=MJ(K)/DIAG_M(I)
1699 ENDDO
1700.AND. IF (I_CHK>0MJ(NC)<EM20) DIAG_M(I)=MJ(NC)
1701C
1702 DEALLOCATE(DIAG_A,MJ)
1703 DEALLOCATE(LT_A,JDIA)
1704 RETURN
1705 END
1706C ------factorized sparse approximate inverse version hybrid SMP inside
1707!||====================================================================
1708!|| imp_fsa_invhp ../engine/source/implicit/imp_fsa_inv.F
1709!||--- called by ------------------------------------------------------
1710!|| lin_solvh1 ../engine/source/implicit/lin_solv.F
1711!||--- calls -----------------------------------------------------
1712!|| fsa_solv ../engine/source/implicit/imp_fsa_inv.F
1713!|| my_barrier ../engine/source/system/machine.F
1714!|| omp_get_thread_num ../engine/source/engine/openmp_stub.F90
1715!|| sp_dim ../engine/source/implicit/imp_fsa_inv.F
1716!||--- uses -----------------------------------------------------
1717!|| message_mod ../engine/share/message_module/message_mod.F
1718!||====================================================================
1719 SUBROUTINE IMP_FSA_INVHP(
1720 1 NDDL ,NNZ ,IADK ,JDIK ,DIAG_K ,
1721 2 LT_K ,DIAG_M ,LT_M ,MAXC ,MAX_A ,
1722 3 NNE ,IDLFT0 ,IDLFT1,MAX_D )
1723C-----------------------------------------------
1724C M o d u l e s
1725C-----------------------------------------------
1726 USE MESSAGE_MOD
1727C-----------------------------------------------
1728C I m p l i c i t T y p e s
1729C-----------------------------------------------
1730#include "implicit_f.inc"
1731C-----------------------------------------------
1732C C o m m o n B l o c k s
1733C-----------------------------------------------
1734#include "task_c.inc"
1735C-----------------------------------------------
1736C D u m m y A r g u m e n t s
1737C-----------------------------------------------
1738 INTEGER NDDL ,NNZ ,IADK(*),JDIK(*),MAXC ,MAX_A ,NNE,
1739 . IDLFT0 ,IDLFT1,MAX_D
1740C REAL
1741 my_real
1742 . DIAG_K(*), LT_K(*),DIAG_M(*), LT_M(*)
1743C-----------------------------------------------
1744C L o c a l V a r i a b l e s
1745C-----------------------------------------------
1746C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
1747 INTEGER I,J,K,M,N,NC,MAX_L,IERR,I_CHK,IER1,
1748 . ITSK,F_DDL,L_DDL,N1
1749 my_real
1750 . DIAG_C(NDDL-IDLFT1+1),LT_C(MAX_D+1)
1751 INTEGER OMP_GET_THREAD_NUM
1752 EXTERNAL OMP_GET_THREAD_NUM
1753C-----------------------------
1754 I_CHK = NNE
1755 NNE = 0
1756 IF ((IDLFT0+1)>NDDL) RETURN
1757C
1758C--------------copy la partie utile------
1759 K=IADK(IDLFT1+1)-1
1760!$OMP PARALLEL PRIVATE(ITSK,F_DDL,L_DDL,NC,MAX_L,I,J,N1)
1761 ITSK = OMP_GET_THREAD_NUM()
1762 N1 = NDDL-IDLFT1
1763 F_DDL = IDLFT1+1+ITSK*N1/ NTHREAD
1764 L_DDL = IDLFT1+(ITSK+1)*N1/ NTHREAD
1765C
1766C----------------------
1767 CALL MY_BARRIER
1768C---------------------
1769 DO I=F_DDL,L_DDL
1770 DIAG_C(I-IDLFT1) = DIAG_M(I)
1771 DO J=IADK(I),IADK(I+1)-1
1772 LT_C(J-K)=LT_M(J)
1773 ENDDO
1774 ENDDO
1775C----------------------
1776 CALL MY_BARRIER
1777C---------------------
1778C
1779C Boucle parallele dynamique SMP
1780C
1781!$OMP DO SCHEDULE(DYNAMIC,1)
1782 DO I=IDLFT0+1,NDDL
1783 CALL SP_DIM(I ,IADK ,JDIK ,NC ,MAX_A ,MAX_L )
1784 CALL FSA_SOLV(
1785 1 NDDL ,NC ,IADK ,JDIK ,DIAG_K ,
1786 2 LT_K ,DIAG_M,LT_M ,DIAG_C,LT_C ,
1787 3 MAX_L ,IDLFT0,IDLFT1,NNE ,I_CHK ,
1788 4 IADK ,JDIK ,I )
1789 ENDDO
1790
1791!$OMP END DO
1792!$OMP END PARALLEL
1793C
1794 RETURN
1795 END
1796!||====================================================================
1797!|| imp_fsa_inv2hp ../engine/source/implicit/imp_fsa_inv.F
1798!||--- called by ------------------------------------------------------
1799!|| lin_solvih2 ../engine/source/implicit/lin_solv.F
1800!||--- calls -----------------------------------------------------
1801!|| fsa_solv ../engine/source/implicit/imp_fsa_inv.F
1802!|| imp_kfiltr ../engine/source/implicit/imp_fsa_inv.F
1803!|| my_barrier ../engine/source/system/machine.F
1804!|| omp_get_thread_num ../engine/source/engine/openmp_stub.F90
1805!|| sp_dim ../engine/source/implicit/imp_fsa_inv.F
1806!||--- uses -----------------------------------------------------
1807!|| message_mod ../engine/share/message_module/message_mod.F
1808!||====================================================================
1809 SUBROUTINE IMP_FSA_INV2HP(
1810 1 NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
1811 2 IADM ,JDIM ,DIAG_M ,LT_M ,MAXC ,
1812 3 MAX_A ,NNE ,IDLFT0 ,IDLFT1 ,MAX_D ,
1813 4 D_TOL ,P_MACH )
1814C-----------------------------------------------
1815C M o d u l e s
1816C-----------------------------------------------
1817 USE MESSAGE_MOD
1818C-----------------------------------------------
1819C I m p l i c i t T y p e s
1820C-----------------------------------------------
1821#include "implicit_f.inc"
1822C-----------------------------------------------
1823C C o m m o n B l o c k s
1824C-----------------------------------------------
1825#include "task_c.inc"
1826C-----------------------------------------------
1827C D u m m y A r g u m e n t s
1828C-----------------------------------------------
1829 INTEGER NDDL ,NNZ ,IADK(*),JDIK(*),MAXC ,MAX_A ,NNE,
1830 . IDLFT0 ,IDLFT1,MAX_D,IADM(*),JDIM(*)
1831C REAL
1832 my_real
1833 . DIAG_K(*), LT_K(*),DIAG_M(*), LT_M(*),D_TOL ,P_MACH
1834C-----------------------------------------------
1835C L o c a l V a r i a b l e s
1836C-----------------------------------------------
1837C--- M->A^-1 LT_M strictly lower in m.c.r.s. format
1838 INTEGER I,J,K,M,N,NC,MAX_L,IERR,I_CHK,IER1,
1839 . ITSK,F_DDL,L_DDL,N1
1840 my_real
1841 . DIAG_C(NDDL-IDLFT1+1),LT_C(MAX_D+1)
1842 INTEGER OMP_GET_THREAD_NUM
1843 EXTERNAL OMP_GET_THREAD_NUM
1844C-----------------------------
1845 I_CHK = NNE
1846 NNE = 0
1847 IF ((IDLFT0+1)>NDDL) RETURN
1848C
1849C--------------copy la partie utile------
1850 K=IADK(IDLFT1+1)-1
1851!$OMP PARALLEL PRIVATE(ITSK,F_DDL,L_DDL,NC,MAX_L,N1,J,I)
1852 ITSK = OMP_GET_THREAD_NUM()
1853 N1 = NDDL-IDLFT1
1854 F_DDL = IDLFT1+1+ITSK*N1/ NTHREAD
1855 L_DDL = IDLFT1+(ITSK+1)*N1/ NTHREAD
1856 DO I=IDLFT1+1,NDDL
1857 DIAG_C(I-IDLFT1) = DIAG_M(I)
1858 DO J=IADK(I),IADK(I+1)-1
1859 LT_C(J-K)=LT_M(J)
1860 ENDDO
1861 ENDDO
1862C----------------------
1863 CALL MY_BARRIER
1864C---------------------
1865C Boucle parallele dynamique SMP
1866C
1867!$OMP DO SCHEDULE(DYNAMIC,1)
1868 DO I=IDLFT0+1,NDDL
1869 CALL SP_DIM(I ,IADM ,JDIM ,NC ,MAX_A ,MAX_L )
1870 CALL FSA_SOLV(
1871 1 NDDL ,NC ,IADK ,JDIK ,DIAG_K ,
1872 2 LT_K ,DIAG_M,LT_M ,DIAG_C,LT_C ,
1873 3 MAX_L ,IDLFT0,IDLFT1,NNE ,I_CHK ,
1874 4 IADM ,JDIM ,I )
1875 ENDDO
1876C
1877!$OMP END DO
1878!$OMP END PARALLEL
1879C
1880Citask0 IF (ITASK == 0 ) THEN
1881 K = IDLFT0+1
1882 IF (D_TOL>ZERO)
1883 . CALL IMP_KFILTR(K ,NDDL ,IADM ,JDIM ,DIAG_M ,
1884 . LT_M ,D_TOL ,P_MACH,DIAG_K)
1885Citask0 END IF
1886C
1887 RETURN
1888 END
1889
1890
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
sp_stat0
subroutine sp_stat0(il, iadk, jdik, nc, jm)
Definition imp_fsa_inv.F:35
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
sp_a2
subroutine sp_a2(nddl, nc, jm, maxc, ifsai)
Definition imp_fsa_inv.F:187
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91
arret
subroutine arret(nn)
Definition arret.F:87