/*
 *  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
 *  spStrip
 *  spDeleteRowAndCol
 *  spPseudoCondition
 *  spCondition
 *  spNorm
 *  spLargestElement
 *  spRoundoff
 *
 *  >>> Other functions contained in this file:
 *  CountTwins
 *  SwapCols
 *  ScaleComplexMatrix
 *  ComplexMatrixMultiply
 *  ComplexCondition
 */


/*
 *  Revision and copyright information.
 *
 *  Copyright (c) 1985,86,87,88
 *  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.
 */

#define spINSIDE_SPARSE


#include "spConfig.h"
#include "spUtils.h"
#include "spmatrix.h"
#include "spDefs.h"
#include "spmalloc.h"
#include "spBuild.h"

static int CountTwins(MatrixPtr Matrix, int Col, ElementPtr *ppTwin1, ElementPtr *ppTwin2 );


static int SwapCols( MatrixPtr Matrix, ElementPtr pTwin1, ElementPtr pTwin2 );

#if spCOMPLEX AND SCALING
static void ScaleComplexMatrix( MatrixPtr Matrix, register  RealVector  RHS_ScaleFactors, register  RealVector  SolutionScaleFactors );
#endif

#if spCOMPLEX AND MULTIPLICATION
static void ComplexMatrixMultiply( MatrixPtr Matrix, RealVector RHS , RealVector Solution IMAG_VECTORS );
#endif

#if spCOMPLEX AND MULTIPLICATION AND TRANSPOSE
static void
ComplexTransposedMatrixMultiply( MatrixPtr Matrix, RealVector RHS, RealVector Solution IMAG_VECTORS );
#endif

#if spCOMPLEX
static RealNumber ComplexCondition( MatrixPtr Matrix, RealNumber NormOfMatrix, int *pError );
#endif

#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:
 *  eMatrix  <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  (SPBOOLEAN)
 *      Flag indicating that at least one zero diagonal with symmetric twins
 *      remain.
 *  StartAt  (int)
 *      Column number of first zero diagonal with symmetric twins.
 *  Swapped  (SPBOOLEAN)
 *      Flag indicating that columns were swapped on this pass.
 *  Twins  (int)
 *      Number of symmetric twins corresponding to current zero diagonal.
 */
#undef SPBOOLEAN
#define  SPBOOLEAN        int

void spMNA_Preorder( char * eMatrix )
{
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    register  int  J, Size;
    ElementPtr  pTwin1, pTwin2;
    int  Twins, StartAt = 1;
    SPBOOLEAN  Swapped, AnotherPassNeeded;

    /* Begin `spMNA_Preorder'. */
    ASSERT( IS_VALID(Matrix) AND NOT 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) AND NOT AnotherPassNeeded)
                {
                    AnotherPassNeeded = YES;
                    StartAt = J;
                }
            }
        }

        /* All lone twins are gone, look for zero diagonals with multiple twins. */
        if (AnotherPassNeeded)
        {
            for (J = StartAt; NOT Swapped AND (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) AND (pTwin2->Row != Col))
            {
                pTwin2 = pTwin2->NextInCol;
            }
            if ((pTwin2 != NULL) AND (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 int
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 = NOT Matrix->NumberOfInterchangesIsOdd;
    return 0;
}
#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 int 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:
 *  eMatrix  <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(char *eMatrix, register  RealVector  RHS_ScaleFactors, register RealVector SolutionScaleFactors )
{
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    register ElementPtr  pElement;
    register int  I, lSize, *pExtOrder;
    RealNumber  ScaleFactor;


    /* Begin `spScale'. */
    ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored );
    if (NOT Matrix->RowsLinked)
    {
        spcLinkRows( Matrix );
    }

#if spCOMPLEX
    if (Matrix->Complex)
    {
        ScaleComplexMatrix( Matrix, RHS_ScaleFactors, SolutionScaleFactors );
        return;
    }
#endif

#if REAL
    lSize = Matrix->Size;

    /* Correct pointers to arrays for ARRAY_OFFSET */
#if NOT ARRAY_OFFSET
    --RHS_ScaleFactors;
    --SolutionScaleFactors;
#endif

    /* 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 /* REAL */
}
#endif /* SCALING */

#if spCOMPLEX AND 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 int 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( MatrixPtr Matrix, register  RealVector  RHS_ScaleFactors, register  RealVector  SolutionScaleFactors )
{
    register ElementPtr  pElement;
    register int  I, lSize, *pExtOrder;
    RealNumber  ScaleFactor;

    /* Begin `ScaleComplexMatrix'. */
    lSize = Matrix->Size;

    /* Correct pointers to arrays for ARRAY_OFFSET */
#if NOT ARRAY_OFFSET
    --RHS_ScaleFactors;
    --SolutionScaleFactors;
#endif

    /* 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 AND spCOMPLEX */

#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:
 *  eMatrix  <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.  This is only necessary if the matrix is
 *      complex and spSEPARATED_COMPLEX_VECTORS is true.
 *  iSolution  <input>  (RealVector)
 *      iSolution is the imaginary portion of the vector being multiplied
 *      by the matrix. This is only necessary if the matrix is
 *      complex and spSEPARATED_COMPLEX_VECTORS is true.
 *
 *  >>> Obscure Macros
 *  IMAG_VECTORS
 *      Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
 *      spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
 *      without a trace.
 */

void spMultiply(char *eMatrix, RealVector RHS, RealVector Solution IMAG_VECTORS )
{
    register  ElementPtr  pElement;
    register  RealVector  Vector;
    register  RealNumber  Sum;
    register  int  I, *pExtOrder;
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;

    /* Begin `spMultiply'. */
    ASSERT( IS_SPARSE( Matrix ) AND NOT Matrix->Factored );
    if (NOT Matrix->RowsLinked)
    {
        spcLinkRows(Matrix);
    }

#if spCOMPLEX
    if (Matrix->Complex)
    {
        ComplexMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS );
        return;
    }
#endif

#if REAL
#if NOT ARRAY_OFFSET
    /* Correct array pointers for ARRAY_OFFSET. */
    --RHS;
    --Solution;
#endif

    /* 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 /* REAL */
}
#endif /* MULTIPLICATION */

#if spCOMPLEX AND 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.
 *      This is only the real portion of the right-hand side if the matrix
 *      is complex and spSEPARATED_COMPLEX_VECTORS is set true.
 *  Solution  <input>  (RealVector)
 *      Solution is the vector being multiplied by the matrix. This is only
 *      the real portion if the matrix is complex and
 *      spSEPARATED_COMPLEX_VECTORS is set true.
 *  iRHS  <output>  (RealVector)
 *      iRHS is the imaginary portion of the right hand side. This is
 *      what is being solved for.  This is only necessary if the matrix is
 *      complex and spSEPARATED_COMPLEX_VECTORS is true.
 *  iSolution  <input>  (RealVector)
 *      iSolution is the imaginary portion of the vector being multiplied
 *      by the matrix. This is only necessary if the matrix is
 *      complex and spSEPARATED_COMPLEX_VECTORS is true.
 *
 *  >>> Obscure Macros
 *  IMAG_VECTORS
 *      Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
 *      spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
 *      without a trace.
 */

static void
ComplexMatrixMultiply( MatrixPtr Matrix, RealVector RHS , RealVector Solution IMAG_VECTORS )
{
    register  ElementPtr  pElement;
    register  ComplexVector  Vector;
    register  ComplexNumber  Sum;
    register  int  I, *pExtOrder;

    /* Begin `ComplexMatrixMultiply'. */

    /* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
#if spSEPARATED_COMPLEX_VECTORS
    --RHS;
    --iRHS;
    --Solution;
    --iSolution;
#else
    RHS -= 2;
    Solution -= 2;
#endif
#endif

    /* Initialize Intermediate vector with reordered Solution vector. */
    Vector = (ComplexVector)Matrix->Intermediate;
    pExtOrder = &Matrix->IntToExtColMap[Matrix->Size];

#if spSEPARATED_COMPLEX_VECTORS
    for (I = Matrix->Size; I > 0; I--)
    {
        Vector[I].Real = Solution[*pExtOrder];
        Vector[I].Imag = iSolution[*(pExtOrder--)];
    }
#else
    for (I = Matrix->Size; I > 0; I--)
    {
        Vector[I] = ((ComplexVector)Solution)[*(pExtOrder--)];
    }
#endif

    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;
        }

#if spSEPARATED_COMPLEX_VECTORS
        RHS[*pExtOrder] = Sum.Real;
        iRHS[*pExtOrder--] = Sum.Imag;
#else
        ((ComplexVector)RHS)[*pExtOrder--] = Sum;
#endif
    }
    return;
}
#endif /* spCOMPLEX AND MULTIPLICATION */

#if MULTIPLICATION AND 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:
 *  eMatrix  <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.  This is only necessary if the matrix is
 *      complex and spSEPARATED_COMPLEX_VECTORS is true.
 *  iSolution  <input>  (RealVector)
 *      iSolution is the imaginary portion of the vector being multiplied
 *      by the matrix. This is only necessary if the matrix is
 *      complex and spSEPARATED_COMPLEX_VECTORS is true.
 *
 *  >>> Obscure Macros
 *  IMAG_VECTORS
 *      Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
 *      spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
 *      without a trace.
 */

void spMultTransposed(char *eMatrix, RealVector RHS, RealVector Solution IMAG_VECTORS )
{
    register  ElementPtr  pElement;
    register  RealVector  Vector;
    register  RealNumber  Sum;
    register  int  I, *pExtOrder;
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;


    /* Begin `spMultTransposed'. */
    ASSERT( IS_SPARSE( Matrix ) AND NOT Matrix->Factored );

#if spCOMPLEX
    if (Matrix->Complex)
    {
        ComplexTransposedMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS );
        return;
    }
#endif

#if REAL
#if NOT ARRAY_OFFSET
    /* Correct array pointers for ARRAY_OFFSET. */
    --RHS;
    --Solution;
#endif

    /* 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 /* REAL */
}
#endif /* MULTIPLICATION AND TRANSPOSE */

#if spCOMPLEX AND MULTIPLICATION AND 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.
 *      This is only the real portion of the right-hand side if the matrix
 *      is complex and spSEPARATED_COMPLEX_VECTORS is set true.
 *  Solution  <input>  (RealVector)
 *      Solution is the vector being multiplied by the matrix. This is only
 *      the real portion if the matrix is complex and
 *      spSEPARATED_COMPLEX_VECTORS is set true.
 *  iRHS  <output>  (RealVector)
 *      iRHS is the imaginary portion of the right hand side. This is
 *      what is being solved for.  This is only necessary if the matrix is
 *      complex and spSEPARATED_COMPLEX_VECTORS is true.
 *  iSolution  <input>  (RealVector)
 *      iSolution is the imaginary portion of the vector being multiplied
 *      by the matrix. This is only necessary if the matrix is
 *      complex and spSEPARATED_COMPLEX_VECTORS is true.
 *
 *  >>> Obscure Macros
 *  IMAG_VECTORS
 *      Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
 *      spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
 *      without a trace.
 */

static void
ComplexTransposedMatrixMultiply( MatrixPtr Matrix, RealVector RHS, RealVector Solution IMAG_VECTORS )
{
    register  ElementPtr  pElement;
    register  ComplexVector  Vector;
    register  ComplexNumber  Sum;
    register  int  I, *pExtOrder;

    /* Begin `ComplexMatrixMultiply'. */

    /* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
#if spSEPARATED_COMPLEX_VECTORS
    --RHS;
    --iRHS;
    --Solution;
    --iSolution;
#else
    RHS -= 2;
    Solution -= 2;
#endif
#endif

    /* Initialize Intermediate vector with reordered Solution vector. */
    Vector = (ComplexVector)Matrix->Intermediate;
    pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size];

#if spSEPARATED_COMPLEX_VECTORS
    for (I = Matrix->Size; I > 0; I--)
    {
        Vector[I].Real = Solution[*pExtOrder];
        Vector[I].Imag = iSolution[*(pExtOrder--)];
    }
#else
    for (I = Matrix->Size; I > 0; I--)
    {
        Vector[I] = ((ComplexVector)Solution)[*(pExtOrder--)];
    }
#endif

    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;
        }

#if spSEPARATED_COMPLEX_VECTORS
        RHS[*pExtOrder] = Sum.Real;
        iRHS[*pExtOrder--] = Sum.Imag;
#else
        ((ComplexVector)RHS)[*pExtOrder--] = Sum;
#endif
    }
    return;
}
#endif /* spCOMPLEX AND MULTIPLICATION AND 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:
 *  eMatrix  <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 register for speed.
 *  Temp  (RealNumber)
 *      Temporary storage for real portion of determinant.
 */

#if spCOMPLEX
void
spDeterminant( char *eMatrix, int *pExponent, register  RealNumber *pDeterminant, RealNumber *piDeterminant )
#else
void
spDeterminant( char *eMatrix, int *pExponent, register  RealNumber *pDeterminant )
#endif
{
    register MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    register 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 ) AND IS_FACTORED(Matrix) );
    *pExponent = 0;

    if (Matrix->Error == spSINGULAR)
    {
        *pDeterminant = 0.0;
#if spCOMPLEX
        if (Matrix->Complex)
        {
            *piDeterminant = 0.0;
        }
#endif
        return;
    }

    Size = Matrix->Size;
    I = 0;

#if spCOMPLEX
    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;
    }
#endif /* spCOMPLEX */
#if REAL AND spCOMPLEX
    else
#endif
#if REAL
    {
        /* 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 /* REAL */
}
#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(char *eMatrix)
{
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    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. */
    {
        register  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. */
    {
        register  ElementPtr pElement, *ppElement;
        register  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;
}
#endif

#if TRANSLATE AND 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:
 *  eMatrix  <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.
 */

extern void spcRowExchange(MatrixPtr Matrix, int  Row1, int Row2);
extern void spcColExchange(MatrixPtr Matrix, int  Col1, int Col2);

void spDeleteRowAndCol( char *eMatrix, int Row, int Col )
{
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    register  ElementPtr  pElement, *ppElement, pLastElement;
    int  Size, ExtRow, ExtCol;


    /* Begin `spDeleteRowAndCol'. */

    ASSERT( IS_SPARSE(Matrix) AND Row > 0 AND Col > 0 );
    ASSERT( Row <= Matrix->ExtSize AND Col <= Matrix->ExtSize );

    Size = Matrix->Size;
    ExtRow = Row;
    ExtCol = Col;
    if (NOT Matrix->RowsLinked)
    {
        spcLinkRows( Matrix );
    }

    Row = Matrix->ExtToIntRowMap[Row];
    Col = Matrix->ExtToIntColMap[Col];
    ASSERT( Row > 0 AND 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:
 *  eMatrix  <input>  (char *)
 *      Pointer to the matrix.
 */

RealNumber spPseudoCondition(char *eMatrix)
{
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    register int I;
    register ArrayOfElementPtrs Diag;
    RealNumber MaxPivot, MinPivot, Mag;

    /* Begin `spPseudoCondition'. */

    ASSERT( IS_SPARSE(Matrix) AND IS_FACTORED(Matrix) );
    if (Matrix->Error == spSINGULAR OR 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:
 *  eMatrix  <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(char *eMatrix, RealNumber NormOfMatrix, int *pError)
{
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    register ElementPtr pElement;
    register RealVector T, Tm;
    register 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;
#define SLACK   1e4

    /* Begin `spCondition'. */

    ASSERT( IS_SPARSE(Matrix) AND IS_FACTORED(Matrix) );
    *pError = Matrix->Error;
    if (Matrix->Error >= spFATAL)
    {
        return 0.0;
    }
    if (NormOfMatrix == 0.0)
    {
        *pError = spSINGULAR;
        return 0.0;
    }

#if spCOMPLEX
    if (Matrix->Complex)
    {
        return ComplexCondition( Matrix, NormOfMatrix, pError );
    }
#endif

#if REAL
    Size = Matrix->Size;
    T = Matrix->Intermediate;
#if spCOMPLEX
    Tm = Matrix->Intermediate + Size;
#else
    Tm = SPALLOC( RealNumber, Size + 1 );
    if (Tm == NULL)
    {
        *pError = spNO_MEMORY;
        return 0.0;
    }
#endif
    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]);
    }

#if NOT spCOMPLEX
    SPFREE( Tm );
#endif

    Linpack = ASy / ASz;
    OLeary = E / MaxY;
    InvNormOfInverse = Min( Linpack, OLeary );
    return InvNormOfInverse / NormOfMatrix;
#endif /* REAL */
}

#if spCOMPLEX
/*
 *  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( MatrixPtr Matrix, RealNumber NormOfMatrix, int *pError )
{
    register ElementPtr pElement;
    register ComplexVector T, Tm;
    register 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 = SPALLOC( 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]);
    }

    SPFREE( Tm );

    Linpack = ASy / ASz;
    OLeary = E / MaxY;
    InvNormOfInverse = Min( Linpack, OLeary );
    return InvNormOfInverse / NormOfMatrix;
}
#endif /* spCOMPLEX */

/*
 *  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:
 *  eMatrix  <input>  (char *)
 *      Pointer to the matrix.
 */

RealNumber spNorm( char *eMatrix)
{
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    register ElementPtr pElement;
    register int I;
    RealNumber Max = 0.0, AbsRowSum;

    /* Begin `spNorm'. */
    ASSERT( IS_SPARSE(Matrix) AND NOT IS_FACTORED(Matrix) );
    if (NOT Matrix->RowsLinked)
    {
        spcLinkRows( Matrix );
    }

    /* Compute row sums. */
#if REAL
    if (NOT 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;
            }
        }
    }
#endif
#if spCOMPLEX
    if (Matrix->Complex)
    {
        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;
            }
        }
    }
#endif
    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:
 *  eMatrix  <input>  (char *)
 *      Pointer to the matrix.
 */

RealNumber spLargestElement( char *eMatrix )
{
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    register int I;
    RealNumber Mag, AbsColSum, Max = 0.0, MaxRow = 0.0, MaxCol = 0.0;
    RealNumber Pivot;
    ComplexNumber cPivot;
    register ElementPtr pElement, pDiag;

    /* Begin `spLargestElement'. */
    ASSERT( IS_SPARSE(Matrix) );

#if REAL
    if (Matrix->Factored AND NOT 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 (NOT 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;
    }
#endif
#if spCOMPLEX
    if (Matrix->Factored AND 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;
    }
#endif
    return MaxRow * MaxCol;
}




/*
 *  MATRIX ROUNDOFF ERROR
 *
 *  >>> Returns:
 *  Returns a bound on the magnitude of the largest element in E = A - LU.
 *
 *  >>> Arguments:
 *  eMatrix  <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(char *eMatrix, RealNumber Rho )
{
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    register ElementPtr pElement;
    register int Count, I, MaxCount = 0;
    RealNumber Reid, Gear;

    /* Begin `spRoundoff'. */
    ASSERT( IS_SPARSE(Matrix) AND IS_FACTORED(Matrix) );

    /* Compute Barlow's bound if it is not given. */
    if (Rho < 0.0)
    {
        Rho = spLargestElement( eMatrix );
    }

    /* 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