sputils.c

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

/* Begin `spMultTransposed'. */
    ASSERT_IS_SPARSE( Matrix );
    ASSERT_IS_NOT_FACTORED( Matrix );
    if (NOT Matrix->InternalVectorsAllocated)
        spcCreateInternalVectors( Matrix );

#if spCOMPLEX
    if (Matrix->Complex)
    {   ComplexTransposedMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS );
        return;
    }
#endif

#if REAL
#if NOT ARRAY_OFFSET
/* Correct array pointers for ARRAY_OFFSET. */
    --RHS;
    --Solution;
#endif

/* 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 /* REAL */
}
#endif /* MULTIPLICATION AND TRANSPOSE */







#if spCOMPLEX AND MULTIPLICATION AND 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.
 *      This is only the real portion of the right-hand side if the matrix
 *      is complex and spSEPARATED_COMPLEX_VECTORS is set true.
 *  Solution  <input>  (RealVector)
 *      Solution is the vector being multiplied by the matrix. This is only
 *      the real portion if the matrix is complex and
 *      spSEPARATED_COMPLEX_VECTORS is set true.
 *  iRHS  <output>  (RealVector)
 *      iRHS is the imaginary portion of the right hand side. This is
 *      what is being solved for.  This is only necessary if the matrix is
 *      complex and spSEPARATED_COMPLEX_VECTORS is true.
 *  iSolution  <input>  (RealVector)
 *      iSolution is the imaginary portion of the vector being multiplied
 *      by the matrix. This is only necessary if the matrix is
 *      complex and spSEPARATED_COMPLEX_VECTORS is true.
 *
 *  >>> Obscure Macros
 *  IMAG_VECTORS
 *      Replaces itself with `, iRHS, iSolution' if the options spCOMPLEX and
 *      spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears
 *      without a trace.
 */

static void
ComplexTransposedMatrixMultiply(
    MatrixPtr  Matrix,
    RealVector RHS,
    RealVector Solution
#if spSEPARATED_COMPLEX_VECTORS
    , RealVector iRHS
    , RealVector iSolution
#endif
)
{
register  ElementPtr  pElement;
register  ComplexVector  Vector;
ComplexNumber  Sum;
register  int  I, *pExtOrder;

/* Begin `ComplexTransposedMatrixMultiply'. */

/* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET
#if spSEPARATED_COMPLEX_VECTORS
    --RHS;              --iRHS;
    --Solution;         --iSolution;
#else
    RHS -= 2;           Solution -= 2;
#endif
#endif

/* Initialize Intermediate vector with reordered Solution vector. */
    Vector = (ComplexVector)Matrix->Intermediate;
    pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size];

#if spSEPARATED_COMPLEX_VECTORS
    for (I = Matrix->Size; I > 0; I--)
    {   Vector[I].Real = Solution[*pExtOrder];
        Vector[I].Imag = iSolution[*(pExtOrder--)];
    }
#else
    for (I = Matrix->Size; I > 0; I--)
        Vector[I] = ((ComplexVector)Solution)[*(pExtOrder--)];
#endif

    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;
        }

#if spSEPARATED_COMPLEX_VECTORS
        RHS[*pExtOrder] = Sum.Real;
        iRHS[*pExtOrder--] = Sum.Imag;
#else
        ((ComplexVector)RHS)[*pExtOrder--] = Sum;
#endif
    }
    return;
}
#endif /* spCOMPLEX AND MULTIPLICATION AND TRANSPOSE */








#if 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
 *  point number.  For this reason the determinant is scaled to a
 *  reasonable value and the logarithm of the scale factor is returned.
 *
 *  \param eMatrix
 *      A pointer to the matrix for which the determinant is desired.
 *  \param pExponent
 *      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.
 *  \param pDeterminant
 *      The real portion of the determinant.   This number is scaled to be
 *      greater than or equal to 1.0 and less than 10.0.
 *  \param piDeterminant
 *      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 register for speed.
 *  Temp  (RealNumber)
 *      Temporary storage for real portion of determinant.
 */

void
spDeterminant(
    spMatrix eMatrix,
    int *pExponent,
    spREAL *pDeterminant
#if spCOMPLEX
    , spREAL *piDeterminant
#endif
)
{
register MatrixPtr  Matrix = (MatrixPtr)eMatrix;
register int I, Size;
#if spCOMPLEX
RealNumber Norm, nr, ni;
ComplexNumber Pivot, cDeterminant;
#endif
#define  NORM(a)     (nr = ABS((a).Real), ni = ABS((a).Imag), MAX (nr,ni))

/* Begin `spDeterminant'. */
    ASSERT_IS_SPARSE( Matrix );
    ASSERT_NO_ERRORS( Matrix );
    ASSERT_IS_FACTORED( Matrix );
    *pExponent = 0;

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

    Size = Matrix->Size;
    I = 0;

#if spCOMPLEX
    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;
    }
#endif /* spCOMPLEX */
#if REAL AND spCOMPLEX
    else
#endif
#if REAL
    {   /* 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 /* REAL */
}
#endif /* DETERMINANT */








#if STRIP

/*!
 *  Strips the matrix of all fill-ins.
 *
 *  \param eMatrix
 *      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( spMatrix 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. */
    {   register  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. */
    {   register  ElementPtr pElement, *ppElement;
        register  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. */
            }