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: bmat.c,v 1.44 2010/08/12 05:18:24 fetk Exp $"
21 * ***************************************************************************
22 */
23
24 /*
25 * ***************************************************************************
26 * File: bmat.c
27 *
28 * Purpose: Class Bmat: methods.
29 *
30 * Author: Michael Holst
31 * ***************************************************************************
32 */
33
34 #include "bmat_p.h"
35
36 VEMBED(rcsid="$Id: bmat.c,v 1.44 2010/08/12 05:18:24 fetk Exp $")
37
38 /*
39 * ***************************************************************************
40 * Class Bmat: Inlineable methods
41 * ***************************************************************************
42 */
43 #if !defined(VINLINE_BAM)
44
45 #endif /* if !defined(VINLINE_BAM) */
46 /*
47 * ***************************************************************************
48 * Class Bmat: Non-inlineable methods
49 * ***************************************************************************
50 */
51
52 /*
53 * ***************************************************************************
54 * Routine: Bmat_ctor
55 *
56 * Purpose: The block sparse matrix constructor.
57 *
58 * Notes: This constructor only does minimal initialization of a Bmat
59 * object after creating the object itself; it doesn't create
60 * any storage for either the integer structure arrays or the
61 * nonzero storage arrays.
62 *
63 * This constructor only fixes the number of blocks and the
64 * numbers of rows and columns in each block; the nonzero
65 * structure is not set yet.
66 *
67 * Author: Michael Holst
68 * ***************************************************************************
69 */
Bmat_ctor(Vmem * vmem,const char * name,int pnumB,int pnumR[MAXV],int pnumC[MAXV],MATmirror pmirror[MAXV][MAXV])70 VPUBLIC Bmat* Bmat_ctor(Vmem *vmem, const char *name,
71 int pnumB, int pnumR[MAXV], int pnumC[MAXV],
72 MATmirror pmirror[MAXV][MAXV])
73 {
74 int p,q;
75 char bname[15];
76 Bmat *thee = VNULL;
77
78 VDEBUGIO("Bmat_ctor: CREATING object..");
79
80 thee = Vmem_malloc( VNULL, 1, sizeof(Bmat) );
81 if (vmem == VNULL) {
82 thee->vmem = Vmem_ctor( "Bmat" );
83 thee->iMadeVmem = 1;
84 } else {
85 thee->vmem = vmem;
86 thee->iMadeVmem = 0;
87 }
88
89 strncpy(thee->name, name, 10);
90 thee->coarse = VNULL;
91 thee->fine = VNULL;
92 thee->numB = pnumB;
93 for (p=0; p<thee->numB; p++) {
94 for (q=0; q<thee->numB; q++) {
95 thee->mirror[p][q] = pmirror[p][q];
96 thee->AD[p][q] = VNULL;
97 }
98 }
99
100 /* create upper-triangular blocks before the lower-triangular blocks */
101 for (p=0; p<thee->numB; p++) {
102 for (q=0; q<thee->numB; q++) {
103 if ( !thee->mirror[p][q] ) {
104 sprintf(bname, "%s_%d%d", thee->name, p, q);
105 thee->AD[p][q] = Mat_ctor(thee->vmem,bname,pnumR[p],pnumC[q]);
106 } else { /* ( thee->mirror[p][q] ) */
107 /* first make sure that the mirror of this guy exists! */
108 VASSERT( !thee->mirror[q][p] );
109 VASSERT( thee->AD[q][p] != VNULL );
110 thee->AD[p][q] = thee->AD[q][p];
111 }
112 }
113 }
114
115 VDEBUGIO("..done.\n");
116
117 return thee;
118 }
119
120 /*
121 * ***************************************************************************
122 * Routine: Bmat_dtor
123 *
124 * Purpose: The block sparse matrix destructor.
125 *
126 * Notes: This destructor does the reverse of Bmat_ctor, and if
127 * necessary first reverses Bmat_initStructure
128 * (or Bmat_copyStructure). I.e., if necessary,
129 * it first frees the large integer and real arrays created by
130 * Bmat_initStructure or Bmat_copyStructure, and then frees the
131 * Bmat object itself at the last moment.
132 *
133 * Author: Michael Holst
134 * ***************************************************************************
135 */
Bmat_dtor(Bmat ** thee)136 VPUBLIC void Bmat_dtor(Bmat **thee)
137 {
138 int p, q;
139
140 /* VASSERT( (*thee) != VNULL ); */
141 if ((*thee) != VNULL) {
142
143 /* first free any coarse guy */
144 if ( (*thee)->coarse != VNULL ) {
145 Bmat_dtor( &((*thee)->coarse) );
146 }
147
148 /* next free the large integer and real nonzero structure */
149 Bmat_killStructure(*thee);
150
151 /* kill the submatrices themselves */
152 Mat_dtor( &((*thee)->AG) );
153 for (p=0; p<(*thee)->numB; p++) {
154 for (q=0; q<(*thee)->numB; q++) {
155 if ( !(*thee)->mirror[p][q] ) {
156 Mat_dtor( &((*thee)->AD[p][q]) );
157 }
158 }
159 }
160
161 /* now kill the object itself */
162 VDEBUGIO("Bmat_dtor: DESTROYING object..");
163 if ((*thee)->iMadeVmem) Vmem_dtor( &((*thee)->vmem) );
164 Vmem_free( VNULL, 1, sizeof(Bmat), (void**)thee );
165 VDEBUGIO("..done.\n");
166
167 (*thee) = VNULL;
168 }
169 }
170
171 /*
172 * ***************************************************************************
173 * Routine: Bmat_initStructure
174 *
175 * Purpose: Initialize the nonzero structure given structure information.
176 *
177 * Notes: This routine actually does the storage creation for both the
178 * integer structure information arrays and the nonzero value
179 * arrays.
180 *
181 * Author: Michael Holst
182 * ***************************************************************************
183 */
Bmat_initStructure(Bmat * thee,MATformat pfrmt[MAXV][MAXV],MATsym psym[MAXV][MAXV],int pnumO[MAXV][MAXV],int * IJA[MAXV][MAXV])184 VPUBLIC void Bmat_initStructure(Bmat *thee,
185 MATformat pfrmt[MAXV][MAXV], MATsym psym[MAXV][MAXV],
186 int pnumO[MAXV][MAXV], int *IJA[MAXV][MAXV])
187 {
188 int p, q;
189
190 for (p=0; p<thee->numB; p++) {
191 for (q=0; q<thee->numB; q++) {
192 if ( !thee->mirror[p][q] ) {
193 Mat_initStructure(thee->AD[p][q],
194 pfrmt[p][q], psym[p][q], pnumO[p][q], IJA[p][q], VNULL);
195 }
196 }
197 }
198 }
199
200 /*
201 * ***************************************************************************
202 * Routine: Bmat_copyStructure
203 *
204 * Purpose: Copy the nonzero structure from an input model.
205 *
206 * Author: Michael Holst
207 * ***************************************************************************
208 */
Bmat_copyStructure(Bmat * thee,Bmat * model)209 VPUBLIC void Bmat_copyStructure(Bmat *thee, Bmat *model)
210 {
211 int p, q;
212
213 for (p=0; p<thee->numB; p++) {
214 for (q=0; q<thee->numB; q++) {
215 VASSERT( thee->mirror[p][q] == model->mirror[p][q] );
216 if ( !thee->mirror[p][q] ) {
217 Mat_copyStructure(thee->AD[p][q], model->AD[p][q]);
218 } else {
219 VASSERT( !thee->mirror[q][p] );
220 thee->AD[p][q] = thee->AD[q][p];
221 }
222 }
223 }
224 }
225
226 /*
227 * ***************************************************************************
228 * Routine: Bmat_killStructure
229 *
230 * Purpose: Kill the nonzero structure.
231 *
232 * Author: Michael Holst
233 * ***************************************************************************
234 */
Bmat_killStructure(Bmat * thee)235 VPUBLIC void Bmat_killStructure(Bmat *thee)
236 {
237 int p, q;
238
239 Mat_killStructure(thee->AG);
240 for (p=0; p<thee->numB; p++) {
241 for (q=0; q<thee->numB; q++) {
242 if ( !thee->mirror[p][q] ) {
243 Mat_killStructure(thee->AD[p][q]);
244 }
245 }
246 }
247 }
248
249 /*
250 * ***************************************************************************
251 * Routine: Bmat_numB
252 *
253 * Purpose: Return the number of blocks.
254 *
255 * Author: Michael Holst
256 * ***************************************************************************
257 */
Bmat_numB(Bmat * thee)258 VPUBLIC int Bmat_numB(Bmat *thee)
259 {
260 return thee->numB;
261 }
262
263 /*
264 * ***************************************************************************
265 * Routine: Bmat_numR
266 *
267 * Purpose: Return the number of rows.
268 *
269 * Notes: For a mirror block, AD[p][q] is a pointer to the mirror,
270 * AD[q][p], and hence, the dimensions must be transposed.
271 *
272 * Author: Michael Holst
273 * ***************************************************************************
274 */
Bmat_numR(Bmat * thee,int p,int q)275 VPUBLIC int Bmat_numR(Bmat *thee, int p, int q)
276 {
277 if( !(thee->mirror[p][q]) ) {
278 return Mat_numR( thee->AD[p][q] );
279 } else {
280 return Mat_numC( thee->AD[q][p] );
281 }
282 }
283
284 /*
285 * ***************************************************************************
286 * Routine: Bmat_numC
287 *
288 * Purpose: Return the number of columns.
289 *
290 * Notes: For a mirror block, AD[p][q] is a pointer to the mirror,
291 * AD[q][p], and hence, the dimensions must be transposed.
292 *
293 * Author: Michael Holst
294 * ***************************************************************************
295 */
Bmat_numC(Bmat * thee,int p,int q)296 VPUBLIC int Bmat_numC(Bmat *thee, int p, int q)
297 {
298 if( !(thee->mirror[p][q]) ) {
299 return Mat_numC( thee->AD[p][q] );
300 } else {
301 return Mat_numR( thee->AD[q][p] );
302 }
303 }
304
305 /*
306 * ***************************************************************************
307 * Routine: Bmat_numA
308 *
309 * Purpose: Return the number of nonzeros we actually store in a block.
310 *
311 * Author: Michael Holst
312 * ***************************************************************************
313 */
Bmat_numA(Bmat * thee,int p,int q)314 VPUBLIC int Bmat_numA(Bmat *thee, int p, int q)
315 {
316 return Mat_numA( thee->AD[p][q] );
317 }
318
319 /*
320 * ***************************************************************************
321 * Routine: Bmat_numO
322 *
323 * Purpose: Return the number of nonzeros we actually store in a block
324 * which are actually in the strict upper-triangle of the
325 * block (DRC only).
326 *
327 * Author: Michael Holst
328 * ***************************************************************************
329 */
Bmat_numO(Bmat * thee,int p,int q)330 VPUBLIC int Bmat_numO(Bmat *thee, int p, int q)
331 {
332 return Mat_numO( thee->AD[p][q] );
333 }
334
335 /*
336 * ***************************************************************************
337 * Routine: Bmat_numZ
338 *
339 * Purpose: Return the number of nonzeros we WOULD actually store in a block
340 * if we ignored symmetry and stored everything.
341 *
342 * Author: Michael Holst
343 * ***************************************************************************
344 */
Bmat_numZ(Bmat * thee,int p,int q)345 VPUBLIC int Bmat_numZ(Bmat *thee, int p, int q)
346 {
347 return Mat_numZ( thee->AD[p][q] );
348 }
349
350 /*
351 * ***************************************************************************
352 * Routine: Bmat_numRT
353 *
354 * Purpose: Return the total number of rows.
355 *
356 * Author: Michael Holst
357 * ***************************************************************************
358 */
Bmat_numRT(Bmat * thee)359 VPUBLIC int Bmat_numRT(Bmat *thee)
360 {
361 int p, numRT;
362
363 numRT = 0;
364 for (p=0; p<thee->numB; p++) {
365 numRT += Bmat_numR(thee,p,p);
366 }
367
368 return numRT;
369 }
370
371 /*
372 * ***************************************************************************
373 * Routine: Bmat_numCT
374 *
375 * Purpose: Return the total number of columns.
376 *
377 * Author: Michael Holst
378 * ***************************************************************************
379 */
Bmat_numCT(Bmat * thee)380 VPUBLIC int Bmat_numCT(Bmat *thee)
381 {
382 int q, numCT;
383
384 numCT = 0;
385 for (q=0; q<thee->numB; q++) {
386 numCT += Bmat_numC(thee,q,q);
387 }
388
389 return numCT;
390 }
391
392 /*
393 * ***************************************************************************
394 * Routine: Bmat_numAT
395 *
396 * Purpose: Return the total number of nonzeros we are actually storing.
397 *
398 * Author: Michael Holst
399 * ***************************************************************************
400 */
Bmat_numAT(Bmat * thee)401 VPUBLIC int Bmat_numAT(Bmat *thee)
402 {
403 int p, q, numAT;
404
405 numAT = 0;
406 for (p=0; p<thee->numB; p++) {
407 for (q=0; q<thee->numB; q++) {
408 numAT += Bmat_numA(thee,p,q);
409 }
410 }
411
412 return numAT;
413 }
414
415 /*
416 * ***************************************************************************
417 * Routine: Bmat_numOT
418 *
419 * Purpose: Return the total number of nonzeros we are actually storing
420 * which are located in the strict upper-triangle.
421 *
422 * Author: Michael Holst
423 * ***************************************************************************
424 */
Bmat_numOT(Bmat * thee)425 VPUBLIC int Bmat_numOT(Bmat *thee)
426 {
427 int p, q, numOT;
428
429 numOT = 0;
430 for (p=0; p<thee->numB; p++) {
431 for (q=0; q<thee->numB; q++) {
432 numOT += Bmat_numO(thee,p,q);
433 }
434 }
435
436 return numOT;
437 }
438
439 /*
440 * ***************************************************************************
441 * Routine: Bmat_numZT
442 *
443 * Purpose: Return the total number of nonzeros we WOULD be storing if we
444 * ignored all symmetry in all blocks and stored everything.
445 *
446 * Author: Michael Holst
447 * ***************************************************************************
448 */
Bmat_numZT(Bmat * thee)449 VPUBLIC int Bmat_numZT(Bmat *thee)
450 {
451 int p, q, numZT;
452
453 numZT = 0;
454 for (p=0; p<thee->numB; p++) {
455 for (q=0; q<thee->numB; q++) {
456 numZT += Bmat_numZ(thee,p,q);
457 }
458 }
459
460 return numZT;
461 }
462
463 /*
464 * ***************************************************************************
465 * Routine: Bmat_format
466 *
467 * Purpose: Return the format of the block.
468 *
469 * Notes: For a mirror block, AD[p][q] is a pointer to the mirror,
470 * AD[q][p], and hence, the format must be transposed.
471 *
472 * Author: Michael Holst and Stephen Bond
473 * ***************************************************************************
474 */
Bmat_format(Bmat * thee,int p,int q)475 VPUBLIC MATformat Bmat_format(Bmat *thee, int p, int q)
476 {
477 if( thee->mirror[p][q] ) { /* Is it a mirror block ? */
478 switch ( Mat_format( thee->AD[q][p] ) ) {
479 case ZERO_FORMAT:
480 return ZERO_FORMAT;
481 break;
482 case ROW_FORMAT:
483 return COL_FORMAT;
484 break;
485 case COL_FORMAT:
486 return ROW_FORMAT;
487 break;
488 default:
489 /* Mirrors of all other cases currently unsupported */
490 VASSERT( 0 );
491 }
492 }
493
494 return Mat_format( thee->AD[p][q] );
495 }
496
497 /*
498 * ***************************************************************************
499 * Routine: Bmat_sym
500 *
501 * Purpose: Return the symmetry of the block.
502 *
503 * Author: Michael Holst
504 * ***************************************************************************
505 */
Bmat_sym(Bmat * thee,int p,int q)506 VPUBLIC MATsym Bmat_sym(Bmat *thee, int p, int q)
507 {
508 return Mat_sym( thee->AD[p][q] );
509 }
510
511 /*
512 * ***************************************************************************
513 * Routine: Bmat_state
514 *
515 * Purpose: Return the state of the block.
516 *
517 * Author: Michael Holst
518 * ***************************************************************************
519 */
Bmat_state(Bmat * thee,int p,int q)520 VPUBLIC MATstate Bmat_state(Bmat *thee, int p, int q)
521 {
522 return Mat_state( thee->AD[p][q] );
523 }
524
525 /*
526 * ***************************************************************************
527 * Routine: Bmat_impl
528 *
529 * Purpose: Return the implicitness of the block.
530 *
531 * Author: Michael Holst
532 * ***************************************************************************
533 */
Bmat_impl(Bmat * thee,int p,int q)534 VPUBLIC MATimpl Bmat_impl(Bmat *thee, int p, int q)
535 {
536 return Mat_impl( thee->AD[p][q] );
537 }
538
539 /*
540 * ***************************************************************************
541 * Routine: Bmat_mirror
542 *
543 * Purpose: Return the mirror-ness of the block.
544 *
545 * Author: Michael Holst
546 * ***************************************************************************
547 */
Bmat_mirror(Bmat * thee,int p,int q)548 VPUBLIC MATmirror Bmat_mirror(Bmat *thee, int p, int q)
549 {
550 return thee->mirror[p][q];
551 }
552
553 /*
554 * ***************************************************************************
555 * Routine: Bmat_sizeA
556 *
557 * Purpose: Return the number of nonzeros in all blocks.
558 *
559 * Author: Michael Holst
560 * ***************************************************************************
561 */
Bmat_sizeA(Bmat * thee)562 VPUBLIC int Bmat_sizeA(Bmat *thee)
563 {
564 int p,q,size = 0;
565 for (p=0; p<thee->numB; p++) {
566 for (q=0; q<thee->numB; q++) {
567 size += Mat_sizeA( thee->AD[p][q] );
568 }
569 }
570 return size;
571 }
572
573 /*
574 * ***************************************************************************
575 * Routine: Bmat_sizeIJA
576 *
577 * Purpose: Return the numer of integer storage locations used.
578 *
579 * Author: Michael Holst
580 * ***************************************************************************
581 */
Bmat_sizeIJA(Bmat * thee)582 VPUBLIC int Bmat_sizeIJA(Bmat *thee)
583 {
584 int p,q,size = 0;
585
586 for (p=0; p<thee->numB; p++) {
587 for (q=0; q<thee->numB; q++) {
588 size += Mat_sizeIJA( thee->AD[p][q] );
589 }
590 }
591
592 return size;
593 }
594
595 /*
596 * ***************************************************************************
597 * Routine: Bmat_IJA
598 *
599 * Purpose: Return the integer structure IJA.
600 *
601 * Author: Michael Holst
602 * ***************************************************************************
603 */
Bmat_IJA(Bmat * thee,int p,int q)604 VPUBLIC int *Bmat_IJA(Bmat *thee, int p, int q)
605 {
606 return Mat_IJA( thee->AD[p][q] );
607 }
608
609 /*
610 * ***************************************************************************
611 * Routine: Bmat_IA
612 *
613 * Purpose: Return the integer structure IA.
614 *
615 * Author: Michael Holst
616 * ***************************************************************************
617 */
Bmat_IA(Bmat * thee,int p,int q)618 VPUBLIC int *Bmat_IA(Bmat *thee, int p, int q)
619 {
620 return Mat_IA( thee->AD[p][q] );
621 }
622
623 /*
624 * ***************************************************************************
625 * Routine: Bmat_JA
626 *
627 * Purpose: Return the integer structure JA.
628 *
629 * Author: Michael Holst
630 * ***************************************************************************
631 */
Bmat_JA(Bmat * thee,int p,int q)632 VPUBLIC int *Bmat_JA(Bmat *thee, int p, int q)
633 {
634 return Mat_JA( thee->AD[p][q] );
635 }
636
637 /*
638 * ***************************************************************************
639 * Routine: Bmat_A
640 *
641 * Purpose: Return the real structure A.
642 *
643 * Author: Michael Holst
644 * ***************************************************************************
645 */
Bmat_A(Bmat * thee,int p,int q)646 VPUBLIC double *Bmat_A(Bmat *thee, int p, int q)
647 {
648 return Mat_A( thee->AD[p][q] );
649 }
650
651 /*
652 * ***************************************************************************
653 * Routine: Bmat_diag
654 *
655 * Purpose: Return the diagonal of A (DRC only).
656 *
657 * Author: Michael Holst
658 * ***************************************************************************
659 */
Bmat_diag(Bmat * thee,int p,int q)660 VPUBLIC double *Bmat_diag(Bmat *thee, int p, int q)
661 {
662 return Mat_diag( thee->AD[p][q] );
663 }
664
665 /*
666 * ***************************************************************************
667 * Routine: Bmat_offU
668 *
669 * Purpose: Return the strict upper triangle of A (DRC only).
670 *
671 * Author: Michael Holst
672 * ***************************************************************************
673 */
Bmat_offU(Bmat * thee,int p,int q)674 VPUBLIC double *Bmat_offU(Bmat *thee, int p, int q)
675 {
676 return Mat_offU( thee->AD[p][q] );
677 }
678
679 /*
680 * ***************************************************************************
681 * Routine: Bmat_offL
682 *
683 * Purpose: Return the strict lower triangle of A (DRC only).
684 *
685 * Author: Michael Holst
686 * ***************************************************************************
687 */
Bmat_offL(Bmat * thee,int p,int q)688 VPUBLIC double *Bmat_offL(Bmat *thee, int p, int q)
689 {
690 return Mat_offL( thee->AD[p][q] );
691 }
692
693 /*
694 * ***************************************************************************
695 * Routine: Bmat_print
696 *
697 * Purpose: Print the prolongation matrix blocks.
698 *
699 * Author: Michael Holst
700 * ***************************************************************************
701 */
Bmat_print(Bmat * thee)702 VPUBLIC void Bmat_print(Bmat *thee)
703 {
704 int p,q;
705
706 for (p=0; p<thee->numB; p++) {
707 for (q=0; q<thee->numB; q++) {
708 Mat_print(thee->AD[p][q]);
709 }
710 }
711 }
712
713 /*
714 * ***************************************************************************
715 * Routine: Bmat_printSp
716 *
717 * Purpose: Print the prolongation matrix in MATLAB sparse form.
718 *
719 * Notes: flag==0 ==> write
720 * flag==1 ==> append
721 *
722 * Author: Michael Holst
723 * ***************************************************************************
724 */
Bmat_printSp(Bmat * thee,char * fname,int pflag)725 VPUBLIC void Bmat_printSp(Bmat *thee, char *fname, int pflag)
726 {
727 char rn[80];
728 FILE *fp;
729 int p, q, numRT, numCT;
730
731 strncpy(rn,"Bmat_printSp:",80);
732
733 numRT = Bmat_numRT( thee );
734 numCT = Bmat_numCT( thee );
735
736 if (pflag == 0) {
737 fp=fopen(fname,"w");
738 } else if (pflag == 1) {
739 fp=fopen(fname,"a");
740 } else {
741 fp=VNULL;
742 }
743
744 if (fp==VNULL) {
745 Vnm_print(2, "%s problem opening file <%s>\n", rn, fname);
746 return;
747 }
748
749 /* print the matrix in matlab sparse format */
750 fprintf(fp,"%% %s matrix <%s> [dim=(%dx%d)]\n",
751 rn, thee->name, numRT, numCT);
752 fprintf(fp,"%% ----------------------------------------\n");
753 fprintf(fp, "fprintf('Defining the blocks of %s..');\n\n", thee->name);
754 fclose(fp);
755
756 for (p=0; p<thee->numB; p++) {
757 for (q=0; q<thee->numB; q++) {
758 if( !(thee->mirror[p][q]) ) {
759 Mat_printSp(thee->AD[p][q],fname,1);
760 }
761 }
762 }
763
764 fp=fopen(fname,"a");
765
766 if (fp==VNULL) {
767 Vnm_print(2, "%s problem opening file <%s>\n", rn, fname);
768 return;
769 }
770
771 fprintf(fp,"%% ----------------------------------------\n");
772 fprintf(fp,"%% Matlab code to generate the block matrix.\n");
773 fprintf(fp,"%% ----------------------------------------\n");
774 fprintf(fp,"%s = [\n",thee->name);
775 for (p=0; p<thee->numB; p++) {
776 for (q=0; q<thee->numB; q++) {
777 if( !thee->mirror[p][q] ) {
778 fprintf(fp," %s ",(thee->AD[p][q])->name);
779 } else {
780 fprintf(fp," %s'",(thee->AD[q][p])->name);
781 }
782 }
783 fprintf(fp,"\n");
784 }
785 fprintf(fp, "];\n\n");
786 fprintf(fp,"%% ----------------------------------------\n");
787
788 /* close file and return */
789 fclose(fp);
790 }
791
792 /*
793 * ***************************************************************************
794 * Routine: Bmat_printNoD
795 *
796 * Purpose: Print the matrix as a DENSE matrix in MATLAB format,
797 * but first zero out any rows/cols corresponding to
798 * Dirichlet boundary points.
799 *
800 * Notes: This routine is useful for e.g. checking that Galerkin
801 * conditions hold for stiffness matrices. Removing the
802 * dirichlet equations is crucial; otherwise the Galerkin
803 * condition cannot hold. Note that the matrix (and the
804 * Galerkin coarse matrix) are then of course singular.
805 *
806 * Author: Michael Holst
807 * ***************************************************************************
808 */
Bmat_printNoD(Bmat * thee)809 VPUBLIC void Bmat_printNoD(Bmat *thee)
810 {
811 int p,q;
812
813 for (p=0; p<thee->numB; p++) {
814 for (q=0; q<thee->numB; q++) {
815 Mat_printNoD(thee->AD[p][q]);
816 }
817 }
818 }
819
820 /*
821 * ***************************************************************************
822 * Routine: Bmat_printSpNoD
823 *
824 * Purpose: Print the matrix as a SPARSE matrix in MATLAB format,
825 * but first zero out any rows/cols corresponding to
826 * Dirichlet boundary points.
827 *
828 * Notes: This routine is useful for e.g. checking that Galerkin
829 * conditions hold for stiffness matrices. Removing the
830 * dirichlet equations is crucial; otherwise the Galerkin
831 * condition cannot hold. Note that the matrix (and the
832 * Galerkin coarse matrix) are then of course singular.
833 *
834 * flag==0 ==> write
835 * flag==1 ==> append
836 *
837 * Author: Michael Holst
838 * ***************************************************************************
839 */
Bmat_printSpNoD(Bmat * thee,char * fname,int pflag)840 VPUBLIC void Bmat_printSpNoD(Bmat *thee, char *fname, int pflag)
841 {
842 char rn[80];
843 FILE *fp;
844 int p, q, numRT, numCT;
845
846 strncpy(rn,"Bmat_printSpNoD:",80);
847
848 numRT = Bmat_numRT( thee );
849 numCT = Bmat_numCT( thee );
850
851 if (pflag == 0) {
852 fp=fopen(fname,"w");
853 } else if (pflag == 1) {
854 fp=fopen(fname,"a");
855 } else {
856 fp=VNULL;
857 }
858
859 if (fp==VNULL) {
860 Vnm_print(2, "%s problem opening file <%s>\n", rn, fname);
861 return;
862 }
863
864 /* print the matrix in matlab sparse format */
865 fprintf(fp,"%% %s matrix <%s> [dim=(%dx%d)]\n",
866 rn, thee->name, numRT, numCT);
867 fprintf(fp,"%% ----------------------------------------\n");
868 fprintf(fp, "fprintf('Defining the blocks of %s..');\n\n", thee->name);
869 fclose(fp);
870
871 for (p=0; p<thee->numB; p++) {
872 for (q=0; q<thee->numB; q++) {
873 Mat_printSpNoD(thee->AD[p][q],fname,1);
874 }
875 }
876
877 fp=fopen(fname,"a");
878
879 if (fp==VNULL) {
880 Vnm_print(2, "%s problem opening file <%s>\n", rn, fname);
881 return;
882 }
883
884 fprintf(fp,"%% ----------------------------------------\n");
885 fprintf(fp,"%% Matlab code to generate the block matrix.\n");
886 fprintf(fp,"%% ----------------------------------------\n");
887 fprintf(fp,"%s = [\n",thee->name);
888 for (p=0; p<thee->numB; p++) {
889 for (q=0; q<thee->numB; q++) {
890 fprintf(fp," %s",(thee->AD[p][q])->name);
891 }
892 fprintf(fp,"\n");
893 }
894 fprintf(fp, "];\n\n");
895 fprintf(fp,"%% ----------------------------------------\n");
896
897 /* close file and return */
898 fclose(fp);
899
900 }
901
902 /*
903 * ***************************************************************************
904 * Routine: Bmat_zero
905 *
906 * Purpose: Clear the floating point storage for the sparse matrix.
907 * Also clear any sparse factorization storage.
908 *
909 * Notes: This is basically done in preparation for an accumulation as
910 * part of a matrix assembly, and before a new sparse factorization.
911 *
912 * Author: Michael Holst
913 * ***************************************************************************
914 */
Bmat_zero(Bmat * thee)915 VPUBLIC void Bmat_zero(Bmat *thee)
916 {
917 int p, q;
918
919 /* zero out all of the floating point storage */
920 for (p=0; p<thee->numB; p++) {
921 for (q=0; q<thee->numB; q++) {
922 if ( !thee->mirror[p][q] ) {
923 Mat_zero(thee->AD[p][q]);
924 }
925 }
926 }
927 }
928
929 /*
930 * ***************************************************************************
931 * Routine: Bmat_diri
932 *
933 * Purpose: Setup the dirichlet equations.
934 *
935 * Author: Michael Holst
936 * ***************************************************************************
937 */
Bmat_diri(Bmat * thee)938 VPUBLIC void Bmat_diri(Bmat *thee)
939 {
940 int p;
941
942 /* stick identity on the boundary rows/cols of diagonal blocks */
943 for (p=0; p<thee->numB; p++) {
944 Mat_diagBRC(thee->AD[p][p]);
945 }
946 }
947
948 /*
949 * ***************************************************************************
950 * Routine: Bmat_set
951 *
952 * Purpose: Set the (i,j)-th entry of the (p,q)-th block to <val>
953 *
954 * Author: Michael Holst
955 * ***************************************************************************
956 */
Bmat_set(Bmat * thee,int p,int q,int i,int j,double val)957 VPUBLIC void Bmat_set(Bmat *thee, int p, int q, int i, int j, double val)
958 {
959 if ( !thee->mirror[p][q] ) {
960 Mat_set(thee->AD[p][q],i,j,val);
961 } else {
962 /* we should NEVER get to here */
963 VASSERT(0);
964 }
965 }
966
967 /*
968 * ***************************************************************************
969 * Routine: Bmat_addTo
970 *
971 * Purpose: Contribute <val> to the (i,j)-th entry of the (p,q)-th block.
972 *
973 * Author: Michael Holst
974 * ***************************************************************************
975 */
Bmat_addTo(Bmat * thee,int p,int q,int i,int j,double val)976 VPUBLIC void Bmat_addTo(Bmat *thee, int p, int q, int i, int j, double val)
977 {
978 if ( !thee->mirror[p][q] ) {
979 Mat_addTo(thee->AD[p][q],i,j,val);
980 } else {
981 /* we should NEVER get to here */
982 VASSERT(0);
983 }
984 }
985
986 /*
987 * ***************************************************************************
988 * Routine: Bmat_galerkin
989 *
990 * Purpose: Enforce the Galerkin conditions algebraically.
991 *
992 * Author: Michael Holst
993 * ***************************************************************************
994 */
Bmat_galerkin(Bmat * thee,Bmat * rmat,Bmat * amat,Bmat * pmat)995 VPUBLIC void Bmat_galerkin(Bmat *thee, Bmat *rmat, Bmat *amat, Bmat *pmat)
996 {
997 int p, q, numR, numC;
998 char bname[15];
999
1000 #if 0
1001 Vnm_tstart(30, "triple matrix product");
1002 Vnm_print(0,"Bmat_galerkin: triple matrix product..");
1003 #endif
1004
1005 /* go through blocks and form B=r*A*p */
1006 for (p=0; p<thee->numB; p++) {
1007 for (q=0; q<thee->numB; q++) {
1008
1009 /* form the product: B=R*A*P */
1010 if ( !Bmat_mirror(thee,p,q) ) {
1011
1012 /* hold onto row/col sizes */
1013 numR = Bmat_numR(thee,p,q);
1014 numC = Bmat_numR(thee,p,q);
1015
1016 /* reinitialize the block */
1017 Mat_dtor(&thee->AD[p][q]);
1018 sprintf(bname, "%s_%d%d", thee->name, p, q);
1019 thee->AD[p][q] = Mat_ctor(thee->vmem, bname, numR, numC);
1020
1021 /* don't forget the possible mirror!!! */
1022 if ( Bmat_mirror(thee,q,p) ) {
1023 thee->AD[q][p] = thee->AD[p][q];
1024 }
1025
1026 /* now form the triple matrix product */
1027 Mat_galerkin(thee->AD[p][q],
1028 rmat->AD[p][p], amat->AD[p][q], pmat->AD[q][q]);
1029 }
1030 }
1031 }
1032
1033 /*
1034 * Now, insert dirichlet identity equations due to singular P.
1035 *
1036 * (The new matrix is singular otherwise, since the prolongation matrix
1037 * has zero COLUMNS at coarse dirichlet points and zero ROWS at fine
1038 * dirichlet points, with the reverse situation for the restriction.)
1039 */
1040 Bmat_diri( thee );
1041
1042 /* Set the coarse and fine pointers within the matrices. */
1043 thee->fine = amat;
1044 amat->coarse = thee;
1045
1046 #if 0
1047 Vnm_print(0,"..done.\n");
1048 Vnm_tstop(30, "triple matrix product");
1049 #endif
1050 }
1051
1052 /*
1053 * ***************************************************************************
1054 * Routine: Bmat_sluDirect
1055 *
1056 * Purpose: Make a decision about whether or not a sparse direct solver
1057 * should be used in place of an iterative solver, based on the
1058 * size of the system.
1059 *
1060 * Notes: This is obviously heuristic in nature; in general the cutoff
1061 * size where iterative methods start to win is smaller in 3D.
1062 *
1063 * Author: Michael Holst
1064 * ***************************************************************************
1065 */
Bmat_sluDirect(Bmat * thee)1066 VPUBLIC int Bmat_sluDirect(Bmat *thee)
1067 {
1068 int ival = 0;
1069 /* the decision about when sparse direct will be the best solver */
1070 if ( Bmat_numRT(thee) < SPARSE_CUTOFF ) ival = 1;
1071 return ival;
1072 }
1073
1074 /*
1075 * ***************************************************************************
1076 * Routine: Bmat_sluCreate
1077 *
1078 * Purpose: Create the global matrix from the blocks in the modified
1079 * ROW or COL storage format. This is useful for preparing a
1080 * single global matrix for input to e.g. a sparse direct solver.
1081 *
1082 * Author: Michael Holst
1083 * ***************************************************************************
1084 */
Bmat_sluCreate(Bmat * thee)1085 VPUBLIC void Bmat_sluCreate(Bmat *thee)
1086 {
1087 int i, j, k, p, q, count, istart, jstart;
1088 Vset *lnk;
1089 LinkA *mt;
1090 Mat *blk, *gmat;
1091 MATmirror mirror;
1092 MATformat format;
1093 double *diag, *offU, *offL;
1094
1095 /* initialize the global matrix datastructure */
1096 gmat = Mat_ctor(thee->vmem, "AG", Bmat_numRT(thee), Bmat_numCT(thee));
1097 Mat_initStructure(gmat, SLU_FORMAT, ISNOT_SYM, Bmat_numZT(thee),
1098 VNULL, VNULL);
1099
1100 /* initialize/clear the dynamic array */
1101 lnk = Vset_ctor( thee->vmem, "lnk", sizeof( LinkA ), VMAX_OBJECTS, 0 );
1102 Vset_reset( lnk );
1103
1104 /* create an empty entry to start each row of global matrix */
1105 for (i=0; i<Bmat_numRT(thee); i++) {
1106 mt=(LinkA*)Vset_create(lnk);
1107 mt->idx = -1;
1108 mt->val = 0.;
1109 mt->next = VNULL;
1110 }
1111
1112 /* now get the COL structure of the matrix */
1113 count = 0;
1114
1115 istart = 0;
1116 for (p=0; p<thee->numB; p++) {
1117
1118 jstart = 0;
1119 for (q=0; q<thee->numB; q++) {
1120
1121 mirror = thee->mirror[p][q];
1122 blk = thee->AD[p][q];
1123 format = Mat_format(blk);
1124 diag = blk->diag;
1125 offU = blk->offU;
1126 offL = blk->offL;
1127 if (mirror) {
1128 if (format == ROW_FORMAT) {
1129 format = COL_FORMAT;
1130 } else if (format == COL_FORMAT) {
1131 format = ROW_FORMAT;
1132 } else {
1133 VASSERT(0);
1134 }
1135 }
1136
1137 for (j=0; j<Mat_numR(blk); j++) {
1138 if (format == DRC_FORMAT) {
1139 i = j;
1140 mContrib(lnk,0,&count,jstart+j,istart+i,diag[i]);
1141 for (k=blk->IA[j]; k<blk->IA[j+1]; k++) {
1142 i = blk->JA[k];
1143 mContrib(lnk,0,&count,jstart+j,istart+i,offL[k]);
1144 /* still in "istart/jstart" block; */
1145 /* but roles of "i" and "j" are reversed */
1146 mContrib(lnk,0,&count,jstart+i,istart+j,offU[k]);
1147 }
1148 } else if (format == ROW_FORMAT) {
1149 for (k=blk->IA[j]; k<blk->IA[j+1]; k++) {
1150 i = blk->JA[k];
1151 mContrib(lnk,0,&count,jstart+i,istart+j,blk->A[k]);
1152 }
1153 } else if (format == COL_FORMAT) {
1154 for (k=blk->IA[j]; k<blk->IA[j+1]; k++) {
1155 i = blk->JA[k];
1156 mContrib(lnk,0,&count,jstart+j,istart+i,blk->A[k]);
1157 }
1158 } else { VASSERT(0); }
1159
1160 }
1161 /* increment index for the column block */
1162 jstart += Bmat_numC(thee,q,q);
1163
1164 }
1165 /* increment index for the row block */
1166 istart += Bmat_numR(thee,p,p);
1167
1168 }
1169
1170 /* now create the matrix from the temporary structure */
1171 count = 0;
1172 gmat->IA[0] = 0;
1173 for (i=0; i<Bmat_numRT(thee); i++) {
1174 for (mt=(LinkA*)Vset_access(lnk,i); mt!=VNULL; mt=mt->next) {
1175 if (mt->idx >= 0) {
1176 gmat->JA[count] = mt->idx;
1177 gmat->A[count] = mt->val;
1178 count++;
1179 }
1180 }
1181 gmat->IA[i+1] = count;
1182 }
1183 thee->AG = gmat;
1184
1185 /* clear/destroy the dynamic array */
1186 Vset_reset( lnk );
1187 Vset_dtor( &lnk );
1188 }
1189
1190 /*
1191 * ***************************************************************************
1192 * Routine: Bmat_sluFactor
1193 *
1194 * Purpose: Create the sparse LU factors for global matrix.
1195 *
1196 * Author: Michael Holst
1197 * ***************************************************************************
1198 */
Bmat_sluFactor(Bmat * thee)1199 VPUBLIC int Bmat_sluFactor(Bmat *thee)
1200 {
1201 if ( thee->AG == VNULL ) {
1202 Bmat_sluCreate(thee);
1203 }
1204
1205 if ( Mat_state(thee->AG) == FACTORED_STATE ) {
1206 return 1;
1207 } else {
1208 return Mat_sluFactor(thee->AG);
1209 }
1210 }
1211
1212 /*
1213 * ***************************************************************************
1214 * Routine: Bmat_sluSolve
1215 *
1216 * Purpose: Forward/backward solve using sparse LU factors of global matrix.
1217 *
1218 * Notes: This requires that Bmat_sluFactor has been previously called.
1219 *
1220 * Author: Michael Holst
1221 * ***************************************************************************
1222 */
Bmat_sluSolve(Bmat * thee,int key,double * f,double * u)1223 VPUBLIC int Bmat_sluSolve(Bmat *thee, int key, double *f, double *u)
1224 {
1225 return Mat_sluSolve(thee->AG, key, f, u);
1226 }
1227
1228 /*
1229 * ***************************************************************************
1230 * Routine: Bmat_sluDestroy
1231 *
1232 * Purpose: Destroy the sparse LU factors for the system matrix.
1233 *
1234 * Notes: This frees our <ia,ja,a> storage, and also the internal
1235 * storage that was malloc'd by the sparse direct library.
1236 *
1237 * Author: Michael Holst
1238 * ***************************************************************************
1239 */
Bmat_sluDestroy(Bmat * thee)1240 VPUBLIC void Bmat_sluDestroy(Bmat *thee)
1241 {
1242 Mat_sluDestroy(thee->AG);
1243 }
1244
1245 /*
1246 * ***************************************************************************
1247 * Routine: Bmat_memChk
1248 *
1249 * Purpose: Print the exact current malloc usage.
1250 *
1251 * Author: Michael Holst
1252 * ***************************************************************************
1253 */
Bmat_memChk(Bmat * thee)1254 VPUBLIC void Bmat_memChk(Bmat *thee)
1255 {
1256 if (thee->iMadeVmem) Vmem_print(thee->vmem);
1257 }
1258
1259 /*
1260 * ***************************************************************************
1261 * Routine: Bmat_clone
1262 *
1263 * Purpose: Construct a clone of a Bmat with the same structure.
1264 *
1265 * Author: Stephen Bond
1266 * ***************************************************************************
1267 */
Bmat_clone(Vmem * vmem,char * name,Bmat * X)1268 VPUBLIC Bmat *Bmat_clone(Vmem *vmem, char *name, Bmat *X)
1269 {
1270
1271 int p, q, numR[MAXV], numC[MAXV];
1272 int numB;
1273 MATmirror pmir[MAXV][MAXV];
1274
1275 numB = Bmat_numB(X);
1276
1277 for (p=0; p<numB; p++) {
1278 numR[p] = Bmat_numR( X, p, p );
1279 numC[p] = Bmat_numC( X, p, p );
1280 for (q=0; q<numB; q++) {
1281 pmir[p][q] = Bmat_mirror( X, p, q );
1282 }
1283 }
1284
1285 return Bmat_ctor(vmem, name, numB, numR, numC, pmir );
1286
1287 }
1288
1289 /*
1290 * ***************************************************************************
1291 * Routine: Bmat_copy
1292 *
1293 * Purpose: Copy a block matrix.
1294 *
1295 * Author: Stephen Bond
1296 * ***************************************************************************
1297 */
Bmat_copy(Bmat * Y,Bmat * X)1298 VPUBLIC void Bmat_copy(Bmat *Y, Bmat *X)
1299 {
1300
1301 int p, q;
1302 int numB;
1303
1304 numB = Bmat_numB(X);
1305
1306 VASSERT( Bmat_numB(Y) == numB );
1307
1308 for (p=0; p<numB; p++) {
1309 for (q=0; q<numB; q++) {
1310 if ( !Y->mirror[p][q] ) {
1311 Mat_copy( Y->AD[p][q], X->AD[p][q] );
1312 } else {
1313 VASSERT( !Y->mirror[q][p] );
1314 VASSERT( Y->AD[q][p] != VNULL );
1315 VASSERT( Y->AD[p][q] == VNULL );
1316 Y->AD[p][q] = Y->AD[q][p];
1317 }
1318 }
1319 }
1320
1321 }
1322
1323 /*
1324 * ***************************************************************************
1325 * Routine: Bmat_squeezeBRC
1326 *
1327 * Purpose: Remove the boundary rows or columns from a block matrix.
1328 *
1329 * Author: Stephen Bond
1330 * ***************************************************************************
1331 */
Bmat_squeezeBRC(Bmat * thee,int key)1332 VPUBLIC void Bmat_squeezeBRC(Bmat *thee, int key)
1333 {
1334
1335 int p, q;
1336 int numB;
1337
1338 numB = Bmat_numB(thee);
1339
1340 /* this only makes sense for prolongation matrices */
1341
1342 for (p=0; p<numB; p++) {
1343 for (q=p+1; q<numB; q++) {
1344 VASSERT( Bmat_mirror(thee,p,q) == ISNOT_MIRROR );
1345 VASSERT( Bmat_mirror(thee,q,p) == ISNOT_MIRROR );
1346 VASSERT( Bmat_format(thee,p,q) == ZERO_FORMAT );
1347 VASSERT( Bmat_format(thee,q,p) == ZERO_FORMAT );
1348 }
1349 Mat_squeezeBRC( thee->AD[p][p], key );
1350 }
1351 }
1352
1353 /*
1354 * ***************************************************************************
1355 * Routine: Bmat_copy2
1356 *
1357 * Purpose: Copy a block matrix.
1358 *
1359 * Author: Stephen Bond
1360 * ***************************************************************************
1361 */
Bmat_copy2(Bmat * Y,Bmat * X)1362 VPUBLIC void Bmat_copy2(Bmat *Y, Bmat *X)
1363 {
1364
1365 int p, q;
1366 int numB;
1367
1368 numB = Bmat_numB(X);
1369
1370 VASSERT( Bmat_numB(Y) == numB );
1371
1372 for (p=0; p<numB; p++) {
1373 for (q=0; q<numB; q++) {
1374 if ( !Y->mirror[p][q] ) {
1375 Mat_copy2( Y->AD[p][q], X->AD[p][q] );
1376 } else {
1377 VASSERT( !Y->mirror[q][p] );
1378 VASSERT( Y->AD[q][p] == Y->AD[q][p] );
1379 }
1380 }
1381 }
1382
1383 }
1384
1385 /*
1386 * ***************************************************************************
1387 * Routine: Bmat_axpy
1388 *
1389 * Purpose: Scalar times a Bmat plus a Bmat: Y += val*X.
1390 *
1391 * Notes: The function of this routine can be controlled using "key"
1392 *
1393 * Author: Stephen Bond
1394 * ***************************************************************************
1395 */
Bmat_axpy(Bmat * Y,Bmat * X,double val,int key)1396 VPUBLIC void Bmat_axpy(Bmat *Y, Bmat *X, double val, int key)
1397 {
1398
1399 int p, q;
1400 int numB;
1401
1402 numB = X->numB;
1403
1404 /* General checks of the block sizes of Y and X */
1405 VASSERT( Bmat_numB(Y) == numB );
1406
1407 for (p=0; p<numB; p++) {
1408 for (q=0; q<numB; q++) {
1409 if ( !Y->mirror[p][q] ) {
1410 Mat_axpy( Y->AD[p][q], X->AD[p][q], val, key );
1411 }
1412 }
1413 }
1414
1415 }
1416
1417