spallocate.c

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

/*
 *  MATRIX ALLOCATION MODULE
 *
 *  Author:                     Advising professor:
 *      Kenneth S. Kundert          Alberto Sangiovanni-Vincentelli
 *      UC Berkeley
 */
/*!\file
 *  This file contains functions for allocating and freeing matrices, configuring them, and for
 *  accessing global information about the matrix (size, error status, etc.).
 *
 *  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:
 *  spCreate
 *  spDestroy
 *  spErrorState
 *  spWhereSingular
 *  spGetSize
 *  spSetReal
 *  spSetComplex
 *  spFillinCount
 *  spElementCount
 *
 *  >>> Other functions contained in this file:
 *  spcGetElement
 *  InitializeElementBlocks
 *  spcGetFillin
 *  RecordAllocation
 *  AllocateBlockOfAllocationList
 *  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"





/*
 * Global strings
 */

char spcMatrixIsNotValid[] = "Matrix passed to Sparse is not valid";
char spcErrorsMustBeCleared[] = "Error not cleared";
char spcMatrixMustBeFactored[] = "Matrix must be factored";
char spcMatrixMustNotBeFactored[] = "Matrix must not be factored";




/*
 *  Function declarations
 */

#if 0 //not used so eliminate warning - JLM
static spError ReserveElements( MatrixPtr, int );
#endif
static void InitializeElementBlocks( MatrixPtr, int, int );
static void RecordAllocation( MatrixPtr, void* );
static void AllocateBlockOfAllocationList( MatrixPtr );



/*!
 *  Allocates and initializes the data structures associated with a matrix.
 *
 *  \return
 *  A pointer to the matrix is returned cast into \a spMatrix (typically a
 *  pointer to a void).  This pointer is then passed and used by the other
 *  matrix routines to refer to a particular matrix.  If an error occurs,
 *  the \a NULL pointer is returned.
 *
 *  \param Size
 *      Size of matrix or estimate of size of matrix if matrix is \a EXPANDABLE.
 *  \param Complex
 *      Type of matrix.  If \a Complex is 0 then the matrix is real, otherwise
 *      the matrix will be complex.  Note that if the routines are not set up
 *      to handle the type of matrix requested, then an \a spPANIC error will occur.
 *      Further note that if a matrix will be both real and complex, it must
 *      be specified here as being complex.
 *  \param pError
 *      Returns error flag, needed because function \a spErrorState() will
 *      not work correctly if \a spCreate() returns \a NULL. Possible errors
 *      include \a  spNO_MEMORY and \a spPANIC.
 */
/*  >>> Local variables:
 *  AllocatedSize  (int)
 *      The size of the matrix being allocated.
 *  Matrix  (MatrixPtr)
 *      A pointer to the matrix frame being created.
 */

spMatrix
spCreate(
    int  Size,
    int  Complex,
    spError *pError
)
{
register  unsigned  SizePlusOne;
register  MatrixPtr  Matrix;
register  int  I;
int  AllocatedSize;

/* Begin `spCreate'. */
/* Clear error flag. */
    *pError = spOKAY;

/* Test for valid size. */
    vASSERT( (Size >= 0) AND (Size != 0 OR EXPANDABLE), "Invalid size" );

/* Test for valid type. */
#if NOT spCOMPLEX
    ASSERT( NOT Complex );
#endif
#if NOT REAL
    ASSERT( Complex );
#endif

/* Create Matrix. */
    AllocatedSize = MAX( Size, MINIMUM_ALLOCATED_SIZE );
    SizePlusOne = (unsigned)(AllocatedSize + 1);

    if ((Matrix = ALLOC(struct MatrixFrame, 1)) == NULL)
    {   *pError = spNO_MEMORY;
        return NULL;
    }

/* Initialize matrix */
    Matrix->ID = SPARSE_ID;
    Matrix->Complex = Complex;
    Matrix->PreviousMatrixWasComplex = Complex;
    Matrix->Factored = NO;
    Matrix->Elements = 0;
    Matrix->Error = *pError;
    Matrix->Fillins = 0;
    Matrix->Reordered = NO;
    Matrix->NeedsOrdering = YES;
    Matrix->NumberOfInterchangesIsOdd = NO;
    Matrix->Partitioned = NO;
    Matrix->RowsLinked = NO;
    Matrix->InternalVectorsAllocated = NO;
    Matrix->SingularCol = 0;
    Matrix->SingularRow = 0;
    Matrix->Size = Size;
    Matrix->AllocatedSize = AllocatedSize;
    Matrix->ExtSize = Size;
    Matrix->AllocatedExtSize = AllocatedSize;
    Matrix->CurrentSize = 0;
    Matrix->ExtToIntColMap = NULL;
    Matrix->ExtToIntRowMap = NULL;
    Matrix->IntToExtColMap = NULL;
    Matrix->IntToExtRowMap = NULL;
    Matrix->MarkowitzRow = NULL;
    Matrix->MarkowitzCol = NULL;
    Matrix->MarkowitzProd = NULL;
    Matrix->DoCmplxDirect = NULL;
    Matrix->DoRealDirect = NULL;
    Matrix->Intermediate = NULL;
    Matrix->RelThreshold = DEFAULT_THRESHOLD;
    Matrix->AbsThreshold = 0.0;

    Matrix->TopOfAllocationList = NULL;
    Matrix->RecordsRemaining = 0;
    Matrix->ElementsRemaining = 0;
    Matrix->FillinsRemaining = 0;

    RecordAllocation( Matrix, (void *)Matrix );
    if (Matrix->Error == spNO_MEMORY) goto MemoryError;

/* Take out the trash. */
    Matrix->TrashCan.Real = 0.0;
#if spCOMPLEX
    Matrix->TrashCan.Imag = 0.0;
#endif
    Matrix->TrashCan.Row = 0;
    Matrix->TrashCan.Col = 0;
    Matrix->TrashCan.NextInRow = NULL;
    Matrix->TrashCan.NextInCol = NULL;
#if INITIALIZE
    Matrix->TrashCan.pInitInfo = NULL;
#endif

/* Allocate space in memory for Diag pointer vector. */
    CALLOC( Matrix->Diag, ElementPtr, SizePlusOne);
    if (Matrix->Diag == NULL)
        goto MemoryError;

/* Allocate space in memory for FirstInCol pointer vector. */
    CALLOC( Matrix->FirstInCol, ElementPtr, SizePlusOne);
    if (Matrix->FirstInCol == NULL)
        goto MemoryError;

/* Allocate space in memory for FirstInRow pointer vector. */
    CALLOC( Matrix->FirstInRow, ElementPtr, SizePlusOne);
    if (Matrix->FirstInRow == NULL)
        goto MemoryError;

/* Allocate space in memory for IntToExtColMap vector. */
    if (( Matrix->IntToExtColMap = ALLOC(int, SizePlusOne)) == NULL)
        goto MemoryError;

/* Allocate space in memory for IntToExtRowMap vector. */
    if (( Matrix->IntToExtRowMap = ALLOC(int, SizePlusOne)) == NULL)
        goto MemoryError;

/* Initialize MapIntToExt vectors. */
    for (I = 1; I <= AllocatedSize; I++)
    {   Matrix->IntToExtRowMap[I] = I;
        Matrix->IntToExtColMap[I] = I;
    }

#if TRANSLATE
/* Allocate space in memory for ExtToIntColMap vector. */
    if (( Matrix->ExtToIntColMap = ALLOC(int, SizePlusOne)) == NULL)
        goto MemoryError;

/* Allocate space in memory for ExtToIntRowMap vector. */
    if (( Matrix->ExtToIntRowMap = ALLOC(int, SizePlusOne)) == NULL)
        goto MemoryError;

/* Initialize MapExtToInt vectors. */
    for (I = 1; I <= AllocatedSize; I++)
    {   Matrix->ExtToIntColMap[I] = -1;
        Matrix->ExtToIntRowMap[I] = -1;
    }
    Matrix->ExtToIntColMap[0] = 0;
    Matrix->ExtToIntRowMap[0] = 0;
#endif

/* Allocate space for fill-ins and initial set of elements. */
    InitializeElementBlocks( Matrix, SPACE_FOR_ELEMENTS*AllocatedSize,
                                     SPACE_FOR_FILL_INS*AllocatedSize );
    if (Matrix->Error == spNO_MEMORY)
        goto MemoryError;

    return (char *)Matrix;

MemoryError:

/* Deallocate matrix and return no pointer to matrix if there is not enough
   memory. */
    *pError = spNO_MEMORY;
    spDestroy( (char *)Matrix);
    return NULL;
}









/*
 *  ELEMENT ALLOCATION
 *
 *  This routine allocates space for matrix elements. It requests large blocks
 *  of storage from the system and doles out individual elements as required.
 *  This technique, as opposed to allocating elements individually, tends to
 *  speed the allocation process.
 *
 *  >>> Returned:
 *  A pointer to an element.
 *
 *  >>> Arguments:
 *  Matrix  <input>  (MatrixPtr)
 *      Pointer to matrix.
 *
 *  >>> Local variables:
 *  pElement  (ElementPtr)
 *      A pointer to the first element in the group of elements being allocated.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 */

ElementPtr
spcGetElement( MatrixPtr Matrix )
{
ElementPtr  pElement;

/* Begin `spcGetElement'. */

/* Allocate block of MatrixElements if necessary. */
    if (Matrix->ElementsRemaining == 0)
    {   pElement = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
        RecordAllocation( Matrix, (void *)pElement );
        if (Matrix->Error == spNO_MEMORY) return NULL;
        Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION;
        Matrix->NextAvailElement = pElement;
    }

/* Update Element counter and return pointer to Element. */
    Matrix->ElementsRemaining--;
    return Matrix->NextAvailElement++;
}








/*
 *  ELEMENT ALLOCATION INITIALIZATION
 *
 *  This routine allocates space for matrix fill-ins and an initial set of
 *  elements.  Besides being faster than allocating space for elements one
 *  at a time, it tends to keep the fill-ins physically close to the other
 *  matrix elements in the computer memory.  This keeps virtual memory paging
 *  to a minimum.
 *
 *  >>> Arguments:
 *  Matrix  <input>    (MatrixPtr)
 *      Pointer to the matrix.
 *  InitialNumberOfElements  <input> (int)
 *      This number is used as the size of the block of memory, in
 *      MatrixElements, reserved for elements. If more than this number of
 *      elements are generated, then more space is allocated later.
 *  NumberOfFillinsExpected  <input> (int)
 *      This number is used as the size of the block of memory, in
 *      MatrixElements, reserved for fill-ins. If more than this number of
 *      fill-ins are generated, then more space is allocated, but they may
 *      not be physically close in computer's memory.
 *
 *  >>> Local variables:
 *  pElement  (ElementPtr)
 *      A pointer to the first element in the group of elements being allocated.
 *
 *  >>> Possible errors:
 *  spNO_MEMORY
 */

static void
InitializeElementBlocks(
    MatrixPtr Matrix,
    int InitialNumberOfElements,
    int NumberOfFillinsExpected
)
{
ElementPtr  pElement;

/* Begin `InitializeElementBlocks'. */

/* Allocate block of MatrixElements for elements. */
    pElement = ALLOC(struct MatrixElement, InitialNumberOfElements);
    RecordAllocation( Matrix, (void *)pElement );
    if (Matrix->Error == spNO_MEMORY) return;
    Matrix->ElementsRemaining = InitialNumberOfElements;
    Matrix->NextAvailElement = pElement;

/* Allocate block of MatrixElements for fill-ins. */
    pElement = ALLOC(struct MatrixElement, NumberOfFillinsExpected);
    RecordAllocation( Matrix, (void *)pElement );
    if (Matrix->Error == spNO_MEMORY) return;
    Matrix->FillinsRemaining = NumberOfFillinsExpected;
    Matrix->NextAvailFillin = pElement;

/* Allocate a fill-in list structure. */
    Matrix->FirstFillinListNode = ALLOC(struct FillinListNodeStruct,1);
    RecordAllocation( Matrix, (void *)Matrix->FirstFillinListNode );
    if (Matrix->Error == spNO_MEMORY) return;
    Matrix->LastFillinListNode = Matrix->FirstFillinListNode;

    Matrix->FirstFillinListNode->pFillinList = pElement;
    Matrix->FirstFillinListNode->NumberOfFillinsInList =NumberOfFillinsExpected;
    Matrix->FirstFillinListNode->Next = NULL;

    return;
}










/*
 *  FILL-IN ALLOCATION
 *
 *  This routine allocates space for matrix fill-ins. It requests large blocks
 *  of storage from the system and doles out individual elements as required.
 *  This technique, as opposed to allocating elements individually, tends to
 *  speed the allocation process.
 *
 *  >>> Returned:
 *  A pointer to the fill-in.