OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dfac_front_aux.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_fac_front_aux_m
15 CONTAINS
16 SUBROUTINE dmumps_fac_h(NFRONT,NASS,IW,LIW,A,LA,
17 & INOPV,NOFFW,
18 & DET_EXPW, DET_MANTW, DET_SIGNW,
19 & IOLDPS,POSELT,UU,SEUIL,KEEP, KEEP8, DKEEP,
20 & PP_FIRST2SWAP_L, PP_LastPanelonDisk_L,
21 & PP_LastPIVRPTRFilled_L,
22 & PP_FIRST2SWAP_U, PP_LastPanelonDisk_U,
23 & PP_LastPIVRPTRFilled_U,MAXFROMN,IS_MAXFROMN_AVAIL,
24 & Inextpiv, OOC_EFFECTIVE_ON_FRONT, NVSCHUR
25 &)
26!$ USE OMP_LIB
27 USE mumps_ooc_common
28 IMPLICIT NONE
29 INTEGER NFRONT,NASS,LIW,INOPV
30 INTEGER(8) :: LA
31 INTEGER :: KEEP(500)
32 INTEGER(8) :: KEEP8(150)
33 DOUBLE PRECISION :: DKEEP(230)
34 DOUBLE PRECISION UU, SEUIL
35 DOUBLE PRECISION A(LA)
36 INTEGER IW(LIW)
37 DOUBLE PRECISION, intent(in) :: MAXFROMN
38 LOGICAL, intent(inout) :: IS_MAXFROMN_AVAIL
39 INTEGER, intent(inout) :: Inextpiv
40 LOGICAL, intent(in) :: OOC_EFFECTIVE_ON_FRONT
41 INTEGER, intent(in) :: NVSCHUR
42 DOUBLE PRECISION AMROW
43 DOUBLE PRECISION RMAX, SEUIL_LOC
44 DOUBLE PRECISION SWOP
45 INTEGER(8) :: APOS, POSELT
46 INTEGER(8) :: J1, J2, J3_8, JJ, IDIAG
47 INTEGER(8) :: J1_ini
48 INTEGER(8) :: NFRONT8
49 INTEGER IOLDPS
50 INTEGER NPIV,IPIV,IPIV_SHIFT
51 INTEGER, intent(inout) :: NOFFW
52 INTEGER, intent(inout) :: DET_EXPW, DET_SIGNW
53 DOUBLE PRECISION, intent(inout) :: DET_MANTW
54 INTEGER J, J3
55 INTEGER NPIVP1,JMAX,ISW,ISWPS1
56 INTEGER ISWPS2,KSW,XSIZE
57 INTEGER I_PIVRPTR_L, I_PIVR_L, NBPANELS_L
58 INTEGER I_PIVRPTR_U, I_PIVR_U, NBPANELS_U
59 INTEGER PP_FIRST2SWAP_L, PP_LastPanelonDisk_L,
60 & pp_lastpivrptrfilled_l,
61 & pp_first2swap_u, pp_lastpanelondisk_u,
62 & pp_lastpivrptrfilled_u
63 INTEGER ISHIFT, K206
64 INTEGER DMUMPS_IXAMAX
65 include 'mumps_headers.h'
66 INTRINSIC max
67 DOUBLE PRECISION, PARAMETER :: RZERO = 0.0d0
68#if defined(_OPENMP)
69 INTEGER :: NOMP, CHUNK, K360
70 K360 = keep(360)
71 nomp = omp_get_max_threads()
72#endif
73 seuil_loc = max(dkeep(1), seuil)
74 nfront8 = int(nfront,8)
75 inopv = 0
76 xsize = keep(ixsz)
77 npiv = iw(ioldps+1+xsize)
78 npivp1 = npiv + 1
79 k206 = keep(206)
80 IF ((keep(50).NE.1).AND.ooc_effective_on_front) THEN
81 CALL dmumps_get_ooc_perm_ptr(typef_l, nbpanels_l,
82 & i_pivrptr_l, i_pivr_l,
83 & ioldps+2*nfront+6+iw(ioldps+5+xsize)
84 & +keep(ixsz),
85 & iw, liw)
86 CALL dmumps_get_ooc_perm_ptr(typef_u, nbpanels_u,
87 & i_pivrptr_u, i_pivr_u,
88 & ioldps+2*nfront+6+iw(ioldps+5+xsize)+xsize,
89 & iw, liw)
90 ENDIF
91 ishift = 0
92 IF (k206.GE.1) THEN
93 IF (inextpiv.GT.npivp1.AND.inextpiv.LE.nass) THEN
94 ishift = inextpiv - npivp1
95 ENDIF
96 IF (ishift.GT.0.AND.is_maxfromn_avail) THEN
97 ipiv = npivp1
98 apos = poselt + nfront8*int(npiv,8) + int(ipiv-1,8)
99 idiag = apos + int(ipiv - npivp1,8)*nfront8
100 IF (abs(a(idiag)) .GE. uu*maxfromn .AND.
101 & abs(a(idiag)) .GT. max(seuil_loc,tiny(rmax))
102 & ) THEN
103 ishift = 0
104 ENDIF
105 ENDIF
106 IF ( ishift .GT. 0) THEN
107 is_maxfromn_avail = .false.
108 ENDIF
109 ENDIF
110 DO 460 ipiv_shift=npivp1+ishift,nass+ishift
111 IF (ipiv_shift .LE. nass) THEN
112 ipiv=ipiv_shift
113 ELSE
114 ipiv=ipiv_shift-nass-1+npivp1
115 ENDIF
116 apos = poselt + nfront8*int(npiv,8) + int(ipiv-1,8)
117 jmax = 1
118 amrow = rzero
119 j1 = apos
120 j3 = nass -npiv
121 jmax = dmumps_ixamax(j3,a(j1),nfront,keep(360))
122 jj = j1 + int(jmax-1,8)*nfront8
123 amrow = abs(a(jj))
124 rmax = amrow
125 j1 = apos + int(nass-npiv,8) * nfront8
126 j3 = nfront - nass - keep(253)-nvschur
127 IF (is_maxfromn_avail) THEN
128 rmax = max(maxfromn,rmax)
129 is_maxfromn_avail = .false.
130 ELSE
131 IF (j3.EQ.0) GOTO 370
132 IF (keep(351).EQ.1) THEN
133 j1_ini = j1
134!$ CHUNK = max(K360/2,(J3+NOMP-1)/NOMP)
135!$OMP PARALLEL DO schedule(static, CHUNK)
136!$OMP& FIRSTPRIVATE(J1_ini,NFRONT8,J3)
137!$omp& reduction(max:rmax) IF (j3.GE.k360)
138 DO j=1,j3
139 rmax = max(abs(a(j1_ini + int(j-1,8) * nfront8)),
140 & rmax)
141 END DO
142!$OMP END PARALLEL DO
143 ELSE
144 DO j=1,j3
145 rmax = max(abs(a(j1)), rmax)
146 j1 = j1 + nfront8
147 END DO
148 ENDIF
149 END IF
150 370 IF (rmax.LE.tiny(rmax)) GO TO 460
151 idiag = apos + int(ipiv - npivp1,8)*nfront8
152 IF (abs(a(idiag)) .GE. uu*rmax .AND.
153 & abs(a(idiag)) .GT. max(seuil_loc,tiny(rmax)) ) THEN
154 jmax = ipiv - npiv
155 GO TO 380
156 ENDIF
157 IF ( .NOT. ( amrow .GE. uu*rmax .AND.
158 & amrow .GT. max(seuil_loc,tiny(rmax))
159 & )
160 & ) GO TO 460
161 noffw = noffw + 1
162 380 CONTINUE
163 IF (k206.GE.1) THEN
164 inextpiv = ipiv + 1
165 ENDIF
166 CALL dmumps_update_minmax_pivot
167 & ( abs(a(apos + int(jmax - 1,8) * nfront8 )),
168 & dkeep, keep, .false.)
169 IF (keep(258) .NE. 0) THEN
170 CALL dmumps_updatedeter(
171 & a(apos + int(jmax - 1,8) * nfront8 ),
172 & det_mantw, det_expw )
173 ENDIF
174 IF (ipiv.EQ.npivp1) GO TO 400
175 IF (keep(405) .EQ.0) THEN
176 keep8(80) = keep8(80)+1
177 ELSE
178!$OMP ATOMIC UPDATE
179 keep8(80) = keep8(80)+1
180!$OMP END ATOMIC
181 ENDIF
182 det_signw = - det_signw
183 j1 = poselt + int(npiv,8)
184 j3_8 = poselt + int(ipiv-1,8)
185 DO j= 1,nfront
186 swop = a(j1)
187 a(j1) = a(j3_8)
188 a(j3_8) = swop
189 j1 = j1 + nfront8
190 j3_8 = j3_8 + nfront8
191 END DO
192 iswps1 = ioldps + 5 + npivp1 + nfront + xsize
193 iswps2 = ioldps + 5 + ipiv + nfront + xsize
194 isw = iw(iswps1)
195 iw(iswps1) = iw(iswps2)
196 iw(iswps2) = isw
197 400 IF (jmax.EQ.1) GO TO 420
198 det_signw = -det_signw
199 j1 = poselt + int(npiv,8) * nfront8
200 j2 = poselt + int(npiv + jmax - 1,8) * nfront8
201 DO ksw=1,nfront
202 swop = a(j1)
203 a(j1) = a(j2)
204 a(j2) = swop
205 j1 = j1 + 1_8
206 j2 = j2 + 1_8
207 END DO
208 iswps1 = ioldps + 5 + npiv + 1 + xsize
209 iswps2 = ioldps + 5 + npiv + jmax + xsize
210 isw = iw(iswps1)
211 iw(iswps1) = iw(iswps2)
212 iw(iswps2) = isw
213 GO TO 420
214 460 CONTINUE
215 inopv = 1
216 GOTO 430
217 420 CONTINUE
218 IF (ooc_effective_on_front) THEN
219 IF (keep(251).EQ.0) THEN
220 CALL dmumps_store_perminfo( iw(i_pivrptr_l),
221 & nbpanels_l,
222 & iw(i_pivr_l), nass, npivp1, npiv+jmax,
223 & pp_lastpanelondisk_l,
224 & pp_lastpivrptrfilled_l)
225 ENDIF
226 CALL dmumps_store_perminfo( iw(i_pivrptr_u),
227 & nbpanels_u,
228 & iw(i_pivr_u), nass, npivp1, ipiv,
229 & pp_lastpanelondisk_u,
230 & pp_lastpivrptrfilled_u)
231 ENDIF
232 430 CONTINUE
233 is_maxfromn_avail = .false.
234 RETURN
235 END SUBROUTINE dmumps_fac_h
236 SUBROUTINE dmumps_fac_n(NFRONT,NASS,IW,LIW,A,LA,
237 & IOLDPS,POSELT,IFINB,XSIZE,
238 & KEEP,MAXFROMN,IS_MAXFROMN_AVAIL,NVSCHUR)
239!$ USE OMP_LIB
240 IMPLICIT NONE
241 include 'mumps_headers.h'
242 INTEGER NFRONT,NASS,LIW,IFINB
243 INTEGER(8) :: LA
244 DOUBLE PRECISION A(LA)
245 INTEGER IW(LIW)
246 DOUBLE PRECISION ALPHA,VALPIV
247 INTEGER(8) :: APOS, POSELT, UUPOS, LPOS, IRWPOS
248 INTEGER(8) :: NFRONT8
249 INTEGER IOLDPS,NPIV,XSIZE
250 INTEGER, intent(in) :: KEEP(500)
251 DOUBLE PRECISION, intent(inout) :: MAXFROMN
252 LOGICAL, intent(inout) :: IS_MAXFROMN_AVAIL
253 INTEGER, intent(in) :: NVSCHUR
254 INTEGER NEL,IROW,NEL2,JCOL,NELMAXM
255 INTEGER NPIVP1
256 DOUBLE PRECISION, PARAMETER :: ONE = 1.0d0
257#if defined(_OPENMP)
258 LOGICAL:: OMP_FLAG
259 INTEGER:: NOMP, K360, CHUNK
260 nomp = omp_get_max_threads()
261 k360 = keep(360)
262#endif
263 nfront8=int(nfront,8)
264 npiv = iw(ioldps+1+xsize)
265 npivp1 = npiv + 1
266 nel = nfront - npivp1
267 nelmaxm= nel -keep(253)-nvschur
268 nel2 = nass - npivp1
269 ifinb = 0
270 IF (npivp1.EQ.nass) ifinb = 1
271 apos = poselt + int(npiv,8)*(nfront8 + 1_8)
272 valpiv = one/a(apos)
273#if defined(_OPENMP)
274 omp_flag = .false.
275 chunk = max(nel,1)
276 IF (nomp.GT.1) THEN
277 IF (nel.LT.k360) THEN
278 IF (nel*nel2.GE.keep(361)) THEN
279 omp_flag = .true.
280 chunk = max(20, (nel+nomp-1)/nomp)
281 ENDIF
282 ELSE
283 omp_flag = .true.
284 chunk = max(k360/2, (nel+nomp-1)/nomp)
285 ENDIF
286 ENDIF
287#endif
288 IF (keep(351).EQ.2) THEN
289 maxfromn = 0.0d0
290 IF (nel2 > 0) THEN
291 is_maxfromn_avail = .true.
292 ENDIF
293!$OMP PARALLEL DO schedule(static, CHUNK)
294!$OMP& PRIVATE(LPOS, UUPOS, IRWPOS, ALPHA, JCOL)
295!$omp& firstprivate(apos,nfront8,valpiv,nel,nel2)
296!$omp& reduction(max:maxfromn) IF(omp_flag)
297 DO irow = 1, nel
298 lpos = apos + nfront8*int(irow,8)
299 a(lpos) = a(lpos)*valpiv
300 alpha = -a(lpos)
301 irwpos = lpos + 1_8
302 uupos = apos + 1_8
303 IF (nel2 > 0) THEN
304 a(irwpos) = a(irwpos) + alpha*a(uupos)
305 IF (irow.LE.nelmaxm)
306 & maxfromn=max(maxfromn, abs(a(irwpos)))
307 irwpos = irwpos+1_8
308 uupos = uupos+1_8
309 DO jcol = 2, nel2
310 a(irwpos) = a(irwpos) + alpha*a(uupos)
311 irwpos = irwpos+1_8
312 uupos = uupos+1_8
313 ENDDO
314 ENDIF
315 END DO
316!$OMP END PARALLEL DO
317 ELSE
318!$OMP PARALLEL DO schedule(static, CHUNK)
319!$OMP& FIRSTPRIVATE(APOS,NFRONT8,VALPIV,NEL,NEL2)
320!$OMP& PRIVATE(LPOS, UUPOS, IRWPOS, ALPHA, JCOL) IF(OMP_FLAG)
321 DO irow = 1, nel
322 lpos = apos + nfront8*int(irow,8)
323 a(lpos) = a(lpos)*valpiv
324 alpha = -a(lpos)
325 irwpos = lpos + 1_8
326 uupos = apos + 1_8
327 DO jcol = 1, nel2
328 a(irwpos) = a(irwpos) + alpha*a(uupos)
329 irwpos = irwpos+1_8
330 uupos = uupos+1_8
331 ENDDO
332 ENDDO
333!$OMP END PARALLEL DO
334 ENDIF
335 RETURN
336 END SUBROUTINE dmumps_fac_n
337 SUBROUTINE dmumps_fac_pt_setlock427( K427_OUT, K427,
338 & K405, K222, NEL1, NASS )
339 INTEGER, INTENT(IN) :: K427, K405, K222, NEL1, NASS
340 INTEGER, INTENT(OUT) :: K427_OUT
341 k427_out = k427
342 IF ( k427_out .GT. 0 ) k427_out = 0
343 IF ( k427_out .LT. 0 ) k427_out = -1
344 IF ( k427_out .GT. 99 ) k427_out = 0
345 IF ( k427_out .LT. -100 ) k427_out = -1
346 RETURN
347 END SUBROUTINE dmumps_fac_pt_setlock427
348 SUBROUTINE dmumps_fac_p(A,LA,NFRONT,
349 & NPIV,NASS,POSELT,CALL_UTRSM, KEEP, INODE,
350 & CALL_OOC, IWFAC, LIWFAC, LAFAC, MonBloc, MYID, KEEP8,
351 & LNextPiv2beWritten, UNextPiv2beWritten,
352 & IFLAG )
353 USE dmumps_ooc, ONLY : io_block, typef_both_lu,
354 & dmumps_ooc_io_lu_panel
355 USE mumps_ooc_common, ONLY : strat_try_write
356 IMPLICIT NONE
357 INTEGER(8) :: LA,POSELT,LAFAC
358 DOUBLE PRECISION A(LA)
359 INTEGER NFRONT, NPIV, NASS
360 LOGICAL, INTENT(IN) :: CALL_UTRSM
361 INTEGER, INTENT(INOUT) :: IFLAG
362 LOGICAL, INTENT(IN) :: CALL_OOC
363 INTEGER LIWFAC, MYID,
364 & lnextpiv2bewritten, unextpiv2bewritten
365 INTEGER IWFAC(LIWFAC)
366 TYPE(IO_BLOCK) :: MonBloc
367 INTEGER :: KEEP(500)
368 INTEGER(8) :: KEEP8(150)
369 INTEGER(8) :: LPOS, LPOS1, LPOS2, UPOS
370 INTEGER NEL1, NEL11, IFLAG_OOC
371 INTEGER :: INODE
372 DOUBLE PRECISION ALPHA, ONE
373 parameter(one = 1.0d0, alpha=-1.0d0)
374 include 'mumps_headers.h'
375 nel1 = nfront - nass
376 nel11 = nfront - npiv
377 lpos2 = poselt + int(nass,8)*int(nfront,8)
378 lpos = lpos2 + int(npiv,8)
379 lpos1 = poselt + int(npiv,8)
380 upos = poselt + int(nass,8)
381 IF ( call_utrsm ) THEN
382 CALL dtrsm('R', 'U', 'N', 'U', nel1, npiv, one,
383 & a(poselt), nfront, a(upos), nfront)
384 ENDIF
385 CALL dtrsm('L','L','N','N',npiv,nel1,one,a(poselt),nfront,
386 & a(lpos2),nfront)
387 IF (call_ooc) THEN
388 CALL dmumps_ooc_io_lu_panel
389 & ( strat_try_write, typef_both_lu,
390 & a(poselt), lafac, monbloc,
391 & lnextpiv2bewritten, unextpiv2bewritten,
392 & iwfac, liwfac,
393 & myid, keep8(31), iflag_ooc,
394 & .false. )
395 IF (iflag_ooc .LT. 0) THEN
396 iflag = iflag_ooc
397 GOTO 500
398 ENDIF
399 ENDIF
400 CALL dgemm('N','N',nel11,nel1,npiv,alpha,a(lpos1),
401 & nfront,a(lpos2),nfront,one,a(lpos),nfront)
402 IF ((call_utrsm).AND.(nass-npiv.GT.0)) THEN
403 lpos2 = poselt + int(npiv,8)*int(nfront,8)
404 lpos = lpos2 + int(nass,8)
405 CALL dgemm('N','N',nel1,nass-npiv,npiv,alpha,a(upos),
406 & nfront,a(lpos2),nfront,one,a(lpos),nfront)
407 ENDIF
408 500 CONTINUE
409 RETURN
410 END SUBROUTINE dmumps_fac_p
411 SUBROUTINE dmumps_fac_t(A,LA,NPIVB,NFRONT,
412 & NPIV,NASS,POSELT)
413 IMPLICIT NONE
414 INTEGER NPIVB,NASS
415 INTEGER(8) :: LA
416 DOUBLE PRECISION A(LA)
417 INTEGER(8) :: APOS, POSELT
418 INTEGER NFRONT, NPIV, NASSL
419 INTEGER(8) :: LPOS, LPOS1, LPOS2
420 INTEGER NEL1, NEL11, NPIVE
421 DOUBLE PRECISION ALPHA, ONE
422 parameter(one = 1.0d0, alpha=-1.0d0)
423 nel1 = nfront - nass
424 nel11 = nfront - npiv
425 npive = npiv - npivb
426 nassl = nass - npivb
427 apos = poselt + int(npivb,8)*int(nfront,8)
428 & + int(npivb,8)
429 lpos2 = apos + int(nassl,8)
430 CALL dtrsm('R','U','N','U',nel1,npive,one,a(apos),nfront,
431 & a(lpos2),nfront)
432 lpos = lpos2 + int(nfront,8)*int(npive,8)
433 lpos1 = apos + int(nfront,8)*int(npive,8)
434 CALL dgemm('N','N',nel1,nel11,npive,alpha,a(lpos2),
435 & nfront,a(lpos1),nfront,one,a(lpos),nfront)
436 RETURN
437 END SUBROUTINE dmumps_fac_t
438 SUBROUTINE dmumps_fac_sq(IBEG_BLOCK, IEND_BLOCK, NPIV,
439 & NFRONT, LAST_ROW, LAST_COL, A, LA, POSELT,
440 & FIRST_COL, CALL_LTRSM, CALL_UTRSM, CALL_GEMM,
441 & WITH_COMM_THREAD, LR_ACTIVATED
442 & )
443!$ USE OMP_LIB
444#if defined(_OPENMP)
445 USE dmumps_buf
446#endif
447 IMPLICIT NONE
448 INTEGER, intent(in) :: IBEG_BLOCK, IEND_BLOCK
449 INTEGER, intent(in) :: NPIV, NFRONT, LAST_ROW, LAST_COL
450 INTEGER, intent(in) :: FIRST_COL
451 INTEGER(8), intent(in) :: LA
452 DOUBLE PRECISION, intent(inout) :: A(LA)
453 INTEGER(8), intent(in) :: POSELT
454 LOGICAL, intent(in) :: CALL_LTRSM, CALL_UTRSM, CALL_GEMM
455 LOGICAL, intent(in) :: WITH_COMM_THREAD, LR_ACTIVATED
456 INTEGER(8) :: NFRONT8, LPOSN, LPOS2N
457 INTEGER(8) :: LPOS, LPOS1, LPOS2, UPOS, POSELT_LOCAL
458 INTEGER :: NELIM, LKJIW, NEL1, NEL11, UTRSM_NCOLS
459 DOUBLE PRECISION ALPHA, ONE
460 PARAMETER (ONE = 1.0d0, alpha=-1.0d0)
461!$ INTEGER :: NOMP
462!$ LOGICAL :: TRSM_GEMM_FINISHED
463!$ LOGICAL :: SAVE_NESTED, SAVE_DYNAMIC
464 nfront8= int(nfront,8)
465 nelim = iend_block - npiv
466 nel1 = last_row - iend_block
467 IF ( nel1 < 0 ) THEN
468 WRITE(*,*)
469 & "Internal error 1 in DMUMPS_FAC_SQ,IEND_BLOCK>LAST_ROW",
470 & iend_block, last_row
471 CALL mumps_abort()
472 ENDIF
473 lkjiw = npiv - ibeg_block + 1
474 nel11 = last_col - npiv
475 lpos2 = poselt + int(iend_block,8)*nfront8 + int(ibeg_block-1,8)
476 utrsm_ncols = last_col - first_col
477 upos = poselt + int(ibeg_block-1,8)*nfront8 + int(first_col,8)
478 poselt_local = poselt + int(ibeg_block-1,8)*nfront8
479 & + int(ibeg_block-1,8)
480 IF ((nel1.NE.0).AND.(lkjiw.NE.0)) THEN
481 IF (with_comm_thread .EQV. .false.) THEN
482 IF (call_ltrsm) THEN
483 CALL dtrsm('L','L','N','N',lkjiw,nel1,one,
484 & a(poselt_local),nfront,
485 & a(lpos2),nfront)
486 ENDIF
487 IF (call_utrsm) THEN
488 CALL dtrsm('R','U','N','U',utrsm_ncols,lkjiw,one,
489 & a(poselt_local),nfront,
490 & a(upos),nfront)
491 lpos2n = poselt + int(npiv,8)*nfront8 + int(ibeg_block-1,8)
492 lposn = poselt + int(npiv,8)*nfront8 + int(first_col,8)
493 CALL dgemm('N','N',utrsm_ncols,nelim,
494 & lkjiw,alpha,a(upos),nfront,a(lpos2n),
495 & nfront,one,a(lposn),nfront)
496 ENDIF
497 IF (call_gemm) THEN
498 lpos = lpos2 + int(lkjiw,8)
499 lpos1 = poselt_local + int(lkjiw,8)
500 CALL dgemm('N','N',nel11,nel1,lkjiw,alpha,a(lpos1),
501 & nfront,a(lpos2),nfront,one,a(lpos),nfront)
502 ENDIF
503 ELSE
504!$ NOMP = OMP_GET_MAX_THREADS()
505!$ CALL OMP_SET_NUM_THREADS(2)
506!$ SAVE_NESTED = OMP_GET_NESTED()
507!$ SAVE_DYNAMIC = OMP_GET_DYNAMIC()
508!$ CALL OMP_SET_NESTED(.TRUE.)
509!$ CALL OMP_SET_DYNAMIC(.FALSE.)
510!$ TRSM_GEMM_FINISHED = .FALSE.
511!$OMP PARALLEL SHARED(TRSM_GEMM_FINISHED)
512!$ IF (OMP_GET_THREAD_NUM() .EQ. 1) THEN
513#if defined(WORKAROUNDINTELILP64OPENMPLIMITATION)
514!$ CALL OMP_SET_NUM_THREADS(int(NOMP,4))
515#else
516!$ CALL OMP_SET_NUM_THREADS(NOMP)
517#endif
518 IF (call_ltrsm) THEN
519 CALL dtrsm('L','L','N','N',lkjiw,nel1,one,
520 & a(poselt_local),nfront,
521 & a(lpos2),nfront)
522 ENDIF
523 IF (call_utrsm) THEN
524 CALL dtrsm('R','U','N','U',utrsm_ncols,lkjiw,one,
525 & a(poselt_local),nfront,
526 & a(upos),nfront)
527 lpos2n = poselt + int(npiv,8)*nfront8 + int(ibeg_block-1,8)
528 lposn = poselt + int(npiv,8)*nfront8 + int(first_col,8)
529 CALL dgemm('N','N',utrsm_ncols,nelim,
530 & lkjiw,alpha,a(upos),nfront,a(lpos2n),
531 & nfront,one,a(lposn),nfront)
532 ENDIF
533 IF (call_gemm) THEN
534 lpos = lpos2 + int(lkjiw,8)
535 lpos1 = poselt_local + int(lkjiw,8)
536 CALL dgemm('N','N',nel11,nel1,lkjiw,alpha,a(lpos1),
537 & nfront,a(lpos2),nfront,one,a(lpos),nfront)
538 END IF
539!$ TRSM_GEMM_FINISHED = .TRUE.
540!$ ELSE
541!$ DO WHILE (.NOT. TRSM_GEMM_FINISHED)
542!$ CALL DMUMPS_BUF_TEST()
543!$ CALL MUMPS_USLEEP(10000)
544!$ END DO
545!$ END IF
546!$omp END parallel
547!$ CALL OMP_SET_NESTED(SAVE_NESTED)
548!$ CALL OMP_SET_DYNAMIC(SAVE_DYNAMIC)
549#if defined(WORKAROUNDINTELILP64OPENMPLIMITATION)
550!$ CALL omp_set_num_threads(int(nomp,4))
551#else
552!$ CALL OMP_SET_NUM_THREADS(NOMP)
553#endif
554 ENDIF
555 ELSE
556 IF (call_utrsm.AND.utrsm_ncols.NE.0) THEN
557 CALL dtrsm('R','U','N','U',utrsm_ncols,lkjiw,one,
558 & a(poselt_local),nfront,
559 & a(upos),nfront)
560 lpos2n = poselt + int(npiv,8)*nfront8 + int(ibeg_block-1,8)
561 lposn = poselt + int(npiv,8)*nfront8 + int(first_col,8)
562 CALL dgemm('N','N',utrsm_ncols,nelim,
563 & lkjiw,alpha,a(upos),nfront,a(lpos2n),
564 & nfront,one,a(lposn),nfront)
565 ENDIF
566 ENDIF
567 RETURN
568 END SUBROUTINE dmumps_fac_sq
569 SUBROUTINE dmumps_fac_mq(IBEG_BLOCK,IEND_BLOCK,
570 & NFRONT, NASS, NPIV, LAST_COL, A, LA, POSELT, IFINB,
571 & LR_ACTIVATED
572 & )
573 IMPLICIT NONE
574 INTEGER, intent(in) :: IBEG_BLOCK, IEND_BLOCK, NFRONT,
575 & NASS, NPIV, LAST_COL
576 INTEGER, intent(out) :: IFINB
577 INTEGER(8), intent(in) :: LA, POSELT
578 DOUBLE PRECISION, intent(inout) :: A(LA)
579 LOGICAL, intent(in) :: LR_ACTIVATED
580 DOUBLE PRECISION :: VALPIV
581 INTEGER(8) :: APOS, UUPOS, LPOS
582 INTEGER(8) :: NFRONT8
583 DOUBLE PRECISION :: ONE, ALPHA
584 INTEGER :: NEL2,NPIVP1,KROW,NEL
585 parameter(one = 1.0d0, alpha=-1.0d0)
586 nfront8= int(nfront,8)
587 npivp1 = npiv + 1
588 nel = last_col - npivp1
589 ifinb = 0
590 nel2 = iend_block - npivp1
591 IF (nel2.EQ.0) THEN
592 IF (iend_block.EQ.nass) THEN
593 ifinb = -1
594 ELSE
595 ifinb = 1
596 ENDIF
597 ELSE
598 apos = poselt + int(npiv,8)*(nfront8 + 1_8)
599 valpiv = one/a(apos)
600 lpos = apos + nfront8
601 DO 541 krow = 1,nel2
602 a(lpos) = a(lpos)*valpiv
603 lpos = lpos + nfront8
604 541 CONTINUE
605 lpos = apos + nfront8
606 uupos = apos + 1_8
607#if defined(MUMPS_USE_BLAS2)
608 CALL dger(nel,nel2,alpha,a(uupos),1,a(lpos),nfront,
609 & a(lpos+1_8),nfront)
610#else
611 CALL dgemm('N','N',nel,nel2,1,alpha,a(uupos),nel,
612 & a(lpos),nfront,one,a(lpos+1_8),nfront)
613#endif
614 ENDIF
615 RETURN
616 END SUBROUTINE dmumps_fac_mq
617 SUBROUTINE dmumps_fac_fr_update_cbrows( INODE, NFRONT, NASS,
618 & CALL_UTRSM, A, LA, LAFAC, POSELT, IW, LIW, IOLDPS,
619 & MonBloc, MYID, NOFFW,
620 & DET_EXPW, DET_MANTW, DET_SIGNW,
621 & LIWFAC,
622 & PP_FIRST2SWAP_L, PP_FIRST2SWAP_U,
623 & LNextPiv2beWritten, UNextPiv2beWritten,
624 & PP_LastPIVRPTRFilled_L, PP_LastPIVRPTRFilled_U,
625 &
626 & XSIZE, SEUIL, UU, DKEEP, KEEP8, KEEP, IFLAG,
627 & OOC_EFFECTIVE_ON_FRONT, NVSCHUR)
628 USE dmumps_ooc, ONLY: io_block
629 IMPLICIT NONE
630 INTEGER, intent(in) :: INODE, NFRONT, NASS,
631 & liw, myid, xsize, ioldps, liwfac
632 INTEGER(8), intent(in) :: LA, POSELT
633 INTEGER, intent(inout) :: NOFFW
634 INTEGER, intent(inout) :: DET_EXPW, DET_SIGNW
635 DOUBLE PRECISION, intent(inout) :: DET_MANTW
636 INTEGER, intent(inout) :: PP_FIRST2SWAP_L, PP_FIRST2SWAP_U,
637 & lnextpiv2bewritten, unextpiv2bewritten,
638 & pp_lastpivrptrfilled_l, pp_lastpivrptrfilled_u,
639 & iflag
640 LOGICAL, intent(in) :: CALL_UTRSM
641 INTEGER, intent(inout) :: IW(LIW)
642 DOUBLE PRECISION, intent(inout) :: A(LA)
643 DOUBLE PRECISION, intent(in) :: SEUIL, UU, DKEEP(230)
644 INTEGER, intent(in) :: KEEP( 500 )
645 INTEGER(8), intent(inout) :: LAFAC
646 INTEGER(8) :: KEEP8(150)
647 INTEGER, intent(in) :: NVSCHUR
648 TYPE(io_block), intent(inout) :: MonBloc
649 LOGICAL, intent(in) :: OOC_EFFECTIVE_ON_FRONT
650 INTEGER :: NPIV, NEL1, IBEG_BLOCK, IFINB, INOPV
651 INTEGER Inextpiv
652 DOUBLE PRECISION :: MAXFROMN
653 LOGICAL :: IS_MAXFROMN_AVAIL
654 NPIV = iw(ioldps+1+xsize)
655 nel1 = nfront - nass
656 IF (keep(206).GE.1) THEN
657 inextpiv = 1
658 ELSE
659 inextpiv = 0
660 ENDIF
661 IF ((npiv.GT.0).AND.(nel1.GT.0)) THEN
662 IF (ooc_effective_on_front) THEN
663 monbloc%LastPiv = npiv
664 ENDIF
665 CALL dmumps_fac_p(a,la,nfront, npiv, nass, poselt,
666 & call_utrsm, keep, inode,
667 & ooc_effective_on_front, iw(ioldps),
668 & liwfac, lafac,
669 & monbloc, myid, keep8,
670 & lnextpiv2bewritten, unextpiv2bewritten,
671 & iflag)
672 ENDIF
673 npiv = iw(ioldps+1+xsize)
674 ibeg_block = npiv
675 IF (nass.EQ.npiv) GOTO 500
676 is_maxfromn_avail = .false.
677 120 CALL dmumps_fac_h(nfront,nass,iw,liw,a,la,
678 & inopv, noffw,
679 & det_expw, det_mantw, det_signw,
680 & ioldps,poselt,uu,seuil,
681 & keep, keep8, dkeep,
682 & pp_first2swap_l, monbloc%LastPanelWritten_L,
683 & pp_lastpivrptrfilled_l,
684 & pp_first2swap_u, monbloc%LastPanelWritten_U,
685 & pp_lastpivrptrfilled_u, maxfromn, is_maxfromn_avail,
686 & inextpiv, ooc_effective_on_front, nvschur
687 & )
688 IF (inopv.NE.1) THEN
689 CALL dmumps_fac_n(nfront,nass,iw,liw,a,la,
690 & ioldps,poselt,ifinb,xsize,
691 & keep, maxfromn, is_maxfromn_avail,
692 & nvschur)
693 iw(ioldps+1+xsize) = iw(ioldps+1+xsize) + 1
694 IF (ifinb.EQ.0) GOTO 120
695 ENDIF
696 npiv = iw(ioldps+1+xsize)
697 nel1 = nfront - nass
698 IF ((npiv.LE.ibeg_block).OR.(nel1.EQ.0)) GO TO 500
699 CALL dmumps_fac_t(a,la,ibeg_block,
700 & nfront,npiv,nass,poselt)
701 500 CONTINUE
702 RETURN
703 END SUBROUTINE dmumps_fac_fr_update_cbrows
704 SUBROUTINE dmumps_fac_i(NFRONT,NASS,LAST_ROW,
705 & IBEG_BLOCK, IEND_BLOCK,
706 & N,INODE,IW,LIW,A,LA,
707 & INOPV,NOFFW,NBTINYW,
708 & DET_EXPW, DET_MANTW, DET_SIGNW,
709 & IFLAG,IOLDPS,POSELT,UU,SEUIL,KEEP,KEEP8,
710 & DKEEP,PIVNUL_LIST,LPN_LIST,
711 &
712 & PP_FIRST2SWAP_L, PP_LastPanelonDisk_L,
713 & PP_LastPIVRPTRFilled_L,
714 & PP_FIRST2SWAP_U, PP_LastPanelonDisk_U,
715 & PP_LastPIVRPTRFilled_U,
716 & PIVOT_OPTION, LR_ACTIVATED, IEND_BLR, Inextpiv,
717 & OOC_EFFECTIVE_ON_FRONT, NVSCHUR, PARPIV_T1,
718 & TIPIV
719 & )
720!$ USE OMP_LIB
721 USE mumps_ooc_common
722 IMPLICIT NONE
723 INTEGER, intent(in) :: IBEG_BLOCK, IEND_BLOCK
724 INTEGER, intent(inout), OPTIONAL :: TIPIV(:)
725 INTEGER(8), intent(in) :: LA
726 DOUBLE PRECISION, intent(inout) :: A(LA)
727 INTEGER, intent(in) :: NFRONT,NASS,N,LIW,INODE,LAST_ROW
728 INTEGER, intent(inout) :: IFLAG,INOPV,NOFFW, NBTINYW
729 INTEGER, intent(inout) :: DET_EXPW, DET_SIGNW
730 DOUBLE PRECISION, intent(inout) :: DET_MANTW
731 DOUBLE PRECISION, intent(in) :: UU, SEUIL
732 INTEGER, intent(inout) :: IW(LIW)
733 INTEGER, intent(in) :: IOLDPS
734 INTEGER(8), intent(in) :: POSELT
735 INTEGER KEEP(500)
736 INTEGER(8) KEEP8(150)
737 INTEGER, intent(in) :: LPN_LIST
738 INTEGER, intent(inout) :: PIVNUL_LIST(LPN_LIST)
739 DOUBLE PRECISION DKEEP(230)
740 INTEGER PP_FIRST2SWAP_L, PP_LastPanelonDisk_L,
741 & PP_LastPIVRPTRFilled_L,
742 & pp_first2swap_u, pp_lastpanelondisk_u,
743 & pp_lastpivrptrfilled_u
744 INTEGER, intent(in) :: PIVOT_OPTION, IEND_BLR
745 LOGICAL, intent(in) :: LR_ACTIVATED
746 INTEGER, intent(inout) :: Inextpiv
747 LOGICAL, intent(in) :: OOC_EFFECTIVE_ON_FRONT
748 INTEGER, intent(in) :: NVSCHUR
749 INTEGER, intent(in) :: PARPIV_T1
750 include 'mumps_headers.h'
751 DOUBLE PRECISION SWOP
752 INTEGER XSIZE
753 INTEGER(8) :: APOS, IDIAG
754 INTEGER(8) :: J1, J2, JJ, J3
755 INTEGER(8) :: NFRONT8
756 INTEGER ILOC
757 DOUBLE PRECISION ZERO
758 PARAMETER( ZERO = 0.0d0 )
759 DOUBLE PRECISION RZERO, RMAX, AMROW, MAX_PREV_in_PARPIV
760 INTEGER(8) :: APOSMAX, APOSROW
761 DOUBLE PRECISION :: RMAX_NORELAX
762 DOUBLE PRECISION PIVNUL, ABS_PIVOT
763 DOUBLE PRECISION FIXA, CSEUIL, PIVOT
764 INTEGER NPIV,IPIV, LRLOC
765 INTEGER NPIVP1,JMAX,J,ISW,ISWPS1
766 INTEGER ISWPS2,KSW, HF, IPIVNUL
767 INTEGER DMUMPS_IXAMAX
768 INTEGER :: ISHIFT, K206
769 INTEGER :: IPIV_SHIFT,IPIV_END
770 INTRINSIC max
771 DATA rzero /0.0d0/
772#if defined(_OPENMP)
773 INTEGER :: NOMP,CHUNK,K361
774#endif
775 INTEGER I_PIVRPTR_L, I_PIVR_L, NBPANELS_L
776 INTEGER I_PIVRPTR_U, I_PIVR_U, NBPANELS_U
777#if defined(_OPENMP)
778 nomp = omp_get_max_threads()
779 k361 = keep(361)
780#endif
781 pivnul = dkeep(1)
782 fixa = dkeep(2)
783 cseuil = seuil
784 nfront8 = int(nfront,8)
785 k206 = keep(206)
786 xsize = keep(ixsz)
787 npiv = iw(ioldps+1+xsize)
788 hf = 6 + iw(ioldps+5+xsize)+xsize
789 npivp1 = npiv + 1
790 aposmax = poselt+nfront8*nfront8-1_8
791 IF (ooc_effective_on_front) THEN
792 CALL dmumps_get_ooc_perm_ptr(typef_l, nbpanels_l,
793 & i_pivrptr_l, i_pivr_l,
794 & ioldps+2*nfront+6+iw(ioldps+5+xsize)+xsize,
795 & iw, liw)
796 CALL dmumps_get_ooc_perm_ptr(typef_u, nbpanels_u,
797 & i_pivrptr_u, i_pivr_u,
798 & ioldps+2*nfront+6+iw(ioldps+5+xsize)+xsize,
799 & iw, liw)
800 ENDIF
801 IF ( present(tipiv) ) THEN
802 iloc = npivp1 - ibeg_block + 1
803 tipiv(iloc) = iloc
804 ENDIF
805 IF (inopv .EQ. -1) THEN
806 apos = poselt + nfront8*int(npivp1-1,8) + int(npiv,8)
807 pivot = a(apos)
808 abs_pivot = abs(pivot)
809 idiag = apos
810 CALL dmumps_update_minmax_pivot
811 & ( abs_pivot, dkeep, keep, .true.)
812 IF(abs_pivot.LT.seuil) THEN
813 IF (dble(pivot) .GE. rzero) THEN
814 a(apos) = cseuil
815 ELSE
816 a(apos) = -cseuil
817 ENDIF
818 nbtinyw = nbtinyw + 1
819 ELSE IF (keep(258) .NE. 0) THEN
820 CALL dmumps_updatedeter(pivot, det_mantw, det_expw )
821 ENDIF
822 IF (ooc_effective_on_front) THEN
823 IF (keep(251).EQ.0) THEN
824 CALL dmumps_store_perminfo( iw(i_pivrptr_l),
825 & nbpanels_l,
826 & iw(i_pivr_l), nass, npivp1, npivp1,
827 & pp_lastpanelondisk_l,
828 & pp_lastpivrptrfilled_l)
829 ENDIF
830 CALL dmumps_store_perminfo( iw(i_pivrptr_u),
831 & nbpanels_u,
832 & iw(i_pivr_u), nass, npivp1, npivp1,
833 & pp_lastpanelondisk_u,
834 & pp_lastpivrptrfilled_u)
835 ENDIF
836 GO TO 420
837 ENDIF
838 inopv = 0
839 ishift = 0
840 ipiv_end = iend_block
841 IF (k206.GE.1) THEN
842 IF (inextpiv.GT.npivp1.AND.inextpiv.LE.iend_block) THEN
843 ishift = inextpiv - npivp1
844 ENDIF
845 IF ( k206.EQ.1
846 & .OR. (k206 .GT.1 .AND. iend_block.EQ.iend_blr) ) THEN
847 ipiv_end = iend_block + ishift
848 ENDIF
849 ENDIF
850 DO 460 ipiv_shift = npivp1+ishift, ipiv_end
851 IF (ipiv_shift .LE. iend_block) THEN
852 ipiv=ipiv_shift
853 ELSE
854 ipiv = ipiv_shift-iend_block-1+npivp1
855 IF (ibeg_block.EQ.npivp1) THEN
856 EXIT
857 ENDIF
858 ENDIF
859 apos = poselt + nfront8*int(ipiv-1,8) + int(npiv,8)
860 jmax = 1
861 IF ((pivot_option.EQ.0).OR.(uu.EQ.rzero)) THEN
862 abs_pivot = abs(a(apos))
863 IF(abs_pivot.LT.seuil) THEN
864 CALL dmumps_update_minmax_pivot
865 & ( abs_pivot, dkeep, keep, .true.)
866 IF (dble(a(apos)) .GE. rzero) THEN
867 a(apos) = cseuil
868 ELSE
869 a(apos) = -cseuil
870 ENDIF
871 nbtinyw = nbtinyw + 1
872 GO TO 420
873 ELSE IF (abs_pivot.EQ.rzero) THEN
874 GO TO 630
875 ENDIF
876 GO TO 380
877 ENDIF
878 amrow = rzero
879 j1 = apos
880 IF (pivot_option.EQ.1 .OR. (lr_activated .AND.
881 & (keep(480).GE.2
882 & ))) THEN
883 j = iend_blr - npiv
884 ELSE
885 j = nass - npiv
886 ENDIF
887 j2 = j1 + j - 1_8
888 jmax = dmumps_ixamax(j,a(j1),1,keep(361))
889 jj = j1 + int(jmax - 1,8)
890 amrow = abs(a(jj))
891 rmax = amrow
892 IF (pivot_option.GE.2) THEN
893 j1 = j2 + 1_8
894 IF (pivot_option.GE.3
895 & ) THEN
896 j2 = apos +
897 & int(- npiv + nfront - 1 - keep(253) - nvschur,8)
898 ELSE
899 j2 = apos +int(- npiv + nass - 1 ,8)
900 ENDIF
901 IF (j2.LT.j1) GO TO 370
902 IF (keep(351).EQ.1) THEN
903!$ CHUNK = max(K361/2,(int(J2-J1)+NOMP-1)/NOMP)
904!$OMP PARALLEL DO schedule(static, CHUNK) PRIVATE(JJ)
905!$OMP& FIRSTPRIVATE(J1,J2)
906!$OMP& REDUCTION(max:RMAX) IF ((J2-J1).GE.K361)
907 DO jj=j1,j2
908 rmax = max(abs(a(jj)),rmax)
909 ENDDO
910!$OMP END PARALLEL DO
911 ELSE
912 DO 360 jj=j1,j2
913 rmax = max(abs(a(jj)),rmax)
914 360 CONTINUE
915 ENDIF
916 370 CONTINUE
917 ENDIF
918 idiag = apos + int(ipiv - npivp1,8)
919 abs_pivot = abs(a(idiag))
920 IF (parpiv_t1.NE.0) THEN
921 rmax_norelax = dble(a(aposmax+int(ipiv,8)))
922 ELSE
923 rmax_norelax = rzero
924 ENDIF
925 rmax = max(rmax,rmax_norelax)
926 IF ( rmax .LE. pivnul ) THEN
927 IF (last_row.EQ.nfront) THEN
928 lrloc = last_row -keep(253)-nvschur
929 ELSE
930 lrloc = last_row
931 ENDIF
932 IF (nfront - keep(253) .EQ. nass) THEN
933 IF (iend_block.NE.nass ) THEN
934 GOTO 460
935 ENDIF
936 j1=poselt+int(ipiv-1,8)+int(npiv,8)*nfront8
937 j2=poselt+int(ipiv-1,8)+int(lrloc-1,8)*nfront8
938 ELSE
939 j1=poselt+int(ipiv-1,8)
940 j2=poselt+int(ipiv-1,8)+int(lrloc-1,8)*nfront8
941 ENDIF
942 DO jj=j1, j2, nfront8
943 IF ( abs(a(jj)) .GT. pivnul ) THEN
944 GOTO 460
945 END IF
946 ENDDO
947 IF ((parpiv_t1.NE.0)
948 & .AND.(parpiv_t1.NE.-1)
949 & .AND.(rmax_norelax.LT.0)
950 & .AND.(ipiv.GT.1)) THEN
951 max_prev_in_parpiv = rzero
952 DO jj=1,ipiv-1
953 max_prev_in_parpiv= max( max_prev_in_parpiv,
954 & dble(a(aposmax+int(jj,8))) )
955 ENDDO
956 IF (max_prev_in_parpiv.GT.pivnul) THEN
957 aposrow = poselt + nfront8*int(ipiv-1,8)
958 DO jj=1,ipiv-1
959 IF (abs(a(aposrow+jj-1)).GT.pivnul) GOTO 460
960 ENDDO
961 ENDIF
962 ENDIF
963 CALL dmumps_update_minmax_pivot
964 & ( abs_pivot, dkeep, keep, .true.)
965!$OMP ATOMIC CAPTURE
966 keep(109) = keep(109)+1
967 ipivnul = keep(109)
968!$OMP END ATOMIC
969 pivnul_list(ipivnul) = iw( ioldps+hf+npiv+ipiv-npivp1 )
970 IF(dble(fixa).GT.rzero) THEN
971 IF(dble(a(idiag)) .GE. rzero) THEN
972 a(idiag) = fixa
973 ELSE
974 a(idiag) = -fixa
975 ENDIF
976 ELSE
977 j1 = apos
978 j2 = apos +
979 & int(- npiv + nfront - 1 - keep(253),8)
980 DO jj=j1,j2
981 a(jj) = zero
982 ENDDO
983 a(idiag) = -fixa
984 ENDIF
985 jmax = ipiv - npiv
986 GOTO 385
987 ENDIF
988 rmax = max(rmax,abs(rmax_norelax))
989 IF (abs_pivot .GE. uu*rmax .AND.
990 & abs_pivot .GT. max(seuil,tiny(rmax))) THEN
991 jmax = ipiv - npiv
992 GO TO 380
993 ENDIF
994 IF ( .NOT. (amrow .GE. uu*rmax .AND.
995 & amrow .GT. max(seuil,tiny(rmax))) ) GO TO 460
996 noffw = noffw + 1
997 380 CONTINUE
998 IF (k206.GE.1) THEN
999 inextpiv = ipiv + 1
1000 ENDIF
1001 CALL dmumps_update_minmax_pivot
1002 & ( abs(a(apos+int(jmax-1,8))),
1003 & dkeep, keep, .false.)
1004 IF (keep(258) .NE. 0) THEN
1005 CALL dmumps_updatedeter( a(apos+int(jmax-1,8)),
1006 & det_mantw,
1007 & det_expw )
1008 ENDIF
1009 385 CONTINUE
1010 IF (ipiv.EQ.npivp1) GO TO 400
1011 IF (keep(405) .EQ. 0) THEN
1012 keep8(80) = keep8(80)+1
1013 ELSE
1014!$OMP ATOMIC UPDATE
1015 keep8(80) = keep8(80)+1
1016!$OMP END ATOMIC
1017 ENDIF
1018 IF (parpiv_t1.NE.0) THEN
1019 swop = a(aposmax+int(npivp1,8))
1020 a(aposmax+int(npivp1,8)) = a(aposmax+int(ipiv,8))
1021 a(aposmax+int(ipiv,8)) = swop
1022 ENDIF
1023 det_signw = - det_signw
1024 j1 = poselt + int(npiv,8)*nfront8
1025 j2 = j1 + nfront8 - 1_8
1026 j3 = poselt + int(ipiv-1,8)*nfront8
1027 DO 390 jj=j1,j2
1028 swop = a(jj)
1029 a(jj) = a(j3)
1030 a(j3) = swop
1031 j3 = j3 + 1_8
1032 390 CONTINUE
1033 iswps1 = ioldps + hf - 1 + npivp1
1034 iswps2 = ioldps + hf - 1 + ipiv
1035 isw = iw(iswps1)
1036 iw(iswps1) = iw(iswps2)
1037 iw(iswps2) = isw
1038 400 IF (jmax.EQ.1) GO TO 420
1039 det_signw = - det_signw
1040 IF ( present(tipiv) ) THEN
1041 tipiv(iloc) = iloc + jmax - 1
1042 ENDIF
1043 j1 = poselt + int(npiv,8)
1044 j2 = poselt + int(npiv + jmax - 1,8)
1045 DO 410 ksw=1,last_row
1046 swop = a(j1)
1047 a(j1) = a(j2)
1048 a(j2) = swop
1049 j1 = j1 + nfront8
1050 j2 = j2 + nfront8
1051 410 CONTINUE
1052 iswps1 = ioldps + hf - 1 + nfront + npiv + 1
1053 iswps2 = ioldps + hf - 1 + nfront + npiv + jmax
1054 isw = iw(iswps1)
1055 iw(iswps1) = iw(iswps2)
1056 iw(iswps2) = isw
1057 GO TO 420
1058 460 CONTINUE
1059 IF (k206 .GE. 1) THEN
1060 inextpiv=iend_block+1
1061 ENDIF
1062 IF (iend_block.EQ.nass) THEN
1063 inopv = 1
1064 ELSE
1065 inopv = 2
1066 ENDIF
1067 GO TO 430
1068 630 CONTINUE
1069 iflag = -10
1070 GOTO 430
1071 420 CONTINUE
1072 IF (ooc_effective_on_front) THEN
1073 IF (keep(251).EQ.0) THEN
1074 CALL dmumps_store_perminfo( iw(i_pivrptr_l),
1075 & nbpanels_l,
1076 & iw(i_pivr_l), nass, npivp1, ipiv,
1077 & pp_lastpanelondisk_l,
1078 & pp_lastpivrptrfilled_l)
1079 ENDIF
1080 CALL dmumps_store_perminfo( iw(i_pivrptr_u),
1081 & nbpanels_u,
1082 & iw(i_pivr_u), nass, npivp1, npiv+jmax,
1083 & pp_lastpanelondisk_u,
1084 & pp_lastpivrptrfilled_u)
1085 ENDIF
1086 430 CONTINUE
1087 RETURN
1088 END SUBROUTINE dmumps_fac_i
1089 SUBROUTINE dmumps_fac_i_ldlt
1090 & ( nfront,nass,inode,ibeg_block,iend_block,
1091 & iw,liw, a,la, inopv,
1092 & nnegw, nb22t1w, nbtinyw,
1093 & det_expw, det_mantw, det_signw,
1094 & iflag,ioldps,poselt,uu, seuil,keep,keep8,pivsiz,
1095 & dkeep,pivnul_list,lpn_list, xsize,
1096 & pp_first2swap_l, pp_lastpanelondisk,
1097 & pp_lastpivrptrindexfilled,maxfromm,is_maxfromm_avail,
1098 & pivot_option, iend_blr, inextpiv,
1099 & ooc_effective_on_front,
1100 & nvschur, parpiv_t1, lr_activated
1101 & )
1102!$ USE OMP_LIB
1103 USE mumps_ooc_common
1104 IMPLICIT NONE
1105 INTEGER(8) :: POSELT, LA
1106 INTEGER NFRONT,NASS,LIW,INODE,IFLAG,INOPV,
1107 & ioldps
1108 INTEGER, intent(inout) :: NNEGW, NB22T1W, NBTINYW
1109 INTEGER, intent(inout) :: DET_EXPW, DET_SIGNW
1110 DOUBLE PRECISION, intent(inout) :: DET_MANTW
1111 INTEGER, intent(in) :: IBEG_BLOCK, IEND_BLOCK
1112 INTEGER, intent(in) :: PIVOT_OPTION,IEND_BLR
1113 INTEGER, intent(inout) :: Inextpiv
1114 LOGICAL, intent(in) :: OOC_EFFECTIVE_ON_FRONT
1115 INTEGER PIVSIZ,LPIV, XSIZE
1116 DOUBLE PRECISION A(LA)
1117 DOUBLE PRECISION UU, UULOC, SEUIL
1118 INTEGER IW(LIW)
1119 INTEGER KEEP(500)
1120 INTEGER(8) KEEP8(150)
1121 INTEGER LPN_LIST
1122 INTEGER PIVNUL_LIST(LPN_LIST)
1123 DOUBLE PRECISION DKEEP(230)
1124 INTEGER PP_FIRST2SWAP_L, PP_LastPanelonDisk
1125 INTEGER PP_LastPIVRPTRIndexFilled
1126 DOUBLE PRECISION, intent(in) :: MAXFROMM
1127 LOGICAL, intent(inout) :: IS_MAXFROMM_AVAIL
1128 INTEGER, intent(in) :: NVSCHUR
1129 INTEGER, intent(in) :: PARPIV_T1
1130 LOGICAL, intent(in) :: LR_ACTIVATED
1131 include 'mpif.h'
1132 INTEGER (8) :: POSPV1,POSPV2,OFFDAG,APOSJ
1133 INTEGER JMAX, LIM, LIM_SWAP
1134 DOUBLE PRECISION RMAX,AMAX,TMAX, MAX_PREV_in_PARPIV, ABS_PIVOT
1135 DOUBLE PRECISION RMAX_NORELAX, TMAX_NORELAX, UULOCM1
1136 INTEGER(8) :: APOSMAX, APOSROW
1137 DOUBLE PRECISION MAXPIV
1138 DOUBLE PRECISION PIVNUL
1139 DOUBLE PRECISION MAXFROMM_UPDATED
1140 DOUBLE PRECISION FIXA, CSEUIL
1141 DOUBLE PRECISION PIVOT,DETPIV
1142 DOUBLE PRECISION ABSDETPIV
1143 INCLUDE 'mumps_headers.h'
1144 INTEGER :: HF, IPIVNUL
1145 INTEGER :: J
1146 INTEGER(8) :: APOS, J1, J2, JJ, NFRONT8, KK, J1_ini, JJ_ini
1147 INTEGER :: LDA
1148 INTEGER(8) :: LDA8
1149 INTEGER NPIV,IPIV
1150 INTEGER NPIVP1,K
1151 INTEGER :: ISHIFT, K206, IPIV_SHIFT, IPIV_END
1152 INTRINSIC max
1153 DOUBLE PRECISION ZERO, ONE
1154 parameter( zero = 0.0d0 )
1155 parameter( one = 1.0d0 )
1156 DOUBLE PRECISION RZERO,RONE
1157 parameter(rzero=0.0d0, rone=1.0d0)
1158#if defined(_OPENMP)
1159 LOGICAL :: OMP_FLAG
1160 INTEGER :: NOMP, CHUNK, J1_end
1161#endif
1162 INTEGER I_PIVRPTR, I_PIVR, NBPANELS_L
1163!$ NOMP = OMP_GET_MAX_THREADS()
1164 pivnul = dkeep(1)
1165 fixa = dkeep(2)
1166 cseuil = seuil
1167 lda = nfront
1168 lda8 = int(lda,8)
1169 nfront8 = int(nfront,8)
1170 k206 = keep(206)
1171 uuloc = uu
1172 IF (uuloc.GT.rzero) THEN
1173 uulocm1 = rone/uuloc
1174 ELSE
1175 uulocm1 = rone
1176 ENDIF
1177 hf = 6 + xsize
1178 IF (keep(50).NE.1 .AND. ooc_effective_on_front) THEN
1179 CALL dmumps_get_ooc_perm_ptr(typef_l, nbpanels_l,
1180 & i_pivrptr, i_pivr, ioldps+2*nfront+6+keep(ixsz),
1181 & iw, liw)
1182 ENDIF
1183 pivsiz = 1
1184 npiv = iw(ioldps+1+xsize)
1185 npivp1 = npiv + 1
1186 aposmax = poselt+lda8*lda8-1_8
1187 IF(inopv .EQ. -1) THEN
1188 apos = poselt + (lda8+1_8) * int(npiv,8)
1189 pospv1 = apos
1190 CALL dmumps_update_minmax_pivot
1191 & ( abs(a(apos)), dkeep, keep, .true.)
1192 IF(abs(a(apos)).LT.seuil) THEN
1193 IF(dble(a(apos)) .GE. rzero) THEN
1194 a(apos) = cseuil
1195 ELSE
1196 a(apos) = -cseuil
1197 nnegw = nnegw+1
1198 ENDIF
1199 nbtinyw = nbtinyw + 1
1200 ELSE
1201 IF (keep(258) .NE. 0) THEN
1202 CALL dmumps_updatedeter( a(apos), det_mantw, det_expw )
1203 ENDIF
1204 ENDIF
1205 IF (keep(50).NE.1 .AND. ooc_effective_on_front) THEN
1206 CALL dmumps_store_perminfo( iw(i_pivrptr), nbpanels_l,
1207 & iw(i_pivr), nass, npivp1, npivp1,
1208 & pp_lastpanelondisk,
1209 & pp_lastpivrptrindexfilled)
1210 ENDIF
1211 GO TO 420
1212 ENDIF
1213 inopv = 0
1214 ishift = 0
1215 ipiv_end = iend_block
1216 IF (k206.GE.1) THEN
1217 IF (inextpiv.GT.npivp1.AND.inextpiv.LE.iend_block) THEN
1218 ishift = inextpiv - npivp1
1219 ENDIF
1220 IF ( k206.EQ.1
1221 & .OR. (k206 .GT.1 .AND. iend_block.EQ.iend_blr) ) THEN
1222 ipiv_end = iend_block + ishift
1223 ENDIF
1224 IF (ishift.GT.0.AND.is_maxfromm_avail) THEN
1225 ipiv = npivp1
1226 apos = poselt + lda8*int(ipiv-1,8) + int(npiv,8)
1227 pospv1 = apos + int(ipiv - npivp1,8)
1228 pivot = a(pospv1)
1229 IF ( maxfromm .GT. pivnul ) THEN
1230 IF (parpiv_t1.NE.0) THEN
1231 maxfromm_updated = max
1232 & ( maxfromm,
1233 & abs(dble(a(aposmax+int(ipiv,8))))
1234 & )
1235 ELSE
1236 maxfromm_updated = maxfromm
1237 ENDIF
1238 IF ( (abs(pivot) .GE. uuloc*maxfromm_updated).AND.
1239 & abs(pivot) .GT. max(seuil,tiny(maxfromm_updated))
1240 & ) THEN
1241 ishift = 0
1242 ENDIF
1243 ENDIF
1244 ENDIF
1245 IF ( ishift .GT. 0) THEN
1246 is_maxfromm_avail = .false.
1247 ENDIF
1248 ENDIF
1249 DO 460 ipiv_shift = npivp1+ishift, ipiv_end
1250 IF (ipiv_shift .LE. iend_block) THEN
1251 ipiv=ipiv_shift
1252 ELSE
1253 ipiv = ipiv_shift-iend_block-1+npivp1
1254 IF (ibeg_block.EQ.npivp1) THEN
1255 EXIT
1256 ENDIF
1257 ENDIF
1258 apos = poselt + lda8*int(ipiv-1,8) + int(npiv,8)
1259 pospv1 = apos + int(ipiv - npivp1,8)
1260 pivot = a(pospv1)
1261 abs_pivot = abs(pivot)
1262 IF (uuloc.EQ.rzero.OR.pivot_option.EQ.0) THEN
1263 IF(abs_pivot.LT.seuil) THEN
1264 CALL dmumps_update_minmax_pivot
1265 & ( abs_pivot, dkeep, keep, .true.)
1266 IF(dble(pivot) .GE. rzero) THEN
1267 a(pospv1) = cseuil
1268 ELSE
1269 a(pospv1) = -cseuil
1270 nnegw = nnegw+1
1271 ENDIF
1272 nbtinyw = nbtinyw + 1
1273 ELSE IF (abs_pivot.EQ.rzero) THEN
1274 GO TO 630
1275 ELSE
1276 IF (pivot.LT.rzero) nnegw = nnegw+1
1277 CALL dmumps_update_minmax_pivot
1278 & ( abs_pivot, dkeep, keep, .false.)
1279 IF (keep(258) .NE. 0) THEN
1280 CALL dmumps_updatedeter(pivot, det_mantw, det_expw )
1281 ENDIF
1282 ENDIF
1283 GO TO 420
1284 ENDIF
1285 IF ( is_maxfromm_avail ) THEN
1286 IF ( maxfromm .GT. pivnul ) THEN
1287 IF (parpiv_t1.NE.0) THEN
1288 maxfromm_updated = max
1289 & ( maxfromm,
1290 & abs(dble(a(aposmax+int(ipiv,8))))
1291 & )
1292 ELSE
1293 maxfromm_updated = maxfromm
1294 ENDIF
1295 IF ( (abs_pivot .GE. uuloc*maxfromm_updated).AND.
1296 & (abs_pivot .GT. max(seuil,tiny(maxfromm_updated)))
1297 & ) THEN
1298 IF (pivot .LT. rzero) nnegw = nnegw+1
1299 CALL dmumps_update_minmax_pivot
1300 & ( abs_pivot,
1301 & dkeep, keep, .false.)
1302 IF (keep(258) .NE. 0) THEN
1303 CALL dmumps_updatedeter(pivot, det_mantw, det_expw )
1304 ENDIF
1305 GOTO 415
1306 ENDIF
1307 ENDIF
1308 is_maxfromm_avail = .false.
1309 ENDIF
1310 amax = -rone
1311 jmax = 0
1312 IF (pivot_option.EQ.3
1313 & ) THEN
1314 lim = nfront - keep(253)-nvschur
1315 ELSEIF (pivot_option.GE.2
1316 & ) THEN
1317 lim = nass
1318 ELSEIF (pivot_option.GE.1) THEN
1319 lim = iend_blr
1320 ELSE
1321 write(*,*) 'Internal error in FAC_I_LDLT 1x1:',
1322 & pivot_option
1323 CALL mumps_abort()
1324 ENDIF
1325 j1 = apos
1326 j2 = pospv1 - 1_8
1327 DO jj=j1,j2
1328 IF(abs(a(jj)) .GT. amax) THEN
1329 amax = abs(a(jj))
1330 jmax = ipiv - int(pospv1-jj)
1331 ENDIF
1332 ENDDO
1333 j1 = pospv1 + lda8
1334 DO j=1, iend_block - ipiv
1335 IF(abs(a(j1)) .GT. amax) THEN
1336 amax = abs(a(j1))
1337 jmax = ipiv + j
1338 ENDIF
1339 j1 = j1 + lda8
1340 ENDDO
1341 rmax = rzero
1342 j1_ini = j1
1343#if defined(_openmp)
1344 j1_end = lim - iend_block
1345 chunk = max(j1_end,1)
1346 IF ( j1_end.GE.keep(360)) THEN
1347 omp_flag = .true.
1348 chunk = max(keep(360)/2,(j1_end+nomp-1)/nomp)
1349 ELSE
1350 omp_flag = .false.
1351 ENDIF
1352#endif
1353!$OMP PARALLEL DO SCHEDULE(static, CHUNK) PRIVATE(J1)
1354!$OMP& REDUCTION(max:RMAX) IF(OMP_FLAG)
1355 DO j=1, lim - iend_block
1356 j1 = j1_ini + int(j-1,8) * lda8
1357 rmax = max(abs(a(j1)),rmax)
1358 ENDDO
1359!$OMP END PARALLEL DO
1360 IF (parpiv_t1.NE.0) THEN
1361 rmax_norelax = dble(a(aposmax+int(ipiv,8)))
1362 ELSE
1363 rmax_norelax = rzero
1364 ENDIF
1365 rmax = max(rmax,rmax_norelax)
1366 IF (max(amax,rmax,abs(pivot)).LE.pivnul) THEN
1367 IF ((parpiv_t1.NE.0)
1368 & .AND.(parpiv_t1.NE.-1)
1369 & .AND.(rmax_norelax.LT.0)
1370 & .AND.(ipiv.GT.1)) THEN
1371 max_prev_in_parpiv = rzero
1372 DO jj=1,ipiv-1
1373 max_prev_in_parpiv= max( max_prev_in_parpiv,
1374 & dble(a(aposmax+int(jj,8))) )
1375 ENDDO
1376 IF (max_prev_in_parpiv.GT.pivnul) THEN
1377 aposrow = poselt + nfront8*int(ipiv-1,8)
1378 DO jj=1,ipiv-1
1379 IF (abs(a(aposrow+jj-1)).GT.pivnul) THEN
1380 GOTO 460
1381 ENDIF
1382 ENDDO
1383 ENDIF
1384 ENDIF
1385 CALL dmumps_update_minmax_pivot
1386 & ( abs(a(pospv1)), dkeep, keep, .true.)
1387!$OMP ATOMIC CAPTURE
1388 keep(109) = keep(109) + 1
1389 ipivnul = keep(109)
1390!$OMP END ATOMIC
1391 pivnul_list(ipivnul) = iw( ioldps+hf+npiv+ipiv-npivp1 )
1392 IF(dble(fixa).GT.rzero) THEN
1393 IF(dble(pivot) .GE. rzero) THEN
1394 a(pospv1) = fixa
1395 ELSE
1396 a(pospv1) = -fixa
1397 ENDIF
1398 ELSE
1399 j1 = apos
1400 j2 = pospv1 - 1_8
1401 DO jj=j1,j2
1402 a(jj) = zero
1403 ENDDO
1404 j1 = pospv1 + lda8
1405 DO j=1, iend_block - ipiv
1406 a(j1) = zero
1407 j1 = j1 + lda8
1408 ENDDO
1409 DO j=1,lim - iend_block
1410 a(j1) = zero
1411 j1 = j1 + lda8
1412 ENDDO
1413 a(pospv1) = one
1414 ENDIF
1415 pivot = a(pospv1)
1416 GO TO 415
1417 ENDIF
1418 rmax = max(rmax,abs(rmax_norelax))
1419 IF ( abs(pivot).GE.uuloc*max(rmax,amax)
1420 & .AND. abs(pivot) .GT. max(seuil,tiny(rmax)) ) THEN
1421 IF (pivot .LT. zero) nnegw = nnegw+1
1422 CALL dmumps_update_minmax_pivot
1423 & ( abs(pivot),
1424 & dkeep, keep, .false.)
1425 IF (keep(258) .NE.0 ) THEN
1426 CALL dmumps_updatedeter(pivot, det_mantw, det_expw )
1427 ENDIF
1428 GO TO 415
1429 END IF
1430 IF (npivp1.EQ.iend_block) THEN
1431 GOTO 460
1432 ELSE IF (jmax.EQ.0) THEN
1433 GOTO 460
1434 ENDIF
1435 IF (max(abs(pivot),rmax,amax).LE.tiny(rmax)) THEN
1436 GOTO 460
1437 ENDIF
1438 IF (
1439 & (keep(19).NE.0).AND.(max(amax,rmax,abs(pivot)).LE.seuil)
1440 & )
1441 & THEN
1442 GO TO 460
1443 ENDIF
1444 IF (rmax.LT.amax) THEN
1445 j1 = apos
1446 j2 = pospv1 - 1_8
1447 DO jj=j1,j2
1448 IF(int(pospv1-jj) .NE. ipiv-jmax) THEN
1449 rmax = max(rmax,abs(a(jj)))
1450 ENDIF
1451 ENDDO
1452 j1 = pospv1 + lda8
1453 DO j=1,iend_block-ipiv
1454 IF(ipiv+j .NE. jmax) THEN
1455 rmax = max(abs(a(j1)),rmax)
1456 ENDIF
1457 j1 = j1 + lda8
1458 ENDDO
1459 ENDIF
1460 aposj = poselt + int(jmax-1,8)*lda8 + int(npiv,8)
1461 pospv2 = aposj + int(jmax - npivp1,8)
1462 IF (ipiv.LT.jmax) THEN
1463 offdag = aposj + int(ipiv - npivp1,8)
1464 ELSE
1465 offdag = apos + int(jmax - npivp1,8)
1466 END IF
1467 tmax = rzero
1468#if defined(_OPENMP)
1469 j1_end = lim-jmax
1470 chunk = max(j1_end,1)
1471 IF (j1_end.GE.keep(360)) THEN
1472 omp_flag = .true.
1473 chunk = max(keep(360)/2,(j1_end+nomp-1)/nomp)
1474 ELSE
1475 omp_flag = .false.
1476 ENDIF
1477#endif
1478 IF (jmax .LT. ipiv) THEN
1479 jj_ini = pospv2
1480!$OMP PARALLEL DO SCHEDULE(static, CHUNK) IF (OMP_FLAG)
1481!$OMP& PRIVATE(JJ) REDUCTION(max:TMAX)
1482 DO k = 1, lim - jmax
1483 jj = jj_ini+ int(k,8)*nfront8
1484 IF (jmax+k.NE.ipiv) THEN
1485 tmax=max(tmax,abs(a(jj)))
1486 ENDIF
1487 ENDDO
1488!$OMP END PARALLEL DO
1489 DO kk = aposj, pospv2-1_8
1490 tmax = max(tmax,abs(a(kk)))
1491 ENDDO
1492 ELSE
1493 jj_ini = pospv2
1494!$OMP PARALLEL DO SCHEDULE(static, CHUNK) PRIVATE(JJ)
1495!$OMP& REDUCTION(max:TMAX) IF(OMP_FLAG)
1496 DO k = 1, lim-jmax
1497 jj = jj_ini + int(k,8)*nfront8
1498 tmax=max(tmax,abs(a(jj)))
1499 ENDDO
1500!$OMP END PARALLEL DO
1501 DO kk = aposj, pospv2 - 1_8
1502 IF (kk.NE.offdag) THEN
1503 tmax = max(tmax,abs(a(kk)))
1504 ENDIF
1505 ENDDO
1506 ENDIF
1507 IF (parpiv_t1.NE.0) THEN
1508 tmax_norelax = max(seuil*uulocm1,
1509 & abs(dble(a(aposmax+int(jmax,8))))
1510 & )
1511 ELSE
1512 tmax_norelax = seuil*uulocm1
1513 ENDIF
1514 tmax = max(tmax,tmax_norelax)
1515 detpiv = a(pospv1)*a(pospv2) - a(offdag)**2
1516 absdetpiv = abs(detpiv)
1517 IF (seuil.GT.rzero) THEN
1518 IF (sqrt(absdetpiv) .LE. seuil ) THEN
1519 GOTO 460
1520 ENDIF
1521 ENDIF
1522 maxpiv = max(abs(a(pospv1)),abs(a(pospv2)))
1523 IF (maxpiv.EQ.rzero) maxpiv = rone
1524 IF ((abs(a(pospv2))*rmax+amax*tmax)*uuloc.GT.
1525 & absdetpiv .OR. (absdetpiv .EQ. rzero) ) THEN
1526 GO TO 460
1527 ENDIF
1528 IF ((abs(a(pospv1))*tmax+amax*rmax)*uuloc.GT.
1529 & absdetpiv .OR. (absdetpiv.EQ. rzero) ) THEN
1530 GO TO 460
1531 ENDIF
1532 CALL dmumps_update_minmax_pivot
1533 & ( sqrt(absdetpiv),
1534 & dkeep, keep, .false.)
1535 IF (keep(258) .NE.0 ) THEN
1536 CALL dmumps_updatedeter(detpiv, det_mantw, det_expw )
1537 ENDIF
1538 pivsiz = 2
1539 nb22t1w = nb22t1w + 1
1540 IF(detpiv .LT. rzero) THEN
1541 nnegw = nnegw+1
1542 ELSE IF(a(pospv2) .LT. rzero) THEN
1543 nnegw = nnegw+2
1544 ENDIF
1545 415 CONTINUE
1546 IF (k206.GE.1) THEN
1547 inextpiv = max(npivp1+pivsiz, ipiv+1)
1548 ENDIF
1549 DO k=1,pivsiz
1550 IF (pivsiz .EQ. 2) THEN
1551 IF (k==1) THEN
1552 lpiv = min(ipiv,jmax)
1553 ELSE
1554 lpiv = max(ipiv,jmax)
1555 ENDIF
1556 ELSE
1557 lpiv = ipiv
1558 ENDIF
1559 IF (lpiv.EQ.npivp1) GOTO 416
1560 IF (keep(405) .EQ. 0) THEN
1561 keep8(80) = keep8(80)+1
1562 ELSE
1563!$OMP ATOMIC UPDATE
1564 keep8(80) = keep8(80)+1
1565!$OMP END ATOMIC
1566 ENDIF
1567 lim_swap = nfront
1568 CALL dmumps_swap_ldlt( a, la, iw, liw,
1569 & ioldps, npivp1, lpiv, poselt, lim_swap,
1570 & lda, nfront, 1, parpiv_t1, keep(50),
1571 & keep(ixsz), -9999)
1572 416 CONTINUE
1573 IF (keep(50).NE.1 .AND. ooc_effective_on_front) THEN
1574 CALL dmumps_store_perminfo( iw(i_pivrptr), nbpanels_l,
1575 & iw(i_pivr), nass, npivp1, lpiv, pp_lastpanelondisk,
1576 & pp_lastpivrptrindexfilled)
1577 ENDIF
1578 npivp1 = npivp1 + 1
1579 ENDDO
1580 IF(pivsiz .EQ. 2) THEN
1581 a(poselt+(lda8+1_8)*int(npiv,8)+1_8) = detpiv
1582 ENDIF
1583 GOTO 420
1584 460 CONTINUE
1585 IF (k206 .GE. 1) THEN
1586 inextpiv=iend_block+1
1587 ENDIF
1588 IF (iend_block.EQ.nass) THEN
1589 inopv = 1
1590 ELSE
1591 inopv = 2
1592 ENDIF
1593 GO TO 420
1594 630 CONTINUE
1595 pivsiz = 0
1596 iflag = -10
1597 420 CONTINUE
1598 is_maxfromm_avail = .false.
1599 RETURN
1600 END SUBROUTINE dmumps_fac_i_ldlt
1601 SUBROUTINE dmumps_fac_mq_ldlt(IEND_BLOCK,
1602 & NFRONT,NASS,NPIV,INODE,
1603 & A,LA,LDA,
1604 & POSELT,IFINB,PIVSIZ,
1605 & MAXFROMM, IS_MAXFROMM_AVAIL, IS_MAX_USEFUL,
1606 & PARPIV_T1, LAST_ROW, IEND_BLR, NVSCHUR_K253,
1607 & LR_ACTIVATED
1608 & )
1609 IMPLICIT NONE
1610 INTEGER, intent(out):: IFINB
1611 INTEGER, intent(in) :: INODE, NFRONT, NASS, NPIV
1612 INTEGER, intent(in) :: IEND_BLOCK
1613 INTEGER, intent(in) :: LDA
1614 INTEGER(8), intent(in) :: LA
1615 DOUBLE PRECISION, intent(inout) :: A(LA)
1616 INTEGER, intent(in) :: LAST_ROW
1617 INTEGER, intent(in) :: IEND_BLR
1618 INTEGER(8) :: POSELT
1619 DOUBLE PRECISION, intent(out) :: MAXFROMM
1620 LOGICAL, intent(out) :: IS_MAXFROMM_AVAIL
1621 LOGICAL, intent(in) :: IS_MAX_USEFUL
1622 INTEGER, intent(in) :: PARPIV_T1
1623 INTEGER, INTENT(in) :: NVSCHUR_K253
1624 LOGICAL, intent(in) :: LR_ACTIVATED
1625 DOUBLE PRECISION VALPIV
1626 DOUBLE PRECISION :: MAXFROMMTMP
1627 INTEGER NCB1
1628 INTEGER(8) :: NFRONT8
1629 INTEGER(8) :: LDA8
1630 INTEGER(8) :: K1POS
1631 INTEGER NEL2
1632 DOUBLE PRECISION ONE, ZERO
1633 DOUBLE PRECISION A11,A22,A12
1634 INTEGER(8) :: APOS, LPOS, LPOS1, LPOS2
1635 INTEGER(8) :: POSPV1, POSPV2
1636 INTEGER :: PIVSIZ,NPIV_NEW,J2,I
1637 INTEGER(8) :: OFFDAG, OFFDAG_OLD, K1, K2, IROW
1638#if defined(__ve__)
1639 INTEGER(8) :: J2_8, KU1, KU2
1640#else
1641 INTEGER(8) :: IBEG, IEND, JJ_LOC, JJ, ROW_SHIFT
1642 INTEGER(8) :: IBEG_LOC, IEND_LOC
1643#endif
1644 DOUBLE PRECISION SWOP,DETPIV,MULT1,MULT2
1645 INTEGER(8) :: APOSMAX
1646!$ LOGICAL :: OMP_FLAG
1647 include 'mumps_headers.h'
1648 parameter(one = 1.0d0,
1649 & zero = 0.0d0)
1650 lda8 = int(lda,8)
1651 nfront8 = int(nfront,8)
1652 npiv_new = npiv + pivsiz
1653 ifinb = 0
1654 is_maxfromm_avail = .false.
1655 ncb1 = last_row - iend_block
1656 nel2 = iend_block - npiv_new
1657 IF (nel2.EQ.0) THEN
1658 IF (iend_block.EQ.nass) THEN
1659 ifinb = -1
1660 ELSE
1661 ifinb = 1
1662 ENDIF
1663 ENDIF
1664 maxfromm = 0.0d0
1665 IF(pivsiz .EQ. 1) THEN
1666 apos = poselt + int(npiv,8)*(nfront8 + 1_8)
1667 valpiv = one/a(apos)
1668 lpos = apos + lda8
1669#if defined(__ve__)
1670 IF (nel2+ncb1.GT.0) THEN
1671!$ OMP_FLAG = (NCB1 + NEL2> 300)
1672!$OMP PARALLEL DO PRIVATE(I,K1POS) IF (OMP_FLAG)
1673!nec$ ivdep
1674 DO i=1, nel2 + ncb1
1675 k1pos = lpos+ int(i-1,8)*lda8
1676 a(apos+int(i,8))=a(k1pos)
1677 ENDDO
1678!$OMP END PARALLEL DO
1679!$OMP PARALLEL DO PRIVATE(I,K1POS) IF (OMP_FLAG)
1680!NEC$ IVDEP
1681 DO i=1, nel2 + ncb1
1682 k1pos = lpos+ int(i-1,8)*lda8
1683 a(k1pos) = a(k1pos) * valpiv
1684 ENDDO
1685!$OMP END PARALLEL DO
1686 IF (.NOT. is_max_useful) THEN
1687!$ OMP_FLAG = (NCB1 > 300).AND.(NEL2.GE.2)
1688!$OMP PARALLEL DO PRIVATE(I,J2,J2_8,K1POS) IF (OMP_FLAG)
1689!NEC$ IVDEP
1690 DO j2 = 1, nel2
1691 j2_8 = int(j2,8)
1692!nec$ ivdep
1693 DO i=j2, nel2 + ncb1
1694 k1pos = lpos+ int(i-1,8)*lda8
1695 a(k1pos+j2_8)=a(k1pos+j2_8)-(a(k1pos)*a(apos+j2_8))
1696 ENDDO
1697 ENDDO
1698!$OMP END PARALLEL DO
1699 ELSE
1700 IF (nel2.GT.0) THEN
1701 maxfrommtmp=0.0d0
1702!$ OMP_FLAG = (NCB1+NEL2 > 300)
1703!$OMP PARALLEL DO PRIVATE(I,K1POS) IF (OMP_FLAG)
1704!$OMP& REDUCTION(max:MAXFROMMTMP)
1705!nec$ ivdep
1706 DO i=1, nel2 + ncb1 - nvschur_k253
1707 k1pos = lpos+ int(i-1,8)*lda8
1708 a(k1pos+1_8)=a(k1pos+1_8) - a(k1pos) * a(apos+1_8)
1709 maxfrommtmp=max(maxfrommtmp, abs(a(k1pos+1_8)))
1710 ENDDO
1711 !$OMP END PARALLEL DO
1712 is_maxfromm_avail = .true.
1713 maxfromm=max(maxfromm, maxfrommtmp)
1714 IF (nvschur_k253.GT.0) THEN
1715 DO i= nel2 + ncb1- nvschur_k253 +1, nel2 + ncb1
1716 k1pos = lpos+ int(i-1,8)*lda8
1717 a(k1pos+1_8)=a(k1pos+1_8) - a(k1pos) * a(apos+1_8)
1718 ENDDO
1719 ENDIF
1720 ENDIF
1721 IF (nel2.GT.1) THEN
1722!$ OMP_FLAG = (NCB1+NEL2 > 300).AND.(NEL2.GE.3)
1723!$OMP PARALLEL DO PRIVATE(I,J2,J2_8,K1POS) IF (OMP_FLAG)
1724!NEC$ IVDEP
1725 DO j2 = 2, nel2
1726 j2_8 = int(j2,8)
1727!NEC$ IVDEP
1728 DO i=j2, nel2 + ncb1
1729 k1pos = lpos+ int(i-1,8)*lda8
1730 a(k1pos+j2_8)=a(k1pos+j2_8)-(a(k1pos)*a(apos+j2_8))
1731 ENDDO
1732 ENDDO
1733!$OMP END PARALLEL DO
1734 ENDIF
1735 ENDIF
1736 ENDIF
1737#else
1738 IF (nel2 > 0) THEN
1739 IF (.NOT. is_max_useful) THEN
1740 DO i=1, nel2
1741 k1pos = lpos + int(i-1,8)*lda8
1742 a(apos+int(i,8))=a(k1pos)
1743 a(k1pos) = a(k1pos) * valpiv
1744 DO jj=1_8, int(i,8)
1745 a(k1pos+jj)=a(k1pos+jj) - a(k1pos) * a(apos+jj)
1746 ENDDO
1747 ENDDO
1748 ELSE
1749 is_maxfromm_avail = .true.
1750 DO i=1, nel2
1751 k1pos = lpos + int(i-1,8)*lda8
1752 a(apos+int(i,8))=a(k1pos)
1753 a(k1pos) = a(k1pos) * valpiv
1754 a(k1pos+1_8)=a(k1pos+1_8) - a(k1pos) * a(apos+1_8)
1755 maxfromm=max( maxfromm,abs(a(k1pos+1_8)) )
1756 DO jj = 2_8, int(i,8)
1757 a(k1pos+jj)=a(k1pos+jj) - a(k1pos) * a(apos+jj)
1758 ENDDO
1759 ENDDO
1760 ENDIF
1761 ENDIF
1762 IF (ncb1.GT.0) THEN
1763 IF (.NOT. is_max_useful) THEN
1764!$OMP PARALLEL DO PRIVATE(JJ,K1POS) IF (NCB1 > 300)
1765 DO i=nel2+1, nel2 + ncb1
1766 k1pos = lpos+ int(i-1,8)*lda8
1767 a(apos+int(i,8))=a(k1pos)
1768 a(k1pos) = a(k1pos) * valpiv
1769 DO jj = 1_8, int(nel2,8)
1770 a(k1pos+jj)=a(k1pos+jj) - a(k1pos) * a(apos+jj)
1771 ENDDO
1772 ENDDO
1773!$OMP END PARALLEL DO
1774 ELSE
1775 maxfrommtmp=0.0d0
1776!$ OMP_FLAG = (NCB1-NVSCHUR_K253>300)
1777!$OMP PARALLEL DO PRIVATE(JJ,K1POS)
1778!$OMP& REDUCTION(max:MAXFROMMTMP) IF (OMP_FLAG)
1779 DO i=nel2+1, nel2 + ncb1 - nvschur_k253
1780 k1pos = lpos+ int(i-1,8)*lda8
1781 a(apos+int(i,8))=a(k1pos)
1782 a(k1pos) = a(k1pos) * valpiv
1783 IF (nel2 > 0) THEN
1784 a(k1pos+1_8) = a(k1pos+1_8) - a(k1pos) * a(apos+1_8)
1785 maxfrommtmp=max(maxfrommtmp, abs(a(k1pos+1_8)))
1786 DO jj = 2_8, int(nel2,8)
1787 a(k1pos+jj)=a(k1pos+jj) - a(k1pos) * a(apos+jj)
1788 ENDDO
1789 ENDIF
1790 ENDDO
1791!$OMP END PARALLEL DO
1792 DO i = nel2 + ncb1 - nvschur_k253 + 1, nel2 + ncb1
1793 k1pos = lpos+ int(i-1,8)*lda8
1794 a(apos+int(i,8))=a(k1pos)
1795 a(k1pos) = a(k1pos) * valpiv
1796 DO jj = 1_8, int(nel2,8)
1797 a(k1pos+jj)=a(k1pos+jj) - a(k1pos) * a(apos+jj)
1798 ENDDO
1799 ENDDO
1800 maxfromm=max(maxfromm, maxfrommtmp)
1801 ENDIF
1802 ENDIF
1803#endif
1804 ELSE
1805 pospv1 = poselt + int(npiv,8)*(nfront8 + 1_8)
1806 pospv2 = pospv1 + nfront8 + 1_8
1807 offdag_old = pospv2 - 1_8
1808 offdag = pospv1 + 1_8
1809 swop = a(pospv2)
1810 detpiv = a(offdag)
1811 a22 = a(pospv1)/detpiv
1812 a11 = swop/detpiv
1813 a12 = -a(offdag_old)/detpiv
1814 a(offdag) = a(offdag_old)
1815 a(offdag_old) = zero
1816 lpos1 = pospv2 + lda8 - 1_8
1817 lpos2 = lpos1 + 1_8
1818#if defined(__ve__)
1819#if defined(check)
1820 IF (lda8.NE.nfront8) CALL mumps_abort()
1821 IF (nel2 + ncb1.NE.last_row-npiv_new) CALL mumps_abort()
1822#endif
1823 CALL dcopy(last_row-npiv_new, a(lpos1), lda, a(pospv1+2_8), 1)
1824 CALL dcopy(last_row-npiv_new, a(lpos2), lda, a(pospv2+1_8), 1)
1825!$ OMP_FLAG = (NEL2+NCB1 > 300)
1826!$OMP PARALLEL DO PRIVATE(J2,J2_8,I,K1,K2,KU1,KU2)
1827!$OMP& IF (OMP_FLAG)
1828!NEC$ IVDEP
1829 DO j2=1, nel2 + ncb1
1830 j2_8 = int(j2,8)
1831 ku1 = pospv1 + 2_8 + (j2_8-1_8)
1832 ku2 = pospv2 + 1_8 + (j2_8-1_8)
1833 k1 = lpos1 + (j2_8-1_8)*nfront8
1834 k2 = k1 + 1_8
1835 a(k1) = a11*a(ku1)+a12*a(ku2)
1836 a(k2) = a12*a(ku1)+a22*a(ku2)
1837 ENDDO
1838 IF (nel2.GT.0) THEN
1839!$ OMP_FLAG = (NCB1+NEL2 > 300).AND.(NEL2.GE.2)
1840!$OMP PARALLEL DO PRIVATE(I,J2,J2_8,K1,K2,MULT1,MULT2,IROW)
1841!$OMP& IF (OMP_FLAG)
1842!NEC$ IVDEP
1843 DO j2 = 1,nel2
1844 j2_8 = int(j2,8)
1845 mult1 = -a(pospv1 + 2_8 + j2_8-1_8)
1846 mult2 = -a(pospv2 + 1_8 + j2_8-1_8)
1847!NEC$ IVDEP
1848 DO i= j2, nel2 + ncb1
1849 k1 = lpos1 + (int(i,8)-1_8)*nfront8
1850 k2 = k1 + 1_8
1851 irow = k2 + j2_8
1852 a(irow) = a(irow) + mult1*a(k1) +
1853 & mult2*a(k2)
1854 ENDDO
1855 ENDDO
1856 ENDIF
1857#else
1858 jj = pospv2 + nfront8-1_8
1859 ibeg = jj + 2_8
1860 iend = ibeg
1861 DO j2 = 1,nel2
1862 k1 = jj
1863 k2 = jj+1_8
1864 mult1 = - (a11*a(k1)+a12*a(k2))
1865 mult2 = - (a12*a(k1)+a22*a(k2))
1866 a(pospv1 + 2_8 + (int(j2,8)-1_8)) = a(k1)
1867 a(pospv2 + 1_8 + (int(j2,8)-1_8)) = a(k2)
1868 k1 = pospv1 + 2_8
1869 k2 = pospv2 + 1_8
1870 DO irow = ibeg, iend
1871 a(irow) = a(irow) + mult1*a(k1) + mult2*a(k2)
1872 k1 = k1 + 1_8
1873 k2 = k2 + 1_8
1874 ENDDO
1875 a( jj ) = -mult1
1876 a( jj + 1_8 ) = -mult2
1877 ibeg = ibeg + nfront8
1878 iend = iend + nfront8 + 1_8
1879 jj = jj+nfront8
1880 ENDDO
1881 iend = iend-1_8
1882!$OMP PARALLEL DO PRIVATE(J2, K1, K2, MULT1, MULT2, IROW, JJ_LOC,
1883!$OMP& ROW_SHIFT, IBEG_LOC, IEND_LOC) IF (LAST_ROW-IEND_BLOCK>300)
1884 DO j2 = 1,last_row-iend_block
1885 row_shift = (j2-1_8)*nfront8
1886 jj_loc = jj + row_shift
1887 ibeg_loc = ibeg + row_shift
1888 iend_loc = iend + row_shift
1889 k1 = jj_loc
1890 k2 = jj_loc+1_8
1891 mult1 = - (a11*a(k1)+a12*a(k2))
1892 mult2 = - (a12*a(k1)+a22*a(k2))
1893 a(pospv1 + 2_8 + nel2 + (j2-1_8)) = a(k1)
1894 a(pospv2 + 1_8 + nel2 + (j2-1_8)) = a(k2)
1895 k1 = pospv1 + 2_8
1896 k2 = pospv2 + 1_8
1897 DO irow = ibeg_loc, iend_loc
1898 a(irow) = a(irow) + mult1*a(k1) + mult2*a(k2)
1899 k1 = k1 + 1_8
1900 k2 = k2 + 1_8
1901 ENDDO
1902 a( jj_loc ) = -mult1
1903 a( jj_loc + 1_8 ) = -mult2
1904 ENDDO
1905!$OMP END PARALLEL DO
1906#endif
1907 ENDIF
1908 IF ((is_maxfromm_avail).AND.(nel2.GT.0)) THEN
1909 IF (parpiv_t1.NE.0) THEN
1910 aposmax = poselt+lda8*lda8-1_8 + int(npiv_new+1,8)
1911 maxfromm = max(maxfromm,
1912 & dble(a(aposmax))
1913 & )
1914 ENDIF
1915 ENDIF
1916 RETURN
1917 END SUBROUTINE dmumps_fac_mq_ldlt
1918 SUBROUTINE dmumps_fac_sq_ldlt(IBEG_BLOCK,IEND_BLOCK,NPIV,
1919 & NFRONT,NASS,INODE,A,LA,
1920 & LDA,
1921 & POSELT,
1922 & KEEP,KEEP8,
1923 & FIRST_ROW_TRSM, LAST_ROW_TRSM,
1924 & LAST_COL_GEMM, LAST_ROW_GEMM,
1925 & CALL_TRSM, CALL_GEMM, LR_ACTIVATED,
1926 & IW, LIW, OFFSET_IW
1927 & )
1928 IMPLICIT NONE
1929 INTEGER, intent(in) :: NPIV
1930 INTEGER, intent(in) :: NFRONT, NASS, IBEG_BLOCK, IEND_BLOCK
1931 INTEGER(8), intent(in) :: LA
1932 DOUBLE PRECISION, intent(inout) :: A(LA)
1933 INTEGER, intent(in) :: INODE
1934 INTEGER :: KEEP(500)
1935 INTEGER(8) :: KEEP8(150)
1936 INTEGER(8), intent(in) :: POSELT
1937 INTEGER, intent(in) :: LDA
1938 INTEGER, intent(in) :: LAST_COL_GEMM
1939 INTEGER, intent(in) :: LAST_ROW_GEMM, LAST_ROW_TRSM,
1940 & first_row_trsm
1941 LOGICAL, intent(in) :: CALL_TRSM, CALL_GEMM, LR_ACTIVATED
1942 INTEGER :: OFFSET_IW, LIW
1943 INTEGER :: IW(LIW)
1944 INTEGER(8) :: LDA8
1945 INTEGER NPIV_BLOCK, NEL1
1946 INTEGER NRHS_TRSM
1947 INTEGER(8) :: LPOS, UPOS, APOS
1948 INTEGER IROW
1949 INTEGER Block
1950 INTEGER BLSIZE
1951 DOUBLE PRECISION ONE, ALPHA
1952 include 'mumps_headers.h'
1953 parameter(one=1.0d0, alpha=-1.0d0)
1954 lda8 = int(lda,8)
1955 nel1 = last_col_gemm - iend_block
1956 nrhs_trsm = last_row_trsm-first_row_trsm
1957 npiv_block = npiv - ibeg_block + 1
1958 IF (npiv_block.EQ.0) GO TO 500
1959 IF (nel1.NE.0) THEN
1960 IF (call_trsm) THEN
1961 apos = poselt + lda8*int(ibeg_block-1,8) + int(ibeg_block-1,8)
1962 lpos = poselt + lda8*int(first_row_trsm,8)+int(ibeg_block-1,8)
1963 upos = poselt + lda8*int(ibeg_block-1,8)+int(first_row_trsm,8)
1964 CALL dtrsm('L', 'U', 'T', 'U', npiv_block, nrhs_trsm,
1965 & one, a(apos), lda, a(lpos), lda)
1966 CALL dmumps_fac_ldlt_copy2u_scalel(nrhs_trsm, 1, keep(424),
1967 & nfront, npiv_block, liw, iw, offset_iw, la, a,
1968 & poselt, lpos, upos, apos, .NOT.lr_activated)
1969 ENDIF
1970 IF (call_gemm) THEN
1971#if defined(GEMMT_AVAILABLE)
1972 IF ( keep(421).EQ. -1) THEN
1973 lpos = poselt + lda8*int(iend_block,8) + int(ibeg_block-1,8)
1974 upos = poselt + lda8*int(ibeg_block-1,8) + int(iend_block,8)
1975 apos = poselt + lda8*int(iend_block,8) + int(iend_block,8)
1976 CALL dgemmt( 'U','N','N', nel1,
1977 & npiv_block,
1978 & alpha, a( upos ), lda,
1979 & a( lpos ), lda, one, a( apos ), lda )
1980 ELSE
1981#endif
1982 IF ( last_col_gemm - iend_block > keep(7) ) THEN
1983 blsize = keep(8)
1984 ELSE
1985 blsize = last_col_gemm - iend_block
1986 END IF
1987 IF ( last_col_gemm - iend_block .GT. 0 ) THEN
1988#if defined(SAK_BYROW)
1989 DO irow = iend_block+1, last_col_gemm, blsize
1990 block = min( blsize, last_col_gemm - irow + 1 )
1991 lpos = poselt + int(irow - 1,8) * lda8 +
1992 & int(ibeg_block - 1,8)
1993 upos = poselt + int(ibeg_block - 1,8) * lda8 +
1994 & int(irow - 1,8)
1995 apos = poselt + int(irow - 1,8) * lda8 +
1996 & int(iend_block,8)
1997 CALL dgemm( 'N','N', irow + block - iend_block - 1,
1998 & block, npiv_block,
1999 & alpha, a( upos ), lda,
2000 & a( lpos ), lda, one, a( apos ), lda )
2001 ENDDO
2002#else
2003 DO irow = iend_block+1, last_col_gemm, blsize
2004 block = min( blsize, last_col_gemm - irow + 1 )
2005 lpos = poselt + int( irow - 1,8) * lda8 +
2006 & int(ibeg_block - 1,8)
2007 upos = poselt + int(ibeg_block - 1,8) * lda8 +
2008 & int( irow - 1,8)
2009 apos = poselt + int( irow - 1,8) * lda8 + int( irow - 1,8)
2010 CALL dgemm( 'N','N', block, last_col_gemm - irow + 1,
2011 & npiv_block, alpha, a( upos ), lda,
2012 & a( lpos ), lda, one, a( apos ), lda )
2013 END DO
2014#endif
2015 END IF
2016#if defined(GEMMT_AVAILABLE)
2017 END IF
2018#endif
2019 lpos = poselt + int(last_col_gemm,8)*lda8 + int(ibeg_block-1,8)
2020 upos = poselt + int(ibeg_block-1,8) * lda8 +
2021 & int(iend_block,8)
2022 apos = poselt + int(last_col_gemm,8)*lda8 + int(iend_block,8)
2023 IF (last_row_gemm .GT. last_col_gemm) THEN
2024 CALL dgemm('N', 'N', nel1, last_row_gemm-last_col_gemm,
2025 & npiv_block, alpha, a(upos), lda, a(lpos), lda,
2026 & one, a(apos), lda)
2027 ENDIF
2028 ENDIF
2029 ENDIF
2030 500 CONTINUE
2031 RETURN
2032 END SUBROUTINE dmumps_fac_sq_ldlt
2033 SUBROUTINE dmumps_swap_ldlt( A, LA, IW, LIW,
2034 & IOLDPS, NPIVP1, IPIV, POSELT, LASTROW2SWAP,
2035 & LDA, NFRONT, LEVEL, PARPIV, K50, XSIZE,
2036 & IBEG_BLOCK_TO_SEND )
2037 IMPLICIT NONE
2038 INTEGER(8) :: POSELT, LA
2039 INTEGER LIW, IOLDPS, NPIVP1, IPIV
2040 INTEGER LDA, NFRONT, LEVEL, PARPIV, K50, XSIZE
2041 INTEGER LASTROW2SWAP
2042 DOUBLE PRECISION A( LA )
2043 INTEGER IW( LIW )
2044 INTEGER, INTENT(IN) :: IBEG_BLOCK_TO_SEND
2045 INCLUDE 'mumps_headers.h'
2046 INTEGER :: IBEG
2047 INTEGER ISW, ISWPS1, ISWPS2, HF
2048 INTEGER(8) :: IDIAG, APOS
2049 INTEGER(8) :: LDA8
2050 DOUBLE PRECISION SWOP
2051 lda8 = int(lda,8)
2052 apos = poselt + lda8*int(ipiv-1,8) + int(npivp1-1,8)
2053 idiag = apos + int(ipiv - npivp1,8)
2054 hf = 6 + iw( ioldps + 5 + xsize) + xsize
2055 iswps1 = ioldps + hf + npivp1 - 1
2056 iswps2 = ioldps + hf + ipiv - 1
2057 isw = iw(iswps1)
2058 iw(iswps1) = iw(iswps2)
2059 iw(iswps2) = isw
2060 isw = iw(iswps1+nfront)
2061 iw(iswps1+nfront) = iw(iswps2+nfront)
2062 iw(iswps2+nfront) = isw
2063 IF ( level .eq. 2 ) THEN
2064 ibeg = ibeg_block_to_send
2065 CALL dswap( npivp1 - 1 - ibeg + 1,
2066 & a( poselt + int(npivp1-1,8) +
2067 & int(ibeg-1,8) * lda8), lda,
2068 & a( poselt + int(ipiv-1,8) +
2069 & int(ibeg-1,8) * lda8), lda )
2070 END IF
2071 CALL dswap( npivp1-1,
2072 & a( poselt+int(npivp1-1,8) * lda8 ), 1,
2073 & a( poselt + int(ipiv-1,8) * lda8 ), 1 )
2074 CALL dswap( ipiv - npivp1 - 1,
2075 & a( poselt+int(npivp1,8) * lda8 + int(npivp1-1,8) ),
2076 & lda, a( apos + 1_8 ), 1 )
2077 swop = a(idiag)
2078 a(idiag) = a( poselt+int(npivp1-1,8)*lda8+int(npivp1-1,8) )
2079 a( poselt + int(npivp1-1,8)*lda8 + int(npivp1-1,8) ) = swop
2080 IF (lastrow2swap - ipiv.GT.0) THEN
2081 CALL dswap( lastrow2swap - ipiv,
2082 & a( apos + lda8 ), lda,
2083 & a( idiag + lda8 ), lda )
2084 ENDIF
2085 IF (parpiv.NE.0 .AND.k50.EQ.2) THEN
2086 IF ( level .eq. 2 .OR. level.eq.1) THEN
2087 apos = poselt+lda8*lda8-1_8
2088 swop = a(apos+int(npivp1,8))
2089 a(apos+int(npivp1,8))= a(apos+int(ipiv,8))
2090 a(apos+int(ipiv,8)) = swop
2091 ENDIF
2092 ENDIF
2093 RETURN
2094 END SUBROUTINE dmumps_swap_ldlt
2095 SUBROUTINE dmumps_fac_ldlt_copy2u_scalel( IROWMAX, IROWMIN,
2096 & SIZECOPY, LDA, NCOLS, LIW, IW, OFFSET_IW,
2097 & LA, A, POSELT, A_LPOS, A_UPOS, A_DPOS,
2098 & COPY_NEEDED )
2099!$ USE OMP_LIB
2100 INTEGER, INTENT(IN) :: IROWMAX, IROWMIN
2101 INTEGER, INTENT(IN) :: SIZECOPY
2102 INTEGER, INTENT(IN) :: LDA, NCOLS
2103 INTEGER, INTENT(IN) :: LIW
2104 INTEGER, INTENT(IN) :: IW(LIW)
2105 INTEGER, INTENT(IN) :: OFFSET_IW
2106 INTEGER(8), INTENT(IN) :: LA
2107 DOUBLE PRECISION, INTENT(INOUT) :: A(LA)
2108 INTEGER(8), INTENT(IN) :: POSELT, A_LPOS, A_UPOS, A_DPOS
2109 LOGICAL, INTENT(IN) :: COPY_NEEDED
2110 INTEGER(8) :: LPOS, UPOS
2111 INTEGER(8) :: DPOS, POSPV1, POSPV2, OFFDAG
2112 INTEGER(8) :: LDA8
2113 INTEGER :: IROWEND, IROW, Block2
2114 INTEGER :: I, J
2115 DOUBLE PRECISION :: MULT1, MULT2, A11, DETPIV, A22, A12
2116 INTEGER :: BLSIZECOPY
2117 DOUBLE PRECISION :: ONE
2118 parameter(one = 1.0d0)
2119 INTEGER(8) :: LPOSI, UPOSI
2120 LOGICAL :: PIVOT_2X2
2121!$ LOGICAL :: OMP_FLAG
2122!$ INTEGER :: NOMP, CHUNK
2123 lda8 = int(lda,8)
2124 IF (sizecopy.NE.0) THEN
2125 blsizecopy = sizecopy
2126 ELSE
2127 blsizecopy = 250
2128 ENDIF
2129!$ NOMP = OMP_GET_MAX_THREADS()
2130!$ omp_flag = .false.
2131!$ CHUNK = (64/4)
2132!$ IF (NOMP .GT. 1 .AND. NCOLS .GE. 4*CHUNK) THEN
2133!$ OMP_FLAG = .TRUE.
2134!$ CHUNK = max(2*CHUNK, NCOLS/NOMP)
2135!$ ENDIF
2136 DO irowend = irowmax, irowmin, -blsizecopy
2137 block2 = min(blsizecopy, irowend)
2138 irow = irowend - block2 + 1
2139 lpos = a_lpos + int(irow-1,8)*lda8
2140 upos = a_upos + int(irow-1,8)
2141!$OMP PARALLEL DO PRIVATE(PIVOT_2X2, A11, DPOS,
2142!$OMP& POSPV1, POSPV2, OFFDAG, A22, A12, DETPIV, J, MULT1, MULT2
2143!$OMP& , LPOSI, UPOSI
2144!$OMP& ) FIRSTPRIVATE(Block2, LDA, LDA8, LPOS, UPOS, A_DPOS)
2145!$OMP& SCHEDULE(STATIC,CHUNK) IF(OMP_FLAG)
2146 DO i=1, ncols
2147 pivot_2x2 = .false.
2148 IF(iw(offset_iw+i-1) .LE. 0) THEN
2149 pivot_2x2 = .true.
2150 ELSE
2151 IF (i .GT. 1) THEN
2152 IF (iw(offset_iw+i-2) .LE. 0) THEN
2153 cycle
2154 ENDIF
2155 ENDIF
2156 ENDIF
2157 dpos = a_dpos + lda8*int(i-1,8) + int(i-1,8)
2158 IF(.not. pivot_2x2) THEN
2159 a11 = one/a(dpos)
2160 lposi = lpos+int(i-1,8)
2161 IF (copy_needed) THEN
2162 uposi = upos+int(i-1,8)*lda8
2163#if defined(__ve__)
2164!NEC$ IVDEP
2165#endif
2166 DO j = 1, block2
2167 a(uposi+int(j-1,8)) = a(lposi+int(j-1,8)*lda8)
2168 END DO
2169 ENDIF
2170#if defined(__ve__)
2171!NEC$ IVDEP
2172#endif
2173 DO j = 1, block2
2174 a(lposi+int(j-1,8)*lda8) = a(lposi+int(j-1,8)*lda8)*a11
2175 END DO
2176 ELSE
2177 IF (copy_needed) THEN
2178 CALL dcopy(block2, a(lpos+int(i-1,8)),
2179 & lda, a(upos+int(i-1,8)*lda8), 1)
2180 CALL dcopy(block2, a(lpos+int(i,8)),
2181 & lda, a(upos+int(i,8)*lda8), 1)
2182 ENDIF
2183 pospv1 = dpos
2184 pospv2 = dpos + int(lda+1,8)
2185 offdag = pospv1+1_8
2186 a11 = a(pospv1)
2187 a22 = a(pospv2)
2188 a12 = a(offdag)
2189 detpiv = a11*a22 - a12**2
2190 a22 = a11/detpiv
2191 a11 = a(pospv2)/detpiv
2192 a12 = -a12/detpiv
2193#if defined(__ve__)
2194!NEC$ IVDEP
2195#endif
2196 DO j = 1,block2
2197 mult1 = a11*a(lpos+int(j-1,8)*lda8+int(i-1,8))
2198 & + a12*a(lpos+int(j-1,8)*lda8+int(i,8))
2199 mult2 = a12*a(lpos+int(j-1,8)*lda8+int(i-1,8))
2200 & + a22*a(lpos+int(j-1,8)*lda8+int(i,8))
2201 a(lpos+int(j-1,8)*lda8+int(i-1,8)) = mult1
2202 a(lpos+int(j-1,8)*lda8+int(i,8)) = mult2
2203 ENDDO
2204 ENDIF
2205 ENDDO
2206!$OMP END PARALLEL DO
2207 ENDDO
2208 END SUBROUTINE dmumps_fac_ldlt_copy2u_scalel
2209 SUBROUTINE dmumps_fac_ldlt_copyscale_u( IROWMAX, IROWMIN,
2210 & SIZECOPY, LDA, NCOLS, LIW, IW, OFFSET_IW,
2211 & LA, A, POSELT, A_LPOS, A_UPOS, A_DPOS )
2212!$ USE OMP_LIB
2213 INTEGER, INTENT(IN) :: IROWMAX, IROWMIN
2214 INTEGER, INTENT(IN) :: SIZECOPY
2215 INTEGER, INTENT(IN) :: LDA, NCOLS
2216 INTEGER, INTENT(IN) :: LIW
2217 INTEGER, INTENT(IN) :: IW(LIW)
2218 INTEGER, INTENT(IN) :: OFFSET_IW
2219 INTEGER(8), INTENT(IN) :: LA
2220 DOUBLE PRECISION, INTENT(INOUT) :: A(LA)
2221 INTEGER(8), INTENT(IN) :: POSELT, A_LPOS, A_UPOS, A_DPOS
2222 INTEGER(8) :: LPOS, UPOS
2223 INTEGER(8) :: DPOS, POSPV1, POSPV2, OFFDAG
2224 INTEGER(8) :: LDA8
2225 INTEGER :: IROWEND, IROW, Block2
2226 INTEGER :: I, J
2227 DOUBLE PRECISION :: MULT1, MULT2, A11, A22, A12
2228 INTEGER :: BLSIZECOPY
2229 DOUBLE PRECISION :: ONE
2230 parameter(one = 1.0d0)
2231 INTEGER(8) :: LPOSI, UPOSI
2232 LOGICAL :: PIVOT_2X2
2233!$ LOGICAL :: OMP_FLAG
2234!$ INTEGER :: NOMP, CHUNK
2235 lda8 = int(lda,8)
2236 IF (sizecopy.NE.0) THEN
2237 blsizecopy = sizecopy
2238 ELSE
2239 blsizecopy = 250
2240 ENDIF
2241!$ NOMP = OMP_GET_MAX_THREADS()
2242!$ OMP_FLAG = .FALSE.
2243!$ CHUNK = (64/4)
2244!$ IF (NOMP .GT. 1 .AND. NCOLS .GE. 4*CHUNK) THEN
2245!$ OMP_FLAG = .TRUE.
2246!$ CHUNK = max(2*CHUNK, NCOLS/NOMP)
2247!$ ENDIF
2248 DO irowend = irowmax, irowmin, -blsizecopy
2249 block2 = min(blsizecopy, irowend)
2250 irow = irowend - block2 + 1
2251 lpos = a_lpos + int(irow-1,8)*lda8
2252 upos = a_upos + int(irow-1,8)
2253!$OMP PARALLEL DO PRIVATE(PIVOT_2X2, A11, DPOS,
2254!$OMP& POSPV1, POSPV2, OFFDAG, A22, A12, J, MULT1, MULT2
2255!$OMP& , LPOSI, UPOSI
2256!$OMP& ) FIRSTPRIVATE(Block2, LDA, LDA8, LPOS, UPOS, POSELT)
2257!$OMP& SCHEDULE(STATIC,CHUNK) IF(OMP_FLAG)
2258 DO i=1, ncols
2259 pivot_2x2 = .false.
2260 IF(iw(offset_iw+i-1) .LE. 0) THEN
2261 pivot_2x2 = .true.
2262 ELSE
2263 IF (i .GT. 1) THEN
2264 IF (iw(offset_iw+i-2) .LE. 0) THEN
2265 cycle
2266 ENDIF
2267 ENDIF
2268 ENDIF
2269 dpos = a_dpos + lda8*int(i-1,8) + int(i-1,8)
2270 IF(.not. pivot_2x2) THEN
2271 a11 = a(dpos)
2272 lposi = lpos+int(i-1,8)
2273 uposi = upos+int(i-1,8)*lda8
2274#if defined(__ve__)
2275!NEC$ IVDEP
2276#endif
2277 DO j = 1, block2
2278 a(uposi+int(j-1,8)) = a(lposi+int(j-1,8)*lda8)*a11
2279 END DO
2280 ELSE
2281 pospv1 = dpos
2282 pospv2 = dpos + int(lda+1,8)
2283 offdag = pospv1+1_8
2284 a11 = a(pospv1)
2285 a22 = a(pospv2)
2286 a12 = a(offdag)
2287#if defined(__ve__)
2288!nec$ ivdep
2289#endif
2290 DO j = 1,block2
2291 mult1 = a11*a(lpos+int(j-1,8)*lda8+int(i-1,8))
2292 & + a12*a(lpos+int(j-1,8)*lda8+int(i,8))
2293 mult2 = a12*a(lpos+int(j-1,8)*lda8+int(i-1,8))
2294 & + a22*a(lpos+int(j-1,8)*lda8+int(i,8))
2295 a(upos+int(i-1,8)*lda8+int(j-1,8)) = mult1
2296 a(upos+int(i,8)*lda8+int(j-1,8)) = mult2
2297 ENDDO
2298 ENDIF
2299 ENDDO
2300!$OMP END PARALLEL DO
2301 ENDDO
2302 RETURN
2303 END SUBROUTINE dmumps_fac_ldlt_copyscale_u
2304 SUBROUTINE dmumps_fac_t_ldlt(NFRONT,NASS,
2305 & IW,LIW,A,LA,
2306 & LDA,
2307 & IOLDPS,POSELT,KEEP,KEEP8,
2308 & POSTPONE_COL_UPDATE, ETATASS,
2309 & TYPEFile, LAFAC, MonBloc, NextPiv2beWritten,
2310 & LIWFAC, MYID, IFLAG, OFFSET_IW, INODE )
2311 USE dmumps_ooc
2312 IMPLICIT NONE
2313 INTEGER NFRONT, NASS,LIW
2314 INTEGER(8) :: LA
2315 DOUBLE PRECISION A(LA)
2316 INTEGER IW(LIW)
2317 INTEGER KEEP(500)
2318 INTEGER(8) KEEP8(150)
2319 INTEGER(8) :: POSELT
2320 INTEGER LDA
2321 INTEGER IOLDPS, ETATASS
2322 LOGICAL POSTPONE_COL_UPDATE
2323 INTEGER(8) :: LAFAC
2324 INTEGER TYPEFile, NextPiv2beWritten
2325 INTEGER LIWFAC, MYID, IFLAG
2326 TYPE(IO_BLOCK):: MonBloc
2327 INTEGER IDUMMY
2328 LOGICAL LAST_CALL
2329 INTEGER :: OFFSET_IW
2330 INTEGER, intent(in):: INODE
2331 include 'mumps_headers.h'
2332 INTEGER(8) :: UPOS, APOS, LPOS
2333 INTEGER(8) :: LDA8
2334 INTEGER BLSIZE, BLSIZE2, Block, IROW, NPIV, IROWEND
2335 INTEGER I2, I2END, Block2
2336 DOUBLE PRECISION ONE, ALPHA, BETA, ZERO
2337 parameter(one = 1.0d0, alpha=-1.0d0)
2338 parameter(zero=0.0d0)
2339 lda8 = int(lda,8)
2340 IF (etatass.EQ.1) THEN
2341 beta = zero
2342 ELSE
2343 beta = one
2344 ENDIF
2345 IF ( nfront - nass > keep(58) ) THEN
2346 IF ( nfront - nass > keep(57) ) THEN
2347 blsize = keep(58)
2348 ELSE
2349 blsize = (nfront - nass)/2
2350 END IF
2351 ELSE
2352 blsize = nfront - nass
2353 END IF
2354 blsize2 = keep(218)
2355 npiv = iw( ioldps + 1 + keep(ixsz))
2356 IF ( nfront - nass .GT. 0 ) THEN
2357 IF ( postpone_col_update ) THEN
2358 lpos = poselt + lda8 * int(nass,8)
2359 CALL dtrsm( 'L', 'U', 'T', 'U',
2360 & npiv, nfront-nass, one,
2361 & a( poselt ), lda,
2362 & a( lpos ), lda )
2363 ENDIF
2364#if defined(GEMMT_AVAILABLE)
2365 IF ( keep(421).EQ. -1) THEN
2366 lpos = poselt + int(nass,8)*lda8
2367 upos = poselt + int(nass,8)
2368 apos = poselt + int(nass,8)*lda8 + int(nass,8)
2369 IF (postpone_col_update) THEN
2370 CALL dmumps_fac_ldlt_copy2u_scalel( nfront - nass, 1,
2371 & keep(424), nfront, npiv,
2372 & liw, iw, offset_iw, la, a, poselt, lpos, upos,
2373 & poselt, .true. )
2374 ENDIF
2375 CALL dgemmt('U', 'N', 'N', nfront-nass, npiv,
2376 & alpha, a( upos ), lda,
2377 & a( lpos ), lda,
2378 & beta,
2379 & a( apos ), lda )
2380 ELSE
2381#endif
2382 DO irowend = nfront - nass, 1, -blsize
2383 block = min( blsize, irowend )
2384 irow = irowend - block + 1
2385 lpos = poselt + int(nass,8)*lda8 + int(irow-1,8) * lda8
2386 apos = poselt + int(nass,8)*lda8 + int(irow-1,8) * lda8 +
2387 & int(nass + irow - 1,8)
2388 upos = poselt + int(nass,8)
2389 IF (.NOT. postpone_col_update) THEN
2390 upos = poselt + int(nass + irow - 1,8)
2391 ENDIF
2392 IF (postpone_col_update) THEN
2393 CALL dmumps_fac_ldlt_copy2u_scalel( block, 1,
2394 & keep(424), nfront, npiv,
2395 & liw, iw, offset_iw, la, a, poselt, lpos, upos,
2396 & poselt, .true. )
2397 ENDIF
2398 DO i2end = block, 1, -blsize2
2399 block2 = min(blsize2, i2end)
2400 i2 = i2end - block2+1
2401 CALL dgemm('N', 'N', block2, block-i2+1, npiv, alpha,
2402 & a(upos+int(i2-1,8)), lda,
2403 & a(lpos+int(i2-1,8)*lda8), lda,
2404 & beta,
2405 & a(apos + int(i2-1,8) + int(i2-1,8)*lda8), lda)
2406 IF (keep(201).EQ.1) THEN
2407 IF (nextpiv2bewritten.LE.npiv) THEN
2408 last_call=.false.
2409 CALL dmumps_ooc_io_lu_panel(
2410 & strat_try_write, typefile,
2411 & a(poselt), lafac, monbloc,
2412 & nextpiv2bewritten, idummy,
2413 & iw(ioldps), liwfac, myid,
2414 & keep8(31),
2415 & iflag,last_call )
2416 IF (iflag .LT. 0 ) RETURN
2417 ENDIF
2418 ENDIF
2419 ENDDO
2420 IF ( nfront - nass - irow + 1 - block > 0 ) THEN
2421 CALL dgemm( 'N', 'N', block, nfront-nass-block-irow+1, npiv,
2422 & alpha, a( upos ), lda,
2423 & a( lpos + lda8 * int(block,8) ), lda,
2424 & beta,
2425 & a( apos + lda8 * int(block,8) ), lda )
2426 ENDIF
2427 END DO
2428#if defined(GEMMT_AVAILABLE)
2429 END IF
2430#endif
2431 IF ( (postpone_col_update).AND.(nass-npiv.GT.0) ) THEN
2432 lpos = poselt + int(npiv,8)*lda8
2433 upos = poselt + int(npiv,8)
2434 CALL dmumps_fac_ldlt_copyscale_u( nass-npiv, 1,
2435 & keep(424), nfront, npiv,
2436 & liw, iw, offset_iw, la, a, poselt, lpos, upos, poselt)
2437 lpos = poselt + lda8 * int(nass,8)
2438 CALL dgemm('N', 'N', nass-npiv, nfront-nass, npiv, alpha,
2439 & a( poselt + int(npiv,8)), lda,
2440 & a( lpos ), lda,
2441 & beta,
2442 & a( lpos + int(npiv,8) ), lda)
2443 ENDIF
2444 END IF
2445 RETURN
2446 END SUBROUTINE dmumps_fac_t_ldlt
2447 SUBROUTINE dmumps_store_perminfo( PIVRPTR, NBPANELS, PIVR, NASS,
2448 & K, P, LastPanelonDisk,
2449 & LastPIVRPTRIndexFilled )
2450 IMPLICIT NONE
2451 INTEGER, intent(in) :: NBPANELS, NASS, K, P
2452 INTEGER, intent(inout) :: PIVRPTR(NBPANELS), PIVR(NASS)
2453 INTEGER LastPanelonDisk, LastPIVRPTRIndexFilled
2454 INTEGER I
2455 IF ( lastpanelondisk+1 > nbpanels ) THEN
2456 WRITE(*,*) "INTERNAL ERROR IN DMUMPS_STORE_PERMINFO!"
2457 WRITE(*,*) "NASS=",nass,"PIVRPTR=",pivrptr(1:nbpanels)
2458 WRITE(*,*) "K=",k, "P=",p, "LastPanelonDisk=",lastpanelondisk
2459 WRITE(*,*) "LastPIVRPTRIndexFilled=", lastpivrptrindexfilled
2460 CALL mumps_abort()
2461 ENDIF
2462 pivrptr(lastpanelondisk+1) = k + 1
2463 IF (lastpanelondisk.NE.0) THEN
2464 pivr(k - pivrptr(1) + 1) = p
2465 DO i = lastpivrptrindexfilled + 1, lastpanelondisk
2466 pivrptr(i)=pivrptr(lastpivrptrindexfilled)
2467 ENDDO
2468 ENDIF
2469 lastpivrptrindexfilled = lastpanelondisk + 1
2470 RETURN
2471 END SUBROUTINE dmumps_store_perminfo
2472 SUBROUTINE dmumps_update_minmax_pivot
2473 & ( diag, dkeep, keep, nullpivot)
2474!$ USE OMP_LIB
2475 IMPLICIT NONE
2476 DOUBLE PRECISION, INTENT(IN) :: DIAG
2477 DOUBLE PRECISION, INTENT(INOUT) :: DKEEP(230)
2478 LOGICAL, INTENT(IN) :: NULLPIVOT
2479 INTEGER, INTENT(IN) :: KEEP(500)
2480 IF (keep(405).EQ.0) THEN
2481 dkeep(21) = max(dkeep(21), diag)
2482 dkeep(19) = min(dkeep(19), diag)
2483 IF (.NOT.nullpivot) THEN
2484 dkeep(20) = min(dkeep(20), diag)
2485 ENDIF
2486 ELSE
2487!$omp atomic update
2488 dkeep(21) = max(dkeep(21), diag)
2489!$OMP END ATOMIC
2490!$OMP ATOMIC UPDATE
2491 dkeep(19) = min(dkeep(19), diag)
2492!$OMP END ATOMIC
2493 IF (.NOT.nullpivot) THEN
2494!$OMP ATOMIC UPDATE
2495 dkeep(20) = min(dkeep(20), diag)
2496!$OMP END ATOMIC
2497 ENDIF
2498 ENDIF
2499 RETURN
2500 END SUBROUTINE dmumps_update_minmax_pivot
2501 SUBROUTINE dmumps_get_size_schur_in_front (
2502 & N, NCB, SIZE_SCHUR, ROW_INDICES, PERM,
2503 & NVSCHUR
2504 & )
2505 IMPLICIT NONE
2506 INTEGER, intent(in) :: N, NCB, SIZE_SCHUR
2507 INTEGER, intent(in) :: ROW_INDICES(NCB), PERM(N)
2508 INTEGER, intent(out):: NVSCHUR
2509 INTEGER :: I, IPOS, IBEG_SCHUR
2510 ibeg_schur = n - size_schur +1
2511 nvschur = 0
2512 ipos = ncb
2513 DO i= ncb,1,-1
2514 IF (abs(row_indices(i)).LE.n) THEN
2515 IF (perm(row_indices(i)).LT.ibeg_schur) EXIT
2516 ENDIF
2517 ipos = ipos -1
2518 ENDDO
2519 nvschur = ncb-ipos
2520 RETURN
2521 END SUBROUTINE dmumps_get_size_schur_in_front
2522 END MODULE dmumps_fac_front_aux_m
mumps_abort
#define mumps_abort
Definition VE_Metis.h:25
dmumps_updatedeter
subroutine dmumps_updatedeter(piv, deter, nexp)
Definition dfac_determinant.F:15
dmumps_get_ooc_perm_ptr
subroutine dmumps_get_ooc_perm_ptr(typef, nbpanels, i_pivptr, i_piv, ipos, iw, liw)
Definition dooc_panel_piv.F:133
alpha
#define alpha
Definition eval.h:35
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dger
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130
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
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
dmumps_buf
Definition dmumps_comm_buffer.F:14
dmumps_fac_front_aux_m
Definition dfac_front_aux.F:14
dmumps_fac_front_aux_m::dmumps_fac_h
subroutine dmumps_fac_h(nfront, nass, iw, liw, a, la, inopv, noffw, det_expw, det_mantw, det_signw, ioldps, poselt, uu, seuil, keep, keep8, dkeep, pp_first2swap_l, pp_lastpanelondisk_l, pp_lastpivrptrfilled_l, pp_first2swap_u, pp_lastpanelondisk_u, pp_lastpivrptrfilled_u, maxfromn, is_maxfromn_avail, inextpiv, ooc_effective_on_front, nvschur)
Definition dfac_front_aux.F:26
dmumps_fac_front_aux_m::dmumps_store_perminfo
subroutine dmumps_store_perminfo(pivrptr, nbpanels, pivr, nass, k, p, lastpanelondisk, lastpivrptrindexfilled)
Definition dfac_front_aux.F:2450
dmumps_fac_front_aux_m::dmumps_fac_ldlt_copyscale_u
subroutine dmumps_fac_ldlt_copyscale_u(irowmax, irowmin, sizecopy, lda, ncols, liw, iw, offset_iw, la, a, poselt, a_lpos, a_upos, a_dpos)
Definition dfac_front_aux.F:2212
dmumps_fac_front_aux_m::dmumps_fac_n
subroutine dmumps_fac_n(nfront, nass, iw, liw, a, la, ioldps, poselt, ifinb, xsize, keep, maxfromn, is_maxfromn_avail, nvschur)
Definition dfac_front_aux.F:239
dmumps_fac_front_aux_m::dmumps_fac_t
subroutine dmumps_fac_t(a, la, npivb, nfront, npiv, nass, poselt)
Definition dfac_front_aux.F:413
dmumps_fac_front_aux_m::dmumps_fac_mq
subroutine dmumps_fac_mq(ibeg_block, iend_block, nfront, nass, npiv, last_col, a, la, poselt, ifinb, lr_activated)
Definition dfac_front_aux.F:573
dmumps_fac_front_aux_m::dmumps_fac_ldlt_copy2u_scalel
subroutine dmumps_fac_ldlt_copy2u_scalel(irowmax, irowmin, sizecopy, lda, ncols, liw, iw, offset_iw, la, a, poselt, a_lpos, a_upos, a_dpos, copy_needed)
Definition dfac_front_aux.F:2099
dmumps_fac_front_aux_m::dmumps_fac_i
subroutine dmumps_fac_i(nfront, nass, last_row, ibeg_block, iend_block, n, inode, iw, liw, a, la, inopv, noffw, nbtinyw, det_expw, det_mantw, det_signw, iflag, ioldps, poselt, uu, seuil, keep, keep8, dkeep, pivnul_list, lpn_list pp_first2swap_l, pp_lastpanelondisk_l, pp_lastpivrptrfilled_l, pp_first2swap_u, pp_lastpanelondisk_u, pp_lastpivrptrfilled_u, pivot_option, lr_activated, iend_blr, inextpiv, ooc_effective_on_front, nvschur, parpiv_t1, tipiv)
Definition dfac_front_aux.F:720
dmumps_fac_front_aux_m::dmumps_fac_p
subroutine dmumps_fac_p(a, la, nfront, npiv, nass, poselt, call_utrsm, keep, inode, call_ooc, iwfac, liwfac, lafac, monbloc, myid, keep8, lnextpiv2bewritten, unextpiv2bewritten, iflag)
Definition dfac_front_aux.F:353
dmumps_fac_front_aux_m::dmumps_get_size_schur_in_front
subroutine dmumps_get_size_schur_in_front(n, ncb, size_schur, row_indices, perm, nvschur)
Definition dfac_front_aux.F:2505
dmumps_fac_front_aux_m::dmumps_fac_mq_ldlt
subroutine dmumps_fac_mq_ldlt(iend_block, nfront, nass, npiv, inode, a, la, lda, poselt, ifinb, pivsiz, maxfromm, is_maxfromm_avail, is_max_useful, parpiv_t1, last_row, iend_blr, nvschur_k253, lr_activated)
Definition dfac_front_aux.F:1609
dmumps_fac_front_aux_m::dmumps_fac_sq
subroutine dmumps_fac_sq(ibeg_block, iend_block, npiv, nfront, last_row, last_col, a, la, poselt, first_col, call_ltrsm, call_utrsm, call_gemm, with_comm_thread, lr_activated)
Definition dfac_front_aux.F:443
dmumps_fac_front_aux_m::dmumps_swap_ldlt
subroutine dmumps_swap_ldlt(a, la, iw, liw, ioldps, npivp1, ipiv, poselt, lastrow2swap, lda, nfront, level, parpiv, k50, xsize, ibeg_block_to_send)
Definition dfac_front_aux.F:2037
dmumps_fac_front_aux_m::dmumps_update_minmax_pivot
subroutine dmumps_update_minmax_pivot(diag, dkeep, keep, nullpivot)
Definition dfac_front_aux.F:2474
dmumps_fac_front_aux_m::dmumps_fac_pt_setlock427
subroutine dmumps_fac_pt_setlock427(k427_out, k427, k405, k222, nel1, nass)
Definition dfac_front_aux.F:339
dmumps_fac_front_aux_m::dmumps_fac_sq_ldlt
subroutine dmumps_fac_sq_ldlt(ibeg_block, iend_block, npiv, nfront, nass, inode, a, la, lda, poselt, keep, keep8, first_row_trsm, last_row_trsm, last_col_gemm, last_row_gemm, call_trsm, call_gemm, lr_activated, iw, liw, offset_iw)
Definition dfac_front_aux.F:1928
dmumps_fac_front_aux_m::dmumps_fac_t_ldlt
subroutine dmumps_fac_t_ldlt(nfront, nass, iw, liw, a, la, lda, ioldps, poselt, keep, keep8, postpone_col_update, etatass, typefile, lafac, monbloc, nextpiv2bewritten, liwfac, myid, iflag, offset_iw, inode)
Definition dfac_front_aux.F:2311
dmumps_fac_front_aux_m::dmumps_fac_i_ldlt
subroutine dmumps_fac_i_ldlt(nfront, nass, inode, ibeg_block, iend_block, iw, liw, a, la, inopv, nnegw, nb22t1w, nbtinyw, det_expw, det_mantw, det_signw, iflag, ioldps, poselt, uu, seuil, keep, keep8, pivsiz, dkeep, pivnul_list, lpn_list, xsize, pp_first2swap_l, pp_lastpanelondisk, pp_lastpivrptrindexfilled, maxfromm, is_maxfromm_avail, pivot_option, iend_blr, inextpiv, ooc_effective_on_front, nvschur, parpiv_t1, lr_activated)
Definition dfac_front_aux.F:1102
dmumps_fac_front_aux_m::dmumps_fac_fr_update_cbrows
subroutine dmumps_fac_fr_update_cbrows(inode, nfront, nass, call_utrsm, a, la, lafac, poselt, iw, liw, ioldps, monbloc, myid, noffw, det_expw, det_mantw, det_signw, liwfac, pp_first2swap_l, pp_first2swap_u, lnextpiv2bewritten, unextpiv2bewritten, pp_lastpivrptrfilled_l, pp_lastpivrptrfilled_u xsize, seuil, uu, dkeep, keep8, keep, iflag, ooc_effective_on_front, nvschur)
Definition dfac_front_aux.F:628
dmumps_ooc
Definition dmumps_ooc.F:14
dmumps_ooc::dmumps_ooc_io_lu_panel
subroutine, public dmumps_ooc_io_lu_panel(strat, typefile, afac, lafac, monbloc, lnextpiv2bewritten, unextpiv2bewritten, iw, liwfac, myid, filesize, ierr, last_call)
Definition dmumps_ooc.F:3032
dmumps_ooc::typef_both_lu
integer, parameter, public typef_both_lu
Definition dmumps_ooc.F:64
mumps_ooc_common
Definition mumps_ooc_common.F:14
mumps_ooc_common::typef_u
integer, public typef_u
Definition mumps_ooc_common.F:18
mumps_ooc_common::strat_try_write
integer, public strat_try_write
Definition mumps_ooc_common.F:45
mumps_ooc_common::typef_l
integer, public typef_l
Definition mumps_ooc_common.F:18