OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dsol_lr.F
Go to the documentation of this file.
1C
2C This file is part of MUMPS 5.5.1, released
3C on Tue Jul 12 13:17:24 UTC 2022
4C
5C
6C Copyright 1991-2022 CERFACS, CNRS, ENS Lyon, INP Toulouse, Inria,
7C Mumps Technologies, University of Bordeaux.
8C
9C This version of MUMPS is provided to you free of charge. It is
10C released under the CeCILL-C license
11C (see doc/CeCILL-C_V1-en.txt, doc/CeCILL-C_V1-fr.txt, and
12C https://cecill.info/licences/Licence_CeCILL-C_V1-en.html)
13C
14 MODULE dmumps_sol_lr
15 USE dmumps_lr_type
16 USE dmumps_lr_core
17 USE dmumps_lr_stats
18 USE dmumps_lr_data_m, only: blr_array
19 IMPLICIT NONE
20 CONTAINS
21 SUBROUTINE dmumps_sol_fwd_lr_su
22 & (inode, n, iwhdlr, npiv_global, nslaves,
23 & iw, ipos_init, liw,
24 & liell, wcb, lwcb,
25 & ld_wcbpiv, ld_wcbcb,
26 & ppiv_init, pcb_init,
27 & rhscomp, lrhscomp, nrhs,
28 & posinrhscomp_fwd, jbdeb, jbfin,
29 & mtype, keep, keep8, oocwrite_compatible_with_blr,
30 & iflag, ierror )
31 INTEGER, INTENT(IN) :: INODE, N, IWHDLR, NPIV_GLOBAL, NSLAVES
32 INTEGER, INTENT(IN) :: MTYPE, LIELL, KEEP(500)
33 INTEGER(8), INTENT(IN) :: KEEP8(150)
34 INTEGER, INTENT(IN) :: LIW, IPOS_INIT, LRHSCOMP
35 INTEGER, INTENT(IN) :: IW(LIW), POSINRHSCOMP_FWD(N)
36 INTEGER(8), INTENT(IN) :: LWCB, PPIV_INIT, PCB_INIT
37 INTEGER, INTENT(IN) :: LD_WCBPIV, LD_WCBCB, NRHS, JBDEB, JBFIN
38 DOUBLE PRECISION, INTENT(INOUT) :: WCB(LWCB)
39 INTEGER, INTENT(INOUT) :: IFLAG, IERROR
40 DOUBLE PRECISION, INTENT(INOUT) :: RHSCOMP(LRHSCOMP, NRHS)
41 LOGICAL, INTENT(IN) :: OOCWRITE_COMPATIBLE_WITH_BLR
42 INTEGER :: I, NPARTSASS, NB_BLR , NELIM, LDADIAG,
43 & DIAGSIZ_DYN, DIAGSIZ_STA, IBEG_BLR, IEND_BLR,
44 & LD_CB, NELIM_GLOBAL, NRHS_B, IPOS, KCB
45 INTEGER(8) :: PPIV, PCB
46 INTEGER :: LAST_BLR
47 DOUBLE PRECISION, POINTER, DIMENSION(:) :: DIAG
48 TYPE(lrb_type), POINTER, DIMENSION(:) :: BLR_PANEL
49 DOUBLE PRECISION :: ONE, MONE, ZERO
50 parameter(one = 1.0d0, mone=-1.0d0)
51 parameter(zero=0.0d0)
52 nrhs_b = jbfin-jbdeb+1
53 IF (mtype.EQ.1) THEN
54 IF (associated(blr_array(iwhdlr)%PANELS_L))
55 & THEN
56 npartsass=size(blr_array(iwhdlr)%PANELS_L)
57 nb_blr = size(blr_array(iwhdlr)%BEGS_BLR_STATIC) -1
58 ELSE
59 WRITE(6,*) " Internal error in DMUMPS_SOL_FWD_SU_MASTER"
60 ENDIF
61 ELSE
62 IF (associated(blr_array(iwhdlr)%PANELS_U))
63 & THEN
64 npartsass=size(blr_array(iwhdlr)%PANELS_U)
65 nb_blr = size(blr_array(iwhdlr)%BEGS_BLR_STATIC) -1
66 ENDIF
67 ENDIF
68 IF (nslaves.EQ.0 .OR. (keep(50).eq.0 .and. mtype .NE.1)) THEN
69 last_blr = nb_blr
70 ELSE
71 last_blr = npartsass
72 ENDIF
73 ipos = ipos_init
74 ppiv = ppiv_init
75 nelim_global =
76 & blr_array(iwhdlr)%BEGS_BLR_STATIC(npartsass+1)
77 & - blr_array(iwhdlr)%BEGS_BLR_DYNAMIC(npartsass+1)
78 DO i=1, npartsass
79 ibeg_blr = blr_array(iwhdlr)%BEGS_BLR_DYNAMIC(i)
80 iend_blr = blr_array(iwhdlr)%BEGS_BLR_DYNAMIC(i+1) -1
81 diagsiz_dyn = blr_array(iwhdlr)%BEGS_BLR_DYNAMIC(i+1) -
82 & ibeg_blr
83 diagsiz_sta = blr_array(iwhdlr)%BEGS_BLR_STATIC(i+1) -
84 & ibeg_blr
85 IF (keep(50).NE.0) THEN
86 ldadiag = diagsiz_dyn
87 ELSE
88 ldadiag = diagsiz_sta
89 ENDIF
90 IF (iend_blr.EQ.npiv_global) THEN
91 pcb = pcb_init
92 ELSE
93 pcb = ppiv + int(diagsiz_dyn,8)
94 ENDIF
95 IF ( diagsiz_dyn.EQ.0) cycle
96 nelim = diagsiz_sta - diagsiz_dyn
97 IF ( mtype .EQ. 1 ) THEN
98 blr_panel => blr_array(iwhdlr)%PANELS_L(i)%LRB_PANEL
99 ELSE
100 blr_panel => blr_array(iwhdlr)%PANELS_U(i)%LRB_PANEL
101 END IF
102 diag => blr_array(iwhdlr)%DIAG_BLOCKS(i)%DIAG_BLOCK
103 CALL dmumps_solve_fwd_trsolve (diag(1), int(size(diag),8), 1_8,
104 & diagsiz_dyn , ldadiag, nrhs_b, wcb, lwcb, npiv_global,
105 & ppiv, mtype, keep)
106 IF (nelim.GT.0) THEN
107 kcb = int(pcb-ppiv_init+1)
108 IF (iend_blr.EQ.npiv_global) THEN
109 ld_cb = ld_wcbcb
110 ELSE
111 ld_cb = ld_wcbpiv
112 ENDIF
113 IF (mtype.EQ.1) THEN
114 IF (kcb.LE.npiv_global .AND.
115 & kcb+nelim-1.GT.npiv_global) THEN
116 CALL dgemm('T', 'N', npiv_global-kcb+1, nrhs_b,
117 & diagsiz_dyn, mone,
118 & diag(1+diagsiz_dyn*ldadiag), diagsiz_dyn,
119 & wcb(ppiv), ld_wcbpiv,
120 & one, wcb(pcb), ld_cb)
121 CALL dgemm('T', 'N', kcb+nelim-npiv_global-1,
122 & nrhs_b, diagsiz_dyn, mone,
123 & diag(1+diagsiz_dyn*ldadiag +
124 & (npiv_global-kcb+1)*diagsiz_dyn),
125 & diagsiz_dyn,
126 & wcb(ppiv), ld_wcbpiv,
127 & one, wcb(pcb_init), ld_wcbcb)
128 ELSE
129 CALL dgemm('T', 'N', nelim, nrhs_b, diagsiz_dyn, mone,
130 & diag(1+diagsiz_dyn*ldadiag), diagsiz_dyn,
131 & wcb(ppiv), ld_wcbpiv,
132 & one, wcb(pcb), ld_cb)
133 ENDIF
134 ELSE
135 IF (kcb.LE.npiv_global .AND.
136 & kcb+nelim-1.GT.npiv_global) THEN
137 CALL dgemm('N', 'N', npiv_global-kcb+1,
138 & nrhs_b, diagsiz_dyn, mone,
139 & diag(1+diagsiz_dyn), diagsiz_sta,
140 & wcb(ppiv), ld_wcbpiv,
141 & one, wcb(pcb), ld_cb)
142 CALL dgemm('N', 'N', kcb+nelim-npiv_global-1,
143 & nrhs_b, diagsiz_dyn, mone,
144 & diag(1+diagsiz_dyn+npiv_global-kcb+1),
145 & diagsiz_sta,
146 & wcb(ppiv), ld_wcbpiv,
147 & one, wcb(pcb_init), ld_wcbcb)
148 ELSE
149 CALL dgemm('N', 'n', NELIM, NRHS_B, DIAGSIZ_DYN, MONE,
150 & DIAG(1+DIAGSIZ_DYN), DIAGSIZ_STA,
151 & WCB(PPIV), LD_WCBPIV,
152 & ONE, WCB(PCB), LD_CB)
153 ENDIF
154 ENDIF
155 ENDIF
156 CALL DMUMPS_SOL_FWD_BLR_UPDATE (
157 & WCB, LWCB, 1, LD_WCBPIV, PPIV_INIT, 1,
158 & WCB, LWCB, LD_WCBCB, PCB_INIT,
159 & PPIV,
160 & NRHS_B, NPIV_GLOBAL,
161 & BLR_PANEL, LAST_BLR, I,
162 & BLR_ARRAY(IWHDLR)%BEGS_BLR_STATIC,
163 & KEEP8, KEEP(34), KEEP(450), .FALSE.,
164 & IFLAG, IERROR)
165.LT. IF (IFLAG0) RETURN
166 CALL DMUMPS_SOL_LD_AND_RELOAD_PANEL (
167 & INODE, N, DIAGSIZ_DYN, LIELL, NELIM, NSLAVES,
168 & PPIV,
169 & IW, IPOS, LIW,
170 & DIAG(1), int(size(DIAG),8), 1_8,
171 & WCB, LWCB, LD_WCBPIV,
172 & RHSCOMP, LRHSCOMP, NRHS,
173 & POSINRHSCOMP_FWD, JBDEB, JBFIN,
174 & MTYPE, KEEP, OOCWRITE_COMPATIBLE_WITH_BLR,
175 & .TRUE.
176 & )
177 PPIV = PPIV + int(DIAGSIZ_DYN,8)
178 IPOS = IPOS + DIAGSIZ_DYN
179 ENDDO
180 RETURN
181 END SUBROUTINE DMUMPS_SOL_FWD_LR_SU
182 SUBROUTINE DMUMPS_SOL_SLAVE_LR_U
183 & (INODE, IWHDLR, NPIV_GLOBAL,
184 & WCB, LWCB,
185 & LDX, LDY,
186 & PTRX_INIT, PTRY_INIT,
187 & JBDEB, JBFIN,
188 & MTYPE, KEEP, KEEP8, IFLAG, IERROR )
189 INTEGER, INTENT(IN) :: INODE, IWHDLR, NPIV_GLOBAL
190 INTEGER, INTENT(IN) :: MTYPE, KEEP(500)
191 INTEGER(8), INTENT(IN) :: KEEP8(150)
192 INTEGER(8), INTENT(IN) :: LWCB, PTRX_INIT, PTRY_INIT
193 INTEGER, INTENT(IN) :: LDX, LDY, JBDEB, JBFIN
194 DOUBLE PRECISION, INTENT(INOUT) :: WCB(LWCB)
195 INTEGER, INTENT(INOUT) :: IFLAG, IERROR
196 INTEGER :: I, NPARTSASS, NB_BLR , NRHS_B
197 INTEGER(8) :: PTRX, PTRY
198 TYPE(LRB_TYPE), POINTER, DIMENSION(:) :: BLR_PANEL
199 DOUBLE PRECISION :: ONE, MONE, ZERO
200 PARAMETER (ONE = 1.0D0, MONE=-1.0D0)
201 PARAMETER (ZERO=0.0D0)
202 NRHS_B = JBFIN-JBDEB+1
203 IF (associated(BLR_ARRAY(IWHDLR)%PANELS_L))
204 & THEN
205 NPARTSASS=size(BLR_ARRAY(IWHDLR)%PANELS_L)
206 NB_BLR = size(BLR_ARRAY(IWHDLR)%BEGS_BLR_STATIC)
207 NB_BLR = NB_BLR - 2
208 ELSE
209 WRITE(6,*) " Internal error 1 in DMUMPS_SOL_SLAVE_LR_U"
210 CALL MUMPS_ABORT()
211 ENDIF
212 PTRX = PTRX_INIT
213 PTRY = PTRY_INIT
214 DO I = 1, NPARTSASS
215 BLR_PANEL => BLR_ARRAY(IWHDLR)%PANELS_L(I)%LRB_PANEL
216 IF (associated(BLR_PANEL)) THEN
217.EQ. IF (MTYPE1) THEN
218 CALL DMUMPS_SOL_FWD_BLR_UPDATE (
219 & WCB, LWCB, 1, LDX, -99999_8, 1,
220 & WCB, LWCB, LDY, PTRY,
221 & PTRX,
222 & NRHS_B, NPIV_GLOBAL,
223 & BLR_PANEL, NB_BLR, 0,
224 & BLR_ARRAY(IWHDLR)%BEGS_BLR_STATIC(2:NB_BLR+2),
225 & KEEP8, KEEP(34), KEEP(450), .TRUE., IFLAG, IERROR )
226 ELSE
227 CALL DMUMPS_SOL_BWD_BLR_UPDATE (
228 & WCB, LWCB, 1, LDY, -99999_8, 1,
229 & WCB, LWCB, LDX, PTRX,
230 & PTRY,
231 & NRHS_B, NPIV_GLOBAL,
232 & BLR_PANEL, NB_BLR, 0,
233 & BLR_ARRAY(IWHDLR)%BEGS_BLR_STATIC(2:NB_BLR+2),
234 & KEEP8, KEEP(34), KEEP(450), .TRUE., IFLAG, IERROR )
235 ENDIF
236.EQ. IF (MTYPE 1) THEN
237 PTRX = PTRX + BLR_PANEL(1)%N
238 ELSE
239 PTRY = PTRY + BLR_PANEL(1)%N
240 ENDIF
241.LT. IF (IFLAG0) RETURN
242 ENDIF
243 ENDDO
244 RETURN
245 END SUBROUTINE DMUMPS_SOL_SLAVE_LR_U
246 SUBROUTINE DMUMPS_SOL_FWD_BLR_UPDATE (
247 & ARRAYPIV, LPIV, LPIVCOL, LDPIV, POSPIV, POSPIVCOL,
248 & ARRAYCB, LCB, LDCB, POSCB,
249 & POSDIAG,
250 & NRHS_B, NPIV,
251 & BLR_PANEL, LAST_BLR,
252 & CURRENT_BLR, BEGS_BLR_STATIC,
253 & KEEP8, K34, K450, IS_T2_SLAVE, IFLAG, IERROR )
254!$ USE OMP_LIB
255 INTEGER(8), INTENT(IN) :: LPIV, LCB, POSPIV, POSCB, POSDIAG
256 INTEGER, INTENT(IN) :: LPIVCOL, POSPIVCOL
257 DOUBLE PRECISION, INTENT(INOUT) :: ARRAYPIV(LPIV,LPIVCOL)
258 DOUBLE PRECISION, INTENT(INOUT) :: ARRAYCB(LCB)
259 INTEGER, INTENT(IN) :: LAST_BLR, NRHS_B, LDPIV, LDCB,
260 & CURRENT_BLR, NPIV, K34, K450
261 TYPE(LRB_TYPE), TARGET,INTENT(IN) ::
262 & BLR_PANEL(:)
263 INTEGER :: BEGS_BLR_STATIC(:)
264 LOGICAL, INTENT(IN) :: IS_T2_SLAVE
265 INTEGER(8), INTENT(IN) :: KEEP8(150)
266 INTEGER, INTENT(INOUT) :: IFLAG, IERROR
267 INTEGER :: I, K, M, N, IBEG_BLOCK, IEND_BLOCK
268 INTEGER :: KMAX
269 INTEGER(8) :: POSBLOCK
270 INTEGER :: allocok
271 TYPE(LRB_TYPE), POINTER :: LRB
272 DOUBLE PRECISION, ALLOCATABLE,DIMENSION(:) :: TEMP_BLOCK
273 DOUBLE PRECISION :: ONE, MONE, ZERO
274 PARAMETER (ONE = 1.0D0, MONE=-1.0D0)
275 PARAMETER (ZERO=0.0D0)
276 KMAX = -1
277 DO I = CURRENT_BLR+1, LAST_BLR
278 KMAX = max(KMAX, BLR_PANEL(I-CURRENT_BLR)%K)
279 ENDDO
280#if defined(BLR_MT)
281!$OMP PARALLEL PRIVATE(TEMP_BLOCK, allocok)
282#endif
283.GT. IF (KMAX0) THEN
284 allocate(TEMP_BLOCK(KMAX*NRHS_B), stat=allocok )
285.GT. IF (allocok 0) THEN
286 IFLAG = -13
287 IERROR = NRHS_B * KMAX
288 write(*,*) 'allocation problem in blr routine
289 & dmumps_sol_fwd_blr_update: ',
290 & 'not enough memory? memory requested = ', IERROR
291 ENDIF
292 ENDIF
293#if defined(BLR_MT)
294!$OMP BARRIER
295#endif
296#if defined(BLR_MT)
297!$OMP DO SCHEDULE(DYNAMIC,1)
298!$OMP& PRIVATE(IBEG_BLOCK, IEND_BLOCK, LRB, K, M, N,
299!$OMP& POSBLOCK)
300#endif
301 DO I = CURRENT_BLR+1, LAST_BLR
302.LT. IF (IFLAG0) CYCLE
303 IBEG_BLOCK = BEGS_BLR_STATIC(I)
304 IEND_BLOCK = BEGS_BLR_STATIC(I+1)-1
305.EQ. IF (IBEG_BLOCK IEND_BLOCK + 1) CYCLE
306 LRB => BLR_PANEL(I-CURRENT_BLR)
307 K = LRB%K
308 M = LRB%M
309 N = LRB%N
310 IF (LRB%ISLR) THEN
311.GT. IF (K0) THEN
312 CALL dgemm('n', 'n', K, NRHS_B, N, ONE,
313 & LRB%R(1,1), K, ARRAYPIV(POSDIAG,POSPIVCOL), LDPIV,
314 & ZERO, TEMP_BLOCK(1), K)
315 IF (IS_T2_SLAVE) THEN
316 POSBLOCK = POSCB+int(IBEG_BLOCK-1,8)
317 CALL dgemm('n', 'n', M, NRHS_B, K, MONE,
318 & LRB%Q(1,1), M, TEMP_BLOCK(1), K,
319 & ONE, ARRAYCB(POSBLOCK), LDCB)
320.LE..AND..GT. ELSEIF (IBEG_BLOCKNPIVIEND_BLOCKNPIV) THEN
321 POSBLOCK = POSPIV+int(IBEG_BLOCK-1,8)
322 CALL dgemm('n', 'n', NPIV-IBEG_BLOCK+1, NRHS_B, K,
323 & MONE, LRB%Q(1,1), M, TEMP_BLOCK(1), K,
324 & ONE, ARRAYPIV(POSBLOCK,POSPIVCOL), LDPIV)
325 CALL dgemm('n', 'n', IBEG_BLOCK+M-NPIV-1, NRHS_B, K,
326 & MONE, LRB%Q(NPIV-IBEG_BLOCK+2,1), M, TEMP_BLOCK(1),
327 & K, ONE, ARRAYCB(POSCB), LDCB)
328.LE. ELSEIF (IBEG_BLOCKNPIV) THEN
329 POSBLOCK = POSPIV+int(IBEG_BLOCK-1,8)
330 CALL dgemm('n', 'n', M, NRHS_B, K, MONE,
331 & LRB%Q(1,1), M, TEMP_BLOCK(1), K,
332 & ONE, ARRAYPIV(POSBLOCK,POSPIVCOL), LDPIV)
333 ELSE
334 POSBLOCK = POSCB+int(IBEG_BLOCK-1-NPIV,8)
335 CALL dgemm('n', 'n', M, NRHS_B, K, MONE,
336 & LRB%Q(1,1), M, TEMP_BLOCK(1), K,
337 & ONE, ARRAYCB(POSBLOCK), LDCB)
338 ENDIF
339 ENDIF
340 ELSE
341 IF (IS_T2_SLAVE) THEN
342 POSBLOCK = POSCB + int(IBEG_BLOCK-1,8)
343 CALL dgemm('n', 'n', M, NRHS_B, N, MONE,
344 & LRB%Q(1,1), M, ARRAYPIV(POSDIAG,POSPIVCOL), LDPIV,
345 & ONE, ARRAYCB(POSBLOCK), LDCB)
346.LE..AND..GT. ELSEIF (IBEG_BLOCKNPIVIEND_BLOCKNPIV) THEN
347 POSBLOCK = POSPIV+int(IBEG_BLOCK-1,8)
348 CALL dgemm('n', 'n', NPIV-IBEG_BLOCK+1, NRHS_B, N, MONE,
349 & LRB%Q(1,1), M, ARRAYPIV(POSDIAG,POSPIVCOL), LDPIV,
350 & ONE, ARRAYPIV(POSBLOCK,POSPIVCOL), LDPIV)
351 CALL dgemm('n', 'n', IBEG_BLOCK+M-NPIV-1, NRHS_B, N, MONE,
352 & LRB%Q(NPIV-IBEG_BLOCK+2,1), M,
353 & ARRAYPIV(POSDIAG,POSPIVCOL),
354 & LDPIV, ONE, ARRAYCB(POSCB), LDCB)
355.LE. ELSEIF (IBEG_BLOCKNPIV) THEN
356 POSBLOCK = POSPIV+int(IBEG_BLOCK-1,8)
357 CALL dgemm('n', 'n', M, NRHS_B, N, MONE,
358 & LRB%Q(1,1), M, ARRAYPIV(POSDIAG,POSPIVCOL), LDPIV,
359 & ONE, ARRAYPIV(POSBLOCK,POSPIVCOL), LDPIV)
360 ELSE
361 POSBLOCK = POSCB + int(IBEG_BLOCK-1-NPIV,8)
362 CALL dgemm('n', 'n', M, NRHS_B, N, MONE,
363 & LRB%Q(1,1), M, ARRAYPIV(POSDIAG,POSPIVCOL), LDPIV,
364 & ONE, ARRAYCB(POSBLOCK), LDCB)
365 ENDIF
366 ENDIF
367 ENDDO
368#if defined(BLR_MT)
369!$OMP END DO
370#endif
371.GT. IF (KMAX0) THEN
372 IF (allocated(TEMP_BLOCK)) deallocate(TEMP_BLOCK)
373 ENDIF
374#if defined(BLR_MT)
375!$OMP END PARALLEL
376#endif
377 RETURN
378 END SUBROUTINE DMUMPS_SOL_FWD_BLR_UPDATE
379 SUBROUTINE DMUMPS_SOL_BWD_LR_SU
380 & ( INODE, IWHDLR, NPIV_GLOBAL, NSLAVES,
381 & LIELL, WCB, LWCB, NRHS_B, PTWCB,
382 & RHSCOMP, LRHSCOMP, NRHS,
383 & IPOSINRHSCOMP, JBDEB,
384 & MTYPE, KEEP, KEEP8,
385 & IFLAG, IERROR )
386 INTEGER, INTENT(IN) :: INODE, IWHDLR, NPIV_GLOBAL, NSLAVES
387 INTEGER, INTENT(IN) :: MTYPE, LIELL, KEEP(500)
388 INTEGER(8), INTENT(IN) :: KEEP8(150)
389 INTEGER, INTENT(IN) :: IPOSINRHSCOMP, JBDEB, LRHSCOMP, NRHS
390 INTEGER(8), INTENT(IN) :: LWCB, PTWCB
391 INTEGER, INTENT(IN) :: NRHS_B
392 INTEGER, INTENT(INOUT) :: IFLAG, IERROR
393 DOUBLE PRECISION, INTENT(INOUT) :: WCB(LWCB)
394 DOUBLE PRECISION RHSCOMP(LRHSCOMP,NRHS)
395 INTEGER :: I, NPARTSASS, NB_BLR, LAST_BLR,
396 & NELIM_PANEL, LD_WCB,
397 & DIAGSIZ_DYN, DIAGSIZ_STA, LDADIAG,
398 & IEND_BLR, IBEG_BLR, PCBINRHSCOMP
399 INTEGER(8) :: PCB_LAST, PWCB
400 INTEGER :: IPIV_PANEL
401 DOUBLE PRECISION, POINTER, DIMENSION(:) :: DIAG
402 TYPE(LRB_TYPE), POINTER, DIMENSION(:) :: BLR_PANEL
403 DOUBLE PRECISION :: ONE, MONE, ZERO
404 PARAMETER (ONE = 1.0D0, MONE=-1.0D0)
405 PARAMETER (ZERO=0.0D0)
406.EQ..AND..EQ. IF ((MTYPE1)(KEEP(50)0)) THEN
407 IF (associated(BLR_ARRAY(IWHDLR)%PANELS_U))
408 & THEN
409 NPARTSASS=size(BLR_ARRAY(IWHDLR)%PANELS_U)
410 NB_BLR = size(BLR_ARRAY(IWHDLR)%BEGS_BLR_STATIC) -1
411 ENDIF
412 ELSE
413 IF (associated(BLR_ARRAY(IWHDLR)%PANELS_L))
414 & THEN
415 NPARTSASS=size(BLR_ARRAY(IWHDLR)%PANELS_L)
416 NB_BLR = size(BLR_ARRAY(IWHDLR)%BEGS_BLR_STATIC) -1
417 ELSE
418 WRITE(6,*) " Internal error in DMUMPS_SOL_FWD_SU_MASTER"
419 ENDIF
420 ENDIF
421 PCBINRHSCOMP= IPOSINRHSCOMP + NPIV_GLOBAL
422 PCB_LAST = PTWCB + int(LIELL ,8)
423 PWCB = PTWCB + int(NPIV_GLOBAL,8)
424 LD_WCB = LIELL
425 DO I=NPARTSASS,1,-1
426 IBEG_BLR = BLR_ARRAY(IWHDLR)%BEGS_BLR_DYNAMIC(I)
427 IEND_BLR = BLR_ARRAY(IWHDLR)%BEGS_BLR_DYNAMIC(I+1) -1
428 DIAGSIZ_DYN = BLR_ARRAY(IWHDLR)%BEGS_BLR_DYNAMIC(I+1) -
429 & IBEG_BLR
430 DIAGSIZ_STA = BLR_ARRAY(IWHDLR)%BEGS_BLR_STATIC(I+1) -
431 & IBEG_BLR
432.NE. IF (KEEP(50)0) THEN
433 LDADIAG = DIAGSIZ_DYN
434 ELSE
435 LDADIAG = DIAGSIZ_STA
436 ENDIF
437.EQ. IF (DIAGSIZ_DYN0) CYCLE
438 NELIM_PANEL = DIAGSIZ_STA - DIAGSIZ_DYN
439 IPIV_PANEL = IPOSINRHSCOMP + IBEG_BLR -1
440.EQ..AND..EQ. IF ( MTYPE 1 KEEP(50)0) THEN
441 BLR_PANEL => BLR_ARRAY(IWHDLR)%PANELS_U(I)%LRB_PANEL
442 ELSE
443 BLR_PANEL => BLR_ARRAY(IWHDLR)%PANELS_L(I)%LRB_PANEL
444 END IF
445.EQ..AND..GT..AND..NE. IF (KEEP(50)0 NSLAVES0 MTYPE1) THEN
446 LAST_BLR = NPARTSASS
447 ELSE
448 LAST_BLR = NB_BLR
449 ENDIF
450 CALL DMUMPS_SOL_BWD_BLR_UPDATE (
451 & RHSCOMP, int(LRHSCOMP,8), NRHS, LRHSCOMP,
452 & int(IPOSINRHSCOMP,8), JBDEB,
453 & WCB, LWCB, LD_WCB, PWCB,
454 & int(IPIV_PANEL,8),
455 & NRHS_B, NPIV_GLOBAL,
456 & BLR_PANEL, LAST_BLR,
457 & I, BLR_ARRAY(IWHDLR)%BEGS_BLR_STATIC,
458 & KEEP8, KEEP(34), KEEP(450), .FALSE., IFLAG, IERROR)
459.LT. IF (IFLAG0) RETURN
460 DIAG => BLR_ARRAY(IWHDLR)%DIAG_BLOCKS(I)%DIAG_BLOCK
461.GT. IF (NELIM_PANEL0) THEN
462.EQ..AND..EQ. IF (MTYPE1KEEP(50)0) THEN
463.EQ. IF (IEND_BLRNPIV_GLOBAL) THEN
464 CALL dgemm('t', 'n', DIAGSIZ_DYN, NRHS_B, NELIM_PANEL,
465 & MONE, DIAG(1+DIAGSIZ_DYN), DIAGSIZ_STA, WCB(PWCB),
466 & LD_WCB, ONE , RHSCOMP(IPIV_PANEL,JBDEB),LRHSCOMP)
467 ELSE
468.LE..AND. IF (IEND_BLR+1NPIV_GLOBAL
469.GT. & IEND_BLR+NELIM_PANELNPIV_GLOBAL) THEN
470 CALL dgemm('t', 'n', diagsiz_dyn, nrhs_b,
471 & npiv_global-iend_blr,
472 & mone, diag(1+diagsiz_dyn), diagsiz_sta,
473 & rhscomp(ipiv_panel+diagsiz_dyn,jbdeb), lrhscomp,
474 & one, rhscomp(ipiv_panel,jbdeb), lrhscomp)
475 CALL dgemm('T', 'N', diagsiz_dyn, nrhs_b,
476 & iend_blr+nelim_panel-npiv_global,
477 & mone, diag(1+diagsiz_dyn+npiv_global-iend_blr),
478 & diagsiz_sta,
479 & wcb(pwcb), ld_wcb,
480 & one, rhscomp(ipiv_panel,jbdeb), lrhscomp)
481 ELSE
482 CALL dgemm('T', 'N', diagsiz_dyn, nrhs_b, nelim_panel,
483 & mone, diag(1+diagsiz_dyn), diagsiz_sta,
484 & rhscomp(ipiv_panel+diagsiz_dyn,jbdeb), lrhscomp,
485 & one, rhscomp(ipiv_panel,jbdeb), lrhscomp)
486 ENDIF
487 ENDIF
488 ELSE
489 IF (iend_blr.EQ.npiv_global) THEN
490 CALL dgemm('N', 'N', diagsiz_dyn, nrhs_b, nelim_panel,
491 & mone, diag(1+diagsiz_dyn*ldadiag), diagsiz_dyn,
492 & wcb(pwcb), ld_wcb, one,
493 & rhscomp(ipiv_panel,jbdeb), lrhscomp)
494 ELSE
495 IF (iend_blr+1.LE.npiv_global .AND.
496 & iend_blr+nelim_panel.GT.npiv_global) THEN
497 CALL dgemm('N', 'n', DIAGSIZ_DYN, NRHS_B,
498 & NPIV_GLOBAL-IEND_BLR,
499 & MONE, DIAG(1+DIAGSIZ_DYN*LDADIAG), DIAGSIZ_DYN,
500 & RHSCOMP(IPIV_PANEL+DIAGSIZ_DYN,JBDEB), LRHSCOMP,
501 & ONE, RHSCOMP(IPIV_PANEL,JBDEB), LRHSCOMP)
502 CALL dgemm('n', 'n', DIAGSIZ_DYN, NRHS_B,
503 & IEND_BLR+NELIM_PANEL-NPIV_GLOBAL,
504 & MONE, DIAG(1+DIAGSIZ_DYN*LDADIAG +
505 & (NPIV_GLOBAL-IEND_BLR)*DIAGSIZ_DYN),
506 & DIAGSIZ_DYN,
507 & WCB(PWCB), LD_WCB,
508 & ONE, RHSCOMP(IPIV_PANEL,JBDEB), LRHSCOMP)
509 ELSE
510 CALL dgemm('n', 'n', DIAGSIZ_DYN, NRHS_B, NELIM_PANEL,
511 & MONE, DIAG(1+DIAGSIZ_DYN*LDADIAG), DIAGSIZ_DYN,
512 & RHSCOMP(IPIV_PANEL+DIAGSIZ_DYN,JBDEB), LRHSCOMP,
513 & ONE, RHSCOMP(IPIV_PANEL,JBDEB), LRHSCOMP)
514 ENDIF
515 ENDIF
516 ENDIF
517 ENDIF
518.LT. IF (IFLAG0) RETURN
519 CALL DMUMPS_SOLVE_BWD_LR_TRSOLVE (
520 & DIAG(1), size(DIAG), DIAGSIZ_DYN, NELIM_PANEL, LIELL,
521 & NRHS_B, WCB, LWCB,
522 & RHSCOMP, LRHSCOMP, NRHS,
523 & IPIV_PANEL, JBDEB,
524 & MTYPE, KEEP )
525 ENDDO
526 RETURN
527 END SUBROUTINE DMUMPS_SOL_BWD_LR_SU
528 SUBROUTINE DMUMPS_SOL_BWD_BLR_UPDATE (
529 & ARRAYPIV, LPIV, LPIVCOL, LDPIV, POSPIV, POSPIVCOL,
530 & ARRAYCB, LCB, LDCB, POSCB,
531 & POSDIAG,
532 & NRHS_B, NPIV,
533 & BLR_PANEL, LAST_BLR, CURRENT_BLR,
534 & BEGS_BLR_STATIC,
535 & KEEP8, K34, K450, IS_T2_SLAVE,
536 & IFLAG, IERROR)
537!$ USE OMP_LIB
538 INTEGER(8), INTENT(IN) :: LPIV, LCB, POSPIV, POSCB, POSDIAG
539 INTEGER,INTENT(IN) :: LPIVCOL, POSPIVCOL
540 DOUBLE PRECISION, INTENT(INOUT) :: ARRAYPIV(LPIV,LPIVCOL)
541 DOUBLE PRECISION, INTENT(INOUT) :: ARRAYCB(LCB)
542 INTEGER, INTENT(IN) :: LAST_BLR, NRHS_B, LDPIV, LDCB,
543 & CURRENT_BLR, NPIV, K34, K450
544 TYPE(LRB_TYPE), TARGET,INTENT(IN) ::
545 & BLR_PANEL(:)
546 INTEGER(8), INTENT(IN) :: KEEP8(150)
547 LOGICAL, INTENT(IN) :: IS_T2_SLAVE
548 INTEGER :: BEGS_BLR_STATIC(:)
549 INTEGER, INTENT(INOUT) :: IFLAG, IERROR
550 INTEGER :: I, K, M, N, IBEG_BLOCK, IEND_BLOCK
551 INTEGER :: KMAX
552 INTEGER(8) :: POSBLOCK
553 TYPE(LRB_TYPE), POINTER :: LRB
554 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: TEMP_BLOCK
555 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: DEST_ARRAY
556 INTEGER :: allocok
557 DOUBLE PRECISION :: ONE, MONE, ZERO
558 PARAMETER (ONE = 1.0D0, MONE=-1.0D0)
559 PARAMETER (ZERO=0.0D0)
560 KMAX = -1
561 DO I = CURRENT_BLR+1, LAST_BLR
562 KMAX = max(KMAX, BLR_PANEL(I-CURRENT_BLR)%K)
563 ENDDO
564.LT. IF (CURRENT_BLRLAST_BLR) THEN
565 N = BLR_PANEL(1)%N
566 ELSE
567 RETURN
568 ENDIF
569 allocate(DEST_ARRAY(N*NRHS_B),stat=allocok)
570.GT. IF (allocok 0) THEN
571 IFLAG = -13
572 IERROR = N * NRHS_B
573 GOTO 100
574 ENDIF
575 DEST_ARRAY = ZERO
576#if defined(BLR_MT)
577!$OMP PARALLEL PRIVATE(TEMP_BLOCK,allocok)
578#endif
579.GT. IF (KMAX0) THEN
580 allocate(TEMP_BLOCK(KMAX*NRHS_B), stat=allocok )
581.GT. IF (allocok 0) THEN
582 IFLAG = -13
583 IERROR = NRHS_B * KMAX
584 write(*,*) 'allocation problem in blr routine
585 & dmumps_sol_bwd_blr_update: ',
586 & 'not enough memory? memory requested = ', ierror
587 ENDIF
588 ENDIF
589#if defined(BLR_MT)
590!$OMP BARRIER
591#endif
592#if defined(BLR_MT)
593!$OMP DO SCHEDULE(DYNAMIC,1)
594!$OMP& PRIVATE(IBEG_BLOCK, IEND_BLOCK, LRB, K, M,
595!$OMP& POSBLOCK)
596!$OMP& REDUCTION(+:DEST_ARRAY)
597#endif
598 DO i = current_blr+1, last_blr
599 IF (iflag.LT.0) cycle
600 ibeg_block = begs_blr_static(i)
601 iend_block = begs_blr_static(i+1)-1
602 lrb => blr_panel(i-current_blr)
603 k = lrb%K
604 m = lrb%M
605 IF (lrb%ISLR) THEN
606 IF (k.GT.0) THEN
607 IF (is_t2_slave) THEN
608 posblock = poscb +int(ibeg_block-1,8)
609 CALL dgemm('T', 'N', k, nrhs_b, m, one,
610 & lrb%Q(1,1), m,
611 & arraycb(posblock), ldcb, zero,
612 & temp_block(1), k)
613 ELSE IF (ibeg_block.LE.npiv.AND.iend_block.GT.npiv) THEN
614 posblock = pospiv+int(ibeg_block-1,8)
615 CALL dgemm('T', 'N', k, nrhs_b, npiv-ibeg_block+1, one,
616 & lrb%Q(1,1), m,
617 & arraypiv(posblock,pospivcol), ldpiv,
618 & zero, temp_block(1), k)
619 CALL dgemm('T', 'N', k, nrhs_b, ibeg_block+m-npiv-1,
620 & one, lrb%Q(npiv-ibeg_block+2,1), m,
621 & arraycb(poscb), ldcb,
622 & one,
623 & temp_block(1), k)
624 ELSEIF (ibeg_block.LE.npiv) THEN
625 posblock = pospiv+int(ibeg_block-1,8)
626 CALL dgemm('T', 'N', k, nrhs_b, m, one,
627 & lrb%Q(1,1), m,
628 & arraypiv(posblock,pospivcol), ldpiv,
629 & zero, temp_block(1), k)
630 ELSE
631 posblock = poscb+int(ibeg_block-1-npiv,8)
632 CALL dgemm('T', 'N', k, nrhs_b, m, one,
633 & lrb%Q(1,1), m,
634 & arraycb(posblock), ldcb, zero,
635 & temp_block(1), k)
636 ENDIF
637 CALL dgemm('T', 'N', n, nrhs_b, k, mone,
638 & lrb%R(1,1), k,
639 & temp_block(1), k, one,
640 & dest_array(1), n)
641 ENDIF
642 ELSE
643 IF (is_t2_slave) THEN
644 posblock = poscb+int(ibeg_block-1,8)
645 CALL dgemm('T', 'N', n, nrhs_b, m, mone,
646 & lrb%Q(1,1), m, arraycb(posblock), ldcb,
647 & one, dest_array(1), n)
648 ELSE IF (ibeg_block.LE.npiv.AND.iend_block.GT.npiv) THEN
649 posblock = pospiv+int(ibeg_block-1,8)
650 CALL dgemm('T', 'N', n, nrhs_b, npiv-ibeg_block+1, mone,
651 & lrb%Q(1,1), m, arraypiv(posblock,pospivcol), ldpiv,
652 & one, dest_array(1), n)
653 CALL dgemm('T', 'N', n, nrhs_b, ibeg_block+m-npiv-1, mone,
654 & lrb%Q(npiv-ibeg_block+2,1), m, arraycb(poscb),
655 & ldcb, one, dest_array(1), n)
656 ELSEIF (ibeg_block.LE.npiv) THEN
657 posblock = pospiv+int(ibeg_block-1,8)
658 CALL dgemm('T', 'N', n, nrhs_b, m, mone,
659 & lrb%Q(1,1), m, arraypiv(posblock,pospivcol), ldpiv,
660 & one, dest_array(1), n)
661 ELSE
662 posblock = poscb+int(ibeg_block-1-npiv,8)
663 CALL dgemm('T', 'N', n, nrhs_b, m, mone,
664 & lrb%Q(1,1), m, arraycb(posblock), ldcb,
665 & one, dest_array(1), n)
666 ENDIF
667 ENDIF
668 ENDDO
669#if defined(BLR_MT)
670!$OMP END DO
671#endif
672 IF (kmax.GT.0) THEN
673 IF (allocated(temp_block)) deallocate(temp_block)
674 ENDIF
675#if defined(BLR_MT)
676!$OMP END PARALLEL
677#endif
678 IF (is_t2_slave) THEN
679 DO i=1,nrhs_b
680 call daxpy(n, one, dest_array((i-1)*n+1), 1,
681 & arraypiv(posdiag+(i-1)*ldpiv,pospivcol), 1)
682 ENDDO
683 ELSE
684 DO i=1,nrhs_b
685 call daxpy(n, one, dest_array((i-1)*n+1), 1,
686 & arraypiv(posdiag,pospivcol+i-1), 1)
687 ENDDO
688 ENDIF
689 100 CONTINUE
690 IF (allocated(dest_array)) DEALLOCATE(dest_array)
691 RETURN
692 END SUBROUTINE dmumps_sol_bwd_blr_update
693 END MODULE dmumps_sol_lr
694 SUBROUTINE dmumps_solve_bwd_lr_trsolve (
695 & DIAG, LDIAG, NPIV, NELIM, LIELL,
696 & NRHS_B, W, LWC,
697 & RHSCOMP, LRHSCOMP, NRHS,
698 & PPIVINRHSCOMP, JBDEB,
699 & MTYPE, KEEP)
700 INTEGER, INTENT(IN) :: MTYPE, LIELL, NPIV, NELIM, KEEP(500)
701 INTEGER, INTENT(IN) :: NRHS_B, LDIAG
702 INTEGER, INTENT(IN) :: PPIVINRHSCOMP, JBDEB, LRHSCOMP, NRHS
703 INTEGER(8), INTENT(IN) :: LWC
704 DOUBLE PRECISION, INTENT(IN) :: DIAG(LDIAG)
705 DOUBLE PRECISION, INTENT(INOUT) :: W(LWC)
706 DOUBLE PRECISION RHSCOMP(LRHSCOMP,NRHS)
707 INTEGER :: LDAJ
708 DOUBLE PRECISION ONE
709 parameter(one = 1.0d0)
710 IF ( mtype .eq. 1 ) THEN
711 ldaj = npiv + nelim
712 CALL dtrsm('L','L','T','N', npiv, nrhs_b, one, diag(1),
713 & ldaj, rhscomp(ppivinrhscomp,jbdeb),
714 & lrhscomp)
715 ELSE
716 IF ( keep(50) .EQ. 0 ) THEN
717 ldaj=npiv+nelim
718 ELSE
719 ldaj=npiv
720 ENDIF
721 CALL dtrsm('L','U','N','U', npiv, nrhs_b, one, diag(1),
722 & ldaj, rhscomp(ppivinrhscomp,jbdeb), lrhscomp)
723 END IF
724 RETURN
725 END SUBROUTINE dmumps_solve_bwd_lr_trsolve
dmumps_solve_fwd_trsolve
subroutine dmumps_solve_fwd_trsolve(a, la, apos, npiv, ldadiag, nrhs_b, wcb, lwcb, lda_wcb, ppiv_courant, mtype, keep)
Definition dsol_aux.F:1147
dmumps_solve_bwd_lr_trsolve
subroutine dmumps_solve_bwd_lr_trsolve(diag, ldiag, npiv, nelim, liell, nrhs_b, w, lwc, rhscomp, lrhscomp, nrhs, ppivinrhscomp, jbdeb, mtype, keep)
Definition dsol_lr.F:700
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187
dtrsm
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181
dmumps_lr_core
Definition dlr_core.F:18
dmumps_lr_data_m
Definition dmumps_lr_data_m.F:14
dmumps_lr_data_m::blr_array
type(blr_struc_t), dimension(:), pointer, save, public blr_array
Definition dmumps_lr_data_m.F:61
dmumps_lr_stats
Definition dlr_stats.F:14
dmumps_lr_type
Definition dlr_type.F:14
dmumps_sol_lr
Definition dsol_lr.F:14
dmumps_sol_lr::dmumps_sol_fwd_blr_update
subroutine dmumps_sol_fwd_blr_update(arraypiv, lpiv, lpivcol, ldpiv, pospiv, pospivcol, arraycb, lcb, ldcb, poscb, posdiag, nrhs_b, npiv, blr_panel, last_blr, current_blr, begs_blr_static, keep8, k34, k450, is_t2_slave, iflag, ierror)
Definition dsol_lr.F:254
dmumps_sol_lr::dmumps_sol_fwd_lr_su
subroutine dmumps_sol_fwd_lr_su(inode, n, iwhdlr, npiv_global, nslaves, iw, ipos_init, liw, liell, wcb, lwcb, ld_wcbpiv, ld_wcbcb, ppiv_init, pcb_init, rhscomp, lrhscomp, nrhs, posinrhscomp_fwd, jbdeb, jbfin, mtype, keep, keep8, oocwrite_compatible_with_blr, iflag, ierror)
Definition dsol_lr.F:31
dmumps_sol_lr::dmumps_sol_bwd_blr_update
subroutine dmumps_sol_bwd_blr_update(arraypiv, lpiv, lpivcol, ldpiv, pospiv, pospivcol, arraycb, lcb, ldcb, poscb, posdiag, nrhs_b, npiv, blr_panel, last_blr, current_blr, begs_blr_static, keep8, k34, k450, is_t2_slave, iflag, ierror)
Definition dsol_lr.F:537
dmumps_lr_type::lrb_type
Definition dlr_type.F:16