OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
clr_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 cmumps_lr_core
19 USE mumps_lr_common
20 USE cmumps_lr_type
21 USE cmumps_lr_stats
22 USE cmumps_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 COMPLEX :: ZERO
118 parameter(zero=(0.0e0,0.0e0))
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 cmumps_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 COMPLEX, 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 COMPLEX, POINTER :: LR_BLOCK_PTR(:,:)
286 COMPLEX :: ONE, MONE, ZERO
287 COMPLEX :: MULT1, MULT2, A11, DETPIV, A22, A12
288 PARAMETER (ONE=(1.0e0,0.0e0), mone=(-1.0e0,0.0e0))
289 parameter(zero=(0.0e0,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 ctrsm('R', 'L', 'T', 'N', m, n, one,
304 & a(poselt_local), nfront,
305 & lr_block_ptr(1,1), m)
306 ELSE
307 CALL ctrsm('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 & 'CMUMPS_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 cscal(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 cmumps_lrtrsm
354 SUBROUTINE cmumps_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 CMUMPS_LRGEMM4)
358 TYPE(lrb_type),INTENT(IN) :: lrb
359 INTEGER(8), intent(in) :: la
360 COMPLEX, intent(inout) :: a(la)
361 COMPLEX, intent(inout), DIMENSION(:,:) :: SCALED
362 INTEGER,INTENT(IN) :: LD_DIAG, NFRONT, IW2(*)
363 INTEGER(8), INTENT(IN) :: POSELTT
364 COMPLEX, INTENT(IN), OPTIONAL :: DIAG(*)
365 INTEGER, INTENT(IN) :: MAXI_CLUSTER
366 COMPLEX, intent(inout) :: BLOCK(MAXI_CLUSTER)
367 INTEGER :: J, NROWS
368 COMPLEX :: 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 cmumps_lrgemm_scaling
393 SUBROUTINE cmumps_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 COMPLEX, 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 COMPLEX, 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 COMPLEX :: ALPHA,BETA
419 LOGICAL, INTENT(OUT) :: BUILDQ
420 COMPLEX, 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 COMPLEX, POINTER, DIMENSION(:,:) :: XY_YZ
427 COMPLEX, ALLOCATABLE, TARGET, DIMENSION(:,:) :: XQ, R_Y
428 COMPLEX, 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 COMPLEX, ALLOCATABLE :: WORK_RRQR(:), TAU_RRQR(:),
435 & y_rrqr(:,:)
436 INTEGER, ALLOCATABLE :: JPVT_RRQR(:)
437 INTEGER :: allocok, MREQ
438 REAL, EXTERNAL ::scnrm2
439 COMPLEX :: ONE, MONE, ZERO
440 parameter(one=(1.0e0,0.0e0), mone=(-1.0e0,0.0e0))
441 parameter(zero=(0.0e0,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 cmumps_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 cgemm('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 cmumps_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 cungqr
540 CALL cungqr
541 & (lrb1%K, rank, rank, y_rrqr(1,1),
542 & lrb1%K, tau_rrqr(1),
543 & work_rrqr(1), lwork, info )
544 CALL cgemm('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.AND..NOT. IF (LRB1%ISLR(LRB2%ISLR)) THEN
575.EQ. IF (LRB1%K0) 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.EQ. IF (SYM 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 CMUMPS_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.NOT..AND. IF ((LRB1%ISLR)LRB2%ISLR) THEN
604.EQ. IF (LRB2%K0) 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.EQ. IF (SYM 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 CMUMPS_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.NOT..AND..NOT. IF ((LRB1%ISLR)(LRB2%ISLR)) THEN
633.EQ. IF (SYM 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 CMUMPS_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.NOT. IF (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.GT. IF (SAVE_K+K_YZMAXI_RANK) THEN
672 write(*,*) 'internal error in cmumps_lrgemm4 1a',
673 & 'k_acc+k_cur>k_max:',SAVE_K,K_YZ,MAXI_RANK
674 CALL MUMPS_ABORT()
675 ENDIF
676.NE. IF (LRB3%MLRB1%M) THEN
677 write(*,*) 'internal error in cmumps_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 cgemm('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.NOT. IF (LUA_ACTIVATED) THEN
691 CALL cgemm('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.NOT. IF (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.GT. IF (SAVE_K+K_XYMAXI_RANK) THEN
706 write(*,*) 'internal error in cmumps_lrgemm4 2a',
707 & 'k_acc+k_cur>k_max:',SAVE_K,K_XY,MAXI_RANK
708 CALL MUMPS_ABORT()
709 ENDIF
710.NE. IF (LRB3%NLRB2%M) THEN
711 write(*,*) 'internal error in cmumps_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 cgemm(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.NOT. IF (LUA_ACTIVATED) THEN
725 CALL cgemm('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 cgemm('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.NOT..AND..NOT. IF ((LRB1%ISLR)(LRB2%ISLR)) THEN
744.NE. IF (SYM 0) deallocate(X)
745.NOT..AND. ELSEIF ((LRB1%ISLR)LRB2%ISLR) THEN
746.NE. IF (SYM 0) deallocate(Y)
747.AND..NOT. ELSEIF (LRB1%ISLR(LRB2%ISLR)) THEN
748.NE. IF (SYM 0) deallocate(Y)
749 ELSE
750.NE. IF (SYM 0) deallocate(Y1)
751.GE..AND. IF ((MIDBLK_COMPRESS1)BUILDQ) THEN
752 deallocate(XQ)
753 deallocate(R_Y)
754 ELSE
755 deallocate(Y)
756 ENDIF
757 ENDIF
758 1200 CONTINUE
759 END SUBROUTINE CMUMPS_LRGEMM4
760 SUBROUTINE CMUMPS_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 COMPLEX, 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 COMPLEX :: ONE, MONE, ZERO
773 PARAMETER (ONE=(1.0E0,0.0E0), MONE=(-1.0E0,0.0E0))
774 PARAMETER (ZERO=(0.0E0,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 cgemm('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 CMUMPS_DECOMPRESS_ACC
785 SUBROUTINE CMUMPS_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 COMPLEX, 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 COMPLEX, 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 COMPLEX :: ONE, MONE, ZERO
804 PARAMETER (ONE=(1.0E0,0.0E0), MONE=(-1.0E0,0.0E0))
805 PARAMETER (ZERO=(0.0E0,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 CMUMPS_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.LT. IF(JRANK) ACC_LRB%R(MIN(RANK,J)+1:
835 & RANK,JPVT_RRQR(J))= ZERO
836 END DO
837 CALL cungqr
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 & cmumps_compress_fr_updates: ',
860 & 'not enough memory? memory requested = ' , MREQ
861 CALL MUMPS_ABORT()
862 RETURN
863 END SUBROUTINE CMUMPS_COMPRESS_FR_UPDATES
864 SUBROUTINE CMUMPS_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 COMPLEX, 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 COMPLEX, ALLOCATABLE :: WORK_RRQR(:), TAU_RRQR(:)
880 COMPLEX, 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 COMPLEX :: ONE, MONE, ZERO
889 PARAMETER (ONE=(1.0E0,0.0E0), MONE=(-1.0E0,0.0E0))
890 PARAMETER (ZERO=(0.0E0,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 CMUMPS_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.EQ. SKIP2 = K0
913 ENDIF
914.AND. IF (SKIP1SKIP2) 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 CMUMPS_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.LT. IF(JRANK1) R1(MIN(RANK1,J)+1:
949 & RANK1,JPVT_RRQR(J))= ZERO
950 END DO
951 CALL cungqr
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 CMUMPS_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.LT. IF(JRANK2) R2(MIN(RANK2,J)+1:
981 & RANK2,JPVT_RRQR(J))= ZERO
982 END DO
983 CALL cungqr
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.OR. IF (BUILDQ1BUILDQ2) 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 CMUMPS_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.LT. IF (IFLAG0) GOTO 100
1019 CALL UPD_FLOP_UPDATE(LRB1, LRB2,
1020 & MIDBLK_COMPRESS-1, RANK, BUILDQ,
1021 & .TRUE., .FALSE., REC_ACC=.TRUE.)
1022 ENDIF
1023.NOT. IF ( SKIP1)
1024 & CALL UPD_FLOP_COMPRESS(LRB1, REC_ACC=.TRUE.)
1025.NOT. IF ( 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.AND..GT. IF (SKIP1(RANK20)) 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 & cmumps_recompress_acc: ',
1043 & 'not enough memory? memory requested = ' , MREQ
1044 CALL MUMPS_ABORT()
1045 RETURN
1046 END SUBROUTINE CMUMPS_RECOMPRESS_ACC
1047 RECURSIVE SUBROUTINE CMUMPS_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 COMPLEX, 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 COMPLEX :: ONE, MONE, ZERO
1072 PARAMETER (ONE=(1.0E0,0.0E0), MONE=(-1.0E0,0.0E0))
1073 PARAMETER (ZERO=(0.0E0,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.NE. IF (NB_NODES_NEW*NARYNB_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 cmumps_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.GE. IF (IMAX2) THEN
1096 DO I=2,IMAX
1097.NE. IF (POS_LIST(IOFF+I)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.NOT..OR..EQ. IF (RESORTLEVEL0) 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.GT. IF (NEW_ACC_RANK0) THEN
1118 CALL CMUMPS_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.GT. IF (NB_NODES_NEW1) 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 cmumps_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.EQ. IF (LEVEL0) 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.GT. IF (LEVEL0) THEN
1161 CALL DEALLOC_LRB(ACC_TMP, KEEP8, 4)
1162 ENDIF
1163 CALL CMUMPS_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 CMUMPS_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.NE. IF (POS_LIST_NEW(1)1) THEN
1178 write(*,*) 'internal error in ',
1179 & 'cmumps_recompress_acc_narytree', POS_LIST_NEW(1)
1180 ENDIF
1181 ACC_LRB%K = RANK_LIST_NEW(1)
1182.AND..GT. IF (RESORTLEVEL0) 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 CMUMPS_RECOMPRESS_ACC_NARYTREE
1196 SUBROUTINE CMUMPS_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 COMPLEX, 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 COMPLEX, ALLOCATABLE :: WORK_RRQR(:), TAU_RRQR(:)
1211 COMPLEX, 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 COMPLEX :: ONE, MONE, ZERO
1219 PARAMETER (ONE=(1.0E0,0.0E0), MONE=(-1.0E0,0.0E0))
1220 PARAMETER (ZERO=(0.0E0,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 cgemm('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 cgemm('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 CMUMPS_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 cgemm('n', 't', k1, n, k, one, proj(1,1), k1,
1263 & q2(1,1), n, one, acc_lrb%R(1,1), maxi_rank)
1264 IF (rank1.GT.0) 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 IF(j.LT.rank1) r1(min(rank1,j)+1:
1274 & rank1,jpvt_rrqr(j))= zero
1275 END DO
1276 CALL cungqr
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 cgemm('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 & CMUMPS_RECOMPRESS_ACC_V2: ',
1299 & 'not enough memory? memory requested = ' , mreq
1300 CALL mumps_abort()
1301 RETURN
1302 END SUBROUTINE cmumps_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 IF (cut(i+1) - cut(i) .GE. maxi_cluster) THEN
1311 maxi_cluster = cut(i+1) - cut(i)
1312 END IF
1313 END DO
1314 END SUBROUTINE max_cluster
1315 SUBROUTINE cmumps_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 IF ((sym.NE.0).AND.(fs_or_cb.EQ.0).AND.(j.NE.0)) THEN
1339 write(6,*) 'Internal error in CMUMPS_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 IF (fs_or_cb.EQ.0) THEN ! FS
1347 IF (j.EQ.0) 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 IF (k474.GE.2) THEN
1361 ind_u = k
1362 ENDIF
1363 ENDIF
1364 CALL cmumps_blr_retrieve_panel_loru(
1365 & iwhandler,
1366 & 0, ! L Panel
1367 & k, blr_l)
1368 IF (sym.EQ.0) THEN
1369 IF (lbandslave.AND.k474.GE.2) THEN
1370 blr_u => blr_u_col
1371 ELSE
1372 CALL cmumps_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 cmumps_get_lua_order
1397 SUBROUTINE cmumps_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 COMPLEX :: 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 COMPLEX, 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 COMPLEX :: ONE, ZERO
1434 parameter(one=(1.0e0,0.0e0))
1435 parameter(zero=(0.0e0,0.0e0))
1436 CALL cmumps_blr_retrieve_begsblr_dyn(iwhandler,
1437 & begs_blr_dynamic)
1438 CALL cmumps_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 IF (sym.EQ.0) 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 IF (sym.EQ.0) 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 IF (i.EQ.nb_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 IF (lrb%ISLR.AND.lrb%K.EQ.0) 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 IF (allocok.GT.0) THEN
1489 write(*,*) 'Not enough memory in CMUMPS_BLR_ASM_NIV1',
1490 & ", Memory requested = ", son_la
1491 CALL mumps_abort()
1492 ENDIF
1493C decompress block
1494 IF (lrb%ISLR) THEN
1495 CALL cgemm('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 IF (i.EQ.j.AND.sym.NE.0) THEN
1501C Diag block and LDLT, copy only lower half
1502 IF (j-nb_inasm.EQ.1.AND.nelim.GT.0) 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 IF (sym.NE.0.AND.j-nb_inasm.EQ.1.AND.nelim.GT.0) 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 IF (son_iw(kk+k1-1).LE.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 IF (i.EQ.j.AND.sym.NE.0) 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 cmumps_blr_free_cb_lrb(iwhandler,
1596C Only CB_LRB structure is left to deallocate
1597 & .true., keep8, keep(34))
1598 IF ((keep(486).EQ.3).OR.keep(486).EQ.0) THEN
1599C Case of FR solve: the BLR structure could not be freed
1600C in CMUMPS_END_FACTO_SLAVE and should be freed here
1601C Not reachable in case of error: set INFO1 to 0
1602 CALL cmumps_blr_end_front(iwhandler, 0, keep8, keep(34),
1603 & mtk405=keep(405))
1604 ENDIF
1605 END SUBROUTINE cmumps_blr_asm_niv1
1606 END MODULE cmumps_lr_core
1607C --------------------------------------------------------------------
1608 SUBROUTINE cmumps_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 COMPLEX :: A(LDA,*), TAU(*)
1685 COMPLEX :: 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 COMPLEX :: AKK
1694 LOGICAL INADMISSIBLE
1695 REAL, PARAMETER :: RZERO=0.0e+0, rone=1.0e+0
1696 COMPLEX :: ZERO
1697 COMPLEX :: ONE
1698 PARAMETER ( ONE = ( 1.0e+0, 0.0e+0 ) )
1699 parameter( zero = ( 0.0e+0, 0.0e+0 ) )
1700 REAL :: slamch
1701 INTEGER :: ilaenv, isamax
1702 EXTERNAL :: isamax, slamch
1703 EXTERNAL cgeqrf, cunmqr, xerbla
1704 EXTERNAL ilaenv
1705 EXTERNAL cgemm, cgemv, clarfg, cswap
1706 REAL, EXTERNAL :: scnrm2
1707 REAL, EXTERNAL :: snrm2
1708 info = 0
1709 islr = .false.
1710 IF( m.LT.0 ) THEN
1711 info = -1
1712 ELSE IF( n.LT.0 ) THEN
1713 info = -2
1714 ELSE IF( lda.LT.max( 1, m ) ) THEN
1715 info = -4
1716 END IF
1717 IF( info.EQ.0 ) THEN
1718 IF( ldw.LT.n ) THEN
1719 info = -8
1720 END IF
1721 END IF
1722 IF( info.NE.0 ) THEN
1723 WRITE(*,999) -info
1724 RETURN
1725 END IF
1726 minmn = min(m,n)
1727 IF( minmn.EQ.0 ) THEN
1728 rank = 0
1729 RETURN
1730 END IF
1731 nb = ilaenv( inb, 'CGEQRF', ' ', m, n, -1, -1 )
1732 SELECT CASE(abs(tol_opt))
1733 CASE(1)
1734 toleps_eff = toleps
1735 CASE(2)
1736C TOLEPS_EFF will be computed at step K=1 below
1737 CASE DEFAULT
1738 write(*,*) 'Internal error in CMUMPS_TRUNCATED_RRQR: TOL_OPT =',
1739 & tol_opt
1740 CALL mumps_abort()
1741 END SELECT
1742 toleps_eff = toleps
1743C
1744C Avoid pointers (and TARGET attribute on RWORK/WORK)
1745C because of implicit interface. An implicit interface
1746C is needed to avoid intermediate array copies
1747C VN1 => RWORK(1:N)
1748C VN2 => RWORK(N+1:2*N)
1749C AUXV => WORK(1:LDW,1:1)
1750C F => WORK(1:LDW,2:NB+1)
1751C LDF = LDW
1752* Initialize partial column norms. The first N elements of work
1753* store the exact column norms.
1754 DO j = 1, n
1755C VN1( J ) = scnrm2( M, A( 1, J ), 1 )
1756 rwork( j ) = scnrm2( m, a( 1, j ), 1 )
1757C VN2( J ) = VN1( J )
1758 rwork( n + j ) = rwork( j )
1759 jpvt(j) = j
1760 END DO
1761 IF (tol_opt.LT.0) THEN
1762C Compute TRUNC_ERR for first step
1763C TRUNC_ERR = snrm2( N, VN1( 1 ), 1 )
1764 trunc_err = snrm2( n, rwork( 1 ), 1 )
1765 ENDIF
1766 offset = 0
1767 tol3z = sqrt(slamch('Epsilon'))
1768 DO
1769 jb = min(nb,minmn-offset)
1770 lsticc = 0
1771 k = 0
1772 DO
1773 IF(k.EQ.jb) EXIT
1774 k = k+1
1775 rk = offset+k
1776C PVT = ( RK-1 ) + ISAMAX( N-RK+1, VN1( RK ), 1 )
1777 pvt = ( rk-1 ) + isamax( n-rk+1, rwork( rk ), 1 )
1778 IF (rk.EQ.1) THEN
1779C IF (abs(TOL_OPT).EQ.2) TOLEPS_EFF = VN1(PVT)*TOLEPS
1780 IF (abs(tol_opt).EQ.2) toleps_eff = rwork(pvt)*toleps
1781 ENDIF
1782 IF (tol_opt.GT.0) THEN
1783C TRUNC_ERR = VN1(PVT)
1784 trunc_err = rwork(pvt)
1785C ELSE
1786C TRUNC_ERR has been already computed at previous step
1787 ENDIF
1788 IF(trunc_err.LT.toleps_eff) THEN
1789 rank = rk-1
1790 islr = .true.
1791 RETURN
1792 ENDIF
1793 inadmissible = (rk.GT.maxrank)
1794 IF (inadmissible) THEN
1795 rank = rk
1796 info = rk
1797 islr = .false.
1798 RETURN
1799 END IF
1800 IF( pvt.NE.rk ) THEN
1801 CALL cswap( m, a( 1, pvt ), 1, a( 1, rk ), 1 )
1802c CALL cswap( K-1, F( PVT-OFFSET, 1 ), LDF,
1803c & F( K, 1 ), LDF )
1804 CALL cswap( k-1, work( pvt-offset, 2 ), ldw,
1805 & work( k, 2 ), ldw )
1806 itemp = jpvt(pvt)
1807 jpvt(pvt) = jpvt(rk)
1808 jpvt(rk) = itemp
1809C VN1(PVT) = VN1(RK)
1810C VN2(PVT) = VN2(RK)
1811 rwork(pvt) = rwork(rk)
1812 rwork(n+pvt) = rwork(n+rk)
1813 END IF
1814* Apply previous Householder reflectors to column K:
1815* A(RK:M,RK) := A(RK:M,RK) - A(RK:M,OFFSET+1:RK-1)*F(K,1:K-1)**H.
1816 IF( k.GT.1 ) THEN
1817 DO j = 1, k-1
1818C F( K, J ) = CONJG( F( K, J ) )
1819 work( k, j+1 ) = conjg( work( k, j+1 ) )
1820 END DO
1821 CALL cgemv( 'No transpose', m-rk+1, k-1, -one,
1822C & A(RK,OFFSET+1), LDA, F(K,1), LDF,
1823 & a(rk,offset+1), lda, work(k,2), ldw,
1824 & one, a(rk,rk), 1 )
1825 DO j = 1, k - 1
1826C F( K, J ) = CONJG( F( K, J ) )
1827 work( k, j + 1 ) = conjg( work( k, j + 1 ) )
1828 END DO
1829 END IF
1830* Generate elementary reflector H(k).
1831 IF( rk.LT.m ) THEN
1832 CALL clarfg( m-rk+1, a(rk,rk), a(rk+1,rk), 1, tau(rk) )
1833 ELSE
1834 CALL clarfg( 1, a(rk,rk), a(rk,rk), 1, tau(rk) )
1835 END IF
1836 akk = a(rk,rk)
1837 a(rk,rk) = one
1838* Compute Kth column of F:
1839* F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**H*A(RK:M,K).
1840 IF( rk.LT.n ) THEN
1841 CALL cgemv( 'Conjugate transpose', m-rk+1, n-rk, tau(rk),
1842 & a(rk,rk+1), lda, a(rk,rk), 1, zero,
1843C & F( K+1, K ), 1 )
1844 & work( k+1, k+1 ), 1 )
1845 END IF
1846* Padding F(1:K,K) with zeros.
1847 DO j = 1, k
1848C F( J, K ) = ZERO
1849 work( j, k+1 ) = zero
1850 END DO
1851* Incremental updating of F:
1852* F(1:N,K) := F(1:N-OFFSET,K) -
1853* tau(RK)*F(1:N,1:K-1)*A(RK:M,OFFSET+1:RK-1)**H*A(RK:M,RK).
1854 IF( k.GT.1 ) THEN
1855 CALL cgemv( 'conjugate transpose', M-RK+1, K-1, -TAU(RK),
1856 & A(RK,OFFSET+1), LDA, A(RK,RK), 1, ZERO,
1857 & WORK(1,1), 1 )
1858C & AUXV(1,1), 1 )
1859 CALL cgemv( 'no transpose', N-OFFSET, K-1, ONE,
1860 & WORK(1,2), LDW, WORK(1,1), 1, ONE, WORK(1,K+1), 1 )
1861C & F(1,1), LDF, AUXV(1,1), 1, ONE, F(1,K), 1 )
1862 END IF
1863* Update the current row of A:
1864* A(RK,RK+1:N) := A(RK,RK+1:N) - A(RK,OFFSET+1:RK)*F(K+1:N,1:K)**H.
1865.LT. IF( RKN ) THEN
1866 CALL cgemm( 'no transpose', 'conjugate transpose',
1867 & 1, N-RK,
1868C & K, -ONE, A( RK, OFFSET+1 ), LDA, F( K+1, 1 ), LDF,
1869 & K, -ONE, A( RK, OFFSET+1 ), LDA, WORK( K+1,2 ), LDW,
1870 & ONE, A( RK, RK+1 ), LDA )
1871 END IF
1872* Update partial column norms.
1873*
1874.LT. IF( RKMINMN ) THEN
1875 DO J = RK + 1, N
1876C IF( VN1( J ).NE.RZERO ) THEN
1877.NE. IF( RWORK( J )RZERO ) THEN
1878*
1879* NOTE: The following 4 lines follow from the analysis in
1880* Lapack Working Note 176.
1881*
1882C TEMP = ABS( A( RK, J ) ) / VN1( J )
1883 TEMP = ABS( A( RK, J ) ) / RWORK( J )
1884 TEMP = MAX( RZERO, ( RONE+TEMP )*( RONE-TEMP ) )
1885C TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
1886 TEMP2 = TEMP*( RWORK( J ) / RWORK( N+J ) )**2
1887.LE. IF( TEMP2 TOL3Z ) THEN
1888C VN2( J ) = REAL( LSTICC )
1889 RWORK( N+J ) = REAL( LSTICC )
1890 LSTICC = J
1891 ELSE
1892C VN1( J ) = VN1( J )*SQRT( TEMP )
1893 RWORK( J ) = RWORK( J )*SQRT( TEMP )
1894 END IF
1895 END IF
1896 END DO
1897 END IF
1898 A( RK, RK ) = AKK
1899.NE. IF (LSTICC0) EXIT
1900.LT. IF (TOL_OPT0) THEN
1901C Compute TRUNC_ERR for next step
1902C TRUNC_ERR = snrm2( N-RK, VN1( RK+1 ), 1 )
1903 TRUNC_ERR = snrm2( N-RK, RWORK( RK+1 ), 1 )
1904 ENDIF
1905 END DO
1906* Apply the block reflector to the rest of the matrix:
1907* A(RK+1:M,RK+1:N) := A(RK+1:M,RK+1:N) -
1908* A(RK+1:M,OFFSET+1:RK)*F(K+1:N-OFFSET,1:K)**H.
1909.LT. IF( RKMIN(N,M) ) THEN
1910 CALL cgemm( 'no transpose', 'conjugate transpose', M-RK,
1911 & N-RK, K, -ONE, A(RK+1,OFFSET+1), LDA,
1912C & F(K+1,1), LDF, ONE, A(RK+1,RK+1), LDA )
1913 & WORK(K+1,2), LDW, ONE, A(RK+1,RK+1), LDA )
1914 END IF
1915* Recomputation of difficult columns.
1916.GT. DO WHILE( LSTICC0 )
1917C ITEMP = NINT( VN2( LSTICC ) )
1918 ITEMP = NINT( RWORK( N + LSTICC ) )
1919C VN1( LSTICC ) = scnrm2( M-RK, A( RK+1, LSTICC ), 1 )
1920 RWORK( LSTICC ) = scnrm2( M-RK, A( RK+1, LSTICC ), 1 )
1921*
1922* NOTE: The computation of RWORK( LSTICC ) relies on the fact that
1923* SNRM2 does not fail on vectors with norm below the value of
1924* SQRT(DLAMCH('S'))
1925*
1926C VN2( LSTICC ) = VN1( LSTICC )
1927 RWORK( N + LSTICC ) = RWORK( LSTICC )
1928 LSTICC = ITEMP
1929 END DO
1930.GE. IF(RKMINMN) EXIT
1931 OFFSET = RK
1932.LT. IF (TOL_OPT0) THEN
1933C Compute TRUNC_ERR for next step
1934C TRUNC_ERR = snrm2( N-RK, VN1( RK+1 ), 1 )
1935 TRUNC_ERR = snrm2( N-RK, RWORK( RK+1 ), 1 )
1936 ENDIF
1937 END DO
1938 RANK = RK
1939.NOT..GT. ISLR = (RKMAXRANK)
1940 RETURN
1941 999 FORMAT ('on entry to cmumps_truncated_rrqr, parameter number',
1942 & I2,' had an illegal value')
1943 END SUBROUTINE CMUMPS_TRUNCATED_RRQR
mumps_abort
#define mumps_abort
Definition VE_Metis.h:25
cmumps_truncated_rrqr
subroutine cmumps_truncated_rrqr(m, n, a, lda, jpvt, tau, work, ldw, rwork, toleps, tol_opt, rank, maxrank, info, islr)
Definition clr_core.F:1611
clarfg
subroutine clarfg(n, alpha, x, incx, tau)
CLARFG generates an elementary reflector (Householder matrix).
Definition clarfg.f:106
cungqr
subroutine cungqr(m, n, k, a, lda, tau, work, lwork, info)
CUNGQR
Definition cungqr.f:128
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
cgemv
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:158
ctrsm
subroutine ctrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRSM
Definition ctrsm.f:180
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:187
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
cmumps_lr_core
Definition clr_core.F:18
cmumps_lr_core::cmumps_recompress_acc_narytree
recursive subroutine cmumps_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 clr_core.F:1053
cmumps_lr_core::cmumps_recompress_acc_v2
subroutine cmumps_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 clr_core.F:1200
cmumps_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 clr_core.F:45
cmumps_lr_core::cmumps_compress_fr_updates
subroutine cmumps_compress_fr_updates(acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, nfront, niv, toleps, tol_opt, kpercent, buildq, loru, cb_compress)
Definition clr_core.F:788
cmumps_lr_core::cmumps_blr_asm_niv1
subroutine cmumps_blr_asm_niv1(a, la, posel1, nfront, nass1, iwhandler, son_iw, liw, lstk, nelim, k1, k2, sym, keep, keep8, opassw)
Definition clr_core.F:1400
cmumps_lr_core::cmumps_recompress_acc
subroutine cmumps_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 clr_core.F:868
cmumps_lr_core::alloc_lrb_from_acc
subroutine alloc_lrb_from_acc(acc_lrb, lrb_out, k, m, n, loru, iflag, ierror, keep8)
Definition clr_core.F:160
cmumps_lr_core::cmumps_get_lua_order
subroutine cmumps_get_lua_order(nb_blocks, order, rank, iwhandler, sym, fs_or_cb, i, j, frfr_updates, lbandslave_in, k474, blr_u_col)
Definition clr_core.F:1318
cmumps_lr_core::max_cluster
subroutine max_cluster(cut, cut_size, maxi_cluster)
Definition clr_core.F:1304
cmumps_lr_core::cmumps_lrtrsm
subroutine cmumps_lrtrsm(a, la, poselt_local, nfront, lda, lrb, niv, sym, loru, iw, offset_iw)
Definition clr_core.F:270
cmumps_lr_core::alloc_lrb
subroutine alloc_lrb(lrb_out, k, m, n, islr, iflag, ierror, keep8)
Definition clr_core.F:111
cmumps_lr_core::init_lrb
subroutine init_lrb(lrb_out, k, m, n, islr)
Definition clr_core.F:27
cmumps_lr_core::cmumps_lrgemm4
subroutine cmumps_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 clr_core.F:406
cmumps_lr_core::cmumps_lrgemm_scaling
subroutine cmumps_lrgemm_scaling(lrb, scaled, a, la, diag, ld_diag, iw2, poseltt, nfront, block, maxi_cluster)
Definition clr_core.F:356
cmumps_lr_core::regrouping2
subroutine regrouping2(cut, npartsass, nass, npartscb, ncb, ibcksz, onlycb, k472)
Definition clr_core.F:184
cmumps_lr_data_m
Definition cmumps_lr_data_m.F:14
cmumps_lr_data_m::cmumps_blr_free_cb_lrb
subroutine, public cmumps_blr_free_cb_lrb(iwhandler, free_only_struct, keep8, k34)
Definition cmumps_lr_data_m.F:970
cmumps_lr_data_m::cmumps_blr_retrieve_cb_lrb
subroutine, public cmumps_blr_retrieve_cb_lrb(iwhandler, thecb)
Definition cmumps_lr_data_m.F:837
cmumps_lr_data_m::cmumps_blr_retrieve_panel_loru
subroutine, public cmumps_blr_retrieve_panel_loru(iwhandler, loru, ipanel, thelrbpanel)
Definition cmumps_lr_data_m.F:750
cmumps_lr_data_m::cmumps_blr_end_front
subroutine, public cmumps_blr_end_front(iwhandler, info1, keep8, k34, lrsolve_act_opt, mtk405)
Definition cmumps_lr_data_m.F:342
cmumps_lr_data_m::cmumps_blr_retrieve_begsblr_dyn
subroutine, public cmumps_blr_retrieve_begsblr_dyn(iwhandler, begs_blr_dynamic)
Definition cmumps_lr_data_m.F:638
cmumps_lr_stats
Definition clr_stats.F:14
cmumps_lr_stats::upd_flop_trsm
subroutine upd_flop_trsm(lrb, loru)
Definition clr_stats.F:199
cmumps_lr_stats::upd_flop_decompress
subroutine upd_flop_decompress(f, cb)
Definition clr_stats.F:140
cmumps_lr_type
Definition clr_type.F:14
cmumps_lr_type::dealloc_lrb
subroutine dealloc_lrb(lrb_out, keep8, k34)
Definition clr_type.F:25
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
cmumps_lr_type::lrb_type
Definition clr_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