1 /* util.c */
2
3 #include "../FrontMtx.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 -----------------------------------------
8 purpose -- produce a map from each column
9 to the front that contains it
10
11 created -- 98may04, cca
12 -----------------------------------------
13 */
14 IV *
FrontMtx_colmapIV(FrontMtx * frontmtx)15 FrontMtx_colmapIV (
16 FrontMtx *frontmtx
17 ) {
18 int ii, J, ncolJ, neqns, nfront, nJ ;
19 int *colindJ, *colmap ;
20 IV *colmapIV ;
21 /*
22 -----------------------------------------
23 get the map from columns to owning fronts
24 -----------------------------------------
25 */
26 neqns = FrontMtx_neqns(frontmtx) ;
27 nfront = FrontMtx_nfront(frontmtx) ;
28 colmapIV = IV_new() ;
29 IV_init(colmapIV, neqns, NULL) ;
30 colmap = IV_entries(colmapIV) ;
31 IVfill(neqns, colmap, -1) ;
32 for ( J = 0 ; J < nfront ; J++ ) {
33 if ( (nJ = FrontMtx_frontSize(frontmtx, J)) > 0 ) {
34 FrontMtx_columnIndices(frontmtx, J, &ncolJ, &colindJ) ;
35 if ( ncolJ > 0 && colindJ != NULL ) {
36 for ( ii = 0 ; ii < nJ ; ii++ ) {
37 colmap[colindJ[ii]] = J ;
38 }
39 }
40 }
41 }
42 return(colmapIV) ; }
43
44 /*--------------------------------------------------------------------*/
45 /*
46 --------------------------------------------------------------------
47 purpose -- produce a map from each row to the front that contains it
48
49 created -- 98may04, cca
50 --------------------------------------------------------------------
51 */
52 IV *
FrontMtx_rowmapIV(FrontMtx * frontmtx)53 FrontMtx_rowmapIV (
54 FrontMtx *frontmtx
55 ) {
56 int ii, J, nrowJ, neqns, nfront, nJ ;
57 int *rowindJ, *rowmap ;
58 IV *rowmapIV ;
59 /*
60 --------------------------------------
61 get the map from rows to owning fronts
62 --------------------------------------
63 */
64 neqns = FrontMtx_neqns(frontmtx) ;
65 nfront = FrontMtx_nfront(frontmtx) ;
66 rowmapIV = IV_new() ;
67 IV_init(rowmapIV, neqns, NULL) ;
68 rowmap = IV_entries(rowmapIV) ;
69 IVfill(neqns, rowmap, -1) ;
70 for ( J = 0 ; J < nfront ; J++ ) {
71 if ( (nJ = FrontMtx_frontSize(frontmtx, J)) > 0 ) {
72 FrontMtx_rowIndices(frontmtx, J, &nrowJ, &rowindJ) ;
73 if ( nrowJ > 0 && rowindJ != NULL ) {
74 for ( ii = 0 ; ii < nJ ; ii++ ) {
75 rowmap[rowindJ[ii]] = J ;
76 }
77 }
78 }
79 }
80 return(rowmapIV) ; }
81
82 /*--------------------------------------------------------------------*/
83 /*
84 -------------------------------------------------------------
85 compute the inertia of a symmetric matrix
86
87 fill *pnnegative with the number of negative eigenvalues of A
88 fill *pnzero with the number of zero eigenvalues of A
89 fill *pnpositive with the number of positive eigenvalues of A
90
91 created -- 98may04, cca
92 -------------------------------------------------------------
93 */
94 void
FrontMtx_inertia(FrontMtx * frontmtx,int * pnnegative,int * pnzero,int * pnpositive)95 FrontMtx_inertia (
96 FrontMtx *frontmtx,
97 int *pnnegative,
98 int *pnzero,
99 int *pnpositive
100 ) {
101 SubMtx *mtx ;
102 double arm, areal, bimag, breal, creal, mid, val ;
103 double *entries ;
104 int ii, ipivot, irow, J, nent, nfront, nJ,
105 nnegative, npositive, nzero ;
106 int *pivotsizes ;
107 /*
108 ---------------
109 check the input
110 ---------------
111 */
112 if ( frontmtx == NULL
113 || pnnegative == NULL || pnzero == NULL || pnpositive == NULL ) {
114 fprintf(stderr, "\n fatal error in FrontMtx_inertia(%p,%p,%p,%p)"
115 "\n bad input\n",
116 frontmtx, pnnegative, pnzero, pnpositive) ;
117 fflush(stdout) ;
118 }
119 if ( FRONTMTX_IS_REAL(frontmtx) && ! FRONTMTX_IS_SYMMETRIC(frontmtx) ) {
120 fprintf(stderr, "\n fatal error in FrontMtx_inertia(%p,%p,%p,%p)"
121 "\n matrix is real and not symmetric \n",
122 frontmtx, pnnegative, pnzero, pnpositive) ;
123 fflush(stdout) ;
124 } else if ( FRONTMTX_IS_COMPLEX(frontmtx)
125 && ! FRONTMTX_IS_HERMITIAN(frontmtx) ) {
126 fprintf(stderr, "\n fatal error in FrontMtx_inertia(%p,%p,%p,%p)"
127 "\n matrix is complex and not hermitian \n",
128 frontmtx, pnnegative, pnzero, pnpositive) ;
129 fflush(stdout) ;
130 }
131 nfront = frontmtx->nfront ;
132 nnegative = nzero = npositive = 0 ;
133 for ( J = 0 ; J < nfront ; J++ ) {
134 mtx = FrontMtx_diagMtx(frontmtx, J) ;
135 if ( mtx != NULL ) {
136 if ( ! FRONTMTX_IS_PIVOTING(frontmtx) ) {
137 /*
138 -----------
139 no pivoting
140 -----------
141 */
142 SubMtx_diagonalInfo(mtx, &nJ, &entries) ;
143 if ( FRONTMTX_IS_REAL(frontmtx) ) {
144 for ( ii = 0 ; ii < nJ ; ii++ ) {
145 if ( entries[ii] < 0.0 ) {
146 nnegative++ ;
147 } else if ( entries[ii] > 0.0 ) {
148 npositive++ ;
149 } else {
150 nzero++ ;
151 }
152 }
153 } else if ( FRONTMTX_IS_COMPLEX(frontmtx) ) {
154 for ( ii = 0 ; ii < nJ ; ii++ ) {
155 if ( entries[2*ii] < 0.0 ) {
156 nnegative++ ;
157 } else if ( entries[2*ii] > 0.0 ) {
158 npositive++ ;
159 } else {
160 nzero++ ;
161 }
162 }
163 }
164 } else {
165 /*
166 --------
167 pivoting
168 --------
169 */
170 SubMtx_blockDiagonalInfo(mtx, &nJ, &nent,
171 &pivotsizes, &entries) ;
172 if ( FRONTMTX_IS_REAL(frontmtx) ) {
173 for ( irow = ipivot = ii = 0 ; irow < nJ ; ipivot++ ) {
174 if ( pivotsizes[ipivot] == 1 ) {
175 val = entries[ii] ;
176 if ( val < 0.0 ) {
177 nnegative++ ;
178 } else if ( val > 0.0 ) {
179 npositive++ ;
180 } else {
181 nzero++ ;
182 }
183 irow++ ; ii++ ;
184 } else {
185 areal = entries[ii] ;
186 breal = entries[ii+1] ;
187 creal = entries[ii+2] ;
188 mid = 0.5*(areal + creal) ;
189 arm = sqrt(0.25*(areal - creal)*(areal - creal)
190 + breal*breal) ;
191 val = mid + arm ;
192 if ( val < 0.0 ) {
193 nnegative++ ;
194 } else if ( val > 0.0 ) {
195 npositive++ ;
196 } else {
197 nzero++ ;
198 }
199 val = mid - arm ;
200 if ( val < 0.0 ) {
201 nnegative++ ;
202 } else if ( val > 0.0 ) {
203 npositive++ ;
204 } else {
205 nzero++ ;
206 }
207 irow += 2 ; ii += 3 ;
208 }
209 }
210 } else if ( FRONTMTX_IS_COMPLEX(frontmtx) ) {
211 for ( irow = ipivot = ii = 0 ; irow < nJ ; ipivot++ ) {
212 if ( pivotsizes[ipivot] == 1 ) {
213 val = entries[2*ii] ;
214 if ( val < 0.0 ) {
215 nnegative++ ;
216 } else if ( val > 0.0 ) {
217 npositive++ ;
218 } else {
219 nzero++ ;
220 }
221 irow++ ; ii++ ;
222 } else {
223 areal = entries[2*ii] ;
224 breal = entries[2*ii+2] ;
225 bimag = entries[2*ii+3] ;
226 creal = entries[2*ii+4] ;
227 mid = 0.5*(areal + creal) ;
228 arm = sqrt(0.25*(areal - creal)*(areal - creal)
229 + breal*breal + bimag*bimag) ;
230 val = mid + arm ;
231 if ( val < 0.0 ) {
232 nnegative++ ;
233 } else if ( val > 0.0 ) {
234 npositive++ ;
235 } else {
236 nzero++ ;
237 }
238 val = mid - arm ;
239 if ( val < 0.0 ) {
240 nnegative++ ;
241 } else if ( val > 0.0 ) {
242 npositive++ ;
243 } else {
244 nzero++ ;
245 }
246 irow += 2 ; ii += 3 ;
247 }
248 }
249 }
250 }
251 }
252 }
253 *pnnegative = nnegative ;
254 *pnzero = nzero ;
255 *pnpositive = npositive ;
256
257 return ; }
258
259 /*--------------------------------------------------------------------*/
260 /*
261 -------------------------------------------------------
262 purpose -- create and return an IV object that contains
263 all the row ids owned by process myid.
264
265 created -- 98jun13, cca
266 -------------------------------------------------------
267 */
268 IV *
FrontMtx_ownedRowsIV(FrontMtx * frontmtx,int myid,IV * ownersIV,int msglvl,FILE * msgFile)269 FrontMtx_ownedRowsIV (
270 FrontMtx *frontmtx,
271 int myid,
272 IV *ownersIV,
273 int msglvl,
274 FILE *msgFile
275 ) {
276 int J, neqns, nfront, nJ, nowned, nrowJ, offset ;
277 int *ownedRows, *owners, *rowindJ ;
278 IV *ownedRowsIV ;
279 /*
280 ---------------
281 check the input
282 ---------------
283 */
284 if ( frontmtx == NULL ) {
285 fprintf(stderr, "\n fatal error in FrontMtx_ownedRowsIV(%p,%d,%p)"
286 "\n bad input\n", frontmtx, myid, ownersIV) ;
287 exit(-1) ;
288 }
289 nfront = frontmtx->nfront ;
290 neqns = frontmtx->neqns ;
291 ownedRowsIV = IV_new() ;
292 if ( ownersIV == NULL ) {
293 IV_init(ownedRowsIV, neqns, NULL) ;
294 IV_ramp(ownedRowsIV, 0, 1) ;
295 } else {
296 owners = IV_entries(ownersIV) ;
297 for ( J = 0, nowned = 0 ; J < nfront ; J++ ) {
298 if ( owners[J] == myid ) {
299 nJ = FrontMtx_frontSize(frontmtx, J) ;
300 nowned += nJ ;
301 }
302 }
303 if ( nowned > 0 ) {
304 IV_init(ownedRowsIV, nowned, NULL) ;
305 ownedRows = IV_entries(ownedRowsIV) ;
306 for ( J = 0, offset = 0 ; J < nfront ; J++ ) {
307 if ( owners[J] == myid ) {
308 nJ = FrontMtx_frontSize(frontmtx, J) ;
309 if ( nJ > 0 ) {
310 FrontMtx_rowIndices(frontmtx, J, &nrowJ, &rowindJ) ;
311 IVcopy(nJ, ownedRows + offset, rowindJ) ;
312 offset += nJ ;
313 }
314 }
315 }
316 }
317 }
318 return(ownedRowsIV) ; }
319
320 /*--------------------------------------------------------------------*/
321 /*
322 -------------------------------------------------------
323 purpose -- create and return an IV object that contains
324 all the column ids owned by process myid.
325
326 created -- 98jun13, cca
327 -------------------------------------------------------
328 */
329 IV *
FrontMtx_ownedColumnsIV(FrontMtx * frontmtx,int myid,IV * ownersIV,int msglvl,FILE * msgFile)330 FrontMtx_ownedColumnsIV (
331 FrontMtx *frontmtx,
332 int myid,
333 IV *ownersIV,
334 int msglvl,
335 FILE *msgFile
336 ) {
337 int J, neqns, nfront, nJ, nowned, ncolJ, offset ;
338 int *ownedColumns, *owners, *colindJ ;
339 IV *ownedColumnsIV ;
340 /*
341 ---------------
342 check the input
343 ---------------
344 */
345 if ( frontmtx == NULL ) {
346 fprintf(stderr, "\n fatal error in FrontMtx_ownedColumnsIV(%p,%d,%p)"
347 "\n bad input\n", frontmtx, myid, ownersIV) ;
348 exit(-1) ;
349 }
350 nfront = frontmtx->nfront ;
351 neqns = frontmtx->neqns ;
352 ownedColumnsIV = IV_new() ;
353 if ( ownersIV == NULL ) {
354 IV_init(ownedColumnsIV, neqns, NULL) ;
355 IV_ramp(ownedColumnsIV, 0, 1) ;
356 } else {
357 owners = IV_entries(ownersIV) ;
358 for ( J = 0, nowned = 0 ; J < nfront ; J++ ) {
359 if ( owners[J] == myid ) {
360 nJ = FrontMtx_frontSize(frontmtx, J) ;
361 nowned += nJ ;
362 }
363 }
364 if ( nowned > 0 ) {
365 IV_init(ownedColumnsIV, nowned, NULL) ;
366 ownedColumns = IV_entries(ownedColumnsIV) ;
367 for ( J = 0, offset = 0 ; J < nfront ; J++ ) {
368 if ( owners[J] == myid ) {
369 nJ = FrontMtx_frontSize(frontmtx, J) ;
370 if ( nJ > 0 ) {
371 FrontMtx_columnIndices(frontmtx, J, &ncolJ, &colindJ) ;
372 IVcopy(nJ, ownedColumns + offset, colindJ) ;
373 offset += nJ ;
374 }
375 }
376 }
377 }
378 }
379 return(ownedColumnsIV) ; }
380
381 /*--------------------------------------------------------------------*/
382 /*
383 ----------------------------------------------------------------
384 purpose -- to create and return an IVL object that holds the
385 submatrix nonzero pattern for the upper triangular factor.
386
387 NOTE: this method supercedes calling IVL_mapEntries() on
388 the column adjacency structure. that gave problems when
389 pivoting forced some fronts to have no eliminated columns.
390 in some cases, solve aggregates were expected when none
391 were forthcoming.
392
393 created -- 98aug20, cca
394 ----------------------------------------------------------------
395 */
396 IVL *
FrontMtx_makeUpperBlockIVL(FrontMtx * frontmtx,IV * colmapIV)397 FrontMtx_makeUpperBlockIVL (
398 FrontMtx *frontmtx,
399 IV *colmapIV
400 ) {
401 int count, ii, J, K, ncol, nfront, nJ ;
402 int *colmap, *colind, *list, *mark ;
403 IVL *upperblockIVL ;
404 /*
405 ---------------
406 check the input
407 ---------------
408 */
409 if ( frontmtx == NULL || colmapIV == NULL ) {
410 fprintf(stderr, "\n fatal error in FrontMtx_makeUpperBlockIVL()"
411 "\n bad input\n") ;
412 exit(-1) ;
413 }
414 nfront = FrontMtx_nfront(frontmtx) ;
415 colmap = IV_entries(colmapIV) ;
416 /*
417 -----------------------------
418 set up the working storage
419 and initialize the IVL object
420 -----------------------------
421 */
422 mark = IVinit(nfront, -1) ;
423 list = IVinit(nfront, -1) ;
424 upperblockIVL = IVL_new() ;
425 IVL_init1(upperblockIVL, IVL_CHUNKED, nfront) ;
426 /*
427 -------------------
428 fill the IVL object
429 -------------------
430 */
431 for ( J = 0 ; J < nfront ; J++ ) {
432 if ( (nJ = FrontMtx_frontSize(frontmtx, J)) > 0 ) {
433 FrontMtx_columnIndices(frontmtx, J, &ncol, &colind) ;
434 if ( ncol > 0 ) {
435 mark[J] = J ;
436 count = 0 ;
437 list[count++] = J ;
438 for ( ii = nJ ; ii < ncol ; ii++ ) {
439 K = colmap[colind[ii]] ;
440 if ( mark[K] != J ) {
441 mark[K] = J ;
442 list[count++] = K ;
443 }
444 }
445 IVL_setList(upperblockIVL, J, count, list) ;
446 }
447 }
448 }
449 /*
450 ------------------------
451 free the working storage
452 ------------------------
453 */
454 IVfree(mark) ;
455 IVfree(list) ;
456
457 return(upperblockIVL) ; }
458
459 /*--------------------------------------------------------------------*/
460 /*
461 ----------------------------------------------------------------
462 purpose -- to create and return an IVL object that holds the
463 submatrix nonzero pattern for the lower triangular factor.
464
465 NOTE: this method supercedes calling IVL_mapEntries() on
466 the row adjacency structure. that gave problems when
467 pivoting forced some fronts to have no eliminated columns.
468 in some cases, solve aggregates were expected when none
469 were forthcoming.
470
471 created -- 98aug20, cca
472 ----------------------------------------------------------------
473 */
474 IVL *
FrontMtx_makeLowerBlockIVL(FrontMtx * frontmtx,IV * rowmapIV)475 FrontMtx_makeLowerBlockIVL (
476 FrontMtx *frontmtx,
477 IV *rowmapIV
478 ) {
479 int count, ii, J, K, nrow, nfront, nJ ;
480 int *rowmap, *rowind, *list, *mark ;
481 IVL *lowerblockIVL ;
482 /*
483 ---------------
484 check the input
485 ---------------
486 */
487 if ( frontmtx == NULL || rowmapIV == NULL ) {
488 fprintf(stderr, "\n fatal error in FrontMtx_makeLowerBlockIVL()"
489 "\n bad input\n") ;
490 exit(-1) ;
491 }
492 nfront = FrontMtx_nfront(frontmtx) ;
493 rowmap = IV_entries(rowmapIV) ;
494 /*
495 -----------------------------
496 set up the working storage
497 and initialize the IVL object
498 -----------------------------
499 */
500 mark = IVinit(nfront, -1) ;
501 list = IVinit(nfront, -1) ;
502 lowerblockIVL = IVL_new() ;
503 IVL_init1(lowerblockIVL, IVL_CHUNKED, nfront) ;
504 /*
505 -------------------
506 fill the IVL object
507 -------------------
508 */
509 for ( J = 0 ; J < nfront ; J++ ) {
510 if ( (nJ = FrontMtx_frontSize(frontmtx, J)) > 0 ) {
511 FrontMtx_rowIndices(frontmtx, J, &nrow, &rowind) ;
512 if ( nrow > 0 ) {
513 mark[J] = J ;
514 count = 0 ;
515 list[count++] = J ;
516 for ( ii = nJ ; ii < nrow ; ii++ ) {
517 K = rowmap[rowind[ii]] ;
518 if ( mark[K] != J ) {
519 mark[K] = J ;
520 list[count++] = K ;
521 }
522 }
523 IVL_setList(lowerblockIVL, J, count, list) ;
524 }
525 }
526 }
527 /*
528 ------------------------
529 free the working storage
530 ------------------------
531 */
532 IVfree(mark) ;
533 IVfree(list) ;
534
535 return(lowerblockIVL) ; }
536
537 /*--------------------------------------------------------------------*/
538 /*
539 ---------------------------------------------------------------
540 purpose -- to compute and return the number of floating point
541 operations to perform a solve with a single right hand side.
542
543 created -- 98sep26, cca
544 ---------------------------------------------------------------
545 */
546 int
FrontMtx_nSolveOps(FrontMtx * frontmtx)547 FrontMtx_nSolveOps (
548 FrontMtx *frontmtx
549 ) {
550 int nsolveops ;
551 /*
552 ---------------
553 check the input
554 ---------------
555 */
556 if ( frontmtx == NULL ) {
557 fprintf(stderr, "\n fatal error in FrontMtx_nSolveOps()"
558 "\n frontmtx is NULL\n") ;
559 exit(-1) ;
560 }
561 switch ( frontmtx->type ) {
562 case SPOOLES_REAL :
563 switch ( frontmtx->symmetryflag ) {
564 case SPOOLES_SYMMETRIC :
565 nsolveops = 4*frontmtx->nentU + frontmtx->nentD ;
566 break ;
567 case SPOOLES_NONSYMMETRIC :
568 nsolveops = 2*frontmtx->nentL + frontmtx->nentD
569 + 2*frontmtx->nentU ;
570 break ;
571 default :
572 fprintf(stderr, "\n fatal error in FrontMtx_nSolveOps()"
573 "\n real type, invalid symmetryflag %d\n",
574 frontmtx->symmetryflag) ;
575 exit(-1) ;
576 break ;
577 }
578 break ;
579 case SPOOLES_COMPLEX :
580 switch ( frontmtx->symmetryflag ) {
581 case SPOOLES_SYMMETRIC :
582 case SPOOLES_HERMITIAN :
583 nsolveops = 16*frontmtx->nentU + 8*frontmtx->nentD ;
584 break ;
585 case SPOOLES_NONSYMMETRIC :
586 nsolveops = 8*frontmtx->nentL + 8*frontmtx->nentD
587 + 8*frontmtx->nentU ;
588 break ;
589 default :
590 fprintf(stderr, "\n fatal error in FrontMtx_nSolveOps()"
591 "\n complex type, invalid symmetryflag %d\n",
592 frontmtx->symmetryflag) ;
593 exit(-1) ;
594 break ;
595 }
596 break ;
597 default :
598 fprintf(stderr, "\n fatal error in FrontMtx_nSolveOps()"
599 "\n invalid type %d\n", frontmtx->type) ;
600 exit(-1) ;
601 break ;
602 }
603 return(nsolveops) ; }
604
605 /*--------------------------------------------------------------------*/
606