sputils.c

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

 *  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 for speed.
 *  Temp  (RealNumber)
 *      Temporary storage for real portion of determinant.
 */

void
spDeterminant(void *eMatrix, int *pExponent, RealNumber *pDeterminant,
	      RealNumber *piDeterminant)
{
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    int I, Size;
    RealNumber Norm, nr, ni;
    ComplexNumber Pivot, cDeterminant;

#define  NORM(a)     (nr = ABS((a).Real), ni = ABS((a).Imag), MAX (nr,ni))

    /* Begin `spDeterminant'. */
    assert( IS_SPARSE( Matrix ) && IS_FACTORED(Matrix) );
    *pExponent = 0;

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

    Size = Matrix->Size;
    I = 0;

    if (Matrix->Complex)        /* Complex Case. */
    {
	cDeterminant.Real = 1.0;
        cDeterminant.Imag = 0.0;

        while (++I <= Size)
        {
	    CMPLX_RECIPROCAL( Pivot, *Matrix->Diag[I] );
            CMPLX_MULT_ASSIGN( cDeterminant, Pivot );

	    /* Scale Determinant. */
            Norm = NORM( cDeterminant );
            if (Norm != 0.0)
            {
		while (Norm >= 1.0e12)
                {
		    cDeterminant.Real *= 1.0e-12;
                    cDeterminant.Imag *= 1.0e-12;
                    *pExponent += 12;
                    Norm = NORM( cDeterminant );
                }
                while (Norm < 1.0e-12)
                {
		    cDeterminant.Real *= 1.0e12;
                    cDeterminant.Imag *= 1.0e12;
                    *pExponent -= 12;
                    Norm = NORM( cDeterminant );
                }
            }
        }

	/* Scale Determinant again, this time to be between 1.0 <= x < 10.0. */
        Norm = NORM( cDeterminant );
        if (Norm != 0.0)
        {
	    while (Norm >= 10.0)
            {
		cDeterminant.Real *= 0.1;
                cDeterminant.Imag *= 0.1;
                (*pExponent)++;
                Norm = NORM( cDeterminant );
            }
            while (Norm < 1.0)
            {
		cDeterminant.Real *= 10.0;
                cDeterminant.Imag *= 10.0;
                (*pExponent)--;
                Norm = NORM( cDeterminant );
            }
        }
        if (Matrix->NumberOfInterchangesIsOdd)
            CMPLX_NEGATE( cDeterminant );
        
        *pDeterminant = cDeterminant.Real;
        *piDeterminant = cDeterminant.Imag;
    }
    else
    {
	/* Real Case. */
        *pDeterminant = 1.0;

        while (++I <= Size)
        {
	    *pDeterminant /= Matrix->Diag[I]->Real;

	    /* Scale Determinant. */
            if (*pDeterminant != 0.0)
            {
		while (ABS(*pDeterminant) >= 1.0e12)
                {
		    *pDeterminant *= 1.0e-12;
                    *pExponent += 12;
                }
                while (ABS(*pDeterminant) < 1.0e-12)
                {
		    *pDeterminant *= 1.0e12;
                    *pExponent -= 12;
                }
            }
        }

	/* Scale Determinant again, this time to be between 1.0 <= x <
           10.0. */
        if (*pDeterminant != 0.0)
        {
	    while (ABS(*pDeterminant) >= 10.0)
            {
		*pDeterminant *= 0.1;
                (*pExponent)++;
            }
            while (ABS(*pDeterminant) < 1.0)
            {
		*pDeterminant *= 10.0;
                (*pExponent)--;
            }
        }
        if (Matrix->NumberOfInterchangesIsOdd)
            *pDeterminant = -*pDeterminant;
    }
}
#endif /* DETERMINANT */









#if STRIP
/*
 *  STRIP FILL-INS FROM MATRIX
 *
 *  Strips the matrix of all fill-ins.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (char *)
 *      Pointer to the matrix to be stripped.
 *
 *  >>> Local variables:
 *  pElement  (ElementPtr)
 *      Pointer that is used to step through the matrix.
 *  ppElement  (ElementPtr *)
 *      Pointer to the location of an ElementPtr.  This location will be
 *      updated if a fill-in is stripped from the matrix.
 *  pFillin  (ElementPtr)
 *      Pointer used to step through the lists of fill-ins while marking them.
 *  pLastFillin  (ElementPtr)
 *      A pointer to the last fill-in in the list.  Used to terminate a loop.
 *  pListNode  (struct  FillinListNodeStruct *)
 *      A pointer to a node in the FillinList linked-list.
 */

void
spStripFills( eMatrix )

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. */
    {
	ElementPtr  pFillin, pLastFillin;

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

        while (pListNode != NULL)
        {
	    pFillin = pListNode->pFillinList;
            pLastFillin = &(pFillin[ pListNode->NumberOfFillinsInList - 1 ]);
            while (pFillin <= pLastFillin)
                (pFillin++)->Row = 0;
            pListNode = pListNode->Next;
        }
    }

    /* Unlink fill-ins by searching for elements marked with Row = 0. */
    {
	ElementPtr pElement, *ppElement;
	int  I, Size = Matrix->Size;

	/* Unlink fill-ins in all columns. */
        for (I = 1; I <= Size; I++)
        {
	    ppElement = &(Matrix->FirstInCol[I]);
            while ((pElement = *ppElement) != NULL)
            {
		if (pElement->Row == 0)
                {
		    *ppElement = pElement->NextInCol;  /* Unlink fill-in. */
                    if (Matrix->Diag[pElement->Col] == pElement)
                        Matrix->Diag[pElement->Col] = NULL;
                }
                else
                    ppElement = &pElement->NextInCol;  /* Skip element. */
            }
        }

	/* Unlink fill-ins in all rows. */
        for (I = 1; I <= Size; I++)
        {
	    ppElement = &(Matrix->FirstInRow[I]);
            while ((pElement = *ppElement) != NULL)
            {
		if (pElement->Row == 0)
                    *ppElement = pElement->NextInRow;  /* Unlink fill-in. */
                else
                    ppElement = &pElement->NextInRow;  /* Skip element. */
            }
        }
    }
    return;
}

/* Same as above, but strips entire matrix without destroying the
 * frame.  This assumes that the matrix will be replaced with one of
 * the same size.  */
void
spStripMatrix( 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. */
    {
	ElementPtr  pElement;
	struct ElementListNodeStruct  *pListNode;

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

    /* Reset the fill-in lists. */
    {
	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 */
    {
	int  I, Size = Matrix->Size;
        for (I = 1; I <= Size; I++)
        {
	    Matrix->FirstInRow[I] = NULL;
            Matrix->FirstInCol[I] = NULL;
	    Matrix->Diag[I] = NULL;
	}
    }
}
#endif







#if TRANSLATE && DELETE
/*
 *  DELETE A ROW AND COLUMN FROM THE MATRIX
 *
 *  Deletes a row and a column from a matrix.
 *
 *  Sparse will abort if an attempt is made to delete a row or column that
 *  doesn't exist.
 *
 *  >>> Arguments:
 *  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(char *eMatrix, int Row, int Col )
{
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    ElementPtr  pElement, *ppElement, pLastElement;
    int  Size, ExtRow, ExtCol;
    ElementPtr  spcFindElementInCol();

    /* Begin `spDeleteRowAndCol'. */

    assert( IS_SPARSE(Matrix) && Row > 0 && Col > 0 );
    assert( Row <= Matrix->ExtSize && Col <= Matrix->ExtSize );

    Size = Matrix->Size;
    ExtRow = Row;
    ExtCol = Col;
    if (!Matrix->RowsLinked)
	spcLinkRows( Matrix );

    Row = Matrix->ExtToIntRowMap[Row];
    Col = Matrix->ExtToIntColMap[Col];
    assert( Row > 0 && Col > 0 );

    /* Move Row so that it is the last row in the matrix. */
    if (Row != Size) spcRowExchange( Matrix, Row, Size );

    /* Move Col so that it is the last column in the matrix. */
    if (Col != Size) spcColExchange( Matrix, Col, Size );

    /* Correct Diag pointers. */
    if (Row == Col)
        SWAP( ElementPtr, Matrix->Diag[Row], Matrix->Diag[Size] );
    else {
	Matrix->Diag[Row] = spcFindElementInCol(Matrix,
						Matrix->FirstInCol+Row,
						Row, Row, NO );
	Matrix->Diag[Col] = spcFindElementInCol(Matrix,
						Matrix->FirstInCol+Col,
						Col, Col, NO );
    }

    /* Delete last row and column of the matrix.  */
    /* Break the column links to every element in the last row. */
    pLastElement = Matrix->FirstInRow[ Size ];
    while (pLastElement != NULL) {
	ppElement = &(Matrix->FirstInCol[ pLastElement->Col ]);
        while ((pElement = *ppElement) != NULL) {
	    if (pElement == pLastElement)
                *ppElement = NULL;  /* Unlink last element in column. */
            else
                ppElement = &pElement->NextInCol;  /* Skip element. */
        }
        pLastElement = pLastElement->NextInRow;
    }

    /* Break the row links to every element in the last column. */
    pLastElement = Matrix->FirstInCol[ Size ];
    while (pLastElement != NULL) {
	ppElement = &(Matrix->FirstInRow[ pLastElement->Row ]);
        while ((pElement = *ppElement) != NULL) {
	    if (pElement == pLastElement)
                *ppElement = NULL; /* Unlink last element in row. */
            else
                ppElement = &pElement->NextInRow; /* Skip element. */
        }
        pLastElement = pLastElement->NextInCol;
    }

    /* Clean up some details. */
    Matrix->Size = Size - 1;
    Matrix->Diag[Size] = NULL;
    Matrix->FirstInRow[Size] = NULL;
    Matrix->FirstInCol[Size] = NULL;
    Matrix->CurrentSize--;
    Matrix->ExtToIntRowMap[ExtRow] = -1;
    Matrix->ExtToIntColMap[ExtCol] = -1;
    Matrix->NeedsOrdering = YES;

    return;
}
#endif