sputils.c

ngspice又一个电子CAD仿真软件代码.功能更全


#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:
 *  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.
 *  iSolution  <input>  (RealVector)
 *      iSolution is the imaginary portion of the vector being multiplied
 *      by the matrix.
 */

void
spMultiply(void *eMatrix, RealVector RHS, RealVector Solution,
	   RealVector iRHS, RealVector iSolution)
{
    ElementPtr  pElement;
    RealVector  Vector;
    RealNumber  Sum;
    int  I, *pExtOrder;
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    extern void ComplexMatrixMultiply();

    /* 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( Matrix, RHS, Solution , iRHS, iSolution )

MatrixPtr  Matrix;
RealVector RHS, Solution , iRHS, 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:
 *  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.
 *  iSolution  <input>  (RealVector)
 *      iSolution is the imaginary portion of the vector being multiplied
 *      by the matrix.
 */

void
spMultTransposed(void *eMatrix, RealVector RHS, RealVector Solution,
		 RealVector iRHS, RealVector iSolution)
{
    ElementPtr  pElement;
    RealVector  Vector;
    RealNumber  Sum;
    int  I, *pExtOrder;
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    extern void ComplexTransposedMatrixMultiply();

    /* 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( Matrix, RHS, Solution , iRHS, iSolution )

MatrixPtr  Matrix;
RealVector RHS, Solution , iRHS, 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