1 /*
2 matrix.c -- MATRIX routines.
3
4 Copyright (C) 1994-97 K. Scott Hunziker.
5 Copyright (C) 1990-94 The Boeing Company.
6
7 See the file COPYING for license, warranty, and permission details.
8 */
9
10 static char rcsid[] =
11 "$Id: matrix.c,v 1.5 1997/04/24 06:01:03 ksh Exp $";
12
13 #include "entity.h"
14 #include "scalar.h"
15 #include "vector.h"
16 #include "matrix.h"
17 #include "table.h"
18 #include "get.h"
19 #include "put.h"
20 #include "file_io.h"
21 #include "binop.h"
22 #include "sparse.h"
23 #include "dense.h"
24
25 /*
26 * This array describes the members of the
27 * MATRIX structure. The fields are `name',
28 * and `id'. The entries must be in
29 * alphabetical order, and there must be
30 * exactly one entry for each member of the
31 * MATRIX_MEMBER enumeration except END_Matrix.
32 */
33
34 MEMBER_ID matrix_member_names[] =
35 {
36 {"cid", MatrixCid},
37 {"class", MatrixClass},
38 {"density", MatrixDensity},
39 {"nc", MatrixNc},
40 {"nn", MatrixNn},
41 {"nr", MatrixNr},
42 {"order", MatrixOrder},
43 {"rid", MatrixRid},
44 {"symmetry", MatrixSymmetry},
45 {"type", MatrixType},
46 };
47
48 MATRIX_MEMBER
matrix_member_search(s)49 matrix_member_search (s)
50 char *s;
51 {
52 MEMBER_ID *m;
53
54 assert (s != NULL);
55
56 m = (MEMBER_ID *) bsearch (s, matrix_member_names, END_Matrix, sizeof (MEMBER_ID), member_cmp);
57
58 return ((m == NULL) ? END_Matrix : m->id);
59 }
60
61 ENTITY *
bi_matrix(n,p)62 bi_matrix (n, p)
63 int n;
64 ENTITY *p;
65 {
66 /*
67 * Convert `p' to matrix. NULL becomes a zero-by-zero real
68 * matrix. Scalars become one-by-one matrices, vectors become
69 * matrices with one row, and matrices are unchanged.
70 */
71
72 return p ? matrix_entity (p) : make_matrix (0, 0, real, dense);
73
74 }
75
76 ENTITY *
matrix_entity(p)77 matrix_entity (p)
78 ENTITY *p;
79 {
80 /* Convert `p' to a matrix. */
81
82 EASSERT (p, 0, 0);
83
84 switch (p->class)
85 {
86 case scalar:
87 return (scalar_to_matrix ((SCALAR *) p));
88 case vector:
89 return (vector_to_matrix ((VECTOR *) p));
90 case matrix:
91 return (p);
92 default:
93 fail ("Can't convert a %s to a matrix.", class_string[p->class]);
94 delete_entity (p);
95 break;
96 }
97 raise_exception ();
98 }
99
100 ENTITY *
scalar_to_matrix(s)101 scalar_to_matrix (s)
102 SCALAR *s;
103 {
104 /* Convert `s' to matrix. */
105
106 MATRIX *m;
107
108 EASSERT (s, scalar, 0);
109
110 m = (MATRIX *) form_matrix (1, 1, s->type, dense);
111 m->symmetry = symmetric;
112
113 switch (s->type)
114 {
115 case integer:
116 m->a.integer[0] = s->v.integer;
117 break;
118 case real:
119 m->a.real[0] = s->v.real;
120 break;
121 case complex:
122 m->a.complex[0] = s->v.complex;
123 break;
124 case character:
125 m->a.character[0] = dup_char (s->v.character);
126 break;
127 default:
128 BAD_TYPE (s->type);
129 delete_scalar (s);
130 delete_matrix (m);
131 raise_exception ();
132 }
133
134 if (s->stuff)
135 m->stuff = (TABLE *) copy_table (s->stuff);
136
137 delete_scalar (s);
138 return (ENT (m));
139 }
140
141 ENTITY *
make_matrix(nr,nc,type,density)142 make_matrix (nr, nc, type, density)
143 int nr;
144 int nc;
145 TYPE type;
146 DENSITY density;
147 {
148 /*
149 * Makes a new matrix with `nr' rows and `nc' columns. It's initialized
150 * to zero (or "" if it's character) and marked as symmetric or hermitian
151 * if it's square.
152 *
153 * NOTE: If the resulting matrix is square, then it is known to be
154 * symmetric, even if you supply your own non-zero values on the
155 * diagonal. But if the matrix is complex, then it is hermitian
156 * only by virtue of the fact that the diagonal values have zero
157 * imaginary parts. An easy way to screw up is to use this
158 * function to form a complex array, then put in complex values.
159 * That should be symmetric, not hermitian.
160 */
161
162 MATRIX *m;
163 int i;
164
165 m = (MATRIX *) form_matrix (nr, nc, type, density);
166
167 m->symmetry = (nr == nc) ?
168 ((type == complex) ? hermitian : symmetric) : general;
169
170 switch (density)
171 {
172 case dense:
173 if (m->nn)
174 memset (m->a.ptr, 0, m->nn * type_size[type]);
175 if (type == character)
176 for (i = 0; i < m->nn; i++)
177 m->a.character[i] = NULL_string;
178 break;
179 case sparse:
180 case sparse_upper:
181 break;
182 default:
183 wipeout ("Bad density.");
184 }
185
186 if (debug_level > 1)
187 inform ("Matrix created: %x.", m);
188
189 return (ENT (m));
190 }
191
192 ENTITY *
form_matrix(nr,nc,type,density)193 form_matrix (nr, nc, type, density)
194 int nr;
195 int nc;
196 TYPE type;
197 DENSITY density;
198 {
199 /*
200 * Forms a new matrix with `nr' rows and `nc' columns. If dense,
201 * memory for its data is allocated but not initialized. (No memory
202 * is allocated for data if it's sparse.) Its symmetry is general.
203 */
204
205 MATRIX *m;
206
207 assert (nr >= 0);
208 assert (nc >= 0);
209
210 m = (MATRIX *) CALLOC (1, sizeof (MATRIX));
211
212 m->entity.ref_count = 1;
213 m->entity.class = matrix;
214
215 m->symmetry = general;
216
217 m->type = type;
218 m->order = ordered;
219 m->density = density;
220
221 m->nr = nr;
222 m->nc = nc;
223
224 switch (density)
225 {
226
227 case dense:
228
229 if ((double) nr * (double) nc * type_size[type] > INT_MAX)
230 {
231 fail ("Out of memory.");
232 raise_exception ();
233 }
234 m->nn = nr * nc;
235 m->a.ptr = (m->nn) ? E_MALLOC (m->nn, type) : NULL;
236
237 /*
238 * For character type, set the first pointer to NULL so that
239 * `free_matrix' knows there's junk there and it shouldn't
240 * call FREE_CHAR on it.
241 */
242
243 if (m->nn && type == character)
244 m->a.character[0] = NULL;
245
246 break;
247
248 case sparse:
249 case sparse_upper:
250 break;
251 default:
252 wipeout ("Bad density.");
253 }
254
255 if (debug_level > 1)
256 inform ("Matrix created: %x.", m);
257
258 return (ENT (m));
259 }
260
261 ENTITY *
matrix_to_scalar(m)262 matrix_to_scalar (m)
263 MATRIX *m;
264 {
265 /* Convert `m' to scalar (only if it's one by one). */
266
267 ENTITY *s;
268
269 EASSERT (m, matrix, 0);
270
271 /* The matrix must have only 1 row and 1 column */
272
273 if (m->nr != 1 || m->nc != 1)
274 {
275 fail ("Can't convert to scalar a matrix with %d row%s and %d column%s.",
276 m->nr, PLURAL (m->nr), m->nc, PLURAL (m->nc));
277 delete_matrix (m);
278 raise_exception ();
279 }
280 else
281 {
282 switch (m->density)
283 {
284 case dense:
285 switch (m->type)
286 {
287 case integer:
288 s = int_to_scalar (m->a.integer[0]);
289 break;
290 case real:
291 s = real_to_scalar (m->a.real[0]);
292 break;
293 case complex:
294 s = complex_to_scalar (m->a.complex[0]);
295 break;
296 case character:
297 s = char_to_scalar (dup_char (m->a.character[0]));
298 break;
299 default:
300 BAD_TYPE (m->type);
301 delete_matrix (m);
302 raise_exception ();
303 }
304 break;
305 case sparse:
306 /*
307 * It's a one-by-one matrix, so it has either 0 or 1 elements. If
308 * it has 1 element, then that's the one we want.
309 */
310 if (m->nn > 0)
311 {
312 switch (m->type)
313 {
314 case integer:
315 s = int_to_scalar (m->a.integer[0]);
316 break;
317 case real:
318 s = real_to_scalar (m->a.real[0]);
319 break;
320 case complex:
321 s = complex_to_scalar (m->a.complex[0]);
322 break;
323 case character:
324 s = char_to_scalar (dup_char (m->a.character[0]));
325 break;
326 default:
327 BAD_TYPE (m->type);
328 delete_matrix (m);
329 raise_exception ();
330 }
331 }
332 else
333 {
334 s = make_scalar (m->type);
335 }
336 break;
337 case sparse_upper:
338 if (m->d.integer != NULL)
339 {
340 switch (m->type)
341 {
342 case integer:
343 s = int_to_scalar (m->d.integer[0]);
344 break;
345 case real:
346 s = real_to_scalar (m->d.real[0]);
347 break;
348 case complex:
349 s = complex_to_scalar (m->d.complex[0]);
350 break;
351 case character:
352 s = char_to_scalar (dup_char (m->d.character[0]));
353 break;
354 default:
355 BAD_TYPE (m->type);
356 delete_matrix (m);
357 raise_exception ();
358 }
359 }
360 else
361 {
362 s = make_scalar (m->type);
363 }
364 break;
365 default:
366 BAD_DENSITY (m->density);
367 delete_matrix (m);
368 raise_exception ();
369 }
370 }
371
372 if (m->stuff)
373 ((SCALAR *) s)->stuff = (TABLE *) copy_table (m->stuff);
374
375 delete_matrix (m);
376 return (s);
377 }
378
379 ENTITY *
dup_matrix(old)380 dup_matrix (old)
381 MATRIX *old;
382 {
383 /*
384 * In Algae, duplicating an entity means that we really make another
385 * copy of it in memory (not just incrementing its reference count).
386 */
387
388 MATRIX *new;
389 int i;
390
391 EASSERT (old, matrix, 0);
392
393 if (old->entity.ref_count == 1)
394 return (ENT (old));
395
396 new = (MATRIX *) form_matrix (old->nr, old->nc, old->type, old->density);
397
398 new->symmetry = old->symmetry;
399
400 if (old->rid)
401 new->rid = copy_entity (old->rid);
402 if (old->cid)
403 new->cid = copy_entity (old->cid);
404
405 new->nn = old->nn;
406
407 if (old->ia)
408 new->ia = (int *) dup_mem (old->ia, (old->nr + 1) * sizeof (int));
409
410 if (old->ja)
411 new->ja = (int *) dup_mem (old->ja, old->nn * sizeof (int));
412
413 if (old->a.integer)
414 {
415 if (new->a.integer == NULL)
416 new->a.integer = (int *) E_MALLOC (old->nn, old->type);
417 memcpy (new->a.integer, old->a.integer,
418 old->nn * type_size[old->type]);
419 if (old->type == character)
420 for (i = 0; i < old->nn; i++)
421 new->a.character[i] = dup_char (new->a.character[i]);
422 }
423
424 if (old->d.integer)
425 {
426 new->d.integer = (int *) dup_mem (old->d.integer,
427 old->nr * type_size[old->type]);
428 if (old->type == character)
429 for (i = 0; i < old->nr; i++)
430 new->d.character[i] = dup_char (new->d.character[i]);
431 }
432
433 if (debug_level > 1)
434 inform ("Matrix created: %x.", new);
435
436 delete_matrix (old);
437 return (ENT (new));
438 }
439
440 ENTITY *
matrix_to_vector(m)441 matrix_to_vector (m)
442 MATRIX *m;
443 {
444 /* Convert matrix `m' to a vector (if it's a single row or column). */
445
446 VECTOR *v;
447 int i, j;
448
449 EASSERT (m, matrix, 0);
450
451 /* Make sure it's a single row or column. */
452
453 if (m->nr != 1 && m->nc != 1)
454 {
455 fail ("Can't convert a %d by %d matrix into a vector.",
456 m->nr, m->nc);
457 delete_matrix (m);
458 raise_exception ();
459 }
460
461 if (m->density == sparse_upper)
462 {
463
464 /* If it's sparse_upper, then it has to be a 1x1 matrix. */
465
466 assert (m->nr == 1 && m->nc == 1);
467 v = (VECTOR *) make_vector (1, m->type, dense);
468 if (m->d.ptr != NULL)
469 memcpy (v->a.ptr, m->d.ptr,
470 type_size[m->type]);
471 if (v->type == character)
472 v->a.character[0] = dup_char (v->a.character[0]);
473
474 }
475 else if (m->entity.ref_count == 1)
476 {
477
478 /* If it has a ref_count of 1, we'll take a shortcut. */
479
480 v = (VECTOR *) CALLOC (1, sizeof (VECTOR));
481 if (debug_level > 1)
482 inform ("Vector created: %x.", v);
483 v->entity.ref_count = 1;
484 v->entity.class = vector;
485 v->type = m->type;
486 v->order = m->order;
487 v->density = m->density;
488 v->nn = m->nn;
489
490 /*
491 * For sparse, if `m' was a row, then `ja' is just what we want.
492 * Otherwise, we have to fill in `ja' with the element numbers.
493 */
494
495 v->ja = m->ja;
496 if (v->ja && m->nr != 1)
497 {
498 for (i = j = 0; i < m->nr; i++)
499 {
500 if (m->ia[i + 1] > m->ia[i])
501 v->ja[j++] = i + 1;
502 }
503 }
504 m->ja = NULL;
505
506 /*
507 * We don't have to worry about cleaning up `m->ia'---it's done
508 * in `delete_matrix' below.
509 */
510
511 v->a = m->a;
512 m->a.ptr = NULL;
513
514 }
515 else
516 {
517
518 /* Here's the general case. */
519
520 v = (VECTOR *) form_vector ((m->nr != 1) ? m->nr :
521 m->nc, m->type, m->density);
522 if (m->nn > 0)
523 {
524 v->nn = m->nn;
525 if (m->density == sparse)
526 {
527 if (m->nr == 1)
528 {
529 v->ja = (int *) dup_mem (m->ja, m->nn * sizeof (int));
530 }
531 else
532 {
533 v->ja = (int *) MALLOC (v->nn * sizeof (int));
534 for (i = j = 0; i < m->nr; i++)
535 {
536 if (m->ia[i + 1] > m->ia[i])
537 v->ja[j++] = i + 1;
538 }
539 }
540 v->a.ptr = dup_mem (m->a.ptr, m->nn * type_size[m->type]);
541 }
542 else
543 {
544 assert (m->density == dense);
545 memcpy (v->a.ptr, m->a.ptr, m->nn * type_size[m->type]);
546 }
547 if (v->type == character)
548 {
549 for (i = 0; i < v->nn; i++)
550 v->a.character[i] = dup_char (v->a.character[i]);
551 }
552 }
553 }
554
555 if (m->nr != 1)
556 {
557 v->ne = m->nr;
558 if (m->rid != NULL)
559 v->eid = copy_entity (m->rid);
560 }
561 else
562 {
563 v->ne = m->nc;
564 if (m->cid != NULL)
565 v->eid = copy_entity (m->cid);
566 }
567
568 if (m->stuff)
569 v->stuff = (TABLE *) copy_table (m->stuff);
570
571 delete_matrix (m);
572 return (apt_vector (v));
573 }
574
575 ENTITY *
vector_to_matrix(v)576 vector_to_matrix (v)
577 VECTOR *v;
578 {
579 /* Convert vector `v' to a matrix (a single row). */
580
581 MATRIX *m;
582 int i;
583
584 EASSERT (v, vector, 0);
585
586 if (v->entity.ref_count == 1)
587 {
588 m = (MATRIX *) CALLOC (1, sizeof (MATRIX));
589 if (debug_level > 1)
590 inform ("Matrix created: %x.", m);
591 m->entity.ref_count = 1;
592 m->entity.class = matrix;
593 m->symmetry = (v->ne == 1) ? symmetric : general;
594 m->type = v->type;
595 m->order = v->order;
596 m->density = v->density;
597 if (v->eid != NULL)
598 m->cid = copy_entity (v->eid);
599 m->nr = 1;
600 m->nc = v->ne;
601 m->nn = v->nn;
602 if (v->density == sparse && v->nn > 0)
603 {
604 m->ia = (int *) CALLOC (2, sizeof (int));
605 m->ia[0] = 1;
606 m->ia[1] = 1 + m->nn;
607 }
608 m->ja = v->ja;
609 v->ja = NULL;
610 m->a = v->a;
611 v->a.ptr = NULL;
612 }
613 else
614 {
615 m = (MATRIX *) make_matrix (1, v->ne, v->type, v->density);
616 m->symmetry = (v->ne == 1) ? symmetric : general;
617 if (v->eid != NULL)
618 m->cid = copy_entity (v->eid);
619 if (v->nn > 0)
620 {
621 m->nn = v->nn;
622 if (v->density == sparse)
623 {
624 m->ia = (int *) CALLOC (2, sizeof (int));
625 m->ia[0] = 1;
626 m->ia[1] = 1 + m->nn;
627 m->ja = (int *) dup_mem (v->ja, m->nn * sizeof (int));
628 m->a.ptr = dup_mem (v->a.ptr, m->nn * type_size[v->type]);
629 }
630 else
631 {
632 memcpy (m->a.ptr, v->a.ptr, m->nn * type_size[v->type]);
633 }
634 if (m->type == character)
635 {
636 for (i = 0; i < m->nn; i++)
637 m->a.character[i] = dup_char (m->a.character[i]);
638 }
639 }
640 }
641
642 if (v->stuff)
643 m->stuff = (TABLE *) copy_table (v->stuff);
644
645 delete_vector (v);
646 return (apt_matrix (m));
647 }
648
649 void
free_matrix(p)650 free_matrix (p)
651 MATRIX *p;
652 {
653 /* Destroy the matrix, and free up the space it occupied. */
654
655 int i;
656
657 assert (p->entity.ref_count == 0);
658
659 delete_table (p->stuff);
660 delete_entity (p->rid);
661 delete_entity (p->cid);
662
663 if (p->ia)
664 FREE (p->ia);
665 if (p->ja)
666 FREE (p->ja);
667
668 if (p->a.ptr)
669 {
670
671 /*
672 * For character type, the understanding is that the pointers are
673 * all valid (and we can call FREE_CHAR on them, provided that
674 * p->a.character[0] is not NULL.
675 */
676
677 if (p->type == character && p->a.character[0])
678 {
679 for (i = 0; i < p->nn; i++)
680 FREE_CHAR (p->a.character[i]);
681 }
682
683 FREE (p->a.ptr);
684 }
685
686 if (p->d.ptr)
687 {
688
689 /* Same as above for character type, but with p->d.character[0]. */
690
691 if (p->type == character && p->d.character[0])
692 {
693 for (i = 0; i < p->nr; i++)
694 FREE_CHAR (p->d.character[i]);
695 }
696
697 FREE (p->d.integer);
698 }
699
700 /* Just to make it harder to use it again inadvertently. */
701
702 p->entity.class = undefined_class;
703 FREE (p);
704 }
705
706 void
DB_delete_matrix(p,file,line)707 DB_delete_matrix (p, file, line)
708 MATRIX *p;
709 char *file;
710 int line;
711 {
712 /*
713 * This is the DEBUG version of `delete_matrix'. It decrements
714 * the matrix's reference count, and frees it if it is unreferenced.
715 */
716
717 if (p)
718 {
719
720 if (--p->entity.ref_count < 0)
721 {
722 wipeout ("A matrix's \"ref_count\" went below zero: %s, %d.",
723 file, line);
724 }
725
726 if (p->entity.ref_count >= 1000 || debug_level > 1)
727 {
728 inform ("matrix \"ref_count\" decrement: %x, %d, %s, %d.",
729 p, p->entity.ref_count, file, line);
730 }
731
732 if (p->entity.ref_count == 0)
733 free_matrix (p);
734
735 }
736 }
737
738 int
put_matrix(s,stream,ent_tree)739 put_matrix (s, stream, ent_tree)
740 MATRIX *s;
741 FILE *stream;
742 struct ent_node *ent_tree;
743 {
744 /* Write matrix `s' out in binary form to file `stream'. */
745
746 int size, i;
747
748 EASSERT (s, matrix, 0);
749
750 assert (s->order == ordered); /* unordered is obsolete */
751
752 if (!WRITE_INT (&s->type, stream) ||
753 !WRITE_INT (&s->density, stream) ||
754 !WRITE_INT (&s->symmetry, stream) ||
755 !WRITE_INT (&s->nr, stream) ||
756 !WRITE_INT (&s->nc, stream) ||
757 !WRITE_INT (&s->nn, stream))
758 goto err;
759
760 if (s->a.ptr)
761 {
762
763 i = 1; /* values follow */
764 if (!WRITE_INT (&i, stream))
765 goto err;
766
767 switch (s->type)
768 {
769
770 case integer:
771 if (!WRITE_INTS (s->a.integer, s->nn, stream))
772 goto err;
773 break;
774
775 case real:
776 if (!WRITE_DOUBLES (s->a.real, s->nn, stream))
777 goto err;
778 break;
779
780 case complex:
781 if (!WRITE_DOUBLES (s->a.real, 2 * s->nn, stream))
782 goto err;
783 break;
784
785 case character:
786 for (i = 0; i < s->nn; i++)
787 {
788 size = strlen (s->a.character[i]);
789 if (!WRITE_INT (&size, stream))
790 goto err;
791 if (size > 0 &&
792 fwrite (s->a.character[i], 1, size, stream) <
793 (size_t) size)
794 {
795 WRITE_WARN (stream);
796 goto err;
797 }
798 }
799 break;
800
801 default:
802 BAD_TYPE (s->type);
803 delete_matrix (s);
804 raise_exception ();
805 }
806
807 }
808 else
809 {
810
811 i = 0; /* no values */
812 if (!WRITE_INT (&i, stream))
813 goto err;
814
815 }
816
817 if (s->d.ptr)
818 {
819
820 i = 1; /* diagonal values follow */
821 if (!WRITE_INT (&i, stream))
822 goto err;
823
824 switch (s->type)
825 {
826
827 case integer:
828 if (!WRITE_INTS (s->d.integer, s->nr, stream))
829 goto err;
830 break;
831
832 case real:
833 if (!WRITE_DOUBLES (s->d.real, s->nr, stream))
834 goto err;
835 break;
836
837 case complex:
838 if (!WRITE_DOUBLES (s->d.real, 2 * s->nr, stream))
839 goto err;
840 break;
841
842 case character:
843 for (i = 0; i < s->nr; i++)
844 {
845 size = strlen (s->d.character[i]);
846 if (!WRITE_INT (&size, stream))
847 goto err;
848 if (size > 0 &&
849 fwrite (s->d.character[i], 1, size, stream) <
850 (size_t) size)
851 {
852 WRITE_WARN (stream);
853 goto err;
854 }
855 }
856 break;
857 }
858
859 }
860 else
861 {
862
863 i = 0; /* no diagonal values */
864 if (!WRITE_INT (&i, stream))
865 goto err;
866
867 }
868
869 if (s->ia)
870 {
871 i = 1; /* ia follows */
872 if (!WRITE_INT (&i, stream) ||
873 !WRITE_INTS (s->ia, s->nr + 1, stream))
874 goto err;
875 }
876 else
877 {
878 i = 0; /* no ia */
879 if (!WRITE_INT (&i, stream))
880 goto err;
881 }
882
883 if (s->ja)
884 {
885 i = 1; /* ja follows */
886 if (!WRITE_INT (&i, stream) ||
887 !WRITE_INTS (s->ja, s->nn, stream))
888 goto err;
889 }
890 else
891 {
892 i = 0; /* no ja */
893 if (!WRITE_INT (&i, stream))
894 goto err;
895 }
896
897 if (s->rid)
898 {
899 i = 1; /* rid follows */
900 if (!WRITE_INT (&i, stream) ||
901 !put_entity (copy_entity (s->rid), stream, ent_tree))
902 goto err;
903 }
904 else
905 {
906 i = 0; /* no rid */
907 if (!WRITE_INT (&i, stream))
908 goto err;
909 }
910
911 if (s->cid)
912 {
913 i = 1; /* cid follows */
914 if (!WRITE_INT (&i, stream) ||
915 !put_entity (copy_entity (s->cid), stream, ent_tree))
916 goto err;
917 }
918 else
919 {
920 i = 0; /* no cid */
921 if (!WRITE_INT (&i, stream))
922 goto err;
923 }
924
925 if (s->stuff)
926 {
927 i = 1; /* stuff follows */
928 if (!WRITE_INT (&i, stream) ||
929 !put_entity (copy_table (s->stuff), stream, ent_tree))
930 goto err;
931 }
932 else
933 {
934 i = 0; /* no stuff */
935 if (!WRITE_INT (&i, stream))
936 goto err;
937 }
938
939 delete_matrix (s);
940 return 1;
941
942 err:
943 delete_matrix (s);
944 return 0;
945 }
946
947 ENTITY *
get_matrix(stream,ver)948 get_matrix (stream, ver)
949 FILE *stream;
950 int ver;
951 {
952 /* Read a vector from binary file `stream'. */
953
954 MATRIX *s;
955 int size, i;
956 static char *warn_msg = "Invalid matrix in file.";
957
958 s = (MATRIX *) CALLOC (1, sizeof (MATRIX));
959 s->entity.ref_count = 1;
960 s->entity.class = matrix;
961 s->order = ordered;
962
963 if (!READ_INT (&s->type, stream) ||
964 !READ_INT (&s->density, stream) ||
965 !READ_INT (&s->symmetry, stream) ||
966 !READ_INT (&s->nr, stream) ||
967 !READ_INT (&s->nc, stream) ||
968 !READ_INT (&s->nn, stream))
969 {
970 FREE (s);
971 return NULL;
972 }
973
974 if (!READ_INT (&i, stream))
975 { /* values follow? */
976 FREE (s);
977 return NULL;
978 }
979
980 if (i)
981 {
982
983 if (s->nn < 1)
984 {
985 warn (warn_msg);
986 FREE (s);
987 return NULL;
988 }
989
990 switch (s->type)
991 {
992
993 case integer:
994 s->a.integer = MALLOC (s->nn * sizeof (int));
995 if (!READ_INTS (s->a.integer, s->nn, stream))
996 goto err;
997 break;
998
999 case real:
1000 s->a.real = MALLOC (s->nn * sizeof (REAL));
1001 if (!READ_DOUBLES (s->a.real, s->nn, stream))
1002 goto err;
1003 break;
1004
1005 case complex:
1006 s->a.real = MALLOC (2 * s->nn * sizeof (REAL));
1007 if (!READ_DOUBLES (s->a.real, 2 * s->nn, stream))
1008 goto err;
1009 break;
1010
1011 case character:
1012 s->a.character = MALLOC (s->nn * sizeof (char *));
1013 for (i = 0; i < s->nn; i++)
1014 s->a.character[i] = NULL_string;
1015 for (i = 0; i < s->nn; i++)
1016 {
1017 if (!READ_INT (&size, stream))
1018 goto err;
1019 if (size > 0)
1020 {
1021 s->a.character[i] = (char *) MALLOC (size + 1);
1022 if (fread (s->a.character[i], 1, size, stream) <
1023 (size_t) size)
1024 {
1025 READ_WARN (stream);
1026 goto err;
1027 }
1028 s->a.character[i][size] = '\0';
1029 }
1030 }
1031 break;
1032
1033 default:
1034 warn (warn_msg);
1035 goto err;
1036
1037 }
1038
1039 }
1040 else
1041 {
1042 s->a.ptr = NULL;
1043 }
1044
1045 if (!READ_INT (&i, stream))
1046 goto err; /* diagonal values follow? */
1047
1048 if (i)
1049 {
1050
1051 if (s->nr < 1)
1052 {
1053 warn (warn_msg);
1054 goto err;
1055 }
1056
1057 switch (s->type)
1058 {
1059
1060 case integer:
1061 s->d.integer = MALLOC (s->nr * sizeof (int));
1062 if (!READ_INTS (s->d.integer, s->nr, stream))
1063 goto err;
1064 break;
1065
1066 case real:
1067 s->d.real = MALLOC (s->nr * sizeof (REAL));
1068 if (!READ_DOUBLES (s->d.real, s->nr, stream))
1069 goto err;
1070 break;
1071
1072 case complex:
1073 s->d.real = MALLOC (2 * s->nr * sizeof (REAL));
1074 if (!READ_DOUBLES (s->d.real, 2 * s->nr, stream))
1075 goto err;
1076 break;
1077
1078 case character:
1079 s->d.character = MALLOC (s->nr * sizeof (char *));
1080 for (i = 0; i < s->nr; i++)
1081 s->d.character[i] = NULL_string;
1082 for (i = 0; i < s->nr; i++)
1083 {
1084 if (!READ_INT (&size, stream))
1085 goto err;
1086 if (size > 0)
1087 {
1088 s->d.character[i] = (char *) MALLOC (size + 1);
1089 if (fread (s->d.character[i], 1, size, stream) <
1090 (size_t) size)
1091 {
1092 READ_WARN (stream);
1093 goto err;
1094 }
1095 s->d.character[i][size] = '\0';
1096 }
1097 }
1098 break;
1099
1100 }
1101
1102 }
1103 else
1104 {
1105 s->d.ptr = NULL;
1106 }
1107
1108 if (!READ_INT (&i, stream))
1109 goto err; /* ia follows? */
1110 if (i)
1111 {
1112 s->ia = (int *) MALLOC ((s->nr + 1) * sizeof (int));
1113 if (!READ_INTS (s->ia, s->nr + 1, stream))
1114 goto err;
1115 }
1116
1117 if (!READ_INT (&i, stream))
1118 goto err; /* ja follows? */
1119 if (i)
1120 {
1121 s->ja = (int *) MALLOC (s->nn * sizeof (int));
1122 if (!READ_INTS (s->ja, s->nn, stream))
1123 goto err;
1124 }
1125
1126 if (!READ_INT (&i, stream))
1127 goto err; /* rid follows */
1128 if (i && !(s->rid = get_entity (stream)))
1129 goto err;
1130
1131 if (!READ_INT (&i, stream))
1132 goto err; /* cid follows */
1133 if (i && !(s->cid = get_entity (stream)))
1134 goto err;
1135
1136 if (!READ_INT (&i, stream))
1137 goto err; /* stuff follows */
1138 if (i && !(s->stuff = (TABLE *)
1139 (ver ? get_entity (stream) : get_table (stream, ver))))
1140 goto err;
1141
1142 if (!ok_entity (ENT (s)))
1143 goto err;
1144
1145 return apt_matrix (s);
1146
1147 err:
1148 delete_matrix (s);
1149 return NULL;
1150 }
1151
1152 ENTITY *
gift_wrap_matrix(nr,nc,type,data)1153 gift_wrap_matrix (nr, nc, type, data)
1154 int nr;
1155 int nc;
1156 TYPE type;
1157 void *data;
1158 {
1159 /*
1160 * This function creates a dense MATRIX structure to point to
1161 * the given array.
1162 */
1163
1164 MATRIX *m;
1165
1166 assert (nr >= 0);
1167 assert (nc >= 0);
1168
1169 m = (MATRIX *) CALLOC (1, sizeof (MATRIX));
1170
1171 m->entity.ref_count = 1;
1172 m->entity.class = matrix;
1173
1174 m->symmetry = general;
1175
1176 m->type = type;
1177 m->order = ordered;
1178 m->density = dense;
1179 m->nr = nr;
1180 m->nc = nc;
1181
1182 if ((double) nr * nc > INT_MAX)
1183 {
1184 fail ("Integer overflow. Matrix dimensions too big.");
1185 raise_exception ();
1186 }
1187 m->nn = nr * nc;
1188 m->a.ptr = (m->nn) ? data : NULL;
1189
1190 if (debug_level > 1)
1191 inform ("Matrix created: %x.", m);
1192
1193 return (ENT (m));
1194 }
1195