1 /*
2 * ***************************************************************************
3 * MC = < Manifold Code >
4 * Copyright (C) 1994-- Michael Holst
5 *
6 * This library is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
10 *
11 * This library is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with this library; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 * rcsid="$Id: mat.c,v 1.45 2010/08/12 05:18:27 fetk Exp $"
21 * ***************************************************************************
22 */
23
24 /*
25 * ***************************************************************************
26 * File: mat.c
27 *
28 * Purpose: Class Mat: methods.
29 *
30 * Author: Michael Holst
31 * ***************************************************************************
32 */
33
34 #include "mat_p.h"
35
36 VEMBED(rcsid="$Id: mat.c,v 1.45 2010/08/12 05:18:27 fetk Exp $")
37
38 /*
39 * ***************************************************************************
40 * Class Mat: Inlineable methods
41 * ***************************************************************************
42 */
43 #if !defined(VINLINE_BAM)
44
45 #endif /* if !defined(VINLINE_BAM) */
46 /*
47 * ***************************************************************************
48 * Class Mat: Non-inlineable methods
49 * ***************************************************************************
50 */
51
52 /*
53 * ***************************************************************************
54 * Routine: Mat_ctor
55 *
56 * Purpose: The sparse matrix constructor.
57 *
58 * Notes: This constructor only fixes the number of rows and columns
59 * in the matrix; the nonzero structure is not set.
60 *
61 * Author: Michael Holst
62 * ***************************************************************************
63 */
Mat_ctor(Vmem * vmem,const char * name,int pnumR,int pnumC)64 VPUBLIC Mat* Mat_ctor(Vmem *vmem, const char *name, int pnumR, int pnumC)
65 {
66 Mat *thee = VNULL;
67
68 VDEBUGIO("Mat_ctor: CREATING object..");
69
70 /* initialize the memory manager and instantiate myself */
71 thee = Vmem_malloc( VNULL, 1, sizeof(Mat) );
72 if (vmem == VNULL) {
73 thee->vmem = Vmem_ctor( "Mat" );
74 thee->iMadeVmem = 1;
75 } else {
76 thee->vmem = vmem;
77 thee->iMadeVmem = 0;
78 }
79
80 thee->format = ZERO_FORMAT;
81 thee->state = NULL_STATE;
82 thee->sym = ISNOT_SYM;
83 thee->impl = ISNOT_IMPL;
84
85 thee->iMallocIJA = 0;
86 thee->iMallocA = 0;
87
88 /* initialize the parameters */
89 strncpy(thee->name, name, 10);
90 thee->numR = pnumR;
91 thee->numC = pnumC;
92 thee->numA = 0;
93 thee->numO = 0;
94 thee->numZ = 0;
95 thee->numBR = 0;
96 thee->numBC = 0;
97
98 /* initialize complex object pointers */
99 thee->slu = VNULL;
100 thee->lnkL = VNULL;
101 thee->lnkU = VNULL;
102 thee->xln = VNULL;
103 thee->xlnt = VNULL;
104
105 /* initialize primary malloc area pointers */
106 thee->IJA = VNULL;
107 thee->A = VNULL;
108 thee->BR = VNULL;
109 thee->BC = VNULL;
110
111 /* initialize alias pointers */
112 thee->IA = VNULL;
113 thee->JA = VNULL;
114 thee->diag = VNULL;
115 thee->offU = VNULL;
116 thee->offL = VNULL;
117
118 VDEBUGIO("..done.\n");
119
120 return thee;
121 }
122
123 /*
124 * ***************************************************************************
125 * Routine: Mat_dtor
126 *
127 * Purpose: The sparse matrix destructor.
128 *
129 * Notes: This destructor does the reverse of Mat_ctor, and if
130 * necessary first reverses Mat_initStructure
131 * (or Mat_copyStructure). I.e., if necessary,
132 * it first frees the large integer and real arrays created
133 * by Mat_initStructure or Mat_copyStructure, and then frees
134 * the Mat object itself at the last moment.
135 *
136 * Author: Michael Holst
137 * ***************************************************************************
138 */
Mat_dtor(Mat ** thee)139 VPUBLIC void Mat_dtor(Mat **thee)
140 {
141 /* VASSERT( (*thee) != VNULL ); */
142 if ((*thee) != VNULL) {
143
144 /* first free the large integer and real nonzero structure */
145 Mat_killStructure(*thee);
146
147 /* now kill the object itself */
148 VDEBUGIO("Mat_dtor: DESTROYING object..");
149 if ((*thee)->iMadeVmem) Vmem_dtor( &((*thee)->vmem) );
150 Vmem_free( VNULL, 1, sizeof(Mat), (void**)thee );
151 VDEBUGIO("..done.\n");
152
153 (*thee) = VNULL;
154 }
155 }
156
157 /*
158 * ***************************************************************************
159 * Routine: Mat_initStructure
160 *
161 * Purpose: Initialize the nonzero structure given structure information.
162 *
163 * Notes: This routine actually does the storage creation for both the
164 * integer structure information arrays and the nonzero value
165 * arrays.
166 *
167 * Author: Michael Holst
168 * ***************************************************************************
169 */
Mat_initStructure(Mat * thee,MATformat frmt,MATsym sym,int numO,int * IJA,double * A)170 VPUBLIC void Mat_initStructure(Mat *thee,
171 MATformat frmt, MATsym sym, int numO, int *IJA, double *A)
172 {
173 /* keep track of malloc actions */
174 thee->iMallocIJA = 0;
175 thee->iMallocA = 0;
176
177 /* build framework from input */
178 thee->format = frmt;
179 thee->sym = sym;
180 thee->numO = numO;
181
182 /* initialize the state */
183 thee->state = ZERO_STATE;
184
185 /* Format-dependent consistency checks and setup */
186 switch (thee->format) {
187
188 /* DRC consistency checks and setup */
189 case DRC_FORMAT:
190
191 /* DRC valid for symmetry OR nonsymmetry */
192 VASSERT( (thee->sym==ISNOT_SYM) || (thee->sym==IS_SYM) );
193
194
195 /* DRC only valid for square matrices */
196 VASSERT( thee->numR==thee->numC );
197
198 /* the setup */
199 thee->IJA = IJA; /* integer structure via hand-off (don't malloc) */
200 thee->IA = thee->IJA;
201 thee->JA = thee->IJA + thee->numR + 1;
202 if (thee->sym == ISNOT_SYM) {
203 thee->numZ = thee->numR + 2*thee->numO;
204 thee->numA = thee->numR + 2*thee->numO;
205 if (thee->numA > 0) {
206 VASSERT( A == VNULL );
207 thee->A = Vmem_malloc(thee->vmem,thee->numA,sizeof(double));
208 thee->iMallocA = 1;
209 }
210 thee->diag = thee->A;
211 thee->offU = thee->A + thee->numR;
212 thee->offL = thee->A + thee->numR + thee->numO;
213
214 } else { /* (thee->sym == IS_SYM) */
215 thee->numZ = thee->numR + 2*thee->numO;
216 thee->numA = thee->numR + thee->numO;
217 if (thee->numA > 0) {
218 VASSERT( A == VNULL );
219 thee->A = Vmem_malloc(thee->vmem,thee->numA,sizeof(double));
220 thee->iMallocA = 1;
221 }
222 thee->diag = thee->A;
223 thee->offU = thee->A + thee->numR;
224 thee->offL = thee->offU;
225 }
226 break;
227
228 /* ROW consistency checks and setup */
229 case ROW_FORMAT:
230
231 /* ROW only valid for nonsymmetry */
232 VASSERT( thee->sym==ISNOT_SYM );
233
234 /* the setup */
235 thee->IJA = IJA; /* integer structure via hand-off (don't malloc) */
236 thee->IA = thee->IJA;
237 thee->JA = thee->IJA + thee->numR + 1;
238 thee->numZ = thee->numO;
239 thee->numA = thee->numO;
240 if ( thee->numA > 0 ) {
241 VASSERT( A == VNULL );
242 thee->A = Vmem_malloc(thee->vmem,thee->numA,sizeof(double));
243 thee->iMallocA = 1;
244 }
245 thee->offU = thee->A;
246 thee->offL = VNULL;
247 thee->diag = VNULL;
248 break;
249
250 /* COL consistency checks */
251 case COL_FORMAT:
252
253 /* COL only valid for nonsymmetry */
254 VASSERT( thee->sym==ISNOT_SYM );
255
256 /* the setup */
257 thee->IJA = IJA; /* integer structure via hand-off (don't malloc) */
258 thee->IA = thee->IJA;
259 thee->JA = thee->IJA + thee->numC + 1;
260 thee->numZ = thee->numO;
261 thee->numA = thee->numO;
262 if ( thee->numA > 0 ) {
263 VASSERT( A == VNULL );
264 thee->A = Vmem_malloc(thee->vmem,thee->numA,sizeof(double));
265 thee->iMallocA = 1;
266 }
267 thee->offL = thee->A;
268 thee->offU = VNULL;
269 thee->diag = VNULL;
270 break;
271
272 /* SLU consistency checks */
273 case SLU_FORMAT:
274
275 /* SLU valid for nonsymmetry */
276 VASSERT( thee->sym==ISNOT_SYM );
277
278 /* the setup (we MUST malloc the integer structure) */
279 thee->IJA = Vmem_malloc( thee->vmem, thee->numR+1+numO, sizeof(int) );
280 thee->iMallocIJA = 1;
281 thee->IA = thee->IJA;
282 thee->JA = thee->IJA + thee->numC + 1;
283 thee->numZ = thee->numO;
284 thee->numA = thee->numO;
285 if ( thee->numA > 0 ) {
286 VASSERT( A == VNULL );
287 thee->A = Vmem_malloc(thee->vmem,thee->numA,sizeof(double));
288 thee->iMallocA = 1;
289 }
290 thee->offL = VNULL;
291 thee->offU = VNULL;
292 thee->diag = VNULL;
293 break;
294
295 /* RLN consistency checks and setup */
296 case RLN_FORMAT:
297
298 Mat_initStructureLN( thee, frmt, sym );
299 break;
300
301 /* CLN consistency checks and setup */
302 case CLN_FORMAT:
303
304 Mat_initStructureLN( thee, frmt, sym );
305 break;
306
307 /* XLN consistency checks and setup */
308 case XLN_FORMAT:
309
310 Mat_initStructureLN( thee, frmt, sym );
311 break;
312
313 /* RFL consistency checks and setup */
314 case RFL_FORMAT:
315
316 /* RFL only valid for nonsymmetry */
317 VASSERT( thee->sym==ISNOT_SYM );
318
319 /* the setup (input numO and IJA are ignored) */
320 thee->numO = thee->numR * thee->numC;
321 thee->numA = thee->numO;
322 thee->numZ = thee->numO;
323 if ( thee->numA > 0 ) {
324 if (A != VNULL) {
325 thee->A = A;
326 } else {
327 VASSERT( A == VNULL );
328 thee->A = Vmem_malloc(thee->vmem,thee->numA,sizeof(double));
329 thee->iMallocA = 1;
330 }
331 }
332 break;
333
334 /* CFL consistency checks and setup */
335 case CFL_FORMAT:
336
337 /* CFL only valid for nonsymmetry */
338 VASSERT( thee->sym==ISNOT_SYM );
339
340 /* the setup (input numO and IJA are ignored) */
341 thee->numO = thee->numR * thee->numC;
342 thee->numA = thee->numO;
343 thee->numZ = thee->numO;
344 if ( thee->numA > 0 ) {
345 if (A != VNULL) {
346 thee->A = A;
347 } else {
348 VASSERT( A == VNULL );
349 thee->A = Vmem_malloc(thee->vmem,thee->numA,sizeof(double));
350 thee->iMallocA = 1;
351 }
352 }
353 break;
354
355 default:
356 VASSERT(0);
357 break;
358 }
359 }
360
361 /*
362 * ***************************************************************************
363 * Routine: Mat_copyStructure
364 *
365 * Purpose: Initialize the nonzero structure given structure information.
366 *
367 * Notes: This routine actually does the storage creation for both the
368 * integer structure information arrays and the nonzero value
369 * arrays.
370 *
371 * Author: Michael Holst
372 * ***************************************************************************
373 */
Mat_copyStructure(Mat * thee,Mat * model)374 VPUBLIC void Mat_copyStructure(Mat *thee, Mat *model)
375 {
376 int i, size;
377
378 VASSERT( thee->numR == model->numR );
379 VASSERT( thee->numC == model->numC );
380
381 /* initialize the state */
382 thee->state = ZERO_STATE;
383
384 /* build framework (numR and numC have already been set in ctor) */
385 thee->format = model->format; /* same format */
386 thee->sym = model->sym; /* keep symmetry */
387 thee->impl = model->impl; /* keep implicitness */
388 thee->numO = model->numO; /* numA and numZ are set below */
389
390 /* keep track of malloc actions */
391 thee->iMallocIJA = 0;
392 thee->iMallocA = 0;
393
394 /* nonzero integer structure is via copy (we malloc) */
395 size = thee->numR + 1 + thee->numO;
396 thee->IJA = Vmem_malloc(thee->vmem,size,sizeof(int));
397 thee->iMallocIJA = 1;
398 for (i=0; i<size; i++) {
399 thee->IJA[i] = model->IJA[i];
400 }
401
402 /* boundary structure is via copy (we malloc) */
403 thee->numBR = model->numBR;
404 if (thee->numBR > 0) {
405 thee->BR = Vmem_malloc(thee->vmem,thee->numBR,sizeof(int));
406 for (i=0; i<thee->numBR; i++) {
407 thee->BR[i] = model->BR[i];
408 }
409 }
410 thee->numBC = model->numBC;
411 if (thee->numBC > 0) {
412 thee->BC = Vmem_malloc(thee->vmem,thee->numBC,sizeof(int));
413 for (i=0; i<thee->numBC; i++) {
414 thee->BC[i] = model->BC[i];
415 }
416 }
417
418 /* Format-dependent consistency checks and setup */
419 switch (thee->format) {
420
421 /* DRC consistency checks and setup */
422 case DRC_FORMAT:
423
424 /* DRC valid for symmetry OR nonsymmetry */
425 VASSERT( (thee->sym==ISNOT_SYM) || (thee->sym==IS_SYM) );
426
427 /* DRC only valid for square matrices */
428 VASSERT( thee->numR==thee->numC );
429
430 /* the setup */
431 thee->IA = thee->IJA;
432 thee->JA = thee->IJA + thee->numR + 1;
433 if (thee->sym == ISNOT_SYM) {
434 thee->numZ = thee->numR + 2*thee->numO;
435 thee->numA = thee->numR + 2*thee->numO;
436 if (thee->numA > 0) {
437 thee->A = Vmem_malloc(thee->vmem,thee->numA,sizeof(double));
438 thee->iMallocA = 1;
439 }
440 thee->diag = thee->A;
441 thee->offU = thee->A + thee->numR;
442 thee->offL = thee->A + thee->numR + thee->numO;
443
444 } else { /* (thee->sym == IS_SYM) */
445 thee->numZ = thee->numR + 2*thee->numO;
446 thee->numA = thee->numR + thee->numO;
447 if (thee->numA > 0) {
448 thee->A = Vmem_malloc(thee->vmem,thee->numA,sizeof(double));
449 thee->iMallocA = 1;
450 }
451 thee->diag = thee->A;
452 thee->offU = thee->A + thee->numR;
453 thee->offL = thee->offU;
454 }
455 break;
456
457 /* ROW consistency checks and setup */
458 case ROW_FORMAT:
459
460 /* ROW only valid for nonsymmetry */
461 VASSERT( thee->sym==ISNOT_SYM );
462
463 /* the setup */
464 thee->IA = thee->IJA;
465 thee->JA = thee->IJA + thee->numR + 1;
466 thee->numZ = thee->numO;
467 thee->numA = thee->numO;
468 if ( thee->numA > 0 ) {
469 thee->A = Vmem_malloc(thee->vmem,thee->numA,sizeof(double));
470 thee->iMallocA = 1;
471 }
472 break;
473
474 /* COL consistency checks */
475 case COL_FORMAT:
476
477 /* COL only valid for nonsymmetry */
478 VASSERT( thee->sym==ISNOT_SYM );
479
480 /* the setup */
481 thee->IA = thee->IJA;
482 thee->JA = thee->IJA + thee->numC + 1;
483 thee->numZ = thee->numO;
484 thee->numA = thee->numO;
485 if ( thee->numA > 0 ) {
486 thee->A = Vmem_malloc(thee->vmem,thee->numA,sizeof(double));
487 thee->iMallocA = 1;
488 }
489 break;
490
491 /* RFL consistency checks and setup */
492 case RFL_FORMAT:
493
494 /* RFL only valid for nonsymmetry */
495 VASSERT( thee->sym==ISNOT_SYM );
496
497 /* the setup (input numO is already copied over) */
498 thee->numA = thee->numO;
499 thee->numZ = thee->numO;
500 if ( thee->numA > 0 ) {
501 thee->A = Vmem_malloc(thee->vmem,thee->numA,sizeof(double));
502 thee->iMallocA = 1;
503 }
504 break;
505
506 /* CFL consistency checks and setup */
507 case CFL_FORMAT:
508
509 /* CFL only valid for nonsymmetry */
510 VASSERT( thee->sym==ISNOT_SYM );
511
512 /* the setup (input numO is already copied over) */
513 thee->numA = thee->numO;
514 thee->numZ = thee->numO;
515 if ( thee->numA > 0 ) {
516 thee->A = Vmem_malloc(thee->vmem,thee->numA,sizeof(double));
517 thee->iMallocA = 1;
518 }
519 break;
520
521 default:
522 VASSERT(0);
523 break;
524 }
525 }
526
527 /*
528 * ***************************************************************************
529 * Routine: Mat_killStructure
530 *
531 * Purpose: Kill the nonzero structure and structure information.
532 *
533 * Notes: This routine does the reverse of Mat_initStructure
534 * (or Mat_copyStructure). It leaves only the information
535 * about the number of blocks, number of rows, and number of
536 * columns per block. I.e., what is left is only what was
537 * present after the initial call to Mat_ctor.
538 *
539 * Author: Michael Holst
540 * ***************************************************************************
541 */
Mat_killStructure(Mat * thee)542 VPUBLIC void Mat_killStructure(Mat *thee)
543 {
544 int size;
545
546 /* VASSERT( (*thee) != VNULL ); */
547 if (thee != VNULL) {
548 if (thee->state != NULL_STATE) {
549
550 /* Format-dependent destruction */
551 switch (thee->format) {
552
553 case DRC_FORMAT:
554
555 /* DRC only valid for symmetry OR nonsymmetry */
556 VASSERT( (thee->sym==ISNOT_SYM)||(thee->sym==IS_SYM) );
557 /* DRC only valid for square matrices */
558 VASSERT( thee->numR==thee->numC );
559
560 /* release of nonzero structure */
561 size = thee->numR + 1 + thee->numO;
562 /* VASSERT( thee->iMallocIJA == 1 ); */
563 Vmem_free( thee->vmem, size, sizeof(int),
564 (void**)&(thee->IJA) );
565
566 /* release of boundary row/column arrays */
567 if ( thee->BR != VNULL ) {
568 Vmem_free( thee->vmem, thee->numBR, sizeof(int),
569 (void**)&(thee->BR) );
570 }
571 if ( thee->BC != VNULL ) {
572 Vmem_free( thee->vmem, thee->numBC, sizeof(int),
573 (void**)&(thee->BC) );
574 }
575 thee->numBR = 0;
576 thee->numBC = 0;
577
578 /* release of nonzeros */
579 if ( thee->numA > 0 ) {
580 /* If VWARN is triggered, it will either be memory leak */
581 /* or bus error */
582 VWARN( thee->iMallocA == 1 );
583 Vmem_free( thee->vmem, thee->numA, sizeof(double),
584 (void**)&(thee->A) );
585 thee->numA = 0;
586 thee->numO = 0;
587 thee->numZ = 0;
588 }
589
590 break;
591
592 case ROW_FORMAT:
593
594 /* ROW only valid for nonsymmetry */
595 VASSERT( (thee->sym==ISNOT_SYM) );
596
597 /* release of nonzero structure */
598 size = thee->numR + 1 + thee->numO;
599 /* VASSERT( thee->iMallocIJA == 1 ); */
600 Vmem_free( thee->vmem, size, sizeof(int),
601 (void**)&(thee->IJA) );
602
603 /* release of boundary row/column arrays */
604 if ( thee->BR != VNULL ) {
605 Vmem_free( thee->vmem, thee->numBR, sizeof(int),
606 (void**)&(thee->BR) );
607 }
608 if ( thee->BC != VNULL ) {
609 Vmem_free( thee->vmem, thee->numBC, sizeof(int),
610 (void**)&(thee->BC) );
611 }
612 thee->numBR = 0;
613 thee->numBC = 0;
614
615 /* release of nonzeros */
616 if ( thee->numA > 0 ) {
617 VASSERT( thee->iMallocA == 1 );
618 Vmem_free( thee->vmem, thee->numA, sizeof(double),
619 (void**)&(thee->A) );
620 thee->numA = 0;
621 thee->numO = 0;
622 thee->numZ = 0;
623 }
624
625 break;
626
627 case COL_FORMAT:
628
629 /* COL only valid for nonsymmetry */
630 VASSERT( (thee->sym==ISNOT_SYM) );
631
632 /* release of nonzero structure */
633 size = thee->numC + 1 + thee->numO;
634 /* VASSERT( thee->iMallocIJA == 1 ); */
635 Vmem_free( thee->vmem, size, sizeof(int),
636 (void**)&(thee->IJA) );
637
638 /* release of boundary row/column arrays */
639 if ( thee->BR != VNULL ) {
640 Vmem_free( thee->vmem, thee->numBR, sizeof(int),
641 (void**)&(thee->BR) );
642 }
643 if ( thee->BC != VNULL ) {
644 Vmem_free( thee->vmem, thee->numBC, sizeof(int),
645 (void**)&(thee->BC) );
646 }
647 thee->numBR = 0;
648 thee->numBC = 0;
649
650 /* release of nonzeros */
651 if ( thee->numA > 0 ) {
652 VASSERT( thee->iMallocA == 1 );
653 Vmem_free( thee->vmem, thee->numA, sizeof(double),
654 (void**)&(thee->A) );
655 thee->numA = 0;
656 thee->numO = 0;
657 thee->numZ = 0;
658 }
659
660 break;
661
662 case SLU_FORMAT:
663
664 /* SLU only valid for nonsymmetry */
665 VASSERT( (thee->sym==ISNOT_SYM) );
666
667 /* release of nonzero structure */
668 size = thee->numC + 1 + thee->numO;
669 /* VASSERT( thee->iMallocIJA == 1 ); */
670 Vmem_free( thee->vmem, size, sizeof(int),
671 (void**)&(thee->IJA) );
672
673 /* release of boundary row/column arrays */
674 if ( thee->BR != VNULL ) {
675 Vmem_free( thee->vmem, thee->numBR, sizeof(int),
676 (void**)&(thee->BR) );
677 }
678 if ( thee->BC != VNULL ) {
679 Vmem_free( thee->vmem, thee->numBC, sizeof(int),
680 (void**)&(thee->BC) );
681 }
682 thee->numBR = 0;
683 thee->numBC = 0;
684
685 /* release of nonzeros */
686 if ( thee->numA > 0 ) {
687 VASSERT( thee->iMallocA == 1 );
688 Vmem_free( thee->vmem, thee->numA, sizeof(double),
689 (void**)&(thee->A) );
690 thee->numA = 0;
691 thee->numO = 0;
692 thee->numZ = 0;
693 }
694
695 /* finally destroy the internal objects */
696 Mat_sluDestroy(thee);
697
698 break;
699
700 case RLN_FORMAT:
701 Mat_killStructureLN( thee );
702 break;
703
704 case CLN_FORMAT:
705 Mat_killStructureLN( thee );
706 break;
707
708 case XLN_FORMAT:
709 Mat_killStructureLN( thee );
710 break;
711
712 case RFL_FORMAT:
713
714 /* RFL only valid for nonsymmetry */
715 VASSERT( (thee->sym==ISNOT_SYM) );
716
717 /* release of nonzeros */
718 if ( thee->numA > 0 ) {
719 if ( thee->iMallocA == 1 ) {
720 Vmem_free( thee->vmem, thee->numA, sizeof(double),
721 (void**)&(thee->A) );
722 }
723 thee->numA = 0;
724 thee->numO = 0;
725 thee->numZ = 0;
726 }
727 break;
728
729 case CFL_FORMAT:
730
731 /* CFL only valid for nonsymmetry */
732 VASSERT( (thee->sym==ISNOT_SYM) );
733
734 /* release of nonzeros */
735 if ( thee->numA > 0 ) {
736 if ( thee->iMallocA == 1 ) {
737 Vmem_free( thee->vmem, thee->numA, sizeof(double),
738 (void**)&(thee->A) );
739 }
740 thee->numA = 0;
741 thee->numO = 0;
742 thee->numZ = 0;
743 }
744 break;
745
746 default:
747 VASSERT(0);
748 break;
749 }
750
751 /* initialize the state */
752 thee->state = NULL_STATE;
753
754 /* keep track of malloc actions */
755 thee->iMallocIJA = 0;
756 thee->iMallocA = 0;
757 }
758 }
759 }
760
761 /*
762 * ***************************************************************************
763 * Routine: Mat_numR
764 *
765 * Purpose: Return number of row in the matrix.
766 *
767 * Author: Michael Holst
768 * ***************************************************************************
769 */
Mat_numR(Mat * thee)770 VPUBLIC int Mat_numR(Mat *thee)
771 {
772 return thee->numR;
773 }
774
775 /*
776 * ***************************************************************************
777 * Routine: Mat_numC
778 *
779 * Purpose: Return number of columns in the matrix.
780 *
781 * Author: Michael Holst
782 * ***************************************************************************
783 */
Mat_numC(Mat * thee)784 VPUBLIC int Mat_numC(Mat *thee)
785 {
786 return thee->numC;
787 }
788
789 /*
790 * ***************************************************************************
791 * Routine: Mat_numA
792 *
793 * Purpose: Return number of nonzeros we are actually storing.
794 *
795 * Author: Michael Holst
796 * ***************************************************************************
797 */
Mat_numA(Mat * thee)798 VPUBLIC int Mat_numA(Mat *thee)
799 {
800 return thee->numA;
801 }
802
803 /*
804 * ***************************************************************************
805 * Routine: Mat_numO
806 *
807 * Purpose: Return number of nonzeros we are actually storing
808 * which are located in upper (or lower) triangle.
809 *
810 * Author: Michael Holst
811 * ***************************************************************************
812 */
Mat_numO(Mat * thee)813 VPUBLIC int Mat_numO(Mat *thee)
814 {
815 return thee->numO;
816 }
817
818 /*
819 * ***************************************************************************
820 * Routine: Mat_numZ
821 *
822 * Purpose: Return number of nonzeros we WOULD be storing if we were
823 * ignoring symmetry and storing all nonzeros.
824 *
825 * Author: Michael Holst
826 * ***************************************************************************
827 */
Mat_numZ(Mat * thee)828 VPUBLIC int Mat_numZ(Mat *thee)
829 {
830 return thee->numZ;
831 }
832
833 /*
834 * ***************************************************************************
835 * Routine: Mat_format
836 *
837 * Purpose: Return the format.
838 *
839 * Author: Michael Holst
840 * ***************************************************************************
841 */
Mat_format(Mat * thee)842 VPUBLIC MATformat Mat_format(Mat *thee)
843 {
844 return thee->format;
845 }
846
847 /*
848 * ***************************************************************************
849 * Routine: Mat_sym
850 *
851 * Purpose: Return the symmetry.
852 *
853 * Author: Michael Holst
854 * ***************************************************************************
855 */
Mat_sym(Mat * thee)856 VPUBLIC MATsym Mat_sym(Mat *thee)
857 {
858 return thee->sym;
859 }
860
861 /*
862 * ***************************************************************************
863 * Routine: Mat_state
864 *
865 * Purpose: Return the state.
866 *
867 * Author: Michael Holst
868 * ***************************************************************************
869 */
Mat_state(Mat * thee)870 VPUBLIC MATstate Mat_state(Mat *thee)
871 {
872 return thee->state;
873 }
874
875 /*
876 * ***************************************************************************
877 * Routine: Mat_impl
878 *
879 * Purpose: Return the impl.
880 *
881 * Author: Michael Holst
882 * ***************************************************************************
883 */
Mat_impl(Mat * thee)884 VPUBLIC MATimpl Mat_impl(Mat *thee)
885 {
886 return thee->impl;
887 }
888
889 /*
890 * ***************************************************************************
891 * Routine: Mat_setFormat
892 *
893 * Purpose: Set the format.
894 *
895 * Author: Michael Holst
896 * ***************************************************************************
897 */
Mat_setFormat(Mat * thee,MATformat format)898 VPUBLIC void Mat_setFormat(Mat *thee, MATformat format)
899 {
900 thee->format = format;
901 }
902
903 /*
904 * ***************************************************************************
905 * Routine: Mat_setSym
906 *
907 * Purpose: Set the symmetry.
908 *
909 * Author: Michael Holst
910 * ***************************************************************************
911 */
Mat_setSym(Mat * thee,MATsym sym)912 VPUBLIC void Mat_setSym(Mat *thee, MATsym sym)
913 {
914 thee->sym = sym;
915 }
916
917 /*
918 * ***************************************************************************
919 * Routine: Mat_setState
920 *
921 * Purpose: Set the state.
922 *
923 * Author: Michael Holst
924 * ***************************************************************************
925 */
Mat_setState(Mat * thee,MATstate state)926 VPUBLIC void Mat_setState(Mat *thee, MATstate state)
927 {
928 thee->state = state;
929 }
930
931 /*
932 * ***************************************************************************
933 * Routine: Mat_setImpl
934 *
935 * Purpose: Set the impl.
936 *
937 * Author: Michael Holst
938 * ***************************************************************************
939 */
Mat_setImpl(Mat * thee,MATimpl impl)940 VPUBLIC void Mat_setImpl(Mat *thee, MATimpl impl)
941 {
942 thee->impl = impl;
943 }
944
945 /*
946 * ***************************************************************************
947 * Routine: Mat_sizeIJA
948 *
949 * Purpose: Return total number of INTEGER STORAGE LOCATIONS in the matrix.
950 *
951 * Author: Michael Holst
952 * ***************************************************************************
953 */
Mat_sizeIJA(Mat * thee)954 VPUBLIC int Mat_sizeIJA(Mat *thee)
955 {
956 return (thee->numR + 1 + thee->numO);
957 }
958
959 /*
960 * ***************************************************************************
961 * Routine: Mat_sizeA
962 *
963 * Purpose: Return total number of REAL STORAGE LOCATIONS in the matrix.
964 *
965 * Author: Michael Holst
966 * ***************************************************************************
967 */
Mat_sizeA(Mat * thee)968 VPUBLIC int Mat_sizeA(Mat *thee)
969 {
970 return Mat_numA(thee);
971 }
972
973 /*
974 * ***************************************************************************
975 * Routine: Mat_IJA
976 *
977 * Purpose: Return the integer structure IJA.
978 *
979 * Author: Michael Holst
980 * ***************************************************************************
981 */
Mat_IJA(Mat * thee)982 VPUBLIC int *Mat_IJA(Mat *thee)
983 {
984 return thee->IJA;
985 }
986
987 /*
988 * ***************************************************************************
989 * Routine: Mat_IA
990 *
991 * Purpose: Return the integer structure IA.
992 *
993 * Author: Michael Holst
994 * ***************************************************************************
995 */
Mat_IA(Mat * thee)996 VPUBLIC int *Mat_IA(Mat *thee)
997 {
998 return thee->IA;
999 }
1000
1001 /*
1002 * ***************************************************************************
1003 * Routine: Mat_JA
1004 *
1005 * Purpose: Return the integer structure JA.
1006 *
1007 * Author: Michael Holst
1008 * ***************************************************************************
1009 */
Mat_JA(Mat * thee)1010 VPUBLIC int *Mat_JA(Mat *thee)
1011 {
1012 return thee->JA;
1013 }
1014
1015 /*
1016 * ***************************************************************************
1017 * Routine: Mat_A
1018 *
1019 * Purpose: Return the real structure A.
1020 *
1021 * Author: Michael Holst
1022 * ***************************************************************************
1023 */
Mat_A(Mat * thee)1024 VPUBLIC double *Mat_A(Mat *thee)
1025 {
1026 return thee->A;
1027 }
1028
1029 /*
1030 * ***************************************************************************
1031 * Routine: Mat_diag
1032 *
1033 * Purpose: Return the diagonal structure A.
1034 *
1035 * Author: Michael Holst
1036 * ***************************************************************************
1037 */
Mat_diag(Mat * thee)1038 VPUBLIC double *Mat_diag(Mat *thee)
1039 {
1040 return thee->diag;
1041 }
1042
1043 /*
1044 * ***************************************************************************
1045 * Routine: Mat_offU
1046 *
1047 * Purpose: Return the upper-triangle structure A.
1048 *
1049 * Author: Michael Holst
1050 * ***************************************************************************
1051 */
Mat_offU(Mat * thee)1052 VPUBLIC double *Mat_offU(Mat *thee)
1053 {
1054 return thee->offU;
1055 }
1056
1057 /*
1058 * ***************************************************************************
1059 * Routine: Mat_offL
1060 *
1061 * Purpose: Return the lower-triangle structure A.
1062 *
1063 * Author: Michael Holst
1064 * ***************************************************************************
1065 */
Mat_offL(Mat * thee)1066 VPUBLIC double *Mat_offL(Mat *thee)
1067 {
1068 return thee->offL;
1069 }
1070
1071 /*
1072 * ***************************************************************************
1073 * Routine: Mat_print
1074 *
1075 * Purpose: Print the matrix as a DENSE matrix in MATLAB format.
1076 *
1077 * Author: Michael Holst
1078 * ***************************************************************************
1079 */
Mat_print(Mat * thee)1080 VPUBLIC void Mat_print(Mat *thee)
1081 {
1082 int i, j, k;
1083 int numR, numC;
1084 char rn[80];
1085 const int MaxRows = 100;
1086 const int MaxCols = 30;
1087 double matrix[100][30];
1088
1089 if ( ( thee->format == RLN_FORMAT ) ||
1090 ( thee->format == CLN_FORMAT ) ||
1091 ( thee->format == XLN_FORMAT ) ) {
1092 Mat_printLN( thee );
1093 return;
1094 }
1095
1096 numR = thee->numR;
1097 numC = thee->numC;
1098
1099 strncpy(rn,"Mat_print:",80);
1100
1101 /* some i/o */
1102 Vnm_print(0, "%s printing <%s>" " [dim=(%dx%d),sym=%d,numA=%d,numO=%d]\n",
1103 rn, thee->name, thee->numR, thee->numC, thee->sym,
1104 thee->numA, thee->numO);
1105
1106 /* size check */
1107 if ((numR > MaxRows) || (numC > MaxCols)) {
1108 Vnm_print(0, "%smatrix too large to view....skipping.\n", rn);
1109 return;
1110 }
1111
1112 /* make a dense matrix first */
1113 for (i=0; i<numR; i++) {
1114 for (j=0; j<numC; j++) {
1115 matrix[i][j] = 0.0;
1116 }
1117 }
1118
1119 if (thee->state != NULL_STATE) {
1120
1121 switch (thee->format) {
1122
1123 case DRC_FORMAT:
1124
1125 for (i=0; i<numR; i++) {
1126 matrix[i][i] = thee->diag[i];
1127 for (k=thee->IA[i]; k<thee->IA[i+1]; k++) {
1128 j = thee->JA[k];
1129 matrix[i][j] = thee->offU[k];
1130 matrix[j][i] = thee->offL[k];
1131 }
1132 }
1133 break;
1134
1135 case ROW_FORMAT:
1136
1137 for (i=0; i<numR; i++) {
1138 for (k=thee->IA[i]; k<thee->IA[i+1]; k++) {
1139 j = thee->JA[k];
1140 matrix[i][j] = thee->A[k];
1141 }
1142 }
1143 break;
1144
1145 case COL_FORMAT:
1146
1147 for (j=0; j<numC; j++) {
1148 for (k=thee->IA[j]; k<thee->IA[j+1]; k++) {
1149 i = thee->JA[k];
1150 matrix[i][j] = thee->A[k];
1151 }
1152 }
1153 break;
1154
1155 case RFL_FORMAT:
1156
1157 for (i=0; i<numR; i++) {
1158 for (j=0; j<numC; j++) {
1159 matrix[i][j] = thee->A[ (i*numC) + j ];
1160 }
1161 }
1162 break;
1163
1164 case CFL_FORMAT:
1165
1166 for (i=0; i<numR; i++) {
1167 for (j=0; j<numC; j++) {
1168 matrix[i][j] = thee->A[ (j*numR) + i ];
1169 }
1170 }
1171 break;
1172
1173 default:
1174 Vnm_print(0,
1175 "%smatrix not in correct format to print....skipping.\n", rn);
1176 break;
1177 }
1178 }
1179
1180 /* print the matrix */
1181 Vnm_print(0, "%s = [\n", thee->name);
1182 for (i=0; i<numR; i++) {
1183 for (j=0; j<numC; j++) {
1184 if (VABS(matrix[i][j]) < 0.0001) {
1185 Vnm_print(0, " 0.0 ");
1186 } else {
1187 Vnm_print(0, "%7.3f", matrix[i][j]);
1188 }
1189 }
1190 Vnm_print(0, "\n");
1191 }
1192 Vnm_print(0, "];\n");
1193 }
1194
1195 /*
1196 * ***************************************************************************
1197 * Routine: Mat_printSp
1198 *
1199 * Purpose: Print the matrix as a SPARSE matrix in MATLAB format.
1200 *
1201 * Notes: flag==0 ==> write
1202 * flag==1 ==> append
1203 *
1204 * Author: Michael Holst
1205 * ***************************************************************************
1206 */
Mat_printSp(Mat * thee,char * fname,int pflag)1207 VPUBLIC void Mat_printSp(Mat *thee, char *fname, int pflag)
1208 {
1209 int i, j, k;
1210 int numR, numC;
1211 char rn[80];
1212 FILE *fp;
1213
1214 if ( ( thee->format == RLN_FORMAT ) ||
1215 ( thee->format == CLN_FORMAT ) ||
1216 ( thee->format == XLN_FORMAT ) ) {
1217 Mat_printLNSp( thee, fname, pflag );
1218 return;
1219 }
1220
1221 numR = thee->numR;
1222 numC = thee->numC;
1223
1224 strncpy(rn,"Mat_printSp:",80);
1225
1226 if (pflag == 0) {
1227 fp=fopen(fname,"w");
1228 } else if (pflag == 1) {
1229 fp=fopen(fname,"a");
1230 } else {
1231 fp=VNULL;
1232 }
1233
1234 if (fp==VNULL) {
1235 Vnm_print(2, "%s problem opening file <%s>\n", rn, fname);
1236 return;
1237 }
1238
1239 Vnm_print(0, "%s printing <%s>" " [dim=(%dx%d),sym=%d,numA=%d,numO=%d]"
1240 " to file <%s>\n",
1241 rn, thee->name, thee->numR, thee->numC, thee->sym,
1242 thee->numA, thee->numO, fname);
1243
1244 /* print the matrix in matlab sparse format */
1245 fprintf(fp,"%% %s matrix <%s> [dim=(%dx%d)]\n",
1246 rn, thee->name, numR, numC);
1247
1248 if (thee->state != NULL_STATE) {
1249
1250 fprintf(fp,"%% ----------------------------------------\n");
1251 fprintf(fp, "fprintf('Defining %s_IJA..');", thee->name);
1252 fprintf(fp, "\n%s_IJA = [\n", thee->name);
1253
1254 switch (thee->format) {
1255
1256 case DRC_FORMAT:
1257 for (i=0; i<numR; i++) {
1258 fprintf(fp, " %d %d %15.8e\n",
1259 i+1, i+1, thee->diag[i]);
1260 for (k=thee->IA[i]; k<thee->IA[i+1]; k++) {
1261 j = thee->JA[k];
1262 fprintf(fp, " %d %d %15.8e\n",
1263 i+1, j+1, thee->offU[k]);
1264 fprintf(fp, " %d %d %15.8e\n",
1265 j+1, i+1, thee->offL[k]);
1266 }
1267 }
1268 break;
1269
1270 case ROW_FORMAT:
1271 for (i=0; i<numR; i++) {
1272 for (k=thee->IA[i]; k<thee->IA[i+1]; k++) {
1273 j = thee->JA[k];
1274 fprintf(fp, " %d %d %15.8e\n",
1275 i+1, j+1, thee->A[k]);
1276 }
1277 }
1278 break;
1279
1280 case COL_FORMAT:
1281 for (j=0; j<numC; j++) {
1282 for (k=thee->IA[j]; k<thee->IA[j+1]; k++) {
1283 i = thee->JA[k];
1284 fprintf(fp, " %d %d %15.8e\n",
1285 i+1, j+1, thee->A[k]);
1286 }
1287 }
1288 break;
1289
1290 default:
1291 Vnm_print(0,
1292 "%smatrix not in correct format to print....skipping.\n", rn);
1293 break;
1294 }
1295
1296 fprintf(fp, "];\n\n");
1297 fprintf(fp, "fprintf('..done.\\n');\n");
1298 }
1299
1300 fprintf(fp,"%% ----------------------------------------\n");
1301 fprintf(fp,"%% Matlab code to generate sparse matrix.\n");
1302 fprintf(fp,"%% ----------------------------------------\n");
1303 if( thee->state != NULL_STATE ) {
1304 fprintf(fp,
1305 "%s = sparse(%s_IJA(:,1),%s_IJA(:,2),%s_IJA(:,3),%d,%d);\n",
1306 thee->name, thee->name, thee->name, thee->name, numR, numC);
1307 } else {
1308 fprintf(fp, "%s = sparse(%d,%d);\n", thee->name, numR, numC);
1309 }
1310 fprintf(fp,"%% ----------------------------------------\n");
1311
1312 Vnm_print(0,"..done.\n");
1313
1314 /* close file and return */
1315 fclose(fp);
1316 }
1317
1318 /*
1319 * ***************************************************************************
1320 * Routine: Mat_printNoD
1321 *
1322 * Purpose: Print the matrix as a DENSE matrix in MATLAB format,
1323 * but first zero out any rows/cols corresponding to
1324 * Dirichlet boundary points.
1325 *
1326 * Notes: This routine is useful for e.g. checking that Galerkin
1327 * conditions hold for stiffness matrices. Removing the
1328 * dirichlet equations is crucial; otherwise the Galerkin
1329 * condition cannot hold. Note that the matrix (and the
1330 * Galerkin coarse matrix) are then of course singular.
1331 *
1332 * Author: Michael Holst
1333 * ***************************************************************************
1334 */
Mat_printNoD(Mat * thee)1335 VPUBLIC void Mat_printNoD(Mat *thee)
1336 {
1337 int i, j, k;
1338 char rn[80];
1339 const int MaxRows = 100;
1340 const int MaxCols = 30;
1341 double matrix[30][30];
1342
1343 strncpy(rn,"Mat_printNoD:",80);
1344
1345 /* some i/o */
1346 Vnm_print(0, "%s printing <%s>" " [dim=(%dx%d),sym=%d,numA=%d,numO=%d]\n",
1347 rn, thee->name, thee->numR, thee->numC, thee->sym,
1348 thee->numA, thee->numO);
1349
1350 /* size check */
1351 if ((thee->numR > MaxRows) || (thee->numC > MaxCols)) {
1352 Vnm_print(0,
1353 "%smatrix too large to view....skipping.\n", rn);
1354 return;
1355 }
1356
1357 /* make a dense matrix first */
1358 for (i=0; i<thee->numR; i++)
1359 for (j=0; j<thee->numC; j++)
1360 matrix[i][j] = 0.0;
1361
1362 switch (thee->format) {
1363
1364 case DRC_FORMAT:
1365 for (i=0; i<thee->numR; i++) {
1366 matrix[i][i] = thee->diag[i];
1367 for (k=thee->IA[i]; k<thee->IA[i+1]; k++) {
1368 j = thee->JA[k];
1369 matrix[i][j] = thee->offU[k];
1370 matrix[j][i] = thee->offL[k];
1371 }
1372 }
1373 break;
1374
1375 case ROW_FORMAT:
1376 for (i=0; i<thee->numR; i++) {
1377 for (k=thee->IA[i]; k<thee->IA[i+1]; k++) {
1378 j = thee->JA[k];
1379 matrix[i][j] = thee->A[k];
1380 }
1381 }
1382 break;
1383
1384 case COL_FORMAT:
1385 for (j=0; j<thee->numC; j++) {
1386 for (k=thee->IA[j]; k<thee->IA[j+1]; k++) {
1387 i = thee->JA[k];
1388 matrix[i][j] = thee->A[k];
1389 }
1390 }
1391 break;
1392
1393 default:
1394 Vnm_print(0,
1395 "%smatrix not in correct format to print....skipping.\n", rn);
1396 break;
1397 }
1398
1399 /* pull out the dirichlet rows and columns */
1400 for (k=0; k<thee->numBR; k++) {
1401 for (j=0; j<thee->numC; j++) {
1402 matrix[ thee->BR[k] ][ j ] = 0.0;
1403 }
1404 }
1405 for (k=0; k<thee->numBC; k++) {
1406 for (i=0; i<thee->numR; i++) {
1407 matrix[ i ][ thee->BC[k] ] = 0.0;
1408 }
1409 }
1410
1411 /* print the matrix */
1412 Vnm_print(0, "%s = [\n", thee->name);
1413 for (i=0; i<thee->numR; i++) {
1414 for (j=0; j<thee->numC; j++) {
1415 if (VABS(matrix[i][j]) < 0.0001) {
1416 Vnm_print(0, " 0.0 ");
1417 } else {
1418 Vnm_print(0, "%7.3f", matrix[i][j]);
1419 }
1420 }
1421 Vnm_print(0, "\n");
1422 }
1423 Vnm_print(0, "];\n");
1424 }
1425
1426 /*
1427 * ***************************************************************************
1428 * Routine: Mat_printSpNoD
1429 *
1430 * Purpose: Print the matrix as a SPARSE matrix in MATLAB format,
1431 * but first zero out any rows/cols corresponding to
1432 * Dirichlet boundary points.
1433 *
1434 * Notes: This routine is useful for e.g. checking that Galerkin
1435 * conditions hold for stiffness matrices. Removing the
1436 * dirichlet equations is crucial; otherwise the Galerkin
1437 * condition cannot hold. Note that the matrix (and the
1438 * Galerkin coarse matrix) are then of course singular.
1439 *
1440 * flag==0 ==> write
1441 * flag==1 ==> append
1442 *
1443 * Author: Michael Holst
1444 * ***************************************************************************
1445 */
Mat_printSpNoD(Mat * thee,char * fname,int pflag)1446 VPUBLIC void Mat_printSpNoD(Mat *thee, char *fname, int pflag)
1447 {
1448 int i, j, k, ii, i_do;
1449 char rn[80];
1450 FILE *fp;
1451
1452 strncpy(rn,"Mat_printSpNoD:",80);
1453
1454 if (pflag == 0) {
1455 fp=fopen(fname,"w");
1456 } else if (pflag == 1) {
1457 fp=fopen(fname,"a");
1458 } else {
1459 fp=VNULL;
1460 }
1461 if (fp==VNULL) {
1462 Vnm_print(2, "%s problem opening file <%s>\n", rn, fname);
1463 return;
1464 }
1465
1466 Vnm_print(0, "%s printing <%s>" " [dim=(%dx%d),sym=%d,numA=%d,numO=%d]"
1467 " to file <%s>\n",
1468 rn, thee->name, thee->numR, thee->numC, thee->sym,
1469 thee->numA, thee->numO, fname);
1470
1471 /* print the matrix in matlab sparse format */
1472 fprintf(fp,"%% ----------------------------------------\n");
1473 fprintf(fp,"%% %s matrix <%s> [dim=(%dx%d), sym=%d]\n",
1474 rn,thee->name,thee->numR,thee->numC,thee->sym);
1475 fprintf(fp,"%% ----------------------------------------\n");
1476 fprintf(fp, "fprintf('Defining %s..');", thee->name);
1477 fprintf(fp, "\n%s_IJA = [\n", thee->name);
1478
1479 switch (thee->format) {
1480
1481 case DRC_FORMAT:
1482 for (i=0; i<thee->numR; i++) {
1483 i_do = 1;
1484 for (ii=0; ii<thee->numBR; ii++)
1485 if (thee->BR[ii] == i) i_do = 0;
1486 if (i_do)
1487 fprintf(fp, " %d %d %15.8e\n",
1488 i+1, i+1, thee->diag[i]);
1489 }
1490
1491 for (i=0; i<thee->numR; i++) {
1492 for (k=thee->IA[i]; k<thee->IA[i+1]; k++) {
1493 j = thee->JA[k];
1494 i_do = 1;
1495 for (ii=0; ii<thee->numBR; ii++)
1496 if ((thee->BR[ii] == i) || (thee->BR[ii] == j)) {
1497 i_do = 0;
1498 }
1499 if (i_do)
1500 fprintf(fp, " %d %d %15.8e\n",
1501 i+1, j+1, thee->offU[k]);
1502 }
1503 }
1504
1505 for (j=0; j<thee->numR; j++) {
1506 for (k=thee->IA[j]; k<thee->IA[j+1]; k++) {
1507 i = thee->JA[k];
1508 i_do = 1;
1509 for (ii=0; ii<thee->numBR; ii++)
1510 if ((thee->BR[ii] == i) || (thee->BR[ii] == j)) {
1511 i_do = 0;
1512 }
1513 if (i_do)
1514 fprintf(fp, " %d %d %15.8e\n",
1515 i+1, j+1, thee->offL[k]);
1516 }
1517 }
1518 break;
1519
1520 case ROW_FORMAT:
1521 for (i=0; i<thee->numR; i++) {
1522 for (k=thee->IA[i]; k<thee->IA[i+1]; k++) {
1523 j = thee->JA[k];
1524 i_do = 1;
1525 for (ii=0; ii<thee->numBR; ii++)
1526 if ((thee->BR[ii] == i) || (thee->BR[ii] == j)) {
1527 i_do = 0;
1528 }
1529 if (i_do)
1530 fprintf(fp, " %d %d %15.8e\n",
1531 i+1, j+1, thee->A[k]);
1532 }
1533 }
1534 break;
1535
1536 case COL_FORMAT:
1537 for (j=0; j<thee->numC; j++) {
1538 for (k=thee->IA[j]; k<thee->IA[j+1]; k++) {
1539 i = thee->JA[k];
1540 i_do = 1;
1541 for (ii=0; ii<thee->numBR; ii++)
1542 if ((thee->BR[ii] == i) || (thee->BR[ii] == j)) {
1543 i_do = 0;
1544 }
1545 if (i_do)
1546 fprintf(fp, " %d %d %15.8e\n",
1547 i+1, j+1, thee->A[k]);
1548 }
1549 }
1550 break;
1551
1552 default:
1553 Vnm_print(0,
1554 "%smatrix not in correct format to print....skipping.\n", rn);
1555 break;
1556 }
1557
1558 fprintf(fp, "];\n");
1559 fprintf(fp, "fprintf('..done.\\n');\n");
1560 fprintf(fp,"%% ----------------------------------------\n");
1561 fprintf(fp,"%% Matlab code to generate sparse matrix.\n");
1562 fprintf(fp,"%% ----------------------------------------\n");
1563 fprintf(fp,
1564 "%s = sparse(%s_IJA(:,1),%s_IJA(:,2),%s_IJA(:,3),%d,%d);\n",
1565 thee->name, thee->name, thee->name, thee->name,
1566 thee->numR, thee->numC);
1567 fprintf(fp,"%% ----------------------------------------\n");
1568 fprintf(fp, "\n\n");
1569
1570 Vnm_print(0,"..done.\n");
1571
1572 /* close file and return */
1573 fclose(fp);
1574 }
1575
1576 /*
1577 * ***************************************************************************
1578 * Routine: Mat_zero
1579 *
1580 * Purpose: Clear the floating point storage for the sparse matrix.
1581 * Also clear any sparse factorization storage.
1582 *
1583 * Notes: This is basically done in preparation for an accumulation as
1584 * part of a matrix assembly, and before a new sparse factorization.
1585 *
1586 * Author: Michael Holst
1587 * ***************************************************************************
1588 */
Mat_zero(Mat * thee)1589 VPUBLIC void Mat_zero(Mat *thee)
1590 {
1591 int i;
1592
1593 /* zero out all of the floating point storage */
1594 if (thee->state != NULL_STATE) {
1595
1596 switch (thee->format) {
1597
1598 case DRC_FORMAT:
1599 for (i=0; i<thee->numA; i++) {
1600 thee->A[i] = 0.;
1601 }
1602 break;
1603
1604 case ROW_FORMAT:
1605 for (i=0; i<thee->numA; i++) {
1606 thee->A[i] = 0.;
1607 }
1608 break;
1609
1610 case COL_FORMAT:
1611 if (thee->sym != 2) {
1612 for (i=0; i<thee->numA; i++) {
1613 thee->A[i] = 0.;
1614 }
1615 }
1616 break;
1617
1618 case SLU_FORMAT:
1619 for (i=0; i<thee->numA; i++) {
1620 thee->A[i] = 0.;
1621 }
1622 Mat_sluDestroy(thee);
1623 break;
1624
1625 case RLN_FORMAT:
1626 if (thee->lnkU != VNULL) {
1627 Vset_reset( thee->lnkU );
1628 }
1629 break;
1630
1631 case CLN_FORMAT:
1632 if (thee->lnkL != VNULL) {
1633 Vset_reset( thee->lnkL );
1634 }
1635 break;
1636
1637 case XLN_FORMAT:
1638 break;
1639
1640 case RFL_FORMAT:
1641 for (i=0; i<thee->numA; i++) {
1642 thee->A[i] = 0.;
1643 }
1644 break;
1645
1646 case CFL_FORMAT:
1647 for (i=0; i<thee->numA; i++) {
1648 thee->A[i] = 0.;
1649 }
1650 break;
1651
1652 default:
1653 VASSERT(0);
1654 break;
1655 }
1656
1657 thee->state = ZERO_STATE;
1658 }
1659 }
1660
1661 /*
1662 * ***************************************************************************
1663 * Routine: Mat_set
1664 *
1665 * Purpose: Set a value in a matrix.
1666 *
1667 * Author: Michael Holst
1668 * ***************************************************************************
1669 */
Mat_set(Mat * thee,int i,int j,double val)1670 VPUBLIC void Mat_set(Mat *thee, int i, int j, double val)
1671 {
1672 VASSERT( thee->state != NULL_STATE );
1673
1674 switch (thee->format) {
1675
1676 case DRC_FORMAT:
1677 thee->state = ASSEMBLED_STATE;
1678 if (i==j) thee->diag[i] = val;
1679 else if (i<j) mPlaceit(thee->IA, thee->JA, thee->offU, 0, i, j, val);
1680 else mPlaceit(thee->IA, thee->JA, thee->offL, 0, j, i, val);
1681 break;
1682
1683 case ROW_FORMAT:
1684 thee->state = ASSEMBLED_STATE;
1685 mPlaceit(thee->IA, thee->JA, thee->A, 0, i, j, val);
1686 break;
1687
1688 case COL_FORMAT:
1689 thee->state = ASSEMBLED_STATE;
1690 mPlaceit(thee->IA, thee->JA, thee->A, 0, j, i, val);
1691 break;
1692
1693 case RFL_FORMAT:
1694 thee->state = ASSEMBLED_STATE;
1695 thee->A[ i*(thee->numC) + j ] = val;
1696 break;
1697
1698 case CFL_FORMAT:
1699 thee->state = ASSEMBLED_STATE;
1700 thee->A[ j*(thee->numR) + i ] = val;
1701 break;
1702
1703 default:
1704 VASSERT(0);
1705 break;
1706 }
1707 }
1708
1709 /*
1710 * ***************************************************************************
1711 * Routine: Mat_addTo
1712 *
1713 * Purpose: Add to a value in a matrix.
1714 *
1715 * Author: Michael Holst
1716 * ***************************************************************************
1717 */
Mat_addTo(Mat * thee,int i,int j,double val)1718 VPUBLIC void Mat_addTo(Mat *thee, int i, int j, double val)
1719 {
1720 VASSERT( thee->state != NULL_STATE );
1721
1722 switch (thee->format) {
1723
1724 case DRC_FORMAT:
1725 thee->state = ASSEMBLED_STATE;
1726 if (i==j) thee->diag[i] += val;
1727 else if (i<j) mPlaceit(thee->IA, thee->JA, thee->offU, 1, i, j, val);
1728 else mPlaceit(thee->IA, thee->JA, thee->offL, 1, j, i, val);
1729 break;
1730
1731 case ROW_FORMAT:
1732 thee->state = ASSEMBLED_STATE;
1733 mPlaceit(thee->IA, thee->JA, thee->A, 1, i, j, val);
1734 break;
1735
1736 case COL_FORMAT:
1737 thee->state = ASSEMBLED_STATE;
1738 mPlaceit(thee->IA, thee->JA, thee->A, 1, j, i, val);
1739 break;
1740
1741 case RFL_FORMAT:
1742 thee->state = ASSEMBLED_STATE;
1743 thee->A[ i*(thee->numC) + j ] += val;
1744 break;
1745
1746 case CFL_FORMAT:
1747 thee->state = ASSEMBLED_STATE;
1748 thee->A[ j*(thee->numR) + i ] += val;
1749 break;
1750
1751 default:
1752 VASSERT(0);
1753 break;
1754 }
1755 }
1756
1757 /*
1758 * ***************************************************************************
1759 * Routine: Mat_buildBRC
1760 *
1761 * Purpose: Set the boundary row and column information.
1762 *
1763 * Notes: key=0 ==> set pointers
1764 * key=1 ==> do a copy
1765 *
1766 * Author: Michael Holst
1767 * ***************************************************************************
1768 */
Mat_buildBRC(Mat * thee,int numBR,int numBC,int * BR,int * BC)1769 VPUBLIC void Mat_buildBRC(Mat *thee, int numBR, int numBC, int *BR, int *BC)
1770 {
1771 int i;
1772
1773 if ( thee->numBR > 0 ) {
1774 VASSERT( thee->BR != VNULL );
1775 Vmem_free( thee->vmem, thee->numBR, sizeof(int),
1776 (void**)&(thee->BR) );
1777 }
1778 if ( thee->numBC > 0 ) {
1779 VASSERT( thee->BC != VNULL );
1780 Vmem_free( thee->vmem, thee->numBC, sizeof(int),
1781 (void**)&(thee->BC) );
1782 }
1783 thee->BR = VNULL;
1784 thee->BC = VNULL;
1785
1786 thee->numBR = numBR;
1787 thee->numBC = numBC;
1788
1789 if (thee->numBR > 0) {
1790 thee->BR = Vmem_malloc(thee->vmem,thee->numBR,sizeof(int));
1791 for (i=0; i<thee->numBR; i++) {
1792 thee->BR[i] = BR[i];
1793 }
1794 }
1795 if (thee->numBC > 0) {
1796 thee->BC = Vmem_malloc(thee->vmem,thee->numBC,sizeof(int));
1797 for (i=0; i<thee->numBC; i++) {
1798 thee->BC[i] = BC[i];
1799 }
1800 }
1801 }
1802
1803 /*
1804 * ***************************************************************************
1805 * Routine: Mat_zeroBRC
1806 *
1807 * Purpose: Apply the boundary row and column information.
1808 *
1809 * Author: Michael Holst
1810 * ***************************************************************************
1811 */
Mat_zeroBRC(Mat * thee)1812 VPUBLIC void Mat_zeroBRC(Mat *thee)
1813 {
1814 int i, j, k, ibase, *BR, *BC;
1815
1816 /* this only makes sense for prolongation matrix at this point */
1817 VASSERT( thee->format == ROW_FORMAT );
1818
1819 /* offset row indices in the case of implicit identity */
1820 ibase = (thee->impl == IS_IMPL) ? thee->numC : 0;
1821
1822 /* create temporary boundary row/col info */
1823 BR = Vmem_malloc( thee->vmem, thee->numR + ibase, sizeof(int) );
1824 BC = Vmem_malloc( thee->vmem, thee->numC, sizeof(int) );
1825
1826 for (i=0; i<thee->numR; i++) {
1827 BR[ i ] = 0;
1828 }
1829 for (i=0; i<thee->numBR; i++) {
1830 BR[ thee->BR[i] ] = 1;
1831 }
1832 for (i=0; i<thee->numC; i++) {
1833 BC[ i ] = 0;
1834 }
1835 for (i=0; i<thee->numBC; i++) {
1836 BC[ thee->BC[i] ] = 1;
1837 }
1838
1839 /* bring in the PDE: zero the dirichlet rows */
1840 for (i=0; i<thee->numR; i++) {
1841 for (k=thee->IA[i]; k<thee->IA[i+1]; k++) {
1842 j = thee->JA[k];
1843 if ( VDIRICHLET( BR[i+ibase] ) ) {
1844 thee->A[ k ] = 0.0;
1845 } else if ( VDIRICHLET( BC[j] ) ) {
1846 thee->A[ k ] = 0.0;
1847 }
1848 }
1849 }
1850
1851 /* destroy temporary boundary row/col info */
1852 Vmem_free( thee->vmem, thee->numR + ibase, sizeof(int), (void**)&BR );
1853 Vmem_free( thee->vmem, thee->numC, sizeof(int), (void**)&BC );
1854 }
1855
1856 /*
1857 * ***************************************************************************
1858 * Routine: Mat_diagBRC
1859 *
1860 * Purpose: Place identity entries on diagonal of boundary row/col.
1861 *
1862 * Author: Michael Holst
1863 * ***************************************************************************
1864 */
Mat_diagBRC(Mat * thee)1865 VPUBLIC void Mat_diagBRC(Mat *thee)
1866 {
1867 int i;
1868
1869 /* this only makes sense for DRC matrices */
1870 VASSERT( thee->format == DRC_FORMAT );
1871 VASSERT( thee->numR == thee->numC );
1872 VASSERT( thee->numBR == thee->numBC );
1873
1874 /* put identity on the boundary rows */
1875 for (i=0; i<thee->numBR; i++) {
1876 thee->diag[ thee->BR[i] ] = 1.0;
1877 }
1878 }
1879
1880 /*
1881 * ***************************************************************************
1882 * Routine: Mat_galerkin
1883 *
1884 * Purpose: Enforce the Galerkin conditions algebraically.
1885 *
1886 * Author: Michael Holst
1887 * ***************************************************************************
1888 */
Mat_galerkin(Mat * thee,Mat * rmat,Mat * amat,Mat * pmat)1889 VPUBLIC void Mat_galerkin(Mat *thee, Mat *rmat, Mat *amat, Mat *pmat)
1890 {
1891 int numBR, numBC, *BR, *BC;
1892
1893 /* are we ready to roll? */
1894 VASSERT( thee->state == NULL_STATE );
1895
1896 /* do dimensions of R,A,P agree for a product? */
1897 VASSERT( amat->numR == Mat_numR(rmat) );
1898 VASSERT( amat->numC == Mat_numR(pmat) );
1899
1900 /* some additional consistency checks */
1901 VASSERT( thee->numR == Mat_numC(rmat) );
1902 VASSERT( thee->numC == Mat_numC(pmat) );
1903
1904 /* form the triple matrix product */
1905 buildG(thee->vmem, amat->format, amat->sym, thee->numR, amat->numR,
1906
1907 &(thee->numO),
1908 &(thee->numA),
1909 &(thee->IJA),
1910 &(thee->A),
1911 &(thee->IA),
1912 &(thee->JA),
1913 &(thee->diag),
1914 &(thee->offU),
1915 &(thee->offL),
1916
1917 Mat_IA(rmat), Mat_JA(rmat), Mat_A(rmat),
1918
1919 Mat_IA(amat), Mat_JA(amat),
1920 Mat_A(amat), Mat_diag(amat),
1921 Mat_offU(amat), Mat_offL(amat),
1922
1923 Mat_IA(pmat), Mat_JA(pmat), Mat_A(pmat)
1924 );
1925
1926 /* we have structure now */
1927 thee->state = ASSEMBLED_STATE;
1928 thee->format = amat->format;
1929 thee->sym = amat->sym;
1930 thee->impl = amat->impl;
1931
1932 /* buildG allocated the space */
1933 thee->iMallocIJA = 1;
1934 if (thee->numA > 0 ) {
1935 thee->iMallocA = 1;
1936 }
1937
1938 /* the numZ parameter */
1939 if (thee->format == DRC_FORMAT) {
1940 if (thee->sym == IS_SYM) {
1941 thee->numZ = thee->numA + thee->numO;
1942 } else { /* (thee->sym == ISNOT_SYM) */
1943 thee->numZ = thee->numA;
1944 }
1945 } else { /* (thee->format != DRC_FORMAT) */
1946 thee->numZ = thee->numA;
1947 }
1948
1949 /* the BRC structure */
1950 numBR = rmat->numBC;
1951 numBC = pmat->numBC;
1952 BR = rmat->BC;
1953 BC = pmat->BC;
1954 Mat_buildBRC(thee, numBR, numBC, BR, BC);
1955 }
1956
1957 /*
1958 * ***************************************************************************
1959 * Routine: Mat_sluDirect
1960 *
1961 * Purpose: Make a decision about whether or not a sparse direct solver
1962 * should be used in place of an iterative solver, based on the
1963 * size of the system.
1964 *
1965 * Notes: This is obviously heuristic in nature; in general the cutoff
1966 * size where iterative methods start to win is smaller in 3D.
1967 *
1968 * Author: Michael Holst
1969 * ***************************************************************************
1970 */
Mat_sluDirect(Mat * thee)1971 VPUBLIC int Mat_sluDirect(Mat *thee)
1972 {
1973 int ival = 0;
1974 /* VASSERT( thee->format == SLU_FORMAT ); */
1975 /* the decision about when sparse direct will be the best solver */
1976 if ( thee->numR < SPARSE_CUTOFF ) ival = 1;
1977 return ival;
1978 }
1979
1980 /*
1981 * ***************************************************************************
1982 * Routine: Mat_sluCreate
1983 *
1984 * Purpose: Setup for a sparse LU factorization of matrix.
1985 *
1986 * Notes: Creates internal <ia,ja,a> storage which is later freed
1987 * by Mat_sluDestroy. Also initializes the sparse direct
1988 * library.
1989 *
1990 * Author: Michael Holst
1991 * ***************************************************************************
1992 */
Mat_sluCreate(Mat * thee)1993 VPUBLIC void Mat_sluCreate(Mat *thee)
1994 {
1995 VASSERT( thee->format == SLU_FORMAT );
1996 }
1997
1998 /*
1999 * ***************************************************************************
2000 * Routine: Mat_sluFactor
2001 *
2002 * Purpose: Create the sparse LU factors for the system matrix.
2003 *
2004 * Author: Michael Holst
2005 * ***************************************************************************
2006 */
Mat_sluFactor(Mat * thee)2007 VPUBLIC int Mat_sluFactor(Mat *thee)
2008 {
2009 int rc = 0; /* 1=factors okay, 0=factors not okay */
2010
2011 VASSERT( thee->format == SLU_FORMAT );
2012 VASSERT( thee->numO == thee->numA );
2013 VASSERT( thee->numO == thee->numZ );
2014
2015 /* matrix has already been sparse factored */
2016 if (thee->state == FACTORED_STATE) {
2017
2018 rc = 1;
2019
2020 /* if it hasn't been factored, then free existing factors and rebuild */
2021 } else if (thee->state == ZERO_STATE) {
2022
2023 VDEBUGIO("Mat_sluFactor: ");
2024
2025 /* first any free space currently taken for global COL matrix */
2026 VDEBUGIO("clearing..");
2027 Mat_sluDestroy(thee);
2028
2029 /* initialize the slu object */
2030 VDEBUGIO("..Slu_ctor..");
2031 VASSERT( thee->slu == VNULL );
2032 thee->slu = Slu_ctor( thee->vmem, 1, thee->numR, thee->numR,
2033 thee->numA, thee->IA, thee->JA, thee->A);
2034
2035 /* attempt to factor the matrix */
2036 VDEBUGIO("..factoring..");
2037 rc = Slu_factor(thee->slu);
2038 VDEBUGIO("..done.\n");
2039
2040 /* matrix has been successfully sparse factored */
2041 if (rc == 1) thee->state = FACTORED_STATE;
2042
2043 /* matrix is not ready to factor */
2044 } else if (thee->state == NULL_STATE) {
2045 VASSERT( 0 );
2046 rc = 0;
2047
2048 /* there is a very serious problem */
2049 } else {
2050 VASSERT( 0 );
2051 rc = 0;
2052 }
2053
2054 return rc;
2055 }
2056
2057 /*
2058 * ***************************************************************************
2059 * Routine: Mat_sluSolve
2060 *
2061 * Purpose: Performs a forward/backward solve using the sparse LU factors.
2062 *
2063 * Notes: This requires that Mat_sluFactor has been previously called.
2064 *
2065 * Author: Michael Holst
2066 * ***************************************************************************
2067 */
Mat_sluSolve(Mat * thee,int key,double * f,double * u)2068 VPUBLIC int Mat_sluSolve(Mat *thee, int key, double *f, double *u)
2069 {
2070 int rc;
2071
2072 VASSERT( thee->format == SLU_FORMAT );
2073
2074 rc = Slu_solve(thee->slu, key, f, u);
2075 return rc;
2076 }
2077
2078 /*
2079 * ***************************************************************************
2080 * Routine: Mat_sluDestroy
2081 *
2082 * Purpose: Destroy the sparse LU factors for the system matrix.
2083 *
2084 * Notes: This frees our <ia,ja,a> storage, and also the internal
2085 * storage that was malloc'd by the sparse direct library.
2086 *
2087 * Author: Michael Holst
2088 * ***************************************************************************
2089 */
Mat_sluDestroy(Mat * thee)2090 VPUBLIC void Mat_sluDestroy(Mat *thee)
2091 {
2092 VASSERT( thee->format == SLU_FORMAT );
2093
2094 if (thee->slu != VNULL) {
2095 if (thee->state != NULL_STATE ) {
2096
2097 /* this frees <ia,ja,a> as well as everything else */
2098 Slu_dtor(&(thee->slu));
2099
2100 /* new state */
2101 thee->state = NULL_STATE;
2102 }
2103 }
2104 }
2105
2106 /*
2107 * ***************************************************************************
2108 * Routine: Mat_memChk
2109 *
2110 * Purpose: Print the exact current malloc usage
2111 *
2112 * Author: Michael Holst
2113 * ***************************************************************************
2114 */
Mat_memChk(Mat * thee)2115 VPUBLIC void Mat_memChk(Mat *thee)
2116 {
2117 if (thee->iMadeVmem) Vmem_print(thee->vmem);
2118 }
2119
2120