OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
slr_core.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
14C
15C Note: the last routine of this file, xMUMPS_TRUNCATED_RRQR is derived from
16C the LAPACK package, for which BSD 3-clause license applies
17C (see header of the routine).
18 MODULE smumps_lr_core
19 USE mumps_lr_common
20 USE smumps_lr_type
21 USE smumps_lr_stats
22 USE smumps_lr_data_m
23!$ USE OMP_LIB
24 IMPLICIT NONE
25 CONTAINS
26 SUBROUTINE init_lrb(LRB_OUT,K,M,N,ISLR)
27C This routine simply initializes a LR block but does NOT allocate it
28C (allocation occurs somewhere else)
29 TYPE(lrb_type), INTENT(OUT) :: LRB_OUT
30 INTEGER,INTENT(IN) :: K,M,N
31 LOGICAL,INTENT(IN) :: ISLR
32 lrb_out%M = m
33 lrb_out%N = n
34 lrb_out%K = k
35 lrb_out%ISLR = islr
36 NULLIFY(lrb_out%Q)
37 NULLIFY(lrb_out%R)
38 END SUBROUTINE init_lrb
39C
40C
41 SUBROUTINE is_front_blr_candidate(INODE, NIV, NFRONT, NASS,
42 & BLRON, K489,
43 & K490, K491, K492, K20, K60, IDAD, K38,
44 & LRSTATUS, N, LRGROUPS)
45 INTEGER,INTENT(IN) :: INODE, NFRONT, NASS, BLRON, K489, K490,
46 & K491, K492, NIV, K20, K60, IDAD, K38
47 INTEGER,INTENT(OUT):: LRSTATUS
48 INTEGER, INTENT(IN):: N
49 INTEGER, INTENT(IN), OPTIONAL :: LRGROUPS(N)
50C
51C Local variables
52 LOGICAL :: COMPRESS_PANEL, COMPRESS_CB
53 LRSTATUS = 0
54 compress_panel = .false.
55 IF ((blron.NE.0).and.(
56 & ((k492.LT.0).and.inode.EQ.abs(k492))
57 & .or.
58 & ( (k492.GT.0).and.(k491.LE.nfront)
59 & .and.(k490.LE.nass)))) THEN
60 compress_panel = .true.
61C Compression for NASS =1 is useless
62 IF (nass.LE.1) THEN
63 compress_panel =.false.
64 ENDIF
65 IF (present(lrgroups)) THEN
66 IF (lrgroups(inode) .LT. 0) compress_panel = .false.
67 ENDIF
68 ENDIF
69 compress_cb = .false.
70 IF ((blron.NE.0).and.
71 & (k489.GT.0.AND.(k489.NE.2.OR.niv.EQ.2))
72 & .and.(
73 & ((k492.LT.0).and.inode.EQ.abs(k492))
74 & .or.
75 & ((k492.GT.0).AND.(nfront-nass.GT.k491))))
76 & THEN
77 compress_cb = .true.
78 ENDIF
79 IF (.NOT.compress_panel) compress_cb=.false.
80 IF (compress_panel.OR.compress_cb) THEN
81 IF (compress_cb.AND.(.NOT.compress_panel)) THEN
82 lrstatus = 1
83 ELSE IF (compress_panel.AND.(.NOT.compress_cb)) THEN
84 lrstatus = 2
85 ELSE
86 lrstatus = 3
87 ENDIF
88 ELSE
89 lrstatus = 0
90 ENDIF
91C
92C Schur complement cannot be BLR for now
93C
94 IF ( inode .EQ. k20 .AND. k60 .NE. 0 ) THEN
95 lrstatus = 0
96 ENDIF
97C
98C Do not compress CB of children of root
99C
100 IF ( idad .EQ. k38 .AND. k38 .NE.0 ) THEN
101 compress_cb = .false.
102 IF (lrstatus.GE.2) THEN
103 lrstatus = 2
104 ELSE
105 lrstatus = 0
106 ENDIF
107 ENDIF
108 RETURN
109 END SUBROUTINE is_front_blr_candidate
110 SUBROUTINE alloc_lrb(LRB_OUT,K,M,N,ISLR,IFLAG,IERROR,KEEP8)
111 TYPE(lrb_type), INTENT(OUT) :: LRB_OUT
112 INTEGER,INTENT(IN) :: K,M,N
113 INTEGER,INTENT(INOUT) :: IFLAG, IERROR
114 LOGICAL,INTENT(IN) :: ISLR
115 INTEGER(8) :: KEEP8(150)
116 INTEGER :: MEM, allocok
117 REAL :: ZERO
118 parameter(zero = 0.0d0)
119 lrb_out%M = m
120 lrb_out%N = n
121 lrb_out%K = k
122 lrb_out%ISLR = islr
123 IF ((m.EQ.0).OR.(n.EQ.0)) THEN
124 nullify(lrb_out%Q)
125 nullify(lrb_out%R)
126 RETURN
127 ENDIF
128 IF (islr) THEN
129 IF (k.EQ.0) THEN
130 nullify(lrb_out%Q)
131 nullify(lrb_out%R)
132 ELSE
133 allocate(lrb_out%Q(m,k),lrb_out%R(k,n),stat=allocok)
134 IF (allocok > 0) THEN
135 iflag = -13
136 ierror = k*(m+n)
137 RETURN
138 ENDIF
139 ENDIF
140 ELSE
141 nullify(lrb_out%R)
142 allocate(lrb_out%Q(m,n),stat=allocok)
143 IF (allocok > 0) THEN
144 iflag = -13
145 ierror = m*n
146 RETURN
147 ENDIF
148 ENDIF
149 IF (islr) THEN
150 mem = m*k + n*k
151 ELSE
152 mem = m*n
153 ENDIF
154 CALL mumps_dm_fac_upd_dyn_memcnts(int(mem,8),
155 & .true., keep8, iflag, ierror, .true., .true.)
156 RETURN
157 END SUBROUTINE alloc_lrb
158 SUBROUTINE alloc_lrb_from_acc(ACC_LRB, LRB_OUT, K, M, N, LorU,
159 & IFLAG, IERROR, KEEP8)
160 TYPE(lrb_type), INTENT(IN) :: ACC_LRB
161 TYPE(lrb_type), INTENT(OUT) :: LRB_OUT
162 INTEGER,INTENT(IN) :: K, M, N, LorU
163 INTEGER,INTENT(INOUT) :: IFLAG, IERROR
164 INTEGER(8) :: KEEP8(150)
165 INTEGER :: I
166 IF (loru.EQ.1) THEN
167 CALL alloc_lrb(lrb_out,k,m,n,.true.,iflag,ierror,keep8)
168 IF (iflag.LT.0) RETURN
169 DO i=1,k
170 lrb_out%Q(1:m,i) = acc_lrb%Q(1:m,i)
171 lrb_out%R(i,1:n) = -acc_lrb%R(i,1:n)
172 ENDDO
173 ELSE
174 CALL alloc_lrb(lrb_out,k,n,m,.true.,iflag,ierror,keep8)
175 IF (iflag.LT.0) RETURN
176 DO i=1,k
177 lrb_out%Q(1:n,i) = acc_lrb%R(i,1:n)
178 lrb_out%R(i,1:m) = -acc_lrb%Q(1:m,i)
179 ENDDO
180 ENDIF
181 END SUBROUTINE alloc_lrb_from_acc
182 SUBROUTINE regrouping2(CUT, NPARTSASS, NASS,
183 & NPARTSCB, NCB, IBCKSZ, ONLYCB, K472)
184 INTEGER, INTENT(IN) :: IBCKSZ, NASS, NCB
185 INTEGER, INTENT(INOUT) :: NPARTSCB, NPARTSASS
186 INTEGER, POINTER, DIMENSION(:) :: CUT
187 INTEGER, POINTER, DIMENSION(:) :: NEW_CUT
188 INTEGER :: I, INEW, MINSIZE, NEW_NPARTSASS, allocok
189 LOGICAL :: ONLYCB, TRACE
190 INTEGER, INTENT(IN) :: K472
191 INTEGER :: IBCKSZ2,IFLAG,IERROR
192 ALLOCATE(new_cut(max(npartsass,1)+npartscb+1),stat=allocok)
193 IF (allocok > 0) THEN
194 iflag = -13
195 ierror = max(npartsass,1)+npartscb+1
196 write(*,*) 'Allocation problem in BLR routine REGROUPING2:',
197 & ' not enough memory? memory requested = ' , ierror
198 RETURN
199 ENDIF
200 CALL compute_blr_vcs(k472, ibcksz2, ibcksz, nass)
201 minsize = int(ibcksz2 / 2)
202 new_npartsass = max(npartsass,1)
203 IF (.NOT. onlycb) THEN
204 new_cut(1) = 1
205 inew = 2
206 i = 2
207 DO WHILE (i .LE. npartsass + 1)
208 new_cut(inew) = cut(i)
209 trace = .false.
210 IF (new_cut(inew) - new_cut(inew-1) .GT. minsize) THEN
211 inew = inew + 1
212 trace = .true.
213 ENDIF
214 i = i + 1
215 END DO
216 IF (trace) THEN
217 inew = inew - 1
218 ELSE
219 IF (inew .NE. 2) THEN
220 new_cut(inew-1) = new_cut(inew)
221 inew = inew - 1
222 ENDIF
223 ENDIF
224 new_npartsass = inew - 1
225 ENDIF
226 IF (onlycb) THEN
227 DO i=1,max(npartsass,1)+1
228 new_cut(i) = cut(i)
229 ENDDO
230 ENDIF
231 IF (ncb .EQ. 0) GO TO 50
232 inew = new_npartsass+2
233 i = max(npartsass,1) + 2
234 DO WHILE (i .LE. max(npartsass,1) + npartscb + 1)
235 new_cut(inew) = cut(i)
236 trace = .false.
237 IF (new_cut(inew) - new_cut(inew-1) .GT. minsize) THEN
238 inew = inew + 1
239 trace = .true.
240 ENDIF
241 i = i + 1
242 END DO
243 IF (trace) THEN
244 inew = inew - 1
245 ELSE
246 IF (inew .NE. new_npartsass+2) THEN
247 new_cut(inew-1) = new_cut(inew)
248 inew = inew - 1
249 ENDIF
250 ENDIF
251 npartscb = inew - 1 - new_npartsass
252 50 CONTINUE
253 npartsass = new_npartsass
254 DEALLOCATE(cut)
255 ALLOCATE(cut(npartsass+npartscb+1),stat=allocok)
256 IF (allocok > 0) THEN
257 iflag = -13
258 ierror = npartsass+npartscb+1
259 write(*,*) 'Allocation problem in BLR routine REGROUPING2:',
260 & ' not enough memory? memory requested = ' , ierror
261 RETURN
262 ENDIF
263 DO i=1,npartsass+npartscb+1
264 cut(i) = new_cut(i)
265 ENDDO
266 DEALLOCATE(new_cut)
267 END SUBROUTINE regrouping2
268 SUBROUTINE smumps_lrtrsm(A, LA, POSELT_LOCAL, NFRONT, LDA, LRB,
269 & NIV, SYM, LorU, IW, OFFSET_IW)
270C -----------
271C Parameters
272C -----------
273 INTEGER(8), intent(in) :: LA
274 INTEGER, intent(in) :: NFRONT, NIV, SYM, LorU, LDA
275 INTEGER(8), intent(in) :: POSELT_LOCAL
276 REAL, intent(inout) :: A(LA)
277 TYPE(LRB_TYPE), intent(inout) :: LRB
278 INTEGER, OPTIONAL:: OFFSET_IW
279 INTEGER, OPTIONAL :: IW(*)
280C -----------
281C Local variables
282C -----------
283 INTEGER(8) :: DPOS, POSPV1, POSPV2, OFFDAG
284 INTEGER :: M, N, I, J
285 REAL, POINTER :: LR_BLOCK_PTR(:,:)
286 REAL :: ONE, MONE, ZERO
287 REAL :: MULT1, MULT2, A11, DETPIV, A22, A12
288 PARAMETER (ONE = 1.0e0, mone=-1.0e0)
289 parameter(zero=0.0e0)
290 n = lrb%N
291 IF (lrb%ISLR) THEN
292 m = lrb%K
293 lr_block_ptr => lrb%R
294 ELSE
295 m = lrb%M
296 lr_block_ptr => lrb%Q
297 END IF
298 IF (m.NE.0) THEN
299C Why is it Right, Lower, Tranpose?
300C Because A is stored by rows
301C but BLR_L is stored by columns
302 IF (sym.EQ.0.AND.loru.EQ.0) THEN
303 CALL strsm('R', 'L', 'T', 'N', m, n, one,
304 & a(poselt_local), nfront,
305 & lr_block_ptr(1,1), m)
306 ELSE
307 CALL strsm('R', 'U', 'N', 'U', m, n, one,
308 & a(poselt_local), lda,
309 & lr_block_ptr(1,1), m)
310 IF (loru.EQ.0) THEN
311C Now apply D scaling
312 IF (.NOT.present(offset_iw)) THEN
313 write(*,*) 'Internal error in ',
314 & 'SMUMPS_LRTRSM'
315 CALL mumps_abort()
316 ENDIF
317 dpos = poselt_local
318 i = 1
319 DO
320 IF(i .GT. n) EXIT
321 IF(iw(offset_iw+i-1) .GT. 0) THEN
322C 1x1 pivot
323 a11 = one/a(dpos)
324 CALL sscal(m, a11, lr_block_ptr(1,i), 1)
325 dpos = dpos + int(lda + 1,8)
326 i = i+1
327 ELSE
328C 2x2 pivot
329 pospv1 = dpos
330 pospv2 = dpos+ int(lda + 1,8)
331 offdag = pospv1+1_8
332 a11 = a(pospv1)
333 a22 = a(pospv2)
334 a12 = a(offdag)
335 detpiv = a11*a22 - a12**2
336 a22 = a11/detpiv
337 a11 = a(pospv2)/detpiv
338 a12 = -a12/detpiv
339 DO j = 1,m
340 mult1 = a11*lr_block_ptr(j,i)+a12*lr_block_ptr(j,i+1)
341 mult2 = a12*lr_block_ptr(j,i)+a22*lr_block_ptr(j,i+1)
342 lr_block_ptr(j,i) = mult1
343 lr_block_ptr(j,i+1) = mult2
344 ENDDO
345 dpos = pospv2 + int(lda + 1,8)
346 i = i+2
347 ENDIF
348 ENDDO
349 ENDIF
350 ENDIF
351 ENDIF
352 CALL upd_flop_trsm(lrb, loru)
353 END SUBROUTINE smumps_lrtrsm
354 SUBROUTINE smumps_lrgemm_scaling(LRB, SCALED, A, LA, DIAG,
355 & LD_DIAG, IW2, POSELTT, NFRONT, BLOCK, MAXI_CLUSTER)
356C This routine does the scaling (for the symmetric case) before
357C computing the LR product (done in SMUMPS_LRGEMM4)
358 TYPE(lrb_type),INTENT(IN) :: lrb
359 INTEGER(8), intent(in) :: la
360 REAL, intent(inout) :: a(la)
361 REAL, intent(inout), DIMENSION(:,:) :: SCALED
362 INTEGER,INTENT(IN) :: LD_DIAG, NFRONT, IW2(*)
363 INTEGER(8), INTENT(IN) :: POSELTT
364 REAL, INTENT(IN), OPTIONAL :: DIAG(*)
365 INTEGER, INTENT(IN) :: MAXI_CLUSTER
366 REAL, intent(inout) :: BLOCK(MAXI_CLUSTER)
367 INTEGER :: J, NROWS
368 REAL :: PIV1, PIV2, OFFDIAG
369 IF (lrb%ISLR) THEN
370 nrows = lrb%K
371 ELSE
372 nrows = lrb%M
373 ENDIF
374 j = 1
375 DO WHILE (j <= lrb%N)
376 IF (iw2(j) > 0) THEN
377 scaled(1:nrows,j) = diag(1+ld_diag*(j-1)+j-1)
378 & * scaled(1:nrows,j)
379 j = j+1
380 ELSE !2x2 pivot
381 piv1 = diag(1+ld_diag*(j-1)+j-1)
382 piv2 = diag(1+ld_diag*j+j)
383 offdiag = diag(1+ld_diag*(j-1)+j)
384 block(1:nrows) = scaled(1:nrows,j)
385 scaled(1:nrows,j) = piv1 * scaled(1:nrows,j)
386 & + offdiag * scaled(1:nrows,j+1)
387 scaled(1:nrows,j+1) = offdiag * block(1:nrows)
388 & + piv2 * scaled(1:nrows,j+1)
389 j=j+2
390 ENDIF
391 END DO
392 END SUBROUTINE smumps_lrgemm_scaling
393 SUBROUTINE smumps_lrgemm4(ALPHA,
394 & LRB1, LRB2, BETA,
395 & A, LA, POSELTT, NFRONT, SYM,
396 & IFLAG, IERROR,
397 & MIDBLK_COMPRESS, TOLEPS, TOL_OPT, KPERCENT,
398 & RANK, BUILDQ,
399 & LUA_ACTIVATED,
400C Start of OPTIONAL arguments
401 & LorU,
402 & LRB3, MAXI_RANK,
403 & MAXI_CLUSTER,
404 & DIAG, LD_DIAG, IW2, BLOCK
405 & )
406C
407CC
408 TYPE(lrb_type),INTENT(IN) :: LRB1,LRB2
409 INTEGER(8), intent(in) :: LA
410 REAL, intent(inout) :: A(LA)
411 INTEGER,INTENT(IN) :: NFRONT, SYM, TOL_OPT
412 INTEGER,INTENT(INOUT) :: IFLAG, IERROR
413 INTEGER(8), INTENT(IN) :: POSELTT
414 REAL, INTENT(IN), OPTIONAL :: DIAG(*)
415 INTEGER,INTENT(IN), OPTIONAL :: LD_DIAG, IW2(*)
416 INTEGER,intent(in) :: MIDBLK_COMPRESS, KPERCENT
417 REAL, intent(in) :: TOLEPS
418 REAL :: ALPHA,BETA
419 LOGICAL, INTENT(OUT) :: BUILDQ
420 REAL, intent(inout), OPTIONAL :: BLOCK(*)
421 INTEGER, INTENT(IN), OPTIONAL :: LorU
422 LOGICAL, INTENT(IN) :: LUA_ACTIVATED
423 INTEGER, INTENT(IN), OPTIONAL :: MAXI_CLUSTER
424 INTEGER, INTENT(IN), OPTIONAL :: MAXI_RANK
425 TYPE(LRB_TYPE), INTENT(INOUT), OPTIONAL :: LRB3
426 REAL, POINTER, DIMENSION(:,:) :: XY_YZ
427 REAL, ALLOCATABLE, TARGET, DIMENSION(:,:) :: XQ, R_Y
428 REAL, POINTER, DIMENSION(:,:) :: X, Y, Y1, Y2, Z
429 CHARACTER(len=1) :: SIDE, TRANSY
430 INTEGER :: K_XY, K_YZ, LDY, LDY1, LDY2, K_Y
431 INTEGER :: LDXY_YZ, SAVE_K
432 INTEGER :: I, J, RANK, MAXRANK, INFO, LWORK
433 REAL, ALLOCATABLE :: RWORK_RRQR(:)
434 REAL, ALLOCATABLE :: WORK_RRQR(:), TAU_RRQR(:),
435 & y_rrqr(:,:)
436 INTEGER, ALLOCATABLE :: JPVT_RRQR(:)
437 INTEGER :: allocok, MREQ
438 REAL, EXTERNAL ::snrm2
439 REAL :: ONE, MONE, ZERO
440 parameter(one = 1.0e0, mone=-1.0e0)
441 parameter(zero=0.0e0)
442 IF (lrb1%M.EQ.0) THEN
443 RETURN
444 ENDIF
445 IF (lrb2%M.EQ.0) THEN
446 ENDIF
447 rank = 0
448 buildq = .false.
449 IF (lrb1%ISLR.AND.lrb2%ISLR) THEN
450 IF ((lrb1%K.EQ.0).OR.(lrb2%K.EQ.0)) THEN
451 GOTO 1200
452 ENDIF
453 allocate(y(lrb1%K,lrb2%K),stat=allocok)
454 IF (allocok > 0) THEN
455 mreq = lrb1%K*lrb2%K
456 GOTO 1570
457 ENDIF
458 x => lrb1%Q
459 k_y = lrb1%N
460 IF (sym .EQ. 0) THEN
461 y1 => lrb1%R
462 ELSE
463 allocate(y1(lrb1%K,lrb1%N),stat=allocok)
464 IF (allocok > 0) THEN
465 mreq = lrb1%K*lrb1%N
466 GOTO 1570
467 ENDIF
468 DO j=1,lrb1%N
469 DO i=1,lrb1%K
470 y1(i,j) = lrb1%R(i,j)
471 ENDDO
472 ENDDO
473 CALL smumps_lrgemm_scaling(lrb1, y1, a, la, diag,
474 & ld_diag, iw2, poseltt, nfront, block,
475 & maxi_cluster)
476 ENDIF
477 ldy1 = lrb1%K
478 z => lrb2%Q
479 y2 => lrb2%R
480 ldy2 = lrb2%K
481 CALL sgemm('N', 'T', lrb1%K, lrb2%K, k_y, one,
482 & y1(1,1), ldy1, y2(1,1), ldy2, zero, y(1,1), lrb1%K )
483 IF (midblk_compress.GE.1) THEN
484 lwork = lrb2%K*(lrb2%K+1)
485 allocate(y_rrqr(lrb1%K,lrb2%K),
486 & work_rrqr(lwork), rwork_rrqr(2*lrb2%K),
487 & tau_rrqr(min(lrb1%K,lrb2%K)),
488 & jpvt_rrqr(lrb2%K),stat=allocok)
489 IF (allocok > 0) THEN
490 mreq = lrb1%K*lrb2%K + lwork + 2*lrb2%K +
491 & min(lrb1%K,lrb2%K) + lrb2%K
492 GOTO 1570
493 ENDIF
494 DO j=1,lrb2%K
495 DO i=1,lrb1%K
496 y_rrqr(i,j) = y(i,j)
497 ENDDO
498 ENDDO
499 maxrank = min(lrb1%K, lrb2%K)-1
500 maxrank = max(1, int((maxrank*kpercent/100)))
501 jpvt_rrqr = 0
502 CALL smumps_truncated_rrqr(lrb1%K, lrb2%K, y_rrqr(1,1),
503 & lrb1%K, jpvt_rrqr, tau_rrqr, work_rrqr,
504 & lrb2%K, rwork_rrqr, toleps, tol_opt, rank,
505 & maxrank, info,
506 & buildq)
507 IF (rank.GT.maxrank) THEN
508 deallocate(y_rrqr, work_rrqr, rwork_rrqr, tau_rrqr,
509 & jpvt_rrqr)
510 buildq = .false.
511 ELSE
512 buildq = .true.
513 ENDIF
514 IF (buildq) THEN
515 IF (rank.EQ.0) THEN
516 deallocate(y_rrqr, work_rrqr, rwork_rrqr, tau_rrqr,
517 & jpvt_rrqr)
518 deallocate(y)
519 nullify(y)
520C GOTO 1580 not ok because BUILDQ .EQV. true
521C would try to free XQ and R_Y that are not allocated
522C in that case. So we free Y1 now if it was allocated.
523 IF (sym .NE. 0) deallocate(y1)
524 GOTO 1200
525 ELSE
526 allocate(xq(lrb1%M,rank), r_y(rank,lrb2%K),
527 & stat=allocok)
528 IF (allocok > 0) THEN
529 mreq = lrb1%M*rank + rank*lrb2%K
530 GOTO 1570
531 ENDIF
532 DO j=1, lrb2%K
533 r_y(1:min(rank,j),jpvt_rrqr(j)) =
534 & y_rrqr(1:min(rank,j),j)
535 IF(j.LT.rank) r_y(min(rank,j)+1:
536 & rank,jpvt_rrqr(j))= zero
537 END DO
538C LWORK=LRB2%K*(LRB2%K+1), with LRB2%K>RANK
539C large enough for sorgqr
540 CALL sorgqr
541 & (lrb1%K, rank, rank, y_rrqr(1,1),
542 & lrb1%K, tau_rrqr(1),
543 & work_rrqr(1), lwork, info )
544 CALL sgemm('N', 'N', lrb1%M, rank, lrb1%K, one,
545 & x(1,1), lrb1%M, y_rrqr(1,1), lrb1%K, zero,
546 & xq(1,1), lrb1%M)
547 deallocate(y_rrqr, work_rrqr, rwork_rrqr, tau_rrqr,
548 & jpvt_rrqr)
549 nullify(x)
550 x => xq
551 k_xy = rank
552 deallocate(y)
553 nullify(y)
554 y => r_y
555 ldy = rank
556 k_yz = lrb2%K
557 transy = 'N'
558 side = 'R'
559 ENDIF
560 ENDIF
561 ENDIF
562 IF (.NOT.buildq) THEN
563 ldy = lrb1%K
564 k_xy = lrb1%K
565 k_yz = lrb2%K
566 transy = 'N'
567 IF (lrb1%K .GE. lrb2%K) THEN
568 side = 'L'
569 ELSE
570 side = 'R'
571 ENDIF
572 ENDIF
573 ENDIF
574 IF (lrb1%ISLR.AND.(.NOT.lrb2%ISLR)) THEN
575 IF (lrb1%K.EQ.0) THEN
576 GOTO 1200
577 ENDIF
578 side = 'R'
579 k_xy = lrb1%K
580 transy = 'N'
581 z => lrb2%Q
582 x => lrb1%Q
583 ldy = lrb1%K
584 IF (sym .EQ. 0) THEN
585 y => lrb1%R
586 ELSE
587 allocate(y(lrb1%K,lrb1%N),stat=allocok)
588 IF (allocok > 0) THEN
589 mreq = lrb1%K*lrb1%N
590 GOTO 1570
591 ENDIF
592 DO j=1,lrb1%N
593 DO i=1,lrb1%K
594 y(i,j) = lrb1%R(i,j)
595 ENDDO
596 ENDDO
597 CALL smumps_lrgemm_scaling(lrb1, y, a, la, diag,
598 & ld_diag, iw2, poseltt, nfront, block,
599 & maxi_cluster)
600 ENDIF
601 k_yz = lrb2%N
602 ENDIF
603 IF ((.NOT.lrb1%ISLR).AND.lrb2%ISLR) THEN
604 IF (lrb2%K.EQ.0) THEN
605 GOTO 1200
606 ENDIF
607 side = 'L'
608 k_yz = lrb2%K
609 x => lrb1%Q
610 transy = 'T'
611 k_xy = lrb1%N
612 IF (sym .EQ. 0) THEN
613 y => lrb2%R
614 ELSE
615 allocate(y(lrb2%K,lrb2%N),stat=allocok)
616 IF (allocok > 0) THEN
617 mreq = lrb2%K*lrb2%N
618 GOTO 1570
619 ENDIF
620 DO j=1,lrb2%N
621 DO i=1,lrb2%K
622 y(i,j) = lrb2%R(i,j)
623 ENDDO
624 ENDDO
625 CALL smumps_lrgemm_scaling(lrb2, y, a, la, diag,
626 & ld_diag, iw2, poseltt, nfront, block,
627 & maxi_cluster)
628 ENDIF
629 ldy = lrb2%K
630 z => lrb2%Q
631 ENDIF
632 IF ((.NOT.lrb1%ISLR).AND.(.NOT.lrb2%ISLR)) THEN
633 IF (sym .EQ. 0) THEN
634 x => lrb1%Q
635 ELSE
636 allocate(x(lrb1%M,lrb1%N),stat=allocok)
637 IF (allocok > 0) THEN
638 mreq = lrb1%M*lrb1%N
639 GOTO 1570
640 ENDIF
641 DO j=1,lrb1%N
642 DO i=1,lrb1%M
643 x(i,j) = lrb1%Q(i,j)
644 ENDDO
645 ENDDO
646 CALL smumps_lrgemm_scaling(lrb1, x, a, la, diag,
647 & ld_diag, iw2, poseltt, nfront, block,
648 & maxi_cluster)
649 ENDIF
650 side = 'N'
651 z => lrb2%Q
652 k_xy = lrb1%N
653 ENDIF
654 IF (lua_activated) THEN
655 save_k = lrb3%K
656 IF (side == 'L') THEN
657 lrb3%K = lrb3%K+k_yz
658 ELSEIF (side == 'R') THEN
659 lrb3%K = lrb3%K+k_xy
660 ENDIF
661 ENDIF
662 IF (side == 'L') THEN ! LEFT: XY_YZ = X*Y; A = XY_YZ*Z
663 IF (.NOT.lua_activated) THEN
664 allocate(xy_yz(lrb1%M,k_yz),stat=allocok)
665 IF (allocok > 0) THEN
666 mreq = lrb1%M*k_yz
667 GOTO 1570
668 ENDIF
669 ldxy_yz = lrb1%M
670 ELSE
671 IF (save_k+k_yz.GT.maxi_rank) THEN
672 write(*,*) 'Internal error in SMUMPS_LRGEMM4 1a',
673 & 'K_ACC+K_CUR>K_MAX:',save_k,k_yz,maxi_rank
674 CALL mumps_abort()
675 ENDIF
676 IF (lrb3%M.NE.lrb1%M) THEN
677 write(*,*) 'Internal error in SMUMPS_LRGEMM4 1b',
678 & 'LRB1%M =/= LRB3%M',lrb1%M,lrb3%M
679 CALL mumps_abort()
680 ENDIF
681 xy_yz => lrb3%Q(1:lrb1%M,save_k+1:save_k+k_yz)
682 ldxy_yz = maxi_cluster
683 DO i=1,k_yz
684 lrb3%R(save_k+i,1:lrb2%M) = z(1:lrb2%M,i)
685 ENDDO
686 ENDIF
687 CALL sgemm('N', transy, lrb1%M, k_yz, k_xy, one,
688 & x(1,1), lrb1%M, y(1,1), ldy, zero, xy_yz(1,1),
689 & ldxy_yz)
690 IF (.NOT.lua_activated) THEN
691 CALL sgemm('N', 'T', lrb1%M, lrb2%M, k_yz, alpha,
692 & xy_yz(1,1), lrb1%M, z(1,1), lrb2%M, beta,
693 & a(poseltt), nfront)
694 deallocate(xy_yz)
695 ENDIF
696 ELSEIF (side == 'R') THEN ! RIGHT: XY_YZ = Y*Z; A = X*XY_YZ
697 IF (.NOT.lua_activated) THEN
698 allocate(xy_yz(k_xy,lrb2%M),stat=allocok)
699 IF (allocok > 0) THEN
700 mreq = k_xy*lrb2%M
701 GOTO 1570
702 ENDIF
703 ldxy_yz = k_xy
704 ELSE
705 IF (save_k+k_xy.GT.maxi_rank) THEN
706 write(*,*) 'Internal error in SMUMPS_LRGEMM4 2a',
707 & 'K_ACC+K_CUR>K_MAX:',save_k,k_xy,maxi_rank
708 CALL mumps_abort()
709 ENDIF
710 IF (lrb3%N.NE.lrb2%M) THEN
711 write(*,*) 'Internal error in SMUMPS_LRGEMM4 2b',
712 & 'LRB2%M =/= LRB3%N',lrb2%M,lrb3%N
713 CALL mumps_abort()
714 ENDIF
715 xy_yz => lrb3%R(save_k+1:save_k+k_xy,1:lrb2%M)
716 ldxy_yz = maxi_rank
717 DO i=1,k_xy
718 lrb3%Q(1:lrb1%M,save_k+i) = x(1:lrb1%M,i)
719 ENDDO
720 ENDIF
721 CALL sgemm(transy, 'T', k_xy, lrb2%M, k_yz, one,
722 & y(1,1), ldy, z(1,1), lrb2%M, zero, xy_yz(1,1),
723 & ldxy_yz)
724 IF (.NOT.lua_activated) THEN
725 CALL sgemm('N', 'N', lrb1%M, lrb2%M, k_xy, alpha,
726 & x(1,1), lrb1%M, xy_yz(1,1), k_xy, beta, a(poseltt),
727 & nfront)
728 deallocate(xy_yz)
729 ENDIF
730 ELSE ! SIDE == 'N' : NONE; A = X*Z
731 CALL sgemm('N', 'T', lrb1%M, lrb2%M, k_xy, alpha,
732 & x(1,1), lrb1%M, z(1,1), lrb2%M, beta, a(poseltt),
733 & nfront)
734 ENDIF
735 GOTO 1580
736 1570 CONTINUE
737C Alloc NOT ok!!
738 iflag = -13
739 ierror = mreq
740 RETURN
741 1580 CONTINUE
742C Alloc ok!!
743 IF ((.NOT.lrb1%ISLR).AND.(.NOT.lrb2%ISLR)) THEN
744 IF (sym .NE. 0) deallocate(x)
745 ELSEIF ((.NOT.lrb1%ISLR).AND.lrb2%ISLR) THEN
746 IF (sym .NE. 0) deallocate(y)
747 ELSEIF (lrb1%ISLR.AND.(.NOT.lrb2%ISLR)) THEN
748 IF (sym .NE. 0) deallocate(y)
749 ELSE
750 IF (sym .NE. 0) deallocate(y1)
751 IF ((midblk_compress.GE.1).AND.buildq) THEN
752 deallocate(xq)
753 deallocate(r_y)
754 ELSE
755 deallocate(y)
756 ENDIF
757 ENDIF
758 1200 CONTINUE
759 END SUBROUTINE smumps_lrgemm4
760 SUBROUTINE smumps_decompress_acc(ACC_LRB, MAXI_CLUSTER,
761 & MAXI_RANK, A, LA, POSELTT, NFRONT, NIV, LorU,
762 & COUNT_FLOPS)
763 TYPE(lrb_type),INTENT(INOUT) :: ACC_LRB
764 INTEGER(8), intent(in) :: LA
765 REAL, intent(inout) :: A(LA)
766 INTEGER,INTENT(IN) :: NFRONT, NIV
767 INTEGER,INTENT(IN) :: MAXI_CLUSTER, MAXI_RANK
768 INTEGER(8), INTENT(IN) :: POSELTT
769 LOGICAL, OPTIONAL :: COUNT_FLOPS
770 LOGICAL :: COUNT_FLOPS_LOC
771 INTEGER :: LorU
772 REAL :: ONE, MONE, ZERO
773 parameter(one = 1.0e0, mone=-1.0e0)
774 parameter(zero=0.0e0)
775 IF (present(count_flops)) THEN
776 count_flops_loc=count_flops
777 ELSE
778 count_flops_loc=.true.
779 ENDIF
780 CALL sgemm('N', 'N', acc_lrb%M, acc_lrb%N, acc_lrb%K,
781 & mone, acc_lrb%Q(1,1), maxi_cluster, acc_lrb%R(1,1),
782 & maxi_rank, one, a(poseltt), nfront)
783 acc_lrb%K = 0
784 END SUBROUTINE smumps_decompress_acc
785 SUBROUTINE smumps_compress_fr_updates(ACC_LRB, MAXI_CLUSTER,
786 & MAXI_RANK, A, LA, POSELTT, NFRONT, NIV,
787 & TOLEPS, TOL_OPT, KPERCENT, BUILDQ, LorU, CB_COMPRESS)
788 TYPE(lrb_type),INTENT(INOUT) :: ACC_LRB
789 INTEGER(8), intent(in) :: LA
790 REAL, intent(inout) :: A(LA)
791 INTEGER,INTENT(IN) :: NFRONT, NIV, LorU, TOL_OPT
792 INTEGER,INTENT(IN) :: MAXI_CLUSTER, MAXI_RANK, KPERCENT
793 INTEGER(8), INTENT(IN) :: POSELTT
794 REAL, intent(in) :: TOLEPS
795 LOGICAL, INTENT(OUT) :: BUILDQ
796 LOGICAL, INTENT(IN) :: CB_COMPRESS
797 REAL, ALLOCATABLE :: RWORK_RRQR(:)
798 REAL, ALLOCATABLE :: WORK_RRQR(:), TAU_RRQR(:)
799 INTEGER, ALLOCATABLE :: JPVT_RRQR(:)
800 INTEGER :: INFO, RANK, MAXRANK, LWORK
801 INTEGER :: I, J, M, N
802 INTEGER :: allocok, MREQ
803 REAL :: ONE, MONE, ZERO
804 PARAMETER (ONE = 1.0e0, mone=-1.0e0)
805 parameter(zero=0.0e0)
806 m = acc_lrb%M
807 n = acc_lrb%N
808 maxrank = floor(real(m*n)/real(m+n))
809 maxrank = max(1, int((maxrank*kpercent/100)))
810 lwork = n*(n+1)
811 allocate(work_rrqr(lwork), rwork_rrqr(2*n),
812 & tau_rrqr(n),
813 & jpvt_rrqr(n), stat=allocok)
814 IF (allocok > 0) THEN
815 mreq = lwork +4 *n
816 GOTO 100
817 ENDIF
818 DO i=1,n
819 acc_lrb%Q(1:m,i)=
820 & - a(poseltt+int(i-1,8)*int(nfront,8) :
821 & poseltt+int(i-1,8)*int(nfront,8) + int(m-1,8) )
822 END DO
823 jpvt_rrqr = 0
824 CALL smumps_truncated_rrqr(m, n, acc_lrb%Q(1,1),
825 & maxi_cluster, jpvt_rrqr(1), tau_rrqr(1),
826 & work_rrqr(1),
827 & n, rwork_rrqr(1), toleps, tol_opt,
828 & rank, maxrank, info,
829 & buildq)
830 IF (buildq) THEN
831 DO j=1, n
832 acc_lrb%R(1:min(rank,j),jpvt_rrqr(j)) =
833 & acc_lrb%Q(1:min(rank,j),j)
834 IF(j.LT.rank) acc_lrb%R(min(rank,j)+1:
835 & rank,jpvt_rrqr(j))= zero
836 END DO
837 CALL sorgqr
838 & (m, rank, rank, acc_lrb%Q(1,1),
839 & maxi_cluster, tau_rrqr(1),
840 & work_rrqr(1), lwork, info )
841 DO i=1,n
842 a( poseltt+int(i-1,8)*int(nfront,8) :
843 & poseltt+int(i-1,8)*int(nfront,8)+int(m-1,8) ) = zero
844 END DO
845 acc_lrb%K = rank
846 CALL upd_flop_compress(acc_lrb, cb_compress=cb_compress)
847 ELSE
848 acc_lrb%K = rank
849 acc_lrb%ISLR = .false.
850 CALL upd_flop_compress(acc_lrb, cb_compress=cb_compress)
851 acc_lrb%ISLR = .true.
852 acc_lrb%K = 0
853 ENDIF
854 deallocate(jpvt_rrqr, tau_rrqr, work_rrqr, rwork_rrqr)
855 RETURN
856 100 CONTINUE
857C Alloc NOT ok!!
858 write(*,*) 'Allocation problem in BLR routine
859 & SMUMPS_COMPRESS_FR_UPDATES: ',
860 & 'not enough memory? memory requested = ' , mreq
861 CALL mumps_abort()
862 RETURN
863 END SUBROUTINE smumps_compress_fr_updates
864 SUBROUTINE smumps_recompress_acc(ACC_LRB, MAXI_CLUSTER,
865 & MAXI_RANK, A, LA, POSELTT, NFRONT, NIV,
866 & MIDBLK_COMPRESS, TOLEPS, TOL_OPT, KPERCENT_RMB,
867 & KPERCENT_LUA, NEW_ACC_RANK)
868 TYPE(lrb_type),INTENT(INOUT) :: acc_lrb
869 INTEGER(8), intent(in) :: la
870 REAL, intent(inout) :: a(la)
871 INTEGER,INTENT(IN) :: NFRONT, NIV, TOL_OPT
872 INTEGER :: IFLAG, IERROR
873 INTEGER,INTENT(IN) :: MAXI_CLUSTER, MAXI_RANK, KPERCENT_LUA
874 INTEGER,INTENT(INOUT) :: NEW_ACC_RANK
875 INTEGER(8), INTENT(IN) :: POSELTT
876 INTEGER,intent(in) :: MIDBLK_COMPRESS, KPERCENT_RMB
877 REAL, intent(in) :: TOLEPS
878 REAL, ALLOCATABLE :: RWORK_RRQR(:)
879 REAL, ALLOCATABLE :: WORK_RRQR(:), TAU_RRQR(:)
880 REAL, ALLOCATABLE, DIMENSION(:,:), TARGET :: Q1, R1,
881 & q2, r2
882 INTEGER, ALLOCATABLE :: JPVT_RRQR(:)
883 TYPE(LRB_TYPE) :: LRB1, LRB2
884 INTEGER :: INFO, RANK1, RANK2, RANK, MAXRANK, LWORK
885 LOGICAL :: BUILDQ, BUILDQ1, BUILDQ2, SKIP1, SKIP2
886 INTEGER :: I, J, M, N, K
887 INTEGER :: allocok, MREQ
888 REAL :: ONE, MONE, ZERO
889 parameter(one = 1.0e0, mone=-1.0e0)
890 parameter(zero=0.0e0)
891 skip1 = .false.
892 skip2 = .false.
893 skip1 = .true.
894 1500 CONTINUE
895 m = acc_lrb%M
896 n = acc_lrb%N
897 k = acc_lrb%K
898 maxrank = k-1
899 maxrank = max(1, int((maxrank*kpercent_lua/100)))
900 lwork = k*(k+1)
901 IF (.false.) THEN
902 CALL smumps_recompress_acc_v2(acc_lrb,
903 & maxi_cluster, maxi_rank, a, la, poseltt,
904 & nfront, niv, midblk_compress, toleps,
905 & tol_opt, kpercent_rmb, kpercent_lua,
906 & new_acc_rank)
907 k = acc_lrb%K
908 maxrank = k-1
909 maxrank = max (1, int((maxrank*kpercent_lua/100)))
910 lwork = k*(k+1)
911 skip1 = .true.
912 skip2 = k.EQ.0
913 ENDIF
914 IF (skip1.AND.skip2) GOTO 1600
915 allocate(q1(m,k), q2(n,k),
916 & work_rrqr(lwork),
917 & rwork_rrqr(2*k),
918 & tau_rrqr(k),
919 & jpvt_rrqr(k), stat=allocok)
920 IF (allocok > 0) THEN
921 mreq = lwork + m*n + n*k+ 4 * k
922 GOTO 100
923 ENDIF
924 IF (skip1) THEN
925 buildq1 = .false.
926 ELSE
927 DO j=1,k
928 DO i=1,m
929 q1(i,j) = acc_lrb%Q(i,j)
930 ENDDO
931 ENDDO
932 jpvt_rrqr = 0
933 CALL smumps_truncated_rrqr(m, k, q1(1,1),
934 & m, jpvt_rrqr(1), tau_rrqr(1), work_rrqr(1),
935 & k, rwork_rrqr(1), toleps, tol_opt, rank1,
936 & maxrank, info,
937 & buildq1)
938 ENDIF
939 IF (buildq1) THEN
940 allocate(r1(rank1,k), stat=allocok)
941 IF (allocok > 0) THEN
942 mreq = rank1*k
943 GOTO 100
944 ENDIF
945 DO j=1, k
946 r1(1:min(rank1,j),jpvt_rrqr(j)) =
947 & q1(1:min(rank1,j),j)
948 IF(j.LT.rank1) r1(min(rank1,j)+1:
949 & rank1,jpvt_rrqr(j))= zero
950 END DO
951 CALL sorgqr
952 & (m, rank1, rank1, q1(1,1),
953 & m, tau_rrqr(1),
954 & work_rrqr(1), lwork, info )
955 ENDIF
956 IF (skip2) THEN
957 buildq2 = .false.
958 ELSE
959 DO j=1,k
960 DO i=1,n
961 q2(i,j) = acc_lrb%R(j,i)
962 ENDDO
963 ENDDO
964 jpvt_rrqr = 0
965 CALL smumps_truncated_rrqr(n, k, q2(1,1),
966 & n, jpvt_rrqr(1), tau_rrqr(1), work_rrqr(1),
967 & k, rwork_rrqr(1), toleps, tol_opt,
968 & rank2, maxrank, info,
969 & buildq2)
970 ENDIF
971 IF (buildq2) THEN
972 allocate(r2(rank2,k), stat=allocok)
973 IF (allocok > 0) THEN
974 mreq = rank2*k
975 GOTO 100
976 ENDIF
977 DO j=1, k
978 r2(1:min(rank2,j),jpvt_rrqr(j)) =
979 & q2(1:min(rank2,j),j)
980 IF(j.LT.rank2) r2(min(rank2,j)+1:
981 & rank2,jpvt_rrqr(j))= zero
982 END DO
983 CALL sorgqr
984 & (n, rank2, rank2, q2(1,1),
985 & n, tau_rrqr(1),
986 & work_rrqr(1), lwork, info )
987 ENDIF
988 CALL init_lrb(lrb1,rank1,m,k,buildq1)
989 CALL init_lrb(lrb2,rank2,n,k,buildq2)
990 IF (buildq1.OR.buildq2) THEN
991 IF (buildq1) THEN
992 lrb1%R => r1
993 ELSE
994 DO j=1,k
995 DO i=1,m
996 q1(i,j) = acc_lrb%Q(i,j)
997 ENDDO
998 ENDDO
999 ENDIF
1000 lrb1%Q => q1
1001 IF (buildq2) THEN
1002 lrb2%R => r2
1003 ELSE
1004 DO j=1,k
1005 DO i=1,n
1006 q2(i,j) = acc_lrb%R(j,i)
1007 ENDDO
1008 ENDDO
1009 ENDIF
1010 lrb2%Q => q2
1011 acc_lrb%K = 0
1012 CALL smumps_lrgemm4(mone, lrb1, lrb2, one,
1013 & a, la, poseltt, nfront, 0, iflag, ierror,
1014 & midblk_compress-1, toleps, tol_opt,
1015 & kpercent_rmb,
1016 & rank, buildq, .true., lrb3=acc_lrb,
1017 & maxi_rank=maxi_rank, maxi_cluster=maxi_cluster)
1018 IF (iflag.LT.0) GOTO 100
1019 CALL upd_flop_update(lrb1, lrb2,
1020 & midblk_compress-1, rank, buildq,
1021 & .true., .false., rec_acc=.true.)
1022 ENDIF
1023 IF (.NOT. skip1)
1024 & CALL upd_flop_compress(lrb1, rec_acc=.true.)
1025 IF (.NOT. skip2)
1026 & CALL upd_flop_compress(lrb2, rec_acc=.true.)
1027 deallocate(q1,q2)
1028 IF (buildq1) deallocate(r1)
1029 IF (buildq2) deallocate(r2)
1030 deallocate(jpvt_rrqr, tau_rrqr, work_rrqr, rwork_rrqr)
1031 IF (skip1.AND.(rank2.GT.0)) THEN
1032 skip1 = .false.
1033 skip2 = .true.
1034 GOTO 1500
1035 ENDIF
1036 1600 CONTINUE
1037 new_acc_rank = 0
1038 RETURN
1039 100 CONTINUE
1040C Alloc NOT ok!!
1041 write(*,*) 'Allocation problem in BLR routine
1042 & SMUMPS_RECOMPRESS_ACC: ',
1043 & 'not enough memory? memory requested = ' , mreq
1044 CALL mumps_abort()
1045 RETURN
1046 END SUBROUTINE smumps_recompress_acc
1047 RECURSIVE SUBROUTINE smumps_recompress_acc_narytree(
1048 & ACC_LRB, MAXI_CLUSTER, MAXI_RANK, A, LA, POSELTT, KEEP8,
1049 & NFRONT, NIV, MIDBLK_COMPRESS, TOLEPS, TOL_OPT,
1050 & KPERCENT_RMB,
1051 & KPERCENT_LUA, K478, RANK_LIST, POS_LIST, NB_NODES,
1052 & LEVEL, ACC_TMP)
1053 TYPE(lrb_type),TARGET,INTENT(INOUT) :: acc_lrb
1054 TYPE(lrb_type),TARGET,INTENT(INOUT),OPTIONAL :: acc_tmp
1055 INTEGER(8), intent(in) :: la
1056 REAL, intent(inout) :: a(la)
1057 INTEGER,INTENT(IN) :: nfront, niv, tol_opt
1058 INTEGER,INTENT(IN) :: maxi_cluster, maxi_rank, kpercent_lua
1059 INTEGER(8), INTENT(IN) :: poseltt
1060 INTEGER(8) :: keep8(150)
1061 INTEGER,intent(in) :: midblk_compress, kpercent_rmb
1062 REAL, intent(in) :: toleps
1063 INTEGER,INTENT(IN) :: k478, nb_nodes, level
1064 INTEGER,INTENT(INOUT) :: rank_list(nb_nodes), pos_list(nb_nodes)
1065 TYPE(lrb_type) :: lrb, acc_new
1066 TYPE(lrb_type), POINTER :: lrb_ptr
1067 LOGICAL :: resort
1068 INTEGER :: i, j, m, n, l, node_rank, nary, ioff, imax, curpos
1069 INTEGER :: nb_nodes_new, ktot, new_acc_rank
1070 INTEGER, ALLOCATABLE :: rank_list_new(:), pos_list_new(:)
1071 REAL :: one, mone, zero
1072 parameter(one = 1.0e0, mone=-1.0e0)
1073 parameter(zero=0.0e0)
1074 INTEGER :: allocok
1075 resort = .false.
1076 m = acc_lrb%M
1077 n = acc_lrb%N
1078 nary = -k478
1079 ioff = 0
1080 nb_nodes_new = nb_nodes/nary
1081 IF (nb_nodes_new*nary.NE.nb_nodes) THEN
1082 nb_nodes_new = nb_nodes_new + 1
1083 ENDIF
1084 ALLOCATE(rank_list_new(nb_nodes_new),pos_list_new(nb_nodes_new),
1085 & stat=allocok)
1086 IF (allocok > 0) THEN
1087 write(*,*) 'Allocation error of RANK_LIST_NEW/POS_LIST_NEW ',
1088 & 'in SMUMPS_RECOMPRESS_ACC_NARYTREE'
1089 call mumps_abort()
1090 ENDIF
1091 DO j=1,nb_nodes_new
1092 node_rank = rank_list(ioff+1)
1093 curpos = pos_list(ioff+1)
1094 imax = min(nary,nb_nodes-ioff)
1095 IF (imax.GE.2) THEN
1096 DO i=2,imax
1097 IF (pos_list(ioff+i).NE.curpos+node_rank) THEN
1098 DO l=0,rank_list(ioff+i)-1
1099 acc_lrb%Q(1:m,curpos+node_rank+l) =
1100 & acc_lrb%Q(1:m,pos_list(ioff+i)+l)
1101 acc_lrb%R(curpos+node_rank+l,1:n) =
1102 & acc_lrb%R(pos_list(ioff+i)+l,1:n)
1103 ENDDO
1104 pos_list(ioff+i) = curpos+node_rank
1105 ENDIF
1106 node_rank = node_rank+rank_list(ioff+i)
1107 ENDDO
1108 CALL init_lrb(lrb,node_rank,m,n,.true.)
1109 IF (.NOT.resort.OR.level.EQ.0) THEN
1110 lrb%Q => acc_lrb%Q(1:m,curpos:curpos+node_rank)
1111 lrb%R => acc_lrb%R(curpos:curpos+node_rank,1:n)
1112 ELSE
1113 lrb%Q => acc_tmp%Q(1:m,curpos:curpos+node_rank)
1114 lrb%R => acc_tmp%R(curpos:curpos+node_rank,1:n)
1115 ENDIF
1116 new_acc_rank = node_rank-rank_list(ioff+1)
1117 IF (new_acc_rank.GT.0) THEN
1118 CALL smumps_recompress_acc(lrb,
1119 & maxi_cluster, maxi_rank, a, la, poseltt,
1120 & nfront, niv, midblk_compress, toleps,
1121 & tol_opt, kpercent_rmb, kpercent_lua, new_acc_rank)
1122 ENDIF
1123 rank_list_new(j) = lrb%K
1124 pos_list_new(j) = curpos
1125 ELSE
1126 rank_list_new(j) = node_rank
1127 pos_list_new(j) = curpos
1128 ENDIF
1129 ioff = ioff+imax
1130 ENDDO
1131 IF (nb_nodes_new.GT.1) THEN
1132 IF (resort) THEN
1133 ktot = sum(rank_list_new)
1134 CALL init_lrb(acc_new,ktot,m,n,.true.)
1135 ALLOCATE(acc_new%Q(maxi_cluster,maxi_rank),
1136 & acc_new%R(maxi_rank,maxi_cluster), stat=allocok)
1137 IF (allocok > 0) THEN
1138 write(*,*) 'Allocation error of ACC_NEW%Q/ACC_NEW%R ',
1139 & 'in SMUMPS_RECOMPRESS_ACC_NARYTREE'
1140 call mumps_abort()
1141 ENDIF
1142 CALL mumps_sort_int(nb_nodes_new, rank_list_new,
1143 & pos_list_new)
1144 curpos = 1
1145 IF (level.EQ.0) THEN
1146 lrb_ptr => acc_lrb
1147 ELSE
1148 lrb_ptr => acc_tmp
1149 ENDIF
1150 DO j=1,nb_nodes_new
1151 DO l=0,rank_list_new(j)-1
1152 acc_new%Q(1:m,curpos+l) =
1153 & lrb_ptr%Q(1:m,pos_list_new(j)+l)
1154 acc_new%R(curpos+l,1:n) =
1155 & lrb_ptr%R(pos_list_new(j)+l,1:n)
1156 ENDDO
1157 pos_list_new(j) = curpos
1158 curpos = curpos + rank_list_new(j)
1159 ENDDO
1160 IF (level.GT.0) THEN
1161 CALL dealloc_lrb(acc_tmp, keep8, 4)
1162 ENDIF
1163 CALL smumps_recompress_acc_narytree(acc_lrb,
1164 & maxi_cluster, maxi_rank, a, la, poseltt, keep8,
1165 & nfront, niv, midblk_compress, toleps, tol_opt,
1166 & kpercent_rmb, kpercent_lua, k478,
1167 & rank_list_new, pos_list_new, nb_nodes_new,
1168 & level+1, acc_new)
1169 ELSE
1170 CALL smumps_recompress_acc_narytree(acc_lrb,
1171 & maxi_cluster, maxi_rank, a, la, poseltt, keep8,
1172 & nfront, niv, midblk_compress, toleps, tol_opt,
1173 & kpercent_rmb, kpercent_lua, k478,
1174 & rank_list_new, pos_list_new, nb_nodes_new, level+1)
1175 ENDIF
1176 ELSE
1177 IF (pos_list_new(1).NE.1) THEN
1178 write(*,*) 'Internal error in ',
1179 & 'SMUMPS_RECOMPRESS_ACC_NARYTREE', pos_list_new(1)
1180 ENDIF
1181 acc_lrb%K = rank_list_new(1)
1182 IF (resort.AND.level.GT.0) THEN
1183 DO l=1,acc_lrb%K
1184 DO i=1,m
1185 acc_lrb%Q(i,l) = acc_tmp%Q(i,l)
1186 ENDDO
1187 DO i=1,n
1188 acc_lrb%R(l,i) = acc_tmp%R(l,i)
1189 ENDDO
1190 ENDDO
1191 CALL dealloc_lrb(acc_tmp, keep8, 4)
1192 ENDIF
1193 ENDIF
1194 DEALLOCATE(rank_list_new, pos_list_new)
1195 END SUBROUTINE smumps_recompress_acc_narytree
1196 SUBROUTINE smumps_recompress_acc_v2(ACC_LRB, MAXI_CLUSTER,
1197 & MAXI_RANK, A, LA, POSELTT, NFRONT, NIV,
1198 & MIDBLK_COMPRESS, TOLEPS, TOL_OPT, KPERCENT_RMB,
1199 & KPERCENT_LUA, NEW_ACC_RANK)
1200 TYPE(lrb_type),INTENT(INOUT) :: ACC_LRB
1201 INTEGER(8), intent(in) :: LA
1202 REAL, intent(inout) :: A(LA)
1203 INTEGER,INTENT(IN) :: NFRONT, NIV, TOL_OPT
1204 INTEGER,INTENT(IN) :: MAXI_CLUSTER, MAXI_RANK, KPERCENT_LUA
1205 INTEGER,INTENT(INOUT) :: NEW_ACC_RANK
1206 INTEGER(8), INTENT(IN) :: POSELTT
1207 INTEGER,intent(in) :: MIDBLK_COMPRESS, KPERCENT_RMB
1208 REAL, intent(in) :: TOLEPS
1209 REAL, ALLOCATABLE :: RWORK_RRQR(:)
1210 REAL, ALLOCATABLE :: WORK_RRQR(:), TAU_RRQR(:)
1211 REAL, ALLOCATABLE, DIMENSION(:,:), TARGET ::
1212 & q1, r1, q2, proj
1213 INTEGER, ALLOCATABLE :: JPVT_RRQR(:)
1214 INTEGER :: INFO, RANK1, MAXRANK, LWORK
1215 LOGICAL :: BUILDQ1
1216 INTEGER :: I, J, M, N, K, K1
1217 INTEGER :: allocok, MREQ
1218 REAL :: ONE, MONE, ZERO
1219 parameter(one = 1.0e0, mone=-1.0e0)
1220 parameter(zero=0.0e0)
1221 m = acc_lrb%M
1222 n = acc_lrb%N
1223 k = new_acc_rank
1224 k1 = acc_lrb%K - k
1225 maxrank = k-1
1226 maxrank = max(1, int((maxrank*kpercent_lua/100)))
1227 lwork = k*(k+1)
1228 allocate(q1(m,k), proj(k1, k),
1229 & work_rrqr(lwork), rwork_rrqr(2*k),
1230 & tau_rrqr(k),
1231 & jpvt_rrqr(k), stat=allocok)
1232 IF (allocok > 0) THEN
1233 mreq = m * k + k1 * k + lwork + 4 * k
1234 GOTO 100
1235 ENDIF
1236 DO j=1,k
1237 DO i=1,m
1238 q1(i,j) = acc_lrb%Q(i,j+k1)
1239 ENDDO
1240 ENDDO
1241 CALL sgemm('t', 'n', K1, K, M, ONE, ACC_LRB%Q(1,1),
1242 & MAXI_CLUSTER, Q1(1,1), M, ZERO, PROJ(1,1), K1)
1243 CALL sgemm('n', 'n', M, K, K1, MONE, ACC_LRB%Q(1,1),
1244 & MAXI_CLUSTER, PROJ(1,1), K1, ONE, Q1(1,1), M)
1245 JPVT_RRQR = 0
1246 CALL SMUMPS_TRUNCATED_RRQR(M, K, Q1(1,1),
1247 & M, JPVT_RRQR(1), TAU_RRQR(1), WORK_RRQR(1),
1248 & K, RWORK_RRQR(1), TOLEPS, TOL_OPT,
1249 & RANK1, MAXRANK, INFO,
1250 & BUILDQ1)
1251 IF (BUILDQ1) THEN
1252 allocate(Q2(N,K), stat=allocok)
1253 IF (allocok > 0) THEN
1254 MREQ = N*K
1255 GOTO 100
1256 ENDIF
1257 DO J=1,K
1258 DO I=1,N
1259 Q2(I,J) = ACC_LRB%R(J+K1,I)
1260 ENDDO
1261 ENDDO
1262 CALL sgemm('n', 't', K1, N, K, ONE, PROJ(1,1), K1,
1263 & Q2(1,1), N, ONE, ACC_LRB%R(1,1), MAXI_RANK)
1264.GT. IF (RANK10) THEN
1265 allocate(R1(RANK1,K), stat=allocok)
1266 IF (allocok > 0) THEN
1267 MREQ = RANK1*K
1268 GOTO 100
1269 ENDIF
1270 DO J=1, K
1271 R1(1:MIN(RANK1,J),JPVT_RRQR(J)) =
1272 & Q1(1:MIN(RANK1,J),J)
1273.LT. IF(JRANK1) R1(MIN(RANK1,J)+1:
1274 & RANK1,JPVT_RRQR(J))= ZERO
1275 END DO
1276 CALL sorgqr
1277 & (M, RANK1, RANK1, Q1(1,1),
1278 & M, TAU_RRQR(1),
1279 & WORK_RRQR(1), LWORK, INFO )
1280 DO J=1,K
1281 DO I=1,M
1282 ACC_LRB%Q(I,J+K1) = Q1(I,J)
1283 ENDDO
1284 ENDDO
1285 CALL sgemm('n', 't', RANK1, N, K, ONE, R1(1,1), RANK1,
1286 & Q2(1,1), N, ZERO, ACC_LRB%R(K1+1,1), MAXI_RANK)
1287 deallocate(R1)
1288 ENDIF
1289 deallocate(Q2)
1290 ACC_LRB%K = K1 + RANK1
1291 ENDIF
1292 deallocate(PROJ)
1293 deallocate(Q1, JPVT_RRQR, TAU_RRQR, WORK_RRQR, RWORK_RRQR)
1294 RETURN
1295 100 CONTINUE
1296C Alloc NOT ok!!
1297 write(*,*) 'allocation problem in blr routine
1298 & smumps_recompress_acc_v2: ',
1299 & 'not enough memory? memory requested = ' , MREQ
1300 CALL MUMPS_ABORT()
1301 RETURN
1302 END SUBROUTINE SMUMPS_RECOMPRESS_ACC_V2
1303 SUBROUTINE MAX_CLUSTER(CUT,CUT_SIZE,MAXI_CLUSTER)
1304 INTEGER, intent(in) :: CUT_SIZE
1305 INTEGER, intent(out) :: MAXI_CLUSTER
1306 INTEGER, POINTER, DIMENSION(:) :: CUT
1307 INTEGER :: I
1308 MAXI_CLUSTER = 0
1309 DO I = 1, CUT_SIZE
1310.GE. IF (CUT(I+1) - CUT(I) MAXI_CLUSTER) THEN
1311 MAXI_CLUSTER = CUT(I+1) - CUT(I)
1312 END IF
1313 END DO
1314 END SUBROUTINE MAX_CLUSTER
1315 SUBROUTINE SMUMPS_GET_LUA_ORDER(NB_BLOCKS, ORDER, RANK, IWHANDLER,
1316 & SYM, FS_OR_CB, I, J, FRFR_UPDATES,
1317 & LBANDSLAVE_IN, K474, BLR_U_COL)
1318C -----------
1319C Parameters
1320C -----------
1321 INTEGER, INTENT(IN) :: NB_BLOCKS, IWHANDLER, SYM, FS_OR_CB, I, J
1322 INTEGER, INTENT(OUT) :: ORDER(NB_BLOCKS), RANK(NB_BLOCKS),
1323 & FRFR_UPDATES
1324 LOGICAL, OPTIONAL, INTENT(IN) :: LBANDSLAVE_IN
1325 INTEGER, OPTIONAL, INTENT(IN) :: K474
1326 TYPE(LRB_TYPE), POINTER, OPTIONAL :: BLR_U_COL(:)
1327C -----------
1328C Local variables
1329C -----------
1330 INTEGER :: K, IND_L, IND_U
1331 LOGICAL :: LBANDSLAVE
1332 TYPE(LRB_TYPE), POINTER :: BLR_L(:), BLR_U(:)
1333 IF (PRESENT(LBANDSLAVE_IN)) THEN
1334 LBANDSLAVE = LBANDSLAVE_IN
1335 ELSE
1336 LBANDSLAVE = .FALSE.
1337 ENDIF
1338.NE..AND..EQ..AND..NE. IF ((SYM0)(FS_OR_CB0)(J0)) THEN
1339 write(6,*) 'internal error in smumps_get_lua_order',
1340 & 'sym, fs_or_cb, j = ',SYM,FS_OR_CB,J
1341 CALL MUMPS_ABORT()
1342 ENDIF
1343 FRFR_UPDATES = 0
1344 DO K = 1, NB_BLOCKS
1345 ORDER(K) = K
1346.EQ. IF (FS_OR_CB0) THEN ! FS
1347.EQ. IF (J0) THEN ! L panel
1348 IND_L = NB_BLOCKS+I-K
1349 IND_U = NB_BLOCKS+1-K
1350 ELSE ! U panel
1351 IND_L = NB_BLOCKS+1-K
1352 IND_U = NB_BLOCKS+I-K
1353 ENDIF
1354 ELSE ! CB
1355 IND_L = I-K
1356 IND_U = J-K
1357 ENDIF
1358 IF (LBANDSLAVE) THEN
1359 IND_L = I
1360.GE. IF (K4742) THEN
1361 IND_U = K
1362 ENDIF
1363 ENDIF
1364 CALL SMUMPS_BLR_RETRIEVE_PANEL_LORU(
1365 & IWHANDLER,
1366 & 0, ! L Panel
1367 & K, BLR_L)
1368.EQ. IF (SYM0) THEN
1369.AND..GE. IF (LBANDSLAVEK4742) THEN
1370 BLR_U => BLR_U_COL
1371 ELSE
1372 CALL SMUMPS_BLR_RETRIEVE_PANEL_LORU(
1373 & IWHANDLER,
1374 & 1, ! L Panel
1375 & K, BLR_U)
1376 ENDIF
1377 ELSE
1378 BLR_U => BLR_L
1379 ENDIF
1380 IF (BLR_L(IND_L)%ISLR) THEN
1381 IF (BLR_U(IND_U)%ISLR) THEN
1382 RANK(K) = min(BLR_L(IND_L)%K, BLR_U(IND_U)%K)
1383 ELSE
1384 RANK(K) = BLR_L(IND_L)%K
1385 ENDIF
1386 ELSE
1387 IF (BLR_U(IND_U)%ISLR) THEN
1388 RANK(K) = BLR_U(IND_U)%K
1389 ELSE
1390 RANK(K) = -1
1391 FRFR_UPDATES = FRFR_UPDATES + 1
1392 ENDIF
1393 ENDIF
1394 ENDDO
1395 CALL MUMPS_SORT_INT(NB_BLOCKS, RANK, ORDER)
1396 END SUBROUTINE SMUMPS_GET_LUA_ORDER
1397 SUBROUTINE SMUMPS_BLR_ASM_NIV1 (A, LA, POSEL1, NFRONT, NASS1,
1398 & IWHANDLER, SON_IW, LIW, LSTK, NELIM, K1, K2, SYM,
1399 & KEEP, KEEP8, OPASSW)
1400C
1401C Purpose
1402C =======
1403C
1404C Called by a level 1 master assembling the contribution
1405C block of a level 1 son that has been BLR-compressed
1406C
1407C
1408C Parameters
1409C ==========
1410C
1411 INTEGER(8) :: LA, POSEL1
1412 INTEGER :: LIW, NFRONT, NASS1, LSTK, NELIM, K1, K2, IWHANDLER
1413 REAL :: A(LA)
1414C INTEGER :: SON_IW(LIW)
1415 INTEGER :: SON_IW(:) ! contiguity information lost but no copy
1416 INTEGER :: KEEP(500)
1417 INTEGER(8) :: KEEP8(150)
1418 INTEGER :: SYM
1419 DOUBLE PRECISION, INTENT(INOUT) :: OPASSW
1420C
1421C Local variables
1422C ===============
1423C
1424 REAL, ALLOCATABLE :: SON_A(:)
1425 INTEGER(8) :: APOS, SON_APOS, IACHK, JJ2, NFRONT8
1426 INTEGER :: KK, KK1, allocok, SON_LA
1427 TYPE(LRB_TYPE), POINTER :: CB_LRB(:,:), LRB
1428 INTEGER, POINTER, DIMENSION(:) :: BEGS_BLR_DYNAMIC
1429 INTEGER :: NB_INCB, NB_INASM, NB_BLR, I, J, M, N, II, NPIV,
1430 & IBIS, IBIS_END, FIRST_ROW, LAST_ROW, FIRST_COL, LAST_COL,
1431 & SON_LDA
1432 DOUBLE PRECISION :: PROMOTE_COST
1433 REAL :: ONE, ZERO
1434 PARAMETER (ONE = 1.0E0)
1435 PARAMETER (ZERO = 0.0D0)
1436 CALL SMUMPS_BLR_RETRIEVE_BEGSBLR_DYN(IWHANDLER,
1437 & BEGS_BLR_DYNAMIC)
1438 CALL SMUMPS_BLR_RETRIEVE_CB_LRB(IWHANDLER, CB_LRB)
1439 NB_BLR = size(BEGS_BLR_DYNAMIC)-1
1440 NB_INCB = size(CB_LRB,1)
1441 NB_INASM = NB_BLR - NB_INCB
1442 NPIV = BEGS_BLR_DYNAMIC(NB_INASM+1)-1
1443 NFRONT8 = int(NFRONT,8)
1444.EQ. IF (SYM0) THEN
1445 IBIS_END = NB_INCB*NB_INCB
1446 ELSE
1447 IBIS_END = NB_INCB*(NB_INCB+1)/2
1448 ENDIF
1449#if defined(BLR_MT)
1450!$OMP PARALLEL
1451!$OMP DO PRIVATE(IBIS, I, J, M, N, SON_LA, SON_LDA, FIRST_ROW,
1452!$OMP& LAST_ROW, FIRST_COL, LAST_COL, LRB, SON_A, II, KK,
1453!$OMP& APOS, IACHK, KK1, JJ2, PROMOTE_COST, allocok, SON_APOS)
1454#endif
1455 DO IBIS = 1,IBIS_END
1456C Determining I,J from IBIS
1457.EQ. IF (SYM0) THEN
1458 I = (IBIS-1)/NB_INCB+1
1459 J = IBIS - (I-1)*NB_INCB
1460 ELSE
1461 I = ceiling((1.0D0+sqrt(1.0D0+8.0D0*dble(IBIS)))/2.0D0)-1
1462 J = IBIS - I*(I-1)/2
1463 ENDIF
1464 I = I+NB_INASM
1465 J = J+NB_INASM
1466.EQ. IF (INB_INASM+1) THEN
1467C first CB block, add NELIM because FIRST_ROW starts at NELIM+1
1468 FIRST_ROW = BEGS_BLR_DYNAMIC(I)-NPIV+NELIM
1469 ELSE
1470 FIRST_ROW = BEGS_BLR_DYNAMIC(I)-NPIV
1471 ENDIF
1472 LAST_ROW = BEGS_BLR_DYNAMIC(I+1)-1-NPIV
1473 M=LAST_ROW-FIRST_ROW+1
1474 FIRST_COL = BEGS_BLR_DYNAMIC(J)-NPIV
1475 LAST_COL = BEGS_BLR_DYNAMIC(J+1)-1-NPIV
1476 N = BEGS_BLR_DYNAMIC(J+1)-BEGS_BLR_DYNAMIC(J)
1477 SON_APOS = 1_8
1478 SON_LA = M*N
1479 SON_LDA = N
1480 LRB => CB_LRB(I-NB_INASM,J-NB_INASM)
1481.AND..EQ. IF (LRB%ISLRLRB%K0) THEN
1482C No need to perform extend-add
1483 CALL DEALLOC_LRB(LRB, KEEP8, KEEP(34))
1484 NULLIFY(LRB)
1485 CYCLE
1486 ENDIF
1487 allocate(SON_A(SON_LA),stat=allocok)
1488.GT. IF (allocok0) THEN
1489 write(*,*) 'not enough memory in smumps_blr_asm_niv1',
1490 & ", Memory requested = ", SON_LA
1491 CALL MUMPS_ABORT()
1492 ENDIF
1493C decompress block
1494 IF (LRB%ISLR) THEN
1495 CALL sgemm('t', 't', N, M, LRB%K, ONE, LRB%R(1,1), LRB%K,
1496 & LRB%Q(1,1), M, ZERO, SON_A(SON_APOS), SON_LDA)
1497 PROMOTE_COST = 2.0D0*M*N*LRB%K
1498 CALL UPD_FLOP_DECOMPRESS(PROMOTE_COST, .TRUE.)
1499 ELSE
1500.EQ..AND..NE. IF (IJSYM0) THEN
1501C Diag block and LDLT, copy only lower half
1502.EQ..AND..GT. IF (J-NB_INASM1NELIM0) THEN
1503C The first diagonal block is rectangular !!
1504C with NELIM more cols than rows
1505 DO II=1,M
1506 DO KK=1,II+NELIM
1507 SON_A(SON_APOS+int(II-1,8)*int(SON_LDA,8) +
1508 & int(KK-1,8))
1509 & = LRB%Q(II,KK)
1510 ENDDO
1511 ENDDO
1512 ELSE
1513 DO II=1,M
1514 DO KK=1,II
1515 SON_A(SON_APOS+int(II-1,8)*int(SON_LDA,8) +
1516 & int(KK-1,8))
1517 & = LRB%Q(II,KK)
1518 ENDDO
1519 ENDDO
1520 ENDIF
1521 ELSE
1522 DO II=1,M
1523 DO KK=1,N
1524 SON_A(SON_APOS+int(II-1,8)*int(SON_LDA,8) +
1525 & int(KK-1,8))
1526 & = LRB%Q(II,KK)
1527 ENDDO
1528 ENDDO
1529 ENDIF
1530 ENDIF
1531C Deallocate block
1532 CALL DEALLOC_LRB(LRB, KEEP8, KEEP(34))
1533 NULLIFY(LRB)
1534C extend add in father
1535.NE..AND..EQ..AND..GT. IF (SYM0J-NB_INASM1NELIM0) THEN
1536C Case of LDLT with NELIM: first-block column is treated
1537C differently as the NELIM are assembled at the end of the
1538C father
1539 DO KK = FIRST_ROW, LAST_ROW
1540 IACHK = 1_8 + int(KK-FIRST_ROW,8)*int(SON_LDA,8)
1541.LE. IF (SON_IW(KK+K1-1)NASS1) THEN
1542C Fully summed row of the father => permute destination in
1543C father, symmetric swap to be done
1544C First NELIM columns
1545 APOS = POSEL1 + int(SON_IW(KK+K1-1),8) - 1_8
1546 DO KK1 = FIRST_COL, FIRST_COL+NELIM-1
1547 JJ2 = APOS + int(SON_IW(K1+KK1-1)-1,8)*NFRONT8
1548 A(JJ2) = A(JJ2) + SON_A(IACHK + int(KK1-FIRST_COL,8))
1549 ENDDO
1550C Remaining columns
1551 APOS = POSEL1 + int(SON_IW(KK+K1-1)-1,8)*NFRONT8
1552C DO KK1 = FIRST_COL+NELIM, LAST_COL
1553C In case I=J and first block, one may have
1554C LAST_COL > KK, but only lower triangular part
1555C should be assembled. We use min(LAST_COL,KK)
1556C below index to cover this case.
1557 DO KK1 = FIRST_COL+NELIM, min(LAST_COL,KK)
1558 JJ2 = APOS + int(SON_IW(K1+KK1-1),8) - 1_8
1559 A(JJ2) = A(JJ2) + SON_A(IACHK + int(KK1-FIRST_COL,8))
1560 ENDDO
1561 ELSE
1562 APOS = POSEL1 + int(SON_IW(KK+K1-1)-1,8)*NFRONT8
1563 DO KK1 = FIRST_COL, min(LAST_COL,KK)
1564 JJ2 = APOS + int(SON_IW(K1+KK1-1),8) - 1_8
1565 A(JJ2) = A(JJ2) + SON_A(IACHK + int(KK1-FIRST_COL,8))
1566 ENDDO
1567 ENDIF
1568 ENDDO
1569 ELSE
1570C Case of LDLT without NELIM or LU: everything is simpler
1571 DO KK = FIRST_ROW, LAST_ROW
1572 APOS = POSEL1 + int(SON_IW(KK+K1-1)-1,8)*NFRONT8
1573 IACHK = 1_8 + int(KK-FIRST_ROW,8)*int(SON_LDA,8)
1574.EQ..AND..NE. IF (IJSYM0) THEN
1575C LDLT diag block: assemble only lower half
1576 DO KK1 = FIRST_COL, KK
1577 JJ2 = APOS + int(SON_IW(K1+KK1-1),8) - 1_8
1578 A(JJ2) = A(JJ2) + SON_A(IACHK + int(KK1-FIRST_COL,8))
1579 ENDDO
1580 ELSE
1581 DO KK1 = FIRST_COL, LAST_COL
1582 JJ2 = APOS + int(SON_IW(K1+KK1-1),8) - 1_8
1583 A(JJ2) = A(JJ2) + SON_A(IACHK + int(KK1-FIRST_COL,8))
1584 ENDDO
1585 ENDIF
1586 ENDDO
1587 ENDIF
1588C Deallocate SON_A
1589 DEALLOCATE(SON_A)
1590 ENDDO
1591#if defined(BLR_MT)
1592!$OMP END DO
1593!$OMP END PARALLEL
1594#endif
1595 CALL SMUMPS_BLR_FREE_CB_LRB(IWHANDLER,
1596C Only CB_LRB structure is left to deallocate
1597 & .TRUE., KEEP8, KEEP(34))
1598.EQ..OR..EQ. IF ((KEEP(486)3)KEEP(486)0) THEN
1599C Case of FR solve: the BLR structure could not be freed
1600C in SMUMPS_END_FACTO_SLAVE and should be freed here
1601C Not reachable in case of error: set INFO1 to 0
1602 CALL SMUMPS_BLR_END_FRONT(IWHANDLER, 0, KEEP8, KEEP(34),
1603 & MTK405=KEEP(405))
1604 ENDIF
1605 END SUBROUTINE SMUMPS_BLR_ASM_NIV1
1606 END MODULE SMUMPS_LR_CORE
1607C --------------------------------------------------------------------
1608 SUBROUTINE SMUMPS_TRUNCATED_RRQR( M, N, A, LDA, JPVT, TAU, WORK,
1609 & LDW, RWORK, TOLEPS, TOL_OPT, RANK, MAXRANK, INFO,
1610 & ISLR)
1611C This routine computes a Rank-Revealing QR factorization of a dense
1612C matrix A. The factorization is truncated when the absolute value of
1613C a diagonal coefficient of the R factor becomes smaller than a
1614C prescribed threshold TOLEPS. The resulting partial Q and R factors
1615C provide a rank-k approximation of the input matrix A with accuracy
1616C TOLEPS.
1617C
1618C This routine is obtained by merging the LAPACK
1619C (http://www.netlib.org/lapack/) CGEQP3 and CLAQPS routines and by
1620C applying a minor modification to the outer factorization loop in
1621C order to stop computations as soon as possible when the required
1622C accuracy is reached.
1623C
1624C Copyright (c) 1992-2017 The University of Tennessee and The
1625C University of Tennessee Research Foundation. All rights reserved.
1626C Copyright (c) 2000-2017 The University of California Berkeley.
1627C All rights reserved.
1628C Copyright (c) 2006-2017 The University of Colorado Denver.
1629C All rights reserved.
1630C
1631C Redistribution and use in source and binary forms, with or without
1632C modification, are permitted provided that the following conditions
1633C are met:
1634C
1635C - Redistributions of source code must retain the above copyright
1636C notice, this list of conditions and the following disclaimer.
1637C
1638C - Redistributions in binary form must reproduce the above
1639C copyright notice, this list of conditions and the following
1640C disclaimer listed in this license in the documentation and/or
1641C other materials provided with the distribution.
1642C
1643C - Neither the name of the copyright holders nor the names of its
1644C contributors may be used to endorse or promote products derived from
1645C this software without specific prior written permission.
1646C
1647C The copyright holders provide no reassurances that the source code
1648C provided does not infringe any patent, copyright, or any other
1649C intellectual property rights of third parties. The copyright holders
1650C disclaim any liability to any recipient for claims brought against
1651C recipient by any third party for infringement of that parties
1652C intellectual property rights.
1653C
1654C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
1655C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
1656C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
1657C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
1658C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
1659C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
1660C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
1661C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
1662C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
1663C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
1664C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
1665C
1666 IMPLICIT NONE
1667C
1668 INTEGER :: INFO, LDA, LDW, M, N, RANK, MAXRANK
1669C TOL_OPT controls the tolerance option used
1670C >0 => use 2-norm (||.||_X = ||.||_2)
1671C <0 => use Frobenius-norm (||.||_X = ||.||_F)
1672C Furthermore, depending on abs(TOL_OPT):
1673C 1 => absolute: ||B_{I(k+1:end),J(k+1:end)}||_X <= TOLEPS
1674C 2 => relative to 2-norm of the compressed block:
1675C ||B_{I(k+1:end),J(k+1:end)}||_X <= TOLEPS*||B_{I,J}||_2
1676C 3 => relative to the max of the 2-norms of the row and column diagonal blocks
1677C ||B_{I(k+1:end),J{k+1:end}}||_X <= TOLEPS*max(||B_{I,I}||_2,||B_{J,J}||_2)
1678C 4 => relative to the sqrt of product of the 2-norms of the row and column diagonal blocks
1679C ||B_{I(k+1:end),J{k+1:end}}||_X <= TOLEPS*sqrt(||B_{I,I}||_2*||B_{J,J}||_2)
1680 INTEGER :: TOL_OPT
1681 REAL :: TOLEPS
1682 INTEGER :: JPVT(*)
1683 REAL :: RWORK(*)
1684 REAL :: A(LDA,*), TAU(*)
1685 REAL :: WORK(LDW,*)
1686 LOGICAL :: ISLR
1687 REAL :: TOLEPS_EFF, TRUNC_ERR
1688 INTEGER, PARAMETER :: INB=1, INBMIN=2
1689 INTEGER :: J, JB, MINMN, NB
1690 INTEGER :: OFFSET, ITEMP
1691 INTEGER :: LSTICC, PVT, K, RK
1692 REAL :: TEMP, TEMP2, TOL3Z
1693 REAL :: AKK
1694 LOGICAL INADMISSIBLE
1695 REAL, PARAMETER :: RZERO=0.0E+0, RONE=1.0E+0
1696 REAL :: ZERO
1697 REAL :: ONE
1698 PARAMETER ( ONE = 1.0E+0 )
1699 PARAMETER ( ZERO = 0.0E+0 )
1700 REAL :: slamch
1701 INTEGER :: ilaenv, isamax
1702 EXTERNAL :: isamax, slamch
1703 EXTERNAL sgeqrf, sormqr, xerbla
1704 EXTERNAL ilaenv
1705 EXTERNAL sgemm, sgemv, slarfg, sswap
1706 REAL, EXTERNAL :: snrm2
1707 INFO = 0
1708 ISLR = .FALSE.
1709.LT. IF( M0 ) THEN
1710 INFO = -1
1711.LT. ELSE IF( N0 ) THEN
1712 INFO = -2
1713.LT. ELSE IF( LDAMAX( 1, M ) ) THEN
1714 INFO = -4
1715 END IF
1716.EQ. IF( INFO0 ) THEN
1717.LT. IF( LDWN ) THEN
1718 INFO = -8
1719 END IF
1720 END IF
1721.NE. IF( INFO0 ) THEN
1722 WRITE(*,999) -INFO
1723 RETURN
1724 END IF
1725 MINMN = MIN(M,N)
1726.EQ. IF( MINMN0 ) THEN
1727 RANK = 0
1728 RETURN
1729 END IF
1730 NB = ilaenv( INB, 'cgeqrf', ' ', M, N, -1, -1 )
1731 SELECT CASE(abs(TOL_OPT))
1732 CASE(1)
1733 TOLEPS_EFF = TOLEPS
1734 CASE(2)
1735C TOLEPS_EFF will be computed at step K=1 below
1736 CASE DEFAULT
1737 write(*,*) 'internal error in smumps_truncated_rrqr: tol_opt =',
1738 & TOL_OPT
1739 CALL MUMPS_ABORT()
1740 END SELECT
1741 TOLEPS_EFF = TOLEPS
1742C
1743C Avoid pointers (and TARGET attribute on RWORK/WORK)
1744C because of implicit interface. An implicit interface
1745C is needed to avoid intermediate array copies
1746C VN1 => RWORK(1:N)
1747C VN2 => RWORK(N+1:2*N)
1748C AUXV => WORK(1:LDW,1:1)
1749C F => WORK(1:LDW,2:NB+1)
1750C LDF = LDW
1751* Initialize partial column norms. The first N elements of work
1752* store the exact column norms.
1753 DO J = 1, N
1754C VN1( J ) = snrm2( M, A( 1, J ), 1 )
1755 RWORK( J ) = snrm2( M, A( 1, J ), 1 )
1756C VN2( J ) = VN1( J )
1757 RWORK( N + J ) = RWORK( J )
1758 JPVT(J) = J
1759 END DO
1760.LT. IF (TOL_OPT0) THEN
1761C Compute TRUNC_ERR for first step
1762C TRUNC_ERR = snrm2( N, VN1( 1 ), 1 )
1763 TRUNC_ERR = snrm2( N, RWORK( 1 ), 1 )
1764 ENDIF
1765 OFFSET = 0
1766 TOL3Z = SQRT(slamch('epsilon'))
1767 DO
1768 JB = MIN(NB,MINMN-OFFSET)
1769 LSTICC = 0
1770 K = 0
1771 DO
1772.EQ. IF(KJB) EXIT
1773 K = K+1
1774 RK = OFFSET+K
1775C PVT = ( RK-1 ) + ISAMAX( N-RK+1, VN1( RK ), 1 )
1776 PVT = ( RK-1 ) + isamax( N-RK+1, RWORK( RK ), 1 )
1777.EQ. IF (RK1) THEN
1778C IF (abs(TOL_OPT).EQ.2) TOLEPS_EFF = VN1(PVT)*TOLEPS
1779.EQ. IF (abs(TOL_OPT)2) TOLEPS_EFF = RWORK(PVT)*TOLEPS
1780 ENDIF
1781.GT. IF (TOL_OPT0) THEN
1782C TRUNC_ERR = VN1(PVT)
1783 TRUNC_ERR = RWORK(PVT)
1784C ELSE
1785C TRUNC_ERR has been already computed at previous step
1786 ENDIF
1787.LT. IF(TRUNC_ERRTOLEPS_EFF) THEN
1788 RANK = RK-1
1789 ISLR = .TRUE.
1790 RETURN
1791 ENDIF
1792.GT. INADMISSIBLE = (RKMAXRANK)
1793 IF (INADMISSIBLE) THEN
1794 RANK = RK
1795 INFO = RK
1796 ISLR = .FALSE.
1797 RETURN
1798 END IF
1799.NE. IF( PVTRK ) THEN
1800 CALL sswap( M, A( 1, PVT ), 1, A( 1, RK ), 1 )
1801c CALL sswap( K-1, F( PVT-OFFSET, 1 ), LDF,
1802c & F( K, 1 ), LDF )
1803 CALL sswap( K-1, WORK( PVT-OFFSET, 2 ), LDW,
1804 & WORK( K, 2 ), LDW )
1805 ITEMP = JPVT(PVT)
1806 JPVT(PVT) = JPVT(RK)
1807 JPVT(RK) = ITEMP
1808C VN1(PVT) = VN1(RK)
1809C VN2(PVT) = VN2(RK)
1810 RWORK(PVT) = RWORK(RK)
1811 RWORK(N+PVT) = RWORK(N+RK)
1812 END IF
1813* Apply previous Householder reflectors to column K:
1814* A(RK:M,RK) := A(RK:M,RK) - A(RK:M,OFFSET+1:RK-1)*F(K,1:K-1)**H.
1815.GT. IF( K1 ) THEN
1816 CALL sgemv( 'no transpose', M-RK+1, K-1, -ONE,
1817C & A(RK,OFFSET+1), LDA, F(K,1), LDF,
1818 & A(RK,OFFSET+1), LDA, WORK(K,2), LDW,
1819 & ONE, A(RK,RK), 1 )
1820 END IF
1821* Generate elementary reflector H(k).
1822.LT. IF( RKM ) THEN
1823 CALL slarfg( M-RK+1, A(RK,RK), A(RK+1,RK), 1, TAU(RK) )
1824 ELSE
1825 CALL slarfg( 1, A(RK,RK), A(RK,RK), 1, TAU(RK) )
1826 END IF
1827 AKK = A(RK,RK)
1828 A(RK,RK) = ONE
1829* Compute Kth column of F:
1830* F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**H*A(RK:M,K).
1831.LT. IF( RKN ) THEN
1832 CALL sgemv( 'transpose', M-RK+1, N-RK, TAU(RK),
1833 & A(RK,RK+1), LDA, A(RK,RK), 1, ZERO,
1834C & F( K+1, K ), 1 )
1835 & WORK( K+1, K+1 ), 1 )
1836 END IF
1837* Padding F(1:K,K) with zeros.
1838 DO J = 1, K
1839C F( J, K ) = ZERO
1840 WORK( J, K+1 ) = ZERO
1841 END DO
1842* Incremental updating of F:
1843* F(1:N,K) := F(1:N-OFFSET,K) -
1844* tau(RK)*F(1:N,1:K-1)*A(RK:M,OFFSET+1:RK-1)**H*A(RK:M,RK).
1845.GT. IF( K1 ) THEN
1846 CALL sgemv( 'transpose', M-RK+1, K-1, -TAU(RK),
1847 & A(RK,OFFSET+1), LDA, A(RK,RK), 1, ZERO,
1848 & WORK(1,1), 1 )
1849C & AUXV(1,1), 1 )
1850 CALL sgemv( 'no transpose', n-offset, k-1, one,
1851 & work(1,2), ldw, work(1,1), 1, one, work(1,k+1), 1 )
1852C & F(1,1), LDF, AUXV(1,1), 1, ONE, F(1,K), 1 )
1853 END IF
1854* Update the current row of A:
1855* A(RK,RK+1:N) := A(RK,RK+1:N) - A(RK,OFFSET+1:RK)*F(K+1:N,1:K)**H.
1856 IF( rk.LT.n ) THEN
1857C CALL sgemv( 'No Transpose', N-RK, K, -ONE, F( K+1, 1 ),
1858 CALL sgemv( 'No Transpose', n-rk, k, -one, work( k+1,2 ),
1859 & ldw,
1860 & a( rk, offset+1 ), lda, one, a( rk, rk+1 ), lda )
1861 END IF
1862* Update partial column norms.
1863*
1864 IF( rk.LT.minmn ) THEN
1865 DO j = rk + 1, n
1866C IF( VN1( J ).NE.RZERO ) THEN
1867 IF( rwork( j ).NE.rzero ) THEN
1868*
1869* NOTE: The following 4 lines follow from the analysis in
1870* Lapack Working Note 176.
1871*
1872C TEMP = ABS( A( RK, J ) ) / VN1( J )
1873 temp = abs( a( rk, j ) ) / rwork( j )
1874 temp = max( rzero, ( rone+temp )*( rone-temp ) )
1875C TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
1876 temp2 = temp*( rwork( j ) / rwork( n+j ) )**2
1877 IF( temp2 .LE. tol3z ) THEN
1878C VN2( J ) = REAL( LSTICC )
1879 rwork( n+j ) = real( lsticc )
1880 lsticc = j
1881 ELSE
1882C VN1( J ) = VN1( J )*SQRT( TEMP )
1883 rwork( j ) = rwork( j )*sqrt( temp )
1884 END IF
1885 END IF
1886 END DO
1887 END IF
1888 a( rk, rk ) = akk
1889 IF (lsticc.NE.0) EXIT
1890 IF (tol_opt.LT.0) THEN
1891C Compute TRUNC_ERR for next step
1892C TRUNC_ERR = snrm2( N-RK, VN1( RK+1 ), 1 )
1893 trunc_err = snrm2( n-rk, rwork( rk+1 ), 1 )
1894 ENDIF
1895 END DO
1896* Apply the block reflector to the rest of the matrix:
1897* A(RK+1:M,RK+1:N) := A(RK+1:M,RK+1:N) -
1898* A(RK+1:M,OFFSET+1:RK)*F(K+1:N-OFFSET,1:K)**H.
1899 IF( rk.LT.min(n,m) ) THEN
1900 CALL sgemm( 'No transpose', 'Transpose', m-rk,
1901 & n-rk, k, -one, a(rk+1,offset+1), lda,
1902C & F(K+1,1), LDF, ONE, A(RK+1,RK+1), LDA )
1903 & work(k+1,2), ldw, one, a(rk+1,rk+1), lda )
1904 END IF
1905* Recomputation of difficult columns.
1906 DO WHILE( lsticc.GT.0 )
1907C ITEMP = NINT( VN2( LSTICC ) )
1908 itemp = nint( rwork( n + lsticc ) )
1909C VN1( LSTICC ) = snrm2( M-RK, A( RK+1, LSTICC ), 1 )
1910 rwork( lsticc ) = snrm2( m-rk, a( rk+1, lsticc ), 1 )
1911*
1912* NOTE: The computation of RWORK( LSTICC ) relies on the fact that
1913* SNRM2 does not fail on vectors with norm below the value of
1914* SQRT(DLAMCH('S'))
1915*
1916C VN2( LSTICC ) = VN1( LSTICC )
1917 rwork( n + lsticc ) = rwork( lsticc )
1918 lsticc = itemp
1919 END DO
1920 IF(rk.GE.minmn) EXIT
1921 offset = rk
1922 IF (tol_opt.LT.0) THEN
1923C Compute TRUNC_ERR for next step
1924C TRUNC_ERR = snrm2( N-RK, VN1( RK+1 ), 1 )
1925 trunc_err = snrm2( n-rk, rwork( rk+1 ), 1 )
1926 ENDIF
1927 END DO
1928 rank = rk
1929 islr = .NOT.(rk.GT.maxrank)
1930 RETURN
1931 999 FORMAT ('On entry to SMUMPS_TRUNCATED_RRQR, parameter number',
1932 & i2,' had an illegal value')
1933 END SUBROUTINE smumps_truncated_rrqr
mumps_abort
#define mumps_abort
Definition VE_Metis.h:25
cgeqrf
subroutine cgeqrf(m, n, a, lda, tau, work, lwork, info)
CGEQRF
Definition cgeqrf.f:146
sorgqr
subroutine sorgqr(m, n, k, a, lda, tau, work, lwork, info)
SORGQR
Definition sorgqr.f:128
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
snrm2
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:156
strsm
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:187
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mumps_lr_common
Definition lr_common.F:14
mumps_lr_common::compute_blr_vcs
subroutine compute_blr_vcs(k472, ibcksz, maxsize, nass)
Definition lr_common.F:18
smumps_lr_core
Definition slr_core.F:18
smumps_lr_core::smumps_recompress_acc_v2
subroutine smumps_recompress_acc_v2(acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, nfront, niv, midblk_compress, toleps, tol_opt, kpercent_rmb, kpercent_lua, new_acc_rank)
Definition slr_core.F:1200
smumps_lr_core::smumps_get_lua_order
subroutine smumps_get_lua_order(nb_blocks, order, rank, iwhandler, sym, fs_or_cb, i, j, frfr_updates, lbandslave_in, k474, blr_u_col)
Definition slr_core.F:1318
smumps_lr_core::alloc_lrb
subroutine alloc_lrb(lrb_out, k, m, n, islr, iflag, ierror, keep8)
Definition slr_core.F:111
smumps_lr_core::smumps_decompress_acc
subroutine smumps_decompress_acc(acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, nfront, niv, loru, count_flops)
Definition slr_core.F:763
smumps_lr_core::smumps_recompress_acc
subroutine smumps_recompress_acc(acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, nfront, niv, midblk_compress, toleps, tol_opt, kpercent_rmb, kpercent_lua, new_acc_rank)
Definition slr_core.F:868
smumps_lr_core::alloc_lrb_from_acc
subroutine alloc_lrb_from_acc(acc_lrb, lrb_out, k, m, n, loru, iflag, ierror, keep8)
Definition slr_core.F:160
smumps_lr_core::smumps_compress_fr_updates
subroutine smumps_compress_fr_updates(acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, nfront, niv, toleps, tol_opt, kpercent, buildq, loru, cb_compress)
Definition slr_core.F:788
smumps_lr_core::smumps_lrtrsm
subroutine smumps_lrtrsm(a, la, poselt_local, nfront, lda, lrb, niv, sym, loru, iw, offset_iw)
Definition slr_core.F:270
smumps_lr_core::smumps_lrgemm_scaling
subroutine smumps_lrgemm_scaling(lrb, scaled, a, la, diag, ld_diag, iw2, poseltt, nfront, block, maxi_cluster)
Definition slr_core.F:356
smumps_lr_core::smumps_recompress_acc_narytree
recursive subroutine smumps_recompress_acc_narytree(acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, keep8, nfront, niv, midblk_compress, toleps, tol_opt, kpercent_rmb, kpercent_lua, k478, rank_list, pos_list, nb_nodes, level, acc_tmp)
Definition slr_core.F:1053
smumps_lr_core::smumps_lrgemm4
subroutine smumps_lrgemm4(alpha, lrb1, lrb2, beta, a, la, poseltt, nfront, sym, iflag, ierror, midblk_compress, toleps, tol_opt, kpercent, rank, buildq, lua_activated, loru, lrb3, maxi_rank, maxi_cluster, diag, ld_diag, iw2, block)
Definition slr_core.F:406
smumps_lr_core::init_lrb
subroutine init_lrb(lrb_out, k, m, n, islr)
Definition slr_core.F:27
smumps_lr_core::smumps_blr_asm_niv1
subroutine smumps_blr_asm_niv1(a, la, posel1, nfront, nass1, iwhandler, son_iw, liw, lstk, nelim, k1, k2, sym, keep, keep8, opassw)
Definition slr_core.F:1400
smumps_lr_core::regrouping2
subroutine regrouping2(cut, npartsass, nass, npartscb, ncb, ibcksz, onlycb, k472)
Definition slr_core.F:184
smumps_lr_core::is_front_blr_candidate
subroutine is_front_blr_candidate(inode, niv, nfront, nass, blron, k489, k490, k491, k492, k20, k60, idad, k38, lrstatus, n, lrgroups)
Definition slr_core.F:45
smumps_lr_data_m
Definition smumps_lr_data_m.F:14
smumps_lr_stats
Definition slr_stats.F:14
smumps_lr_stats::upd_flop_compress
subroutine upd_flop_compress(lr_b, rec_acc, cb_compress, frswap)
Definition slr_stats.F:154
smumps_lr_stats::upd_flop_update
subroutine upd_flop_update(lrb1, lrb2, midblk_compress, rank_in, buildq, is_symdiag, lua_activated, rec_acc)
Definition slr_stats.F:226
smumps_lr_stats::upd_flop_trsm
subroutine upd_flop_trsm(lrb, loru)
Definition slr_stats.F:199
smumps_lr_type
Definition slr_type.F:14
smumps_lr_type::dealloc_lrb
subroutine dealloc_lrb(lrb_out, keep8, k34)
Definition slr_type.F:25
smumps_truncated_rrqr
subroutine smumps_truncated_rrqr(m, n, a, lda, jpvt, tau, work, ldw, rwork, toleps, tol_opt, rank, maxrank, info, islr)
Definition slr_core.F:1611
smumps_lr_type::lrb_type
Definition slr_type.F:16
mumps_sort_int
subroutine mumps_sort_int(n, val, id)
Definition tools_common.F:709
mumps_dm_fac_upd_dyn_memcnts
subroutine mumps_dm_fac_upd_dyn_memcnts(mem_count_allocated, atomic_updates, keep8, iflag, ierror, k69upd, k71upd)
Definition tools_common.F:1476