1 /*
2 * MATRIX UTILITY MODULE
3 *
4 * Author: Advising professor:
5 * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
6 * UC Berkeley
7 */
8 /*! \file
9 * This file contains various optional utility routines.
10 *
11 * Objects that begin with the \a spc prefix are considered private
12 * and should not be used.
13 *
14 * \author
15 * Kenneth S. Kundert <kundert@users.sourceforge.net>
16 */
17 /* >>> User accessible functions contained in this file:
18 * spMNA_Preorder
19 * spScale
20 * spMultiply
21 * spMultTransposed
22 * spDeterminant
23 * spStrip
24 * spDeleteRowAndCol
25 * spPseudoCondition
26 * spCondition
27 * spNorm
28 * spLargestElement
29 * spRoundoff
30 * spErrorMessage
31 *
32 * >>> Other functions contained in this file:
33 * CountTwins
34 * SwapCols
35 * ScaleComplexMatrix
36 * ComplexMatrixMultiply
37 * ComplexCondition
38 */
39
40
41 /*
42 * Revision and copyright information.
43 *
44 * Copyright (c) 1985-2003 by Kenneth S. Kundert
45 */
46
47 /*
48 * IMPORTS
49 *
50 * >>> Import descriptions:
51 * spConfig.h
52 * Macros that customize the sparse matrix routines.
53 * spMatrix.h
54 * Macros and declarations to be imported by the user.
55 * spDefs.h
56 * Matrix type and macro definitions for the sparse matrix routines.
57 */
58
59 #define spINSIDE_SPARSE
60 #include <stdio.h>
61 #include "spConfig.h"
62 #include "spMatrix.h"
63 #include "spDefs.h"
64
65
66
67
68
69 /*
70 * Function declarations
71 */
72
73 static int CountTwins( MatrixPtr, int, ElementPtr*, ElementPtr* );
74 static void SwapCols( MatrixPtr, ElementPtr, ElementPtr );
75 #if spCOMPLEX
76 static void ScaleComplexMatrix( MatrixPtr, RealVector, RealVector );
77 #endif
78 #if spSEPARATED_COMPLEX_VECTORS
79 static void ComplexMatrixMultiply( MatrixPtr,
80 RealVector, RealVector, RealVector, RealVector );
81 static void ComplexTransposedMatrixMultiply( MatrixPtr,
82 RealVector, RealVector, RealVector, RealVector );
83 #elseif spCOMPLEX
84 static void ComplexMatrixMultiply( MatrixPtr,
85 RealVector, RealVector );
86 static void ComplexTransposedMatrixMultiply( MatrixPtr,
87 RealVector, RealVector );
88 #endif
89 #if spCOMPLEX
90 static RealNumber ComplexCondition( MatrixPtr, RealNumber, int* );
91 #endif
92
93
94
95
96
97 #if MODIFIED_NODAL
98 /*!
99 * This routine massages modified node admittance matrices to remove
100 * zeros from the diagonal. It takes advantage of the fact that the
101 * row and column associated with a zero diagonal usually have
102 * structural ones placed symmetricly. This routine should be used
103 * only on modified node admittance matrices and should be executed
104 * after the matrix has been built but before the factorization
105 * begins. It should be executed for the initial factorization only
106 * and should be executed before the rows have been linked. Thus it
107 * should be run before using spScale(), spMultiply(),
108 * spDeleteRowAndCol(), or spNorm().
109 *
110 * This routine exploits the fact that the structural ones are placed
111 * in the matrix in symmetric twins. For example, the stamps for
112 * grounded and a floating voltage sources are \code
113 * grounded: floating:
114 * [ x x 1 ] [ x x 1 ]
115 * [ x x ] [ x x -1 ]
116 * [ 1 ] [ 1 -1 ]
117 * \endcode
118 * Notice for the grounded source, there is one set of twins, and for
119 * the floating, there are two sets. We remove the zero from the diagonal
120 * by swapping the rows associated with a set of twins. For example: \code
121 * grounded: floating 1: floating 2:
122 * [ 1 ] [ 1 -1 ] [ x x 1 ]
123 * [ x x ] [ x x -1 ] [ 1 -1 ]
124 * [ x x 1 ] [ x x 1 ] [ x x -1 ]
125 * \endcode
126 *
127 * It is important to deal with any zero diagonals that only have one
128 * set of twins before dealing with those that have more than one because
129 * swapping row destroys the symmetry of any twins in the rows being
130 * swapped, which may limit future moves. Consider \code
131 * [ x x 1 ]
132 * [ x x -1 1 ]
133 * [ 1 -1 ]
134 * [ 1 ]
135 * \endcode
136 * There is one set of twins for diagonal 4 and two for diagonal 3.
137 * Dealing with diagonal 4 first requires swapping rows 2 and 4. \code
138 * [ x x 1 ]
139 * [ 1 ]
140 * [ 1 -1 ]
141 * [ x x -1 1 ]
142 * \endcode
143 * We can now deal with diagonal 3 by swapping rows 1 and 3. \code
144 * [ 1 -1 ]
145 * [ 1 ]
146 * [ x x 1 ]
147 * [ x x -1 1 ]
148 * \endcode
149 * And we are done, there are no zeros left on the diagonal. However, if
150 * we originally dealt with diagonal 3 first, we could swap rows 2 and 3 \code
151 * [ x x 1 ]
152 * [ 1 -1 ]
153 * [ x x -1 1 ]
154 * [ 1 ]
155 * \endcode
156 * Diagonal 4 no longer has a symmetric twin and we cannot continue.
157 *
158 * So we always take care of lone twins first. When none remain, we
159 * choose arbitrarily a set of twins for a diagonal with more than one set
160 * and swap the rows corresponding to that twin. We then deal with any
161 * lone twins that were created and repeat the procedure until no
162 * zero diagonals with symmetric twins remain.
163 *
164 * In this particular implementation, columns are swapped rather than rows.
165 * The algorithm used in this function was developed by Ken Kundert and
166 * Tom Quarles.
167 *
168 * \param * eMatrix
169 * Pointer to the matrix to be preordered.
170 */
171 /* >>> Local variables;
172 * J (int)
173 * Column with zero diagonal being currently considered.
174 * pTwin1 (ElementPtr)
175 * Pointer to the twin found in the column belonging to the zero diagonal.
176 * pTwin2 (ElementPtr)
177 * Pointer to the twin found in the row belonging to the zero diagonal.
178 * belonging to the zero diagonal.
179 * AnotherPassNeeded (BOOLEAN)
180 * Flag indicating that at least one zero diagonal with symmetric twins
181 * remain.
182 * StartAt (int)
183 * Column number of first zero diagonal with symmetric twins.
184 * Swapped (BOOLEAN)
185 * Flag indicating that columns were swapped on this pass.
186 * Twins (int)
187 * Number of symmetric twins corresponding to current zero diagonal.
188 */
189
190 void
spMNA_Preorder(spMatrix eMatrix)191 spMNA_Preorder( spMatrix eMatrix )
192 {
193 MatrixPtr Matrix = (MatrixPtr)eMatrix;
194 register int J, Size;
195 ElementPtr pTwin1=0, pTwin2=0;
196 int Twins, StartAt = 1;
197 BOOLEAN Swapped, AnotherPassNeeded;
198
199 /* Begin `spMNA_Preorder'. */
200 ASSERT_IS_SPARSE( Matrix );
201 ASSERT_NO_ERRORS( Matrix );
202 ASSERT_IS_NOT_FACTORED( Matrix );
203
204 if (Matrix->RowsLinked) return;
205 Size = Matrix->Size;
206 Matrix->Reordered = YES;
207
208 do
209 { AnotherPassNeeded = Swapped = NO;
210
211 /* Search for zero diagonals with lone twins. */
212 for (J = StartAt; J <= Size; J++)
213 { if (Matrix->Diag[J] == NULL)
214 { Twins = CountTwins( Matrix, J, &pTwin1, &pTwin2 );
215 if (Twins == 1)
216 { /* Lone twins found, swap rows. */
217 SwapCols( Matrix, pTwin1, pTwin2 );
218 Swapped = YES;
219 }
220 else if ((Twins > 1) AND NOT AnotherPassNeeded)
221 { AnotherPassNeeded = YES;
222 StartAt = J;
223 }
224 }
225 }
226
227 /* All lone twins are gone, look for zero diagonals with multiple twins. */
228 if (AnotherPassNeeded)
229 { for (J = StartAt; NOT Swapped AND (J <= Size); J++)
230 { if (Matrix->Diag[J] == NULL)
231 { // Twins = CountTwins( Matrix, J, &pTwin1, &pTwin2 );
232 CountTwins( Matrix, J, &pTwin1, &pTwin2 );
233 SwapCols( Matrix, pTwin1, pTwin2 );
234 Swapped = YES;
235 }
236 }
237 }
238 } while (AnotherPassNeeded);
239 return;
240 }
241
242
243
244
245 /*
246 * COUNT TWINS
247 *
248 * This function counts the number of symmetric twins associated with
249 * a zero diagonal and returns one set of twins if any exist. The
250 * count is terminated early at two.
251 */
252
253 static int
CountTwins(MatrixPtr Matrix,int Col,ElementPtr * ppTwin1,ElementPtr * ppTwin2)254 CountTwins(
255 MatrixPtr Matrix,
256 int Col,
257 ElementPtr *ppTwin1,
258 ElementPtr *ppTwin2
259 )
260 {
261 int Row, Twins = 0;
262 ElementPtr pTwin1, pTwin2;
263
264 /* Begin `CountTwins'. */
265
266 pTwin1 = Matrix->FirstInCol[Col];
267 while (pTwin1 != NULL)
268 { if (ABS(pTwin1->Real) == 1.0)
269 { Row = pTwin1->Row;
270 pTwin2 = Matrix->FirstInCol[Row];
271 while ((pTwin2 != NULL) AND (pTwin2->Row != Col))
272 pTwin2 = pTwin2->NextInCol;
273 if ((pTwin2 != NULL) AND (ABS(pTwin2->Real) == 1.0))
274 { /* Found symmetric twins. */
275 if (++Twins >= 2) return Twins;
276 (*ppTwin1 = pTwin1)->Col = Col;
277 (*ppTwin2 = pTwin2)->Col = Row;
278 }
279 }
280 pTwin1 = pTwin1->NextInCol;
281 }
282 return Twins;
283 }
284
285
286
287
288 /*
289 * SWAP COLUMNS
290 *
291 * This function swaps two columns and is applicable before the rows are
292 * linked.
293 */
294
295 static void
SwapCols(MatrixPtr Matrix,ElementPtr pTwin1,ElementPtr pTwin2)296 SwapCols(
297 MatrixPtr Matrix,
298 ElementPtr pTwin1,
299 ElementPtr pTwin2
300 )
301 {
302 int Col1 = pTwin1->Col, Col2 = pTwin2->Col;
303
304 /* Begin `SwapCols'. */
305
306 SWAP (ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]);
307 SWAP (int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]);
308 #if TRANSLATE
309 Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col2]]=Col2;
310 Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col1]]=Col1;
311 #endif
312
313 Matrix->Diag[Col1] = pTwin2;
314 Matrix->Diag[Col2] = pTwin1;
315 Matrix->NumberOfInterchangesIsOdd = NOT Matrix->NumberOfInterchangesIsOdd;
316 return;
317 }
318 #endif /* MODIFIED_NODAL */
319
320
321
322
323
324
325
326
327
328 #if SCALING
329 /*!
330 * This function scales the matrix to enhance the possibility of
331 * finding a good pivoting order. Note that scaling enhances accuracy
332 * of the solution only if it affects the pivoting order, so it makes
333 * no sense to scale the matrix before spFactor(). If scaling is
334 * desired it should be done before spOrderAndFactor(). There
335 * are several things to take into account when choosing the scale
336 * factors. First, the scale factors are directly multiplied against
337 * the elements in the matrix. To prevent roundoff, each scale factor
338 * should be equal to an integer power of the number base of the
339 * machine. Since most machines operate in base two, scale factors
340 * should be a power of two. Second, the matrix should be scaled such
341 * that the matrix of element uncertainties is equilibrated. Third,
342 * this function multiplies the scale factors by the elements, so if
343 * one row tends to have uncertainties 1000 times smaller than the
344 * other rows, then its scale factor should be 1024, not 1/1024.
345 * Fourth, to save time, this function does not scale rows or columns
346 * if their scale factors are equal to one. Thus, the scale factors
347 * should be normalized to the most common scale factor. Rows and
348 * columns should be normalized separately. For example, if the size
349 * of the matrix is 100 and 10 rows tend to have uncertainties near
350 * 1e-6 and the remaining 90 have uncertainties near 1e-12, then the
351 * scale factor for the 10 should be 1/1,048,576 and the scale factors
352 * for the remaining 90 should be 1. Fifth, since this routine
353 * directly operates on the matrix, it is necessary to apply the scale
354 * factors to the RHS and Solution vectors. It may be easier to
355 * simply use spOrderAndFactor() on a scaled matrix to choose the
356 * pivoting order, and then throw away the matrix. Subsequent
357 * factorizations, performed with spFactor(), will not need to have
358 * the RHS and Solution vectors descaled. Lastly, this function
359 * should not be executed before the function spMNA_Preorder().
360 *
361 * \param eMatrix
362 * Pointer to the matrix to be scaled.
363 * \param SolutionScaleFactors
364 * The array of Solution scale factors. These factors scale the columns.
365 * All scale factors are real valued.
366 * \param RHS_ScaleFactors
367 * The array of RHS scale factors. These factors scale the rows.
368 * All scale factors are real valued.
369 */
370 /* >>> Local variables:
371 * lSize (int)
372 * Local version of the size of the matrix.
373 * pElement (ElementPtr)
374 * Pointer to an element in the matrix.
375 * pExtOrder (int *)
376 * Pointer into either IntToExtRowMap or IntToExtColMap vector. Used to
377 * compensate for any row or column swaps that have been performed.
378 * ScaleFactor (RealNumber)
379 * The scale factor being used on the current row or column.
380 */
381
382 void
spScale(spMatrix eMatrix,spREAL RHS_ScaleFactors[],spREAL SolutionScaleFactors[])383 spScale(
384 spMatrix eMatrix,
385 spREAL RHS_ScaleFactors[],
386 spREAL SolutionScaleFactors[]
387 )
388 {
389 MatrixPtr Matrix = (MatrixPtr)eMatrix;
390 register ElementPtr pElement;
391 register int I, lSize, *pExtOrder;
392 RealNumber ScaleFactor;
393 void ScaleComplexMatrix();
394
395 /* Begin `spScale'. */
396 ASSERT_IS_SPARSE( Matrix );
397 ASSERT_NO_ERRORS( Matrix );
398 ASSERT_IS_NOT_FACTORED( Matrix );
399 if (NOT Matrix->RowsLinked) spcLinkRows( Matrix );
400
401 #if spCOMPLEX
402 if (Matrix->Complex)
403 { ScaleComplexMatrix( Matrix, RHS_ScaleFactors, SolutionScaleFactors );
404 return;
405 }
406 #endif
407
408 #if REAL
409 lSize = Matrix->Size;
410
411 /* Correct pointers to arrays for ARRAY_OFFSET */
412 #if NOT ARRAY_OFFSET
413 --RHS_ScaleFactors;
414 --SolutionScaleFactors;
415 #endif
416
417 /* Scale Rows */
418 pExtOrder = &Matrix->IntToExtRowMap[1];
419 for (I = 1; I <= lSize; I++)
420 { if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0)
421 { pElement = Matrix->FirstInRow[I];
422
423 while (pElement != NULL)
424 { pElement->Real *= ScaleFactor;
425 pElement = pElement->NextInRow;
426 }
427 }
428 }
429
430 /* Scale Columns */
431 pExtOrder = &Matrix->IntToExtColMap[1];
432 for (I = 1; I <= lSize; I++)
433 { if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0)
434 { pElement = Matrix->FirstInCol[I];
435
436 while (pElement != NULL)
437 { pElement->Real *= ScaleFactor;
438 pElement = pElement->NextInCol;
439 }
440 }
441 }
442 return;
443
444 #endif /* REAL */
445 }
446 #endif /* SCALING */
447
448
449
450
451
452
453
454
455
456 #if spCOMPLEX AND SCALING
457 /*
458 * SCALE COMPLEX MATRIX
459 *
460 * This function scales the matrix to enhance the possibility of
461 * finding a good pivoting order. Note that scaling enhances accuracy
462 * of the solution only if it affects the pivoting order, so it makes
463 * no sense to scale the matrix before spFactor(). If scaling is
464 * desired it should be done before spOrderAndFactor(). There
465 * are several things to take into account when choosing the scale
466 * factors. First, the scale factors are directly multiplied against
467 * the elements in the matrix. To prevent roundoff, each scale factor
468 * should be equal to an integer power of the number base of the
469 * machine. Since most machines operate in base two, scale factors
470 * should be a power of two. Second, the matrix should be scaled such
471 * that the matrix of element uncertainties is equilibrated. Third,
472 * this function multiplies the scale factors by the elements, so if
473 * one row tends to have uncertainties 1000 times smaller than the
474 * other rows, then its scale factor should be 1024, not 1/1024.
475 * Fourth, to save time, this function does not scale rows or columns
476 * if their scale factors are equal to one. Thus, the scale factors
477 * should be normalized to the most common scale factor. Rows and
478 * columns should be normalized separately. For example, if the size
479 * of the matrix is 100 and 10 rows tend to have uncertainties near
480 * 1e-6 and the remaining 90 have uncertainties near 1e-12, then the
481 * scale factor for the 10 should be 1/1,048,576 and the scale factors
482 * for the remaining 90 should be 1. Fifth, since this routine
483 * directly operates on the matrix, it is necessary to apply the scale
484 * factors to the RHS and Solution vectors. It may be easier to
485 * simply use spOrderAndFactor() on a scaled matrix to choose the
486 * pivoting order, and then throw away the matrix. Subsequent
487 * factorizations, performed with spFactor(), will not need to have
488 * the RHS and Solution vectors descaled. Lastly, this function
489 * should not be executed before the function spMNA_Preorder.
490 *
491 * >>> Arguments:
492 * Matrix <input> (char *)
493 * Pointer to the matrix to be scaled.
494 * SolutionScaleFactors <input> (RealVector)
495 * The array of Solution scale factors. These factors scale the columns.
496 * All scale factors are real valued.
497 * RHS_ScaleFactors <input> (RealVector)
498 * The array of RHS scale factors. These factors scale the rows.
499 * All scale factors are real valued.
500 *
501 * >>> Local variables:
502 * lSize (int)
503 * Local version of the size of the matrix.
504 * pElement (ElementPtr)
505 * Pointer to an element in the matrix.
506 * pExtOrder (int *)
507 * Pointer into either IntToExtRowMap or IntToExtColMap vector. Used to
508 * compensate for any row or column swaps that have been performed.
509 * ScaleFactor (RealNumber)
510 * The scale factor being used on the current row or column.
511 */
512
513 static void
ScaleComplexMatrix(MatrixPtr Matrix,register RealVector RHS_ScaleFactors,register RealVector SolutionScaleFactors)514 ScaleComplexMatrix(
515 MatrixPtr Matrix,
516 register RealVector RHS_ScaleFactors,
517 register RealVector SolutionScaleFactors
518 )
519 {
520 register ElementPtr pElement;
521 register int I, lSize, *pExtOrder;
522 RealNumber ScaleFactor;
523
524 /* Begin `ScaleComplexMatrix'. */
525 lSize = Matrix->Size;
526
527 /* Correct pointers to arrays for ARRAY_OFFSET */
528 #if NOT ARRAY_OFFSET
529 --RHS_ScaleFactors;
530 --SolutionScaleFactors;
531 #endif
532
533 /* Scale Rows */
534 pExtOrder = &Matrix->IntToExtRowMap[1];
535 for (I = 1; I <= lSize; I++)
536 { if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0)
537 { pElement = Matrix->FirstInRow[I];
538
539 while (pElement != NULL)
540 { pElement->Real *= ScaleFactor;
541 pElement->Imag *= ScaleFactor;
542 pElement = pElement->NextInRow;
543 }
544 }
545 }
546
547 /* Scale Columns */
548 pExtOrder = &Matrix->IntToExtColMap[1];
549 for (I = 1; I <= lSize; I++)
550 { if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0)
551 { pElement = Matrix->FirstInCol[I];
552
553 while (pElement != NULL)
554 { pElement->Real *= ScaleFactor;
555 pElement->Imag *= ScaleFactor;
556 pElement = pElement->NextInCol;
557 }
558 }
559 }
560 return;
561 }
562 #endif /* SCALING AND spCOMPLEX */
563
564
565
566
567
568
569
570
571 #if MULTIPLICATION
572 /*!
573 * Multiplies matrix by solution vector to find source vector.
574 * Assumes matrix has not been factored. This routine can be used
575 * as a test to see if solutions are correct. It should not be used
576 * before spMNA_Preorder().
577 *
578 * \param eMatrix
579 * Pointer to the matrix.
580 * \param RHS
581 * RHS is the right hand side. This is what is being solved for.
582 * \param Solution
583 * Solution is the vector being multiplied by the matrix.
584 * \param iRHS
585 * iRHS is the imaginary portion of the right hand side. This is
586 * what is being solved for. This is only necessary if the matrix is
587 * complex and \a spSEPARATED_COMPLEX_VECTORS is true.
588 * \param iSolution
589 * iSolution is the imaginary portion of the vector being multiplied
590 * by the matrix. This is only necessary if the matrix is
591 * complex and \a spSEPARATED_COMPLEX_VECTORS is true.
592 */
593
594 void
spMultiply(spMatrix eMatrix,spREAL RHS[],spREAL Solution[],spREAL iRHS[],spREAL iSolution[])595 spMultiply(
596 spMatrix eMatrix,
597 spREAL RHS[],
598 spREAL Solution[]
599 #if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS
600 , spREAL iRHS[]
601 , spREAL iSolution[]
602 #endif
603 )
604 {
605 register ElementPtr pElement;
606 register RealVector Vector;
607 register RealNumber Sum;
608 register int I, *pExtOrder;
609 MatrixPtr Matrix = (MatrixPtr)eMatrix;
610 extern void ComplexMatrixMultiply();
611
612 /* Begin `spMultiply'. */
613 ASSERT_IS_SPARSE( Matrix );
614 ASSERT_IS_NOT_FACTORED( Matrix );
615 if (NOT Matrix->RowsLinked)
616 spcLinkRows(Matrix);
617 if (NOT Matrix->InternalVectorsAllocated)
618 spcCreateInternalVectors( Matrix );
619
620 #if spCOMPLEX
621 if (Matrix->Complex)
622 { ComplexMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS );
623 return;
624 }
625 #endif
626
627 #if REAL
628 #if NOT ARRAY_OFFSET
629 /* Correct array pointers for ARRAY_OFFSET. */
630 --RHS;
631 --Solution;
632 #endif
633
634 /* Initialize Intermediate vector with reordered Solution vector. */
635 Vector = Matrix->Intermediate;
636 pExtOrder = &Matrix->IntToExtColMap[Matrix->Size];
637 for (I = Matrix->Size; I > 0; I--)
638 Vector[I] = Solution[*(pExtOrder--)];
639
640 pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size];
641 for (I = Matrix->Size; I > 0; I--)
642 { pElement = Matrix->FirstInRow[I];
643 Sum = 0.0;
644
645 while (pElement != NULL)
646 { Sum += pElement->Real * Vector[pElement->Col];
647 pElement = pElement->NextInRow;
648 }
649 RHS[*pExtOrder--] = Sum;
650 }
651 return;
652 #endif /* REAL */
653 }
654 #endif /* MULTIPLICATION */
655
656
657
658
659
660
661
662 #if spCOMPLEX AND MULTIPLICATION
663 /*
664 * COMPLEX MATRIX MULTIPLICATION
665 *
666 * Multiplies matrix by solution vector to find source vector.
667 * Assumes matrix has not been factored. This routine can be used
668 * as a test to see if solutions are correct.
669 *
670 * >>> Arguments:
671 * Matrix <input> (char *)
672 * Pointer to the matrix.
673 * RHS <output> (RealVector)
674 * RHS is the right hand side. This is what is being solved for.
675 * This is only the real portion of the right-hand side if the matrix
676 * is complex and spSEPARATED_COMPLEX_VECTORS is set true.
677 * Solution <input> (RealVector)
678 * Solution is the vector being multiplied by the matrix. This is only
679 * the real portion if the matrix is complex and
680 * spSEPARATED_COMPLEX_VECTORS is set true.
681 * iRHS <output> (RealVector)
682 * iRHS is the imaginary portion of the right hand side. This is
683 * what is being solved for. This is only necessary if the matrix is
684 * complex and spSEPARATED_COMPLEX_VECTORS is true.
685 * iSolution <input> (RealVector)
686 * iSolution is the imaginary portion of the vector being multiplied
687 * by the matrix. This is only necessary if the matrix is
688 * complex and spSEPARATED_COMPLEX_VECTORS is true.
689 */
690
691 static void
ComplexMatrixMultiply(MatrixPtr Matrix,RealVector RHS,RealVector Solution,RealVector iRHS,RealVector iSolution)692 ComplexMatrixMultiply(
693 MatrixPtr Matrix,
694 RealVector RHS,
695 RealVector Solution
696 #if spSEPARATED_COMPLEX_VECTORS
697 , RealVector iRHS
698 , RealVector iSolution
699 #endif
700 )
701 {
702 register ElementPtr pElement;
703 register ComplexVector Vector;
704 ComplexNumber Sum;
705 register int I, *pExtOrder;
706
707 /* Begin `ComplexMatrixMultiply'. */
708
709 /* Correct array pointers for ARRAY_OFFSET. */
710 #if NOT ARRAY_OFFSET
711 #if spSEPARATED_COMPLEX_VECTORS
712 --RHS; --iRHS;
713 --Solution; --iSolution;
714 #else
715 RHS -= 2; Solution -= 2;
716 #endif
717 #endif
718
719 /* Initialize Intermediate vector with reordered Solution vector. */
720 Vector = (ComplexVector)Matrix->Intermediate;
721 pExtOrder = &Matrix->IntToExtColMap[Matrix->Size];
722
723 #if spSEPARATED_COMPLEX_VECTORS
724 for (I = Matrix->Size; I > 0; I--)
725 { Vector[I].Real = Solution[*pExtOrder];
726 Vector[I].Imag = iSolution[*(pExtOrder--)];
727 }
728 #else
729 for (I = Matrix->Size; I > 0; I--)
730 Vector[I] = ((ComplexVector)Solution)[*(pExtOrder--)];
731 #endif
732
733 pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size];
734 for (I = Matrix->Size; I > 0; I--)
735 { pElement = Matrix->FirstInRow[I];
736 Sum.Real = Sum.Imag = 0.0;
737
738 while (pElement != NULL)
739 { /* Cmplx expression : Sum += Element * Vector[Col] */
740 CMPLX_MULT_ADD_ASSIGN( Sum, *pElement, Vector[pElement->Col] );
741 pElement = pElement->NextInRow;
742 }
743
744 #if spSEPARATED_COMPLEX_VECTORS
745 RHS[*pExtOrder] = Sum.Real;
746 iRHS[*pExtOrder--] = Sum.Imag;
747 #else
748 ((ComplexVector)RHS)[*pExtOrder--] = Sum;
749 #endif
750 }
751 return;
752 }
753 #endif /* spCOMPLEX AND MULTIPLICATION */
754
755
756
757
758
759
760
761
762 #if MULTIPLICATION AND TRANSPOSE
763 /*!
764 * Multiplies transposed matrix by solution vector to find source vector.
765 * Assumes matrix has not been factored. This routine can be used
766 * as a test to see if solutions are correct. It should not be used
767 * before spMNA_Preorder().
768 *
769 * \param eMatrix
770 * Pointer to the matrix.
771 * \param RHS
772 * RHS is the right hand side. This is what is being solved for.
773 * \param Solution
774 * Solution is the vector being multiplied by the matrix.
775 * \param iRHS
776 * iRHS is the imaginary portion of the right hand side. This is
777 * what is being solved for. This is only necessary if the matrix is
778 * complex and \a spSEPARATED_COMPLEX_VECTORS is true.
779 * \param iSolution
780 * iSolution is the imaginary portion of the vector being multiplied
781 * by the matrix. This is only necessary if the matrix is
782 * complex and \a spSEPARATED_COMPLEX_VECTORS is true.
783 */
784
785 void
spMultTransposed(spMatrix eMatrix,spREAL RHS[],spREAL Solution[],spREAL iRHS[],spREAL iSolution[])786 spMultTransposed(
787 spMatrix eMatrix,
788 spREAL RHS[],
789 spREAL Solution[]
790 #if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS
791 , spREAL iRHS[]
792 , spREAL iSolution[]
793 #endif
794 )
795 {
796 register ElementPtr pElement;
797 register RealVector Vector;
798 register RealNumber Sum;
799 register int I, *pExtOrder;
800 MatrixPtr Matrix = (MatrixPtr)eMatrix;
801 extern void ComplexTransposedMatrixMultiply();
802
803 /* Begin `spMultTransposed'. */
804 ASSERT_IS_SPARSE( Matrix );
805 ASSERT_IS_NOT_FACTORED( Matrix );
806 if (NOT Matrix->InternalVectorsAllocated)
807 spcCreateInternalVectors( Matrix );
808
809 #if spCOMPLEX
810 if (Matrix->Complex)
811 { ComplexTransposedMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS );
812 return;
813 }
814 #endif
815
816 #if REAL
817 #if NOT ARRAY_OFFSET
818 /* Correct array pointers for ARRAY_OFFSET. */
819 --RHS;
820 --Solution;
821 #endif
822
823 /* Initialize Intermediate vector with reordered Solution vector. */
824 Vector = Matrix->Intermediate;
825 pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size];
826 for (I = Matrix->Size; I > 0; I--)
827 Vector[I] = Solution[*(pExtOrder--)];
828
829 pExtOrder = &Matrix->IntToExtColMap[Matrix->Size];
830 for (I = Matrix->Size; I > 0; I--)
831 { pElement = Matrix->FirstInCol[I];
832 Sum = 0.0;
833
834 while (pElement != NULL)
835 { Sum += pElement->Real * Vector[pElement->Row];
836 pElement = pElement->NextInCol;
837 }
838 RHS[*pExtOrder--] = Sum;
839 }
840 return;
841 #endif /* REAL */
842 }
843 #endif /* MULTIPLICATION AND TRANSPOSE */
844
845
846
847
848
849
850
851 #if spCOMPLEX AND MULTIPLICATION AND TRANSPOSE
852 /*
853 * COMPLEX TRANSPOSED MATRIX MULTIPLICATION
854 *
855 * Multiplies transposed matrix by solution vector to find source vector.
856 * Assumes matrix has not been factored. This routine can be used
857 * as a test to see if solutions are correct.
858 *
859 * >>> Arguments:
860 * Matrix <input> (char *)
861 * Pointer to the matrix.
862 * RHS <output> (RealVector)
863 * RHS is the right hand side. This is what is being solved for.
864 * This is only the real portion of the right-hand side if the matrix
865 * is complex and spSEPARATED_COMPLEX_VECTORS is set true.
866 * Solution <input> (RealVector)
867 * Solution is the vector being multiplied by the matrix. This is only
868 * the real portion if the matrix is complex and
869 * spSEPARATED_COMPLEX_VECTORS is set true.
870 * iRHS <output> (RealVector)
871 * iRHS is the imaginary portion of the right hand side. This is
872 * what is being solved for. This is only necessary if the matrix is
873 * complex and spSEPARATED_COMPLEX_VECTORS is true.
874 * iSolution <input> (RealVector)
875 * iSolution is the imaginary portion of the vector being multiplied
876 * by the matrix. This is only necessary if the matrix is
877 * complex and spSEPARATED_COMPLEX_VECTORS is true.
878 *
879 * >>> Obscure Macros
880 * IMAG_VECTORS
881 * Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
882 * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
883 * without a trace.
884 */
885
886 static void
ComplexTransposedMatrixMultiply(MatrixPtr Matrix,RealVector RHS,RealVector Solution,RealVector iRHS,RealVector iSolution)887 ComplexTransposedMatrixMultiply(
888 MatrixPtr Matrix,
889 RealVector RHS,
890 RealVector Solution
891 #if spSEPARATED_COMPLEX_VECTORS
892 , RealVector iRHS
893 , RealVector iSolution
894 #endif
895 )
896 {
897 register ElementPtr pElement;
898 register ComplexVector Vector;
899 ComplexNumber Sum;
900 register int I, *pExtOrder;
901
902 /* Begin `ComplexTransposedMatrixMultiply'. */
903
904 /* Correct array pointers for ARRAY_OFFSET. */
905 #if NOT ARRAY_OFFSET
906 #if spSEPARATED_COMPLEX_VECTORS
907 --RHS; --iRHS;
908 --Solution; --iSolution;
909 #else
910 RHS -= 2; Solution -= 2;
911 #endif
912 #endif
913
914 /* Initialize Intermediate vector with reordered Solution vector. */
915 Vector = (ComplexVector)Matrix->Intermediate;
916 pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size];
917
918 #if spSEPARATED_COMPLEX_VECTORS
919 for (I = Matrix->Size; I > 0; I--)
920 { Vector[I].Real = Solution[*pExtOrder];
921 Vector[I].Imag = iSolution[*(pExtOrder--)];
922 }
923 #else
924 for (I = Matrix->Size; I > 0; I--)
925 Vector[I] = ((ComplexVector)Solution)[*(pExtOrder--)];
926 #endif
927
928 pExtOrder = &Matrix->IntToExtColMap[Matrix->Size];
929 for (I = Matrix->Size; I > 0; I--)
930 { pElement = Matrix->FirstInCol[I];
931 Sum.Real = Sum.Imag = 0.0;
932
933 while (pElement != NULL)
934 { /* Cmplx expression : Sum += Element * Vector[Row] */
935 CMPLX_MULT_ADD_ASSIGN( Sum, *pElement, Vector[pElement->Row] );
936 pElement = pElement->NextInCol;
937 }
938
939 #if spSEPARATED_COMPLEX_VECTORS
940 RHS[*pExtOrder] = Sum.Real;
941 iRHS[*pExtOrder--] = Sum.Imag;
942 #else
943 ((ComplexVector)RHS)[*pExtOrder--] = Sum;
944 #endif
945 }
946 return;
947 }
948 #endif /* spCOMPLEX AND MULTIPLICATION AND TRANSPOSE */
949
950
951
952
953
954
955
956
957 #if DETERMINANT
958 /*!
959 * This routine in capable of calculating the determinant of the
960 * matrix once the LU factorization has been performed. Hence, only
961 * use this routine after spFactor() and before spClear().
962 * The determinant equals the product of all the diagonal elements of
963 * the lower triangular matrix L, except that this product may need
964 * negating. Whether the product or the negative product equals the
965 * determinant is determined by the number of row and column
966 * interchanges performed. Note that the determinants of matrices can
967 * be very large or very small. On large matrices, the determinant
968 * can be far larger or smaller than can be represented by a floating
969 * point number. For this reason the determinant is scaled to a
970 * reasonable value and the logarithm of the scale factor is returned.
971 *
972 * \param eMatrix
973 * A pointer to the matrix for which the determinant is desired.
974 * \param pExponent
975 * The logarithm base 10 of the scale factor for the determinant. To find
976 * the actual determinant, Exponent should be added to the exponent of
977 * Determinant.
978 * \param pDeterminant
979 * The real portion of the determinant. This number is scaled to be
980 * greater than or equal to 1.0 and less than 10.0.
981 * \param piDeterminant
982 * The imaginary portion of the determinant. When the matrix is real
983 * this pointer need not be supplied, nothing will be returned. This
984 * number is scaled to be greater than or equal to 1.0 and less than 10.0.
985 */
986 /* >>> Local variables:
987 * Norm (RealNumber)
988 * L-infinity norm of a complex number.
989 * Size (int)
990 * Local storage for Matrix->Size. Placed in a register for speed.
991 * Temp (RealNumber)
992 * Temporary storage for real portion of determinant.
993 */
994
995 void
spDeterminant(spMatrix eMatrix,int * pExponent,spREAL * pDeterminant,spREAL * piDeterminant)996 spDeterminant(
997 spMatrix eMatrix,
998 int *pExponent,
999 spREAL *pDeterminant
1000 #if spCOMPLEX
1001 , spREAL *piDeterminant
1002 #endif
1003 )
1004 {
1005 register MatrixPtr Matrix = (MatrixPtr)eMatrix;
1006 register int I, Size;
1007 #if spCOMPLEX
1008 RealNumber Norm, nr, ni;
1009 ComplexNumber Pivot, cDeterminant;
1010 #endif
1011 #define NORM(a) (nr = ABS((a).Real), ni = ABS((a).Imag), MAX (nr,ni))
1012
1013 /* Begin `spDeterminant'. */
1014 ASSERT_IS_SPARSE( Matrix );
1015 ASSERT_NO_ERRORS( Matrix );
1016 ASSERT_IS_FACTORED( Matrix );
1017 *pExponent = 0;
1018
1019 if (Matrix->Error == spSINGULAR)
1020 { *pDeterminant = 0.0;
1021 #if spCOMPLEX
1022 if (Matrix->Complex) *piDeterminant = 0.0;
1023 #endif
1024 return;
1025 }
1026
1027 Size = Matrix->Size;
1028 I = 0;
1029
1030 #if spCOMPLEX
1031 if (Matrix->Complex) /* Complex Case. */
1032 { cDeterminant.Real = 1.0;
1033 cDeterminant.Imag = 0.0;
1034
1035 while (++I <= Size)
1036 { CMPLX_RECIPROCAL( Pivot, *Matrix->Diag[I] );
1037 CMPLX_MULT_ASSIGN( cDeterminant, Pivot );
1038
1039 /* Scale Determinant. */
1040 Norm = NORM( cDeterminant );
1041 if (Norm != 0.0)
1042 { while (Norm >= 1.0e12)
1043 { cDeterminant.Real *= 1.0e-12;
1044 cDeterminant.Imag *= 1.0e-12;
1045 *pExponent += 12;
1046 Norm = NORM( cDeterminant );
1047 }
1048 while (Norm < 1.0e-12)
1049 { cDeterminant.Real *= 1.0e12;
1050 cDeterminant.Imag *= 1.0e12;
1051 *pExponent -= 12;
1052 Norm = NORM( cDeterminant );
1053 }
1054 }
1055 }
1056
1057 /* Scale Determinant again, this time to be between 1.0 <= x < 10.0. */
1058 Norm = NORM( cDeterminant );
1059 if (Norm != 0.0)
1060 { while (Norm >= 10.0)
1061 { cDeterminant.Real *= 0.1;
1062 cDeterminant.Imag *= 0.1;
1063 (*pExponent)++;
1064 Norm = NORM( cDeterminant );
1065 }
1066 while (Norm < 1.0)
1067 { cDeterminant.Real *= 10.0;
1068 cDeterminant.Imag *= 10.0;
1069 (*pExponent)--;
1070 Norm = NORM( cDeterminant );
1071 }
1072 }
1073 if (Matrix->NumberOfInterchangesIsOdd)
1074 CMPLX_NEGATE( cDeterminant );
1075
1076 *pDeterminant = cDeterminant.Real;
1077 *piDeterminant = cDeterminant.Imag;
1078 }
1079 #endif /* spCOMPLEX */
1080 #if REAL AND spCOMPLEX
1081 else
1082 #endif
1083 #if REAL
1084 { /* Real Case. */
1085 *pDeterminant = 1.0;
1086
1087 while (++I <= Size)
1088 { *pDeterminant /= Matrix->Diag[I]->Real;
1089
1090 /* Scale Determinant. */
1091 if (*pDeterminant != 0.0)
1092 { while (ABS(*pDeterminant) >= 1.0e12)
1093 { *pDeterminant *= 1.0e-12;
1094 *pExponent += 12;
1095 }
1096 while (ABS(*pDeterminant) < 1.0e-12)
1097 { *pDeterminant *= 1.0e12;
1098 *pExponent -= 12;
1099 }
1100 }
1101 }
1102
1103 /* Scale Determinant again, this time to be between 1.0 <= x < 10.0. */
1104 if (*pDeterminant != 0.0)
1105 { while (ABS(*pDeterminant) >= 10.0)
1106 { *pDeterminant *= 0.1;
1107 (*pExponent)++;
1108 }
1109 while (ABS(*pDeterminant) < 1.0)
1110 { *pDeterminant *= 10.0;
1111 (*pExponent)--;
1112 }
1113 }
1114 if (Matrix->NumberOfInterchangesIsOdd)
1115 *pDeterminant = -*pDeterminant;
1116 }
1117 #endif /* REAL */
1118 }
1119 #endif /* DETERMINANT */
1120
1121
1122
1123
1124
1125
1126
1127
1128 #if STRIP
1129
1130 /*!
1131 * Strips the matrix of all fill-ins.
1132 *
1133 * \param eMatrix
1134 * Pointer to the matrix to be stripped.
1135 */
1136 /* >>> Local variables:
1137 * pElement (ElementPtr)
1138 * Pointer that is used to step through the matrix.
1139 * ppElement (ElementPtr *)
1140 * Pointer to the location of an ElementPtr. This location will be
1141 * updated if a fill-in is stripped from the matrix.
1142 * pFillin (ElementPtr)
1143 * Pointer used to step through the lists of fill-ins while marking them.
1144 * pLastFillin (ElementPtr)
1145 * A pointer to the last fill-in in the list. Used to terminate a loop.
1146 * pListNode (struct FillinListNodeStruct *)
1147 * A pointer to a node in the FillinList linked-list.
1148 */
1149
1150 void
spStripFills(spMatrix eMatrix)1151 spStripFills( spMatrix eMatrix )
1152 {
1153 MatrixPtr Matrix = (MatrixPtr)eMatrix;
1154 struct FillinListNodeStruct *pListNode;
1155
1156 /* Begin `spStripFills'. */
1157 ASSERT_IS_SPARSE( Matrix );
1158 if (Matrix->Fillins == 0) return;
1159 Matrix->NeedsOrdering = YES;
1160 Matrix->Elements -= Matrix->Fillins;
1161 Matrix->Fillins = 0;
1162
1163 /* Mark the fill-ins. */
1164 { register ElementPtr pFillin, pLastFillin;
1165
1166 pListNode = Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
1167 Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList;
1168 Matrix->NextAvailFillin = pListNode->pFillinList;
1169
1170 while (pListNode != NULL)
1171 { pFillin = pListNode->pFillinList;
1172 pLastFillin = &(pFillin[ pListNode->NumberOfFillinsInList - 1 ]);
1173 while (pFillin <= pLastFillin)
1174 (pFillin++)->Row = 0;
1175 pListNode = pListNode->Next;
1176 }
1177 }
1178
1179 /* Unlink fill-ins by searching for elements marked with Row = 0. */
1180 { register ElementPtr pElement, *ppElement;
1181 register int I, Size = Matrix->Size;
1182
1183 /* Unlink fill-ins in all columns. */
1184 for (I = 1; I <= Size; I++)
1185 { ppElement = &(Matrix->FirstInCol[I]);
1186 while ((pElement = *ppElement) != NULL)
1187 { if (pElement->Row == 0)
1188 { *ppElement = pElement->NextInCol; /* Unlink fill-in. */
1189 if (Matrix->Diag[pElement->Col] == pElement)
1190 Matrix->Diag[pElement->Col] = NULL;
1191 }
1192 else
1193 ppElement = &pElement->NextInCol; /* Skip element. */
1194 }
1195 }
1196
1197 /* Unlink fill-ins in all rows. */
1198 for (I = 1; I <= Size; I++)
1199 { ppElement = &(Matrix->FirstInRow[I]);
1200 while ((pElement = *ppElement) != NULL)
1201 { if (pElement->Row == 0)
1202 *ppElement = pElement->NextInRow; /* Unlink fill-in. */
1203 else
1204 ppElement = &pElement->NextInRow; /* Skip element. */
1205 }
1206 }
1207 }
1208 return;
1209 }
1210 #endif
1211
1212
1213
1214
1215
1216
1217
1218 #if TRANSLATE AND DELETE
1219 /*!
1220 * Deletes a row and a column from a matrix.
1221 *
1222 * Sparse will abort if an attempt is made to delete a row or column that
1223 * doesn't exist.
1224 *
1225 * \param eMatrix
1226 * Pointer to the matrix in which the row and column are to be deleted.
1227 * \param Row
1228 * Row to be deleted.
1229 * \param Col
1230 * Column to be deleted.
1231 */
1232 /* >>> Local variables:
1233 * ExtCol (int)
1234 * The external column that is being deleted.
1235 * ExtRow (int)
1236 * The external row that is being deleted.
1237 * pElement (ElementPtr)
1238 * Pointer to an element in the matrix. Used when scanning rows and
1239 * columns in order to eliminate elements from the last row or column.
1240 * ppElement (ElementPtr *)
1241 * Pointer to the location of an ElementPtr. This location will be
1242 * filled with a NULL pointer if it is the new last element in its row
1243 * or column.
1244 * pElement (ElementPtr)
1245 * Pointer to an element in the last row or column of the matrix.
1246 * Size (int)
1247 * The local version Matrix->Size, the size of the matrix.
1248 */
1249
1250 void
spDeleteRowAndCol(spMatrix eMatrix,int Row,int Col)1251 spDeleteRowAndCol(
1252 spMatrix eMatrix,
1253 int Row,
1254 int Col
1255 )
1256 {
1257 MatrixPtr Matrix = (MatrixPtr)eMatrix;
1258 register ElementPtr pElement, *ppElement, pLastElement;
1259 int Size, ExtRow, ExtCol;
1260
1261 /* Begin `spDeleteRowAndCol'. */
1262 ASSERT_IS_SPARSE( Matrix );
1263 vASSERT( (Row > 0) AND (Col > 0), "Nonpositive row or column number" );
1264 vASSERT( (Row <= Matrix->ExtSize) AND (Col <= Matrix->ExtSize),
1265 "Row or column number too large" );
1266
1267 Size = Matrix->Size;
1268 ExtRow = Row;
1269 ExtCol = Col;
1270 if (NOT Matrix->RowsLinked) spcLinkRows( Matrix );
1271
1272 Row = Matrix->ExtToIntRowMap[Row];
1273 Col = Matrix->ExtToIntColMap[Col];
1274 ASSERT( Row > 0 AND Col > 0 );
1275
1276 /* Move Row so that it is the last row in the matrix. */
1277 if (Row != Size) spcRowExchange( Matrix, Row, Size );
1278
1279 /* Move Col so that it is the last column in the matrix. */
1280 if (Col != Size) spcColExchange( Matrix, Col, Size );
1281
1282 /* Correct Diag pointers. */
1283 if (Row == Col)
1284 SWAP( ElementPtr, Matrix->Diag[Row], Matrix->Diag[Size] )
1285 else
1286 { Matrix->Diag[Row] = spcFindDiag( Matrix, Row );
1287 Matrix->Diag[Col] = spcFindDiag( Matrix, Col );
1288 }
1289
1290 /*
1291 * Delete last row and column of the matrix.
1292 */
1293 /* Break the column links to every element in the last row. */
1294 pLastElement = Matrix->FirstInRow[ Size ];
1295 while (pLastElement != NULL)
1296 { ppElement = &(Matrix->FirstInCol[ pLastElement->Col ]);
1297 while ((pElement = *ppElement) != NULL)
1298 { if (pElement == pLastElement)
1299 *ppElement = NULL; /* Unlink last element in column. */
1300 else
1301 ppElement = &pElement->NextInCol; /* Skip element. */
1302 }
1303 pLastElement = pLastElement->NextInRow;
1304 }
1305
1306 /* Break the row links to every element in the last column. */
1307 pLastElement = Matrix->FirstInCol[ Size ];
1308 while (pLastElement != NULL)
1309 { ppElement = &(Matrix->FirstInRow[ pLastElement->Row ]);
1310 while ((pElement = *ppElement) != NULL)
1311 { if (pElement == pLastElement)
1312 *ppElement = NULL; /* Unlink last element in row. */
1313 else
1314 ppElement = &pElement->NextInRow; /* Skip element. */
1315 }
1316 pLastElement = pLastElement->NextInCol;
1317 }
1318
1319 /* Clean up some details. */
1320 Matrix->Size = Size - 1;
1321 Matrix->Diag[Size] = NULL;
1322 Matrix->FirstInRow[Size] = NULL;
1323 Matrix->FirstInCol[Size] = NULL;
1324 Matrix->CurrentSize--;
1325 Matrix->ExtToIntRowMap[ExtRow] = -1;
1326 Matrix->ExtToIntColMap[ExtCol] = -1;
1327 Matrix->NeedsOrdering = YES;
1328
1329 return;
1330 }
1331 #endif
1332
1333
1334
1335
1336
1337
1338
1339
1340 #if PSEUDOCONDITION
1341 /*!
1342 * Computes the magnitude of the ratio of the largest to the smallest
1343 * pivots. This quantity is an indicator of ill-conditioning in the
1344 * matrix. If this ratio is large, and if the matrix is scaled such
1345 * that uncertainties in the RHS and the matrix entries are
1346 * equilibrated, then the matrix is ill-conditioned. However, a small
1347 * ratio does not necessarily imply that the matrix is
1348 * well-conditioned. This routine must only be used after a matrix has
1349 * been factored by spOrderAndFactor() or spFactor() and before it is
1350 * cleared by spClear() or spInitialize(). The pseudocondition is
1351 * faster to compute than the condition number calculated by
1352 * spCondition(), but is not as informative.
1353 *
1354 * \return
1355 * The magnitude of the ratio of the largest to smallest pivot used during
1356 * previous factorization. If the matrix was singular, zero is returned.
1357 *
1358 * \param eMatrix
1359 * Pointer to the matrix.
1360 */
1361
1362 spREAL
spPseudoCondition(spMatrix eMatrix)1363 spPseudoCondition( spMatrix eMatrix )
1364 {
1365 MatrixPtr Matrix = (MatrixPtr)eMatrix;
1366 register int I;
1367 register ArrayOfElementPtrs Diag;
1368 RealNumber MaxPivot, MinPivot, Mag;
1369
1370 /* Begin `spPseudoCondition'. */
1371 ASSERT_IS_SPARSE( Matrix );
1372 ASSERT_NO_ERRORS( Matrix );
1373 ASSERT_IS_FACTORED( Matrix );
1374 if (Matrix->Error == spSINGULAR OR Matrix->Error == spZERO_DIAG)
1375 return 0.0;
1376
1377 Diag = Matrix->Diag;
1378 MaxPivot = MinPivot = ELEMENT_MAG( Diag[1] );
1379 for (I = 2; I <= Matrix->Size; I++)
1380 { Mag = ELEMENT_MAG( Diag[I] );
1381 if (Mag > MaxPivot)
1382 MaxPivot = Mag;
1383 else if (Mag < MinPivot)
1384 MinPivot = Mag;
1385 }
1386 ASSERT( MaxPivot > 0.0 );
1387 return MaxPivot / MinPivot;
1388 }
1389 #endif
1390
1391
1392
1393
1394
1395
1396
1397
1398 #if CONDITION
1399 /*!
1400 * Computes an estimate of the condition number using a variation on
1401 * the LINPACK condition number estimation algorithm. This quantity is
1402 * an indicator of ill-conditioning in the matrix. To avoid problems
1403 * with overflow, the reciprocal of the condition number is returned.
1404 * If this number is small, and if the matrix is scaled such that
1405 * uncertainties in the RHS and the matrix entries are equilibrated,
1406 * then the matrix is ill-conditioned. If the this number is near
1407 * one, the matrix is well conditioned. This routine must only be
1408 * used after a matrix has been factored by spOrderAndFactor() or
1409 * spFactor() and before it is cleared by spClear() or spInitialize().
1410 *
1411 * Unlike the LINPACK condition number estimator, this routines
1412 * returns the L infinity condition number. This is an artifact of
1413 * Sparse placing ones on the diagonal of the upper triangular matrix
1414 * rather than the lower. This difference should be of no importance.
1415 *
1416 * \b References:
1417 *
1418 * A.K. Cline, C.B. Moler, G.W. Stewart, J.H. Wilkinson. An estimate
1419 * for the condition number of a matrix. SIAM Journal on Numerical
1420 * Analysis. Vol. 16, No. 2, pages 368-375, April 1979.
1421 *
1422 * J.J. Dongarra, C.B. Moler, J.R. Bunch, G.W. Stewart. LINPACK
1423 * User's Guide. SIAM, 1979.
1424 *
1425 * Roger G. Grimes, John G. Lewis. Condition number estimation for
1426 * sparse matrices. SIAM Journal on Scientific and Statistical
1427 * Computing. Vol. 2, No. 4, pages 384-388, December 1981.
1428 *
1429 * Dianne Prost O'Leary. Estimating matrix condition numbers. SIAM
1430 * Journal on Scientific and Statistical Computing. Vol. 1, No. 2,
1431 * pages 205-209, June 1980.
1432 *
1433 * \return
1434 * The reciprocal of the condition number. If the matrix was singular,
1435 * zero is returned.
1436 *
1437 * \param eMatrix
1438 * Pointer to the matrix.
1439 * \param NormOfMatrix
1440 * The L-infinity norm of the unfactored matrix as computed by
1441 * spNorm().
1442 * \param pError
1443 * Used to return error code. Possible errors include \a spSINGULAR
1444 * or \a spNO_MEMORY.
1445 */
1446
1447 spREAL
spCondition(spMatrix eMatrix,spREAL NormOfMatrix,int * pError)1448 spCondition(
1449 spMatrix eMatrix,
1450 spREAL NormOfMatrix,
1451 int *pError
1452 )
1453 {
1454 MatrixPtr Matrix = (MatrixPtr)eMatrix;
1455 register ElementPtr pElement;
1456 register RealVector T, Tm;
1457 register int I, K, Row;
1458 ElementPtr pPivot;
1459 int Size;
1460 RealNumber E, Em, Wp, Wm, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor;
1461 RealNumber Linpack, OLeary, InvNormOfInverse, ComplexCondition();
1462 #define SLACK 1e4
1463
1464 /* Begin `spCondition'. */
1465 ASSERT_IS_SPARSE( Matrix );
1466 ASSERT_NO_ERRORS( Matrix );
1467 ASSERT_IS_FACTORED( Matrix );
1468 *pError = Matrix->Error;
1469 if (Matrix->Error >= spFATAL) return 0.0;
1470 if (NormOfMatrix == 0.0)
1471 { *pError = spSINGULAR;
1472 return 0.0;
1473 }
1474
1475 #if spCOMPLEX
1476 if (Matrix->Complex)
1477 return ComplexCondition( Matrix, NormOfMatrix, pError );
1478 #endif
1479
1480 #if REAL
1481 Size = Matrix->Size;
1482 T = Matrix->Intermediate;
1483 #if spCOMPLEX
1484 Tm = Matrix->Intermediate + Size;
1485 #else
1486 Tm = ALLOC( RealNumber, Size+1 );
1487 if (Tm == NULL)
1488 { *pError = spNO_MEMORY;
1489 return 0.0;
1490 }
1491 #endif
1492 for (I = Size; I > 0; I--) T[I] = 0.0;
1493
1494 /*
1495 * Part 1. Ay = e.
1496 * Solve Ay = LUy = e where e consists of +1 and -1 terms with the sign
1497 * chosen to maximize the size of w in Lw = e. Since the terms in w can
1498 * get very large, scaling is used to avoid overflow.
1499 */
1500
1501 /* Forward elimination. Solves Lw = e while choosing e. */
1502 E = 1.0;
1503 for (I = 1; I <= Size; I++)
1504 { pPivot = Matrix->Diag[I];
1505 if (T[I] < 0.0) Em = -E; else Em = E;
1506 Wm = (Em + T[I]) * pPivot->Real;
1507 if (ABS(Wm) > SLACK)
1508 { ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(Wm) );
1509 for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
1510 E *= ScaleFactor;
1511 Em *= ScaleFactor;
1512 Wm = (Em + T[I]) * pPivot->Real;
1513 }
1514 Wp = (T[I] - Em) * pPivot->Real;
1515 ASp = ABS(T[I] - Em);
1516 ASm = ABS(Em + T[I]);
1517
1518 /* Update T for both values of W, minus value is placed in Tm. */
1519 pElement = pPivot->NextInCol;
1520 while (pElement != NULL)
1521 { Row = pElement->Row;
1522 Tm[Row] = T[Row] - (Wm * pElement->Real);
1523 T[Row] -= (Wp * pElement->Real);
1524 ASp += ABS(T[Row]);
1525 ASm += ABS(Tm[Row]);
1526 pElement = pElement->NextInCol;
1527 }
1528
1529 /* If minus value causes more growth, overwrite T with its values. */
1530 if (ASm > ASp)
1531 { T[I] = Wm;
1532 pElement = pPivot->NextInCol;
1533 while (pElement != NULL)
1534 { T[pElement->Row] = Tm[pElement->Row];
1535 pElement = pElement->NextInCol;
1536 }
1537 }
1538 else T[I] = Wp;
1539 }
1540
1541 /* Compute 1-norm of T, which now contains w, and scale ||T|| to 1/SLACK. */
1542 for (ASw = 0.0, I = Size; I > 0; I--) ASw += ABS(T[I]);
1543 ScaleFactor = 1.0 / (SLACK * ASw);
1544 if (ScaleFactor < 0.5)
1545 { for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
1546 E *= ScaleFactor;
1547 }
1548
1549 /* Backward Substitution. Solves Uy = w.*/
1550 for (I = Size; I >= 1; I--)
1551 { pElement = Matrix->Diag[I]->NextInRow;
1552 while (pElement != NULL)
1553 { T[I] -= pElement->Real * T[pElement->Col];
1554 pElement = pElement->NextInRow;
1555 }
1556 if (ABS(T[I]) > SLACK)
1557 { ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) );
1558 for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
1559 E *= ScaleFactor;
1560 }
1561 }
1562
1563 /* Compute 1-norm of T, which now contains y, and scale ||T|| to 1/SLACK. */
1564 for (ASy = 0.0, I = Size; I > 0; I--) ASy += ABS(T[I]);
1565 ScaleFactor = 1.0 / (SLACK * ASy);
1566 if (ScaleFactor < 0.5)
1567 { for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
1568 ASy = 1.0 / SLACK;
1569 E *= ScaleFactor;
1570 }
1571
1572 /* Compute infinity-norm of T for O'Leary's estimate. */
1573 for (MaxY = 0.0, I = Size; I > 0; I--)
1574 if (MaxY < ABS(T[I])) MaxY = ABS(T[I]);
1575
1576 /*
1577 * Part 2. A* z = y where the * represents the transpose.
1578 * Recall that A = LU implies A* = U* L*.
1579 */
1580
1581 /* Forward elimination, U* v = y. */
1582 for (I = 1; I <= Size; I++)
1583 { pElement = Matrix->Diag[I]->NextInRow;
1584 while (pElement != NULL)
1585 { T[pElement->Col] -= T[I] * pElement->Real;
1586 pElement = pElement->NextInRow;
1587 }
1588 if (ABS(T[I]) > SLACK)
1589 { ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) );
1590 for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
1591 ASy *= ScaleFactor;
1592 }
1593 }
1594
1595 /* Compute 1-norm of T, which now contains v, and scale ||T|| to 1/SLACK. */
1596 for (ASv = 0.0, I = Size; I > 0; I--) ASv += ABS(T[I]);
1597 ScaleFactor = 1.0 / (SLACK * ASv);
1598 if (ScaleFactor < 0.5)
1599 { for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
1600 ASy *= ScaleFactor;
1601 }
1602
1603 /* Backward Substitution, L* z = v. */
1604 for (I = Size; I >= 1; I--)
1605 { pPivot = Matrix->Diag[I];
1606 pElement = pPivot->NextInCol;
1607 while (pElement != NULL)
1608 { T[I] -= pElement->Real * T[pElement->Row];
1609 pElement = pElement->NextInCol;
1610 }
1611 T[I] *= pPivot->Real;
1612 if (ABS(T[I]) > SLACK)
1613 { ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) );
1614 for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
1615 ASy *= ScaleFactor;
1616 }
1617 }
1618
1619 /* Compute 1-norm of T, which now contains z. */
1620 for (ASz = 0.0, I = Size; I > 0; I--) ASz += ABS(T[I]);
1621
1622 #if NOT spCOMPLEX
1623 FREE( Tm );
1624 #endif
1625
1626 Linpack = ASy / ASz;
1627 OLeary = E / MaxY;
1628 InvNormOfInverse = MIN( Linpack, OLeary );
1629 return InvNormOfInverse / NormOfMatrix;
1630 #endif /* REAL */
1631 }
1632
1633
1634
1635
1636
1637 #if spCOMPLEX
1638 /*
1639 * ESTIMATE CONDITION NUMBER
1640 *
1641 * Complex version of spCondition().
1642 *
1643 * >>> Returns:
1644 * The reciprocal of the condition number.
1645 *
1646 * >>> Arguments:
1647 * Matrix <input> (MatrixPtr)
1648 * Pointer to the matrix.
1649 * NormOfMatrix <input> (RealNumber)
1650 * The L-infinity norm of the unfactored matrix as computed by
1651 * spNorm().
1652 * pError <output> (int *)
1653 * Used to return error code.
1654 *
1655 * >>> Possible errors:
1656 * spNO_MEMORY
1657 */
1658
1659 static RealNumber
ComplexCondition(MatrixPtr Matrix,RealNumber NormOfMatrix,int * pError)1660 ComplexCondition(
1661 MatrixPtr Matrix,
1662 RealNumber NormOfMatrix,
1663 int *pError
1664 )
1665 {
1666 register ElementPtr pElement;
1667 register ComplexVector T, Tm;
1668 register int I, K, Row;
1669 ElementPtr pPivot;
1670 int Size;
1671 RealNumber E, Em, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor;
1672 RealNumber Linpack, OLeary, InvNormOfInverse;
1673 ComplexNumber Wp, Wm;
1674
1675 /* Begin `ComplexCondition'. */
1676
1677 Size = Matrix->Size;
1678 T = (ComplexVector)Matrix->Intermediate;
1679 Tm = ALLOC( ComplexNumber, Size+1 );
1680 if (Tm == NULL)
1681 { *pError = spNO_MEMORY;
1682 return 0.0;
1683 }
1684 for (I = Size; I > 0; I--) T[I].Real = T[I].Imag = 0.0;
1685
1686 /*
1687 * Part 1. Ay = e.
1688 * Solve Ay = LUy = e where e consists of +1 and -1 terms with the sign
1689 * chosen to maximize the size of w in Lw = e. Since the terms in w can
1690 * get very large, scaling is used to avoid overflow.
1691 */
1692
1693 /* Forward elimination. Solves Lw = e while choosing e. */
1694 E = 1.0;
1695 for (I = 1; I <= Size; I++)
1696 { pPivot = Matrix->Diag[I];
1697 if (T[I].Real < 0.0) Em = -E; else Em = E;
1698 Wm = T[I];
1699 Wm.Real += Em;
1700 ASm = CMPLX_1_NORM( Wm );
1701 CMPLX_MULT_ASSIGN( Wm, *pPivot );
1702 if (CMPLX_1_NORM(Wm) > SLACK)
1703 { ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(Wm) );
1704 for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
1705 E *= ScaleFactor;
1706 Em *= ScaleFactor;
1707 ASm *= ScaleFactor;
1708 SCLR_MULT_ASSIGN( Wm, ScaleFactor );
1709 }
1710 Wp = T[I];
1711 Wp.Real -= Em;
1712 ASp = CMPLX_1_NORM( Wp );
1713 CMPLX_MULT_ASSIGN( Wp, *pPivot );
1714
1715 /* Update T for both values of W, minus value is placed in Tm. */
1716 pElement = pPivot->NextInCol;
1717 while (pElement != NULL)
1718 { Row = pElement->Row;
1719 /* Cmplx expr: Tm[Row] = T[Row] - (Wp * *pElement). */
1720 CMPLX_MULT_SUBT( Tm[Row], Wm, *pElement, T[Row] );
1721 /* Cmplx expr: T[Row] -= Wp * *pElement. */
1722 CMPLX_MULT_SUBT_ASSIGN( T[Row], Wm, *pElement );
1723 ASp += CMPLX_1_NORM(T[Row]);
1724 ASm += CMPLX_1_NORM(Tm[Row]);
1725 pElement = pElement->NextInCol;
1726 }
1727
1728 /* If minus value causes more growth, overwrite T with its values. */
1729 if (ASm > ASp)
1730 { T[I] = Wm;
1731 pElement = pPivot->NextInCol;
1732 while (pElement != NULL)
1733 { T[pElement->Row] = Tm[pElement->Row];
1734 pElement = pElement->NextInCol;
1735 }
1736 }
1737 else T[I] = Wp;
1738 }
1739
1740 /* Compute 1-norm of T, which now contains w, and scale ||T|| to 1/SLACK. */
1741 for (ASw = 0.0, I = Size; I > 0; I--) ASw += CMPLX_1_NORM(T[I]);
1742 ScaleFactor = 1.0 / (SLACK * ASw);
1743 if (ScaleFactor < 0.5)
1744 { for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
1745 E *= ScaleFactor;
1746 }
1747
1748 /* Backward Substitution. Solves Uy = w.*/
1749 for (I = Size; I >= 1; I--)
1750 { pElement = Matrix->Diag[I]->NextInRow;
1751 while (pElement != NULL)
1752 { /* Cmplx expr: T[I] -= T[pElement->Col] * *pElement. */
1753 CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Col], *pElement );
1754 pElement = pElement->NextInRow;
1755 }
1756 if (CMPLX_1_NORM(T[I]) > SLACK)
1757 { ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
1758 for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
1759 E *= ScaleFactor;
1760 }
1761 }
1762
1763 /* Compute 1-norm of T, which now contains y, and scale ||T|| to 1/SLACK. */
1764 for (ASy = 0.0, I = Size; I > 0; I--) ASy += CMPLX_1_NORM(T[I]);
1765 ScaleFactor = 1.0 / (SLACK * ASy);
1766 if (ScaleFactor < 0.5)
1767 { for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
1768 ASy = 1.0 / SLACK;
1769 E *= ScaleFactor;
1770 }
1771
1772 /* Compute infinity-norm of T for O'Leary's estimate. */
1773 for (MaxY = 0.0, I = Size; I > 0; I--)
1774 if (MaxY < CMPLX_1_NORM(T[I])) MaxY = CMPLX_1_NORM(T[I]);
1775
1776 /*
1777 * Part 2. A* z = y where the * represents the transpose.
1778 * Recall that A = LU implies A* = U* L*.
1779 */
1780
1781 /* Forward elimination, U* v = y. */
1782 for (I = 1; I <= Size; I++)
1783 { pElement = Matrix->Diag[I]->NextInRow;
1784 while (pElement != NULL)
1785 { /* Cmplx expr: T[pElement->Col] -= T[I] * *pElement. */
1786 CMPLX_MULT_SUBT_ASSIGN( T[pElement->Col], T[I], *pElement );
1787 pElement = pElement->NextInRow;
1788 }
1789 if (CMPLX_1_NORM(T[I]) > SLACK)
1790 { ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
1791 for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
1792 ASy *= ScaleFactor;
1793 }
1794 }
1795
1796 /* Compute 1-norm of T, which now contains v, and scale ||T|| to 1/SLACK. */
1797 for (ASv = 0.0, I = Size; I > 0; I--) ASv += CMPLX_1_NORM(T[I]);
1798 ScaleFactor = 1.0 / (SLACK * ASv);
1799 if (ScaleFactor < 0.5)
1800 { for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
1801 ASy *= ScaleFactor;
1802 }
1803
1804 /* Backward Substitution, L* z = v. */
1805 for (I = Size; I >= 1; I--)
1806 { pPivot = Matrix->Diag[I];
1807 pElement = pPivot->NextInCol;
1808 while (pElement != NULL)
1809 { /* Cmplx expr: T[I] -= T[pElement->Row] * *pElement. */
1810 CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Row], *pElement );
1811 pElement = pElement->NextInCol;
1812 }
1813 CMPLX_MULT_ASSIGN( T[I], *pPivot );
1814 if (CMPLX_1_NORM(T[I]) > SLACK)
1815 { ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
1816 for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
1817 ASy *= ScaleFactor;
1818 }
1819 }
1820
1821 /* Compute 1-norm of T, which now contains z. */
1822 for (ASz = 0.0, I = Size; I > 0; I--) ASz += CMPLX_1_NORM(T[I]);
1823
1824 FREE( Tm );
1825
1826 Linpack = ASy / ASz;
1827 OLeary = E / MaxY;
1828 InvNormOfInverse = MIN( Linpack, OLeary );
1829 return InvNormOfInverse / NormOfMatrix;
1830 }
1831 #endif /* spCOMPLEX */
1832
1833
1834
1835
1836
1837 /*!
1838 * Computes the L-infinity norm of an unfactored matrix. It is a fatal
1839 * error to pass this routine a factored matrix.
1840 *
1841 * \return
1842 * The largest absolute row sum of matrix.
1843 *
1844 * \param eMatrix
1845 * Pointer to the matrix.
1846 */
1847
1848 spREAL
spNorm(spMatrix eMatrix)1849 spNorm( spMatrix eMatrix )
1850 {
1851 MatrixPtr Matrix = (MatrixPtr)eMatrix;
1852 register ElementPtr pElement;
1853 register int I;
1854 RealNumber Max = 0.0, AbsRowSum;
1855
1856 /* Begin `spNorm'. */
1857 ASSERT_IS_SPARSE( Matrix );
1858 ASSERT_NO_ERRORS( Matrix );
1859 ASSERT_IS_NOT_FACTORED( Matrix );
1860 if (NOT Matrix->RowsLinked) spcLinkRows( Matrix );
1861
1862 /* Compute row sums. */
1863 #if REAL
1864 if (NOT Matrix->Complex)
1865 { for (I = Matrix->Size; I > 0; I--)
1866 { pElement = Matrix->FirstInRow[I];
1867 AbsRowSum = 0.0;
1868 while (pElement != NULL)
1869 { AbsRowSum += ABS( pElement->Real );
1870 pElement = pElement->NextInRow;
1871 }
1872 if (Max < AbsRowSum) Max = AbsRowSum;
1873 }
1874 }
1875 #endif
1876 #if spCOMPLEX
1877 if (Matrix->Complex)
1878 { for (I = Matrix->Size; I > 0; I--)
1879 { pElement = Matrix->FirstInRow[I];
1880 AbsRowSum = 0.0;
1881 while (pElement != NULL)
1882 { AbsRowSum += CMPLX_1_NORM( *pElement );
1883 pElement = pElement->NextInRow;
1884 }
1885 if (Max < AbsRowSum) Max = AbsRowSum;
1886 }
1887 }
1888 #endif
1889 return Max;
1890 }
1891 #endif /* CONDITION */
1892
1893
1894
1895
1896
1897
1898 #if STABILITY
1899 /*!
1900 * This routine, along with spRoundoff(), are used to gauge the stability of a
1901 * factorization. If the factorization is determined to be too unstable,
1902 * then the matrix should be reordered. The routines compute quantities
1903 * that are needed in the computation of a bound on the error attributed
1904 * to any one element in the matrix during the factorization. In other
1905 * words, there is a matrix \f$ E = [e_{ij}] \f$ of error terms such that
1906 * \f$ A+E = LU \f$. This routine finds a bound on \f$ |e_{ij}| \f$.
1907 * Erisman & Reid [1] showed that \f$ |e_{ij}| < 3.01 u \rho m_{ij} \f$,
1908 * where \f$ u \f$ is the machine rounding unit,
1909 * \f$ \rho = \max a_{ij} \f$ where the max is taken over every row \f$ i \f$,
1910 * column \f$ j \f$, and step \f$ k \f$, and \f$ m_{ij} \f$ is the number
1911 * of multiplications required in the computation of \f$ l_{ij} \f$ if
1912 * \f$ i > j \f$ or \f$ u_{ij} \f$ otherwise. Barlow [2] showed that
1913 * \f$ \rho < \max_i || l_i ||_p \max_j || u_j ||_q \f$ where
1914 * \f$ 1/p + 1/q = 1 \f$.
1915 *
1916 * spLargestElement() finds the magnitude on the largest element in the
1917 * matrix. If the matrix has not yet been factored, the largest
1918 * element is found by direct search. If the matrix is factored, a
1919 * bound on the largest element in any of the reduced submatrices is
1920 * computed using Barlow with \f$ p = \infty \f$ and \f$ q = 1 \f$.
1921 * The ratio of these
1922 * two numbers is the growth, which can be used to determine if the
1923 * pivoting order is adequate. A large growth implies that
1924 * considerable error has been made in the factorization and that it
1925 * is probably a good idea to reorder the matrix. If a large growth
1926 * in encountered after using spFactor(), reconstruct the matrix and
1927 * refactor using spOrderAndFactor(). If a large growth is
1928 * encountered after using spOrderAndFactor(), refactor using
1929 * spOrderAndFactor() with the pivot threshold increased, say to 0.1.
1930 *
1931 * Using only the size of the matrix as an upper bound on \f$ m_{ij} \f$ and
1932 * Barlow's bound, the user can estimate the size of the matrix error
1933 * terms \f$ e_{ij} \f$ using the bound of Erisman and Reid. spRoundoff()
1934 * computes a tighter bound (with more work) based on work by Gear
1935 * [3], \f$ |e_{ij}| < 1.01 u \rho (t c^3 + (1 + t)c^2) \f$ where
1936 * \f$ t \f$ is the threshold and \f$ c \f$ is the maximum number of
1937 * off-diagonal elements in any row of \f$ L \f$. The expensive part
1938 * of computing this bound is determining the maximum number of
1939 * off-diagonals in \f$ L \f$, which changes
1940 * only when the order of the matrix changes. This number is computed
1941 * and saved, and only recomputed if the matrix is reordered.
1942 *
1943 * [1] A. M. Erisman, J. K. Reid. Monitoring the stability of the
1944 * triangular factorization of a sparse matrix. Numerische
1945 * Mathematik. Vol. 22, No. 3, 1974, pp 183-186.
1946 *
1947 * [2] J. L. Barlow. A note on monitoring the stability of triangular
1948 * decomposition of sparse matrices. "SIAM Journal of Scientific
1949 * and Statistical Computing." Vol. 7, No. 1, January 1986, pp 166-168.
1950 *
1951 * [3] I. S. Duff, A. M. Erisman, J. K. Reid. "Direct Methods for Sparse
1952 * Matrices." Oxford 1986. pp 99.
1953 *
1954 * \return
1955 * If matrix is not factored, returns the magnitude of the largest element in
1956 * the matrix. If the matrix is factored, a bound on the magnitude of the
1957 * largest element in any of the reduced submatrices is returned.
1958 *
1959 * \param eMatrix
1960 * Pointer to the matrix.
1961 */
1962
1963 spREAL
spLargestElement(spMatrix eMatrix)1964 spLargestElement( spMatrix eMatrix )
1965 {
1966 MatrixPtr Matrix = (MatrixPtr)eMatrix;
1967 register int I;
1968 RealNumber Mag, AbsColSum, Max = 0.0, MaxRow = 0.0, MaxCol = 0.0;
1969 RealNumber Pivot;
1970 #if spCOMPLEX
1971 ComplexNumber cPivot;
1972 #endif
1973 register ElementPtr pElement, pDiag;
1974
1975 /* Begin `spLargestElement'. */
1976 ASSERT_IS_SPARSE( Matrix );
1977
1978 #if REAL
1979 if (Matrix->Factored AND NOT Matrix->Complex)
1980 { if (Matrix->Error == spSINGULAR) return 0.0;
1981
1982 /* Find the bound on the size of the largest element over all factorization. */
1983 for (I = 1; I <= Matrix->Size; I++)
1984 { pDiag = Matrix->Diag[I];
1985
1986 /* Lower triangular matrix. */
1987 Pivot = 1.0 / pDiag->Real;
1988 Mag = ABS( Pivot );
1989 if (Mag > MaxRow) MaxRow = Mag;
1990 pElement = Matrix->FirstInRow[I];
1991 while (pElement != pDiag)
1992 { Mag = ABS( pElement->Real );
1993 if (Mag > MaxRow) MaxRow = Mag;
1994 pElement = pElement->NextInRow;
1995 }
1996
1997 /* Upper triangular matrix. */
1998 pElement = Matrix->FirstInCol[I];
1999 AbsColSum = 1.0; /* Diagonal of U is unity. */
2000 while (pElement != pDiag)
2001 { AbsColSum += ABS( pElement->Real );
2002 pElement = pElement->NextInCol;
2003 }
2004 if (AbsColSum > MaxCol) MaxCol = AbsColSum;
2005 }
2006 }
2007 else if (NOT Matrix->Complex)
2008 { for (I = 1; I <= Matrix->Size; I++)
2009 { pElement = Matrix->FirstInCol[I];
2010 while (pElement != NULL)
2011 { Mag = ABS( pElement->Real );
2012 if (Mag > Max) Max = Mag;
2013 pElement = pElement->NextInCol;
2014 }
2015 }
2016 return Max;
2017 }
2018 #endif
2019 #if spCOMPLEX
2020 if (Matrix->Factored AND Matrix->Complex)
2021 { if (Matrix->Error == spSINGULAR) return 0.0;
2022
2023 /* Find the bound on the size of the largest element over all factorization. */
2024 for (I = 1; I <= Matrix->Size; I++)
2025 { pDiag = Matrix->Diag[I];
2026
2027 /* Lower triangular matrix. */
2028 CMPLX_RECIPROCAL( cPivot, *pDiag );
2029 Mag = CMPLX_INF_NORM( cPivot );
2030 if (Mag > MaxRow) MaxRow = Mag;
2031 pElement = Matrix->FirstInRow[I];
2032 while (pElement != pDiag)
2033 { Mag = CMPLX_INF_NORM( *pElement );
2034 if (Mag > MaxRow) MaxRow = Mag;
2035 pElement = pElement->NextInRow;
2036 }
2037
2038 /* Upper triangular matrix. */
2039 pElement = Matrix->FirstInCol[I];
2040 AbsColSum = 1.0; /* Diagonal of U is unity. */
2041 while (pElement != pDiag)
2042 { AbsColSum += CMPLX_INF_NORM( *pElement );
2043 pElement = pElement->NextInCol;
2044 }
2045 if (AbsColSum > MaxCol) MaxCol = AbsColSum;
2046 }
2047 }
2048 else if (Matrix->Complex)
2049 { for (I = 1; I <= Matrix->Size; I++)
2050 { pElement = Matrix->FirstInCol[I];
2051 while (pElement != NULL)
2052 { Mag = CMPLX_INF_NORM( *pElement );
2053 if (Mag > Max) Max = Mag;
2054 pElement = pElement->NextInCol;
2055 }
2056 }
2057 return Max;
2058 }
2059 #endif
2060 return MaxRow * MaxCol;
2061 }
2062
2063
2064
2065
2066 /*!
2067 * This routine, along with spLargestElement(), are used to gauge the
2068 * stability of a factorization. See description of spLargestElement()
2069 * for more information.
2070 *
2071 * \return
2072 * Returns a bound on the magnitude of the largest element in
2073 * \f$ E = A - LU \f$.
2074 *
2075 * \param eMatrix
2076 * Pointer to the matrix.
2077 * \param Rho
2078 * The bound on the magnitude of the largest element in any of the
2079 * reduced submatrices. This is the number computed by the function
2080 * spLargestElement() when given a factored matrix. If this number is
2081 * negative, the bound will be computed automatically.
2082 */
2083
2084 spREAL
spRoundoff(spMatrix eMatrix,spREAL Rho)2085 spRoundoff(
2086 spMatrix eMatrix,
2087 spREAL Rho
2088 )
2089 {
2090 MatrixPtr Matrix = (MatrixPtr)eMatrix;
2091 register ElementPtr pElement;
2092 register int Count, I, MaxCount = 0;
2093 RealNumber Reid, Gear;
2094
2095 /* Begin `spRoundoff'. */
2096 ASSERT_IS_SPARSE( Matrix );
2097 ASSERT_NO_ERRORS( Matrix );
2098 ASSERT_IS_FACTORED( Matrix );
2099
2100 /* Compute Barlow's bound if it is not given. */
2101 if (Rho < 0.0) Rho = spLargestElement( eMatrix );
2102
2103 /* Find the maximum number of off-diagonals in L if not previously computed. */
2104 if (Matrix->MaxRowCountInLowerTri < 0)
2105 { for (I = Matrix->Size; I > 0; I--)
2106 { pElement = Matrix->FirstInRow[I];
2107 Count = 0;
2108 while (pElement->Col < I)
2109 { Count++;
2110 pElement = pElement->NextInRow;
2111 }
2112 if (Count > MaxCount) MaxCount = Count;
2113 }
2114 Matrix->MaxRowCountInLowerTri = MaxCount;
2115 }
2116 else MaxCount = Matrix->MaxRowCountInLowerTri;
2117
2118 /* Compute error bound. */
2119 Gear = 1.01*((MaxCount + 1) * Matrix->RelThreshold + 1.0) * SQR(MaxCount);
2120 Reid = 3.01 * Matrix->Size;
2121
2122 if (Gear < Reid)
2123 return (MACHINE_RESOLUTION * Rho * Gear);
2124 else
2125 return (MACHINE_RESOLUTION * Rho * Reid);
2126 }
2127 #endif
2128
2129
2130
2131
2132
2133
2134
2135 #if DOCUMENTATION
2136 /*!
2137 * This routine prints a short message describing the error error state
2138 * of sparse. No message is produced if there is no error.
2139 * The error state is cleared.
2140 *
2141 * \param eMatrix
2142 * Matrix for which the error message is to be printed.
2143 * \param Stream
2144 * Stream to which the error message is to be printed.
2145 * \param Originator
2146 * Name of originator of error message. If NULL, `sparse' is used.
2147 * If zero-length string, no originator is printed.
2148 */
2149 #endif /* DOCUMENTATION */
2150