1 /* split.c */
2
3 #include "../FrontMtx.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 ----------------------------------------------------------------
8 purpose -- for each U_{J,bnd{J}} matrix, remove from hash table,
9 split into their U_{J,K} submatrices and insert
10 into the hash table.
11
12 created -- 98may04, cca
13 ----------------------------------------------------------------
14 */
15 void
FrontMtx_splitUpperMatrices(FrontMtx * frontmtx,int msglvl,FILE * msgFile)16 FrontMtx_splitUpperMatrices (
17 FrontMtx *frontmtx,
18 int msglvl,
19 FILE *msgFile
20 ) {
21 SubMtx *mtxUJ, *mtxUJJ, *mtxUJK ;
22 SubMtxManager *manager ;
23 double *entUJ, *entUJK ;
24 int count, first, ii, inc1, inc2, jcol, jj, J, K, nbytes,
25 ncolJ, ncolUJ, ncolUJK, nentUJ, nentUJK, neqns, nfront,
26 nJ, nrowUJ, nrowUJK, offset, v ;
27 int *colindJ, *colindUJ, *colindUJK, *colmap, *indicesUJ,
28 *indicesUJK, *locmap, *rowindUJ, *rowindUJK, *sizesUJ,
29 *sizesUJK ;
30 I2Ohash *upperhash ;
31 /*
32 ---------------
33 check the input
34 ---------------
35 */
36 if ( frontmtx == NULL || (msglvl > 0 && msgFile == NULL) ) {
37 fprintf(stderr,
38 "\n fatal error in FrontMtx_splitUpperMatrices(%p,%d,%p)"
39 "\n bad input\n", frontmtx, msglvl, msgFile) ;
40 exit(-1) ;
41 }
42 nfront = FrontMtx_nfront(frontmtx) ;
43 neqns = FrontMtx_neqns(frontmtx) ;
44 upperhash = frontmtx->upperhash ;
45 manager = frontmtx->manager ;
46 /*
47 -----------------------------------
48 construct the column and local maps
49 -----------------------------------
50 */
51 colmap = IVinit(neqns, -1) ;
52 locmap = IVinit(neqns, -1) ;
53 for ( J = 0 ; J < nfront ; J++ ) {
54 if ( (nJ = FrontMtx_frontSize(frontmtx, J)) > 0 ) {
55 FrontMtx_columnIndices(frontmtx, J, &ncolJ, &colindJ) ;
56 if ( ncolJ > 0 && colindJ != NULL ) {
57 for ( ii = 0 ; ii < nJ ; ii++ ) {
58 v = colindJ[ii] ;
59 colmap[v] = J ;
60 locmap[v] = ii ;
61 }
62 }
63 }
64 }
65 if ( msglvl > 2 ) {
66 fprintf(msgFile, "\n\n colmap[]") ;
67 IVfprintf(msgFile, neqns, colmap) ;
68 fprintf(msgFile, "\n\n locmap[]") ;
69 IVfprintf(msgFile, neqns, locmap) ;
70 fflush(msgFile) ;
71 }
72 /*
73 ---------------------------------------------
74 move the U_{J,J} matrices into the hash table
75 ---------------------------------------------
76 */
77 for ( J = 0 ; J < nfront ; J++ ) {
78 if ( (mtxUJJ = FrontMtx_upperMtx(frontmtx, J, J)) != NULL ) {
79 I2Ohash_insert(frontmtx->upperhash, J, J, mtxUJJ) ;
80 }
81 }
82 /*
83 ------------------------------------------------------------
84 now split the U_{J,bnd{J}} matrices into U_{J,K} matrices.
85 note: columns of U_{J,bnd{J}} are assumed to be in ascending
86 order with respect to the column ordering of the matrix.
87 ------------------------------------------------------------
88 */
89 for ( J = 0 ; J < nfront ; J++ ) {
90 mtxUJ = FrontMtx_upperMtx(frontmtx, J, nfront) ;
91 if ( msglvl > 2 ) {
92 fprintf(msgFile, "\n\n ### J = %d, mtxUJ = %p", J, mtxUJ) ;
93 fflush(msgFile) ;
94 }
95 if ( mtxUJ != NULL ) {
96 if ( msglvl > 2 ) {
97 SubMtx_writeForHumanEye(mtxUJ, msgFile) ;
98 fflush(msgFile) ;
99 }
100 SubMtx_columnIndices(mtxUJ, &ncolUJ, &colindUJ) ;
101 SubMtx_rowIndices(mtxUJ, &nrowUJ, &rowindUJ) ;
102 if ( msglvl > 2 ) {
103 fprintf(msgFile, "\n column indices for J") ;
104 IVfprintf(msgFile, ncolUJ, colindUJ) ;
105 fprintf(msgFile, "\n row indices for UJ") ;
106 IVfprintf(msgFile, nrowUJ, rowindUJ) ;
107 fflush(msgFile) ;
108 }
109 if ( (K = colmap[colindUJ[0]]) == colmap[colindUJ[ncolUJ-1]] ) {
110 if ( msglvl > 2 ) {
111 fprintf(msgFile, "\n front %d supports only %d", J, K) ;
112 fflush(msgFile) ;
113 }
114 /*
115 -------------------------------------------------
116 U_{J,bnd{J}} is one submatrix, bnd{J} \subseteq K
117 set row and column indices and change column id
118 -------------------------------------------------
119 */
120 IVramp(nrowUJ, rowindUJ, 0, 1) ;
121 for ( ii = 0 ; ii < ncolUJ ; ii++ ) {
122 colindUJ[ii] = locmap[colindUJ[ii]] ;
123 }
124 SubMtx_setFields(mtxUJ, mtxUJ->type, mtxUJ->mode, J, K,
125 mtxUJ->nrow, mtxUJ->ncol, mtxUJ->nent) ;
126 /*
127 mtxUJ->colid = K ;
128 */
129 if ( msglvl > 2 ) {
130 fprintf(msgFile, "\n\n ## inserting U(%d,%d) ", J, K) ;
131 SubMtx_writeForHumanEye(mtxUJ, msgFile) ;
132 fflush(msgFile) ;
133 }
134 I2Ohash_insert(upperhash, J, K, (void *) mtxUJ) ;
135 } else {
136 /*
137 -----------------------------------
138 split U_{J,bnd{J}} into submatrices
139 -----------------------------------
140 */
141 nJ = FrontMtx_frontSize(frontmtx, J) ;
142 if ( SUBMTX_IS_DENSE_COLUMNS(mtxUJ) ) {
143 SubMtx_denseInfo(mtxUJ,
144 &nrowUJ, &ncolUJ, &inc1, &inc2, &entUJ) ;
145 } else if ( SUBMTX_IS_SPARSE_COLUMNS(mtxUJ) ) {
146 SubMtx_sparseColumnsInfo(mtxUJ, &ncolUJ, &nentUJ,
147 &sizesUJ, &indicesUJ, &entUJ) ;
148 offset = 0 ;
149 count = sizesUJ[0] ;
150 }
151 first = 0 ;
152 K = colmap[colindUJ[0]] ;
153 for ( jcol = 1 ; jcol <= ncolUJ ; jcol++ ) {
154 if ( msglvl > 2 ) {
155 fprintf(msgFile, "\n jcol = %d", jcol) ;
156 if ( jcol < ncolUJ ) {
157 fprintf(msgFile, ", colmap[%d] = %d",
158 colindUJ[jcol], colmap[colindUJ[jcol]]);
159 }
160 fflush(msgFile) ;
161 }
162 if ( jcol == ncolUJ || K != colmap[colindUJ[jcol]] ) {
163 ncolUJK = jcol - first ;
164 if ( SUBMTX_IS_DENSE_COLUMNS(mtxUJ) ) {
165 nentUJK = nJ*ncolUJK ;
166 } else if ( SUBMTX_IS_SPARSE_COLUMNS(mtxUJ) ) {
167 if ( count == 0 ) {
168 goto no_entries ;
169 }
170 nentUJK = count ;
171 }
172 nbytes = SubMtx_nbytesNeeded(mtxUJ->type, mtxUJ->mode,
173 nJ, ncolUJK, nentUJK) ;
174 if ( msglvl > 2 ) {
175 fprintf(msgFile,
176 "\n ncolUJK %d, nentUJK %d, nbytes %d",
177 ncolUJK, nentUJK, nbytes) ;
178 fflush(msgFile) ;
179 }
180 mtxUJK = SubMtxManager_newObjectOfSizeNbytes(manager,
181 nbytes) ;
182 SubMtx_init(mtxUJK, mtxUJ->type, mtxUJ->mode, J, K,
183 nJ, ncolUJK, nentUJK) ;
184 if ( SUBMTX_IS_DENSE_COLUMNS(mtxUJ) ) {
185 SubMtx_denseInfo(mtxUJK,
186 &nrowUJK, &ncolUJK, &inc1, &inc2, &entUJK) ;
187 if ( FRONTMTX_IS_REAL(frontmtx) ) {
188 DVcopy(nentUJK, entUJK, entUJ + first*nJ) ;
189 } else if ( FRONTMTX_IS_COMPLEX(frontmtx) ) {
190 DVcopy(2*nentUJK, entUJK, entUJ + 2*first*nJ) ;
191 }
192 } else if ( SUBMTX_IS_SPARSE_COLUMNS(mtxUJ) ) {
193 SubMtx_sparseColumnsInfo(mtxUJK, &ncolUJK, &nentUJK,
194 &sizesUJK, &indicesUJK, &entUJK) ;
195 IVcopy(ncolUJK, sizesUJK, sizesUJ + first) ;
196 IVcopy(nentUJK, indicesUJK, indicesUJ + offset) ;
197 if ( FRONTMTX_IS_REAL(frontmtx) ) {
198 DVcopy(nentUJK, entUJK, entUJ + offset) ;
199 } else if ( FRONTMTX_IS_COMPLEX(frontmtx) ) {
200 DVcopy(2*nentUJK, entUJK, entUJ + 2*offset) ;
201 }
202 count = 0 ;
203 offset += nentUJK ;
204 }
205 /*
206 -------------------------------------
207 initialize the row and column indices
208 -------------------------------------
209 */
210 if ( msglvl > 2 ) {
211 fprintf(msgFile, "\n setting row and column indices");
212 fflush(msgFile) ;
213 }
214 SubMtx_rowIndices(mtxUJK, &nrowUJK, &rowindUJK) ;
215 IVramp(nJ, rowindUJK, 0, 1) ;
216 SubMtx_columnIndices(mtxUJK, &ncolUJK, &colindUJK) ;
217 for ( ii = 0, jj = first ; ii < ncolUJK ; ii++, jj++ ) {
218 colindUJK[ii] = locmap[colindUJ[jj]] ;
219 }
220 /*
221 ----------------------------------
222 insert U_{J,K} into the hash table
223 ----------------------------------
224 */
225 if ( msglvl > 2 ) {
226 fprintf(msgFile,
227 "\n\n ## inserting U(%d,%d) ", J, K) ;
228 SubMtx_writeForHumanEye(mtxUJK, msgFile) ;
229 fflush(msgFile) ;
230 }
231 I2Ohash_insert(upperhash, J, K, (void *) mtxUJK) ;
232 /*
233 -----------------------------------
234 we jump to here if there were no
235 entries to be stored in the matrix.
236 -----------------------------------
237 */
238 no_entries :
239 /*
240 ----------------------------------------------------
241 reset first and K to new first location and front id
242 ----------------------------------------------------
243 */
244 first = jcol ;
245 if ( jcol < ncolUJ ) {
246 K = colmap[colindUJ[jcol]] ;
247 }
248 }
249 if ( jcol < ncolUJ && SUBMTX_IS_SPARSE_COLUMNS(mtxUJ) ) {
250 count += sizesUJ[jcol] ;
251 }
252 }
253 /*
254 --------------------------------------------
255 give U_{J,bnd{J}} back to the matrix manager
256 --------------------------------------------
257 */
258 SubMtxManager_releaseObject(manager, mtxUJ) ;
259 }
260 }
261 }
262 /*
263 ------------------------
264 free the working storage
265 ------------------------
266 */
267 IVfree(colmap) ;
268 IVfree(locmap) ;
269
270 return ; }
271
272 /*--------------------------------------------------------------------*/
273 /*
274 ----------------------------------------------------------------
275 purpose -- for each L_{bnd{J},J} matrix, remove from hash table,
276 split into their L_{K,J} submatrices and insert
277 into the hash table.
278
279 created -- 98may04, cca
280 ----------------------------------------------------------------
281 */
282 void
FrontMtx_splitLowerMatrices(FrontMtx * frontmtx,int msglvl,FILE * msgFile)283 FrontMtx_splitLowerMatrices (
284 FrontMtx *frontmtx,
285 int msglvl,
286 FILE *msgFile
287 ) {
288 SubMtx *mtxLJ, *mtxLJJ, *mtxLKJ ;
289 SubMtxManager *manager ;
290 double *entLJ, *entLKJ ;
291 int count, first, ii, inc1, inc2, irow, jj, J, K, nbytes,
292 ncolLJ, ncolLKJ, nentLJ, nentLKJ, neqns, nfront, nJ,
293 nrowJ, nrowLJ, nrowLKJ, offset, v ;
294 int *colindLJ, *colindLKJ, *rowmap, *indicesLJ, *indicesLKJ,
295 *locmap, *rowindJ, *rowindLJ, *rowindLKJ, *sizesLJ,
296 *sizesLKJ ;
297 I2Ohash *lowerhash ;
298 /*
299 ---------------
300 check the input
301 ---------------
302 */
303 if ( frontmtx == NULL || (msglvl > 0 && msgFile == NULL) ) {
304 fprintf(stderr,
305 "\n fatal error in FrontMtx_splitLowerMatrices(%p,%d,%p)"
306 "\n bad input\n", frontmtx, msglvl, msgFile) ;
307 exit(-1) ;
308 }
309 nfront = FrontMtx_nfront(frontmtx) ;
310 neqns = FrontMtx_neqns(frontmtx) ;
311 lowerhash = frontmtx->lowerhash ;
312 manager = frontmtx->manager ;
313 /*
314 --------------------------------
315 construct the row and local maps
316 --------------------------------
317 */
318 rowmap = IVinit(neqns, -1) ;
319 locmap = IVinit(neqns, -1) ;
320 for ( J = 0 ; J < nfront ; J++ ) {
321 if ( (nJ = FrontMtx_frontSize(frontmtx, J)) > 0 ) {
322 FrontMtx_rowIndices(frontmtx, J, &nrowJ, &rowindJ) ;
323 if ( nrowJ > 0 && rowindJ != NULL ) {
324 for ( ii = 0 ; ii < nJ ; ii++ ) {
325 v = rowindJ[ii] ;
326 rowmap[v] = J ;
327 locmap[v] = ii ;
328 }
329 }
330 }
331 }
332 if ( msglvl > 2 ) {
333 fprintf(msgFile, "\n\n rowmap[]") ;
334 IVfprintf(msgFile, neqns, rowmap) ;
335 fprintf(msgFile, "\n\n locmap[]") ;
336 IVfprintf(msgFile, neqns, locmap) ;
337 fflush(msgFile) ;
338 }
339 /*
340 ---------------------------------------------
341 move the L_{J,J} matrices into the hash table
342 ---------------------------------------------
343 */
344 for ( J = 0 ; J < nfront ; J++ ) {
345 if ( (mtxLJJ = FrontMtx_lowerMtx(frontmtx, J, J)) != NULL ) {
346 I2Ohash_insert(frontmtx->lowerhash, J, J, mtxLJJ) ;
347 }
348 }
349 /*
350 ------------------------------------------------------------
351 now split the L_{bnd{J},J} matrices into L_{K,J} matrices.
352 note: columns of L_{bnd{J},J} are assumed to be in ascending
353 order with respect to the column ordering of the matrix.
354 ------------------------------------------------------------
355 */
356 for ( J = 0 ; J < nfront ; J++ ) {
357 mtxLJ = FrontMtx_lowerMtx(frontmtx, nfront, J) ;
358 if ( msglvl > 2 ) {
359 fprintf(msgFile, "\n\n ### J = %d, mtxLJ = %p", J, mtxLJ) ;
360 fflush(msgFile) ;
361 }
362 if ( mtxLJ != NULL ) {
363 if ( msglvl > 2 ) {
364 SubMtx_writeForHumanEye(mtxLJ, msgFile) ;
365 fflush(msgFile) ;
366 }
367 SubMtx_columnIndices(mtxLJ, &ncolLJ, &colindLJ) ;
368 SubMtx_rowIndices(mtxLJ, &nrowLJ, &rowindLJ) ;
369 if ( msglvl > 2 ) {
370 fprintf(msgFile, "\n column indices for J") ;
371 IVfprintf(msgFile, ncolLJ, colindLJ) ;
372 fprintf(msgFile, "\n row indices for LJ") ;
373 IVfprintf(msgFile, nrowLJ, rowindLJ) ;
374 fflush(msgFile) ;
375 }
376 if ( (K = rowmap[rowindLJ[0]]) == rowmap[rowindLJ[nrowLJ-1]] ) {
377 if ( msglvl > 2 ) {
378 fprintf(msgFile, "\n front %d supports only %d", J, K) ;
379 fflush(msgFile) ;
380 }
381 /*
382 -------------------------------------------------
383 L_{bnd{J},J} is one submatrix, bnd{J} \subseteq K
384 set row and column indices and change column id
385 -------------------------------------------------
386 */
387 IVramp(ncolLJ, colindLJ, 0, 1) ;
388 for ( ii = 0 ; ii < nrowLJ ; ii++ ) {
389 rowindLJ[ii] = locmap[rowindLJ[ii]] ;
390 }
391 /*
392 mtxLJ->rowid = K ;
393 */
394 SubMtx_setFields(mtxLJ, mtxLJ->type, mtxLJ->mode, K, J,
395 mtxLJ->nrow, mtxLJ->ncol, mtxLJ->nent) ;
396 if ( msglvl > 2 ) {
397 fprintf(msgFile, "\n\n ## inserting L(%d,%d) ", K, J) ;
398 SubMtx_writeForHumanEye(mtxLJ, msgFile) ;
399 fflush(msgFile) ;
400 }
401 I2Ohash_insert(lowerhash, K, J, (void *) mtxLJ) ;
402 } else {
403 /*
404 -----------------------------------
405 split L_{bnd{J},J} into submatrices
406 -----------------------------------
407 */
408 nJ = FrontMtx_frontSize(frontmtx, J) ;
409 if ( SUBMTX_IS_DENSE_ROWS(mtxLJ) ) {
410 SubMtx_denseInfo(mtxLJ,
411 &nrowLJ, &ncolLJ, &inc1, &inc2, &entLJ) ;
412 } else if ( SUBMTX_IS_SPARSE_ROWS(mtxLJ) ) {
413 SubMtx_sparseRowsInfo(mtxLJ, &nrowLJ, &nentLJ,
414 &sizesLJ, &indicesLJ, &entLJ) ;
415 offset = 0 ;
416 count = sizesLJ[0] ;
417 }
418 first = 0 ;
419 K = rowmap[rowindLJ[0]] ;
420 for ( irow = 1 ; irow <= nrowLJ ; irow++ ) {
421 if ( msglvl > 2 ) {
422 fprintf(msgFile, "\n irow = %d", irow) ;
423 if ( irow < nrowLJ ) {
424 fprintf(msgFile, ", rowmap[%d] = %d",
425 rowindLJ[irow], rowmap[rowindLJ[irow]]);
426 }
427 fflush(msgFile) ;
428 }
429 if ( irow == nrowLJ || K != rowmap[rowindLJ[irow]] ) {
430 nrowLKJ = irow - first ;
431 if ( SUBMTX_IS_DENSE_ROWS(mtxLJ) ) {
432 nentLKJ = nJ*nrowLKJ ;
433 } else if ( SUBMTX_IS_SPARSE_ROWS(mtxLJ) ) {
434 if ( count == 0 ) {
435 goto no_entries ;
436 }
437 nentLKJ = count ;
438 }
439 nbytes = SubMtx_nbytesNeeded(mtxLJ->type, mtxLJ->mode,
440 nrowLKJ, nJ, nentLKJ) ;
441 mtxLKJ = SubMtxManager_newObjectOfSizeNbytes(manager,
442 nbytes) ;
443 SubMtx_init(mtxLKJ, mtxLJ->type, mtxLJ->mode, K, J,
444 nrowLKJ, nJ, nentLKJ) ;
445 if ( SUBMTX_IS_DENSE_ROWS(mtxLJ) ) {
446 SubMtx_denseInfo(mtxLKJ,
447 &nrowLKJ, &ncolLKJ, &inc1, &inc2, &entLKJ) ;
448 if ( FRONTMTX_IS_REAL(frontmtx) ) {
449 DVcopy(nentLKJ, entLKJ, entLJ + first*nJ) ;
450 } else if ( FRONTMTX_IS_COMPLEX(frontmtx) ) {
451 DVcopy(2*nentLKJ, entLKJ, entLJ + 2*first*nJ) ;
452 }
453 } else if ( SUBMTX_IS_SPARSE_ROWS(mtxLJ) ) {
454 SubMtx_sparseRowsInfo(mtxLKJ, &nrowLKJ, &nentLKJ,
455 &sizesLKJ, &indicesLKJ, &entLKJ) ;
456 IVcopy(nrowLKJ, sizesLKJ, sizesLJ + first) ;
457 IVcopy(nentLKJ, indicesLKJ, indicesLJ + offset) ;
458 if ( FRONTMTX_IS_REAL(frontmtx) ) {
459 DVcopy(nentLKJ, entLKJ, entLJ + offset) ;
460 } else if ( FRONTMTX_IS_COMPLEX(frontmtx) ) {
461 DVcopy(2*nentLKJ, entLKJ, entLJ + 2*offset) ;
462 }
463 count = 0 ;
464 offset += nentLKJ ;
465 }
466 /*
467 -------------------------------------
468 initialize the row and column indices
469 -------------------------------------
470 */
471 SubMtx_rowIndices(mtxLKJ, &nrowLKJ, &rowindLKJ) ;
472 for ( ii = 0, jj = first ; ii < nrowLKJ ; ii++, jj++ ) {
473 rowindLKJ[ii] = locmap[rowindLJ[jj]] ;
474 }
475 SubMtx_columnIndices(mtxLKJ, &ncolLKJ, &colindLKJ) ;
476 IVramp(ncolLKJ, colindLKJ, 0, 1) ;
477 /*
478 ----------------------------------
479 insert L_{K,J} into the hash table
480 ----------------------------------
481 */
482 if ( msglvl > 2 ) {
483 fprintf(msgFile,
484 "\n\n ## inserting L(%d,%d) ", K, J) ;
485 SubMtx_writeForHumanEye(mtxLKJ, msgFile) ;
486 fflush(msgFile) ;
487 }
488 I2Ohash_insert(lowerhash, K, J, (void *) mtxLKJ) ;
489 /*
490 -----------------------------------
491 we jump to here if there were no
492 entries to be stored in the matrix.
493 -----------------------------------
494 */
495 no_entries :
496 /*
497 ----------------------------------------------------
498 reset first and K to new first location and front id
499 ----------------------------------------------------
500 */
501 first = irow ;
502 if ( irow < nrowLJ ) {
503 K = rowmap[rowindLJ[irow]] ;
504 }
505 }
506 if ( irow < nrowLJ && SUBMTX_IS_SPARSE_ROWS(mtxLJ) ) {
507 count += sizesLJ[irow] ;
508 }
509 }
510 /*
511 --------------------------------------------
512 give L_{bnd{J},J} back to the matrix manager
513 --------------------------------------------
514 */
515 SubMtxManager_releaseObject(manager, mtxLJ) ;
516 }
517 }
518 }
519 /*
520 ------------------------
521 free the working storage
522 ------------------------
523 */
524 IVfree(rowmap) ;
525 IVfree(locmap) ;
526
527 return ; }
528
529 /*--------------------------------------------------------------------*/
530