epilectrik
/
pyNgSpice
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
5787 Commits
1 Branch
0 Tags
36 MiB
 
 
 
 
 
 
Tree: ce2e704f56
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'ce2e704f56'
${ noResults }
pyNgSpice/src/maths/sparse/sputils.c

2180 lines
65 KiB
Raw Blame History

/*
* MATRIX UTILITY MODULE
*
* Author: Advising professor:
* Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
* UC Berkeley
*
* This file contains various optional utility routines.
*
* >>> User accessible functions contained in this file:
* spMNA_Preorder
* spScale
* spMultiply
* spMultTransposed
* spDeterminant
* spStripFills
* spDeleteRowAndCol
* spPseudoCondition
* spCondition
* spNorm
* spLargestElement
* spRoundoff
* spErrorMessage
*
* >>> Other functions contained in this file:
* CountTwins
* SwapCols
* ScaleComplexMatrix
* ComplexMatrixMultiply
* ComplexCondition
*/
/*
* Revision and copyright information.
*
* Copyright (c) 1985,86,87,88,89,90
* by Kenneth S. Kundert and the University of California.
*
* Permission to use, copy, modify, and distribute this software and
* its documentation for any purpose and without fee is hereby granted,
* provided that the copyright notices appear in all copies and
* supporting documentation and that the authors and the University of
* California are properly credited. The authors and the University of
* California make no representations as to the suitability of this
* software for any purpose. It is provided `as is', without express
* or implied warranty.
*/
/*
* IMPORTS
*
* >>> Import descriptions:
* spConfig.h
* Macros that customize the sparse matrix routines.
* spMatrix.h
* Macros and declarations to be imported by the user.
* spDefs.h
* Matrix type and macro definitions for the sparse matrix routines.
*/
#include <assert.h>
#include <stdlib.h>
#define spINSIDE_SPARSE
#include "spconfig.h"
#include "ngspice/spmatrix.h"
#include "spdefs.h"
/* Function declarations */
static int CountTwins( MatrixPtr, int, ElementPtr*, ElementPtr* );
static void SwapCols( MatrixPtr, ElementPtr, ElementPtr );
static void ComplexMatrixMultiply( MatrixPtr, RealVector, RealVector,
RealVector, RealVector );
static void ComplexTransposedMatrixMultiply( MatrixPtr, RealVector, RealVector,
RealVector, RealVector);
#if MODIFIED_NODAL
/*
* PREORDER MODIFIED NODE ADMITTANCE MATRIX TO REMOVE ZEROS FROM DIAGONAL
*
* This routine massages modified node admittance matrices to remove
* zeros from the diagonal. It takes advantage of the fact that the
* row and column associated with a zero diagonal usually have
* structural ones placed symmetricly. This routine should be used
* only on modified node admittance matrices and should be executed
* after the matrix has been built but before the factorization
* begins. It should be executed for the initial factorization only
* and should be executed before the rows have been linked. Thus it
* should be run before using spScale(), spMultiply(),
* spDeleteRowAndCol(), or spNorm().
*
* This routine exploits the fact that the structural ones are placed
* in the matrix in symmetric twins. For example, the stamps for
* grounded and a floating voltage sources are
* grounded: floating:
* [ x x 1 ] [ x x 1 ]
* [ x x ] [ x x -1 ]
* [ 1 ] [ 1 -1 ]
* Notice for the grounded source, there is one set of twins, and for
* the floating, there are two sets. We remove the zero from the diagonal
* by swapping the rows associated with a set of twins. For example:
* grounded: floating 1: floating 2:
* [ 1 ] [ 1 -1 ] [ x x 1 ]
* [ x x ] [ x x -1 ] [ 1 -1 ]
* [ x x 1 ] [ x x 1 ] [ x x -1 ]
*
* It is important to deal with any zero diagonals that only have one
* set of twins before dealing with those that have more than one because
* swapping row destroys the symmetry of any twins in the rows being
* swapped, which may limit future moves. Consider
* [ x x 1 ]
* [ x x -1 1 ]
* [ 1 -1 ]
* [ 1 ]
* There is one set of twins for diagonal 4 and two for diagonal 3.
* Dealing with diagonal 4 first requires swapping rows 2 and 4.
* [ x x 1 ]
* [ 1 ]
* [ 1 -1 ]
* [ x x -1 1 ]
* We can now deal with diagonal 3 by swapping rows 1 and 3.
* [ 1 -1 ]
* [ 1 ]
* [ x x 1 ]
* [ x x -1 1 ]
* And we are done, there are no zeros left on the diagonal. However, if
* we originally dealt with diagonal 3 first, we could swap rows 2 and 3
* [ x x 1 ]
* [ 1 -1 ]
* [ x x -1 1 ]
* [ 1 ]
* Diagonal 4 no longer has a symmetric twin and we cannot continue.
*
* So we always take care of lone twins first. When none remain, we
* choose arbitrarily a set of twins for a diagonal with more than one set
* and swap the rows corresponding to that twin. We then deal with any
* lone twins that were created and repeat the procedure until no
* zero diagonals with symmetric twins remain.
*
* In this particular implementation, columns are swapped rather than rows.
* The algorithm used in this function was developed by Ken Kundert and
* Tom Quarles.
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix to be preordered.
*
* >>> Local variables;
* J (int)
* Column with zero diagonal being currently considered.
* pTwin1 (ElementPtr)
* Pointer to the twin found in the column belonging to the zero diagonal.
* pTwin2 (ElementPtr)
* Pointer to the twin found in the row belonging to the zero diagonal.
* belonging to the zero diagonal.
* AnotherPassNeeded (int)
* Flag indicating that at least one zero diagonal with symmetric twins
* remain.
* StartAt (int)
* Column number of first zero diagonal with symmetric twins.
* Swapped (int)
* Flag indicating that columns were swapped on this pass.
* Twins (int)
* Number of symmetric twins corresponding to current zero diagonal.
*/
void
spMNA_Preorder(MatrixPtr Matrix)
{
int J, Size;
ElementPtr pTwin1, pTwin2;
int Twins, StartAt = 1;
int Swapped, AnotherPassNeeded;
/* Begin `spMNA_Preorder'. */
assert( IS_VALID(Matrix) && !Matrix->Factored );
if (Matrix->RowsLinked) return;
Size = Matrix->Size;
Matrix->Reordered = YES;
do
{
AnotherPassNeeded = Swapped = NO;
/* Search for zero diagonals with lone twins. */
for (J = StartAt; J <= Size; J++)
{
if (Matrix->Diag[J] == NULL)
{
Twins = CountTwins( Matrix, J, &pTwin1, &pTwin2 );
if (Twins == 1)
{
/* Lone twins found, swap rows. */
SwapCols( Matrix, pTwin1, pTwin2 );
Swapped = YES;
}
else if ((Twins > 1) && !AnotherPassNeeded)
{
AnotherPassNeeded = YES;
StartAt = J;
}
}
}
/* All lone twins are gone, look for zero diagonals with multiple twins. */
if (AnotherPassNeeded)
{
for (J = StartAt; !Swapped && (J <= Size); J++)
{
if (Matrix->Diag[J] == NULL)
{
Twins = CountTwins( Matrix, J, &pTwin1, &pTwin2 );
SwapCols( Matrix, pTwin1, pTwin2 );
Swapped = YES;
}
}
}
} while (AnotherPassNeeded);
return;
}
/*
* COUNT TWINS
*
* This function counts the number of symmetric twins associated with
* a zero diagonal and returns one set of twins if any exist. The
* count is terminated early at two.
*/
static int
CountTwins( MatrixPtr Matrix, int Col, ElementPtr *ppTwin1, ElementPtr *ppTwin2 )
{
int Row, Twins = 0;
ElementPtr pTwin1, pTwin2;
/* Begin `CountTwins'. */
pTwin1 = Matrix->FirstInCol[Col];
while (pTwin1 != NULL)
{
if (ABS(pTwin1->Real) == 1.0)
{
Row = pTwin1->Row;
pTwin2 = Matrix->FirstInCol[Row];
while ((pTwin2 != NULL) && (pTwin2->Row != Col))
pTwin2 = pTwin2->NextInCol;
if ((pTwin2 != NULL) && (ABS(pTwin2->Real) == 1.0))
{
/* Found symmetric twins. */
if (++Twins >= 2) return Twins;
(*ppTwin1 = pTwin1)->Col = Col;
(*ppTwin2 = pTwin2)->Col = Row;
}
}
pTwin1 = pTwin1->NextInCol;
}
return Twins;
}
/*
* SWAP COLUMNS
*
* This function swaps two columns and is applicable before the rows are
* linked.
*/
static void
SwapCols( MatrixPtr Matrix, ElementPtr pTwin1, ElementPtr pTwin2 )
{
int Col1 = pTwin1->Col, Col2 = pTwin2->Col;
/* Begin `SwapCols'. */
SWAP (ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]);
SWAP (int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]);
#if TRANSLATE
Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col2]]=Col2;
Matrix->ExtToIntColMap[Matrix->IntToExtColMap[Col1]]=Col1;
#endif
Matrix->Diag[Col1] = pTwin2;
Matrix->Diag[Col2] = pTwin1;
Matrix->NumberOfInterchangesIsOdd = !Matrix->NumberOfInterchangesIsOdd;
return;
}
#endif /* MODIFIED_NODAL */
#if SCALING
/*
* SCALE MATRIX
*
* This function scales the matrix to enhance the possibility of
* finding a good pivoting order. Note that scaling enhances accuracy
* of the solution only if it affects the pivoting order, so it makes
* no sense to scale the matrix before spFactor(). If scaling is
* desired it should be done before spOrderAndFactor(). There
* are several things to take into account when choosing the scale
* factors. First, the scale factors are directly multiplied against
* the elements in the matrix. To prevent roundoff, each scale factor
* should be equal to an integer power of the number base of the
* machine. Since most machines operate in base two, scale factors
* should be a power of two. Second, the matrix should be scaled such
* that the matrix of element uncertainties is equilibrated. Third,
* this function multiplies the scale factors by the elements, so if
* one row tends to have uncertainties 1000 times smaller than the
* other rows, then its scale factor should be 1024, not 1/1024.
* Fourth, to save time, this function does not scale rows or columns
* if their scale factors are equal to one. Thus, the scale factors
* should be normalized to the most common scale factor. Rows and
* columns should be normalized separately. For example, if the size
* of the matrix is 100 and 10 rows tend to have uncertainties near
* 1e-6 and the remaining 90 have uncertainties near 1e-12, then the
* scale factor for the 10 should be 1/1,048,576 and the scale factors
* for the remaining 90 should be 1. Fifth, since this routine
* directly operates on the matrix, it is necessary to apply the scale
* factors to the RHS and Solution vectors. It may be easier to
* simply use spOrderAndFactor() on a scaled matrix to choose the
* pivoting order, and then throw away the matrix. Subsequent
* factorizations, performed with spFactor(), will not need to have
* the RHS and Solution vectors descaled. Lastly, this function
* should not be executed before the function spMNA_Preorder.
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix to be scaled.
* SolutionScaleFactors <input> (RealVector)
* The array of Solution scale factors. These factors scale the columns.
* All scale factors are real valued.
* RHS_ScaleFactors <input> (RealVector)
* The array of RHS scale factors. These factors scale the rows.
* All scale factors are real valued.
*
* >>> Local variables:
* lSize (int)
* Local version of the size of the matrix.
* pElement (ElementPtr)
* Pointer to an element in the matrix.
* pExtOrder (int *)
* Pointer into either IntToExtRowMap or IntToExtColMap vector. Used to
* compensate for any row or column swaps that have been performed.
* ScaleFactor (RealNumber)
* The scale factor being used on the current row or column.
*/
void
spScale(MatrixPtr Matrix, RealVector RHS_ScaleFactors, RealVector SolutionScaleFactors)
{
ElementPtr pElement;
int I, lSize, *pExtOrder;
RealNumber ScaleFactor;
void ScaleComplexMatrix();
/* Begin `spScale'. */
assert( IS_VALID(Matrix) && !Matrix->Factored );
if (!Matrix->RowsLinked) spcLinkRows( Matrix );
if (Matrix->Complex)
{
ScaleComplexMatrix( Matrix, RHS_ScaleFactors, SolutionScaleFactors );
return;
}
lSize = Matrix->Size;
/* Scale Rows */
pExtOrder = &Matrix->IntToExtRowMap[1];
for (I = 1; I <= lSize; I++)
{
if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0)
{
pElement = Matrix->FirstInRow[I];
while (pElement != NULL)
{
pElement->Real *= ScaleFactor;
pElement = pElement->NextInRow;
}
}
}
/* Scale Columns */
pExtOrder = &Matrix->IntToExtColMap[1];
for (I = 1; I <= lSize; I++)
{
if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0)
{
pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{
pElement->Real *= ScaleFactor;
pElement = pElement->NextInCol;
}
}
}
return;
}
#endif /* SCALING */
#if SCALING
/*
* SCALE COMPLEX MATRIX
*
* This function scales the matrix to enhance the possibility of
* finding a good pivoting order. Note that scaling enhances accuracy
* of the solution only if it affects the pivoting order, so it makes
* no sense to scale the matrix before spFactor(). If scaling is
* desired it should be done before spOrderAndFactor(). There
* are several things to take into account when choosing the scale
* factors. First, the scale factors are directly multiplied against
* the elements in the matrix. To prevent roundoff, each scale factor
* should be equal to an integer power of the number base of the
* machine. Since most machines operate in base two, scale factors
* should be a power of two. Second, the matrix should be scaled such
* that the matrix of element uncertainties is equilibrated. Third,
* this function multiplies the scale factors by the elements, so if
* one row tends to have uncertainties 1000 times smaller than the
* other rows, then its scale factor should be 1024, not 1/1024.
* Fourth, to save time, this function does not scale rows or columns
* if their scale factors are equal to one. Thus, the scale factors
* should be normalized to the most common scale factor. Rows and
* columns should be normalized separately. For example, if the size
* of the matrix is 100 and 10 rows tend to have uncertainties near
* 1e-6 and the remaining 90 have uncertainties near 1e-12, then the
* scale factor for the 10 should be 1/1,048,576 and the scale factors
* for the remaining 90 should be 1. Fifth, since this routine
* directly operates on the matrix, it is necessary to apply the scale
* factors to the RHS and Solution vectors. It may be easier to
* simply use spOrderAndFactor() on a scaled matrix to choose the
* pivoting order, and then throw away the matrix. Subsequent
* factorizations, performed with spFactor(), will not need to have
* the RHS and Solution vectors descaled. Lastly, this function
* should not be executed before the function spMNA_Preorder.
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix to be scaled.
* SolutionScaleFactors <input> (RealVector)
* The array of Solution scale factors. These factors scale the columns.
* All scale factors are real valued.
* RHS_ScaleFactors <input> (RealVector)
* The array of RHS scale factors. These factors scale the rows.
* All scale factors are real valued.
*
* >>> Local variables:
* lSize (int)
* Local version of the size of the matrix.
* pElement (ElementPtr)
* Pointer to an element in the matrix.
* pExtOrder (int *)
* Pointer into either IntToExtRowMap or IntToExtColMap vector. Used to
* compensate for any row or column swaps that have been performed.
* ScaleFactor (RealNumber)
* The scale factor being used on the current row or column.
*/
static void
ScaleComplexMatrix( Matrix, RHS_ScaleFactors, SolutionScaleFactors )
MatrixPtr Matrix;
RealVector RHS_ScaleFactors, SolutionScaleFactors;
{
ElementPtr pElement;
int I, lSize, *pExtOrder;
RealNumber ScaleFactor;
/* Begin `ScaleComplexMatrix'. */
lSize = Matrix->Size;
/* Scale Rows */
pExtOrder = &Matrix->IntToExtRowMap[1];
for (I = 1; I <= lSize; I++)
{
if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0)
{
pElement = Matrix->FirstInRow[I];
while (pElement != NULL)
{
pElement->Real *= ScaleFactor;
pElement->Imag *= ScaleFactor;
pElement = pElement->NextInRow;
}
}
}
/* Scale Columns */
pExtOrder = &Matrix->IntToExtColMap[1];
for (I = 1; I <= lSize; I++)
{
if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0)
{
pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{
pElement->Real *= ScaleFactor;
pElement->Imag *= ScaleFactor;
pElement = pElement->NextInCol;
}
}
}
return;
}
#endif /* SCALING */
#if MULTIPLICATION
/*
* MATRIX MULTIPLICATION
*
* Multiplies matrix by solution vector to find source vector.
* Assumes matrix has not been factored. This routine can be used
* as a test to see if solutions are correct. It should not be used
* before spMNA_Preorder().
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix.
* RHS <output> (RealVector)
* RHS is the right hand side. This is what is being solved for.
* Solution <input> (RealVector)
* Solution is the vector being multiplied by the matrix.
* iRHS <output> (RealVector)
* iRHS is the imaginary portion of the right hand side. This is
* what is being solved for.
* iSolution <input> (RealVector)
* iSolution is the imaginary portion of the vector being multiplied
* by the matrix.
*/
void
spMultiply(MatrixPtr Matrix, RealVector RHS, RealVector Solution,
RealVector iRHS, RealVector iSolution)
{
ElementPtr pElement;
RealVector Vector;
RealNumber Sum;
int I, *pExtOrder;
/* Begin `spMultiply'. */
assert( IS_SPARSE( Matrix ) && !Matrix->Factored );
if (!Matrix->RowsLinked)
spcLinkRows(Matrix);
if (!Matrix->InternalVectorsAllocated)
spcCreateInternalVectors( Matrix );
if (Matrix->Complex)
{
ComplexMatrixMultiply( Matrix, RHS, Solution , iRHS, iSolution );
return;
}
/* Initialize Intermediate vector with reordered Solution vector. */
Vector = Matrix->Intermediate;
pExtOrder = &Matrix->IntToExtColMap[Matrix->Size];
for (I = Matrix->Size; I > 0; I--)
Vector[I] = Solution[*(pExtOrder--)];
pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size];
for (I = Matrix->Size; I > 0; I--)
{
pElement = Matrix->FirstInRow[I];
Sum = 0.0;
while (pElement != NULL)
{
Sum += pElement->Real * Vector[pElement->Col];
pElement = pElement->NextInRow;
}
RHS[*pExtOrder--] = Sum;
}
return;
}
#endif /* MULTIPLICATION */
#if MULTIPLICATION
/*
* COMPLEX MATRIX MULTIPLICATION
*
* Multiplies matrix by solution vector to find source vector.
* Assumes matrix has not been factored. This routine can be used
* as a test to see if solutions are correct.
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix.
* RHS <output> (RealVector)
* RHS is the right hand side. This is what is being solved for.
* Solution <input> (RealVector)
* Solution is the vector being multiplied by the matrix.
* iRHS <output> (RealVector)
* iRHS is the imaginary portion of the right hand side. This is
* what is being solved for.
* iSolution <input> (RealVector)
* iSolution is the imaginary portion of the vector being multiplied
* by the matrix.
*/
static void
ComplexMatrixMultiply( MatrixPtr Matrix, RealVector RHS, RealVector Solution , RealVector iRHS, RealVector iSolution )
{
ElementPtr pElement;
ComplexVector Vector;
ComplexNumber Sum;
int I, *pExtOrder;
/* Begin `ComplexMatrixMultiply'. */
/* Initialize Intermediate vector with reordered Solution vector. */
Vector = (ComplexVector)Matrix->Intermediate;
pExtOrder = &Matrix->IntToExtColMap[Matrix->Size];
for (I = Matrix->Size; I > 0; I--)
{
Vector[I].Real = Solution[*pExtOrder];
Vector[I].Imag = iSolution[*(pExtOrder--)];
}
pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size];
for (I = Matrix->Size; I > 0; I--)
{
pElement = Matrix->FirstInRow[I];
Sum.Real = Sum.Imag = 0.0;
while (pElement != NULL)
{
/* Cmplx expression : Sum += Element * Vector[Col] */
CMPLX_MULT_ADD_ASSIGN( Sum, *pElement, Vector[pElement->Col] );
pElement = pElement->NextInRow;
}
RHS[*pExtOrder] = Sum.Real;
iRHS[*pExtOrder--] = Sum.Imag;
}
return;
}
#endif /* MULTIPLICATION */
#if MULTIPLICATION && TRANSPOSE
/*
* TRANSPOSED MATRIX MULTIPLICATION
*
* Multiplies transposed matrix by solution vector to find source vector.
* Assumes matrix has not been factored. This routine can be used
* as a test to see if solutions are correct. It should not be used
* before spMNA_Preorder().
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix.
* RHS <output> (RealVector)
* RHS is the right hand side. This is what is being solved for.
* Solution <input> (RealVector)
* Solution is the vector being multiplied by the matrix.
* iRHS <output> (RealVector)
* iRHS is the imaginary portion of the right hand side. This is
* what is being solved for.
* iSolution <input> (RealVector)
* iSolution is the imaginary portion of the vector being multiplied
* by the matrix.
*/
void
spMultTransposed(MatrixPtr Matrix, RealVector RHS, RealVector Solution,
RealVector iRHS, RealVector iSolution)
{
ElementPtr pElement;
RealVector Vector;
RealNumber Sum;
int I, *pExtOrder;
/* Begin `spMultTransposed'. */
assert( IS_SPARSE( Matrix ) && !Matrix->Factored );
if (!Matrix->InternalVectorsAllocated)
spcCreateInternalVectors( Matrix );
if (Matrix->Complex)
{
ComplexTransposedMatrixMultiply( Matrix, RHS, Solution , iRHS, iSolution );
return;
}
/* Initialize Intermediate vector with reordered Solution vector. */
Vector = Matrix->Intermediate;
pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size];
for (I = Matrix->Size; I > 0; I--)
Vector[I] = Solution[*(pExtOrder--)];
pExtOrder = &Matrix->IntToExtColMap[Matrix->Size];
for (I = Matrix->Size; I > 0; I--)
{
pElement = Matrix->FirstInCol[I];
Sum = 0.0;
while (pElement != NULL)
{
Sum += pElement->Real * Vector[pElement->Row];
pElement = pElement->NextInCol;
}
RHS[*pExtOrder--] = Sum;
}
return;
}
#endif /* MULTIPLICATION && TRANSPOSE */
#if MULTIPLICATION && TRANSPOSE
/*
* COMPLEX TRANSPOSED MATRIX MULTIPLICATION
*
* Multiplies transposed matrix by solution vector to find source vector.
* Assumes matrix has not been factored. This routine can be used
* as a test to see if solutions are correct.
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix.
* RHS <output> (RealVector)
* RHS is the right hand side. This is what is being solved for.
* Solution <input> (RealVector)
* Solution is the vector being multiplied by the matrix.
* iRHS <output> (RealVector)
* iRHS is the imaginary portion of the right hand side. This is
* what is being solved for.
* iSolution <input> (RealVector)
* iSolution is the imaginary portion of the vector being multiplied
* by the matrix.
*/
static void
ComplexTransposedMatrixMultiply( MatrixPtr Matrix, RealVector RHS, RealVector Solution , RealVector iRHS, RealVector iSolution )
{
ElementPtr pElement;
ComplexVector Vector;
ComplexNumber Sum;
int I, *pExtOrder;
/* Begin `ComplexMatrixMultiply'. */
/* Initialize Intermediate vector with reordered Solution vector. */
Vector = (ComplexVector)Matrix->Intermediate;
pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size];
for (I = Matrix->Size; I > 0; I--)
{
Vector[I].Real = Solution[*pExtOrder];
Vector[I].Imag = iSolution[*(pExtOrder--)];
}
pExtOrder = &Matrix->IntToExtColMap[Matrix->Size];
for (I = Matrix->Size; I > 0; I--)
{
pElement = Matrix->FirstInCol[I];
Sum.Real = Sum.Imag = 0.0;
while (pElement != NULL)
{
/* Cmplx expression : Sum += Element * Vector[Row] */
CMPLX_MULT_ADD_ASSIGN( Sum, *pElement, Vector[pElement->Row] );
pElement = pElement->NextInCol;
}
RHS[*pExtOrder] = Sum.Real;
iRHS[*pExtOrder--] = Sum.Imag;
}
return;
}
#endif /* MULTIPLICATION && TRANSPOSE */
#if DETERMINANT
/*
* CALCULATE DETERMINANT
*
* This routine in capable of calculating the determinant of the
* matrix once the LU factorization has been performed. Hence, only
* use this routine after spFactor() and before spClear().
* The determinant equals the product of all the diagonal elements of
* the lower triangular matrix L, except that this product may need
* negating. Whether the product or the negative product equals the
* determinant is determined by the number of row and column
* interchanges performed. Note that the determinants of matrices can
* be very large or very small. On large matrices, the determinant
* can be far larger or smaller than can be represented by a floating
* point number. For this reason the determinant is scaled to a
* reasonable value and the logarithm of the scale factor is returned.
*
* >>> Arguments:
* Matrix <input> (char *)
* A pointer to the matrix for which the determinant is desired.
* pExponent <output> (int *)
* The logarithm base 10 of the scale factor for the determinant. To find
* the actual determinant, Exponent should be added to the exponent of
* Determinant.
* pDeterminant <output> (RealNumber *)
* The real portion of the determinant. This number is scaled to be
* greater than or equal to 1.0 and less than 10.0.
* piDeterminant <output> (RealNumber *)
* The imaginary portion of the determinant. When the matrix is real
* this pointer need not be supplied, nothing will be returned. This
* number is scaled to be greater than or equal to 1.0 and less than 10.0.
*
* >>> Local variables:
* Norm (RealNumber)
* L-infinity norm of a complex number.
* Size (int)
* Local storage for Matrix->Size. Placed in a for speed.
* Temp (RealNumber)
* Temporary storage for real portion of determinant.
*/
void
spDeterminant(MatrixPtr Matrix, int *pExponent, RealNumber *pDeterminant,
RealNumber *piDeterminant)
{
int I, Size;
RealNumber Norm, nr, ni;
ComplexNumber Pivot, cDeterminant;
#define NORM(a) (nr = ABS((a).Real), ni = ABS((a).Imag), MAX (nr,ni))
/* Begin `spDeterminant'. */
assert( IS_SPARSE( Matrix ) && IS_FACTORED(Matrix) );
*pExponent = 0;
if (Matrix->Error == spSINGULAR)
{
*pDeterminant = 0.0;
if (Matrix->Complex) *piDeterminant = 0.0;
return;
}
Size = Matrix->Size;
I = 0;
if (Matrix->Complex) /* Complex Case. */
{
cDeterminant.Real = 1.0;
cDeterminant.Imag = 0.0;
while (++I <= Size)
{
CMPLX_RECIPROCAL( Pivot, *Matrix->Diag[I] );
CMPLX_MULT_ASSIGN( cDeterminant, Pivot );
/* Scale Determinant. */
Norm = NORM( cDeterminant );
if (Norm != 0.0)
{
while (Norm >= 1.0e12)
{
cDeterminant.Real *= 1.0e-12;
cDeterminant.Imag *= 1.0e-12;
*pExponent += 12;
Norm = NORM( cDeterminant );
}
while (Norm < 1.0e-12)
{
cDeterminant.Real *= 1.0e12;
cDeterminant.Imag *= 1.0e12;
*pExponent -= 12;
Norm = NORM( cDeterminant );
}
}
}
/* Scale Determinant again, this time to be between 1.0 <= x < 10.0. */
Norm = NORM( cDeterminant );
if (Norm != 0.0)
{
while (Norm >= 10.0)
{
cDeterminant.Real *= 0.1;
cDeterminant.Imag *= 0.1;
(*pExponent)++;
Norm = NORM( cDeterminant );
}
while (Norm < 1.0)
{
cDeterminant.Real *= 10.0;
cDeterminant.Imag *= 10.0;
(*pExponent)--;
Norm = NORM( cDeterminant );
}
}
if (Matrix->NumberOfInterchangesIsOdd)
CMPLX_NEGATE( cDeterminant );
*pDeterminant = cDeterminant.Real;
*piDeterminant = cDeterminant.Imag;
}
else
{
/* Real Case. */
*pDeterminant = 1.0;
while (++I <= Size)
{
*pDeterminant /= Matrix->Diag[I]->Real;
/* Scale Determinant. */
if (*pDeterminant != 0.0)
{
while (ABS(*pDeterminant) >= 1.0e12)
{
*pDeterminant *= 1.0e-12;
*pExponent += 12;
}
while (ABS(*pDeterminant) < 1.0e-12)
{
*pDeterminant *= 1.0e12;
*pExponent -= 12;
}
}
}
/* Scale Determinant again, this time to be between 1.0 <= x <
10.0. */
if (*pDeterminant != 0.0)
{
while (ABS(*pDeterminant) >= 10.0)
{
*pDeterminant *= 0.1;
(*pExponent)++;
}
while (ABS(*pDeterminant) < 1.0)
{
*pDeterminant *= 10.0;
(*pExponent)--;
}
}
if (Matrix->NumberOfInterchangesIsOdd)
*pDeterminant = -*pDeterminant;
}
}
#endif /* DETERMINANT */
#if STRIP
/*
* STRIP FILL-INS FROM MATRIX
*
* Strips the matrix of all fill-ins.
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix to be stripped.
*
* >>> Local variables:
* pElement (ElementPtr)
* Pointer that is used to step through the matrix.
* ppElement (ElementPtr *)
* Pointer to the location of an ElementPtr. This location will be
* updated if a fill-in is stripped from the matrix.
* pFillin (ElementPtr)
* Pointer used to step through the lists of fill-ins while marking them.
* pLastFillin (ElementPtr)
* A pointer to the last fill-in in the list. Used to terminate a loop.
* pListNode (struct FillinListNodeStruct *)
* A pointer to a node in the FillinList linked-list.
*/
void
spStripFills(MatrixPtr Matrix)
{
struct FillinListNodeStruct *pListNode;
/* Begin `spStripFills'. */
assert( IS_SPARSE( Matrix ) );
if (Matrix->Fillins == 0) return;
Matrix->NeedsOrdering = YES;
Matrix->Elements -= Matrix->Fillins;
Matrix->Fillins = 0;
/* Mark the fill-ins. */
{
ElementPtr pFillin, pLastFillin;
pListNode = Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList;
Matrix->NextAvailFillin = pListNode->pFillinList;
while (pListNode != NULL)
{
pFillin = pListNode->pFillinList;
pLastFillin = &(pFillin[ pListNode->NumberOfFillinsInList - 1 ]);
while (pFillin <= pLastFillin)
(pFillin++)->Row = 0;
pListNode = pListNode->Next;
}
}
/* Unlink fill-ins by searching for elements marked with Row = 0. */
{
ElementPtr pElement, *ppElement;
int I, Size = Matrix->Size;
/* Unlink fill-ins in all columns. */
for (I = 1; I <= Size; I++)
{
ppElement = &(Matrix->FirstInCol[I]);
while ((pElement = *ppElement) != NULL)
{
if (pElement->Row == 0)
{
*ppElement = pElement->NextInCol; /* Unlink fill-in. */
if (Matrix->Diag[pElement->Col] == pElement)
Matrix->Diag[pElement->Col] = NULL;
}
else
ppElement = &pElement->NextInCol; /* Skip element. */
}
}
/* Unlink fill-ins in all rows. */
for (I = 1; I <= Size; I++)
{
ppElement = &(Matrix->FirstInRow[I]);
while ((pElement = *ppElement) != NULL)
{
if (pElement->Row == 0)
*ppElement = pElement->NextInRow; /* Unlink fill-in. */
else
ppElement = &pElement->NextInRow; /* Skip element. */
}
}
}
return;
}
/* Same as above, but strips entire matrix without destroying the
* frame. This assumes that the matrix will be replaced with one of
* the same size. */
void
spStripMatrix(MatrixPtr Matrix)
{
/* Begin `spStripMatrix'. */
assert( IS_SPARSE( Matrix ) );
if (Matrix->Elements == 0) return;
Matrix->RowsLinked = NO;
Matrix->NeedsOrdering = YES;
Matrix->Elements = 0;
Matrix->Originals = 0;
Matrix->Fillins = 0;
/* Reset the element lists. */
{
struct ElementListNodeStruct *pListNode;
pListNode = Matrix->LastElementListNode = Matrix->FirstElementListNode;
Matrix->ElementsRemaining = pListNode->NumberOfElementsInList;
Matrix->NextAvailElement = pListNode->pElementList;
}
/* Reset the fill-in lists. */
{
struct FillinListNodeStruct *pListNode;
pListNode = Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList;
Matrix->NextAvailFillin = pListNode->pFillinList;
}
/* Reset the Row, Column and Diag pointers */
{
int I, Size = Matrix->Size;
for (I = 1; I <= Size; I++)
{
Matrix->FirstInRow[I] = NULL;
Matrix->FirstInCol[I] = NULL;
Matrix->Diag[I] = NULL;
}
}
}
#endif
#if TRANSLATE && DELETE
/*
* DELETE A ROW AND COLUMN FROM THE MATRIX
*
* Deletes a row and a column from a matrix.
*
* Sparse will abort if an attempt is made to delete a row or column that
* doesn't exist.
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix in which the row and column are to be deleted.
* Row <input> (int)
* Row to be deleted.
* Col <input> (int)
* Column to be deleted.
*
* >>> Local variables:
* ExtCol (int)
* The external column that is being deleted.
* ExtRow (int)
* The external row that is being deleted.
* pElement (ElementPtr)
* Pointer to an element in the matrix. Used when scanning rows and
* columns in order to eliminate elements from the last row or column.
* ppElement (ElementPtr *)
* Pointer to the location of an ElementPtr. This location will be
* filled with a NULL pointer if it is the new last element in its row
* or column.
* pElement (ElementPtr)
* Pointer to an element in the last row or column of the matrix.
* Size (int)
* The local version Matrix->Size, the size of the matrix.
*/
void
spDeleteRowAndCol(MatrixPtr Matrix, int Row, int Col)
{
ElementPtr pElement, *ppElement, pLastElement;
int Size, ExtRow, ExtCol;
ElementPtr spcFindElementInCol();
/* Begin `spDeleteRowAndCol'. */
assert( IS_SPARSE(Matrix) && Row > 0 && Col > 0 );
assert( Row <= Matrix->ExtSize && Col <= Matrix->ExtSize );
Size = Matrix->Size;
ExtRow = Row;
ExtCol = Col;
if (!Matrix->RowsLinked)
spcLinkRows( Matrix );
Row = Matrix->ExtToIntRowMap[Row];
Col = Matrix->ExtToIntColMap[Col];
assert( Row > 0 && Col > 0 );
/* Move Row so that it is the last row in the matrix. */
if (Row != Size) spcRowExchange( Matrix, Row, Size );
/* Move Col so that it is the last column in the matrix. */
if (Col != Size) spcColExchange( Matrix, Col, Size );
/* Correct Diag pointers. */
if (Row == Col)
SWAP( ElementPtr, Matrix->Diag[Row], Matrix->Diag[Size] );
else {
Matrix->Diag[Row] = spcFindElementInCol(Matrix,
Matrix->FirstInCol+Row,
Row, Row, NO );
Matrix->Diag[Col] = spcFindElementInCol(Matrix,
Matrix->FirstInCol+Col,
Col, Col, NO );
}
/* Delete last row and column of the matrix. */
/* Break the column links to every element in the last row. */
pLastElement = Matrix->FirstInRow[ Size ];
while (pLastElement != NULL) {
ppElement = &(Matrix->FirstInCol[ pLastElement->Col ]);
while ((pElement = *ppElement) != NULL) {
if (pElement == pLastElement)
*ppElement = NULL; /* Unlink last element in column. */
else
ppElement = &pElement->NextInCol; /* Skip element. */
}
pLastElement = pLastElement->NextInRow;
}
/* Break the row links to every element in the last column. */
pLastElement = Matrix->FirstInCol[ Size ];
while (pLastElement != NULL) {
ppElement = &(Matrix->FirstInRow[ pLastElement->Row ]);
while ((pElement = *ppElement) != NULL) {
if (pElement == pLastElement)
*ppElement = NULL; /* Unlink last element in row. */
else
ppElement = &pElement->NextInRow; /* Skip element. */
}
pLastElement = pLastElement->NextInCol;
}
/* Clean up some details. */
Matrix->Size = Size - 1;
Matrix->Diag[Size] = NULL;
Matrix->FirstInRow[Size] = NULL;
Matrix->FirstInCol[Size] = NULL;
Matrix->CurrentSize--;
Matrix->ExtToIntRowMap[ExtRow] = -1;
Matrix->ExtToIntColMap[ExtCol] = -1;
Matrix->NeedsOrdering = YES;
return;
}
#endif
#if PSEUDOCONDITION
/*
* CALCULATE PSEUDOCONDITION
*
* Computes the magnitude of the ratio of the largest to the smallest
* pivots. This quantity is an indicator of ill-conditioning in the
* matrix. If this ratio is large, and if the matrix is scaled such
* that uncertainties in the RHS and the matrix entries are
* equilibrated, then the matrix is ill-conditioned. However, a
* small ratio does not necessarily imply that the matrix is
* well-conditioned. This routine must only be used after a matrix
* has been factored by spOrderAndFactor() or spFactor() and before
* it is cleared by spClear() or spInitialize(). The pseudocondition
* is faster to compute than the condition number calculated by
* spCondition(), but is not as informative.
*
* >>> Returns:
* The magnitude of the ratio of the largest to smallest pivot used during
* previous factorization. If the matrix was singular, zero is returned.
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix. */
RealNumber
spPseudoCondition(MatrixPtr Matrix)
{
int I;
ArrayOfElementPtrs Diag;
RealNumber MaxPivot, MinPivot, Mag;
/* Begin `spPseudoCondition'. */
assert( IS_SPARSE(Matrix) && IS_FACTORED(Matrix) );
if (Matrix->Error == spSINGULAR || Matrix->Error == spZERO_DIAG)
return 0.0;
Diag = Matrix->Diag;
MaxPivot = MinPivot = ELEMENT_MAG( Diag[1] );
for (I = 2; I <= Matrix->Size; I++)
{
Mag = ELEMENT_MAG( Diag[I] );
if (Mag > MaxPivot)
MaxPivot = Mag;
else if (Mag < MinPivot)
MinPivot = Mag;
}
assert( MaxPivot > 0.0);
return MaxPivot / MinPivot;
}
#endif
#if CONDITION
/*
* ESTIMATE CONDITION NUMBER
*
* Computes an estimate of the condition number using a variation on
* the LINPACK condition number estimation algorithm. This quantity
* is an indicator of ill-conditioning in the matrix. To avoid
* problems with overflow, the reciprocal of the condition number is
* returned. If this number is small, and if the matrix is scaled
* such that uncertainties in the RHS and the matrix entries are
* equilibrated, then the matrix is ill-conditioned. If the this
* number is near one, the matrix is well conditioned. This routine
* must only be used after a matrix has been factored by
* spOrderAndFactor() or spFactor() and before it is cleared by
* spClear() or spInitialize().
*
* Unlike the LINPACK condition number estimator, this routines
* returns the L infinity condition number. This is an artifact of
* Sparse placing ones on the diagonal of the upper triangular matrix
* rather than the lower. This difference should be of no
* importance.
*
* References:
* A.K. Cline, C.B. Moler, G.W. Stewart, J.H. Wilkinson. An estimate
* for the condition number of a matrix. SIAM Journal on Numerical
* Analysis. Vol. 16, No. 2, pages 368-375, April 1979.
*
* J.J. Dongarra, C.B. Moler, J.R. Bunch, G.W. Stewart. LINPACK
* User's Guide. SIAM, 1979.
*
* Roger G. Grimes, John G. Lewis. Condition number estimation for
* sparse matrices. SIAM Journal on Scientific and Statistical
* Computing. Vol. 2, No. 4, pages 384-388, December 1981.
*
* Dianne Prost O'Leary. Estimating matrix condition numbers. SIAM
* Journal on Scientific and Statistical Computing. Vol. 1, No. 2,
* pages 205-209, June 1980.
*
* >>> Returns:
* The reciprocal of the condition number. If the matrix was singular,
* zero is returned.
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix.
* NormOfMatrix <input> (RealNumber)
* The L-infinity norm of the unfactored matrix as computed by
* spNorm().
* pError <output> (int *)
* Used to return error code.
*
* >>> Possible errors:
* spSINGULAR
* spNO_MEMORY */
RealNumber
spCondition(MatrixPtr Matrix, RealNumber NormOfMatrix, int *pError)
{
ElementPtr pElement;
RealVector T, Tm;
int I, K, Row;
ElementPtr pPivot;
int Size;
RealNumber E, Em, Wp, Wm, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor;
RealNumber Linpack, OLeary, InvNormOfInverse, ComplexCondition();
#define SLACK 1e4
/* Begin `spCondition'. */
assert( IS_SPARSE(Matrix) && IS_FACTORED(Matrix) );
*pError = Matrix->Error;
if (Matrix->Error >= spFATAL) return 0.0;
if (NormOfMatrix == 0.0)
{
*pError = spSINGULAR;
return 0.0;
}
if (Matrix->Complex)
return ComplexCondition( Matrix, NormOfMatrix, pError );
Size = Matrix->Size;
T = Matrix->Intermediate;
Tm = Matrix->Intermediate + Size;
for (I = Size; I > 0; I--) T[I] = 0.0;
/*
* Part 1. Ay = e.
*
* Solve Ay = LUy = e where e consists of +1 and -1 terms with the
* sign chosen to maximize the size of w in Lw = e. Since the
* terms in w can get very large, scaling is used to avoid
* overflow. */
/* Forward elimination. Solves Lw = e while choosing e. */
E = 1.0;
for (I = 1; I <= Size; I++)
{
pPivot = Matrix->Diag[I];
if (T[I] < 0.0) Em = -E; else Em = E;
Wm = (Em + T[I]) * pPivot->Real;
if (ABS(Wm) > SLACK)
{
ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(Wm) );
for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
E *= ScaleFactor;
Em *= ScaleFactor;
Wm = (Em + T[I]) * pPivot->Real;
}
Wp = (T[I] - Em) * pPivot->Real;
ASp = ABS(T[I] - Em);
ASm = ABS(Em + T[I]);
/* Update T for both values of W, minus value is placed in
Tm. */
pElement = pPivot->NextInCol;
while (pElement != NULL)
{
Row = pElement->Row;
Tm[Row] = T[Row] - (Wm * pElement->Real);
T[Row] -= (Wp * pElement->Real);
ASp += ABS(T[Row]);
ASm += ABS(Tm[Row]);
pElement = pElement->NextInCol;
}
/* If minus value causes more growth, overwrite T with its
values. */
if (ASm > ASp)
{
T[I] = Wm;
pElement = pPivot->NextInCol;
while (pElement != NULL)
{
T[pElement->Row] = Tm[pElement->Row];
pElement = pElement->NextInCol;
}
}
else T[I] = Wp;
}
/* Compute 1-norm of T, which now contains w, and scale ||T|| to
1/SLACK. */
for (ASw = 0.0, I = Size; I > 0; I--) ASw += ABS(T[I]);
ScaleFactor = 1.0 / (SLACK * ASw);
if (ScaleFactor < 0.5)
{
for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
E *= ScaleFactor;
}
/* Backward Substitution. Solves Uy = w.*/
for (I = Size; I >= 1; I--)
{
pElement = Matrix->Diag[I]->NextInRow;
while (pElement != NULL)
{
T[I] -= pElement->Real * T[pElement->Col];
pElement = pElement->NextInRow;
}
if (ABS(T[I]) > SLACK)
{
ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) );
for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
E *= ScaleFactor;
}
}
/* Compute 1-norm of T, which now contains y, and scale ||T|| to
1/SLACK. */
for (ASy = 0.0, I = Size; I > 0; I--) ASy += ABS(T[I]);
ScaleFactor = 1.0 / (SLACK * ASy);
if (ScaleFactor < 0.5)
{
for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
ASy = 1.0 / SLACK;
E *= ScaleFactor;
}
/* Compute infinity-norm of T for O'Leary's estimate. */
for (MaxY = 0.0, I = Size; I > 0; I--)
if (MaxY < ABS(T[I])) MaxY = ABS(T[I]);
/*
* Part 2.
*
* A* z = y where the * represents the transpose. Recall that A =
* LU implies A* = U* L*. */
/* Forward elimination, U* v = y. */
for (I = 1; I <= Size; I++)
{
pElement = Matrix->Diag[I]->NextInRow;
while (pElement != NULL)
{
T[pElement->Col] -= T[I] * pElement->Real;
pElement = pElement->NextInRow;
}
if (ABS(T[I]) > SLACK)
{
ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) );
for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
ASy *= ScaleFactor;
}
}
/* Compute 1-norm of T, which now contains v, and scale ||T|| to 1/SLACK. */
for (ASv = 0.0, I = Size; I > 0; I--) ASv += ABS(T[I]);
ScaleFactor = 1.0 / (SLACK * ASv);
if (ScaleFactor < 0.5)
{
for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
ASy *= ScaleFactor;
}
/* Backward Substitution, L* z = v. */
for (I = Size; I >= 1; I--)
{
pPivot = Matrix->Diag[I];
pElement = pPivot->NextInCol;
while (pElement != NULL)
{
T[I] -= pElement->Real * T[pElement->Row];
pElement = pElement->NextInCol;
}
T[I] *= pPivot->Real;
if (ABS(T[I]) > SLACK)
{
ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) );
for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
ASy *= ScaleFactor;
}
}
/* Compute 1-norm of T, which now contains z. */
for (ASz = 0.0, I = Size; I > 0; I--) ASz += ABS(T[I]);
Linpack = ASy / ASz;
OLeary = E / MaxY;
InvNormOfInverse = MIN( Linpack, OLeary );
return InvNormOfInverse / NormOfMatrix;
}
/*
* ESTIMATE CONDITION NUMBER
*
* Complex version of spCondition().
*
* >>> Returns:
* The reciprocal of the condition number.
*
* >>> Arguments:
* Matrix <input> (MatrixPtr)
* Pointer to the matrix.
* NormOfMatrix <input> (RealNumber)
* The L-infinity norm of the unfactored matrix as computed by
* spNorm().
* pError <output> (int *)
* Used to return error code.
*
* >>> Possible errors:
* spNO_MEMORY
*/
static RealNumber
ComplexCondition( Matrix, NormOfMatrix, pError )
MatrixPtr Matrix;
RealNumber NormOfMatrix;
int *pError;
{
ElementPtr pElement;
ComplexVector T, Tm;
int I, K, Row;
ElementPtr pPivot;
int Size;
RealNumber E, Em, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor;
RealNumber Linpack, OLeary, InvNormOfInverse;
ComplexNumber Wp, Wm;
/* Begin `ComplexCondition'. */
Size = Matrix->Size;
T = (ComplexVector)Matrix->Intermediate;
Tm = SP_MALLOC( ComplexNumber, Size+1 );
if (Tm == NULL)
{
*pError = spNO_MEMORY;
return 0.0;
}
for (I = Size; I > 0; I--) T[I].Real = T[I].Imag = 0.0;
/*
* Part 1. Ay = e.
*
* Solve Ay = LUy = e where e consists of +1 and -1 terms with the
* sign chosen to maximize the size of w in Lw = e. Since the
* terms in w can get very large, scaling is used to avoid
* overflow. */
/* Forward elimination. Solves Lw = e while choosing e. */
E = 1.0;
for (I = 1; I <= Size; I++)
{
pPivot = Matrix->Diag[I];
if (T[I].Real < 0.0) Em = -E; else Em = E;
Wm = T[I];
Wm.Real += Em;
ASm = CMPLX_1_NORM( Wm );
CMPLX_MULT_ASSIGN( Wm, *pPivot );
if (CMPLX_1_NORM(Wm) > SLACK)
{
ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(Wm) );
for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
E *= ScaleFactor;
Em *= ScaleFactor;
ASm *= ScaleFactor;
SCLR_MULT_ASSIGN( Wm, ScaleFactor );
}
Wp = T[I];
Wp.Real -= Em;
ASp = CMPLX_1_NORM( Wp );
CMPLX_MULT_ASSIGN( Wp, *pPivot );
/* Update T for both values of W, minus value is placed in Tm. */
pElement = pPivot->NextInCol;
while (pElement != NULL)
{
Row = pElement->Row;
/* Cmplx expr: Tm[Row] = T[Row] - (Wp * *pElement). */
CMPLX_MULT_SUBT( Tm[Row], Wm, *pElement, T[Row] );
/* Cmplx expr: T[Row] -= Wp * *pElement. */
CMPLX_MULT_SUBT_ASSIGN( T[Row], Wm, *pElement );
ASp += CMPLX_1_NORM(T[Row]);
ASm += CMPLX_1_NORM(Tm[Row]);
pElement = pElement->NextInCol;
}
/* If minus value causes more growth, overwrite T with its
values. */
if (ASm > ASp)
{
T[I] = Wm;
pElement = pPivot->NextInCol;
while (pElement != NULL)
{
T[pElement->Row] = Tm[pElement->Row];
pElement = pElement->NextInCol;
}
}
else T[I] = Wp;
}
/* Compute 1-norm of T, which now contains w, and scale ||T|| to
1/SLACK. */
for (ASw = 0.0, I = Size; I > 0; I--) ASw += CMPLX_1_NORM(T[I]);
ScaleFactor = 1.0 / (SLACK * ASw);
if (ScaleFactor < 0.5)
{
for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
E *= ScaleFactor;
}
/* Backward Substitution. Solves Uy = w.*/
for (I = Size; I >= 1; I--)
{
pElement = Matrix->Diag[I]->NextInRow;
while (pElement != NULL)
{
/* Cmplx expr: T[I] -= T[pElement->Col] * *pElement. */
CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Col], *pElement );
pElement = pElement->NextInRow;
}
if (CMPLX_1_NORM(T[I]) > SLACK)
{
ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
E *= ScaleFactor;
}
}
/* Compute 1-norm of T, which now contains y, and scale ||T|| to
1/SLACK. */
for (ASy = 0.0, I = Size; I > 0; I--) ASy += CMPLX_1_NORM(T[I]);
ScaleFactor = 1.0 / (SLACK * ASy);
if (ScaleFactor < 0.5)
{
for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
ASy = 1.0 / SLACK;
E *= ScaleFactor;
}
/* Compute infinity-norm of T for O'Leary's estimate. */
for (MaxY = 0.0, I = Size; I > 0; I--)
if (MaxY < CMPLX_1_NORM(T[I])) MaxY = CMPLX_1_NORM(T[I]);
/* Part 2. A* z = y where the * represents the transpose. Recall
* that A = LU implies A* = U* L*. */
/* Forward elimination, U* v = y. */
for (I = 1; I <= Size; I++)
{
pElement = Matrix->Diag[I]->NextInRow;
while (pElement != NULL)
{
/* Cmplx expr: T[pElement->Col] -= T[I] * *pElement. */
CMPLX_MULT_SUBT_ASSIGN( T[pElement->Col], T[I], *pElement );
pElement = pElement->NextInRow;
}
if (CMPLX_1_NORM(T[I]) > SLACK)
{
ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
ASy *= ScaleFactor;
}
}
/* Compute 1-norm of T, which now contains v, and scale ||T|| to
1/SLACK. */
for (ASv = 0.0, I = Size; I > 0; I--) ASv += CMPLX_1_NORM(T[I]);
ScaleFactor = 1.0 / (SLACK * ASv);
if (ScaleFactor < 0.5)
{
for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
ASy *= ScaleFactor;
}
/* Backward Substitution, L* z = v. */
for (I = Size; I >= 1; I--)
{
pPivot = Matrix->Diag[I];
pElement = pPivot->NextInCol;
while (pElement != NULL)
{
/* Cmplx expr: T[I] -= T[pElement->Row] * *pElement. */
CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Row], *pElement );
pElement = pElement->NextInCol;
}
CMPLX_MULT_ASSIGN( T[I], *pPivot );
if (CMPLX_1_NORM(T[I]) > SLACK)
{
ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
ASy *= ScaleFactor;
}
}
/* Compute 1-norm of T, which now contains z. */
for (ASz = 0.0, I = Size; I > 0; I--) ASz += CMPLX_1_NORM(T[I]);
SP_FREE( Tm );
Linpack = ASy / ASz;
OLeary = E / MaxY;
InvNormOfInverse = MIN( Linpack, OLeary );
return InvNormOfInverse / NormOfMatrix;
}
/*
* L-INFINITY MATRIX NORM
*
* Computes the L-infinity norm of an unfactored matrix. It is a fatal
* error to pass this routine a factored matrix.
*
* One difficulty is that the rows may not be linked.
*
* >>> Returns:
* The largest absolute row sum of matrix.
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix.
*/
RealNumber
spNorm(MatrixPtr Matrix)
{
ElementPtr pElement;
int I;
RealNumber Max = 0.0, AbsRowSum;
/* Begin `spNorm'. */
assert( IS_SPARSE(Matrix) && !IS_FACTORED(Matrix) );
if (!Matrix->RowsLinked) spcLinkRows( Matrix );
/* Compute row sums. */
if (!Matrix->Complex)
{
for (I = Matrix->Size; I > 0; I--)
{
pElement = Matrix->FirstInRow[I];
AbsRowSum = 0.0;
while (pElement != NULL)
{
AbsRowSum += ABS( pElement->Real );
pElement = pElement->NextInRow;
}
if (Max < AbsRowSum) Max = AbsRowSum;
}
}
else
{
for (I = Matrix->Size; I > 0; I--)
{
pElement = Matrix->FirstInRow[I];
AbsRowSum = 0.0;
while (pElement != NULL)
{
AbsRowSum += CMPLX_1_NORM( *pElement );
pElement = pElement->NextInRow;
}
if (Max < AbsRowSum) Max = AbsRowSum;
}
}
return Max;
}
#endif /* CONDITION */
#if STABILITY
/*
* STABILITY OF FACTORIZATION
*
* The following routines are used to gauge the stability of a
* factorization. If the factorization is determined to be too
* unstable, then the matrix should be reordered. The routines
* compute quantities that are needed in the computation of a bound
* on the error attributed to any one element in the matrix during
* the factorization. In other words, there is a matrix E = [e_ij]
* of error terms such that A+E = LU. This routine finds a bound on
* |e_ij|. Erisman & Reid [1] showed that |e_ij| < 3.01 u rho m_ij,
* where u is the machine rounding unit, rho = max a_ij where the max
* is taken over every row i, column j, and step k, and m_ij is the
* number of multiplications required in the computation of l_ij if i
* > j or u_ij otherwise. Barlow [2] showed that rho < max_i || l_i
* ||_p max_j || u_j ||_q where 1/p + 1/q = 1.
*
* The first routine finds the magnitude on the largest element in
* the matrix. If the matrix has not yet been factored, the largest
* element is found by direct search. If the matrix is factored, a
* bound on the largest element in any of the reduced submatrices is
* computed using Barlow with p = oo and q = 1. The ratio of these
* two numbers is the growth, which can be used to determine if the
* pivoting order is adequate. A large growth implies that
* considerable error has been made in the factorization and that it
* is probably a good idea to reorder the matrix. If a large growth
* in encountered after using spFactor(), reconstruct the matrix and
* refactor using spOrderAndFactor(). If a large growth is
* encountered after using spOrderAndFactor(), refactor using
* spOrderAndFactor() with the pivot threshold increased, say to 0.1.
*
* Using only the size of the matrix as an upper bound on m_ij and
* Barlow's bound, the user can estimate the size of the matrix error
* terms e_ij using the bound of Erisman and Reid. The second
* routine computes a tighter bound (with more work) based on work by
* Gear [3], |e_ij| < 1.01 u rho (t c^3 + (1 + t)c^2) where t is the
* threshold and c is the maximum number of off-diagonal elements in
* any row of L. The expensive part of computing this bound is
* determining the maximum number of off-diagonals in L, which
* changes only when the order of the matrix changes. This number is
* computed and saved, and only recomputed if the matrix is
* reordered.
*
* [1] A. M. Erisman, J. K. Reid. Monitoring the stability of the
* triangular factorization of a sparse matrix. Numerische
* Mathematik. Vol. 22, No. 3, 1974, pp 183-186.
*
* [2] J. L. Barlow. A note on monitoring the stability of triangular
* decomposition of sparse matrices. "SIAM Journal of Scientific
* and Statistical Computing." Vol. 7, No. 1, January 1986, pp 166-168.
*
* [3] I. S. Duff, A. M. Erisman, J. K. Reid. "Direct Methods for Sparse
* Matrices." Oxford 1986. pp 99. */
/*
* LARGEST ELEMENT IN MATRIX
*
* >>> Returns:
* If matrix is not factored, returns the magnitude of the largest element in
* the matrix. If the matrix is factored, a bound on the magnitude of the
* largest element in any of the reduced submatrices is returned.
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix.
*/
RealNumber
spLargestElement(MatrixPtr Matrix)
{
int I;
RealNumber Mag, AbsColSum, Max = 0.0, MaxRow = 0.0, MaxCol = 0.0;
RealNumber Pivot;
ComplexNumber cPivot;
ElementPtr pElement, pDiag;
/* Begin `spLargestElement'. */
assert( IS_SPARSE(Matrix) );
if (Matrix->Factored && !Matrix->Complex)
{
if (Matrix->Error == spSINGULAR) return 0.0;
/* Find the bound on the size of the largest element over all
factorization. */
for (I = 1; I <= Matrix->Size; I++)
{
pDiag = Matrix->Diag[I];
/* Lower triangular matrix. */
Pivot = 1.0 / pDiag->Real;
Mag = ABS( Pivot );
if (Mag > MaxRow) MaxRow = Mag;
pElement = Matrix->FirstInRow[I];
while (pElement != pDiag)
{
Mag = ABS( pElement->Real );
if (Mag > MaxRow) MaxRow = Mag;
pElement = pElement->NextInRow;
}
/* Upper triangular matrix. */
pElement = Matrix->FirstInCol[I];
AbsColSum = 1.0; /* Diagonal of U is unity. */
while (pElement != pDiag)
{
AbsColSum += ABS( pElement->Real );
pElement = pElement->NextInCol;
}
if (AbsColSum > MaxCol) MaxCol = AbsColSum;
}
}
else if (!Matrix->Complex)
{
for (I = 1; I <= Matrix->Size; I++)
{
pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{
Mag = ABS( pElement->Real );
if (Mag > Max) Max = Mag;
pElement = pElement->NextInCol;
}
}
return Max;
}
if (Matrix->Factored && Matrix->Complex)
{
if (Matrix->Error == spSINGULAR) return 0.0;
/* Find the bound on the size of the largest element over all
factorization. */
for (I = 1; I <= Matrix->Size; I++)
{
pDiag = Matrix->Diag[I];
/* Lower triangular matrix. */
CMPLX_RECIPROCAL( cPivot, *pDiag );
Mag = CMPLX_1_NORM( cPivot );
if (Mag > MaxRow) MaxRow = Mag;
pElement = Matrix->FirstInRow[I];
while (pElement != pDiag)
{
Mag = CMPLX_1_NORM( *pElement );
if (Mag > MaxRow) MaxRow = Mag;
pElement = pElement->NextInRow;
}
/* Upper triangular matrix. */
pElement = Matrix->FirstInCol[I];
AbsColSum = 1.0; /* Diagonal of U is unity. */
while (pElement != pDiag)
{
AbsColSum += CMPLX_1_NORM( *pElement );
pElement = pElement->NextInCol;
}
if (AbsColSum > MaxCol) MaxCol = AbsColSum;
}
}
else if (Matrix->Complex)
{
for (I = 1; I <= Matrix->Size; I++)
{
pElement = Matrix->FirstInCol[I];
while (pElement != NULL)
{
Mag = CMPLX_1_NORM( *pElement );
if (Mag > Max) Max = Mag;
pElement = pElement->NextInCol;
}
}
return Max;
}
return MaxRow * MaxCol;
}
/*
* MATRIX ROUNDOFF ERROR
*
* >>> Returns:
* Returns a bound on the magnitude of the largest element in E = A - LU.
*
* >>> Arguments:
* Matrix <input> (char *)
* Pointer to the matrix.
* Rho <input> (RealNumber)
* The bound on the magnitude of the largest element in any of the
* reduced submatrices. This is the number computed by the function
* spLargestElement() when given a factored matrix. If this number is
* negative, the bound will be computed automatically.
*/
RealNumber
spRoundoff(MatrixPtr Matrix, RealNumber Rho)
{
ElementPtr pElement;
int Count, I, MaxCount = 0;
RealNumber Reid, Gear;
/* Begin `spRoundoff'. */
assert( IS_SPARSE(Matrix) && IS_FACTORED(Matrix) );
/* Compute Barlow's bound if it is not given. */
if (Rho < 0.0) Rho = spLargestElement( Matrix );
/* Find the maximum number of off-diagonals in L if not previously computed. */
if (Matrix->MaxRowCountInLowerTri < 0)
{
for (I = Matrix->Size; I > 0; I--)
{
pElement = Matrix->FirstInRow[I];
Count = 0;
while (pElement->Col < I)
{
Count++;
pElement = pElement->NextInRow;
}
if (Count > MaxCount) MaxCount = Count;
}
Matrix->MaxRowCountInLowerTri = MaxCount;
}
else MaxCount = Matrix->MaxRowCountInLowerTri;
/* Compute error bound. */
Gear = 1.01*((MaxCount + 1) * Matrix->RelThreshold + 1.0) * SQR(MaxCount);
Reid = 3.01 * Matrix->Size;
if (Gear < Reid)
return (MACHINE_RESOLUTION * Rho * Gear);
else
return (MACHINE_RESOLUTION * Rho * Reid);
}
#endif
#if DOCUMENTATION
/*
* SPARSE ERROR MESSAGE
*
* This routine prints a short message to a stream describing the error
* error state of sparse. No message is produced if there is no error.
*
* >>> Arguments:
* Matrix <input> (char *)
* Matrix for which the error message is to be printed.
* Stream <input> (FILE *)
* Stream to which the error message is to be printed.
* Originator <input> (char *)
* Name of originator of error message. If NULL, `sparse' is used.
* If zero-length string, no originator is printed.
*/
void
spErrorMessage(MatrixPtr Matrix, FILE *Stream, char *Originator)
{
int Row, Col, Error;
/* Begin `spErrorMessage'. */
if (Matrix == NULL)
Error = spNO_MEMORY;
else
{
assert(Matrix->ID == SPARSE_ID);
Error = Matrix->Error;
}
if (Error == spOKAY) return;
if (Originator == NULL) Originator = "sparse";
if (Originator[0] != '\0') fprintf( Stream, "%s: ", Originator);
if (Error >= spFATAL)
fprintf( Stream, "fatal error, ");
else
fprintf( Stream, "warning, ");
/* Print particular error message. Do not use switch statement
* because error codes may not be unique. */
if (Error == spPANIC)
fprintf( Stream, "Sparse called improperly.\n");
else if (Error == spNO_MEMORY)
fprintf( Stream, "insufficient memory available.\n");
else if (Error == spSINGULAR)
{
spWhereSingular( Matrix, &Row, &Col );
fprintf( Stream, "singular matrix detected at row %d and column %d.\n",
Row, Col);
}
else if (Error == spZERO_DIAG)
{
spWhereSingular( Matrix, &Row, &Col );
fprintf( Stream, "zero diagonal detected at row %d and column %d.\n",
Row, Col);
}
else if (Error == spSMALL_PIVOT)
{
fprintf( Stream,
"unable to find a pivot that is larger than absolute threshold.\n");
}
else
{
abort();
}
return;
}
#endif /* DOCUMENTATION */