1 /* storeFront.c */
2
3 #include "../FrontMtx.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 ----------------------------------------------------------------
8 purpose -- to store the factor indicies and entries from front J
9
10 pivotsizesIV -- used for symmetric or hermitian and pivoting
11 droptol -- used for drop tolerance factorization,
12 an entry is stored if its magnitude > droptol
13
14 created -- 98may04, cca
15 ----------------------------------------------------------------
16 */
17 void
FrontMtx_storeFront(FrontMtx * frontmtx,Chv * frontJ,IV * pivotsizesIV,double droptol,int msglvl,FILE * msgFile)18 FrontMtx_storeFront (
19 FrontMtx *frontmtx,
20 Chv *frontJ,
21 IV *pivotsizesIV,
22 double droptol,
23 int msglvl,
24 FILE *msgFile
25 ) {
26 SubMtx *mtx ;
27 double *entries ;
28 int inc1, inc2, ipivot, irow, jcol, J, m, nbytes,
29 ncol, nD, nent, nfront, nL, npivot, nrow, nU ;
30 int *colind, *colindJ, *firstlocs, *indices, *pivots,
31 *pivotsizes, *rowind, *rowindJ, *sizes ;
32 SubMtxManager *manager ;
33 /*
34 -------------------------------------
35 extract information from front matrix
36 -------------------------------------
37 */
38 nfront = frontmtx->nfront ;
39 manager = frontmtx->manager ;
40 if ( ( FRONTMTX_IS_SYMMETRIC(frontmtx)
41 || FRONTMTX_IS_HERMITIAN(frontmtx) )
42 && FRONTMTX_IS_PIVOTING(frontmtx) ) {
43 IV_sizeAndEntries(pivotsizesIV, &npivot, &pivotsizes) ;
44 } else {
45 npivot = 0, pivotsizes = NULL ;
46 }
47 /*
48 --------------------------------
49 extract information from chevron
50 --------------------------------
51 */
52 J = frontJ->id ;
53 Chv_dimensions(frontJ, &nD, &nL, &nU) ;
54 Chv_columnIndices(frontJ, &ncol, &colindJ) ;
55 Chv_rowIndices(frontJ, &nrow, &rowindJ) ;
56 if ( msglvl > 2 ) {
57 fprintf(msgFile, "\n inside FrontMtx_storeFront, J = %d", J) ;
58 fflush(msgFile) ;
59 }
60 if ( FRONTMTX_IS_PIVOTING(frontmtx) ) {
61 if ( frontmtx->lock != NULL ) {
62 if ( msglvl > 2 ) {
63 fprintf(msgFile, "\n locking lock") ;
64 fflush(msgFile) ;
65 }
66 Lock_lock(frontmtx->lock) ;
67 frontmtx->nlocks++ ;
68 }
69 /*
70 --------------------
71 store the front size
72 --------------------
73 */
74 FrontMtx_setFrontSize(frontmtx, J, nD) ;
75 /*
76 ------------------------
77 store the column indices
78 ------------------------
79 */
80 if ( msglvl > 2 ) {
81 fprintf(msgFile, "\n setting column indices, J = %d", J) ;
82 fflush(msgFile) ;
83 }
84 IVL_setList(frontmtx->coladjIVL, J, ncol, colindJ) ;
85 if ( FRONTMTX_IS_NONSYMMETRIC(frontmtx) ) {
86 /*
87 ---------------------
88 store the row indices
89 ---------------------
90 */
91 if ( msglvl > 2 ) {
92 fprintf(msgFile, "\n setting row indices, J = %d", J) ;
93 fflush(msgFile) ;
94 }
95 IVL_setList(frontmtx->rowadjIVL, J, nrow, rowindJ) ;
96 }
97 if ( frontmtx->lock != NULL ) {
98 if ( msglvl > 2 ) {
99 fprintf(msgFile, "\n unlocking lock") ;
100 fflush(msgFile) ;
101 }
102 Lock_unlock(frontmtx->lock) ;
103 }
104 }
105 if ( nD == 0 ) {
106 return ;
107 }
108 /*
109 ------------------------
110 store the D_{J,J} matrix
111 ------------------------
112 */
113 if ( pivotsizes != NULL ) {
114 /*
115 ---------------------------------------------------
116 symmetric and pivoting, store block diagonal matrix
117 ---------------------------------------------------
118 */
119 nent = Chv_countEntries(frontJ, npivot, pivotsizes, CHV_DIAGONAL) ;
120 nbytes = SubMtx_nbytesNeeded(frontJ->type, SUBMTX_BLOCK_DIAGONAL_SYM,
121 nD, nD, nent) ;
122 mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
123 if ( FRONTMTX_IS_SYMMETRIC(frontmtx) ) {
124 SubMtx_init(mtx, frontJ->type, SUBMTX_BLOCK_DIAGONAL_SYM,
125 J, J, nD, nD, nent) ;
126 } else if ( FRONTMTX_IS_HERMITIAN(frontmtx) ) {
127 SubMtx_init(mtx, frontJ->type, SUBMTX_BLOCK_DIAGONAL_HERM,
128 J, J, nD, nD, nent) ;
129 }
130 SubMtx_blockDiagonalInfo(mtx, &nrow, &nent, &pivots, &entries) ;
131 IVzero(nD, pivots) ;
132 IVcopy(npivot, pivots, pivotsizes) ;
133 Chv_copyEntriesToVector(frontJ, npivot, pivotsizes, nent, entries,
134 CHV_DIAGONAL, CHV_BY_ROWS) ;
135 if ( msglvl > 2 ) {
136 fprintf(msgFile, "\n\n block diagonal matrix") ;
137 SubMtx_writeForHumanEye(mtx, msgFile) ;
138 fflush(msgFile) ;
139 }
140 frontmtx->nentD += nent ;
141 } else {
142 /*
143 ---------------
144 diagonal matrix
145 ---------------
146 */
147 mtx = frontmtx->p_mtxDJJ[J] ;
148 if ( mtx == NULL ) {
149 nbytes = SubMtx_nbytesNeeded(frontJ->type, SUBMTX_DIAGONAL,
150 nD, nD, nD) ;
151 mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
152 SubMtx_init(mtx, frontJ->type, SUBMTX_DIAGONAL, J, J, nD, nD, nD);
153 }
154 SubMtx_diagonalInfo(mtx, &nent, &entries) ;
155 Chv_copyEntriesToVector(frontJ, npivot, pivotsizes, nent, entries,
156 CHV_DIAGONAL, CHV_BY_ROWS) ;
157 frontmtx->nentD += nD ;
158 if ( msglvl > 2 ) {
159 fprintf(msgFile, "\n\n entries in SubMtx") ;
160 DVfprintf(msgFile, nent, entries) ;
161 fprintf(msgFile, "\n\n block diagonal matrix") ;
162 SubMtx_writeForHumanEye(mtx, msgFile) ;
163 fflush(msgFile) ;
164 }
165 }
166 frontmtx->p_mtxDJJ[J] = mtx ;
167 SubMtx_columnIndices(mtx, &ncol, &colind) ;
168 IVcopy(ncol, colind, colindJ) ;
169 SubMtx_rowIndices(mtx, &nrow, &rowind) ;
170 IVcopy(nrow, rowind, rowindJ) ;
171 /*
172 ------------------------------
173 store the upper U_{J,J} matrix
174 ------------------------------
175 */
176 if ( FRONTMTX_IS_DENSE_FRONTS(frontmtx) ) {
177 /*
178 -----------
179 dense front
180 -----------
181 */
182 nent = Chv_countEntries(frontJ, npivot, pivotsizes,
183 CHV_STRICT_UPPER_11) ;
184 if ( nent > 0 ) {
185 mtx = frontmtx->p_mtxUJJ[J] ;
186 if ( mtx == NULL ) {
187 nbytes = SubMtx_nbytesNeeded(frontJ->type,
188 SUBMTX_DENSE_SUBCOLUMNS, nD, nD,nent);
189 mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
190 SubMtx_init(mtx, frontJ->type, SUBMTX_DENSE_SUBCOLUMNS,
191 J, J, nD, nD, nent) ;
192 }
193 SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
194 &firstlocs, &sizes, &entries) ;
195 if ( pivotsizes != NULL ) {
196 for ( jcol = ipivot = 0 ; jcol < nD ; ipivot++ ) {
197 m = pivotsizes[ipivot] ;
198 if ( m == 1 ) {
199 firstlocs[jcol] = 0 ;
200 sizes[jcol] = jcol ;
201 jcol++ ;
202 } else if ( m == 2 ) {
203 firstlocs[jcol] = firstlocs[jcol+1] = 0 ;
204 sizes[jcol] = sizes[jcol+1] = jcol ;
205 jcol += 2 ;
206 }
207 }
208 } else {
209 for ( jcol = 0 ; jcol < nD ; jcol++ ) {
210 firstlocs[jcol] = 0 ;
211 sizes[jcol] = jcol ;
212 }
213 }
214 Chv_copyEntriesToVector(frontJ, npivot, pivotsizes, nent, entries,
215 CHV_STRICT_UPPER_11, CHV_BY_COLUMNS) ;
216 frontmtx->p_mtxUJJ[J] = mtx ;
217 SubMtx_columnIndices(mtx, &ncol, &colind) ;
218 IVcopy(ncol, colind, colindJ) ;
219 SubMtx_rowIndices(mtx, &nrow, &rowind) ;
220 IVcopy(nrow, rowind, rowindJ) ;
221 if ( FRONTMTX_IS_PIVOTING(frontmtx) ) {
222 frontmtx->nentU += nent ;
223 }
224 if ( msglvl > 2 ) {
225 fprintf(msgFile, "\n\n UJJ matrix") ;
226 SubMtx_writeForHumanEye(mtx, msgFile) ;
227 fflush(msgFile) ;
228 }
229 }
230 } else {
231 /*
232 ------------
233 sparse front
234 ------------
235 */
236 nent = Chv_countBigEntries(frontJ, npivot, pivotsizes,
237 CHV_STRICT_UPPER_11, droptol) ;
238 if ( nent > 0 ) {
239 nbytes = SubMtx_nbytesNeeded(frontJ->type, SUBMTX_SPARSE_COLUMNS,
240 nD, nD, nent) ;
241 if ( msglvl > 2 ) {
242 fprintf(msgFile,
243 "\n U_{%d,%d}, nD %d, nent %d, nbytes %d",
244 J, J, nD, nent, nbytes) ;
245 fflush(msgFile) ;
246 }
247 mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
248 SubMtx_init(mtx, frontJ->type, SUBMTX_SPARSE_COLUMNS,
249 J, J, nD, nD, nent) ;
250 SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
251 &sizes, &indices, &entries) ;
252 Chv_copyBigEntriesToVector(frontJ, npivot, pivotsizes, sizes,
253 indices, entries, CHV_STRICT_UPPER_11,
254 CHV_BY_COLUMNS, droptol) ;
255 frontmtx->p_mtxUJJ[J] = mtx ;
256 SubMtx_columnIndices(mtx, &ncol, &colind) ;
257 IVcopy(ncol, colind, colindJ) ;
258 SubMtx_rowIndices(mtx, &nrow, &rowind) ;
259 IVcopy(nrow, rowind, rowindJ) ;
260 frontmtx->nentU += nent ;
261 if ( msglvl > 2 ) {
262 fprintf(msgFile, "\n\n UJJ matrix") ;
263 SubMtx_writeForHumanEye(mtx, msgFile) ;
264 fflush(msgFile) ;
265 }
266 }
267 }
268 /*
269 -----------------------------------
270 store the upper U_{J,bnd{J}} matrix
271 -----------------------------------
272 */
273 if ( FRONTMTX_IS_DENSE_FRONTS(frontmtx) ) {
274 /*
275 -----------
276 dense front
277 -----------
278 */
279 nent = Chv_countEntries(frontJ, npivot, pivotsizes, CHV_UPPER_12) ;
280 if ( nent > 0 ) {
281 mtx = frontmtx->p_mtxUJN[J] ;
282 if ( mtx == NULL ) {
283 nbytes = SubMtx_nbytesNeeded(frontJ->type,
284 SUBMTX_DENSE_COLUMNS, nD, nU, nent) ;
285 mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
286 SubMtx_init(mtx, frontJ->type, SUBMTX_DENSE_COLUMNS,
287 J, nfront, nD, nU, nent) ;
288 }
289 SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
290 Chv_copyEntriesToVector(frontJ, npivot, pivotsizes, nent, entries,
291 CHV_UPPER_12, CHV_BY_COLUMNS) ;
292 frontmtx->p_mtxUJN[J] = mtx ;
293 SubMtx_columnIndices(mtx, &ncol, &colind) ;
294 IVcopy(ncol, colind, colindJ + nD) ;
295 SubMtx_rowIndices(mtx, &nrow, &rowind) ;
296 IVcopy(nrow, rowind, rowindJ) ;
297 if ( FRONTMTX_IS_PIVOTING(frontmtx) ) {
298 frontmtx->nentU += nent ;
299 }
300 if ( msglvl > 2 ) {
301 fprintf(msgFile, "\n\n UJJN matrix") ;
302 SubMtx_writeForHumanEye(mtx, msgFile) ;
303 fflush(msgFile) ;
304 }
305 }
306 } else {
307 /*
308 ------------
309 sparse front
310 ------------
311 */
312 nent = Chv_countBigEntries(frontJ, npivot, pivotsizes,
313 CHV_UPPER_12, droptol) ;
314 if ( nent > 0 ) {
315 nbytes = SubMtx_nbytesNeeded(frontJ->type, SUBMTX_SPARSE_COLUMNS,
316 nD, nU, nent) ;
317 if ( msglvl > 2 ) {
318 fprintf(msgFile,
319 "\n U_{%d,%d}, nD %d, nU %d, nent %d, nbytes %d",
320 J, nfront, nD, nU, nent, nbytes) ;
321 fflush(msgFile) ;
322 }
323 mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
324 SubMtx_init(mtx, frontJ->type, SUBMTX_SPARSE_COLUMNS,
325 J, nfront, nD, nU, nent) ;
326 SubMtx_sparseColumnsInfo(mtx,
327 &nrow, &nent, &sizes, &indices, &entries) ;
328 Chv_copyBigEntriesToVector(frontJ, npivot, pivotsizes, sizes,
329 indices, entries, CHV_UPPER_12,
330 CHV_BY_COLUMNS, droptol) ;
331 frontmtx->p_mtxUJN[J] = mtx ;
332 SubMtx_columnIndices(mtx, &ncol, &colind) ;
333 IVcopy(ncol, colind, colindJ + nD) ;
334 SubMtx_rowIndices(mtx, &nrow, &rowind) ;
335 IVcopy(nrow, rowind, rowindJ) ;
336 frontmtx->nentU += nent ;
337 if ( msglvl > 2 ) {
338 fprintf(msgFile, "\n\n UJJN matrix") ;
339 SubMtx_writeForHumanEye(mtx, msgFile) ;
340 fflush(msgFile) ;
341 }
342 }
343 }
344 if ( FRONTMTX_IS_NONSYMMETRIC(frontmtx) ) {
345 /*
346 ------------------------------
347 store the lower L_{J,J} matrix
348 ------------------------------
349 */
350 if ( FRONTMTX_IS_DENSE_FRONTS(frontmtx) ) {
351 /*
352 -----------
353 dense front
354 -----------
355 */
356 nent = Chv_countEntries(frontJ, npivot, pivotsizes,
357 CHV_STRICT_LOWER_11) ;
358 if ( nent > 0 ) {
359 mtx = frontmtx->p_mtxLJJ[J] ;
360 if ( mtx == NULL ) {
361 nbytes = SubMtx_nbytesNeeded(frontJ->type,
362 SUBMTX_DENSE_SUBROWS, nD, nD,nent);
363 mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
364 SubMtx_init(mtx, frontJ->type, SUBMTX_DENSE_SUBROWS,
365 J, J, nD, nD, nent) ;
366 }
367 SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
368 &firstlocs, &sizes, &entries) ;
369 for ( irow = 0 ; irow < nD ; irow++ ) {
370 firstlocs[irow] = 0 ;
371 sizes[irow] = irow ;
372 }
373 Chv_copyEntriesToVector(frontJ, npivot, pivotsizes, nent,
374 entries, CHV_STRICT_LOWER_11, CHV_BY_ROWS) ;
375 frontmtx->p_mtxLJJ[J] = mtx ;
376 SubMtx_columnIndices(mtx, &ncol, &colind) ;
377 IVcopy(ncol, colind, colindJ) ;
378 SubMtx_rowIndices(mtx, &nrow, &rowind) ;
379 IVcopy(nrow, rowind, rowindJ) ;
380 if ( FRONTMTX_IS_PIVOTING(frontmtx) ) {
381 frontmtx->nentL += nent ;
382 }
383 if ( msglvl > 2 ) {
384 fprintf(msgFile, "\n\n LJJ matrix") ;
385 SubMtx_writeForHumanEye(mtx, msgFile) ;
386 fflush(msgFile) ;
387 }
388 }
389 } else {
390 /*
391 ------------
392 sparse front
393 ------------
394 */
395 nent = Chv_countBigEntries(frontJ, npivot, pivotsizes,
396 CHV_STRICT_LOWER_11, droptol);
397 if ( nent > 0 ) {
398 nbytes = SubMtx_nbytesNeeded(frontJ->type, SUBMTX_SPARSE_ROWS,
399 nD, nD, nent) ;
400 if ( msglvl > 2 ) {
401 fprintf(msgFile,
402 "\n L_{%d,%d}, nD %d, nent %d, nbytes %d",
403 J, J, nD, nent, nbytes) ;
404 fflush(msgFile) ;
405 }
406 mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
407 SubMtx_init(mtx, frontJ->type, SUBMTX_SPARSE_ROWS,
408 J, J, nD, nD, nent) ;
409 SubMtx_sparseRowsInfo(mtx,
410 &nrow, &nent, &sizes, &indices, &entries);
411 Chv_copyBigEntriesToVector(frontJ, npivot, pivotsizes, sizes,
412 indices, entries, CHV_STRICT_LOWER_11,
413 CHV_BY_ROWS, droptol) ;
414 frontmtx->p_mtxLJJ[J] = mtx ;
415 SubMtx_columnIndices(mtx, &ncol, &colind) ;
416 IVcopy(ncol, colind, colindJ) ;
417 SubMtx_rowIndices(mtx, &nrow, &rowind) ;
418 IVcopy(nrow, rowind, rowindJ) ;
419 frontmtx->nentL += nent ;
420 if ( msglvl > 2 ) {
421 fprintf(msgFile, "\n\n LJJ matrix") ;
422 SubMtx_writeForHumanEye(mtx, msgFile) ;
423 fflush(msgFile) ;
424 }
425 }
426 }
427 /*
428 -----------------------------------
429 store the lower L_{bnd{J},J} matrix
430 -----------------------------------
431 */
432 if ( FRONTMTX_IS_DENSE_FRONTS(frontmtx) ) {
433 /*
434 -----------
435 dense front
436 -----------
437 */
438 nent = Chv_countEntries(frontJ, npivot, pivotsizes, CHV_LOWER_21);
439 if ( nent > 0 ) {
440 mtx = frontmtx->p_mtxLNJ[J] ;
441 if ( mtx == NULL ) {
442 nbytes = SubMtx_nbytesNeeded(frontJ->type,
443 SUBMTX_DENSE_ROWS, nL, nD, nent) ;
444 mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
445 SubMtx_init(mtx, frontJ->type, SUBMTX_DENSE_ROWS,
446 nfront, J, nL, nD, nent) ;
447 }
448 SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
449 Chv_copyEntriesToVector(frontJ, npivot, pivotsizes, nent,
450 entries, CHV_LOWER_21, CHV_BY_ROWS) ;
451 frontmtx->p_mtxLNJ[J] = mtx ;
452 SubMtx_columnIndices(mtx, &ncol, &colind) ;
453 IVcopy(ncol, colind, colindJ) ;
454 SubMtx_rowIndices(mtx, &nrow, &rowind) ;
455 IVcopy(nrow, rowind, rowindJ + nD) ;
456 if ( FRONTMTX_IS_PIVOTING(frontmtx) ) {
457 frontmtx->nentL += nent ;
458 }
459 if ( msglvl > 2 ) {
460 fprintf(msgFile, "\n\n LNJ matrix") ;
461 SubMtx_writeForHumanEye(mtx, msgFile) ;
462 fflush(msgFile) ;
463 }
464 }
465 } else {
466 /*
467 ------------
468 sparse front
469 ------------
470 */
471 nent = Chv_countBigEntries(frontJ, npivot, pivotsizes,
472 CHV_LOWER_21, droptol) ;
473 if ( nent > 0 ) {
474 nbytes = SubMtx_nbytesNeeded(frontJ->type, SUBMTX_SPARSE_ROWS,
475 nL, nD, nent) ;
476 if ( msglvl > 2 ) {
477 fprintf(msgFile,
478 "\n L_{%d,%d}, nL %d, nD %d, nent %d, nbytes %d",
479 nfront, J, nL, nD, nent, nbytes) ;
480 fflush(msgFile) ;
481 }
482 mtx = SubMtxManager_newObjectOfSizeNbytes(manager, nbytes) ;
483 SubMtx_init(mtx, frontJ->type, SUBMTX_SPARSE_ROWS,
484 nfront, J, nL, nD, nent) ;
485 SubMtx_sparseRowsInfo(mtx,
486 &nrow, &nent, &sizes, &indices, &entries) ;
487 Chv_copyBigEntriesToVector(frontJ, npivot, pivotsizes, sizes,
488 indices, entries, CHV_LOWER_21,
489 CHV_BY_ROWS, droptol) ;
490 frontmtx->p_mtxLNJ[J] = mtx ;
491 SubMtx_columnIndices(mtx, &ncol, &colind) ;
492 IVcopy(ncol, colind, colindJ) ;
493 SubMtx_rowIndices(mtx, &nrow, &rowind) ;
494 IVcopy(nrow, rowind, rowindJ + nD) ;
495 frontmtx->nentL += nent ;
496 if ( msglvl > 2 ) {
497 fprintf(msgFile, "\n\n LNJ matrix") ;
498 SubMtx_writeForHumanEye(mtx, msgFile) ;
499 fflush(msgFile) ;
500 }
501 }
502 }
503 }
504 return ; }
505
506 /*--------------------------------------------------------------------*/
507