OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
smumps_lr_core Module Reference

Functions/Subroutines

subroutine init_lrb (lrb_out, k, m, n, islr)
subroutine is_front_blr_candidate (inode, niv, nfront, nass, blron, k489, k490, k491, k492, k20, k60, idad, k38, lrstatus, n, lrgroups)
subroutine alloc_lrb (lrb_out, k, m, n, islr, iflag, ierror, keep8)
subroutine alloc_lrb_from_acc (acc_lrb, lrb_out, k, m, n, loru, iflag, ierror, keep8)
subroutine regrouping2 (cut, npartsass, nass, npartscb, ncb, ibcksz, onlycb, k472)
subroutine smumps_lrtrsm (a, la, poselt_local, nfront, lda, lrb, niv, sym, loru, iw, offset_iw)
subroutine smumps_lrgemm_scaling (lrb, scaled, a, la, diag, ld_diag, iw2, poseltt, nfront, block, maxi_cluster)
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)
subroutine smumps_decompress_acc (acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, nfront, niv, loru, count_flops)
subroutine smumps_compress_fr_updates (acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, nfront, niv, toleps, tol_opt, kpercent, buildq, loru, cb_compress)
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)
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)
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)
subroutine max_cluster (cut, cut_size, maxi_cluster)
subroutine smumps_get_lua_order (nb_blocks, order, rank, iwhandler, sym, fs_or_cb, i, j, frfr_updates, lbandslave_in, k474, blr_u_col)
subroutine smumps_blr_asm_niv1 (a, la, posel1, nfront, nass1, iwhandler, son_iw, liw, lstk, nelim, k1, k2, sym, keep, keep8, opassw)

Function/Subroutine Documentation

◆ alloc_lrb()

subroutine smumps_lr_core::alloc_lrb ( type(lrb_type), intent(out) lrb_out,
integer, intent(in) k,
integer, intent(in) m,
integer, intent(in) n,
logical, intent(in) islr,
integer, intent(inout) iflag,
integer, intent(inout) ierror,
integer(8), dimension(150) keep8 )

Definition at line 110 of file slr_core.F.

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
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

◆ alloc_lrb_from_acc()

subroutine smumps_lr_core::alloc_lrb_from_acc ( type(lrb_type), intent(in) acc_lrb,
type(lrb_type), intent(out) lrb_out,
integer, intent(in) k,
integer, intent(in) m,
integer, intent(in) n,
integer, intent(in) loru,
integer, intent(inout) iflag,
integer, intent(inout) ierror,
integer(8), dimension(150) keep8 )

Definition at line 158 of file slr_core.F.

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

◆ init_lrb()

subroutine smumps_lr_core::init_lrb ( type(lrb_type), intent(out) lrb_out,
integer, intent(in) k,
integer, intent(in) m,
integer, intent(in) n,
logical, intent(in) islr )

Definition at line 26 of file slr_core.F.

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)

◆ is_front_blr_candidate()

subroutine smumps_lr_core::is_front_blr_candidate ( integer, intent(in) inode,
integer, intent(in) niv,
integer, intent(in) nfront,
integer, intent(in) nass,
integer, intent(in) blron,
integer, intent(in) k489,
integer, intent(in) k490,
integer, intent(in) k491,
integer, intent(in) k492,
integer, intent(in) k20,
integer, intent(in) k60,
integer, intent(in) idad,
integer, intent(in) k38,
integer, intent(out) lrstatus,
integer, intent(in) n,
integer, dimension(n), intent(in), optional lrgroups )

Definition at line 41 of file slr_core.F.

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

◆ max_cluster()

subroutine smumps_lr_core::max_cluster ( integer, dimension(:), pointer cut,
integer, intent(in) cut_size,
integer, intent(out) maxi_cluster )

Definition at line 1303 of file slr_core.F.

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

◆ regrouping2()

subroutine smumps_lr_core::regrouping2 ( integer, dimension(:), pointer cut,
integer, intent(inout) npartsass,
integer, intent(in) nass,
integer, intent(inout) npartscb,
integer, intent(in) ncb,
integer, intent(in) ibcksz,
logical onlycb,
integer, intent(in) k472 )

Definition at line 182 of file slr_core.F.

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)
max
#define max(a, b)
Definition macros.h:21

◆ smumps_blr_asm_niv1()

subroutine smumps_lr_core::smumps_blr_asm_niv1 ( real, dimension(la) a,
integer(8) la,
integer(8) posel1,
integer nfront,
integer nass1,
integer iwhandler,
integer, dimension(:) son_iw,
integer liw,
integer lstk,
integer nelim,
integer k1,
integer k2,
integer sym,
integer, dimension(500) keep,
integer(8), dimension(150) keep8,
double precision, intent(inout) opassw )

Definition at line 1397 of file slr_core.F.

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 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 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 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 smumps_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 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
mumps_abort
#define mumps_abort
Definition VE_Metis.h:25
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:187
min
#define min(a, b)
Definition macros.h:20

◆ smumps_compress_fr_updates()

subroutine smumps_lr_core::smumps_compress_fr_updates ( type(lrb_type), intent(inout) acc_lrb,
integer, intent(in) maxi_cluster,
integer, intent(in) maxi_rank,
real, dimension(la), intent(inout) a,
integer(8), intent(in) la,
integer(8), intent(in) poseltt,
integer, intent(in) nfront,
integer, intent(in) niv,
real, intent(in) toleps,
integer, intent(in) tol_opt,
integer, intent(in) kpercent,
logical, intent(out) buildq,
integer, intent(in) loru,
logical, intent(in) cb_compress )

Definition at line 785 of file slr_core.F.

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
sorgqr
subroutine sorgqr(m, n, k, a, lda, tau, work, lwork, info)
SORGQR
Definition sorgqr.f:128
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_decompress_acc()

subroutine smumps_lr_core::smumps_decompress_acc ( type(lrb_type), intent(inout) acc_lrb,
integer, intent(in) maxi_cluster,
integer, intent(in) maxi_rank,
real, dimension(la), intent(inout) a,
integer(8), intent(in) la,
integer(8), intent(in) poseltt,
integer, intent(in) nfront,
integer, intent(in) niv,
integer loru,
logical, optional count_flops )

Definition at line 760 of file slr_core.F.

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

◆ smumps_get_lua_order()

subroutine smumps_lr_core::smumps_get_lua_order ( integer, intent(in) nb_blocks,
integer, dimension(nb_blocks), intent(out) order,
integer, dimension(nb_blocks), intent(out) rank,
integer, intent(in) iwhandler,
integer, intent(in) sym,
integer, intent(in) fs_or_cb,
integer, intent(in) i,
integer, intent(in) j,
integer, intent(out) frfr_updates,
logical, intent(in), optional lbandslave_in,
integer, intent(in), optional k474,
type(lrb_type), dimension(:), optional, pointer blr_u_col )

Definition at line 1315 of file slr_core.F.

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 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 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 smumps_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 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)
mumps_sort_int
subroutine mumps_sort_int(n, val, id)
Definition tools_common.F:709

◆ smumps_lrgemm4()

subroutine smumps_lr_core::smumps_lrgemm4 ( real alpha,
type(lrb_type), intent(in) lrb1,
type(lrb_type), intent(in) lrb2,
real beta,
real, dimension(la), intent(inout) a,
integer(8), intent(in) la,
integer(8), intent(in) poseltt,
integer, intent(in) nfront,
integer, intent(in) sym,
integer, intent(inout) iflag,
integer, intent(inout) ierror,
integer, intent(in) midblk_compress,
real, intent(in) toleps,
integer, intent(in) tol_opt,
integer, intent(in) kpercent,
integer rank,
logical, intent(out) buildq,
logical, intent(in) lua_activated,
integer, intent(in), optional loru,
type(lrb_type), intent(inout), optional lrb3,
integer, intent(in), optional maxi_rank,
integer, intent(in), optional maxi_cluster,
real, dimension(*), intent(in), optional diag,
integer, intent(in), optional ld_diag,
integer, dimension(*), intent(in), optional iw2,
real, dimension(*), intent(inout), optional block )

Definition at line 393 of file slr_core.F.

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
alpha
#define alpha
Definition eval.h:35

◆ smumps_lrgemm_scaling()

subroutine smumps_lr_core::smumps_lrgemm_scaling ( type(lrb_type), intent(in) lrb,
real, dimension(:,:), intent(inout) scaled,
real, dimension(la), intent(inout) a,
integer(8), intent(in) la,
real, dimension(*), intent(in), optional diag,
integer, intent(in) ld_diag,
integer, dimension(*), intent(in) iw2,
integer(8), intent(in) poseltt,
integer, intent(in) nfront,
real, dimension(maxi_cluster), intent(inout) block,
integer, intent(in) maxi_cluster )

Definition at line 354 of file slr_core.F.

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

◆ smumps_lrtrsm()

subroutine smumps_lr_core::smumps_lrtrsm ( real, dimension(la), intent(inout) a,
integer(8), intent(in) la,
integer(8), intent(in) poselt_local,
integer, intent(in) nfront,
integer, intent(in) lda,
type(lrb_type), intent(inout) lrb,
integer, intent(in) niv,
integer, intent(in) sym,
integer, intent(in) loru,
integer, dimension(*), optional iw,
integer, optional offset_iw )

Definition at line 268 of file slr_core.F.

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)
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
strsm
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181

◆ smumps_recompress_acc()

subroutine smumps_lr_core::smumps_recompress_acc ( type(lrb_type), intent(inout) acc_lrb,
integer, intent(in) maxi_cluster,
integer, intent(in) maxi_rank,
real, dimension(la), intent(inout) a,
integer(8), intent(in) la,
integer(8), intent(in) poseltt,
integer, intent(in) nfront,
integer, intent(in) niv,
integer, intent(in) midblk_compress,
real, intent(in) toleps,
integer, intent(in) tol_opt,
integer, intent(in) kpercent_rmb,
integer, intent(in) kpercent_lua,
integer, intent(inout) new_acc_rank )

Definition at line 864 of file slr_core.F.

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

◆ smumps_recompress_acc_narytree()

recursive subroutine smumps_lr_core::smumps_recompress_acc_narytree ( type(lrb_type), intent(inout), target acc_lrb,
integer, intent(in) maxi_cluster,
integer, intent(in) maxi_rank,
real, dimension(la), intent(inout) a,
integer(8), intent(in) la,
integer(8), intent(in) poseltt,
integer(8), dimension(150) keep8,
integer, intent(in) nfront,
integer, intent(in) niv,
integer, intent(in) midblk_compress,
real, intent(in) toleps,
integer, intent(in) tol_opt,
integer, intent(in) kpercent_rmb,
integer, intent(in) kpercent_lua,
integer, intent(in) k478,
integer, dimension(nb_nodes), intent(inout) rank_list,
integer, dimension(nb_nodes), intent(inout) pos_list,
integer, intent(in) nb_nodes,
integer, intent(in) level,
type(lrb_type), intent(inout), optional, target acc_tmp )

Definition at line 1047 of file slr_core.F.

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)

◆ smumps_recompress_acc_v2()

subroutine smumps_lr_core::smumps_recompress_acc_v2 ( type(lrb_type), intent(inout) acc_lrb,
integer, intent(in) maxi_cluster,
integer, intent(in) maxi_rank,
real, dimension(la), intent(inout) a,
integer(8), intent(in) la,
integer(8), intent(in) poseltt,
integer, intent(in) nfront,
integer, intent(in) niv,
integer, intent(in) midblk_compress,
real, intent(in) toleps,
integer, intent(in) tol_opt,
integer, intent(in) kpercent_rmb,
integer, intent(in) kpercent_lua,
integer, intent(inout) new_acc_rank )

Definition at line 1196 of file slr_core.F.

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 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 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