1 /* instance.c */
2
3 #include "../SubMtx.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 ------------------------------------------
8 purpose -- fill *prowid with the row id
9 and *pcolid with the column id
10
11 created -- 98may01, cca
12 ------------------------------------------
13 */
14 void
SubMtx_ids(SubMtx * mtx,int * prowid,int * pcolid)15 SubMtx_ids (
16 SubMtx *mtx,
17 int *prowid,
18 int *pcolid
19 ) {
20 /*
21 ---------------
22 check the input
23 ---------------
24 */
25 if ( mtx == NULL || prowid == NULL || pcolid == NULL ) {
26 fprintf(stderr,
27 "\n fatal error in SubMtx_ids(%p,%p,%p)"
28 "\n bad input\n", mtx, prowid, pcolid) ;
29 exit(-1) ;
30 }
31 *prowid = mtx->rowid ;
32 *pcolid = mtx->colid ;
33
34 return ; }
35
36 /*--------------------------------------------------------------------*/
37 /*
38 -------------------------------------
39 purpose -- set the row and column ids
40
41 created -- 98may01, cca
42 -------------------------------------
43 */
44 void
SubMtx_setIds(SubMtx * mtx,int rowid,int colid)45 SubMtx_setIds (
46 SubMtx *mtx,
47 int rowid,
48 int colid
49 ) {
50 int *buffer ;
51 /*
52 ---------------
53 check the input
54 ---------------
55 */
56 if ( mtx == NULL ) {
57 fprintf(stderr,
58 "\n fatal error in SubMtx_ids(%p,%d,%d)"
59 "\n bad input\n", mtx, rowid, colid) ;
60 exit(-1) ;
61 }
62 buffer = (int *) mtx->wrkDV.vec ;
63 buffer[2] = mtx->rowid = rowid ;
64 buffer[3] = mtx->colid = colid ;
65
66 return ; }
67
68 /*--------------------------------------------------------------------*/
69 /*
70 ---------------------------------------------------
71 purpose -- fill *pnrow with the # of rows
72 and *pncol with the # of columns
73 and *pnent with the # of matrix entries
74
75 created -- 98may01, cca
76 ---------------------------------------------------
77 */
78 void
SubMtx_dimensions(SubMtx * mtx,int * pnrow,int * pncol,int * pnent)79 SubMtx_dimensions (
80 SubMtx *mtx,
81 int *pnrow,
82 int *pncol,
83 int *pnent
84 ) {
85 /*
86 ---------------
87 check the input
88 ---------------
89 */
90 if ( mtx == NULL || pnrow == NULL || pncol == NULL || pnent == NULL ) {
91 fprintf(stderr,
92 "\n fatal error in SubMtx_ids(%p,%p,%p,%p)"
93 "\n bad input\n", mtx, pnrow, pncol, pnent) ;
94 exit(-1) ;
95 }
96 *pnrow = mtx->nrow ;
97 *pncol = mtx->ncol ;
98 *pnent = mtx->nent ;
99
100 return ; }
101
102 /*--------------------------------------------------------------------*/
103 /*
104 ------------------------------------------------------------
105 purpose --
106 fill *pnrow with the number of rows
107 if prowind != NULL then
108 fill *prowind with the base location of the row indices
109 endif
110
111 created -- 98may01, cca
112 ------------------------------------------------------------
113 */
114 void
SubMtx_rowIndices(SubMtx * mtx,int * pnrow,int ** prowind)115 SubMtx_rowIndices (
116 SubMtx *mtx,
117 int *pnrow,
118 int **prowind
119 ) {
120 *pnrow = mtx->nrow ;
121 if ( prowind != NULL ) {
122 *prowind = (int *) mtx->wrkDV.vec + 7 ;
123 }
124 return ; }
125
126 /*--------------------------------------------------------------------*/
127 /*
128 ---------------------------------------------------------------
129 purpose --
130 fill *pncol with the number of column
131 if pcolind != NULL then
132 fill *prowind with the base location of the column indices
133 endif
134
135 created -- 98may01, cca
136 ---------------------------------------------------------------
137 */
138 void
SubMtx_columnIndices(SubMtx * mtx,int * pncol,int ** pcolind)139 SubMtx_columnIndices (
140 SubMtx *mtx,
141 int *pncol,
142 int **pcolind
143 ) {
144 *pncol = mtx->ncol ;
145 if ( pcolind != NULL ) {
146 *pcolind = (int *) mtx->wrkDV.vec + 7 + mtx->nrow ;
147 }
148 return ; }
149
150 /*--------------------------------------------------------------------*/
151 /*
152 ---------------------------------
153 purpose -- for dense storage
154 *pnrow with mtx->nrow
155 *pncol with mtx->ncol
156 *pinc1 with row increment
157 *pinc2 with column increment
158 *pentries with mtx->entries
159
160 created -- 98may01, cca
161 ---------------------------------
162 */
163 void
SubMtx_denseInfo(SubMtx * mtx,int * pnrow,int * pncol,int * pinc1,int * pinc2,double ** pentries)164 SubMtx_denseInfo (
165 SubMtx *mtx,
166 int *pnrow,
167 int *pncol,
168 int *pinc1,
169 int *pinc2,
170 double **pentries
171 ) {
172 double *dbuffer ;
173 int nint ;
174 /*
175 ---------------
176 check the input
177 ---------------
178 */
179 if ( mtx == NULL || pnrow == NULL || pncol == NULL
180 || pinc1 == NULL || pinc2 == NULL || pentries == NULL ) {
181 fprintf(stderr,
182 "\n fatal error in SubMtx_denseInfo(%p,%p,%p,%p,%p,%p)"
183 "\n bad input\n",
184 mtx, pnrow, pncol, pinc1, pinc2, pentries) ;
185 exit(-1) ;
186 }
187 if ( ! (SUBMTX_IS_REAL(mtx) || SUBMTX_IS_COMPLEX(mtx)) ) {
188 fprintf(stderr,
189 "\n fatal error in SubMtx_denseInfo(%p,%p,%p,%p,%p,%p)"
190 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
191 mtx, pnrow, pncol, pinc1, pinc2, pentries, mtx->type) ;
192 exit(-1) ;
193 }
194 if ( ! (SUBMTX_IS_DENSE_ROWS(mtx) || SUBMTX_IS_DENSE_COLUMNS(mtx)) ) {
195 fprintf(stderr,
196 "\n fatal error in SubMtx_denseInfo(%p,%p,%p,%p,%p,%p)"
197 "\n bad mode %d"
198 "\n must be SUBMTX_DENSE_ROWS or SUBMTX_DENSE_COLUMNS\n",
199 mtx, pnrow, pncol, pinc1, pinc2, pentries, mtx->mode) ;
200 exit(-1) ;
201 }
202 *pnrow = mtx->nrow ;
203 *pncol = mtx->ncol ;
204 if ( SUBMTX_IS_DENSE_ROWS(mtx) ) {
205 *pinc1 = mtx->ncol ;
206 *pinc2 = 1 ;
207 } else {
208 *pinc1 = 1 ;
209 *pinc2 = mtx->nrow ;
210 }
211 dbuffer = mtx->wrkDV.vec ;
212 nint = 7 + mtx->nrow + mtx->ncol ;
213 if ( sizeof(int) == sizeof(double) ) {
214 *pentries = dbuffer + nint ;
215 } else if ( 2*sizeof(int) == sizeof(double) ) {
216 *pentries = dbuffer + (nint+1)/2 ;
217 }
218 return ; }
219
220 /*--------------------------------------------------------------------*/
221 /*
222 ---------------------------------------------
223 purpose -- for sparse rows, fill
224 *pnrow with # of rows
225 *pnent with # of indices
226 *psizes with sizes[] of rows
227 *pindices with indices[] for column indices
228 *pentries with entries[] for matrix entries
229
230 created -- 98may01, cca
231 ---------------------------------------------
232 */
233 void
SubMtx_sparseRowsInfo(SubMtx * mtx,int * pnrow,int * pnent,int ** psizes,int ** pindices,double ** pentries)234 SubMtx_sparseRowsInfo (
235 SubMtx *mtx,
236 int *pnrow,
237 int *pnent,
238 int **psizes,
239 int **pindices,
240 double **pentries
241 ) {
242 double *dbuffer ;
243 int nint ;
244 int *ibuffer ;
245 /*
246 ---------------
247 check the input
248 ---------------
249 */
250 if ( mtx == NULL || pnrow == NULL || pnent == NULL
251 || psizes == NULL || pindices == NULL || pentries == NULL ) {
252 fprintf(stderr,
253 "\n fatal error in SubMtx_sparseRowsInfo(%p,%p,%p,%p,%p,%p)"
254 "\n bad input\n",
255 mtx, pnrow, pnent, psizes, pindices, pentries) ;
256 exit(-1) ;
257 }
258 if ( ! (SUBMTX_IS_REAL(mtx) || SUBMTX_IS_COMPLEX(mtx)) ) {
259 fprintf(stderr,
260 "\n fatal error in SubMtx_sparseRowsInfo(%p,%p,%p,%p,%p,%p)"
261 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
262 mtx, pnrow, pnent, psizes, pindices, pentries, mtx->type) ;
263 exit(-1) ;
264 }
265 if ( ! SUBMTX_IS_SPARSE_ROWS(mtx) ) {
266 fprintf(stderr,
267 "\n fatal error in SubMtx_sparseRowsInfo(%p,%p,%p,%p,%p,%p)"
268 "\n bad mode %d, must be SUBMTX_SPARSE_ROWS\n",
269 mtx, pnrow, pnent, psizes, pindices, pentries, mtx->mode) ;
270 exit(-1) ;
271 }
272 *pnrow = mtx->nrow ;
273 *pnent = mtx->nent ;
274 dbuffer = mtx->wrkDV.vec ;
275 ibuffer = (int *) dbuffer ;
276 nint = 7 + mtx->nrow + mtx->ncol ;
277 *psizes = ibuffer + nint ;
278 nint += mtx->nrow ;
279 *pindices = ibuffer + nint ;
280 nint += mtx->nent ;
281 if ( sizeof(int) == sizeof(double) ) {
282 *pentries = dbuffer + nint ;
283 } else if ( 2*sizeof(int) == sizeof(double) ) {
284 *pentries = dbuffer + (nint+1)/2 ;
285 }
286 return ; }
287
288 /*--------------------------------------------------------------------*/
289 /*
290 ----------------------------------------------
291 purpose -- for sparse columns, fill
292 *pncol with # of columns
293 *pnent with # of matrix entries
294 *psizes with sizes[ncol], column sizes
295 *pindices with indices[nent], matrix row ids
296 *pentries with entries[nent], matrix entries
297
298 created -- 98may01, cca
299 ----------------------------------------------
300 */
301 void
SubMtx_sparseColumnsInfo(SubMtx * mtx,int * pncol,int * pnent,int ** psizes,int ** pindices,double ** pentries)302 SubMtx_sparseColumnsInfo (
303 SubMtx *mtx,
304 int *pncol,
305 int *pnent,
306 int **psizes,
307 int **pindices,
308 double **pentries
309 ) {
310 double *dbuffer ;
311 int nint ;
312 int *ibuffer ;
313 /*
314 ---------------
315 check the input
316 ---------------
317 */
318 if ( mtx == NULL || pncol == NULL || pnent == NULL
319 || psizes == NULL || pindices == NULL || pentries == NULL ) {
320 fprintf(stderr,
321 "\n fatal error in SubMtx_sparseColumnsInfo(%p,%p,%p,%p,%p,%p)"
322 "\n bad input\n",
323 mtx, pncol, pnent, psizes, pindices, pentries) ;
324 exit(-1) ;
325 }
326 if ( ! (SUBMTX_IS_REAL(mtx) || SUBMTX_IS_COMPLEX(mtx)) ) {
327 fprintf(stderr,
328 "\n fatal error in SubMtx_sparseColumnsInfo(%p,%p,%p,%p,%p,%p)"
329 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
330 mtx, pncol, pnent, psizes, pindices, pentries, mtx->type) ;
331 exit(-1) ;
332 }
333 if ( ! SUBMTX_IS_SPARSE_COLUMNS(mtx) ) {
334 fprintf(stderr,
335 "\n fatal error in SubMtx_sparseColumnsInfo(%p,%p,%p,%p,%p,%p)"
336 "\n bad mode %d"
337 "\n must be SUBMTX_SPARSE_COLUMNS\n",
338 mtx, pncol, pnent, psizes, pindices, pentries, mtx->mode) ;
339 exit(-1) ;
340 }
341 *pncol = mtx->ncol ;
342 *pnent = mtx->nent ;
343 dbuffer = mtx->wrkDV.vec ;
344 ibuffer = (int *) dbuffer ;
345 nint = 7 + mtx->nrow + mtx->ncol ;
346 *psizes = ibuffer + nint ;
347 nint += mtx->ncol ;
348 *pindices = ibuffer + nint ;
349 nint += mtx->nent ;
350 if ( sizeof(int) == sizeof(double) ) {
351 *pentries = dbuffer + nint ;
352 } else if ( 2*sizeof(int) == sizeof(double) ) {
353 *pentries = dbuffer + (nint+1)/2 ;
354 }
355 return ; }
356
357 /*--------------------------------------------------------------------*/
358 /*
359 ----------------------------------------------------
360 purpose -- for sparse triples, fill
361 *pnent with # of matrix entries
362 *prowids with rowids[nent], row ids of entries
363 *prowids with colids[nent], column ids of entries
364 *pentries with entries[nent], matrix entries
365
366 created -- 98may01, cca
367 ----------------------------------------------------
368 */
369 void
SubMtx_sparseTriplesInfo(SubMtx * mtx,int * pnent,int ** prowids,int ** pcolids,double ** pentries)370 SubMtx_sparseTriplesInfo (
371 SubMtx *mtx,
372 int *pnent,
373 int **prowids,
374 int **pcolids,
375 double **pentries
376 ) {
377 double *dbuffer ;
378 int nint ;
379 int *ibuffer ;
380 /*
381 ---------------
382 check the input
383 ---------------
384 */
385 if ( mtx == NULL || pnent == NULL
386 || prowids == NULL || pcolids == NULL || pentries == NULL ) {
387 fprintf(stderr,
388 "\n fatal error in SubMtx_sparseTriplesInfo(%p,%p,%p,%p,%p)"
389 "\n bad input\n",
390 mtx, pnent, prowids, pcolids, pentries) ;
391 exit(-1) ;
392 }
393 if ( ! (SUBMTX_IS_REAL(mtx) || SUBMTX_IS_COMPLEX(mtx)) ) {
394 fprintf(stderr,
395 "\n fatal error in SubMtx_sparseTriplesInfo(%p,%p,%p,%p,%p)"
396 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
397 mtx, pnent, prowids, pcolids, pentries, mtx->type) ;
398 exit(-1) ;
399 }
400 if ( ! SUBMTX_IS_SPARSE_TRIPLES(mtx) ) {
401 fprintf(stderr,
402 "\n fatal error in SubMtx_sparseTriplesInfo(%p,%p,%p,%p,%p)"
403 "\n bad mode %d"
404 "\n must be SUBMTX_SPARSE_TRIPLES\n",
405 mtx, pnent, prowids, pcolids, pentries, mtx->mode) ;
406 exit(-1) ;
407 }
408 *pnent = mtx->nent ;
409 dbuffer = mtx->wrkDV.vec ;
410 ibuffer = (int *) dbuffer ;
411 nint = 7 + mtx->nrow + mtx->ncol ;
412 *prowids = ibuffer + nint ;
413 nint += mtx->nent ;
414 *pcolids = ibuffer + nint ;
415 nint += mtx->nent ;
416 if ( sizeof(int) == sizeof(double) ) {
417 *pentries = dbuffer + nint ;
418 } else if ( 2*sizeof(int) == sizeof(double) ) {
419 *pentries = dbuffer + (nint+1)/2 ;
420 }
421 return ; }
422
423 /*--------------------------------------------------------------------*/
424 /*
425 -----------------------------------------------------------
426 purpose -- for dense subrows, fill
427 *pnrow with # of rows
428 *pnent with # of matrix entries
429 *pfirstlocs with firstlocs[nrow], column of first nonzero
430 *psizes with sizes[nrow], number of nonzero columns
431 *pentries with entries[nent], matrix entries
432
433 created -- 98may01, cca
434 -----------------------------------------------------------
435 */
436 void
SubMtx_denseSubrowsInfo(SubMtx * mtx,int * pnrow,int * pnent,int ** pfirstlocs,int ** psizes,double ** pentries)437 SubMtx_denseSubrowsInfo (
438 SubMtx *mtx,
439 int *pnrow,
440 int *pnent,
441 int **pfirstlocs,
442 int **psizes,
443 double **pentries
444 ) {
445 double *dbuffer ;
446 int nint ;
447 int *ibuffer ;
448 /*
449 ---------------
450 check the input
451 ---------------
452 */
453 if ( mtx == NULL || pnrow == NULL || pnent == NULL
454 || pfirstlocs == NULL || psizes == NULL || pentries == NULL ) {
455 fprintf(stderr,
456 "\n fatal error in SubMtx_denseSubrowsInfo(%p,%p,%p,%p,%p,%p)"
457 "\n bad input\n",
458 mtx, pnrow, pnent, pfirstlocs, psizes, pentries) ;
459 if ( mtx != NULL ) {
460 SubMtx_writeForHumanEye(mtx, stderr) ;
461 }
462 exit(-1) ;
463 }
464 if ( ! (SUBMTX_IS_REAL(mtx) || SUBMTX_IS_COMPLEX(mtx)) ) {
465 fprintf(stderr,
466 "\n fatal error in SubMtx_denseSubrowsInfo(%p,%p,%p,%p,%p,%p)"
467 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
468 mtx, pnrow, pnent, pfirstlocs, psizes, pentries, mtx->type) ;
469 exit(-1) ;
470 }
471 if ( ! SUBMTX_IS_DENSE_SUBROWS(mtx) ) {
472 fprintf(stderr,
473 "\n fatal error in SubMtx_denseSubrowsInfo(%p,%p,%p,%p,%p,%p)"
474 "\n bad mode %d"
475 "\n must be SUBMTX_DENSE_SUBROWS\n",
476 mtx, pnrow, pnent, pfirstlocs, psizes, pentries, mtx->mode) ;
477 exit(-1) ;
478 }
479 *pnrow = mtx->nrow ;
480 *pnent = mtx->nent ;
481 dbuffer = mtx->wrkDV.vec ;
482 ibuffer = (int *) dbuffer ;
483 nint = 7 + mtx->nrow + mtx->ncol ;
484 *pfirstlocs = ibuffer + nint ;
485 nint += mtx->nrow ;
486 *psizes = ibuffer + nint ;
487 nint += mtx->nrow ;
488 if ( sizeof(int) == sizeof(double) ) {
489 *pentries = dbuffer + nint ;
490 } else if ( 2*sizeof(int) == sizeof(double) ) {
491 *pentries = dbuffer + (nint+1)/2 ;
492 }
493
494 return ; }
495
496 /*--------------------------------------------------------------------*/
497 /*
498 -----------------------------------------------------------
499 purpose -- for dense subcolumns, fill
500 *pncol with # of columns
501 *pnent with # of matrix entries
502 *pfirstlocs with firstlocs[ncol], row of first nonzero
503 *psizes with sizes[ncol], number of nonzero rows
504 *pentries with entries[nent], matrix entries
505
506 created -- 98may01, cca
507 -----------------------------------------------------------
508 */
509 void
SubMtx_denseSubcolumnsInfo(SubMtx * mtx,int * pncol,int * pnent,int ** pfirstlocs,int ** psizes,double ** pentries)510 SubMtx_denseSubcolumnsInfo (
511 SubMtx *mtx,
512 int *pncol,
513 int *pnent,
514 int **pfirstlocs,
515 int **psizes,
516 double **pentries
517 ) {
518 double *dbuffer ;
519 int nint ;
520 int *ibuffer ;
521 /*
522 ---------------
523 check the input
524 ---------------
525 */
526 if ( mtx == NULL
527 || pfirstlocs == NULL || psizes == NULL || pentries == NULL ) {
528 fprintf(stderr,
529 "\n fatal error in SubMtx_denseSubcolumnsInfo(%p,%p,%p,%p,%p,%p)"
530 "\n bad input\n",
531 mtx, pncol, pnent, pfirstlocs, psizes, pentries) ;
532 exit(-1) ;
533 }
534 if ( ! (SUBMTX_IS_REAL(mtx) || SUBMTX_IS_COMPLEX(mtx)) ) {
535 fprintf(stderr,
536 "\n fatal error in SubMtx_denseSubcolumsInfo(%p,%p,%p,%p,%p,%p)"
537 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
538 mtx, pncol, pnent, pfirstlocs, psizes, pentries, mtx->type) ;
539 exit(-1) ;
540 }
541 if ( ! SUBMTX_IS_DENSE_SUBCOLUMNS(mtx) ) {
542 fprintf(stderr,
543 "\n fatal error in SubMtx_denseSubcolumnsInfo(%p,%p,%p,%p,%p,%p)"
544 "\n bad mode %d"
545 "\n must be SUBMTX_DENSE_SUBCOLUMNS\n",
546 mtx, pncol, pnent, pfirstlocs, psizes, pentries, mtx->mode) ;
547 exit(-1) ;
548 }
549 *pncol = mtx->ncol ;
550 *pnent = mtx->nent ;
551 dbuffer = mtx->wrkDV.vec ;
552 ibuffer = (int *) dbuffer ;
553 nint = 7 + mtx->nrow + mtx->ncol ;
554 *pfirstlocs = ibuffer + nint ;
555 nint += mtx->ncol ;
556 *psizes = ibuffer + nint ;
557 nint += mtx->ncol ;
558 if ( sizeof(int) == sizeof(double) ) {
559 *pentries = dbuffer + nint ;
560 } else if ( 2*sizeof(int) == sizeof(double) ) {
561 *pentries = dbuffer + (nint+1)/2 ;
562 }
563
564 return ; }
565
566 /*--------------------------------------------------------------------*/
567 /*
568 ------------------------------------------------
569 purpose -- for a diagonal matrix, fill
570 *pncol with # of columns
571 *pentries with entries[nent], matrix entries
572
573 created -- 98may01, cca
574 ------------------------------------------------
575 */
576 void
SubMtx_diagonalInfo(SubMtx * mtx,int * pncol,double ** pentries)577 SubMtx_diagonalInfo (
578 SubMtx *mtx,
579 int *pncol,
580 double **pentries
581 ) {
582 double *dbuffer ;
583 int nint ;
584 int *ibuffer ;
585 /*
586 ---------------
587 check the input
588 ---------------
589 */
590 if ( mtx == NULL || pncol == NULL || pentries == NULL ) {
591 fprintf(stderr,
592 "\n fatal error in SubMtx_diagonalInfo(%p,%p,%p)"
593 "\n bad input\n",
594 mtx, pncol, pentries) ;
595 exit(-1) ;
596 }
597 if ( ! (SUBMTX_IS_REAL(mtx) || SUBMTX_IS_COMPLEX(mtx)) ) {
598 fprintf(stderr,
599 "\n fatal error in SubMtx_diagonalInfo(%p,%p,%p)"
600 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
601 mtx, pncol, pentries, mtx->type) ;
602 exit(-1) ;
603 }
604 if ( ! SUBMTX_IS_DIAGONAL(mtx) ) {
605 fprintf(stderr,
606 "\n fatal error in SubMtx_diagonalInfo(%p,%p,%p)"
607 "\n bad mode %d"
608 "\n must be SUBMTX_DIAGONAL\n",
609 mtx, pncol, pentries, mtx->mode) ;
610 exit(-1) ;
611 }
612 *pncol = mtx->ncol ;
613 dbuffer = mtx->wrkDV.vec ;
614 ibuffer = (int *) dbuffer ;
615 nint = 7 + mtx->nrow + mtx->ncol ;
616 if ( sizeof(int) == sizeof(double) ) {
617 *pentries = dbuffer + nint ;
618 } else if ( 2*sizeof(int) == sizeof(double) ) {
619 *pentries = dbuffer + (nint+1)/2 ;
620 }
621 return ; }
622
623 /*--------------------------------------------------------------------*/
624 /*
625 ------------------------------------------------------
626 purpose -- for a block diagonal symmetric matrix, fill
627 *pncol with # of columns
628 *pnent with # of entries
629 *ppivotsizes with pivotsizes[ncol]
630 *pentries with entries[nent], matrix entries
631
632 created -- 98may01, cca
633 ------------------------------------------------------
634 */
635 void
SubMtx_blockDiagonalInfo(SubMtx * mtx,int * pncol,int * pnent,int ** ppivotsizes,double ** pentries)636 SubMtx_blockDiagonalInfo (
637 SubMtx *mtx,
638 int *pncol,
639 int *pnent,
640 int **ppivotsizes,
641 double **pentries
642 ) {
643 double *dbuffer ;
644 int nint ;
645 int *ibuffer ;
646 /*
647 ---------------
648 check the input
649 ---------------
650 */
651 if ( mtx == NULL
652 || pncol == NULL || pnent == NULL
653 || ppivotsizes == NULL || pentries == NULL ) {
654 fprintf(stderr,
655 "\n fatal error in SubMtx_blockDiagonalInfo(%p,%p,%p,%p,%p)"
656 "\n bad input\n",
657 mtx, pncol, pnent, ppivotsizes, pentries) ;
658 exit(-1) ;
659 }
660 if ( ! (SUBMTX_IS_REAL(mtx) || SUBMTX_IS_COMPLEX(mtx)) ) {
661 fprintf(stderr,
662 "\n fatal error in SubMtx_blockDiagonalInfo(%p,%p,%p,%p,%p)"
663 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
664 mtx, pncol, pnent, ppivotsizes, pentries, mtx->type) ;
665 exit(-1) ;
666 }
667 if ( ! (SUBMTX_IS_BLOCK_DIAGONAL_SYM(mtx)
668 || SUBMTX_IS_BLOCK_DIAGONAL_HERM(mtx)) ) {
669 fprintf(stderr,
670 "\n fatal error in SubMtx_blockDiagonalInfo(%p,%p,%p,%p,%p)"
671 "\n bad mode %d"
672 "\n must be SUBMTX_BLOCK_DIAGONAL_SYM or SUBMTX_BLOCK_DIAGONAL_HERM \n",
673 mtx, pncol, pnent, ppivotsizes, pentries, mtx->mode) ;
674 exit(-1) ;
675 }
676 *pncol = mtx->ncol ;
677 *pnent = mtx->nent ;
678 dbuffer = mtx->wrkDV.vec ;
679 ibuffer = (int *) dbuffer ;
680 nint = 7 + 2*mtx->nrow ;
681 *ppivotsizes = ibuffer + nint ;
682 nint += mtx->nrow ;
683 if ( sizeof(int) == sizeof(double) ) {
684 *pentries = dbuffer + nint ;
685 } else if ( 2*sizeof(int) == sizeof(double) ) {
686 *pentries = dbuffer + (nint+1)/2 ;
687 }
688 return ; }
689
690 /*--------------------------------------------------------------------*/
691 /*
692 -------------------------------------------------------
693 purpose -- to find matrix entry (irow,jcol) if present.
694
695 return value --
696 if entry (irow,jcol) is not present then
697 *pValue is 0.0
698 return value is -1
699 else entry (irow,jcol) is present then
700 *pValue is the matrix entry
701 return value is offset into entries array
702 endif
703
704 created -- 98may01, cca
705 -------------------------------------------------------
706 */
707 int
SubMtx_realEntry(SubMtx * mtx,int irow,int jcol,double * pValue)708 SubMtx_realEntry (
709 SubMtx *mtx,
710 int irow,
711 int jcol,
712 double *pValue
713 ) {
714 /*
715 ---------------
716 check the input
717 ---------------
718 */
719 if ( mtx == NULL || irow < 0 || irow >= mtx->nrow || jcol < 0
720 || jcol >= mtx->ncol || pValue == NULL ) {
721 fprintf(stderr,
722 "\n fatal error in SubMtx_realEntry(%p,%d,%d,%p)"
723 "\n bad input\n", mtx, irow, jcol, pValue) ;
724 exit(-1) ;
725 }
726 if ( ! SUBMTX_IS_REAL(mtx) ) {
727 fprintf(stderr,
728 "\n fatal error in SubMtx_realEntry(%p,%d,%d,%p)"
729 "\n bad type %d, must be SPOOLES_REAL\n",
730 mtx, irow, jcol, pValue, mtx->type) ;
731 exit(-1) ;
732 }
733 *pValue = 0 ;
734 switch ( mtx->mode ) {
735 case SUBMTX_DENSE_ROWS :
736 case SUBMTX_DENSE_COLUMNS : {
737 double *entries ;
738 int inc1, inc2, ncol, nrow, offset ;
739
740 SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
741 if ( irow < 0 || irow >= nrow || jcol < 0 || jcol >= ncol ) {
742 return(-1) ;
743 }
744 offset = irow*inc1 + jcol*inc2 ;
745 *pValue = entries[offset] ;
746 return(offset) ;
747 } break ;
748 case SUBMTX_SPARSE_ROWS : {
749 double *entries ;
750 int ii, jj, nent, nrow, offset, *indices, *sizes ;
751
752 SubMtx_sparseRowsInfo(mtx, &nrow, &nent, &sizes, &indices, &entries) ;
753 if ( irow < 0 || irow >= nrow ) {
754 return(-1) ;
755 }
756 for ( ii = offset = 0 ; ii < irow ; ii++ ) {
757 offset += sizes[ii] ;
758 }
759 for ( ii = 0, jj = offset ; ii < sizes[irow] ; ii++, jj++ ) {
760 if ( indices[jj] == jcol ) {
761 *pValue = entries[jj] ;
762 return(jj) ;
763 }
764 }
765 return(-1) ;
766 } break ;
767 case SUBMTX_SPARSE_COLUMNS : {
768 double *entries ;
769 int ii, jj, nent, ncol, offset, *indices, *sizes ;
770
771 SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
772 &sizes, &indices, &entries) ;
773 if ( jcol < 0 || jcol >= ncol ) {
774 return(-1) ;
775 }
776 for ( ii = offset = 0 ; ii < jcol ; ii++ ) {
777 offset += sizes[ii] ;
778 }
779 for ( ii = 0, jj = offset ; ii < sizes[jcol] ; ii++, jj++ ) {
780 if ( indices[jj] == irow ) {
781 *pValue = entries[jj] ;
782 return(jj) ;
783 }
784 }
785 return(-1) ;
786 } break ;
787 case SUBMTX_SPARSE_TRIPLES : {
788 double *entries ;
789 int ii, nent, *colids, *rowids ;
790
791 SubMtx_sparseTriplesInfo(mtx, &nent, &rowids, &colids, &entries) ;
792 for ( ii = 0 ; ii < nent ; ii++ ) {
793 if ( irow == rowids[ii] && jcol == colids[ii] ) {
794 *pValue = entries[ii] ;
795 return(ii) ;
796 }
797 }
798 return(-1) ;
799 } break ;
800 case SUBMTX_DENSE_SUBROWS : {
801 double *entries ;
802 int ii, joff, nent, nrow, offset, *firstlocs, *sizes ;
803
804 SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
805 &firstlocs, &sizes, &entries) ;
806 if ( irow < 0 || irow >= nrow || sizes[irow] == 0 ) {
807 return(-1) ;
808 }
809 for ( ii = offset = 0 ; ii < irow ; ii++ ) {
810 offset += sizes[ii] ;
811 }
812 if ( 0 <= (joff = jcol - firstlocs[irow]) && joff < sizes[irow] ) {
813 offset += joff ;
814 *pValue = entries[offset] ;
815 return(offset) ;
816 }
817 return(-1) ;
818 } break ;
819 case SUBMTX_DENSE_SUBCOLUMNS : {
820 double *entries ;
821 int ii, ioff, nent, ncol, offset, *firstlocs, *sizes ;
822
823 SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
824 &firstlocs, &sizes, &entries) ;
825 if ( jcol < 0 || jcol >= ncol || sizes[jcol] == 0 ) {
826 return(-1) ;
827 }
828 for ( ii = offset = 0 ; ii < jcol ; ii++ ) {
829 offset += sizes[ii] ;
830 }
831 if ( 0 <= (ioff = irow - firstlocs[jcol]) && ioff < sizes[jcol] ) {
832 offset += ioff ;
833 *pValue = entries[offset] ;
834 return(offset) ;
835 }
836 return(-1) ;
837 } break ;
838 case SUBMTX_DIAGONAL : {
839 double *entries ;
840 int ncol ;
841
842 if ( irow < 0 || jcol < 0 || irow != jcol ) {
843 return(-1) ;
844 }
845 SubMtx_diagonalInfo(mtx, &ncol, &entries) ;
846 if ( irow >= ncol || jcol >= ncol ) {
847 return(-1) ;
848 }
849 *pValue = entries[irow] ;
850 return(irow) ;
851 } break ;
852 case SUBMTX_BLOCK_DIAGONAL_SYM : {
853 double *entries ;
854 int ii, ipivot, jrow, kk, m, ncol, nent, size ;
855 int *pivotsizes ;
856
857 if ( irow < 0 || jcol < 0 ) {
858 return(-1) ;
859 }
860 if ( irow > jcol ) {
861 ii = irow ;
862 irow = jcol ;
863 jcol = ii ;
864 }
865 SubMtx_blockDiagonalInfo(mtx, &ncol, &nent, &pivotsizes, &entries) ;
866 if ( irow >= ncol || jcol >= ncol ) {
867 return(-1) ;
868 }
869 for ( jrow = ipivot = kk = 0 ; jrow <= irow ; ipivot++ ) {
870 size = m = pivotsizes[ipivot] ;
871 for ( ii = 0 ; ii < m ; ii++, jrow++ ) {
872 if ( jrow == irow ) {
873 if ( jcol - irow > m - ii - 1 ) {
874 return(-1) ;
875 } else {
876 kk += jcol - irow ;
877 *pValue = entries[kk] ;
878 return(kk) ;
879 }
880 } else {
881 kk += size-- ;
882 }
883 }
884 }
885 return(kk) ;
886 } break ;
887 case SUBMTX_BLOCK_DIAGONAL_HERM : {
888 double sign ;
889 double *entries ;
890 int ii, ipivot, jrow, kk, m, ncol, nent, size ;
891 int *pivotsizes ;
892
893 if ( irow < 0 || jcol < 0 ) {
894 return(-1) ;
895 }
896 if ( irow > jcol ) {
897 ii = irow ;
898 irow = jcol ;
899 jcol = ii ;
900 sign = -1.0 ;
901 } else {
902 sign = 1.0 ;
903 }
904 SubMtx_blockDiagonalInfo(mtx, &ncol, &nent, &pivotsizes, &entries) ;
905 if ( irow >= ncol || jcol >= ncol ) {
906 return(-1) ;
907 }
908 for ( jrow = ipivot = kk = 0 ; jrow <= irow ; ipivot++ ) {
909 size = m = pivotsizes[ipivot] ;
910 for ( ii = 0 ; ii < m ; ii++, jrow++ ) {
911 if ( jrow == irow ) {
912 if ( jcol - irow > m - ii - 1 ) {
913 return(-1) ;
914 } else {
915 kk += jcol - irow ;
916 *pValue = entries[kk] ;
917 return(kk) ;
918 }
919 } else {
920 kk += size-- ;
921 }
922 }
923 }
924 return(kk) ;
925 } break ;
926 default :
927 fprintf(stderr,
928 "\n fatal error in SubMtx_realEntry(%p,%d,%d,%p)"
929 "\n bad mode %d", mtx, irow, jcol, pValue, mtx->mode) ;
930 exit(-1) ;
931 break ;
932 }
933 return(-1) ; }
934
935 /*--------------------------------------------------------------------*/
936 /*
937 -------------------------------------------------------
938 purpose -- to find matrix entry (irow,jcol) if present.
939
940 return value --
941 if entry (irow,jcol) is not present then
942 *pReal and *pImag are 0.0
943 return value is -1
944 else entry (irow,jcol) is present then
945 (*pReal,*pImag) is the matrix entry
946 return value is offset into entries array
947 endif
948
949 created -- 98may01, cca
950 -------------------------------------------------------
951 */
952 int
SubMtx_complexEntry(SubMtx * mtx,int irow,int jcol,double * pReal,double * pImag)953 SubMtx_complexEntry (
954 SubMtx *mtx,
955 int irow,
956 int jcol,
957 double *pReal,
958 double *pImag
959 ) {
960 /*
961 ---------------
962 check the input
963 ---------------
964 */
965 if ( mtx == NULL || irow < 0 || irow >= mtx->nrow || jcol < 0
966 || jcol >= mtx->ncol || pReal == NULL || pImag == NULL ) {
967 fprintf(stderr,
968 "\n fatal error in SubMtx_complexEntry(%p,%d,%d,%p,%p)"
969 "\n bad input\n", mtx, irow, jcol, pReal, pImag) ;
970 exit(-1) ;
971 }
972 if ( ! SUBMTX_IS_COMPLEX(mtx) ) {
973 fprintf(stderr,
974 "\n fatal error in SubMtx_complexEntry(%p,%d,%d,%p,%p)"
975 "\n bad type %d, must be SPOOLES_COMPLEX\n",
976 mtx, irow, jcol, pReal, pImag, mtx->type) ;
977 exit(-1) ;
978 }
979 *pReal = *pImag = 0 ;
980 switch ( mtx->mode ) {
981 case SUBMTX_DENSE_ROWS :
982 case SUBMTX_DENSE_COLUMNS : {
983 double *entries ;
984 int inc1, inc2, ncol, nrow, offset ;
985
986 SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
987 if ( irow < 0 || irow >= nrow || jcol < 0 || jcol >= ncol ) {
988 return(-1) ;
989 }
990 offset = irow*inc1 + jcol*inc2 ;
991 *pReal = entries[2*offset] ;
992 *pImag = entries[2*offset+1] ;
993 return(offset) ;
994 } break ;
995 case SUBMTX_SPARSE_ROWS : {
996 double *entries ;
997 int ii, jj, nent, nrow, offset, *indices, *sizes ;
998
999 SubMtx_sparseRowsInfo(mtx, &nrow, &nent, &sizes, &indices, &entries) ;
1000 if ( irow < 0 || irow >= nrow ) {
1001 return(-1) ;
1002 }
1003 for ( ii = offset = 0 ; ii < irow ; ii++ ) {
1004 offset += sizes[ii] ;
1005 }
1006 for ( ii = 0, jj = offset ; ii < sizes[irow] ; ii++, jj++ ) {
1007 if ( indices[jj] == jcol ) {
1008 *pReal = entries[2*jj] ;
1009 *pImag = entries[2*jj+1] ;
1010 return(jj) ;
1011 }
1012 }
1013 return(-1) ;
1014 } break ;
1015 case SUBMTX_SPARSE_COLUMNS : {
1016 double *entries ;
1017 int ii, jj, nent, ncol, offset, *indices, *sizes ;
1018
1019 SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
1020 &sizes, &indices, &entries) ;
1021 if ( jcol < 0 || jcol >= ncol ) {
1022 return(-1) ;
1023 }
1024 for ( ii = offset = 0 ; ii < jcol ; ii++ ) {
1025 offset += sizes[ii] ;
1026 }
1027 for ( ii = 0, jj = offset ; ii < sizes[jcol] ; ii++, jj++ ) {
1028 if ( indices[jj] == irow ) {
1029 *pReal = entries[2*jj] ;
1030 *pImag = entries[2*jj+1] ;
1031 return(jj) ;
1032 }
1033 }
1034 return(-1) ;
1035 } break ;
1036 case SUBMTX_SPARSE_TRIPLES : {
1037 double *entries ;
1038 int ii, nent, *colids, *rowids ;
1039
1040 SubMtx_sparseTriplesInfo(mtx, &nent, &rowids, &colids, &entries) ;
1041 for ( ii = 0 ; ii < nent ; ii++ ) {
1042 if ( irow == rowids[ii] && jcol == colids[ii] ) {
1043 *pReal = entries[2*ii] ;
1044 *pImag = entries[2*ii+1] ;
1045 return(ii) ;
1046 }
1047 }
1048 return(-1) ;
1049 } break ;
1050 case SUBMTX_DENSE_SUBROWS : {
1051 double *entries ;
1052 int ii, joff, nent, nrow, offset, *firstlocs, *sizes ;
1053
1054 SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
1055 &firstlocs, &sizes, &entries) ;
1056 if ( irow < 0 || irow >= nrow || sizes[irow] == 0 ) {
1057 return(-1) ;
1058 }
1059 for ( ii = offset = 0 ; ii < irow ; ii++ ) {
1060 offset += sizes[ii] ;
1061 }
1062 if ( 0 <= (joff = jcol - firstlocs[irow]) && joff < sizes[irow] ) {
1063 offset += joff ;
1064 *pReal = entries[2*offset] ;
1065 *pImag = entries[2*offset+1] ;
1066 return(offset) ;
1067 }
1068 return(-1) ;
1069 } break ;
1070 case SUBMTX_DENSE_SUBCOLUMNS : {
1071 double *entries ;
1072 int ii, ioff, nent, ncol, offset, *firstlocs, *sizes ;
1073
1074 SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
1075 &firstlocs, &sizes, &entries) ;
1076 if ( jcol < 0 || jcol >= ncol || sizes[jcol] == 0 ) {
1077 return(-1) ;
1078 }
1079 for ( ii = offset = 0 ; ii < jcol ; ii++ ) {
1080 offset += sizes[ii] ;
1081 }
1082 if ( 0 <= (ioff = irow - firstlocs[jcol]) && ioff < sizes[jcol] ) {
1083 offset += ioff ;
1084 *pReal = entries[2*offset] ;
1085 *pImag = entries[2*offset+1] ;
1086 return(offset) ;
1087 }
1088 return(-1) ;
1089 } break ;
1090 case SUBMTX_DIAGONAL : {
1091 double *entries ;
1092 int ncol ;
1093
1094 if ( irow < 0 || jcol < 0 || irow != jcol ) {
1095 return(-1) ;
1096 }
1097 SubMtx_diagonalInfo(mtx, &ncol, &entries) ;
1098 if ( irow >= ncol || jcol >= ncol ) {
1099 return(-1) ;
1100 }
1101 *pReal = entries[2*irow] ;
1102 *pImag = entries[2*irow+1] ;
1103 return(irow) ;
1104 } break ;
1105 case SUBMTX_BLOCK_DIAGONAL_SYM : {
1106 double *entries ;
1107 int ii, ipivot, jrow, kk, m, ncol, nent, size ;
1108 int *pivotsizes ;
1109
1110 if ( irow < 0 || jcol < 0 ) {
1111 return(-1) ;
1112 }
1113 if ( irow > jcol ) {
1114 ii = irow ;
1115 irow = jcol ;
1116 jcol = ii ;
1117 }
1118 SubMtx_blockDiagonalInfo(mtx, &ncol, &nent, &pivotsizes, &entries) ;
1119 if ( irow >= ncol || jcol >= ncol ) {
1120 return(-1) ;
1121 }
1122 for ( jrow = ipivot = kk = 0 ; jrow <= irow ; ipivot++ ) {
1123 size = m = pivotsizes[ipivot] ;
1124 for ( ii = 0 ; ii < m ; ii++, jrow++ ) {
1125 if ( jrow == irow ) {
1126 if ( jcol - irow > m - ii - 1 ) {
1127 return(-1) ;
1128 } else {
1129 kk += jcol - irow ;
1130 *pReal = entries[2*kk] ;
1131 *pImag = entries[2*kk+1] ;
1132 return(kk) ;
1133 }
1134 } else {
1135 kk += size-- ;
1136 }
1137 }
1138 }
1139 return(kk) ;
1140 } break ;
1141 case SUBMTX_BLOCK_DIAGONAL_HERM : {
1142 double sign ;
1143 double *entries ;
1144 int ii, ipivot, jrow, kk, m, ncol, nent, size ;
1145 int *pivotsizes ;
1146
1147 if ( irow < 0 || jcol < 0 ) {
1148 return(-1) ;
1149 }
1150 if ( irow > jcol ) {
1151 ii = irow ;
1152 irow = jcol ;
1153 jcol = ii ;
1154 sign = -1.0 ;
1155 } else {
1156 sign = 1.0 ;
1157 }
1158 SubMtx_blockDiagonalInfo(mtx, &ncol, &nent, &pivotsizes, &entries) ;
1159 if ( irow >= ncol || jcol >= ncol ) {
1160 return(-1) ;
1161 }
1162 for ( jrow = ipivot = kk = 0 ; jrow <= irow ; ipivot++ ) {
1163 size = m = pivotsizes[ipivot] ;
1164 for ( ii = 0 ; ii < m ; ii++, jrow++ ) {
1165 if ( jrow == irow ) {
1166 if ( jcol - irow > m - ii - 1 ) {
1167 return(-1) ;
1168 } else {
1169 kk += jcol - irow ;
1170 *pReal = entries[2*kk] ;
1171 *pImag = sign*entries[2*kk+1] ;
1172 return(kk) ;
1173 }
1174 } else {
1175 kk += size-- ;
1176 }
1177 }
1178 }
1179 return(kk) ;
1180 } break ;
1181 default :
1182 fprintf(stderr,
1183 "\n fatal error in SubMtx_complexEntry(%p,%d,%d,%p,%p)"
1184 "\n bad mode %d", mtx, irow, jcol, pReal, pImag, mtx->mode) ;
1185 exit(-1) ;
1186 break ;
1187 }
1188 return(-1) ; }
1189
1190 /*--------------------------------------------------------------------*/
1191 /*
1192 -------------------------------------------------
1193 purpose -- to return a pointer to the location of
1194 matrix entry (irow,jcol) if present.
1195
1196 if entry (irow,jcol) is not present then
1197 *ppValue is NULL
1198 else entry (irow,jcol) is present then
1199 *ppValue is the location of the matrix entry
1200 endif
1201
1202 created -- 98may01, cca
1203 -------------------------------------------------
1204 */
1205 void
SubMtx_locationOfRealEntry(SubMtx * mtx,int irow,int jcol,double ** ppValue)1206 SubMtx_locationOfRealEntry (
1207 SubMtx *mtx,
1208 int irow,
1209 int jcol,
1210 double **ppValue
1211 ) {
1212 /*
1213 ---------------
1214 check the input
1215 ---------------
1216 */
1217 if ( mtx == NULL || irow < 0 || irow >= mtx->nrow || jcol < 0
1218 || jcol >= mtx->ncol || ppValue == NULL ) {
1219 fprintf(stderr,
1220 "\n fatal error in SubMtx_locationOfRealEntry(%p,%d,%d,%p)"
1221 "\n bad input\n", mtx, irow, jcol, ppValue) ;
1222 exit(-1) ;
1223 }
1224 if ( ! SUBMTX_IS_REAL(mtx) ) {
1225 fprintf(stderr,
1226 "\n fatal error in SubMtx_locationOfRealEntry(%p,%d,%d,%p)"
1227 "\n bad type %d, must be SPOOLES_REAL\n",
1228 mtx, irow, jcol, ppValue, mtx->type) ;
1229 exit(-1) ;
1230 }
1231 *ppValue = NULL ;
1232 switch ( mtx->mode ) {
1233 case SUBMTX_DENSE_ROWS :
1234 case SUBMTX_DENSE_COLUMNS : {
1235 double *entries ;
1236 int inc1, inc2, ncol, nrow, offset ;
1237
1238 SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
1239 if ( irow >= 0 && irow < nrow && jcol >= 0 && jcol < ncol ) {
1240 offset = irow*inc1 + jcol*inc2 ;
1241 *ppValue = entries + offset ;
1242 }
1243 } break ;
1244 case SUBMTX_SPARSE_ROWS : {
1245 double *entries ;
1246 int ii, jj, nent, nrow, offset, *indices, *sizes ;
1247
1248 SubMtx_sparseRowsInfo(mtx, &nrow, &nent, &sizes, &indices, &entries);
1249 if ( irow >= 0 && irow < nrow ) {
1250 for ( ii = offset = 0 ; ii < irow ; ii++ ) {
1251 offset += sizes[ii] ;
1252 }
1253 for ( ii = 0, jj = offset ; ii < sizes[irow] ; ii++, jj++ ) {
1254 if ( indices[jj] == jcol ) {
1255 *ppValue = entries + jj ;
1256 break ;
1257 }
1258 }
1259 }
1260 } break ;
1261 case SUBMTX_SPARSE_COLUMNS : {
1262 double *entries ;
1263 int ii, jj, nent, ncol, offset, *indices, *sizes ;
1264
1265 SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
1266 &sizes, &indices, &entries) ;
1267 if ( jcol >= 0 && jcol < ncol ) {
1268 for ( ii = offset = 0 ; ii < jcol ; ii++ ) {
1269 offset += sizes[ii] ;
1270 }
1271 for ( ii = 0, jj = offset ; ii < sizes[jcol] ; ii++, jj++ ) {
1272 if ( indices[jj] == irow ) {
1273 *ppValue = entries + jj ;
1274 break ;
1275 }
1276 }
1277 }
1278 } break ;
1279 case SUBMTX_SPARSE_TRIPLES : {
1280 double *entries ;
1281 int ii, nent, *colids, *rowids ;
1282
1283 SubMtx_sparseTriplesInfo(mtx, &nent, &rowids, &colids, &entries) ;
1284 for ( ii = 0 ; ii < nent ; ii++ ) {
1285 if ( irow == rowids[ii] && jcol == colids[ii] ) {
1286 *ppValue = entries + ii ;
1287 break ;
1288 }
1289 }
1290 } break ;
1291 case SUBMTX_DENSE_SUBROWS : {
1292 double *entries ;
1293 int ii, joff, nent, nrow, offset, *firstlocs, *sizes ;
1294
1295 SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
1296 &firstlocs, &sizes, &entries) ;
1297 if ( irow >= 0 && irow < nrow && sizes[irow] != 0 ) {
1298 for ( ii = offset = 0 ; ii < irow ; ii++ ) {
1299 offset += sizes[ii] ;
1300 }
1301 if ( 0 <= (joff = jcol - firstlocs[irow])
1302 && joff < sizes[irow] ) {
1303 offset += joff ;
1304 *ppValue = entries + offset ;
1305 break ;
1306 }
1307 }
1308 } break ;
1309 case SUBMTX_DENSE_SUBCOLUMNS : {
1310 double *entries ;
1311 int ii, ioff, nent, ncol, offset, *firstlocs, *sizes ;
1312
1313 SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
1314 &firstlocs, &sizes, &entries) ;
1315 if ( jcol >= 0 && jcol < ncol && sizes[jcol] != 0 ) {
1316 for ( ii = offset = 0 ; ii < jcol ; ii++ ) {
1317 offset += sizes[jcol] ;
1318 }
1319 if ( 0 <= (ioff = irow - firstlocs[jcol])
1320 && ioff < sizes[jcol] ) {
1321 offset += ioff ;
1322 *ppValue = entries + offset ;
1323 break ;
1324 }
1325 }
1326 } break ;
1327 case SUBMTX_DIAGONAL : {
1328 double *entries ;
1329 int ncol ;
1330
1331 if ( irow >= 0 && jcol >= 0 && irow == jcol ) {
1332 SubMtx_diagonalInfo(mtx, &ncol, &entries) ;
1333 if ( irow < ncol && jcol < ncol ) {
1334 *ppValue = entries + irow ;
1335 }
1336 }
1337 } break ;
1338 case SUBMTX_BLOCK_DIAGONAL_SYM :
1339 case SUBMTX_BLOCK_DIAGONAL_HERM : {
1340 double *entries ;
1341 int ii, ipivot, jrow, kk, m, ncol, nent, size ;
1342 int *pivotsizes ;
1343
1344 if ( irow >= 0 && jcol >= 0 ) {
1345 SubMtx_blockDiagonalInfo(mtx, &ncol, &nent,
1346 &pivotsizes, &entries) ;
1347 if ( irow < ncol && jcol < ncol ) {
1348 for ( jrow = ipivot = kk = 0 ; jrow <= irow ; ipivot++ ) {
1349 size = m = pivotsizes[ipivot] ;
1350 for ( ii = 0 ; ii < m ; ii++, jrow++ ) {
1351 if ( jrow == irow ) {
1352 if ( jrow - irow > m - ii ) {
1353 kk = -1 ;
1354 } else {
1355 kk += jrow - irow ;
1356 }
1357 } else {
1358 kk += size-- ;
1359 }
1360 }
1361 }
1362 if ( kk != -1 ) {
1363 *ppValue = entries + kk ;
1364 }
1365 }
1366 }
1367 } break ;
1368 default :
1369 fprintf(stderr,
1370 "\n fatal error in SubMtx_locationOfRealEntry(%p,%d,%d,%p)"
1371 "\n bad mode %d", mtx, irow, jcol, ppValue, mtx->mode) ;
1372 exit(-1) ;
1373 break ;
1374 }
1375 return ; }
1376
1377 /*--------------------------------------------------------------------*/
1378 /*
1379 --------------------------------------------------------
1380 purpose -- to return a pointer to the location of
1381 matrix entry (irow,jcol) if present.
1382
1383 if entry (irow,jcol) is not present then
1384 (*ppReal,*ppImag) is (NULL,NULL)
1385 else entry (irow,jcol) is present then
1386 (*ppReal,*ppImag) is the location of the matrix entry
1387 endif
1388
1389 created -- 98may01, cca
1390 --------------------------------------------------------
1391 */
1392 void
SubMtx_locationOfComplexEntry(SubMtx * mtx,int irow,int jcol,double ** ppReal,double ** ppImag)1393 SubMtx_locationOfComplexEntry (
1394 SubMtx *mtx,
1395 int irow,
1396 int jcol,
1397 double **ppReal,
1398 double **ppImag
1399 ) {
1400 /*
1401 ---------------
1402 check the input
1403 ---------------
1404 */
1405 if ( mtx == NULL || irow < 0 || irow >= mtx->nrow || jcol < 0
1406 || jcol >= mtx->ncol || ppReal == NULL || ppImag == NULL ) {
1407 fprintf(stderr,
1408 "\n fatal error in SubMtx_locationOfComplexEntry(%p,%d,%d,%p,%p)"
1409 "\n bad input\n", mtx, irow, jcol, ppReal, ppImag) ;
1410 exit(-1) ;
1411 }
1412 if ( ! SUBMTX_IS_COMPLEX(mtx) ) {
1413 fprintf(stderr,
1414 "\n fatal error in SubMtx_locationOfComplexEntry(%p,%d,%d,%p,%p)"
1415 "\n bad type %d, must be SPOOLES_COMPLEX\n",
1416 mtx, irow, jcol, ppReal, ppImag, mtx->type) ;
1417 exit(-1) ;
1418 }
1419 *ppReal = NULL ;
1420 *ppImag = NULL ;
1421 switch ( mtx->mode ) {
1422 case SUBMTX_DENSE_ROWS :
1423 case SUBMTX_DENSE_COLUMNS : {
1424 double *entries ;
1425 int inc1, inc2, ncol, nrow, offset ;
1426
1427 SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
1428 if ( irow >= 0 && irow < nrow && jcol >= 0 && jcol < ncol ) {
1429 offset = irow*inc1 + jcol*inc2 ;
1430 *ppReal = entries + 2*offset ;
1431 *ppImag = entries + 2*offset + 1 ;
1432 }
1433 } break ;
1434 case SUBMTX_SPARSE_ROWS : {
1435 double *entries ;
1436 int ii, jj, nent, nrow, offset, *indices, *sizes ;
1437
1438 SubMtx_sparseRowsInfo(mtx, &nrow, &nent, &sizes, &indices, &entries);
1439 if ( irow >= 0 && irow < nrow ) {
1440 for ( ii = offset = 0 ; ii < irow ; ii++ ) {
1441 offset += sizes[ii] ;
1442 }
1443 for ( ii = 0, jj = offset ; ii < sizes[irow] ; ii++, jj++ ) {
1444 if ( indices[jj] == jcol ) {
1445 *ppReal = entries + 2*jj ;
1446 *ppImag = entries + 2*jj + 1 ;
1447 break ;
1448 }
1449 }
1450 }
1451 } break ;
1452 case SUBMTX_SPARSE_COLUMNS : {
1453 double *entries ;
1454 int ii, jj, nent, ncol, offset, *indices, *sizes ;
1455
1456 SubMtx_sparseColumnsInfo(mtx, &ncol, &nent,
1457 &sizes, &indices, &entries) ;
1458 if ( jcol >= 0 && jcol < ncol ) {
1459 for ( ii = offset = 0 ; ii < jcol ; ii++ ) {
1460 offset += sizes[ii] ;
1461 }
1462 for ( ii = 0, jj = offset ; ii < sizes[jcol] ; ii++, jj++ ) {
1463 if ( indices[jj] == irow ) {
1464 *ppReal = entries + 2*jj ;
1465 *ppImag = entries + 2*jj + 1 ;
1466 break ;
1467 }
1468 }
1469 }
1470 } break ;
1471 case SUBMTX_SPARSE_TRIPLES : {
1472 double *entries ;
1473 int ii, nent, *colids, *rowids ;
1474
1475 SubMtx_sparseTriplesInfo(mtx, &nent, &rowids, &colids, &entries) ;
1476 for ( ii = 0 ; ii < nent ; ii++ ) {
1477 if ( irow == rowids[ii] && jcol == colids[ii] ) {
1478 *ppReal = entries + 2*ii ;
1479 *ppImag = entries + 2*ii + 1 ;
1480 break ;
1481 }
1482 }
1483 } break ;
1484 case SUBMTX_DENSE_SUBROWS : {
1485 double *entries ;
1486 int ii, joff, nent, nrow, offset, *firstlocs, *sizes ;
1487
1488 SubMtx_denseSubrowsInfo(mtx, &nrow, &nent,
1489 &firstlocs, &sizes, &entries) ;
1490 if ( irow >= 0 && irow < nrow && sizes[irow] != 0 ) {
1491 for ( ii = offset = 0 ; ii < irow ; ii++ ) {
1492 offset += sizes[ii] ;
1493 }
1494 if ( 0 <= (joff = jcol - firstlocs[irow])
1495 && joff < sizes[irow] ) {
1496 offset += joff ;
1497 *ppReal = entries + 2*offset ;
1498 *ppImag = entries + 2*offset + 1 ;
1499 break ;
1500 }
1501 }
1502 } break ;
1503 case SUBMTX_DENSE_SUBCOLUMNS : {
1504 double *entries ;
1505 int ii, ioff, nent, ncol, offset, *firstlocs, *sizes ;
1506
1507 SubMtx_denseSubcolumnsInfo(mtx, &ncol, &nent,
1508 &firstlocs, &sizes, &entries) ;
1509 if ( jcol >= 0 && jcol < ncol && sizes[jcol] != 0 ) {
1510 for ( ii = offset = 0 ; ii < jcol ; ii++ ) {
1511 offset += sizes[jcol] ;
1512 }
1513 if ( 0 <= (ioff = irow - firstlocs[jcol])
1514 && ioff < sizes[jcol] ) {
1515 offset += ioff ;
1516 *ppReal = entries + 2*offset ;
1517 *ppImag = entries + 2*offset + 1 ;
1518 break ;
1519 }
1520 }
1521 } break ;
1522 case SUBMTX_DIAGONAL : {
1523 double *entries ;
1524 int ncol ;
1525
1526 if ( irow >= 0 && jcol >= 0 && irow == jcol ) {
1527 SubMtx_diagonalInfo(mtx, &ncol, &entries) ;
1528 if ( irow < ncol && jcol < ncol ) {
1529 *ppReal = entries + 2*irow ;
1530 *ppImag = entries + 2*irow + 1 ;
1531 }
1532 }
1533 } break ;
1534 case SUBMTX_BLOCK_DIAGONAL_SYM :
1535 case SUBMTX_BLOCK_DIAGONAL_HERM : {
1536 double *entries ;
1537 int ii, ipivot, jrow, kk, m, ncol, nent, size ;
1538 int *pivotsizes ;
1539
1540 if ( irow >= 0 && jcol >= 0 ) {
1541 SubMtx_blockDiagonalInfo(mtx, &ncol, &nent,
1542 &pivotsizes, &entries) ;
1543 if ( irow < ncol && jcol < ncol ) {
1544 for ( jrow = ipivot = kk = 0 ; jrow <= irow ; ipivot++ ) {
1545 size = m = pivotsizes[ipivot] ;
1546 for ( ii = 0 ; ii < m ; ii++, jrow++ ) {
1547 if ( jrow == irow ) {
1548 if ( jrow - irow > m - ii ) {
1549 kk = -1 ;
1550 } else {
1551 kk += jrow - irow ;
1552 }
1553 } else {
1554 kk += size-- ;
1555 }
1556 }
1557 }
1558 if ( kk != -1 ) {
1559 *ppReal = entries + 2*kk ;
1560 *ppImag = entries + 2*kk + 1 ;
1561 }
1562 }
1563 }
1564 } break ;
1565 default :
1566 fprintf(stderr,
1567 "\n fatal error in SubMtx_locationOfComplexEntry(%p,%d,%d,%p,%p)"
1568 "\n bad mode %d", mtx, irow, jcol, ppReal, ppImag, mtx->mode) ;
1569 exit(-1) ;
1570 break ;
1571 }
1572 return ; }
1573
1574 /*--------------------------------------------------------------------*/
1575