sputils.c

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

    /* 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]);

    FREE( Tm );

    Linpack = ASy / ASz;
    OLeary = E / MaxY;
    InvNormOfInverse = MIN( Linpack, OLeary );
    return InvNormOfInverse / NormOfMatrix;
}





/*
 *  L-INFINITY MATRIX NORM 
 *
 *  Computes the L-infinity norm of an unfactored matrix.  It is a fatal
 *  error to pass this routine a factored matrix.
 *
 *  One difficulty is that the rows may not be linked.
 *
 *  >>> Returns:
 *  The largest absolute row sum of matrix.
 *
 *  >>> Arguments:
 *  eMatrix  <input>  (char *)
 *      Pointer to the matrix.
 */

RealNumber
spNorm( eMatrix )

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

    /* Begin `spNorm'. */
    assert( IS_SPARSE(Matrix) && !IS_FACTORED(Matrix) );
    if (!Matrix->RowsLinked) spcLinkRows( Matrix );

    /* Compute row sums. */
    if (!Matrix->Complex)
    {
	for (I = Matrix->Size; I > 0; I--)
        {
	    pElement = Matrix->FirstInRow[I];
            AbsRowSum = 0.0;
            while (pElement != NULL)
            {
		AbsRowSum += ABS( pElement->Real );
                pElement = pElement->NextInRow;
            }
            if (Max < AbsRowSum) Max = AbsRowSum;
        }
    }
    else
    {
	for (I = Matrix->Size; I > 0; I--)
        {
	    pElement = Matrix->FirstInRow[I];
            AbsRowSum = 0.0;
            while (pElement != NULL)
            {
		AbsRowSum += CMPLX_1_NORM( *pElement );
                pElement = pElement->NextInRow;
            }
            if (Max < AbsRowSum) Max = AbsRowSum;
        }
    }
    return Max;
}
#endif /* CONDITION */






#if STABILITY
/*
 *  STABILITY OF FACTORIZATION
 *
 *  The following routines are used to gauge the stability of a
 *  factorization.  If the factorization is determined to be too
 *  unstable, then the matrix should be reordered.  The routines
 *  compute quantities that are needed in the computation of a bound
 *  on the error attributed to any one element in the matrix during
 *  the factorization.  In other words, there is a matrix E = [e_ij]
 *  of error terms such that A+E = LU.  This routine finds a bound on
 *  |e_ij|.  Erisman & Reid [1] showed that |e_ij| < 3.01 u rho m_ij,
 *  where u is the machine rounding unit, rho = max a_ij where the max
 *  is taken over every row i, column j, and step k, and m_ij is the
 *  number of multiplications required in the computation of l_ij if i
 *  > j or u_ij otherwise.  Barlow [2] showed that rho < max_i || l_i
 *  ||_p max_j || u_j ||_q where 1/p + 1/q = 1.
 *
 *  The first routine finds the magnitude on the largest element in
 *  the matrix.  If the matrix has not yet been factored, the largest
 *  element is found by direct search.  If the matrix is factored, a
 *  bound on the largest element in any of the reduced submatrices is
 *  computed using Barlow with p = oo and q = 1.  The ratio of these
 *  two numbers is the growth, which can be used to determine if the
 *  pivoting order is adequate.  A large growth implies that
 *  considerable error has been made in the factorization and that it
 *  is probably a good idea to reorder the matrix.  If a large growth
 *  in encountered after using spFactor(), reconstruct the matrix and
 *  refactor using spOrderAndFactor().  If a large growth is
 *  encountered after using spOrderAndFactor(), refactor using
 *  spOrderAndFactor() with the pivot threshold increased, say to 0.1.
 *
 *  Using only the size of the matrix as an upper bound on m_ij and
 *  Barlow's bound, the user can estimate the size of the matrix error
 *  terms e_ij using the bound of Erisman and Reid.  The second
 *  routine computes a tighter bound (with more work) based on work by
 *  Gear [3], |e_ij| < 1.01 u rho (t c^3 + (1 + t)c^2) where t is the
 *  threshold and c is the maximum number of off-diagonal elements in
 *  any row of L.  The expensive part of computing this bound is
 *  determining the maximum number of off-diagonals in L, which
 *  changes only when the order of the matrix changes.  This number is
 *  computed and saved, and only recomputed if the matrix is
 *  reordered.
 *
 *  [1] A. M. Erisman, J. K. Reid.  Monitoring the stability of the
 *      triangular factorization of a sparse matrix.  Numerische
 *      Mathematik.  Vol. 22, No. 3, 1974, pp 183-186.
 *
 *  [2] J. L. Barlow.  A note on monitoring the stability of triangular
 *      decomposition of sparse matrices.  "SIAM Journal of Scientific
 *      and Statistical Computing."  Vol. 7, No. 1, January 1986, pp 166-168.
 *
 *  [3] I. S. Duff, A. M. Erisman, J. K. Reid.  "Direct Methods for Sparse
 *      Matrices."  Oxford 1986. pp 99.  */

/*
 *  LARGEST ELEMENT IN MATRIX
 *
 *  >>> Returns:
 *  If matrix is not factored, returns the magnitude of the largest element in
 *  the matrix.  If the matrix is factored, a bound on the magnitude of the
 *  largest element in any of the reduced submatrices is returned.
 *
 *  >>> Arguments:
 *  eMatrix  <input>  (char *)
 *      Pointer to the matrix.
 */

RealNumber
spLargestElement( eMatrix )

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

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

    if (Matrix->Factored && !Matrix->Complex)
    {
	if (Matrix->Error == spSINGULAR) return 0.0;

	/* Find the bound on the size of the largest element over all
           factorization. */
        for (I = 1; I <= Matrix->Size; I++)
        {
	    pDiag = Matrix->Diag[I];

	    /* Lower triangular matrix. */
            Pivot = 1.0 / pDiag->Real;
            Mag = ABS( Pivot );
            if (Mag > MaxRow) MaxRow = Mag;
            pElement = Matrix->FirstInRow[I];
            while (pElement != pDiag)
            {
		Mag = ABS( pElement->Real );
                if (Mag > MaxRow) MaxRow = Mag;
                pElement = pElement->NextInRow;
            }

	    /* Upper triangular matrix. */
            pElement = Matrix->FirstInCol[I];
            AbsColSum = 1.0;  /* Diagonal of U is unity. */
            while (pElement != pDiag)
            {
		AbsColSum += ABS( pElement->Real );
                pElement = pElement->NextInCol;
            }
            if (AbsColSum > MaxCol) MaxCol = AbsColSum;
        }
    }
    else if (!Matrix->Complex)
    {
	for (I = 1; I <= Matrix->Size; I++)
        {
	    pElement = Matrix->FirstInCol[I];
            while (pElement != NULL)
            {
		Mag = ABS( pElement->Real );
                if (Mag > Max) Max = Mag;
                pElement = pElement->NextInCol;
            }
        }
        return Max;
    }
    if (Matrix->Factored && Matrix->Complex)
    {
	if (Matrix->Error == spSINGULAR) return 0.0;

	/* Find the bound on the size of the largest element over all
           factorization. */
        for (I = 1; I <= Matrix->Size; I++)
        {
	    pDiag = Matrix->Diag[I];

	    /* Lower triangular matrix. */
            CMPLX_RECIPROCAL( cPivot, *pDiag );
            Mag = CMPLX_1_NORM( cPivot );
            if (Mag > MaxRow) MaxRow = Mag;
            pElement = Matrix->FirstInRow[I];
            while (pElement != pDiag)
            {
		Mag = CMPLX_1_NORM( *pElement );
                if (Mag > MaxRow) MaxRow = Mag;
                pElement = pElement->NextInRow;
            }

	    /* Upper triangular matrix. */
            pElement = Matrix->FirstInCol[I];
            AbsColSum = 1.0;  /* Diagonal of U is unity. */
            while (pElement != pDiag)
            {
		AbsColSum += CMPLX_1_NORM( *pElement );
                pElement = pElement->NextInCol;
            }
            if (AbsColSum > MaxCol) MaxCol = AbsColSum;
        }
    }
    else if (Matrix->Complex)
    {
	for (I = 1; I <= Matrix->Size; I++)
        {
	    pElement = Matrix->FirstInCol[I];
            while (pElement != NULL)
            {
		Mag = CMPLX_1_NORM( *pElement );
                if (Mag > Max) Max = Mag;
                pElement = pElement->NextInCol;
            }
        }
        return Max;
    }
    return MaxRow * MaxCol;
}




/*
 *  MATRIX ROUNDOFF ERROR
 *
 *  >>> Returns:
 *  Returns a bound on the magnitude of the largest element in E = A - LU.
 *
 *  >>> Arguments:
 *  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( eMatrix, Rho )

char *eMatrix;
RealNumber Rho;
{
    MatrixPtr  Matrix = (MatrixPtr)eMatrix;
    ElementPtr pElement;
    int Count, I, MaxCount = 0;
    RealNumber Reid, Gear;

    /* Begin `spRoundoff'. */
    assert( IS_SPARSE(Matrix) && IS_FACTORED(Matrix) );

    /* Compute Barlow's bound if it is not given. */
    if (Rho < 0.0) Rho = spLargestElement( 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








#if DOCUMENTATION
/*
 *  SPARSE ERROR MESSAGE
 *
 *  This routine prints a short message to a stream describing the error
 *  error state of sparse.  No message is produced if there is no error.
 *
 *  >>> Arguments:
 *  eMatrix  <input>  (char *)
 *	Matrix for which the error message is to be printed.
 *  Stream  <input>  (FILE *)
 *	Stream to which