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: matln.c,v 1.6 2010/08/12 05:18:28 fetk Exp $"
21 * ***************************************************************************
22 */
23
24 /*
25 * ***************************************************************************
26 * File: matln.c
27 *
28 * Purpose: Class Mat: linked list data structures.
29 *
30 * Author: Stephen Bond and Michael Holst
31 * ***************************************************************************
32 */
33
34 #include "mat_p.h"
35
36 VEMBED(rcsid="$Id: matln.c,v 1.6 2010/08/12 05:18:28 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_initStructureLN
55 *
56 * Purpose: Initialize the nonzero structure given structure information.
57 *
58 * Notes: This routine actually does the storage creation for all
59 * internal Vset, link arrays, and link pointer arrays.
60 *
61 * Author: Stephen Bond and Michael Holst
62 * ***************************************************************************
63 */
Mat_initStructureLN(Mat * thee,MATformat frmt,MATsym sym)64 VPUBLIC void Mat_initStructureLN(Mat *thee, MATformat frmt, MATsym sym)
65 {
66 int i, sizeA;
67 LinkA *mt;
68 LinkRCS *mtX;
69
70 /* use an upperbound for the maximum size of the Vset */
71 sizeA = (thee->numR)*LN_MAX_ENTRIES_PER_ROW;
72
73 /* build framework from input */
74 thee->format = frmt;
75 thee->sym = sym;
76
77 /* initialize the state */
78 thee->state = ZERO_STATE;
79
80 /* Format-dependent consistency checks and setup */
81 switch (thee->format) {
82
83 /* RLN consistency checks and setup */
84 case RLN_FORMAT:
85
86 /* RLN valid for nonsymmetry */
87 VASSERT( thee->sym==ISNOT_SYM );
88
89 thee->lnkU = Vset_ctor( thee->vmem, "U", sizeof(LinkA), sizeA, 0 );
90
91 /* create an empty entry to start each row of the matrix */
92 for (i=0; i<(thee->numR); i++) {
93 mt=(LinkA*)Vset_create(thee->lnkU);
94 mt->idx = -1;
95 mt->val = 0.;
96 mt->next = VNULL;
97 }
98
99 thee->numO = 0;
100
101 break;
102
103 /* CLN consistency checks and setup */
104 case CLN_FORMAT:
105
106 /* RLN valid for nonsymmetry */
107 VASSERT( thee->sym==ISNOT_SYM );
108
109 thee->lnkL = Vset_ctor( thee->vmem, "L", sizeof(LinkA), sizeA, 0 );
110
111 /* create an empty entry to start each column of the matrix */
112 for (i=0; i<(thee->numC); i++) {
113 mt=(LinkA*)Vset_create(thee->lnkL);
114 mt->idx = -1;
115 mt->val = 0.;
116 mt->next = VNULL;
117 }
118
119 thee->numO = 0;
120
121 break;
122
123 /* XLN consistency checks and setup */
124 case XLN_FORMAT:
125
126 switch( sym ) {
127 case IS_SYM:
128 thee->lnkU = Vset_ctor( thee->vmem,"X",sizeof(LinkRCS),sizeA,0 );
129 thee->numA = thee->numR;
130 thee->xln = Vmem_malloc( thee->vmem,thee->numA,sizeof(LinkRCS) );
131
132 /* create a zero entries to start the diagonal of the matrix */
133 for (i=0; i<(thee->numR); i++) {
134 mtX = &( ((LinkRCS*) (thee->xln))[i] );
135 mtX->idxT = i;
136 mtX->idx = i;
137 mtX->nxtT = VNULL;
138 mtX->next = VNULL;
139 mtX->val = 0.;
140 }
141 break;
142 case STRUC_SYM:
143 thee->lnkU = Vset_ctor( thee->vmem,"X",sizeof(LinkRCN),sizeA,0 );
144 thee->numA = thee->numR;
145 thee->xln = Vmem_malloc( thee->vmem,thee->numA,sizeof(LinkRCS) );
146
147 /* create a zero entries to start the diagonal of the matrix */
148 for (i=0; i<(thee->numR); i++) {
149 mtX = &( ((LinkRCS*) (thee->xln))[i] );
150 mtX->idxT = i;
151 mtX->idx = i;
152 mtX->nxtT = VNULL;
153 mtX->next = VNULL;
154 mtX->val = 0.;
155 }
156 break;
157 case ISNOT_SYM:
158 thee->lnkU = Vset_ctor( thee->vmem,"X",sizeof(LinkRCS),sizeA,0 );
159 thee->numA = thee->numR + thee->numC;
160 thee->xln = Vmem_malloc( thee->vmem,thee->numA,sizeof(LinkRC*) );
161 thee->xlnt = (void*) &( ((LinkRC**) (thee->xln))[thee->numR] );
162 break;
163 default:
164 VASSERT(0);
165 break;
166 }
167
168 thee->numO = 0;
169
170 break;
171 default:
172 VASSERT(0);
173 break;
174 }
175 }
176
177 /*
178 * ***************************************************************************
179 * Routine: Mat_killStructureLN
180 *
181 * Purpose: Kill the nonzero structure and structure information.
182 *
183 * Notes: This routine does the reverse of Mat_initStructureLN
184 * (or Mat_copyStructureLN). It leaves only the information
185 * about the number of blocks, number of rows, and number of
186 * columns per block. I.e., what is left is only what was
187 * present after the initial call to Mat_ctor.
188 *
189 * Author: Stephen Bond and Michael Holst
190 * ***************************************************************************
191 */
Mat_killStructureLN(Mat * thee)192 VPUBLIC void Mat_killStructureLN(Mat *thee)
193 {
194 if (thee != VNULL) {
195 if (thee->state != NULL_STATE) {
196
197 /* Format-dependent destruction */
198 switch (thee->format) {
199
200 case RLN_FORMAT:
201 if (thee->lnkU != VNULL) {
202 Vset_dtor( &(thee->lnkU) );
203 }
204 break;
205
206 case CLN_FORMAT:
207 if (thee->lnkL != VNULL) {
208 Vset_dtor( &(thee->lnkL) );
209 }
210 break;
211
212 case XLN_FORMAT:
213 if (thee->lnkU != VNULL)
214 Vset_dtor( &(thee->lnkU) );
215
216 if (thee->xln != VNULL) {
217 if ( Mat_sym(thee) == ISNOT_SYM ) {
218 Vmem_free( thee->vmem, thee->numA, sizeof(LinkRC*),
219 (void**)&(thee->xln));
220 } else {
221 Vmem_free( thee->vmem, thee->numA, sizeof(LinkRCS),
222 (void**)&(thee->xln));
223 }
224 }
225
226 thee->numA = 0;
227 break;
228
229 default:
230 VASSERT(0);
231 break;
232 }
233
234 /* initialize the state */
235 thee->state = NULL_STATE;
236
237 }
238 }
239 }
240
241 /*
242 * ***************************************************************************
243 * Routine: Mat_accessXLN
244 *
245 * Purpose: Access the first element in the ROW or COL of an XLN.
246 *
247 * Notes: key==0 ==> return a ROW pointer
248 * key==1 ==> return a COL pointer
249 *
250 * In the symmetric cases, we are just returning a pointer
251 * to the diagonal. In the nonsymmetric case, we are returning
252 * a pointer to the first element in the row or column. In
253 * the symmetric cases the columns are linked in reverse index
254 * ordering to save the storage of an additional pointer array.
255 *
256 * Author: Stephen Bond
257 * ***************************************************************************
258 */
Mat_accessXLN(Mat * thee,int idx,int key)259 VPUBLIC LinkRC* Mat_accessXLN(Mat *thee, int idx, int key)
260 {
261 if ( thee->sym != ISNOT_SYM /* && ( key == 0 || key == 1 ) */ ) {
262
263 return (LinkRC*) &(((LinkRCS*) (thee->xln))[idx]);
264
265 } else if ( key == 0 /* && thee->sym == ISNOT_SYM */ ) {
266
267 return ((LinkRC**) (thee->xln))[idx];
268
269 } else /* if (key == 1 && thee->sym == ISNOT_SYM ) */ {
270
271 return ((LinkRC**) (thee->xlnt))[idx];
272 }
273 }
274
275 /*
276 * ***************************************************************************
277 * Routine: Mat_contribXLN
278 *
279 * Purpose: Set or add a value to a doubly linked matrix entry array.
280 *
281 * Notes: key==0 ==> Set the value
282 * key==1 ==> Add the value
283 *
284 * This is a doubly linked variant of mContrib.
285 *
286 * Author: Stephen Bond
287 * ***************************************************************************
288 */
Mat_contribXLN(Mat * thee,int key,int i,int j,double val)289 VPUBLIC void Mat_contribXLN(Mat *thee, int key, int i, int j, double val)
290 {
291 switch (thee->sym) {
292 case IS_SYM:
293 Mat_contribSYMXLN( thee, key, i, j, val );
294 break;
295 case STRUC_SYM:
296 Mat_contribSSYMXLN( thee, key, i, j, val );
297 break;
298 case ISNOT_SYM:
299 Mat_contribNSYMXLN( thee, key, i, j, val );
300 break;
301 default:
302 VASSERT( 0 );
303 }
304 }
305
306 /*
307 * ***************************************************************************
308 * Routine: Mat_contribNSYMXLN
309 *
310 * Purpose: Set or add a value to a NOTSYM XLN matrix.
311 *
312 * Notes: key==0 ==> Set the value
313 * key==1 ==> Add the value
314 *
315 * This is a doubly linked variant of mContrib.
316 *
317 * Author: Stephen Bond
318 * ***************************************************************************
319 */
Mat_contribNSYMXLN(Mat * thee,int key,int i,int j,double val)320 VPUBLIC void Mat_contribNSYMXLN(Mat *thee, int key, int i, int j, double val)
321 {
322 int flag = 0, done = 0;
323 LinkRC *curr = VNULL, *prev = VNULL, *mnew = VNULL;
324
325 VASSERT( thee->sym == ISNOT_SYM );
326
327 /* STEP 1: Traverse rowwise to the location for insert/set/contrib */
328 curr = ((LinkRC**) (thee->xln))[i];
329 prev = VNULL;
330 done = 0;
331 flag = 0;
332 while (!done) {
333 if (curr == VNULL) {
334 /* first guy in row is "NULL", fill him with the info */
335 mnew = (LinkRC*)Vset_create(thee->lnkU);
336 mnew->idxT = i;
337 mnew->idx = j;
338 mnew->nextT = VNULL;
339 mnew->next = VNULL;
340 ((LinkRCS*) mnew)->val = val;
341 ((LinkRC**) (thee->xln))[i] = mnew;
342 thee->numO++;
343 done = 1;
344 flag = 1;
345 } else if (curr->idx > j) {
346 /* not in the list; insert new link between curr and prev */
347 mnew = (LinkRC*)Vset_create(thee->lnkU);
348 mnew->idxT = i;
349 mnew->idx = j;
350 mnew->nextT = VNULL;
351 mnew->next = VNULL;
352 ((LinkRCS*) mnew)->val = val;
353 if (prev == VNULL) {
354 ((LinkRC**) (thee->xln))[i] = mnew;
355 mnew->next = curr;
356 } else {
357 prev->next = mnew;
358 mnew->next = curr;
359 }
360 thee->numO++;
361 done = 1;
362 flag = 1;
363 } else if (curr->idx == j) {
364 /* found the position; just add the contribution */
365 mnew = (LinkRC*)curr;
366 if (key == 0)
367 ((LinkRCS*) mnew)->val = val;
368 else
369 ((LinkRCS*) mnew)->val += val;
370 done = 1;
371 flag = 0;
372 } else if (curr->next == VNULL) {
373 /* not found; no one left, append on the end */
374 mnew = (LinkRC*)Vset_create(thee->lnkU);
375 mnew->idxT = i;
376 mnew->idx = j;
377 ((LinkRCS*) mnew)->val = val;
378 mnew->next = VNULL;
379 curr->next = (LinkRC*)mnew;
380 thee->numO++;
381 done = 1;
382 flag = 1;
383 /* not found yet; still hope */
384 } else {
385 prev = curr;
386 curr = curr->next;
387 }
388 }
389
390 if (flag == 0) {
391 /* exit now if no new links were inserted */
392 return;
393 }
394
395 /* STEP 2: Traverse colwise to find the location for insertion */
396 curr = ((LinkRC**) (thee->xlnt))[j];
397 prev = VNULL;
398 done = 0;
399 while (!done) {
400 if (curr == VNULL) {
401 /* first guy in col is "NULL", set to point to new link */
402 ((LinkRC**) (thee->xlnt))[j] = mnew;
403 done = 1;
404 } else if (curr->idxT > i) {
405 /* not in the list; insert new link between curr and prev */
406 if (prev == VNULL) {
407 ((LinkRC**) (thee->xlnt))[j] = mnew;
408 mnew->nextT = curr;
409 } else {
410 prev->nextT = mnew;
411 mnew->nextT = curr;
412 }
413 done = 1;
414 } else if (curr->nextT == VNULL) {
415 /* not found; no one left; set the last pointer to new link */
416 curr->nextT = mnew;
417 done = 1;
418 /* not found yet; still hope */
419 } else {
420 prev = curr;
421 curr = curr->nextT;
422 }
423 }
424 }
425
426 /*
427 * ***************************************************************************
428 * Routine: Mat_contribSSYMXLN
429 *
430 * Purpose: Set or add a value to a STRUC_SYM XLN matrix.
431 *
432 * Notes: key==0 ==> Set the value
433 * key==1 ==> Add the value
434 *
435 * This is a doubly linked variant of mContrib.
436 *
437 * Author: Stephen Bond
438 * ***************************************************************************
439 */
Mat_contribSSYMXLN(Mat * thee,int key,int i,int j,double val)440 VPUBLIC void Mat_contribSSYMXLN(Mat *thee, int key, int i, int j, double val)
441 {
442 int ii = i, jj = j, flag = 0, done = 0, contribTYPE = 0;
443 LinkRC *curr = VNULL, *prev = VNULL, *mnew = VNULL;
444
445 VASSERT( thee->sym == STRUC_SYM );
446
447 if ( i > j ) {
448 ii = j;
449 jj = i;
450 contribTYPE = 1;
451 }
452
453 /* STEP 1: Traverse rowwise to the location for insert/set/contrib */
454 curr = (LinkRC*) &(((LinkRCS*) (thee->xln))[ii]);
455 prev = VNULL;
456 done = 0;
457 flag = 0;
458 while (!done) {
459 if (curr->idx > jj) {
460 /* not in the list; insert new link between curr and prev */
461 mnew = (LinkRC*)Vset_create(thee->lnkU);
462 mnew->idxT = ii;
463 mnew->idx = jj;
464 if ( contribTYPE == 0) {
465 ((LinkRCN*) mnew)->val = val;
466 ((LinkRCN*) mnew)->valT = 0.0;
467 } else {
468 ((LinkRCN*) mnew)->val = 0.0;
469 ((LinkRCN*) mnew)->valT = val;
470 }
471 mnew->nextT = VNULL;
472 mnew->next = VNULL;
473 prev->next = mnew;
474 mnew->next = curr;
475 thee->numO++;
476 done = 1;
477 flag = 1;
478 } else if (curr->idx == jj) {
479 /* found the position; just add the contribution */
480 mnew = (LinkRC*)curr;
481 if ( contribTYPE == 0 ) {
482 if (key == 0)
483 ((LinkRCS*) mnew)->val = val;
484 else
485 ((LinkRCS*) mnew)->val += val;
486 } else {
487 if (key == 0)
488 ((LinkRCN*) mnew)->valT = val;
489 else
490 ((LinkRCN*) mnew)->valT += val;
491 }
492 done = 1;
493 flag = 0;
494 } else if (curr->next == VNULL) {
495 /* not found; no one left, append on the end */
496 mnew = (LinkRC*)Vset_create(thee->lnkU);
497 mnew->idxT = ii;
498 mnew->idx = jj;
499 if ( contribTYPE == 0) {
500 ((LinkRCN*) mnew)->val = val;
501 ((LinkRCN*) mnew)->valT = 0.0;
502 } else {
503 ((LinkRCN*) mnew)->val = 0.0;
504 ((LinkRCN*) mnew)->valT = val;
505 }
506 mnew->next = VNULL;
507 curr->next = (LinkRC*)mnew;
508 thee->numO++;
509 done = 1;
510 flag = 1;
511 /* not found yet; still hope */
512 } else {
513 prev = curr;
514 curr = curr->next;
515 }
516 }
517
518 if (flag == 0) {
519 /* exit now if no new links were inserted */
520 return;
521 }
522
523 /* STEP 2: Traverse colwise to find the location for insertion */
524 curr = (LinkRC*) &(((LinkRCS*) (thee->xln))[jj]);
525 prev = VNULL;
526 done = 0;
527 while (!done) {
528 if (curr->idxT < ii) {
529 /* not in the list; insert new link between curr and prev */
530 prev->nextT = mnew;
531 mnew->nextT = curr;
532 done = 1;
533 } else if (curr->nextT == VNULL) {
534 /* not found; no one left; set the last pointer to new link */
535 curr->nextT = mnew;
536 done = 1;
537 /* not found yet; still hope */
538 } else {
539 prev = curr;
540 curr = curr->nextT;
541 }
542 }
543 }
544
545 /*
546 * ***************************************************************************
547 * Routine: Mat_contribSYMXLN
548 *
549 * Purpose: Set or add a value to a SYM XLN Matrix.
550 *
551 * Notes: key==0 ==> Set the value
552 * key==1 ==> Add the value
553 *
554 * This is a doubly linked variant of mContrib.
555 *
556 * Author: Stephen Bond
557 * ***************************************************************************
558 */
Mat_contribSYMXLN(Mat * thee,int key,int i,int j,double val)559 VPUBLIC void Mat_contribSYMXLN(Mat *thee, int key, int i, int j, double val)
560 {
561 int flag = 0, done = 0;
562 LinkRC *curr = VNULL, *prev = VNULL, *mnew = VNULL;
563
564 VASSERT( thee->sym == IS_SYM );
565 VASSERT( i <= j );
566
567 /* STEP 1: Traverse rowwise to the location for insert/set/contrib */
568 curr = (LinkRC*) &(((LinkRCS*) (thee->xln))[i]);
569 prev = VNULL;
570 done = 0;
571 flag = 0;
572 while (!done) {
573 if (curr->idx > j) {
574 /* not in the list; insert new link between curr and prev */
575 mnew = (LinkRC*)Vset_create(thee->lnkU);
576 mnew->idxT = i;
577 mnew->idx = j;
578 mnew->nextT = VNULL;
579 mnew->next = VNULL;
580 ((LinkRCS*) mnew)->val = val;
581 VASSERT(prev != VNULL);
582 prev->next = mnew;
583 mnew->next = curr;
584 thee->numO++;
585 done = 1;
586 flag = 1;
587 } else if (curr->idx == j) {
588 /* found the position; just add the contribution */
589 mnew = (LinkRC*)curr;
590 if (key == 0)
591 ((LinkRCS*) mnew)->val = val;
592 else
593 ((LinkRCS*) mnew)->val += val;
594 done = 1;
595 flag = 0;
596 } else if (curr->next == VNULL) {
597 /* not found; no one left, append on the end */
598 mnew = (LinkRC*)Vset_create(thee->lnkU);
599 mnew->idxT = i;
600 mnew->idx = j;
601 ((LinkRCS*) mnew)->val = val;
602 mnew->next = VNULL;
603 curr->next = (LinkRC*)mnew;
604 thee->numO++;
605 done = 1;
606 flag = 1;
607 /* not found yet; still hope */
608 } else {
609 prev = curr;
610 curr = curr->next;
611 }
612 }
613
614 if (flag == 0) {
615 /* exit now if no new links were inserted */
616 return;
617 }
618
619 /* STEP 2: Traverse colwise to find the location for insertion */
620 curr = (LinkRC*) &(((LinkRCS*) (thee->xln))[j]);
621 prev = VNULL;
622 done = 0;
623 while (!done) {
624 if (curr->idxT < i) {
625 /* not in the list; insert new link between curr and prev */
626 prev->nextT = mnew;
627 mnew->nextT = curr;
628 done = 1;
629 } else if (curr->nextT == VNULL) {
630 /* not found; no one left; set the last pointer to new link */
631 curr->nextT = mnew;
632 done = 1;
633 /* not found yet; still hope */
634 } else {
635 prev = curr;
636 curr = curr->nextT;
637 }
638 }
639 }
640
641 /*
642 * ***************************************************************************
643 * Routine: Mat_printLN
644 *
645 * Purpose: Print an LN format matrix as a DENSE matrix in MATLAB format.
646 *
647 * Author: Stephen Bond and Michael Holst
648 * ***************************************************************************
649 */
Mat_printLN(Mat * thee)650 VPUBLIC void Mat_printLN(Mat *thee)
651 {
652 int i, j;
653 int numR, numC;
654 char rn[80];
655 const int MaxRows = 30;
656 const int MaxCols = 30;
657 double matrix[30][30];
658 LinkA *mt;
659 LinkRC *mtX;
660
661 numR = thee->numR;
662 numC = thee->numC;
663
664 strncpy(rn,"Mat_printLN:",80);
665
666 /* some i/o */
667 Vnm_print(0, "%s printing <%s>" " [dim=(%dx%d),sym=%d,numA=%d]\n",
668 rn, thee->name, numR, numC, thee->sym, thee->numA);
669
670 /* size check */
671 if ((numR > MaxRows) || (numC > MaxCols)) {
672 Vnm_print(0, "%smatrix too large to view....skipping.\n", rn);
673 return;
674 }
675
676 /* make a dense matrix first */
677 for (i=0; i<numR; i++)
678 for (j=0; j<numC; j++)
679 matrix[i][j] = 0.0;
680
681 if (thee->state != NULL_STATE) {
682
683 switch (thee->format) {
684
685 case RLN_FORMAT:
686 for (i=0; i<numR; i++) {
687 for (mt=(LinkA*)Vset_access(thee->lnkU,i);
688 mt!=VNULL; mt=mt->next) {
689 if (mt->idx >= 0) {
690 j = mt->idx;
691 matrix[i][j] = mt->val;
692 }
693 }
694 }
695 break;
696
697 case CLN_FORMAT:
698 for (j=0; j<numC; j++) {
699 for (mt=(LinkA*)Vset_access(thee->lnkL,j);
700 mt!=VNULL; mt=mt->next) {
701 if (mt->idx >= 0) {
702 i = mt->idx;
703 matrix[i][j] = mt->val;
704 }
705 }
706 }
707 break;
708
709 case XLN_FORMAT:
710 for (i=0; i<numR; i++) {
711 if ( thee->sym == ISNOT_SYM ) {
712 mtX = ((LinkRC**) thee->xln)[i];
713 } else {
714 mtX = (LinkRC*) &( ((LinkRCS*) thee->xln)[i] );
715 matrix[i][i] = ((LinkRCS*) mtX)->val;
716 mtX = mtX->next;
717 }
718 for ( /* no-op */ ; mtX!=VNULL; mtX=mtX->next) {
719 j = mtX->idx;
720 if (j < numC) {
721 if ( thee->sym == ISNOT_SYM ) {
722 matrix[i][j] = ((LinkRCS*) mtX)->val;
723 } else if ( thee->sym == IS_SYM ) {
724 matrix[i][j] = ((LinkRCS*) mtX)->val;
725 matrix[j][i] = ((LinkRCS*) mtX)->val;
726 } else {
727 matrix[i][j] = ((LinkRCN*) mtX)->val;
728 matrix[j][i] = ((LinkRCN*) mtX)->valT;
729 }
730 }
731 }
732 }
733 break;
734 default:
735 Vnm_print(0,
736 "%smatrix not in correct format to print....skipping.\n", rn);
737 break;
738 }
739 }
740
741 /* print the matrix */
742 Vnm_print(0, "%s = [\n", thee->name);
743 for (i=0; i<numR; i++) {
744 for (j=0; j<numC; j++) {
745 if (VABS(matrix[i][j]) < 0.0001) {
746 Vnm_print(0, " 0.0 ");
747 } else {
748 Vnm_print(0, "%7.3f", matrix[i][j]);
749 }
750 }
751 Vnm_print(0, "\n");
752 }
753 Vnm_print(0, "];\n");
754 }
755
756 /*
757 * ***************************************************************************
758 * Routine: Mat_printLNSp
759 *
760 * Purpose: Print an LN format matrix as a SPARSE matrix in MATLAB format.
761 *
762 * Notes: flag==0 ==> write
763 * flag==1 ==> append
764 *
765 * Author: Stephen Bond and Michael Holst
766 * ***************************************************************************
767 */
Mat_printLNSp(Mat * thee,char * fname,int pflag)768 VPUBLIC void Mat_printLNSp(Mat *thee, char *fname, int pflag)
769 {
770 int i, j;
771 int numR, numC;
772 char rn[80];
773 FILE *fp;
774 LinkA *mt;
775 LinkRC *mtX;
776
777 numR = thee->numR;
778 numC = thee->numC;
779
780 strncpy(rn,"Mat_printSp:",80);
781
782 if (pflag == 0) {
783 fp=fopen(fname,"w");
784 } else if (pflag == 1) {
785 fp=fopen(fname,"a");
786 } else {
787 fp=VNULL;
788 }
789
790 if (fp==VNULL) {
791 Vnm_print(2, "%s problem opening file <%s>\n", rn, fname);
792 return;
793 }
794
795 Vnm_print(0,"%s printing <%s> [dim=(%dx%d)] to file <%s>..",
796 rn,thee->name,numR,numC,fname);
797
798 /* print the matrix in matlab sparse format */
799 fprintf(fp,"%% %s matrix <%s> [dim=(%dx%d)]\n",
800 rn, thee->name, numR, numC);
801 fprintf(fp,"%% ----------------------------------------\n");
802 fprintf(fp, "fprintf('Defining %s_IJA..');", thee->name);
803 fprintf(fp, "\n%s_IJA = [\n", thee->name);
804
805 if (thee->state != NULL_STATE) {
806
807 switch (thee->format) {
808
809 case RLN_FORMAT:
810 for (i=0; i<numR; i++) {
811 for (mt=(LinkA*)Vset_access(thee->lnkU,i);
812 mt!=VNULL; mt=mt->next) {
813 if (mt->idx >= 0) {
814 j = mt->idx;
815 fprintf(fp, " %d %d %15.8e\n",
816 i+1, j+1, mt->val);
817 }
818 }
819 }
820 break;
821
822 case CLN_FORMAT:
823 for (j=0; j<numC; j++) {
824 for (mt=(LinkA*)Vset_access(thee->lnkL,j);
825 mt!=VNULL; mt=mt->next) {
826 if (mt->idx >= 0) {
827 i = mt->idx;
828 fprintf(fp, " %d %d %15.8e\n",
829 i+1, j+1, mt->val);
830 }
831 }
832 }
833 break;
834
835 case XLN_FORMAT:
836 for (i=0; i<numR; i++) {
837 if ( thee->sym == ISNOT_SYM ) {
838 mtX = ((LinkRC**) thee->xln)[i];
839 } else {
840 mtX = (LinkRC*) &( ((LinkRCS*) thee->xln)[i] );
841 fprintf(fp, " %d %d %15.8e\n",
842 i+1, i+1, ((LinkRCS*) mtX)->val );
843 mtX = mtX->next;
844 }
845 for ( /* no-op */ ; mtX!=VNULL; mtX=mtX->next) {
846 j = mtX->idx;
847 if (j < numC) {
848 if ( thee->sym == ISNOT_SYM ) {
849 fprintf(fp, " %d %d %15.8e\n",
850 i+1, j+1, ((LinkRCS*) mtX)->val );
851 } else if ( thee->sym == IS_SYM ) {
852 fprintf(fp, " %d %d %15.8e\n",
853 i+1, j+1, ((LinkRCS*) mtX)->val );
854 fprintf(fp, " %d %d %15.8e\n",
855 j+1, i+1, ((LinkRCS*) mtX)->val );
856 } else {
857 fprintf(fp, " %d %d %15.8e\n",
858 i+1, j+1, ((LinkRCN*) mtX)->val );
859 fprintf(fp, " %d %d %15.8e\n",
860 j+1, i+1, ((LinkRCN*) mtX)->valT );
861 }
862 }
863 }
864 }
865 break;
866 default:
867 Vnm_print(0,
868 "%smatrix not in correct format to print....skipping.\n", rn);
869 break;
870 }
871 }
872
873 fprintf(fp, "];\n\n");
874 fprintf(fp, "fprintf('..done.\\n');\n");
875 fprintf(fp,"%% ----------------------------------------\n");
876 fprintf(fp,"%% Matlab code to generate sparse matrix.\n");
877 fprintf(fp,"%% ----------------------------------------\n");
878 fprintf(fp,
879 "%s = sparse(%s_IJA(:,1),%s_IJA(:,2),%s_IJA(:,3),%d,%d);\n",
880 thee->name, thee->name, thee->name, thee->name, numR, numC);
881 fprintf(fp,"%% ----------------------------------------\n");
882
883 Vnm_print(0,"..done.\n");
884
885 /* close file and return */
886 fclose(fp);
887 }
888
889
890