sputils.c

spice中支持多层次元件模型仿真的可单独运行的插件源码

            }
        }
    }
    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( eMatrix )

char *eMatrix;
{
MatrixPtr  Matrix = (MatrixPtr)eMatrix;

/* 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. */
    {   register  ElementPtr  pElement;
	struct ElementListNodeStruct  *pListNode;

        pListNode = Matrix->LastElementListNode = Matrix->FirstElementListNode;
        Matrix->ElementsRemaining = pListNode->NumberOfElementsInList;
        Matrix->NextAvailElement = pListNode->pElementList;
    }

/* Reset the fill-in lists. */
    {   register  ElementPtr  pFillin;
	struct FillinListNodeStruct  *pListNode;

        pListNode = Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
        Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList;
        Matrix->NextAvailFillin = pListNode->pFillinList;
    }

/* Reset the Row, Column and Diag pointers */
    {   register  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 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.
 */

void
spDeleteRowAndCol( eMatrix, Row, Col )

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

/* 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( eMatrix )

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( eMatrix, NormOfMatrix, pError )

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, ComplexCondition();
#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 = ALLOC( 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]) );