OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
imp_pc_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/.
23C ------approximate inverse (minimisation:least squares by colonne)------
24!||====================================================================
25!|| imp_pc_inv ../engine/source/implicit/imp_pc_inv.F
26!||--- calls -----------------------------------------------------
27!|| arret ../engine/source/system/arret.F
28!|| get_suba ../engine/source/implicit/imp_pc_inv.F
29!|| imp_saic ../engine/source/implicit/imp_pc_inv.F
30!|| sp_static ../engine/source/implicit/imp_fsa_inv.F
31!||====================================================================
32 SUBROUTINE imp_pc_inv(
33 1 NDDL ,NNZ ,IADK ,JDIK ,DIAG_K ,
34 2 LT_K ,IADM ,JDIM ,DIAG_M, LT_M ,
35 3 PSI ,NNZM ,MAXC ,MAXA ,MAX_L ,
36 4 IOPT ,NNE )
37C-----------------------------------------------
38C I m p l i c i t T y p e s
39C-----------------------------------------------
40#include "implicit_f.inc"
41C-----------------------------------------------
42C D u m m y A r g u m e n t s
43C-----------------------------------------------
44 INTEGER NDDL ,NNZ ,IADK(*),JDIK(*),NNZM ,IADM(*),JDIM(*),
45 . MAXC ,MAXA ,MAX_L,IOPT,NNE
46C REAL
47 my_real
48 . diag_k(*), diag_m(*), lt_k(*) ,lt_m(*) ,psi
49C-----------------------------------------------
50C L o c a l V a r i a b l e s
51C-----------------------------------------------
52C--- DIAG_M->[D]^-1, LT_M strictly upper in c.r.s. format
53 INTEGER I,J,K,M,N,NC(NDDL),JM(MAXC,NDDL),I1
54 my_real
55 . PSR,TOL,MJ(MAXA),A(MAXA,MAXC)
56C-----------------------------
57 CALL sp_static(nddl ,iadk ,jdik ,diag_k ,lt_k ,
58 . iadm ,jdim ,nnzm ,nc ,jm ,
59 . maxc ,psi ,iopt )
60C------------parallelisation ici----
61 nnzm = 0
62 iadm(1)=1
63 nne=0
64 tol = psi*psi
65 DO i=1,nddl
66 CALL get_suba(nddl ,iadk ,jdik ,diag_k ,lt_k ,
67 . nc ,jm ,a ,maxc ,maxa ,
68 . i ,m ,mj )
69 IF (m>maxa) WRITE(*,*)'M>MAXB',m,maxa
70C write(*,*)'m,n,i=',m,nc(i),i
71 CALL imp_saic(m ,nc(i) ,a ,mj ,maxa )
72C------------post-filtrage----
73 DO k =1,nc(i)
74 j=jm(k,i)
75 IF (j==i) THEN
76 diag_m(i)=mj(k)
77 IF (diag_m(i)<=em20) nne=nne+1
78 ELSE
79 i1=min(i,j)
80 psr = tol*abs(diag_m(i1))
81 IF (abs(mj(k))>=psr) THEN
82 nnzm = nnzm+1
83 IF (nnzm>max_l) THEN
84 WRITE(*,*)'NNZM>MAX_L',nnzm,maxc,i
85 CALL arret(2)
86 ENDIF
87 jdim(nnzm)=j
88 lt_m(nnzm)=mj(k)
89 ENDIF
90 ENDIF
91 ENDDO
92 iadm(i+1)=nnzm+1
93 ENDDO
94C------------eventuellement deuxieme niveau----
95C
96 RETURN
97 END
98C-------------resol A(M,N).MJ-I=0 ----
99!||====================================================================
100!|| imp_saic ../engine/source/implicit/imp_pc_inv.F
101!||--- called by ------------------------------------------------------
102!|| imp_pc_inv ../engine/source/implicit/imp_pc_inv.F
103!||--- calls -----------------------------------------------------
104!|| imp_qrf ../engine/source/implicit/imp_pc_inv.F
105!|| lt_solv ../engine/source/implicit/imp_pc_inv.F
106!|| mav_qt ../engine/source/implicit/imp_pc_inv.f
107!||====================================================================
108 SUBROUTINE imp_saic(M ,N ,A ,MJ ,MAXC )
109C-----------------------------------------------
110C I m p l i c i t T y p e s
111C-----------------------------------------------
112#include "implicit_f.inc"
113C-----------------------------------------------
114C D u m m y A r g u m e n t s
115C-----------------------------------------------
116 INTEGER M,N ,MAXC
117C REAL
118 my_real
119 . A(MAXC,*),MJ(*)
120C-----------------------------------------------
121C L o c a l V a r i a b l e s
122C-----------------------------------------------
123 INTEGER I
124 my_real
125 . D_R(N),TAU(N)
126C-----------------------------------------------
127 CALL imp_qrf(m ,n ,a ,maxc ,d_r ,
128 . tau )
129C
130 CALL mav_qt(m ,n ,a ,maxc ,tau ,
131 . mj )
132C
133 CALL lt_solv(n ,a ,maxc ,d_r ,mj )
134C--------------------------------------------
135 RETURN
136 END
137C-------------set submatrix A(M,N) for SAI ----
138!||====================================================================
139!|| get_suba ../engine/source/implicit/imp_pc_inv.F
140!||--- called by ------------------------------------------------------
141!|| imp_pc_inv ../engine/source/implicit/imp_pc_inv.F
142!||--- calls -----------------------------------------------------
143!|| get_kijs ../engine/source/implicit/imp_pc_inv.F
144!|| intab2 ../engine/source/implicit/imp_fsa_inv.F
145!||====================================================================
146 SUBROUTINE get_suba(NDDL ,IADK ,JDIK ,DIAG_K ,LT_K ,
147 . NC ,JM ,A ,MAXC ,MAXA ,
148 . IM ,M ,MJ )
149C-----------------------------------------------
150C I m p l i c i t T y p e s
151C-----------------------------------------------
152#include "implicit_f.inc"
153C-----------------------------------------------
154C D u m m y A r g u m e n t s
155C-----------------------------------------------
156 INTEGER NDDL ,MAXC,MAXA,IADK(*) ,JDIK(*),IM,M
157 INTEGER NC(*),JM(MAXC,*)
158C REAL
159 my_real
160 . A(MAXA,*),LT_K(*),DIAG_K(*),MJ(*)
161C-----------------------------------------------
162C External function
163C-----------------------------------------------
164 INTEGER INTAB2
165 EXTERNAL intab2
166C-----------------------------------------------
167C L o c a l V a r i a b l e s
168C-----------------------------------------------
169 LOGICAL EXTN(NDDL)
170 INTEGER I,I1,J,K,K0,NL,N,JN,IUN
171C--------------------------------------------
172 m=0
173 DO i=1,nddl
174 extn(i)=.true.
175 ENDDO
176C------d'abord NxN-----------
177 DO i1=1,nc(im)
178 i=jm(i1,im)
179 extn(i)=.false.
180 m=m+1
181c write(*,*)'I1,IM,M=',I1,IM,M
182 IF (i==im) THEN
183 n=1
184 iun=m
185 100 j=jm(n,im)
186 IF (j<i) THEN
187 CALL get_kijs(j ,i ,iadk,jdik,lt_k ,a(m,n))
188 n=n+1
189 GOTO 100
190 ENDIF
191C------A(M,N)=DIAG_K(I)------
192 a(m,n)=diag_k(i)
193 n=n+1
194 DO k=n,nc(im)
195 j=jm(k,im)
196 CALL get_kijs(i ,j ,iadk,jdik,lt_k ,a(m,k))
197 ENDDO
198 ELSE
199 DO k=1,nc(im)
200 j=jm(k,im)
201 IF (j<i) THEN
202 CALL get_kijs(j ,i ,iadk,jdik,lt_k ,a(m,k))
203 ELSEIF (j==i) THEN
204 a(m,k)=diag_k(i)
205 ELSE
206 CALL get_kijs(i ,j ,iadk,jdik,lt_k ,a(m,k))
207 ENDIF
208 ENDDO
209 ENDIF
210 ENDDO
211C------ajoute autres lignes-----------
212 DO i=1,nddl
213 IF (extn(i)) THEN
214 n=intab2(nc(im),jm(1,im),nc(i),jm(1,i))
215 IF (n>0) THEN
216 IF(m==maxa) write(*,*)'mem',n,i
217 m=m+1
218 DO k=1,n-1
219 a(m,k)=zero
220 ENDDO
221 DO k=n,nc(im)
222 j=jm(k,im)
223 IF (j<i) THEN
224 CALL get_kijs(j ,i ,iadk,jdik,lt_k ,a(m,k))
225 ELSEIF (j==i) THEN
226 a(m,k)=diag_k(i)
227 ELSE
228 CALL get_kijs(i ,j ,iadk,jdik,lt_k ,a(m,k))
229 ENDIF
230 ENDDO
231 ENDIF
232 ENDIF
233 ENDDO
234C------INITIAL ej------
235 DO i=1,m
236 mj(i)=zero
237 ENDDO
238 mj(iun)=one
239C
240 RETURN
241 END
242C-------------QR factorisation ----
243!||====================================================================
244!|| imp_qrf ../engine/source/implicit/imp_pc_inv.F
245!||--- called by ------------------------------------------------------
246!|| imp_saic ../engine/source/implicit/imp_pc_inv.F
247!||--- calls -----------------------------------------------------
248!|| produt_v ../engine/source/implicit/produt_v.F
249!||====================================================================
250 SUBROUTINE imp_qrf(M ,N ,A ,MAXC ,D_R ,
251 . TAU )
252C-----------------------------------------------
253C I m p l i c i t T y p e s
254C-----------------------------------------------
255#include "implicit_f.inc"
256C-----------------------------------------------
257C D u m m y A r g u m e n t s
258C-----------------------------------------------
259 INTEGER MAXC,M,N
260C REAL
261 my_real
262 . A(MAXC,*),D_R(*),TAU(*)
263C-----------------------------------------------
264C L o c a l V a r i a b l e s
265C-----------------------------------------------
266C--A(M,N), a la sortie,Qi=1-uu^t/tau est stocke dans triag inf de A
267C----------R est stoke dans triag sup,+ diag dans D-------
268 INTEGER I,J,K,L
269 my_real
270 . S,NORM2,SCAL
271C------------Qi=1-u(i)u(i)^t/tau(i)--------------------------------
272 DO J=1,n-1
273 scal=zero
274 DO i=j,m
275 scal=max(scal,abs(a(i,j)))
276 ENDDO
277 IF (scal==zero) THEN
278 tau(j)=zero
279 d_r(j)=zero
280 WRITE(*,*)'SIGNULAR A'
281 ELSE
282 scal=one/scal
283 DO i=j,m
284 a(i,j)=a(i,j)*scal
285 ENDDO
286 l=m-j+1
287 CALL produt_v( l ,a(j,j) ,a(j,j) ,norm2)
288 s =sign(sqrt(norm2),a(j,j))
289 a(j,j)=a(j,j)+s
290 tau(j)=s*a(j,j)
291 d_r(j)=-s/scal
292C------------R=QA--------------------------------
293 DO k=j+1,n
294 CALL produt_v( l ,a(j,j) ,a(j,k) ,norm2)
295 s = norm2/tau(j)
296 DO i=j,n
297 a(i,k)=a(i,k)-s*a(i,j)
298 ENDDO
299 ENDDO
300 ENDIF
301 ENDDO
302 d_r(n)=a(n,n)
303C
304 RETURN
305 END
306C-------------b=Q^tb ---------------------------------
307!||====================================================================
308!|| mav_qt ../engine/source/implicit/imp_pc_inv.F
309!||--- called by ------------------------------------------------------
310!|| imp_saic ../engine/source/implicit/imp_pc_inv.F
311!||--- calls -----------------------------------------------------
312!|| produt_v ../engine/source/implicit/produt_v.F
313!||====================================================================
314 SUBROUTINE mav_qt(M ,N ,A ,MAXC ,TAU ,
315 . B )
316C-----------------------------------------------
317C I m p l i c i t T y p e s
318C-----------------------------------------------
319#include "implicit_f.inc"
320C-----------------------------------------------
321C D u m m y A r g u m e n t s
322C-----------------------------------------------
323 INTEGER MAXC,M,N
324C REAL
325 my_real
326 . A(MAXC,*),TAU(*),B(*)
327C-----------------------------------------------
328C L o c a l V a r i a b l e s
329C-----------------------------------------------
330 INTEGER I,J,K
331 my_real
332 . s,norm2
333C--------------------------------------------
334 DO j=1,n-1
335 k=m-j+1
336 CALL produt_v( k ,a(j,j) ,b(j) ,norm2)
337 s = norm2/tau(j)
338 DO i=j,n
339 b(i)=b(i)-s*a(i,j)
340 ENDDO
341 ENDDO
342C
343 RETURN
344 END
345C-------------Rx=b ---------------------------------
346!||====================================================================
347!|| lt_solv ../engine/source/implicit/imp_pc_inv.F
348!||--- called by ------------------------------------------------------
349!|| imp_saic ../engine/source/implicit/imp_pc_inv.F
350!||====================================================================
351 SUBROUTINE lt_solv(N ,A ,MAXC ,D_R ,B )
352C----------R est stoke dans triag_sup de A,+ diag dans D_R-------
353C----------X est sortie dans B-------
354C-----------------------------------------------
355C I m p l i c i t T y p e s
356C-----------------------------------------------
357#include "implicit_f.inc"
358C-----------------------------------------------
359C D u m m y A r g u m e n t s
360C-----------------------------------------------
361 INTEGER MAXC,N
362C REAL
363 my_real
364 . A(MAXC,*),D_R(*),B(*)
365C-----------------------------------------------
366C L o c a l V a r i a b l e s
367C-----------------------------------------------
368 INTEGER I,J,K
369 my_real
370 . s
371C--------------------------------------------
372 b(n)=b(n)/d_r(n)
373 DO i=n-1,1,-1
374 s = zero
375 DO j=i+1,n
376 s=s+a(i,j)*b(j)
377 ENDDO
378 b(i)=(b(i)-s)/d_r(i)
379 ENDDO
380C
381 RETURN
382 END
383!||====================================================================
384!|| get_kijs ../engine/source/implicit/imp_pc_inv.F
385!||--- called by ------------------------------------------------------
386!|| get_suba ../engine/source/implicit/imp_pc_inv.F
387!||--- calls -----------------------------------------------------
388!|| intab0 ../engine/source/implicit/imp_fsa_inv.F
389!||====================================================================
390 SUBROUTINE get_kijs(I ,J ,IADK,JDIK,K_LT ,KIJ)
391C-----------------------------------------------
392C I m p l i c i t T y p e s
393C-----------------------------------------------
394#include "implicit_f.inc"
395C-----------------------------------------------
396C D u m m y A r g u m e n t s
397C-----------------------------------------------
398 INTEGER I,J
399 INTEGER IADK(*) , JDIK(*)
400C REAL
401 my_real
402 . K_LT(*) ,KIJ
403C-----------------------------------------------
404C External function
405C-----------------------------------------------
406 INTEGER INTAB0
407 EXTERNAL intab0
408C-----------------------------------------------
409C L o c a l V a r i a b l e s
410C-----------------------------------------------
411 INTEGER K,K0,NL,N
412C----6---------------------------------------------------------------7---------8
413 K0=iadk(i)
414 nl=iadk(i+1)-k0
415 n=intab0(nl,jdik(k0),j)
416 IF (n>0) THEN
417 kij=k_lt(n+k0-1)
418 ELSE
419 kij=zero
420 ENDIF
421C
422 RETURN
423 END
424
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
imp_pc_inv
subroutine imp_pc_inv(nddl, nnz, iadk, jdik, diag_k, lt_k, iadm, jdim, diag_m, lt_m, psi, nnzm, maxc, maxa, max_l, iopt, nne)
Definition imp_pc_inv.F:37
lt_solv
subroutine lt_solv(n, a, maxc, d_r, b)
Definition imp_pc_inv.F:352
get_kijs
subroutine get_kijs(i, j, iadk, jdik, k_lt, kij)
Definition imp_pc_inv.F:391
mav_qt
subroutine mav_qt(m, n, a, maxc, tau, b)
Definition imp_pc_inv.F:316
get_suba
subroutine get_suba(nddl, iadk, jdik, diag_k, lt_k, nc, jm, a, maxc, maxa, im, m, mj)
Definition imp_pc_inv.F:149
imp_saic
subroutine imp_saic(m, n, a, mj, maxc)
Definition imp_pc_inv.F:109
imp_qrf
subroutine imp_qrf(m, n, a, maxc, d_r, tau)
Definition imp_pc_inv.F:252
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
produt_v
subroutine produt_v(nddl, x, y, r)
Definition produt_v.F:33
arret
subroutine arret(nn)
Definition arret.F:87