spbuild.c

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

/*
 *  MATRIX BUILD MODULE
 *
 *  Author:                     Advising professor:
 *     Kenneth S. Kundert           Alberto Sangiovanni-Vincentelli
 *     UC Berkeley
 */
/*!\file
 *  This file contains the routines associated with clearing, loading and
 *  preprocessing the matrix.
 *
 *  Objects that begin with the \a spc prefix are considered private
 *  and should not be used.
 *
 *  \author
 *  Kenneth S. Kundert <[e-mail]>
 */
/*  >>> User accessible functions contained in this file:
 *  spClear
 *  spFindElement
 *  spGetElement
 *  spGetAdmittance
 *  spGetQuad
 *  spGetOnes
 *  spInstallInitInfo
 *  spGetInitInfo
 *  spInitialize
 *
 *  >>> Other functions contained in this file:
 *  Translate
 *  spcFindDiag
 *  spcCreateElement
 *  spcLinkRows
 *  EnlargeMatrix
 *  ExpandTranslationArrays
 */


/*
 *  Revision and copyright information.
 *
 *  Copyright (c) 1985-2003 by Kenneth S. Kundert
 */

#if 0
static char copyright[] =
    "Sparse1.4: Copyright (c) 1985-2003 by Kenneth S. Kundert";
#endif




/*
 *  IMPORTS
 *
 *  >>> Import descriptions:
 *  spConfig.h
 *     Macros that customize the sparse matrix routines.
 *  spMatrix.h
 *     Macros and declarations to be imported by the user.
 *  spDefs.h
 *     Matrix type and macro definitions for the sparse matrix routines.
 */

#define spINSIDE_SPARSE
#include <stdio.h>
#include "spConfig.h"
#include "spMatrix.h"
#include "spDefs.h"





/*
 *  Function declarations
 */

static void Translate( MatrixPtr, int*, int* );
static void EnlargeMatrix( MatrixPtr, int );
static void ExpandTranslationArrays( MatrixPtr, int );






/*!
 *  Sets every element of the matrix to zero and clears the error flag.
 *
 *  \param eMatrix
 *     Pointer to matrix that is to be cleared.
 */
/*  >>> Local variables:
 *  pElement  (ElementPtr)
 *     A pointer to the element being cleared.
 */

void
spClear( spMatrix eMatrix )
{
MatrixPtr  Matrix = (MatrixPtr)eMatrix;
register  ElementPtr  pElement;
register  int  I;

/* Begin `spClear'. */
    ASSERT_IS_SPARSE( Matrix );

/* Clear matrix. */
#if spCOMPLEX
    if (Matrix->PreviousMatrixWasComplex OR Matrix->Complex)
    {   for (I = Matrix->Size; I > 0; I--)
        {   pElement = Matrix->FirstInCol[I];
            while (pElement != NULL)
            {   pElement->Real = 0.0;
                pElement->Imag = 0.0;
                pElement = pElement->NextInCol;
            }
        }
    }
    else
#endif
    {   for (I = Matrix->Size; I > 0; I--)
        {   pElement = Matrix->FirstInCol[I];
            while (pElement != NULL)
            {   pElement->Real = 0.0;
                pElement = pElement->NextInCol;
            }
        }
    }

/* Empty the trash. */
    Matrix->TrashCan.Real = 0.0;
#if spCOMPLEX
    Matrix->TrashCan.Imag = 0.0;
#endif

    Matrix->Error = spOKAY;
    Matrix->Factored = NO;
    Matrix->SingularCol = 0;
    Matrix->SingularRow = 0;
    Matrix->PreviousMatrixWasComplex = Matrix->Complex;
    return;
}










/*!
 *  This routine is used to find an element given its indices.  It will not
 *  create it if it does not exist.
 *
 *  \return
 *  A pointer to the desired element, or \a NULL if it does not exist.
 *
 *  \param eMatrix
 *      Pointer to matrix.
 *  \param Row
 *      Row index for element.
 *  \param Col
 *      Column index for element.
 *
 *  \see spGetElement()
 */
/*  >>> Local variables:
 *  pElement  (ElementPtr)
 *      Pointer to an element in the matrix.
 */

spElement *
spFindElement(
    spMatrix eMatrix,
    int Row,
    int Col
)
{
MatrixPtr  Matrix = (MatrixPtr)eMatrix;
register ElementPtr  pElement;
int StartAt=0, Min = LARGEST_INTEGER;
#define BorderRight 0   /* Start at left border, move right. */
#define BorderDown  1   /* Start at top border, move down. */
#define DiagRight   2   /* Start at diagonal, move right. */
#define DiagDown    3   /* Start at diagonal, move down. */

/* Begin `spFindElement'. */
    if (Row == Col) return &Matrix->Diag[Row]->Real;

/* Determine where to start the search. */
    if (Matrix->RowsLinked)
    {   if ((Col >= Row) AND Matrix->Diag[Row])
        {   Min = Col - Row;
            StartAt = DiagRight;
        }
        else
        {   Min = Col;
            StartAt = BorderRight;
        }
    }
    if ((Row >= Col) AND Matrix->Diag[Col])
    {   if (Row - Col < Min)
            StartAt = DiagDown;
    }
    else if (Row < Min)
        StartAt = BorderDown;
    
/* Search column for element. */
    if ((StartAt == BorderDown) OR (StartAt == DiagDown))
    {   if (StartAt == BorderDown)
            pElement = Matrix->FirstInCol[Col];
        else
            pElement = Matrix->Diag[Col];

        while ((pElement != NULL) AND (pElement->Row < Row))
            pElement = pElement->NextInCol;
        if (pElement AND (pElement->Row == Row))
            return &pElement->Real;
        else
            return NULL;
    }

/* Search row for element. */
    if (StartAt == BorderRight)
        pElement = Matrix->FirstInRow[Row];
    else
        pElement = Matrix->Diag[Row];

    while ((pElement != NULL) AND (pElement->Col < Col))
        pElement = pElement->NextInRow;
    if (pElement AND (pElement->Col == Col))
        return &pElement->Real;
    else
        return NULL;
}









/*!
 *  Finds element [Row,Col] and returns a pointer to it.  If element is
 *  not found then it is created and spliced into matrix.  This routine
 *  is only to be used after spCreate() and before spMNA_Preorder(),
 *  spFactor() or spOrderAndFactor().  Returns a pointer to the
 *  real portion of an \a spElement.  This pointer is later used by
 *  \a spADD_xxx_ELEMENT to directly access element.
 *
 *  \return
 *  Returns a pointer to the element.  This pointer is then used to directly
 *  access the element during successive builds.
 *
 *  \param eMatrix
 *     Pointer to the matrix that the element is to be added to.
 *  \param Row
 *     Row index for element.  Must be in the range of [0..Size] unless
 *     the options \a EXPANDABLE or \a TRANSLATE are used. Elements placed in
 *     row zero are discarded.  In no case may \a Row be less than zero.
 *  \param Col
 *     Column index for element.  Must be in the range of [0..Size] unless
 *     the options \a EXPANDABLE or \a TRANSLATE are used. Elements placed in
 *     column zero are discarded.  In no case may \a Col be less than zero.

 *  \see spFindElement()
 */
/*  >>> Local variables:
 *  pElement  (RealNumber *)
 *     Pointer to the element.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 *  Error is not cleared in this routine.
 */

spElement *
spGetElement(
    spMatrix eMatrix,
    int Row,
    int Col
)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
ElementPtr pElement;

/* Begin `spGetElement'. */
    ASSERT_IS_SPARSE( Matrix );
    vASSERT( Row >= 0 AND Col >= 0, "Negative row or column number" );

    if ((Row == 0) OR (Col == 0))
        return &Matrix->TrashCan.Real;

#if NOT TRANSLATE
    vASSERT( NOT Matrix->Reordered,
             "Set TRANSLATE to add elements to a reordered matrix" );
#endif

#if TRANSLATE
    Translate( Matrix, &Row, &Col );
    if (Matrix->Error == spNO_MEMORY) return NULL;
#endif

#if NOT TRANSLATE
#if NOT EXPANDABLE
    vASSERT( (Row <= Matrix->Size) AND (Col <= Matrix->Size),
             "Row or column number too large" );
#endif

#if EXPANDABLE
/* Re-size Matrix if necessary. */
    if ((Row > Matrix->Size) OR (Col > Matrix->Size))
        EnlargeMatrix( Matrix, MAX(Row, Col) );
    if (Matrix->Error == spNO_MEMORY) return NULL;
#endif
#endif

    if ((Row != Col) OR ((pElement = Matrix->Diag[Row]) == NULL))
    {   /*
         * Element does not exist or does not reside along diagonal.  Search
         * for element and if it does not exist, create it.
         */
        pElement = spcCreateElement( Matrix, Row, Col,
                                     &(Matrix->FirstInRow[Row]),
                                     &(Matrix->FirstInCol[Col]), NO );
    }
/*
 * Cast pointer into a pointer to a RealNumber.  This requires that Real
 * be the first record in the MatrixElement structure.
 */
    return &pElement->Real;
}







#if TRANSLATE

/*
 *  TRANSLATE EXTERNAL INDICES TO INTERNAL
 *
 *  Convert internal row and column numbers to internal row and column numbers.
 *  Also updates Ext/Int maps.
 *
 *
 *  >>> Arguments:
 *  Matrix  <input>    (MatrixPtr)
 *      Pointer to the matrix.
 *  Row  <input/output>  (int *)
 *     Upon entry Row is either a external row number of an external node
 *     number.  Upon entry, the internal equivalent is supplied.
 *  Col  <input/output>  (int *)
 *     Upon entry Column is either a external column number of an external node
 *     number.  Upon entry, the internal equivalent is supplied.
 *
 *  >>> Local variables:
 *  ExtCol  (int)
 *     Temporary variable used to hold the external column or node number
 *     during the external to internal column number translation.
 *  ExtRow  (int)
 *     Temporary variable used to hold the external row or node number during
 *     the external to internal row number translation.
 *  IntCol  (int)
 *     Temporary variable used to hold the internal column or node number
 *     during the external to internal column number translation.
 *  IntRow  (int)
 *     Temporary variable used to hold the internal row or node number during
 *     the external to internal row number translation.
 */

static void
Translate( 
    MatrixPtr Matrix,
    int *Row,
    int *Col
)
{
register int IntRow, IntCol, ExtRow, ExtCol;

/* Begin `Translate'. */
    ExtRow = *Row;
    ExtCol = *Col;

/* Expand translation arrays if necessary. */
    if ((ExtRow > Matrix->AllocatedExtSize) OR
        (ExtCol > Matrix->AllocatedExtSize))
    {
        ExpandTranslationArrays( Matrix, MAX(ExtRow, ExtCol) );
        if (Matrix->Error == spNO_MEMORY) return;
    }

/* Set ExtSize if necessary. */
    if ((ExtRow > Matrix->ExtSize) OR (ExtCol > Matrix->ExtSize))
        Matrix->ExtSize = MAX(ExtRow, ExtCol);

/* Translate external row or node number to internal row or node number. */
    if ((IntRow = Matrix->ExtToIntRowMap[ExtRow]) == -1)
    {   Matrix->ExtToIntRowMap[ExtRow] = ++Matrix->CurrentSize;
        Matrix->ExtToIntColMap[ExtRow] = Matrix->CurrentSize;
        IntRow = Matrix->CurrentSize;