OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ssol_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 smumps_sol_lr
15 USE smumps_lr_type
16 USE smumps_lr_core
17 USE smumps_lr_stats
18 USE smumps_lr_data_m, only: blr_array
19 IMPLICIT NONE
20 CONTAINS
21 SUBROUTINE smumps_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 REAL, INTENT(INOUT) :: WCB(LWCB)
39 INTEGER, INTENT(INOUT) :: IFLAG, IERROR
40 REAL, 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 REAL, POINTER, DIMENSION(:) :: DIAG
48 TYPE(lrb_type), POINTER, DIMENSION(:) :: BLR_PANEL
49 REAL :: ONE, MONE, ZERO
50 parameter(one = 1.0e0, mone=-1.0e0)
51 parameter(zero=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 SMUMPS_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 smumps_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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 smumps_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 IF (iflag.LT.0) RETURN
166 CALL smumps_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 smumps_sol_fwd_lr_su
182 SUBROUTINE smumps_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 REAL, 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 REAL :: ONE, MONE, ZERO
200 parameter(one = 1.0e0, mone=-1.0e0)
201 parameter(zero=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 SMUMPS_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 IF (mtype.EQ.1) THEN
218 CALL smumps_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 smumps_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 IF (mtype .EQ. 1) THEN
237 ptrx = ptrx + blr_panel(1)%N
238 ELSE
239 ptry = ptry + blr_panel(1)%N
240 ENDIF
241 IF (iflag.LT.0) RETURN
242 ENDIF
243 ENDDO
244 RETURN
245 END SUBROUTINE smumps_sol_slave_lr_u
246 SUBROUTINE smumps_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 REAL, INTENT(INOUT) :: ARRAYPIV(LPIV,LPIVCOL)
258 REAL, 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 REAL, ALLOCATABLE,DIMENSION(:) :: TEMP_BLOCK
273 REAL :: ONE, MONE, ZERO
274 parameter(one = 1.0e0, mone=-1.0e0)
275 parameter(zero=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 IF (kmax.GT.0) THEN
284 allocate(temp_block(kmax*nrhs_b), stat=allocok )
285 IF (allocok .GT. 0) THEN
286 iflag = -13
287 ierror = nrhs_b * kmax
288 write(*,*) 'Allocation problem in BLR routine
289 & SMUMPS_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 IF (iflag.LT.0) cycle
303 ibeg_block = begs_blr_static(i)
304 iend_block = begs_blr_static(i+1)-1
305 IF (ibeg_block .EQ. 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 IF (k.GT.0) THEN
312 CALL sgemm('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 sgemm('N', 'N', m, nrhs_b, k, mone,
318 & lrb%Q(1,1), m, temp_block(1), k,
319 & one, arraycb(posblock), ldcb)
320 ELSEIF (ibeg_block.LE.npiv.AND.iend_block.GT.npiv) THEN
321 posblock = pospiv+int(ibeg_block-1,8)
322 CALL sgemm('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 sgemm('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 ELSEIF (ibeg_block.LE.npiv) THEN
329 posblock = pospiv+int(ibeg_block-1,8)
330 CALL sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 smumps_sol_fwd_blr_update
379 SUBROUTINE smumps_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 REAL, INTENT(INOUT) :: WCB(LWCB)
394 REAL 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 REAL, POINTER, DIMENSION(:) :: DIAG
402 TYPE(lrb_type), POINTER, DIMENSION(:) :: BLR_PANEL
403 REAL :: ONE, MONE, ZERO
404 PARAMETER (ONE = 1.0e0, mone=-1.0e0)
405 parameter(zero=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 SMUMPS_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 smumps_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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 IF (iflag.LT.0) RETURN
519 CALL smumps_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 smumps_sol_bwd_lr_su
528 SUBROUTINE smumps_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 REAL, INTENT(INOUT) :: ARRAYPIV(LPIV,LPIVCOL)
541 REAL, 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 REAL, ALLOCATABLE, DIMENSION(:) :: TEMP_BLOCK
555 REAL, ALLOCATABLE, DIMENSION(:) :: DEST_ARRAY
556 INTEGER :: allocok
557 REAL :: ONE, MONE, ZERO
558 PARAMETER (ONE = 1.0e0, mone=-1.0e0)
559 parameter(zero=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 IF (current_blr.LT.last_blr) THEN
565 n = blr_panel(1)%N
566 ELSE
567 RETURN
568 ENDIF
569 allocate(dest_array(n*nrhs_b),stat=allocok)
570 IF (allocok .GT. 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 IF (kmax.GT.0) THEN
580 allocate(temp_block(kmax*nrhs_b), stat=allocok )
581 IF (allocok .GT. 0) THEN
582 iflag = -13
583 ierror = nrhs_b * kmax
584 write(*,*) 'Allocation problem in BLR routine
585 & SMUMPS_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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 sgemm('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 saxpy(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 saxpy(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 SMUMPS_SOL_BWD_BLR_UPDATE
693 END MODULE SMUMPS_SOL_LR
694 SUBROUTINE SMUMPS_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 REAL, INTENT(IN) :: DIAG(LDIAG)
705 REAL, INTENT(INOUT) :: W(LWC)
706 REAL RHSCOMP(LRHSCOMP,NRHS)
707 INTEGER :: LDAJ
708 REAL ONE
709 PARAMETER (ONE = 1.0E0)
710.eq. IF ( MTYPE 1 ) THEN
711 LDAJ = NPIV + NELIM
712 CALL strsm('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 strsm('l','u','n','u', NPIV, NRHS_B, ONE, DIAG(1),
722 & LDAJ, RHSCOMP(PPIVINRHSCOMP,JBDEB), LRHSCOMP)
723 END IF
724 RETURN
725 END SUBROUTINE SMUMPS_SOLVE_BWD_LR_TRSOLVE
mumps_abort
#define mumps_abort
Definition VE_Metis.h:25
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:187
max
#define max(a, b)
Definition macros.h:21
smumps_lr_core
Definition slr_core.F:18
smumps_lr_data_m
Definition smumps_lr_data_m.F:14
smumps_lr_data_m::blr_array
type(blr_struc_t), dimension(:), pointer, save, public blr_array
Definition smumps_lr_data_m.F:61
smumps_lr_stats
Definition slr_stats.F:14
smumps_lr_type
Definition slr_type.F:14
smumps_sol_lr
Definition ssol_lr.F:14
smumps_sol_lr::smumps_sol_bwd_blr_update
subroutine smumps_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 ssol_lr.F:537
smumps_sol_lr::smumps_sol_fwd_lr_su
subroutine smumps_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 ssol_lr.F:31
smumps_sol_lr::smumps_sol_bwd_lr_su
subroutine smumps_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 ssol_lr.F:386
smumps_sol_lr::smumps_sol_fwd_blr_update
subroutine smumps_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 ssol_lr.F:254
smumps_sol_lr::smumps_sol_slave_lr_u
subroutine smumps_sol_slave_lr_u(inode, iwhdlr, npiv_global, wcb, lwcb, ldx, ldy, ptrx_init, ptry_init, jbdeb, jbfin, mtype, keep, keep8, iflag, ierror)
Definition ssol_lr.F:189
smumps_sol_ld_and_reload_panel
subroutine smumps_sol_ld_and_reload_panel(inode, n, npiv, liell, nelim, nslaves, ppiv_courant, iw, ipos, liw, a, la, apos, wcb, lwcb, ld_wcbpiv, rhscomp, lrhscomp, nrhs, posinrhscomp_fwd, jbdeb, jbfin, mtype, keep, oocwrite_compatible_with_blr, ignore_k459)
Definition ssol_aux.F:1379
smumps_solve_fwd_trsolve
subroutine smumps_solve_fwd_trsolve(a, la, apos, npiv, ldadiag, nrhs_b, wcb, lwcb, lda_wcb, ppiv_courant, mtype, keep)
Definition ssol_aux.F:1147
smumps_solve_bwd_lr_trsolve
subroutine smumps_solve_bwd_lr_trsolve(diag, ldiag, npiv, nelim, liell, nrhs_b, w, lwc, rhscomp, lrhscomp, nrhs, ppivinrhscomp, jbdeb, mtype, keep)
Definition ssol_lr.F:700
smumps_lr_type::lrb_type
Definition slr_type.F:16