1 /*
2 * MATRIX BUILD MODULE
3 *
4 * Author: Advising professor:
5 * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
6 * UC Berkeley
7 *
8 * This file contains the routines associated with clearing, loading and
9 * preprocessing the matrix for the sparse matrix routines.
10 *
11 * >>> User accessible functions contained in this file:
12 * spClear
13 * spGetElement
14 * spGetAdmittance
15 * spGetQuad
16 * spGetOnes
17 * spInstallInitInfo
18 * spGetInitInfo
19 * spInitialize
20 *
21 * >>> Other functions contained in this file:
22 * spcFindElementInCol
23 * Translate
24 * spcCreateElement
25 * spcLinkRows
26 * EnlargeMatrix
27 * ExpandTranslationArrays
28 */
29
30
31 /*
32 * Revision and copyright information.
33 *
34 * Copyright (c) 1985,86,87,88,89,90
35 * by Kenneth S. Kundert and the University of California.
36 *
37 * Permission to use, copy, modify, and distribute this software and
38 * its documentation for any purpose and without fee is hereby granted,
39 * provided that the copyright notices appear in all copies and
40 * supporting documentation and that the authors and the University of
41 * California are properly credited. The authors and the University of
42 * California make no representations as to the suitability of this
43 * software for any purpose. It is provided `as is', without express
44 * or implied warranty.
45 */
46
47 #ifdef notdef
48 static char copyright[] =
49 "Sparse1.3: Copyright (c) 1985,86,87,88,89,90 by Kenneth S. Kundert";
50 static char RCSid[] =
51 "@(#)$Header: spBuild.c,v 1.3 88/06/24 05:00:31 kundert Exp $";
52 #endif
53
54
55
56
57 /*
58 * IMPORTS
59 *
60 * >>> Import descriptions:
61 * spConfig.h
62 * Macros that customize the sparse matrix routines.
63 * spMatrix.h
64 * Macros and declarations to be imported by the user.
65 * spDefs.h
66 * Matrix type and macro definitions for the sparse matrix routines.
67 */
68
69 #define spINSIDE_SPARSE
70 #include "spconfig.h"
71 #include "spmatrix.h"
72 #include "spdefs.h"
73
74
75
76
77
78 /*
79 * Function declarations
80 */
81
82 #ifdef __STDC__
83 static void Translate( MatrixPtr, int*, int* );
84 static void EnlargeMatrix( MatrixPtr, int );
85 static void ExpandTranslationArrays( MatrixPtr, int );
86 #else /* __STDC__ */
87 static void Translate();
88 static void EnlargeMatrix();
89 static void ExpandTranslationArrays();
90 #endif /* __STDC__ */
91
92
93
94
95
96
97
98 /* Routines for JSPICE3 support
99 * S. Whiteley 1993
100 */
101
102 /*
103 * COPY INIT TO REAL PART OF MATRIX
104 *
105 * Sets the real component to the value of the init component.
106 * This is useful for initializing the real matrix with values
107 * stored in the init part, and more efficient than the
108 * INITIALIZE option routines.
109 *
110 * >>> Arguments:
111 * Matrix <input> (char *)
112 * Pointer to matrix that is to be cleared.
113 *
114 * >>> Local variables:
115 * pElement (ElementPtr)
116 * A pointer to the element being cleared.
117 */
118
119 void
spItoR(eMatrix)120 spItoR( eMatrix )
121
122 char *eMatrix;
123 {
124 MatrixPtr Matrix = (MatrixPtr)eMatrix;
125 register ElementPtr pElement;
126 register int I;
127
128 /* Begin `spClear'. */
129 ASSERT( IS_SPARSE( Matrix ) );
130
131 #if spCOMPLEX
132 { for (I = Matrix->Size; I > 0; I--)
133 { pElement = Matrix->FirstInCol[I];
134 while (pElement != NULL)
135 { pElement->Real = pElement->Init;
136 pElement = pElement->NextInCol;
137 }
138 }
139 }
140 #endif
141
142 /* Empty the trash. */
143 Matrix->TrashCan.Real = 0.0;
144 #if spCOMPLEX
145 Matrix->TrashCan.Imag = 0.0;
146 #endif
147
148 Matrix->Error = spOKAY;
149 Matrix->Factored = NO;
150 Matrix->SingularCol = 0;
151 Matrix->SingularRow = 0;
152 Matrix->PreviousMatrixWasComplex = Matrix->Complex;
153 return;
154 }
155
156
157
158 /*
159 * COPY REAL PART TO INIT PART OF MATRIX
160 *
161 * Sets the init component to the value of the real component.
162 * This is useful for initializing the real matrix with values
163 * stored in the init part, and more efficient than the
164 * INITIALIZE option routines.
165 *
166 * >>> Arguments:
167 * Matrix <input> (char *)
168 * Pointer to matrix that is to be cleared.
169 *
170 * >>> Local variables:
171 * pElement (ElementPtr)
172 * A pointer to the element being cleared.
173 */
174
175 void
spRtoI(eMatrix)176 spRtoI( eMatrix )
177
178 char *eMatrix;
179 {
180 MatrixPtr Matrix = (MatrixPtr)eMatrix;
181 register ElementPtr pElement;
182 register int I;
183
184 /* Begin `spClear'. */
185 ASSERT( IS_SPARSE( Matrix ) );
186
187 #if spCOMPLEX
188 { for (I = Matrix->Size; I > 0; I--)
189 { pElement = Matrix->FirstInCol[I];
190 while (pElement != NULL)
191 { pElement->Init = pElement->Real;
192 pElement = pElement->NextInCol;
193 }
194 }
195 }
196 #endif
197
198 /* Empty the trash. */
199 Matrix->TrashCan.Real = 0.0;
200 #if spCOMPLEX
201 Matrix->TrashCan.Imag = 0.0;
202 #endif
203
204 Matrix->Error = spOKAY;
205 Matrix->Factored = NO;
206 Matrix->SingularCol = 0;
207 Matrix->SingularRow = 0;
208 Matrix->PreviousMatrixWasComplex = Matrix->Complex;
209 return;
210 }
211
212 /*
213 * LOAD GMIN
214 *
215 * This routine adds Gmin to each diagonal element. Because Gmin is
216 * added to the current diagonal, which may bear little relation to
217 * what the outside world thinks is a diagonal, and because the
218 * elements that are diagonals may change after calling spOrderAndFactor,
219 * use of this routine is not recommended. It is included here simply
220 * for compatibility with Spice3.
221 */
222
223 void
spLoadGmin(Matrixp,Gmin)224 spLoadGmin( Matrixp, Gmin )
225
226 char *Matrixp;
227 double Gmin;
228 {
229 MatrixPtr Matrix = (MatrixPtr)Matrixp;
230 int I;
231 ArrayOfElementPtrs Diag;
232 ElementPtr diag;
233
234 /* Begin `LoadGmin'. */
235 ASSERT( IS_SPARSE( Matrix ) );
236
237 if (Gmin != 0.0) {
238 Diag = Matrix->Diag;
239 for (I = Matrix->Size; I > 0; I--) {
240 if (diag = Diag[I])
241 diag->Real += Gmin;
242 }
243 }
244 return;
245 }
246
247
248 /*
249 * Return some interesting things
250 */
251 void
spGetStat(Matrixp,Size,NonZero,FillIns)252 spGetStat( Matrixp, Size, NonZero, FillIns )
253
254 char *Matrixp;
255 int *Size, *NonZero, *FillIns;
256 {
257 MatrixPtr Matrix = (MatrixPtr)Matrixp;
258
259 if (Matrix == (MatrixPtr)0) return;
260 *Size = Matrix->Size;
261 *NonZero = Matrix->Elements;
262 *FillIns = Matrix->Fillins;
263 }
264
265
266 /*
267 * spAddCol()
268 */
269 int
spAddCol(Matrixp,Accum_Col,Addend_Col)270 spAddCol( Matrixp, Accum_Col, Addend_Col )
271
272 char *Matrixp;
273 int Accum_Col, Addend_Col;
274 {
275 MatrixPtr Matrix = (MatrixPtr)Matrixp;
276 ElementPtr Accum, Addend, Prev;
277
278 Accum_Col = Matrix->ExtToIntColMap[Accum_Col];
279 Addend_Col = Matrix->ExtToIntColMap[Accum_Col];
280
281 Accum = Matrix->FirstInCol[Accum_Col];
282 Addend = Matrix->FirstInCol[Addend_Col];
283 Prev = NULL;
284 while (Addend != NULL) {
285 while (Accum && Accum->Row < Addend->Row) {
286 Prev = Accum;
287 Accum = Accum->NextInCol;
288 }
289 if (!Accum || Accum->Row > Addend->Row) {
290 Accum =
291 spcCreateElement(Matrix, Addend->Row, Accum_Col, &Prev, 0);
292 }
293 Accum->Real += Addend->Real;
294 Accum->Imag += Addend->Imag;
295 Addend = Addend->NextInCol;
296 }
297
298 return spError( (char *)Matrix );
299 }
300
301
302 /*
303 * spZeroCol()
304 */
305 int
spZeroCol(Matrixp,Col)306 spZeroCol( Matrixp, Col )
307
308 char *Matrixp;
309 int Col;
310 {
311 MatrixPtr Matrix = (MatrixPtr)Matrixp;
312 ElementPtr Element;
313
314 Col = Matrix->ExtToIntColMap[Col];
315
316 for (Element = Matrix->FirstInCol[Col]; Element != NULL;
317 Element = Element->NextInCol) {
318 Element->Real = 0.0;
319 Element->Imag = 0.0;
320 }
321
322 return spError( (char *)Matrix );
323 }
324
325
326 #ifndef M_LN2
327 #define M_LN2 0.69314718055994530942
328 #endif
329 #ifndef M_LN10
330 #define M_LN10 2.30258509299404568402
331 #endif
332
333 /* Avoid messing with checking for the presence of these functions
334 * in the library, and whether logb() returns int or double.
335 */
336 static double
my_logb(x)337 my_logb(x)
338
339 double x;
340 {
341 double y = 0.0;
342
343 if (x != 0.0) {
344 if (x < 0.0)
345 x = - x;
346 while (x > 2.0) {
347 y += 1.0;
348 x /= 2.0;
349 }
350 while (x < 1.0) {
351 y -= 1.0;
352 x *= 2.0;
353 }
354 }
355 else
356 y = 0.0;
357
358 return (y);
359 }
360
361
362 static double
my_scalb(x,n)363 my_scalb(x, n)
364
365 double x;
366 int n;
367 {
368 double y, z = 1.0, k = 2.0;
369
370 if (n < 0) {
371 n = -n;
372 k = 0.5;
373 }
374
375 if (x != 0.0)
376 for (y = 1.0; n; n >>= 1) {
377 y *= k;
378 if (n & 1)
379 z *= y;
380 }
381
382 return (x * z);
383 }
384
385
386 /*
387 * spDProd()
388 */
389 int
spDProd(Matrix,pMantissaR,pMantissaI,pExponent)390 spDProd( Matrix, pMantissaR, pMantissaI, pExponent)
391
392 char *Matrix;
393 double *pMantissaR;
394 double *pMantissaI;
395 int *pExponent;
396 {
397 double re, im, x, y, z;
398 int p;
399
400 spDeterminant( (char *)Matrix, &p, &re, &im);
401
402 #ifdef debug_print
403 printf("Determinant 10: (%20g,%20g)^%d\n", re, im, p);
404 #endif
405
406 /* Convert base 10 numbers to base 2 numbers, for comparison */
407 y = p * M_LN10 / M_LN2;
408 x = (int) y;
409 y -= x;
410
411 /* ASSERT
412 * x = integral part of exponent, y = fraction part of exponent
413 */
414
415 /* Fold in the fractional part */
416 #ifdef debug_print
417 printf(" ** base10 -> base2 int = %g, frac = %20g\n", x, y);
418 #endif
419 z = pow(2.0, y);
420 re *= z;
421 im *= z;
422 #ifdef debug_print
423 printf(" ** multiplier = %20g\n", z);
424 #endif
425
426 /* Re-normalize (re or im may be > 2.0 or both < 1.0 */
427 if (re != 0.0) {
428 y = my_logb(re);
429 if (im != 0.0)
430 z = my_logb(im);
431 else
432 z = 0;
433 }
434 else if (im != 0.0) {
435 z = my_logb(im);
436 y = 0;
437 }
438 else {
439 /* Singular */
440 /*printf("10 -> singular\n");*/
441 y = 0;
442 z = 0;
443 }
444
445 #ifdef debug_print
446 printf(" ** renormalize changes = %g,%g\n", y, z);
447 #endif
448 if (y < z)
449 y = z;
450
451 *pExponent = x + y;
452 x = my_scalb(re, (int) -y);
453 z = my_scalb(im, (int) -y);
454 #ifdef debug_print
455 printf(" ** values are: re %g, im %g, y %g, re' %g, im' %g\n",
456 re, im, y, x, z);
457 #endif
458 *pMantissaR = my_scalb(re, (int) -y);
459 *pMantissaI = my_scalb(im, (int) -y);
460
461 #ifdef debug_print
462 printf("Determinant 10->2: (%20g,%20g)^%d\n", *pMantissaR,
463 *pMantissaI, *pExponent);
464 #endif
465 return spError( (char *)Matrix );
466 }
467
468
469 /* end of added routines */
470
471
472
473
474
475
476
477 /*
478 * CLEAR MATRIX
479 *
480 * Sets every element of the matrix to zero and clears the error flag.
481 *
482 * >>> Arguments:
483 * Matrix <input> (char *)
484 * Pointer to matrix that is to be cleared.
485 *
486 * >>> Local variables:
487 * pElement (ElementPtr)
488 * A pointer to the element being cleared.
489 */
490
491 void
spClear(eMatrix)492 spClear( eMatrix )
493
494 char *eMatrix;
495 {
496 MatrixPtr Matrix = (MatrixPtr)eMatrix;
497 register ElementPtr pElement;
498 register int I;
499
500 /* Begin `spClear'. */
501 ASSERT( IS_SPARSE( Matrix ) );
502
503 /* Clear matrix. */
504 #if spCOMPLEX
505 if (Matrix->PreviousMatrixWasComplex OR Matrix->Complex)
506 { for (I = Matrix->Size; I > 0; I--)
507 { pElement = Matrix->FirstInCol[I];
508 while (pElement != NULL)
509 { pElement->Real = 0.0;
510 pElement->Imag = 0.0;
511 pElement = pElement->NextInCol;
512 }
513 }
514 }
515 else
516 #endif
517 { for (I = Matrix->Size; I > 0; I--)
518 { pElement = Matrix->FirstInCol[I];
519 while (pElement != NULL)
520 { pElement->Real = 0.0;
521 pElement = pElement->NextInCol;
522 }
523 }
524 }
525
526 /* Empty the trash. */
527 Matrix->TrashCan.Real = 0.0;
528 #if spCOMPLEX
529 Matrix->TrashCan.Imag = 0.0;
530 #endif
531
532 Matrix->Error = spOKAY;
533 Matrix->Factored = NO;
534 Matrix->SingularCol = 0;
535 Matrix->SingularRow = 0;
536 Matrix->PreviousMatrixWasComplex = Matrix->Complex;
537 return;
538 }
539
540
541
542
543
544
545
546
547
548
549
550 /*
551 * SINGLE ELEMENT ADDITION TO MATRIX BY INDEX
552 *
553 * Finds element [Row,Col] and returns a pointer to it. If element is
554 * not found then it is created and spliced into matrix. This routine
555 * is only to be used after spCreate() and before spMNA_Preorder(),
556 * spFactor() or spOrderAndFactor(). Returns a pointer to the
557 * Real portion of a MatrixElement. This pointer is later used by
558 * spADD_xxx_ELEMENT to directly access element.
559 *
560 * >>> Returns:
561 * Returns a pointer to the element. This pointer is then used to directly
562 * access the element during successive builds.
563 *
564 * >>> Arguments:
565 * Matrix <input> (char *)
566 * Pointer to the matrix that the element is to be added to.
567 * Row <input> (int)
568 * Row index for element. Must be in the range of [0..Size] unless
569 * the options EXPANDABLE or TRANSLATE are used. Elements placed in
570 * row zero are discarded. In no case may Row be less than zero.
571 * Col <input> (int)
572 * Column index for element. Must be in the range of [0..Size] unless
573 * the options EXPANDABLE or TRANSLATE are used. Elements placed in
574 * column zero are discarded. In no case may Col be less than zero.
575 *
576 * >>> Local variables:
577 * pElement (RealNumber *)
578 * Pointer to the element.
579 *
580 * >>> Possible errors:
581 * spNO_MEMORY
582 * Error is not cleared in this routine.
583 */
584
585 RealNumber *
spGetElement(eMatrix,Row,Col)586 spGetElement( eMatrix, Row, Col )
587
588 char *eMatrix;
589 int Row, Col;
590 {
591 MatrixPtr Matrix = (MatrixPtr)eMatrix;
592 RealNumber *pElement;
593 ElementPtr spcFindElementInCol();
594 void Translate();
595
596 /* Begin `spGetElement'. */
597 ASSERT( IS_SPARSE( Matrix ) AND Row >= 0 AND Col >= 0 );
598
599 if ((Row == 0) OR (Col == 0))
600 return &Matrix->TrashCan.Real;
601
602 #if NOT TRANSLATE
603 ASSERT(Matrix->NeedsOrdering);
604 #endif
605
606 #if TRANSLATE
607 Translate( Matrix, &Row, &Col );
608 if (Matrix->Error == spNO_MEMORY) return NULL;
609 #endif
610
611 #if NOT TRANSLATE
612 #if NOT EXPANDABLE
613 ASSERT(Row <= Matrix->Size AND Col <= Matrix->Size);
614 #endif
615
616 #if EXPANDABLE
617 /* Re-size Matrix if necessary. */
618 if ((Row > Matrix->Size) OR (Col > Matrix->Size))
619 EnlargeMatrix( Matrix, MAX(Row, Col) );
620 if (Matrix->Error == spNO_MEMORY) return NULL;
621 #endif
622 #endif
623
624 /*
625 * The condition part of the following if statement tests to see if the
626 * element resides along the diagonal, if it does then it tests to see
627 * if the element has been created yet (Diag pointer not NULL). The
628 * pointer to the element is then assigned to Element after it is cast
629 * into a pointer to a RealNumber. This casting makes the pointer into
630 * a pointer to Real. This statement depends on the fact that Real
631 * is the first record in the MatrixElement structure.
632 */
633
634 if ((Row != Col) OR ((pElement = (RealNumber *)Matrix->Diag[Row]) == NULL))
635 {
636 /*
637 * Element does not exist or does not reside along diagonal. Search
638 * column for element. As in the if statement above, the pointer to the
639 * element which is returned by spcFindElementInCol is cast into a
640 * pointer to Real, a RealNumber.
641 */
642 pElement = (RealNumber*)spcFindElementInCol( Matrix,
643 &(Matrix->FirstInCol[Col]),
644 Row, Col, YES );
645 }
646 return pElement;
647 }
648
649
650
651
652
653
654
655
656
657
658
659 /*
660 * FIND ELEMENT BY SEARCHING COLUMN
661 *
662 * Searches column starting at element specified at PtrAddr and finds element
663 * in Row. If Element does not exists, it is created. The pointer to the
664 * element is returned.
665 *
666 * >>> Returned:
667 * A pointer to the desired element:
668 *
669 * >>> Arguments:
670 * Matrix <input> (MatrixPtr)
671 * Pointer to Matrix.
672 * LastAddr <input-output> (ElementPtr *)
673 * Address of pointer that initially points to the element in Col at which
674 * the search is started. The pointer in this location may be changed if
675 * a fill-in is required in and adjacent element. For this reason it is
676 * important that LastAddr be the address of a FirstInCol or a NextInCol
677 * rather than a temporary variable.
678 * Row <input> (int)
679 * Row being searched for.
680 * Col (int)
681 * Column being searched.
682 * CreateIfMissing <input> (BOOLEAN)
683 * Indicates what to do if element is not found, create one or return a
684 * NULL pointer.
685 *
686 * Local variables:
687 * pElement (ElementPtr)
688 * Pointer used to search through matrix.
689 */
690
691 ElementPtr
spcFindElementInCol(Matrix,LastAddr,Row,Col,CreateIfMissing)692 spcFindElementInCol( Matrix, LastAddr, Row, Col, CreateIfMissing )
693
694 MatrixPtr Matrix;
695 register ElementPtr *LastAddr;
696 register int Row;
697 int Col;
698 BOOLEAN CreateIfMissing;
699 {
700 register ElementPtr pElement;
701 ElementPtr spcCreateElement();
702
703 /* Begin `spcFindElementInCol'. */
704 pElement = *LastAddr;
705
706 /* Search for element. */
707 while (pElement != NULL)
708 { if (pElement->Row < Row)
709 {
710 /* Have not reached element yet. */
711 LastAddr = &(pElement->NextInCol);
712 pElement = pElement->NextInCol;
713 }
714 else if (pElement->Row == Row)
715 {
716 /* Reached element. */
717 return pElement;
718 }
719 else break; /* while loop */
720 }
721
722 /* Element does not exist and must be created. */
723 if (CreateIfMissing)
724 return spcCreateElement( Matrix, Row, Col, LastAddr, NO );
725 else return NULL;
726 }
727
728
729
730
731
732
733
734
735 #if TRANSLATE
736
737 /*
738 * TRANSLATE EXTERNAL INDICES TO INTERNAL
739 *
740 * Convert external row and column numbers to internal row and column numbers.
741 * Also updates Ext/Int maps.
742 *
743 *
744 * >>> Arguments:
745 * Matrix <input> (MatrixPtr)
746 * Pointer to the matrix.
747 * Row <input/output> (int *)
748 * Upon entry Row is either a external row number or an external node
749 * number. Upon return, the internal equivalent is supplied.
750 * Col <input/output> (int *)
751 * Upon entry Column is either a external column number or an external node
752 * number. Upon return, the internal equivalent is supplied.
753 *
754 * >>> Local variables:
755 * ExtCol (int)
756 * Temporary variable used to hold the external column or node number
757 * during the external to internal column number translation.
758 * ExtRow (int)
759 * Temporary variable used to hold the external row or node number during
760 * the external to internal row number translation.
761 * IntCol (int)
762 * Temporary variable used to hold the internal column or node number
763 * during the external to internal column number translation.
764 * IntRow (int)
765 * Temporary variable used to hold the internal row or node number during
766 * the external to internal row number translation.
767 */
768
769 static void
Translate(Matrix,Row,Col)770 Translate( Matrix, Row, Col )
771
772 MatrixPtr Matrix;
773 int *Row, *Col;
774 {
775 register int IntRow, IntCol, ExtRow, ExtCol;
776
777 /* Begin `Translate'. */
778 ExtRow = *Row;
779 ExtCol = *Col;
780
781 /* Expand translation arrays if necessary. */
782 if ((ExtRow > Matrix->AllocatedExtSize) OR
783 (ExtCol > Matrix->AllocatedExtSize))
784 {
785 ExpandTranslationArrays( Matrix, MAX(ExtRow, ExtCol) );
786 if (Matrix->Error == spNO_MEMORY) return;
787 }
788
789 /* Set ExtSize if necessary. */
790 if ((ExtRow > Matrix->ExtSize) OR (ExtCol > Matrix->ExtSize))
791 Matrix->ExtSize = MAX(ExtRow, ExtCol);
792
793 /* Translate external row or node number to internal row or node number. */
794 if ((IntRow = Matrix->ExtToIntRowMap[ExtRow]) == -1)
795 { Matrix->ExtToIntRowMap[ExtRow] = ++Matrix->CurrentSize;
796 Matrix->ExtToIntColMap[ExtRow] = Matrix->CurrentSize;
797 IntRow = Matrix->CurrentSize;
798
799 #if NOT EXPANDABLE
800 ASSERT(IntRow <= Matrix->Size);
801 #endif
802
803 #if EXPANDABLE
804 /* Re-size Matrix if necessary. */
805 if (IntRow > Matrix->Size)
806 EnlargeMatrix( Matrix, IntRow );
807 if (Matrix->Error == spNO_MEMORY) return;
808 #endif
809
810 Matrix->IntToExtRowMap[IntRow] = ExtRow;
811 Matrix->IntToExtColMap[IntRow] = ExtRow;
812 }
813
814 /* Translate external column or node number to internal column or node number.*/
815 if ((IntCol = Matrix->ExtToIntColMap[ExtCol]) == -1)
816 { Matrix->ExtToIntRowMap[ExtCol] = ++Matrix->CurrentSize;
817 Matrix->ExtToIntColMap[ExtCol] = Matrix->CurrentSize;
818 IntCol = Matrix->CurrentSize;
819
820 #if NOT EXPANDABLE
821 ASSERT(IntCol <= Matrix->Size);
822 #endif
823
824 #if EXPANDABLE
825 /* Re-size Matrix if necessary. */
826 if (IntCol > Matrix->Size)
827 EnlargeMatrix( Matrix, IntCol );
828 if (Matrix->Error == spNO_MEMORY) return;
829 #endif
830
831 Matrix->IntToExtRowMap[IntCol] = ExtCol;
832 Matrix->IntToExtColMap[IntCol] = ExtCol;
833 }
834
835 *Row = IntRow;
836 *Col = IntCol;
837 return;
838 }
839 #endif
840
841
842
843
844
845
846 #if QUAD_ELEMENT
847 /*
848 * ADDITION OF ADMITTANCE TO MATRIX BY INDEX
849 *
850 * Performs same function as spGetElement except rather than one
851 * element, all four Matrix elements for a floating component are
852 * added. This routine also works if component is grounded. Positive
853 * elements are placed at [Node1,Node2] and [Node2,Node1]. This
854 * routine is only to be used after spCreate() and before
855 * spMNA_Preorder(), spFactor() or spOrderAndFactor().
856 *
857 * >>> Returns:
858 * Error code.
859 *
860 * >>> Arguments:
861 * Matrix <input> (char *)
862 * Pointer to the matrix that component is to be entered in.
863 * Node1 <input> (int)
864 * Row and column indices for elements. Must be in the range of [0..Size]
865 * unless the options EXPANDABLE or TRANSLATE are used. Node zero is the
866 * ground node. In no case may Node1 be less than zero.
867 * Node2 <input> (int)
868 * Row and column indices for elements. Must be in the range of [0..Size]
869 * unless the options EXPANDABLE or TRANSLATE are used. Node zero is the
870 * ground node. In no case may Node2 be less than zero.
871 * Template <output> (struct spTemplate *)
872 * Collection of pointers to four elements that are later used to directly
873 * address elements. User must supply the template, this routine will
874 * fill it.
875 *
876 * Possible errors:
877 * spNO_MEMORY
878 * Error is not cleared in this routine.
879 */
880
881 int
spGetAdmittance(Matrix,Node1,Node2,Template)882 spGetAdmittance( Matrix, Node1, Node2, Template )
883
884 char *Matrix;
885 int Node1, Node2;
886 struct spTemplate *Template;
887 {
888
889 /* Begin `spGetAdmittance'. */
890 Template->Element1 = spGetElement(Matrix, Node1, Node1 );
891 Template->Element2 = spGetElement(Matrix, Node2, Node2 );
892 Template->Element3Negated = spGetElement( Matrix, Node2, Node1 );
893 Template->Element4Negated = spGetElement( Matrix, Node1, Node2 );
894 if
895 ( (Template->Element1 == NULL)
896 OR (Template->Element2 == NULL)
897 OR (Template->Element3Negated == NULL)
898 OR (Template->Element4Negated == NULL)
899 ) return spNO_MEMORY;
900
901 if (Node1 == 0)
902 SWAP( RealNumber*, Template->Element1, Template->Element2 );
903
904 return spOKAY;
905 }
906 #endif /* QUAD_ELEMENT */
907
908
909
910
911
912
913
914
915
916 #if QUAD_ELEMENT
917 /*
918 * ADDITION OF FOUR ELEMENTS TO MATRIX BY INDEX
919 *
920 * Similar to spGetAdmittance, except that spGetAdmittance only
921 * handles 2-terminal components, whereas spGetQuad handles simple
922 * 4-terminals as well. These 4-terminals are simply generalized
923 * 2-terminals with the option of having the sense terminals different
924 * from the source and sink terminals. spGetQuad adds four
925 * elements to the matrix. Positive elements occur at Row1,Col1
926 * Row2,Col2 while negative elements occur at Row1,Col2 and Row2,Col1.
927 * The routine works fine if any of the rows and columns are zero.
928 * This routine is only to be used after spCreate() and before
929 * spMNA_Preorder(), spFactor() or spOrderAndFactor()
930 * unless TRANSLATE is set true.
931 *
932 * >>> Returns:
933 * Error code.
934 *
935 * >>> Arguments:
936 * Matrix <input> (char *)
937 * Pointer to the matrix that component is to be entered in.
938 * Row1 <input> (int)
939 * First row index for elements. Must be in the range of [0..Size]
940 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the
941 * ground row. In no case may Row1 be less than zero.
942 * Row2 <input> (int)
943 * Second row index for elements. Must be in the range of [0..Size]
944 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the
945 * ground row. In no case may Row2 be less than zero.
946 * Col1 <input> (int)
947 * First column index for elements. Must be in the range of [0..Size]
948 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the
949 * ground column. In no case may Col1 be less than zero.
950 * Col2 <input> (int)
951 * Second column index for elements. Must be in the range of [0..Size]
952 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the
953 * ground column. In no case may Col2 be less than zero.
954 * Template <output> (struct spTemplate *)
955 * Collection of pointers to four elements that are later used to directly
956 * address elements. User must supply the template, this routine will
957 * fill it.
958 * Real <input> (RealNumber)
959 * Real data to be added to elements.
960 * Imag <input> (RealNumber)
961 * Imag data to be added to elements. If matrix is real, this argument
962 * may be deleted.
963 *
964 * Possible errors:
965 * spNO_MEMORY
966 * Error is not cleared in this routine.
967 */
968
969 int
spGetQuad(Matrix,Row1,Row2,Col1,Col2,Template)970 spGetQuad( Matrix, Row1, Row2, Col1, Col2, Template )
971
972 char *Matrix;
973 int Row1, Row2, Col1, Col2;
974 struct spTemplate *Template;
975 {
976 /* Begin `spGetQuad'. */
977 Template->Element1 = spGetElement( Matrix, Row1, Col1);
978 Template->Element2 = spGetElement( Matrix, Row2, Col2 );
979 Template->Element3Negated = spGetElement( Matrix, Row2, Col1 );
980 Template->Element4Negated = spGetElement( Matrix, Row1, Col2 );
981 if
982 ( (Template->Element1 == NULL)
983 OR (Template->Element2 == NULL)
984 OR (Template->Element3Negated == NULL)
985 OR (Template->Element4Negated == NULL)
986 ) return spNO_MEMORY;
987
988 if (Template->Element1 == &((MatrixPtr)Matrix)->TrashCan.Real)
989 SWAP( RealNumber *, Template->Element1, Template->Element2 );
990
991 return spOKAY;
992 }
993 #endif /* QUAD_ELEMENT */
994
995
996
997
998
999
1000
1001
1002
1003 #if QUAD_ELEMENT
1004 /*
1005 * ADDITION OF FOUR STRUCTURAL ONES TO MATRIX BY INDEX
1006 *
1007 * Performs similar function to spGetQuad() except this routine is
1008 * meant for components that do not have an admittance representation.
1009 *
1010 * The following stamp is used:
1011 * Pos Neg Eqn
1012 * Pos [ . . 1 ]
1013 * Neg [ . . -1 ]
1014 * Eqn [ 1 -1 . ]
1015 *
1016 * >>> Returns:
1017 * Error code.
1018 *
1019 * >>> Arguments:
1020 * Matrix <input> (char *)
1021 * Pointer to the matrix that component is to be entered in.
1022 * Pos <input> (int)
1023 * See stamp above. Must be in the range of [0..Size]
1024 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the
1025 * ground row. In no case may Pos be less than zero.
1026 * Neg <input> (int)
1027 * See stamp above. Must be in the range of [0..Size]
1028 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the
1029 * ground row. In no case may Neg be less than zero.
1030 * Eqn <input> (int)
1031 * See stamp above. Must be in the range of [0..Size]
1032 * unless the options EXPANDABLE or TRANSLATE are used. Zero is the
1033 * ground row. In no case may Eqn be less than zero.
1034 * Template <output> (struct spTemplate *)
1035 * Collection of pointers to four elements that are later used to directly
1036 * address elements. User must supply the template, this routine will
1037 * fill it.
1038 *
1039 * Possible errors:
1040 * spNO_MEMORY
1041 * Error is not cleared in this routine.
1042 */
1043
1044 int
spGetOnes(Matrix,Pos,Neg,Eqn,Template)1045 spGetOnes(Matrix, Pos, Neg, Eqn, Template)
1046
1047 char *Matrix;
1048 int Pos, Neg, Eqn;
1049 struct spTemplate *Template;
1050 {
1051 /* Begin `spGetOnes'. */
1052 Template->Element4Negated = spGetElement( Matrix, Neg, Eqn );
1053 Template->Element3Negated = spGetElement( Matrix, Eqn, Neg );
1054 Template->Element2 = spGetElement( Matrix, Pos, Eqn );
1055 Template->Element1 = spGetElement( Matrix, Eqn, Pos );
1056 if
1057 ( (Template->Element1 == NULL)
1058 OR (Template->Element2 == NULL)
1059 OR (Template->Element3Negated == NULL)
1060 OR (Template->Element4Negated == NULL)
1061 ) return spNO_MEMORY;
1062
1063 spADD_REAL_QUAD( *Template, 1.0 );
1064 return spOKAY;
1065 }
1066 #endif /* QUAD_ELEMENT */
1067
1068
1069
1070
1071
1072
1073
1074 /*
1075 *
1076 * CREATE AND SPLICE ELEMENT INTO MATRIX
1077 *
1078 * This routine is used to create new matrix elements and splice them into the
1079 * matrix.
1080 *
1081 * >>> Returned:
1082 * A pointer to the element that was created is returned.
1083 *
1084 * >>> Arguments:
1085 * Matrix <input> (MatrixPtr)
1086 * Pointer to matrix.
1087 * Row <input> (int)
1088 * Row index for element.
1089 * Col <input> (int)
1090 * Column index for element.
1091 * LastAddr <input-output> (ElementPtr *)
1092 * This contains the address of the pointer to the element just above the
1093 * one being created. It is used to speed the search and it is updated with
1094 * address of the created element.
1095 * Fillin <input> (BOOLEAN)
1096 * Flag that indicates if created element is to be a fill-in.
1097 *
1098 * >>> Local variables:
1099 * pElement (ElementPtr)
1100 * Pointer to an element in the matrix. It is used to refer to the newly
1101 * created element and to restring the pointers of the element's row and
1102 * column.
1103 * pLastElement (ElementPtr)
1104 * Pointer to the element in the matrix that was just previously pointed
1105 * to by pElement. It is used to restring the pointers of the element's
1106 * row and column.
1107 * pCreatedElement (ElementPtr)
1108 * Pointer to the desired element, the one that was just created.
1109 *
1110 * >>> Possible errors:
1111 * spNO_MEMORY
1112 */
1113
1114 ElementPtr
spcCreateElement(Matrix,Row,Col,LastAddr,Fillin)1115 spcCreateElement( Matrix, Row, Col, LastAddr, Fillin )
1116
1117 MatrixPtr Matrix;
1118 int Row;
1119 register int Col;
1120 register ElementPtr *LastAddr;
1121 BOOLEAN Fillin;
1122 {
1123 register ElementPtr pElement, pLastElement;
1124 ElementPtr pCreatedElement, spcGetElement(), spcGetFillin();
1125
1126 /* Begin `spcCreateElement'. */
1127
1128 if (Matrix->RowsLinked)
1129 {
1130 /* Row pointers cannot be ignored. */
1131 if (Fillin)
1132 { pElement = spcGetFillin( Matrix );
1133 Matrix->Fillins++;
1134 }
1135 else
1136 { pElement = spcGetElement( Matrix );
1137 Matrix->NeedsOrdering = YES;
1138 }
1139 if (pElement == NULL) return NULL;
1140
1141 /* If element is on diagonal, store pointer in Diag. */
1142 if (Row == Col) Matrix->Diag[Row] = pElement;
1143
1144 /* Initialize Element. */
1145 pCreatedElement = pElement;
1146 pElement->Row = Row;
1147 pElement->Col = Col;
1148 pElement->Real = 0.0;
1149 #if spCOMPLEX
1150 pElement->Imag = 0.0;
1151 #endif
1152 #if INITIALIZE
1153 pElement->pInitInfo = NULL;
1154 #endif
1155
1156 /* Splice element into column. */
1157 pElement->NextInCol = *LastAddr;
1158 *LastAddr = pElement;
1159
1160 /* Search row for proper element position. */
1161 pElement = Matrix->FirstInRow[Row];
1162 pLastElement = NULL;
1163 while (pElement != NULL)
1164 {
1165 /* Search for element row position. */
1166 if (pElement->Col < Col)
1167 {
1168 /* Have not reached desired element. */
1169 pLastElement = pElement;
1170 pElement = pElement->NextInRow;
1171 }
1172 else pElement = NULL;
1173 }
1174
1175 /* Splice element into row. */
1176 pElement = pCreatedElement;
1177 if (pLastElement == NULL)
1178 {
1179 /* Element is first in row. */
1180 pElement->NextInRow = Matrix->FirstInRow[Row];
1181 Matrix->FirstInRow[Row] = pElement;
1182 }
1183 else
1184 /* Element is not first in row. */
1185 {
1186 pElement->NextInRow = pLastElement->NextInRow;
1187 pLastElement->NextInRow = pElement;
1188 }
1189
1190 }
1191 else
1192 {
1193 /*
1194 * Matrix has not been factored yet. Thus get element rather than fill-in.
1195 * Also, row pointers can be ignored.
1196 */
1197
1198 /* Allocate memory for Element. */
1199 pElement = spcGetElement( Matrix );
1200 if (pElement == NULL) return NULL;
1201
1202 /* If element is on diagonal, store pointer in Diag. */
1203 if (Row == Col) Matrix->Diag[Row] = pElement;
1204
1205 /* Initialize Element. */
1206 pCreatedElement = pElement;
1207 pElement->Row = Row;
1208 #if DEBUG
1209 pElement->Col = Col;
1210 #endif
1211 pElement->Real = 0.0;
1212 #if spCOMPLEX
1213 pElement->Imag = 0.0;
1214 #endif
1215 #if INITIALIZE
1216 pElement->pInitInfo = NULL;
1217 #endif
1218
1219 /* Splice element into column. */
1220 pElement->NextInCol = *LastAddr;
1221 *LastAddr = pElement;
1222 }
1223
1224 Matrix->Elements++;
1225 return pCreatedElement;
1226 }
1227
1228
1229
1230
1231
1232
1233
1234
1235 /*
1236 *
1237 * LINK ROWS
1238 *
1239 * This routine is used to generate the row links. The spGetElement()
1240 * routines do not create row links, which are needed by the spFactor()
1241 * routines.
1242 *
1243 * >>> Arguments:
1244 * Matrix <input> (MatrixPtr)
1245 * Pointer to the matrix.
1246 *
1247 * >>> Local variables:
1248 * pElement (ElementPtr)
1249 * Pointer to an element in the matrix.
1250 * FirstInRowEntry (ElementPtr *)
1251 * A pointer into the FirstInRow array. Points to the FirstInRow entry
1252 * currently being operated upon.
1253 * FirstInRowArray (ArrayOfElementPtrs)
1254 * A pointer to the FirstInRow array. Same as Matrix->FirstInRow but
1255 * resides in a register and requires less indirection so is faster to
1256 * use.
1257 * Col (int)
1258 * Column currently being operated upon.
1259 */
1260
1261 void
spcLinkRows(Matrix)1262 spcLinkRows( Matrix )
1263
1264 MatrixPtr Matrix;
1265 {
1266 register ElementPtr pElement, *FirstInRowEntry;
1267 register ArrayOfElementPtrs FirstInRowArray;
1268 register int Col;
1269
1270 /* Begin `spcLinkRows'. */
1271 FirstInRowArray = Matrix->FirstInRow;
1272 for (Col = Matrix->Size; Col >= 1; Col--)
1273 {
1274 /* Generate row links for the elements in the Col'th column. */
1275 pElement = Matrix->FirstInCol[Col];
1276
1277 while (pElement != NULL)
1278 { pElement->Col = Col;
1279 FirstInRowEntry = &FirstInRowArray[pElement->Row];
1280 pElement->NextInRow = *FirstInRowEntry;
1281 *FirstInRowEntry = pElement;
1282 pElement = pElement->NextInCol;
1283 }
1284 }
1285 Matrix->RowsLinked = YES;
1286 return;
1287 }
1288
1289
1290
1291
1292
1293
1294
1295
1296 /*
1297 * ENLARGE MATRIX
1298 *
1299 * Increases the size of the matrix.
1300 *
1301 * >>> Arguments:
1302 * Matrix <input> (MatrixPtr)
1303 * Pointer to the matrix.
1304 * NewSize <input> (int)
1305 * The new size of the matrix.
1306 *
1307 * >>> Local variables:
1308 * OldAllocatedSize (int)
1309 * The allocated size of the matrix before it is expanded.
1310 */
1311
1312 static void
EnlargeMatrix(Matrix,NewSize)1313 EnlargeMatrix( Matrix, NewSize )
1314
1315 MatrixPtr Matrix;
1316 register int NewSize;
1317 {
1318 register int I, OldAllocatedSize = Matrix->AllocatedSize;
1319
1320 /* Begin `EnlargeMatrix'. */
1321 Matrix->Size = NewSize;
1322
1323 if (NewSize <= OldAllocatedSize)
1324 return;
1325
1326 /* Expand the matrix frame. */
1327 NewSize = MAX( NewSize, EXPANSION_FACTOR * OldAllocatedSize );
1328 Matrix->AllocatedSize = NewSize;
1329
1330 if (( REALLOC(Matrix->IntToExtColMap, int, NewSize+1)) == NULL)
1331 { Matrix->Error = spNO_MEMORY;
1332 return;
1333 }
1334 if (( REALLOC(Matrix->IntToExtRowMap, int, NewSize+1)) == NULL)
1335 { Matrix->Error = spNO_MEMORY;
1336 return;
1337 }
1338 if (( REALLOC(Matrix->Diag, ElementPtr, NewSize+1)) == NULL)
1339 { Matrix->Error = spNO_MEMORY;
1340 return;
1341 }
1342 if (( REALLOC(Matrix->FirstInCol, ElementPtr, NewSize+1)) == NULL)
1343 { Matrix->Error = spNO_MEMORY;
1344 return;
1345 }
1346 if (( REALLOC(Matrix->FirstInRow, ElementPtr, NewSize+1)) == NULL)
1347 { Matrix->Error = spNO_MEMORY;
1348 return;
1349 }
1350
1351 /*
1352 * Destroy the Markowitz and Intermediate vectors, they will be recreated
1353 * in spOrderAndFactor().
1354 */
1355 FREE( Matrix->MarkowitzRow );
1356 FREE( Matrix->MarkowitzCol );
1357 FREE( Matrix->MarkowitzProd );
1358 FREE( Matrix->DoRealDirect );
1359 FREE( Matrix->DoCmplxDirect );
1360 FREE( Matrix->Intermediate );
1361 Matrix->InternalVectorsAllocated = NO;
1362
1363 /* Initialize the new portion of the vectors. */
1364 for (I = OldAllocatedSize+1; I <= NewSize; I++)
1365 { Matrix->IntToExtColMap[I] = I;
1366 Matrix->IntToExtRowMap[I] = I;
1367 Matrix->Diag[I] = NULL;
1368 Matrix->FirstInRow[I] = NULL;
1369 Matrix->FirstInCol[I] = NULL;
1370 }
1371
1372 return;
1373 }
1374
1375
1376
1377
1378
1379
1380
1381
1382 #if TRANSLATE
1383
1384 /*
1385 * EXPAND TRANSLATION ARRAYS
1386 *
1387 * Increases the size arrays that are used to translate external to internal
1388 * row and column numbers.
1389 *
1390 * >>> Arguments:
1391 * Matrix <input> (MatrixPtr)
1392 * Pointer to the matrix.
1393 * NewSize <input> (int)
1394 * The new size of the translation arrays.
1395 *
1396 * >>> Local variables:
1397 * OldAllocatedSize (int)
1398 * The allocated size of the translation arrays before being expanded.
1399 */
1400
1401 static void
ExpandTranslationArrays(Matrix,NewSize)1402 ExpandTranslationArrays( Matrix, NewSize )
1403
1404 MatrixPtr Matrix;
1405 register int NewSize;
1406 {
1407 register int I, OldAllocatedSize = Matrix->AllocatedExtSize;
1408
1409 /* Begin `ExpandTranslationArrays'. */
1410 Matrix->ExtSize = NewSize;
1411
1412 if (NewSize <= OldAllocatedSize)
1413 return;
1414
1415 /* Expand the translation arrays ExtToIntRowMap and ExtToIntColMap. */
1416 NewSize = MAX( NewSize, EXPANSION_FACTOR * OldAllocatedSize );
1417 Matrix->AllocatedExtSize = NewSize;
1418
1419 if (( REALLOC(Matrix->ExtToIntRowMap, int, NewSize+1)) == NULL)
1420 { Matrix->Error = spNO_MEMORY;
1421 return;
1422 }
1423 if (( REALLOC(Matrix->ExtToIntColMap, int, NewSize+1)) == NULL)
1424 { Matrix->Error = spNO_MEMORY;
1425 return;
1426 }
1427
1428 /* Initialize the new portion of the vectors. */
1429 for (I = OldAllocatedSize+1; I <= NewSize; I++)
1430 { Matrix->ExtToIntRowMap[I] = -1;
1431 Matrix->ExtToIntColMap[I] = -1;
1432 }
1433
1434 return;
1435 }
1436 #endif
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446 #if INITIALIZE
1447 /*
1448 * INITIALIZE MATRIX
1449 *
1450 * With the INITIALIZE compiler option (see spConfig.h) set true,
1451 * Sparse allows the user to keep initialization information with each
1452 * structurally nonzero matrix element. Each element has a pointer
1453 * that is set and used by the user. The user can set this pointer
1454 * using spInstallInitInfo and may be read using spGetInitInfo. Both
1455 * may be used only after the element exists. The function
1456 * spInitialize() is a user customizable way to initialize the matrix.
1457 * Passed to this routine is a function pointer. spInitialize() sweeps
1458 * through every element in the matrix and checks the pInitInfo
1459 * pointer (the user supplied pointer). If the pInitInfo is NULL,
1460 * which is true unless the user changes it (almost always true for
1461 * fill-ins), then the element is zeroed. Otherwise, the function
1462 * pointer is called and passed the pInitInfo pointer as well as the
1463 * element pointer and the external row and column numbers. If the
1464 * user sets the value of each element, then spInitialize() replaces
1465 * spClear().
1466 *
1467 * The user function is expected to return a nonzero integer if there
1468 * is a fatal error and zero otherwise. Upon encountering a nonzero
1469 * return code, spInitialize() terminates and returns the error code.
1470 *
1471 * >>> Arguments:
1472 * Matrix <input> (char *)
1473 * Pointer to matrix.
1474 *
1475 * >>> Possible Errors:
1476 * Returns nonzero if error, zero otherwise.
1477 */
1478
1479 void
spInstallInitInfo(pElement,pInitInfo)1480 spInstallInitInfo( pElement, pInitInfo )
1481
1482 RealNumber *pElement;
1483 char *pInitInfo;
1484 {
1485 /* Begin `spInstallInitInfo'. */
1486 ASSERT(pElement != NULL);
1487
1488 ((ElementPtr)pElement)->pInitInfo = pInitInfo;
1489 }
1490
1491
1492 char *
spGetInitInfo(pElement)1493 spGetInitInfo( pElement )
1494
1495 RealNumber *pElement;
1496 {
1497 /* Begin `spGetInitInfo'. */
1498 ASSERT(pElement != NULL);
1499
1500 return (char *)((ElementPtr)pElement)->pInitInfo;
1501 }
1502
1503
1504 int
spInitialize(eMatrix,pInit)1505 spInitialize( eMatrix, pInit )
1506
1507 char *eMatrix;
1508 int (*pInit)();
1509 {
1510 MatrixPtr Matrix = (MatrixPtr)eMatrix;
1511 register ElementPtr pElement;
1512 int J, Error, Col;
1513
1514 /* Begin `spInitialize'. */
1515 ASSERT( IS_SPARSE( Matrix ) );
1516
1517 #if spCOMPLEX
1518 /* Clear imaginary part of matrix if matrix is real but was complex. */
1519 if (Matrix->PreviousMatrixWasComplex AND NOT Matrix->Complex)
1520 { for (J = Matrix->Size; J > 0; J--)
1521 { pElement = Matrix->FirstInCol[J];
1522 while (pElement != NULL)
1523 { pElement->Imag = 0.0;
1524 pElement = pElement->NextInCol;
1525 }
1526 }
1527 }
1528 #endif /* spCOMPLEX */
1529
1530 /* Initialize the matrix. */
1531 for (J = Matrix->Size; J > 0; J--)
1532 { pElement = Matrix->FirstInCol[J];
1533 Col = Matrix->IntToExtColMap[J];
1534 while (pElement != NULL)
1535 { if (pElement->pInitInfo == NULL)
1536 { pElement->Real = 0.0;
1537 # if spCOMPLEX
1538 pElement->Imag = 0.0;
1539 # endif
1540 }
1541 else
1542 { Error = (*pInit)((RealNumber *)pElement, pElement->pInitInfo,
1543 Matrix->IntToExtRowMap[pElement->Row], Col);
1544 if (Error)
1545 { Matrix->Error = spFATAL;
1546 return Error;
1547 }
1548
1549 }
1550 pElement = pElement->NextInCol;
1551 }
1552 }
1553
1554 /* Empty the trash. */
1555 Matrix->TrashCan.Real = 0.0;
1556 #if spCOMPLEX
1557 Matrix->TrashCan.Imag = 0.0;
1558 #endif
1559
1560 Matrix->Error = spOKAY;
1561 Matrix->Factored = NO;
1562 Matrix->SingularCol = 0;
1563 Matrix->SingularRow = 0;
1564 Matrix->PreviousMatrixWasComplex = Matrix->Complex;
1565 return 0;
1566 }
1567 #endif /* INITIALIZE */
1568