sputils.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

        }
    }
    return;
}
#endif







#if TRANSLATE AND DELETE
/*!
 *  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.
 *
 *  \param eMatrix
 *      Pointer to the matrix in which the row and column are to be deleted.
 *  \param Row
 *      Row to be deleted.
 *  \param Col
 *      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(
    spMatrix eMatrix,
    int  Row,
    int  Col
)
{
MatrixPtr  Matrix = (MatrixPtr)eMatrix;
register  ElementPtr  pElement, *ppElement, pLastElement;
int  Size, ExtRow, ExtCol;

/* Begin `spDeleteRowAndCol'. */
    ASSERT_IS_SPARSE( Matrix );
    vASSERT( (Row > 0) AND (Col > 0), "Nonpositive row or column number" );
    vASSERT( (Row <= Matrix->ExtSize) AND (Col <= Matrix->ExtSize),
             "Row or column number too large" );

    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] = spcFindDiag( Matrix, Row );
        Matrix->Diag[Col] = spcFindDiag( Matrix, Col );
    }

/*
 * 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
/*!
 *  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.
 *
 *  \return
 *  The magnitude of the ratio of the largest to smallest pivot used during
 *  previous factorization.  If the matrix was singular, zero is returned.
 *
 *  \param eMatrix
 *      Pointer to the matrix.
 */

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

    /* Begin `spPseudoCondition'. */
    ASSERT_IS_SPARSE( Matrix );
    ASSERT_NO_ERRORS( Matrix );
    ASSERT_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
/*!
 *  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.
 *
 *  \b 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.
 *
 *  \return
 *  The reciprocal of the condition number.  If the matrix was singular,
 *  zero is returned.
 *
 *  \param eMatrix
 *      Pointer to the matrix.
 *  \param NormOfMatrix
 *      The L-infinity norm of the unfactored matrix as computed by
 *      spNorm().
 *  \param pError
 *      Used to return error code.  Possible errors include \a spSINGULAR
 *      or \a spNO_MEMORY.
 */

spREAL
spCondition( 
    spMatrix eMatrix,
    spREAL 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 );
    ASSERT_NO_ERRORS( Matrix );
    ASSERT_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]) );
            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];