OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
csol_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 cmumps_sol_lr
15 USE cmumps_lr_type
16 USE cmumps_lr_core
17 USE cmumps_lr_stats
18 USE cmumps_lr_data_m, only: blr_array
19 IMPLICIT NONE
20 CONTAINS
21 SUBROUTINE cmumps_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 COMPLEX, INTENT(INOUT) :: WCB(LWCB)
39 INTEGER, INTENT(INOUT) :: IFLAG, IERROR
40 COMPLEX, 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 COMPLEX, POINTER, DIMENSION(:) :: DIAG
48 TYPE(lrb_type), POINTER, DIMENSION(:) :: BLR_PANEL
49 COMPLEX :: ONE, MONE, ZERO
50 parameter(one=(1.0e0,0.0e0), mone=(-1.0e0,0.0e0))
51 parameter(zero=(0.0e0,0.0e0))
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 CMUMPS_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 cmumps_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 cgemm('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 cgemm('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 cgemm('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 cgemm('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 cgemm('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 cgemm('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 CMUMPS_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 CMUMPS_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 CMUMPS_SOL_FWD_LR_SU
182 SUBROUTINE CMUMPS_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 COMPLEX, 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 COMPLEX :: ONE, MONE, ZERO
200 PARAMETER (ONE=(1.0E0,0.0E0), MONE=(-1.0E0,0.0E0))
201 PARAMETER (ZERO=(0.0E0,0.0E0))
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 CMUMPS_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 CMUMPS_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 CMUMPS_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 CMUMPS_SOL_SLAVE_LR_U
246 SUBROUTINE CMUMPS_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 COMPLEX, INTENT(INOUT) :: ARRAYPIV(LPIV,LPIVCOL)
258 COMPLEX, 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 COMPLEX, ALLOCATABLE,DIMENSION(:) :: TEMP_BLOCK
273 COMPLEX :: ONE, MONE, ZERO
274 PARAMETER (ONE=(1.0E0,0.0E0), MONE=(-1.0E0,0.0E0))
275 PARAMETER (ZERO=(0.0E0,0.0E0))
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 & cmumps_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 cgemm('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 cgemm('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 cgemm('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 cgemm('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 cgemm('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 cgemm('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 cgemm('n', 'n', m, nrhs_b, n, mone,
344 & lrb%Q(1,1), m, arraypiv(posdiag,pospivcol), ldpiv,
345 & one, arraycb(posblock), ldcb)
346 ELSEIF (ibeg_block.LE.npiv.AND.iend_block.GT.npiv) THEN
347 posblock = pospiv+int(ibeg_block-1,8)
348 CALL cgemm('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 cgemm('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 ELSEIF (ibeg_block.LE.npiv) THEN
356 posblock = pospiv+int(ibeg_block-1,8)
357 CALL cgemm('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 cgemm('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 IF (kmax.GT.0) 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 cmumps_sol_fwd_blr_update
379 SUBROUTINE cmumps_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 COMPLEX, INTENT(INOUT) :: WCB(LWCB)
394 COMPLEX 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 COMPLEX, POINTER, DIMENSION(:) :: DIAG
402 TYPE(lrb_type), POINTER, DIMENSION(:) :: BLR_PANEL
403 COMPLEX :: ONE, MONE, ZERO
404 parameter(one=(1.0e0,0.0e0), mone=(-1.0e0,0.0e0))
405 parameter(zero=(0.0e0,0.0e0))
406 IF ((mtype.EQ.1).AND.(keep(50).EQ.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 CMUMPS_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 IF (keep(50).NE.0) THEN
433 ldadiag = diagsiz_dyn
434 ELSE
435 ldadiag = diagsiz_sta
436 ENDIF
437 IF (diagsiz_dyn.EQ.0) cycle
438 nelim_panel = diagsiz_sta - diagsiz_dyn
439 ipiv_panel = iposinrhscomp + ibeg_blr -1
440 IF ( mtype .EQ. 1 .AND. keep(50).EQ.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 IF (keep(50).EQ.0 .AND. nslaves.GT.0 .AND. mtype.NE.1) THEN
446 last_blr = npartsass
447 ELSE
448 last_blr = nb_blr
449 ENDIF
450 CALL cmumps_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 IF (iflag.LT.0) RETURN
460 diag => blr_array(iwhdlr)%DIAG_BLOCKS(i)%DIAG_BLOCK
461 IF (nelim_panel.GT.0) THEN
462 IF (mtype.EQ.1.AND.keep(50).EQ.0) THEN
463 IF (iend_blr.EQ.npiv_global) THEN
464 CALL cgemm('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 IF (iend_blr+1.LE.npiv_global .AND.
469 & iend_blr+nelim_panel.GT.npiv_global) THEN
470 CALL cgemm('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 cgemm('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 cgemm('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 cgemm('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.LE..AND. IF (IEND_BLR+1NPIV_GLOBAL
496.GT. & IEND_BLR+NELIM_PANELNPIV_GLOBAL) THEN
497 CALL cgemm('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 cgemm('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 cgemm('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 CMUMPS_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 CMUMPS_SOL_BWD_LR_SU
528 SUBROUTINE CMUMPS_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 COMPLEX, INTENT(INOUT) :: ARRAYPIV(LPIV,LPIVCOL)
541 COMPLEX, 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 COMPLEX, ALLOCATABLE, DIMENSION(:) :: TEMP_BLOCK
555 COMPLEX, ALLOCATABLE, DIMENSION(:) :: DEST_ARRAY
556 INTEGER :: allocok
557 COMPLEX :: ONE, MONE, ZERO
558 PARAMETER (ONE=(1.0E0,0.0E0), MONE=(-1.0E0,0.0E0))
559 PARAMETER (ZERO=(0.0E0,0.0E0))
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 & cmumps_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.LT. IF (IFLAG0) 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.GT. IF (K0) THEN
607 IF (IS_T2_SLAVE) THEN
608 POSBLOCK = POSCB +int(IBEG_BLOCK-1,8)
609 CALL cgemm('t', 'n', K, NRHS_B, M, ONE,
610 & LRB%Q(1,1), M,
611 & ARRAYCB(POSBLOCK), LDCB, ZERO,
612 & TEMP_BLOCK(1), K)
613.LE..AND..GT. ELSE IF (IBEG_BLOCKNPIVIEND_BLOCKNPIV) THEN
614 POSBLOCK = POSPIV+int(IBEG_BLOCK-1,8)
615 CALL cgemm('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 cgemm('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.LE. ELSEIF (IBEG_BLOCKNPIV) THEN
625 POSBLOCK = POSPIV+int(IBEG_BLOCK-1,8)
626 CALL cgemm('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 cgemm('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 cgemm('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 cgemm('t', 'n', N, NRHS_B, M, MONE,
646 & LRB%Q(1,1), M, ARRAYCB(POSBLOCK), LDCB,
647 & ONE, DEST_ARRAY(1), N)
648.LE..AND..GT. ELSE IF (IBEG_BLOCKNPIVIEND_BLOCKNPIV) THEN
649 POSBLOCK = POSPIV+int(IBEG_BLOCK-1,8)
650 CALL cgemm('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 cgemm('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.LE. ELSEIF (IBEG_BLOCKNPIV) THEN
657 POSBLOCK = POSPIV+int(IBEG_BLOCK-1,8)
658 CALL cgemm('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 cgemm('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.GT. IF (KMAX0) 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 caxpy(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 caxpy(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 CMUMPS_SOL_BWD_BLR_UPDATE
693 END MODULE CMUMPS_SOL_LR
694 SUBROUTINE CMUMPS_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 COMPLEX, INTENT(IN) :: DIAG(LDIAG)
705 COMPLEX, INTENT(INOUT) :: W(LWC)
706 COMPLEX RHSCOMP(LRHSCOMP,NRHS)
707 INTEGER :: LDAJ
708 COMPLEX ONE
709 PARAMETER ( ONE=(1.0E0,0.0E0) )
710.eq. IF ( MTYPE 1 ) THEN
711 LDAJ = NPIV + NELIM
712 CALL ctrsm('l','l','t','n', NPIV, NRHS_B, ONE, DIAG(1),
713 & LDAJ, RHSCOMP(PPIVINRHSCOMP,JBDEB),
714 & LRHSCOMP)
715 ELSE
716.EQ. IF ( KEEP(50) 0 ) THEN
717 LDAJ=NPIV+NELIM
718 ELSE
719 LDAJ=NPIV
720 ENDIF
721 CALL ctrsm('l','u','n','u', NPIV, NRHS_B, ONE, DIAG(1),
722 & LDAJ, RHSCOMP(PPIVINRHSCOMP,JBDEB), LRHSCOMP)
723 END IF
724 RETURN
725 END SUBROUTINE CMUMPS_SOLVE_BWD_LR_TRSOLVE
cmumps_solve_fwd_trsolve
subroutine cmumps_solve_fwd_trsolve(a, la, apos, npiv, ldadiag, nrhs_b, wcb, lwcb, lda_wcb, ppiv_courant, mtype, keep)
Definition csol_aux.F:1148
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:187
cmumps_lr_core
Definition clr_core.F:18
cmumps_lr_data_m
Definition cmumps_lr_data_m.F:14
cmumps_lr_data_m::blr_array
type(blr_struc_t), dimension(:), pointer, save, public blr_array
Definition cmumps_lr_data_m.F:61
cmumps_lr_stats
Definition clr_stats.F:14
cmumps_lr_type
Definition clr_type.F:14
cmumps_sol_lr
Definition csol_lr.F:14
cmumps_sol_lr::cmumps_sol_bwd_blr_update
subroutine cmumps_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 csol_lr.F:537
cmumps_sol_lr::cmumps_sol_bwd_lr_su
subroutine cmumps_sol_bwd_lr_su(inode, iwhdlr, npiv_global, nslaves, liell, wcb, lwcb, nrhs_b, ptwcb, rhscomp, lrhscomp, nrhs, iposinrhscomp, jbdeb, mtype, keep, keep8, iflag, ierror)
Definition csol_lr.F:386
cmumps_sol_lr::cmumps_sol_fwd_lr_su
subroutine cmumps_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 csol_lr.F:31
cmumps_sol_lr::cmumps_sol_fwd_blr_update
subroutine cmumps_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 csol_lr.F:254
cmumps_lr_type::lrb_type
Definition clr_type.F:16