/*
* MATRIX SPALLOCATION MODULE
*
* Author: Advising professor:
* Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
* UC Berkeley
*
* This file contains the allocation and deallocation routines for the
* sparse matrix routines.
*
* >>> User accessible functions contained in this file:
* spCreate
* spDestroy
* spError
* 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,86,87,88
* by Kenneth S. Kundert and the University of California.
*
* Permission to use, copy, modify, and distribute this software and
* its documentation for any purpose and without fee is hereby granted,
* provided that the copyright notices appear in all copies and
* supporting documentation and that the authors and the University of
* California are properly credited. The authors and the University of
* California make no representations as to the suitability of this
* software for any purpose. It is provided `as is', without express
* or implied warranty.
*/
/*
* 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 "spConfig.h"
#include "spmatrix.h"
#include "spDefs.h"
#include "spmalloc.h"
#include "spAllocate.h"
#include "MALLOC.h"
static int InitializeElementBlocks( MatrixPtr Matrix, int InitialNumberOfElements, int NumberOfFillinsExpected );
static int RecordAllocation( MatrixPtr Matrix, char *AllocatedPtr );
static int AllocateBlockOfAllocationList( MatrixPtr Matrix );
/*
* MATRIX SPALLOCATION
*
* Allocates and initializes the data structures associated with a matrix.
*
* >>> Returned:
* A pointer to the matrix is returned cast into the form of a pointer to
* a character. This pointer is then passed and used by the other matrix
* routines to refer to a particular matrix. If an error occurs, the NULL
* pointer is returned.
*
* >>> Arguments:
* Size (int)
* Size of matrix or estimate of size of matrix if matrix is EXPANDABLE.
* Complex (int)
* Type of matrix. If 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 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.
* pError (int *)
* Returns error flag, needed because function spError() will not work
* correctly if spCreate() returns NULL.
*
* >>> Local variables:
* AllocatedSize (int)
* The size of the matrix being allocated.
* Matrix (MatrixPtr)
* A pointer to the matrix frame being created.
*
* >>> Possible errors:
* spNO_MEMORY
* spPANIC
* Error is cleared in this routine.
*/
char *spCreate(int Size, SPBOOLEAN Complex, int *pError )
{
register unsigned SizePlusOne;
register MatrixPtr Matrix;
register int I;
int AllocatedSize;
/* Begin `spCreate'. */
/* Clear error flag. */
*pError = spOKAY;
/* Test for valid size. */
if ((Size < 0) OR (Size == 0 AND NOT EXPANDABLE))
{
*pError = spPANIC;
return NULL;
}
/* Test for valid type. */
#if NOT spCOMPLEX
if (Complex)
{
*pError = spPANIC;
return NULL;
}
#endif
#if NOT REAL
if (NOT Complex)
{
*pError = spPANIC;
return NULL;
}
#endif
/* Create Matrix. */
AllocatedSize = Max( Size, MINIMUM_ALLOCATED_SIZE );
SizePlusOne = (unsigned)(AllocatedSize + 1);
if ((Matrix = SPALLOC(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, (char *)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. */
SPCALLOC( Matrix->Diag, ElementPtr, SizePlusOne);
if (Matrix->Diag == NULL)
{
goto MemoryError;
}
/* Allocate space in memory for FirstInCol pointer vector. */
SPCALLOC( Matrix->FirstInCol, ElementPtr, SizePlusOne);
if (Matrix->FirstInCol == NULL)
{
goto MemoryError;
}
/* Allocate space in memory for FirstInRow pointer vector. */
SPCALLOC( Matrix->FirstInRow, ElementPtr, SizePlusOne);
if (Matrix->FirstInRow == NULL)
{
goto MemoryError;
}
/* Allocate space in memory for IntToExtColMap vector. */
if (( Matrix->IntToExtColMap = SPALLOC(int, SizePlusOne)) == NULL)
{
goto MemoryError;
}
/* Allocate space in memory for IntToExtRowMap vector. */
if (( Matrix->IntToExtRowMap = SPALLOC(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 = SPALLOC(int, SizePlusOne)) == NULL)
{
goto MemoryError;
}
/* Allocate space in memory for ExtToIntRowMap vector. */
if (( Matrix->ExtToIntRowMap = SPALLOC(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 SPALLOCATION
*
* 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 (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 = SPALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, (char *)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 SPALLOCATION 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 (MatrixPtr)
* Pointer to the matrix.
* InitialNumberOfElements (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 (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 int
InitializeElementBlocks( MatrixPtr Matrix, int InitialNumberOfElements,
int NumberOfFillinsExpected )
{
ElementPtr pElement;
/* Begin `InitializeElementBlocks'. */
/* Allocate block of MatrixElements for elements. */
pElement = SPALLOC(struct MatrixElement, InitialNumberOfElements);
RecordAllocation( Matrix, (char *)pElement );
if (Matrix->Error == spNO_MEMORY)
{
return 0;
}
Matrix->ElementsRemaining = InitialNumberOfElements;
Matrix->NextAvailElement = pElement;
/* Allocate block of MatrixElements for fill-ins. */
pElement = SPALLOC(struct MatrixElement, NumberOfFillinsExpected);
RecordAllocation( Matrix, (char *)pElement );
if (Matrix->Error == spNO_MEMORY)
{
return 0;
}
Matrix->FillinsRemaining = NumberOfFillinsExpected;
Matrix->NextAvailFillin = pElement;
/* Allocate a fill-in list structure. */
Matrix->FirstFillinListNode = SPALLOC(struct FillinListNodeStruct, 1);
RecordAllocation( Matrix, (char *)Matrix->FirstFillinListNode );
if (Matrix->Error == spNO_MEMORY)
{
return 0;
}
Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
Matrix->FirstFillinListNode->pFillinList = pElement;
Matrix->FirstFillinListNode->NumberOfFillinsInList = NumberOfFillinsExpected;
Matrix->FirstFillinListNode->Next = NULL;
return 0;
}
/*
* FILL-IN SPALLOCATION
*
* 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.
*
* >>> Arguments:
* Matrix (MatrixPtr)
* Pointer to matrix.
*
* >>> Possible errors:
* spNO_MEMORY
*/
ElementPtr spcGetFillin( MatrixPtr Matrix )
{
struct FillinListNodeStruct *pListNode;
ElementPtr pFillins;
/* Begin `spcGetFillin'. */
#if NOT STRIP OR LINT
if (Matrix->FillinsRemaining == 0)
{
return spcGetElement( Matrix );
}
#endif
#if STRIP OR LINT
if (Matrix->FillinsRemaining == 0)
{
pListNode = Matrix->LastFillinListNode;
/* First see if there are any stripped fill-ins left. */
if (pListNode->Next != NULL)
{
Matrix->LastFillinListNode = pListNode = pListNode->Next;
Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList;
Matrix->NextAvailFillin = pListNode->pFillinList;
}
else
{
/* Allocate block of fill-ins. */
pFillins = SPALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, (char *)pFillins );
if (Matrix->Error == spNO_MEMORY)
{
return NULL;
}
Matrix->FillinsRemaining = ELEMENTS_PER_ALLOCATION;
Matrix->NextAvailFillin = pFillins;
/* Allocate a fill-in list structure. */
pListNode->Next = SPALLOC(struct FillinListNodeStruct, 1);
RecordAllocation( Matrix, (char *)pListNode->Next );
if (Matrix->Error == spNO_MEMORY)
{
return NULL;
}
Matrix->LastFillinListNode = pListNode = pListNode->Next;
pListNode->pFillinList = pFillins;
pListNode->NumberOfFillinsInList = ELEMENTS_PER_ALLOCATION;
pListNode->Next = NULL;
}
}
#endif
/* Update Fill-in counter and return pointer to Fill-in. */
Matrix->FillinsRemaining--;
return Matrix->NextAvailFillin++;
}
/*
* RECORD A MEMORY SPALLOCATION
*
* This routine is used to record all memory allocations so that the memory
* can be freed later.
*
* >>> Arguments:
* Matrix (MatrixPtr)
* Pointer to the matrix.
* AllocatedPtr (char *)
* The pointer returned by malloc or calloc. These pointers are saved in
* a list so that they can be easily freed.
*
* >>> Possible errors:
* spNO_MEMORY
*/
static int
RecordAllocation( MatrixPtr Matrix, char *AllocatedPtr )
{
/* Begin `RecordAllocation'. */
/*
* If Allocated pointer is NULL, assume that malloc returned a NULL pointer,
* which indicates a spNO_MEMORY error.
*/
if (AllocatedPtr == NULL)
{
Matrix->Error = spNO_MEMORY;
return 0;
}
/* Allocate block of MatrixElements if necessary. */
if (Matrix->RecordsRemaining == 0)
{
AllocateBlockOfAllocationList( Matrix );
if (Matrix->Error == spNO_MEMORY)
{
SPFREE(AllocatedPtr);
return 0;
}
}
/* Add Allocated pointer to Allocation List. */
(++Matrix->TopOfAllocationList)->AllocatedPtr = AllocatedPtr;
Matrix->RecordsRemaining--;
return 0;
}
/*
* ADD A BLOCK OF SLOTS TO SPALLOCATION LIST
*
* This routine increases the size of the allocation list.
*
* >>> Arguments:
* Matrix (MatrixPtr)
* Pointer to the matrix.
*
* >>> Local variables:
* ListPtr (AllocationListPtr)
* Pointer to the list that contains the pointers to segments of memory
* that were allocated by the operating system for the current matrix.
*
* >>> Possible errors:
* spNO_MEMORY
*/
static int
AllocateBlockOfAllocationList( MatrixPtr Matrix )
{
register int I;
register AllocationListPtr ListPtr;
/* Begin `AllocateBlockOfAllocationList'. */
/* Allocate block of records for allocation list. */
ListPtr = SPALLOC(struct AllocationRecord, (ELEMENTS_PER_ALLOCATION + 1));
if (ListPtr == NULL)
{
Matrix->Error = spNO_MEMORY;
return 0;
}
/* String entries of allocation list into singly linked list. List is linked
such that any record points to the one before it. */
ListPtr->NextRecord = Matrix->TopOfAllocationList;
Matrix->TopOfAllocationList = ListPtr;
ListPtr += ELEMENTS_PER_ALLOCATION;
for (I = ELEMENTS_PER_ALLOCATION; I > 0; I--)
{
ListPtr->NextRecord = ListPtr - 1;
ListPtr--;
}
/* Record allocation of space for allocation list on allocation list. */
Matrix->TopOfAllocationList->AllocatedPtr = (char *)ListPtr;
Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION;
return 0;
}
/*
* MATRIX DEALLOCATION
*
* Deallocates pointers and elements of Matrix.
*
* >>> Arguments:
* Matrix (char *)
* Pointer to the matrix frame which is to be removed from memory.
*
* >>> Local variables:
* ListPtr (AllocationListPtr)
* Pointer into the linked list of pointers to allocated data structures.
* Points to pointer to structure to be freed.
* NextListPtr (AllocationListPtr)
* Pointer into the linked list of pointers to allocated data structures.
* Points to the next pointer to structure to be freed. This is needed
* because the data structure to be freed could include the current node
* in the allocation list.
*/
void spDestroy(register char *eMatrix)
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
register AllocationListPtr ListPtr, NextListPtr;
/* Begin `spDestroy'. */
ASSERT( IS_SPARSE( Matrix ) );
/* Deallocate the vectors that are located in the matrix frame. */
SPFREE( Matrix->IntToExtColMap );
SPFREE( Matrix->IntToExtRowMap );
SPFREE( Matrix->ExtToIntColMap );
SPFREE( Matrix->ExtToIntRowMap );
SPFREE( Matrix->Diag );
SPFREE( Matrix->FirstInRow );
SPFREE( Matrix->FirstInCol );
SPFREE( Matrix->MarkowitzRow );
SPFREE( Matrix->MarkowitzCol );
SPFREE( Matrix->MarkowitzProd );
SPFREE( Matrix->DoCmplxDirect );
SPFREE( Matrix->DoRealDirect );
SPFREE( Matrix->Intermediate );
/* Sequentially step through the list of allocated pointers freeing pointers
* along the way. */
ListPtr = Matrix->TopOfAllocationList;
while (ListPtr != NULL )
{
char *LocPtr;
/* dans certain cas le pointeur ds la zone a desalouer
se trouve lui meme ds la dite zone en fait quand
( ListPtr == ListPtr->AllocatedPtr )
donc un free(x) suivit de x=0
fait que l'on essaye d'ecrire ds une zone que l'on vient de desalouer
ce qui plante sur linux
fprintf(stderr,"Warning bad SPFREE\n");
je regle le probleme en mettant a zero avant le free !
*/
NextListPtr = ListPtr->NextRecord;
/* BUGUED : SPFREE( ListPtr->AllocatedPtr) */
LocPtr = ListPtr->AllocatedPtr;
ListPtr->AllocatedPtr = NULL;
if ( LocPtr != NULL)
{
FREE(LocPtr);
}
ListPtr = NextListPtr;
}
return;
}
/*
* RETURN MATRIX ERROR STATUS
*
* This function is used to determine the error status of the given matrix.
*
* >>> Returned:
* The error status of the given matrix.
*
* >>> Arguments:
* eMatrix (char *)
* The matrix for which the error status is desired.
*/
int spError( char *eMatrix )
{
/* Begin `spError'. */
if (eMatrix != NULL)
{
ASSERT(((MatrixPtr)eMatrix)->ID == SPARSE_ID);
return ((MatrixPtr)eMatrix)->Error;
}
else
{
return spNO_MEMORY;
} /* This error may actually be spPANIC,
* no way to tell. */
}
/*
* WHERE IS MATRIX SINGULAR
*
* This function returns the row and column number where the matrix was
* detected as singular or where a zero was detected on the diagonal.
*
* >>> Arguments:
* eMatrix (char *)
* The matrix for which the error status is desired.
* pRow (int *)
* The row number.
* pCol (int *)
* The column number.
*/
void spWhereSingular( char *eMatrix, int *pRow, int *pCol )
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
/* Begin `spWhereSingular'. */
ASSERT( IS_SPARSE( Matrix ) );
if (Matrix->Error == spSINGULAR OR Matrix->Error == spZERO_DIAG)
{
*pRow = Matrix->SingularRow;
*pCol = Matrix->SingularCol;
}
else
{
*pRow = *pCol = 0;
}
return;
}
/*
* MATRIX SIZE
*
* Returns the size of the matrix. Either the internal or external size of
* the matrix is returned.
*
* >>> Arguments:
* eMatrix (char *)
* Pointer to matrix.
* External (SPBOOLEAN)
* If External is set true, the external size , i.e., the value of the
* largest external row or column number encountered is returned.
* Otherwise the true size of the matrix is returned. These two sizes
* may differ if the TRANSLATE option is set true.
*/
int spGetSize(char *eMatrix, SPBOOLEAN External )
{
MatrixPtr Matrix = (MatrixPtr)eMatrix;
/* Begin `spGetSize'. */
ASSERT( IS_SPARSE( Matrix ) );
#if TRANSLATE
if (External)
{
return Matrix->ExtSize;
}
else
{
return Matrix->Size;
}
#else
return Matrix->Size;
#endif
}
/*
* SET MATRIX COMPLEX OR REAL
*
* Forces matrix to be either real or complex.
*
* >>> Arguments:
* eMatrix (char *)
* Pointer to matrix.
*/
void spSetReal( char *eMatrix )
{
/* Begin `spSetReal'. */
ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) AND REAL);
((MatrixPtr)eMatrix)->Complex = NO;
return;
}
void spSetComplex( char *eMatrix )
{
/* Begin `spSetComplex'. */
ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) AND spCOMPLEX);
((MatrixPtr)eMatrix)->Complex = YES;
return;
}
/*
* ELEMENT OR FILL-IN COUNT
*
* Two functions used to return simple statistics. Either the number
* of total elements, or the number of fill-ins can be returned.
*
* >>> Arguments:
* eMatrix (char *)
* Pointer to matrix.
*/
int spFillinCount( char *eMatrix )
{
/* Begin `spFillinCount'. */
ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) );
return ((MatrixPtr)eMatrix)->Fillins;
}
int spElementCount( char *eMatrix )
{
/* Begin `spElementCount'. */
ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) );
return ((MatrixPtr)eMatrix)->Elements;
}