From 7ac9278389afcac531554de6eecfed14d89bfaa2 Mon Sep 17 00:00:00 2001 From: arno Date: Mon, 3 Jul 2000 15:28:50 +0000 Subject: [PATCH] * src/include/spconfig.h: Removed spCOMPLEX, spSEPARATED_COMPLEX_VECTORS and spCOMPATIBILITY defines. This made including this file from src/include/spmatrix.h unnecessary. Moved this file to src/maths/sparse/spconfig.h. * src/include/spmatrix.h: Removed include of src/include/spconfig.h. * src/maths/sparse/spalloc.c, src/maths/sparse/spbuild.c, src/maths/sparse/spcombin.c, src/maths/sparse/spdefs.h, src/maths/sparse/spfactor.c, src/maths/sparse/spoutput.c, src/maths/sparse/spsmp.c, src/maths/sparse/spsolve.c, src/maths/sparse/sputils.c: The other files affected by the removal of spCOMPLEX, spSEPARATED_COMPLEX_VECTORS and spCOMPATIBILITY defines. Also: assertions are enabled by default. * src/include/smpdefs.h, src/maths/sparse/spsmp.c: SMPmatrix is now a typedef for void, instead of char. Updated all function declarations to match this. Also added function prototypes not previously mentioned in src/include/smpdefs.h. * src/include/complex.h: Updates of cast from char * to void * --- src/include/Makefile.am | 1 - src/include/complex.h | 11 +- src/include/smpdefs.h | 66 +- src/include/spmatrix.h | 172 +-- src/maths/sparse/Makefile.am | 21 +- src/maths/sparse/spalloc.c | 452 ++++---- src/maths/sparse/spbuild.c | 496 ++++----- src/maths/sparse/spcombin.c | 44 +- src/{include => maths/sparse}/spconfig.h | 54 - src/maths/sparse/spdefs.h | 264 +---- src/maths/sparse/spfactor.c | 1227 ++++++++++----------- src/maths/sparse/spoutput.c | 414 ++++---- src/maths/sparse/spsmp.c | 237 ++--- src/maths/sparse/spsolve.c | 325 +++--- src/maths/sparse/sputils.c | 1229 +++++++++++----------- 15 files changed, 2199 insertions(+), 2814 deletions(-) rename src/{include => maths/sparse}/spconfig.h (86%) diff --git a/src/include/Makefile.am b/src/include/Makefile.am index 4b3f6460e..846fddf20 100644 --- a/src/include/Makefile.am +++ b/src/include/Makefile.am @@ -47,7 +47,6 @@ noinst_HEADERS = \ sensgen.h \ sim.h \ smpdefs.h \ - spconfig.h \ sperror.h \ spmatrix.h \ suffix.h \ diff --git a/src/include/complex.h b/src/include/complex.h index c9c76b482..6c6d8b014 100644 --- a/src/include/complex.h +++ b/src/include/complex.h @@ -269,12 +269,9 @@ typedef struct (A).imag = (B.imag) - (C.imag); \ } -/* Macro function that returns the approx absolute value of a complex number. */ -#if spCOMPLEX +/* Macro function that returns the approx absolute value of a complex + number. */ #define ELEMENT_MAG(ptr) (ABS((ptr)->Real) + ABS((ptr)->Imag)) -#else -#define ELEMENT_MAG(ptr) ((ptr)->Real < 0.0 ? -(ptr)->Real : (ptr)->Real) -#endif /* Complex assignment statements. */ #define CMPLX_ASSIGN(to,from) \ @@ -481,8 +478,8 @@ typedef struct /* Complex reciprocation: to = 1.0 / den */ #define CMPLX_RECIPROCAL(to,den) \ { RealNumber r_; \ - if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \ - ((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \ + if (((den).Real >= (den).Imag && (den).Real > -(den).Imag) || \ + ((den).Real < (den).Imag && (den).Real <= -(den).Imag)) \ { r_ = (den).Imag / (den).Real; \ (to).Imag = -r_*((to).Real = 1.0/((den).Real + r_*(den).Imag)); \ } \ diff --git a/src/include/smpdefs.h b/src/include/smpdefs.h index bb3b6e8f5..38694cf71 100644 --- a/src/include/smpdefs.h +++ b/src/include/smpdefs.h @@ -1,7 +1,7 @@ #ifndef SMP #define SMP -typedef char SMPmatrix; +typedef void SMPmatrix; typedef struct MatrixElement *SMPelement; /********** @@ -12,52 +12,34 @@ Author: 1985 Thomas L. Quarles #include "complex.h" #include -#ifdef __STDC__ int SMPaddElt( SMPmatrix *, int , int , double ); -void SMPcClear( SMPmatrix *); -int SMPcLUfac( SMPmatrix *, double ); -int SMPcProdDiag( SMPmatrix *, SPcomplex *, int *); -int SMPcReorder( SMPmatrix * , double , double , int *); -void SMPcSolve( SMPmatrix *, double [], double [], double [], double []); -void SMPclear( SMPmatrix *); -void SMPcolSwap( SMPmatrix * , int , int ); -void SMPdestroy( SMPmatrix *); -int SMPfillup( SMPmatrix * ); -SMPelement * SMPfindElt( SMPmatrix *, int , int , int ); -void SMPgetError( SMPmatrix *, int *, int *); -int SMPluFac( SMPmatrix *, double , double ); double * SMPmakeElt( SMPmatrix * , int , int ); +void SMPcClear( SMPmatrix *); +void SMPclear( SMPmatrix *); +int SMPcLUfac( SMPmatrix *, double ); +int SMPluFac( SMPmatrix *, double , double ); +int SMPcReorder( SMPmatrix * , double , double , int *); +int SMPreorder( SMPmatrix * , double , double , double ); +void SMPcaSolve(SMPmatrix *Matrix, double RHS[], double iRHS[], + double Spare[], double iSpare[]); +void SMPcSolve( SMPmatrix *, double [], double [], double [], double []); +void SMPsolve( SMPmatrix *, double [], double []); int SMPmatSize( SMPmatrix *); int SMPnewMatrix( SMPmatrix ** ); -int SMPnewNode( int , SMPmatrix *); +void SMPdestroy( SMPmatrix *); int SMPpreOrder( SMPmatrix *); void SMPprint( SMPmatrix * , FILE *); -int SMPreorder( SMPmatrix * , double , double , double ); -void SMProwSwap( SMPmatrix * , int , int ); -void SMPsolve( SMPmatrix *, double [], double []); -#else /* stdc */ -int SMPaddElt(); -void SMPcClear(); -int SMPcLUfac(); -int SMPcProdDiag(); -int SMPcReorder(); -void SMPcSolve(); -void SMPclear(); -void SMPcolSwap(); -void SMPdestroy(); -int SMPfillup(); -SMPelement * SMPfindElt(); -void SMPgetError(); -int SMPluFac(); -double * SMPmakeElt(); -int SMPmatSize(); -int SMPnewMatrix(); -int SMPnewNode(); -int SMPpreOrder(); -void SMPprint(); -int SMPreorder(); -void SMProwSwap(); -void SMPsolve(); -#endif /* stdc */ +void SMPgetError( SMPmatrix *, int *, int *); +int SMPcProdDiag( SMPmatrix *, SPcomplex *, int *); +int SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent); +SMPelement * SMPfindElt( SMPmatrix *, int , int , int ); +int SMPcZeroCol(SMPmatrix *eMatrix, int Col); +int SMPcAddCol(SMPmatrix *eMatrix, int Accum_Col, int Addend_Col); +int SMPzeroRow(SMPmatrix *eMatrix, int Row); +#ifdef PARALLEL_ARCH +void SMPcombine(SMPmatrix *Matrix, double RHS[], double Spare[]); +void SMPcCombine(SMPmatrix *Matrix, double RHS[], double Spare[], + double iRHS[], double iSpare[]); +#endif #endif /*SMP*/ diff --git a/src/include/spmatrix.h b/src/include/spmatrix.h index 1ea24fb50..225264284 100644 --- a/src/include/spmatrix.h +++ b/src/include/spmatrix.h @@ -42,17 +42,6 @@ #ifndef spOKAY -/* - * IMPORTS - * - * >>> Import descriptions: - * spConfig.h - * Macros that customize the sparse matrix routines. - */ - -#include "spconfig.h" - - @@ -76,8 +65,7 @@ * Fatal error. A zero was encountered on the diagonal the matrix. This * does not necessarily imply that the matrix is singular. When this * error occurs, the matrix should be reconstructed and factored using - * spOrderAndFactor(). In spCOMPATIBILITY mode, spZERO_DIAG is - * indistinguishable from spSINGULAR. + * spOrderAndFactor(). * spSINGULAR * Fatal error. Matrix is singular, so no unique solution exists. * spNO_MEMORY @@ -108,21 +96,6 @@ #define spFATAL E_BADMATRIX -#if spCOMPATIBILITY -#define NO_ERROR spOKAY -#define UNDER_FLOW spOKAY -#define OVER_FLOW spOKAY -#define ILL_CONDITIONED spSMALL_PIVOT -#define SINGULAR spSINGULAR -#define NO_MEMORY spNO_MEMORY -#define RANGE spPANIC - -#define FATAL spFATAL - -#undef spZERO_DIAG -#define spZERO_DIAG spSINGULAR -#endif /* spCOMPATIBILITY */ - @@ -145,10 +118,6 @@ #define spREAL double /* #define void int */ -#if spCOMPATIBILITY -#define SPARSE_REAL spREAL -#endif - /* @@ -242,15 +211,6 @@ *((template).Element4Negated+1) -= imag; \ } -#if spCOMPATIBILITY -#define ADD_REAL_ELEMENT_TO_MATRIX spADD_REAL_ELEMENT -#define ADD_IMAG_ELEMENT_TO_MATRIX spADD_IMAG_ELEMENT -#define ADD_COMPLEX_ELEMENT_TO_MATRIX spADD_COMPLEX_ELEMENT -#define ADD_REAL_QUAD_ELEMENT_TO_MATRIX spADD_REAL_ELEMENT -#define ADD_IMAG_QUAD_ELEMENT_TO_MATRIX spADD_IMAG_ELEMENT -#define ADD_COMPLEX_QUAD_ELEMENT_TO_MATRIX spADD_COMPLEX_ELEMENT -#endif - @@ -272,10 +232,6 @@ * before being placed in Element3 and Element4. */ -#if spCOMPATIBILITY -#define spTemplate TemplateStruct -#endif - /* Begin `spTemplate'. */ struct spTemplate { spREAL *Element1 ; @@ -296,97 +252,47 @@ struct spTemplate /* Begin function declarations. */ -extern void spClear( char* ); -extern spREAL spCondition( char*, spREAL, int* ); -extern char *spCreate( int, int, int* ); -extern void spDeleteRowAndCol( char*, int, int ); -extern void spDestroy( char* ); -extern int spElementCount( char* ); -extern int spError( char* ); -extern int spFactor( char* ); -extern int spFileMatrix( char*, char*, char*, int, int, int ); -extern int spFileStats( char*, char*, char* ); -extern int spFillinCount( char* ); -extern int spGetAdmittance( char*, int, int, struct spTemplate* ); -extern spREAL *spGetElement( char*, int, int ); +extern void spClear( void * ); +extern spREAL spCondition( void *, spREAL, int* ); +extern void *spCreate( int, int, int* ); +extern void spDeleteRowAndCol( void *, int, int ); +extern void spDestroy(void *); +extern int spElementCount( void * ); +extern int spError( void * ); +extern int spFactor( void * ); +extern int spFileMatrix( void *, char *, char *, int, int, int ); +extern int spFileStats( void *, char *, char * ); +extern int spFillinCount( void * ); +extern int spGetAdmittance( void *, int, int, struct spTemplate* ); +extern spREAL *spGetElement( void *, int, int ); extern char *spGetInitInfo( spREAL* ); -extern int spGetOnes( char*, int, int, int, struct spTemplate* ); -extern int spGetQuad( char*, int, int, int, int, struct spTemplate* ); -extern int spGetSize( char*, int ); -extern int spInitialize( char*, int (*)() ); -extern void spInstallInitInfo( spREAL*, char* ); -extern spREAL spLargestElement( char* ); -extern void spMNA_Preorder( char* ); -extern spREAL spNorm( char* ); -extern int spOrderAndFactor( char*, spREAL*, spREAL, spREAL, int ); -extern int spOriginalCount( char *); -extern void spPartition( char*, int ); -extern void spPrint( char*, int, int, int ); -extern spREAL spPseudoCondition( char* ); -extern spREAL spRoundoff( char*, spREAL ); -extern void spScale( char*, spREAL*, spREAL* ); -extern void spSetComplex( char* ); -extern void spSetReal( char* ); -extern void spStripFills( char* ); -extern void spWhereSingular( char*, int*, int* ); +extern int spGetOnes( void *, int, int, int, struct spTemplate* ); +extern int spGetQuad( void *, int, int, int, int, struct spTemplate* ); +extern int spGetSize( void *, int ); +extern int spInitialize( void *, int (*)() ); +extern void spInstallInitInfo( spREAL*, void * ); +extern spREAL spLargestElement( void * ); +extern void spMNA_Preorder( void * ); +extern spREAL spNorm( void * ); +extern int spOrderAndFactor( void *, spREAL*, spREAL, spREAL, int ); +extern int spOriginalCount( void *); +extern void spPartition( void *, int ); +extern void spPrint( void *, int, int, int ); +extern spREAL spPseudoCondition( void * ); +extern spREAL spRoundoff( void *, spREAL ); +extern void spScale( void *, spREAL*, spREAL* ); +extern void spSetComplex( void * ); +extern void spSetReal( void * ); +extern void spStripFills( void * ); +extern void spWhereSingular( void *, int*, int* ); /* Functions with argument lists that are dependent on options. */ -#if spCOMPLEX -extern void spDeterminant ( char*, int*, spREAL*, spREAL* ); -#else /* NOT spCOMPLEX */ -extern void spDeterminant ( char*, int*, spREAL* ); -#endif /* NOT spCOMPLEX */ -#if spCOMPLEX && spSEPARATED_COMPLEX_VECTORS -extern int spFileVector( char*, char* , spREAL*, spREAL*); -extern void spMultiply( char*, spREAL*, spREAL*, spREAL*, spREAL* ); -extern void spMultTransposed(char*,spREAL*,spREAL*,spREAL*,spREAL*); -extern void spSolve( char*, spREAL*, spREAL*, spREAL*, spREAL* ); -extern void spSolveTransposed(char*,spREAL*,spREAL*,spREAL*,spREAL*); -#else /* NOT (spCOMPLEX && spSEPARATED_COMPLEX_VECTORS) */ -extern int spFileVector( char*, char* , spREAL* ); -extern void spMultiply( char*, spREAL*, spREAL* ); -extern void spMultTransposed( char*, spREAL*, spREAL* ); -extern void spSolve( char*, spREAL*, spREAL* ); -extern void spSolveTransposed( char*, spREAL*, spREAL* ); -#endif /* NOT (spCOMPLEX && spSEPARATED_COMPLEX_VECTORS) */ - -#if spCOMPATIBILITY -extern char *AllocateMatrix(); -extern spREAL *AddElementToMatrix(); -extern void AddRealElementToMatrix(); -extern void AddImagElementToMatrix(); -extern void AddComplexElementToMatrix(); -extern void AddAdmittanceToMatrix(); -extern void AddOnesToMatrix(); -extern void AddQuadToMatrix(); -extern void AddRealQuadElementToMatrix(); -extern void AddImagQuadElementToMatrix(); -extern void AddComplexQuadElementToMatrix(); -extern void CleanMatrix(); -extern void ClearMatrix(); -extern int ClearMatrixError(); -extern void DeallocateMatrix(); -extern void DeleteRowAndColFromMatrix(); -extern void Determinant(); -extern int DecomposeMatrix(); -extern int GetMatrixSize(); -extern int MatrixElementCount(); -extern int MatrixFillinCount(); -extern void MatrixMultiply(); -extern spREAL MatrixRoundoffError(); -extern int MatrixError(); -extern int OrderAndDecomposeMatrix(); -extern void OutputMatrixToFile(); -extern void OutputStatisticsToFile(); -extern void OutputVectorToFile(); -extern void PreorderForModifiedNodal(); -extern void PrintMatrix(); -extern void SetMatrixComplex(); -extern void SetMatrixReal(); -extern void SolveMatrix(); -extern void SolveTransposedMatrix(); -extern void ScaleMatrix(); -#endif /* spCOMPATIBILITY */ +extern void spDeterminant ( void *, int*, spREAL*, spREAL* ); +extern int spFileVector( void *, char * , spREAL*, spREAL*); +extern void spMultiply( void *, spREAL*, spREAL*, spREAL*, spREAL* ); +extern void spMultTransposed(void *,spREAL*,spREAL*,spREAL*,spREAL*); +extern void spSolve( void *, spREAL*, spREAL*, spREAL*, spREAL* ); +extern void spSolveTransposed(void *,spREAL*,spREAL*,spREAL*,spREAL*); #endif /* spOKAY */ diff --git a/src/maths/sparse/Makefile.am b/src/maths/sparse/Makefile.am index a9fb22371..ddedb9686 100644 --- a/src/maths/sparse/Makefile.am +++ b/src/maths/sparse/Makefile.am @@ -3,16 +3,17 @@ noinst_LIBRARIES = libsparse.a libsparse_a_SOURCES = \ - spalloc.c \ - spbuild.c \ - spcombin.c \ - spdefs.h \ - spextra.c \ - spfactor.c \ - spoutput.c \ - spsmp.c \ - spsolve.c \ - sputils.c + spalloc.c \ + spbuild.c \ + spcombin.c \ + spconfig.h \ + spdefs.h \ + spextra.c \ + spfactor.c \ + spoutput.c \ + spsmp.c \ + spsolve.c \ + sputils.c diff --git a/src/maths/sparse/spalloc.c b/src/maths/sparse/spalloc.c index 3f6827096..9011902f5 100644 --- a/src/maths/sparse/spalloc.c +++ b/src/maths/sparse/spalloc.c @@ -58,6 +58,7 @@ * spDefs.h * Matrix type and macro definitions for the sparse matrix routines. */ +#include #define spINSIDE_SPARSE #include "spconfig.h" @@ -75,15 +76,9 @@ * Function declarations */ -#ifdef __STDC__ static void InitializeElementBlocks( MatrixPtr, int, int ); -static void RecordAllocation( MatrixPtr, char* ); +static void RecordAllocation( MatrixPtr, void *); static void AllocateBlockOfAllocationList( MatrixPtr ); -#else /* __STDC__ */ -static void InitializeElementBlocks(); -static void RecordAllocation(); -static void AllocateBlockOfAllocationList(); -#endif /* __STDC__ */ @@ -124,51 +119,42 @@ static void AllocateBlockOfAllocationList(); * Error is cleared in this routine. */ -char * -spCreate( Size, Complex, pError ) - -int Size, *pError; -BOOLEAN Complex; +void * +spCreate(int Size, int Complex, int *pError) { -register unsigned SizePlusOne; -register MatrixPtr Matrix; -register int I; -int AllocatedSize; + unsigned SizePlusOne; + MatrixPtr Matrix; + int I; + int AllocatedSize; -/* Begin `spCreate'. */ -/* Clear error flag. */ + /* Begin `spCreate'. */ + /* Clear error flag. */ *pError = spOKAY; -/* Test for valid size. */ - if ((Size < 0) OR (Size == 0 AND NOT EXPANDABLE)) - { *pError = spPANIC; + /* Test for valid size. */ + if ((Size < 0) || (Size == 0 && !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; + /* Test for valid type. */ +#if !REAL + if (!Complex) { + *pError = spPANIC; + return NULL; } #endif -/* Create Matrix. */ + /* Create Matrix. */ AllocatedSize = MAX( Size, MINIMUM_ALLOCATED_SIZE ); SizePlusOne = (unsigned)(AllocatedSize + 1); - if ((Matrix = ALLOC(struct MatrixFrame, 1)) == NULL) - { *pError = spNO_MEMORY; - return NULL; + if ((Matrix = ALLOC(struct MatrixFrame, 1)) == NULL) { + *pError = spNO_MEMORY; + return NULL; } -/* Initialize matrix */ + /* Initialize matrix */ Matrix->ID = SPARSE_ID; Matrix->Complex = Complex; Matrix->PreviousMatrixWasComplex = Complex; @@ -208,14 +194,12 @@ int AllocatedSize; Matrix->ElementsRemaining = 0; Matrix->FillinsRemaining = 0; - RecordAllocation( Matrix, (char *)Matrix ); + RecordAllocation( Matrix, (void *)Matrix ); if (Matrix->Error == spNO_MEMORY) goto MemoryError; -/* Take out the trash. */ + /* 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; @@ -224,67 +208,68 @@ int AllocatedSize; Matrix->TrashCan.pInitInfo = NULL; #endif -/* Allocate space in memory for Diag pointer vector. */ + /* 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. */ + /* 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. */ + /* 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. */ + /* Allocate space in memory for IntToExtColMap vector. */ if (( Matrix->IntToExtColMap = ALLOC(int, SizePlusOne)) == NULL) goto MemoryError; -/* Allocate space in memory for IntToExtRowMap vector. */ + /* Allocate space in memory for IntToExtRowMap vector. */ if (( Matrix->IntToExtRowMap = ALLOC(int, SizePlusOne)) == NULL) goto MemoryError; -/* Initialize MapIntToExt vectors. */ + /* Initialize MapIntToExt vectors. */ for (I = 1; I <= AllocatedSize; I++) - { Matrix->IntToExtRowMap[I] = I; + { + Matrix->IntToExtRowMap[I] = I; Matrix->IntToExtColMap[I] = I; } #if TRANSLATE -/* Allocate space in memory for ExtToIntColMap vector. */ + /* Allocate space in memory for ExtToIntColMap vector. */ if (( Matrix->ExtToIntColMap = ALLOC(int, SizePlusOne)) == NULL) goto MemoryError; -/* Allocate space in memory for ExtToIntRowMap vector. */ + /* 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; + /* 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. */ + /* Allocate space for fill-ins and initial set of elements. */ InitializeElementBlocks( Matrix, SPACE_FOR_ELEMENTS*AllocatedSize, - SPACE_FOR_FILL_INS*AllocatedSize ); + SPACE_FOR_FILL_INS*AllocatedSize ); if (Matrix->Error == spNO_MEMORY) goto MemoryError; - return (char *)Matrix; + return (void *)Matrix; -MemoryError: + MemoryError: -/* Deallocate matrix and return no pointer to matrix if there is not enough - memory. */ + /* Deallocate matrix and return no pointer to matrix if there is not enough + memory. */ *pError = spNO_MEMORY; - spDestroy( (char *)Matrix); + spDestroy( (void *)Matrix); return NULL; } @@ -320,48 +305,46 @@ MemoryError: */ ElementPtr -spcGetElement( Matrix ) - -MatrixPtr Matrix; +spcGetElement(MatrixPtr Matrix) { -ElementPtr pElements; + ElementPtr pElements; -/* Begin `spcGetElement'. */ + /* Begin `spcGetElement'. */ -#if NOT COMBINE OR STRIP OR LINT -/* Allocate block of MatrixElements if necessary. */ - if (Matrix->ElementsRemaining == 0) - { pElements = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION); - RecordAllocation( Matrix, (char *)pElements ); +#if !COMBINE || STRIP || LINT + /* Allocate block of MatrixElements if necessary. */ + if (Matrix->ElementsRemaining == 0) { + pElements = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION); + RecordAllocation( Matrix, (void *)pElements ); if (Matrix->Error == spNO_MEMORY) return NULL; Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION; Matrix->NextAvailElement = pElements; } #endif -#if COMBINE OR STRIP OR LINT +#if COMBINE || STRIP || LINT if (Matrix->ElementsRemaining == 0) - { pListNode = Matrix->LastElementListNode; + { + pListNode = Matrix->LastElementListNode; -/* First see if there are any stripped elements left. */ - if (pListNode->Next != NULL) - { Matrix->LastElementListNode = pListNode = pListNode->Next; + /* First see if there are any stripped elements left. */ + if (pListNode->Next != NULL) { + Matrix->LastElementListNode = pListNode = pListNode->Next; Matrix->ElementsRemaining = pListNode->NumberOfElementsInList; Matrix->NextAvailElement = pListNode->pElementList; - } - else - { -/* Allocate block of elements. */ + } else { + /* Allocate block of elements. */ pElements = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION); - RecordAllocation( Matrix, (char *)pElements ); + RecordAllocation( Matrix, (void *)pElements ); if (Matrix->Error == spNO_MEMORY) return NULL; Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION; Matrix->NextAvailElement = pElements; -/* Allocate an element list structure. */ + /* Allocate an element list structure. */ pListNode->Next = ALLOC(struct ElementListNodeStruct,1); - RecordAllocation( Matrix, (char *)pListNode->Next ); - if (Matrix->Error == spNO_MEMORY) return NULL; + RecordAllocation( Matrix, (void *)pListNode->Next ); + if (Matrix->Error == spNO_MEMORY) + return NULL; Matrix->LastElementListNode = pListNode = pListNode->Next; pListNode->pElementList = pElements; @@ -371,7 +354,7 @@ ElementPtr pElements; } #endif -/* Update Element counter and return pointer to Element. */ + /* Update Element counter and return pointer to Element. */ Matrix->ElementsRemaining--; return Matrix->NextAvailElement++; @@ -387,11 +370,11 @@ ElementPtr pElements; /* * 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. + * 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) @@ -411,30 +394,26 @@ ElementPtr pElements; * A pointer to the first element in the group of elements being allocated. * * >>> Possible errors: - * spNO_MEMORY - */ + * spNO_MEMORY */ static void -InitializeElementBlocks( Matrix, InitialNumberOfElements, - NumberOfFillinsExpected ) - -MatrixPtr Matrix; -int InitialNumberOfElements, NumberOfFillinsExpected; +InitializeElementBlocks(MatrixPtr Matrix, int InitialNumberOfElements, + int NumberOfFillinsExpected) { -ElementPtr pElement; + ElementPtr pElement; -/* Begin `InitializeElementBlocks'. */ + /* Begin `InitializeElementBlocks'. */ -/* Allocate block of MatrixElements for elements. */ + /* Allocate block of MatrixElements for elements. */ pElement = ALLOC(struct MatrixElement, InitialNumberOfElements); - RecordAllocation( Matrix, (char *)pElement ); + RecordAllocation( Matrix, (void *)pElement ); if (Matrix->Error == spNO_MEMORY) return; Matrix->ElementsRemaining = InitialNumberOfElements; Matrix->NextAvailElement = pElement; -/* Allocate an element list structure. */ + /* Allocate an element list structure. */ Matrix->FirstElementListNode = ALLOC(struct ElementListNodeStruct,1); - RecordAllocation( Matrix, (char *)Matrix->FirstElementListNode ); + RecordAllocation( Matrix, (void *)Matrix->FirstElementListNode ); if (Matrix->Error == spNO_MEMORY) return; Matrix->LastElementListNode = Matrix->FirstElementListNode; @@ -443,16 +422,16 @@ ElementPtr pElement; InitialNumberOfElements; Matrix->FirstElementListNode->Next = NULL; -/* Allocate block of MatrixElements for fill-ins. */ + /* Allocate block of MatrixElements for fill-ins. */ pElement = ALLOC(struct MatrixElement, NumberOfFillinsExpected); - RecordAllocation( Matrix, (char *)pElement ); + RecordAllocation( Matrix, (void *)pElement ); if (Matrix->Error == spNO_MEMORY) return; Matrix->FillinsRemaining = NumberOfFillinsExpected; Matrix->NextAvailFillin = pElement; -/* Allocate a fill-in list structure. */ + /* Allocate a fill-in list structure. */ Matrix->FirstFillinListNode = ALLOC(struct FillinListNodeStruct,1); - RecordAllocation( Matrix, (char *)Matrix->FirstFillinListNode ); + RecordAllocation( Matrix, (void *)Matrix->FirstFillinListNode ); if (Matrix->Error == spNO_MEMORY) return; Matrix->LastFillinListNode = Matrix->FirstFillinListNode; @@ -475,10 +454,10 @@ ElementPtr pElement; /* * 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. + * 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. @@ -488,43 +467,38 @@ ElementPtr pElement; * Pointer to matrix. * * >>> Possible errors: - * spNO_MEMORY - */ + * spNO_MEMORY */ ElementPtr -spcGetFillin( Matrix ) - -MatrixPtr Matrix; +spcGetFillin(MatrixPtr Matrix) { -/* Begin `spcGetFillin'. */ + /* Begin `spcGetFillin'. */ -#if NOT STRIP OR LINT +#if !STRIP || LINT if (Matrix->FillinsRemaining == 0) return spcGetElement( Matrix ); #endif -#if STRIP OR LINT +#if STRIP || LINT - if (Matrix->FillinsRemaining == 0) - { pListNode = Matrix->LastFillinListNode; + 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; + /* 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. */ + } else { + /* Allocate block of fill-ins. */ pFillins = ALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION); - RecordAllocation( Matrix, (char *)pFillins ); + RecordAllocation( Matrix, (void *)pFillins ); if (Matrix->Error == spNO_MEMORY) return NULL; Matrix->FillinsRemaining = ELEMENTS_PER_ALLOCATION; Matrix->NextAvailFillin = pFillins; -/* Allocate a fill-in list structure. */ + /* Allocate a fill-in list structure. */ pListNode->Next = ALLOC(struct FillinListNodeStruct,1); - RecordAllocation( Matrix, (char *)pListNode->Next ); + RecordAllocation( Matrix, (void *)pListNode->Next ); if (Matrix->Error == spNO_MEMORY) return NULL; Matrix->LastFillinListNode = pListNode = pListNode->Next; @@ -535,7 +509,7 @@ MatrixPtr Matrix; } #endif -/* Update Fill-in counter and return pointer to Fill-in. */ + /* Update Fill-in counter and return pointer to Fill-in. */ Matrix->FillinsRemaining--; return Matrix->NextAvailFillin++; } @@ -551,50 +525,43 @@ MatrixPtr Matrix; /* * RECORD A MEMORY ALLOCATION * - * This routine is used to record all memory allocations so that the memory - * can be freed later. + * 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. + * AllocatedPtr (void *) + * 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 - */ + * spNO_MEMORY */ static void -RecordAllocation( Matrix, AllocatedPtr ) - -MatrixPtr Matrix; -char *AllocatedPtr; +RecordAllocation(MatrixPtr Matrix, void *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; + /* 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; } -/* Allocate block of MatrixElements if necessary. */ - if (Matrix->RecordsRemaining == 0) - { AllocateBlockOfAllocationList( Matrix ); - if (Matrix->Error == spNO_MEMORY) - { FREE(AllocatedPtr); + /* Allocate block of MatrixElements if necessary. */ + if (Matrix->RecordsRemaining == 0) { + AllocateBlockOfAllocationList( Matrix ); + if (Matrix->Error == spNO_MEMORY) { + FREE(AllocatedPtr); return; } } -/* Add Allocated pointer to Allocation List. */ + /* Add Allocated pointer to Allocation List. */ (++Matrix->TopOfAllocationList)->AllocatedPtr = AllocatedPtr; Matrix->RecordsRemaining--; return; - } @@ -615,42 +582,41 @@ char *AllocatedPtr; * * >>> 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. + * 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 - */ + * spNO_MEMORY */ static void -AllocateBlockOfAllocationList( Matrix ) - -MatrixPtr Matrix; +AllocateBlockOfAllocationList(MatrixPtr Matrix) { -register int I; -register AllocationListPtr ListPtr; + int I; + AllocationListPtr ListPtr; -/* Begin `AllocateBlockOfAllocationList'. */ -/* Allocate block of records for allocation list. */ + /* Begin `AllocateBlockOfAllocationList'. */ + /* Allocate block of records for allocation list. */ ListPtr = ALLOC(struct AllocationRecord, (ELEMENTS_PER_ALLOCATION+1)); - if (ListPtr == NULL) - { Matrix->Error = spNO_MEMORY; + if (ListPtr == NULL) { + Matrix->Error = spNO_MEMORY; return; } -/* String entries of allocation list into singly linked list. List is linked - such that any record points to the one before it. */ + /* 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--; + 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; + /* Record allocation of space for allocation list on allocation list. */ + Matrix->TopOfAllocationList->AllocatedPtr = (void *)ListPtr; Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION; return; @@ -669,7 +635,7 @@ register AllocationListPtr ListPtr; * Deallocates pointers and elements of Matrix. * * >>> Arguments: - * Matrix (char *) + * Matrix (void *) * Pointer to the matrix frame which is to be removed from memory. * * >>> Local variables: @@ -684,18 +650,16 @@ register AllocationListPtr ListPtr; */ void -spDestroy( eMatrix ) - -register char *eMatrix; +spDestroy(void *eMatrix) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register AllocationListPtr ListPtr, NextListPtr; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + AllocationListPtr ListPtr, NextListPtr; -/* Begin `spDestroy'. */ - ASSERT( IS_SPARSE( Matrix ) ); + /* Begin `spDestroy'. */ + assert( IS_SPARSE( Matrix ) ); -/* Deallocate the vectors that are located in the matrix frame. */ + /* Deallocate the vectors that are located in the matrix frame. */ FREE( Matrix->IntToExtColMap ); FREE( Matrix->IntToExtRowMap ); FREE( Matrix->ExtToIntColMap ); @@ -710,17 +674,14 @@ register AllocationListPtr ListPtr, NextListPtr; FREE( Matrix->DoRealDirect ); FREE( Matrix->Intermediate ); -/* Sequentially step through the list of allocated pointers freeing pointers - * along the way. */ + /* Sequentially step through the list of allocated pointers + * freeing pointers along the way. */ ListPtr = Matrix->TopOfAllocationList; - while (ListPtr != NULL) - { NextListPtr = ListPtr->NextRecord; - if ((char *) ListPtr == ListPtr->AllocatedPtr) - { + while (ListPtr != NULL) { + NextListPtr = ListPtr->NextRecord; + if ((void *) ListPtr == ListPtr->AllocatedPtr) { FREE( ListPtr ); - } - else - { + } else { FREE( ListPtr->AllocatedPtr ); } ListPtr = NextListPtr; @@ -737,29 +698,27 @@ register AllocationListPtr ListPtr, NextListPtr; /* * RETURN MATRIX ERROR STATUS * - * This function is used to determine the error status of the given matrix. + * 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. - */ - + * eMatrix (void *) + * The matrix for which the error status is desired. */ int -spError( eMatrix ) - -char *eMatrix; +spError(void *eMatrix ) { -/* Begin `spError'. */ + /* Begin `spError'. */ - if (eMatrix != NULL) - { ASSERT(((MatrixPtr)eMatrix)->ID == SPARSE_ID); + if (eMatrix != NULL) { + assert(((MatrixPtr)eMatrix)->ID == SPARSE_ID); return ((MatrixPtr)eMatrix)->Error; + } else { + /* This error may actually be spPANIC, no way to tell. */ + return spNO_MEMORY; } - else return spNO_MEMORY; /* This error may actually be spPANIC, - * no way to tell. */ } @@ -777,7 +736,7 @@ char *eMatrix; * detected as singular or where a zero was detected on the diagonal. * * >>> Arguments: - * eMatrix (char *) + * eMatrix (void *) * The matrix for which the error status is desired. * pRow (int *) * The row number. @@ -786,18 +745,16 @@ char *eMatrix; */ void -spWhereSingular( eMatrix, pRow, pCol ) - -char *eMatrix; -int *pRow, *pCol; +spWhereSingular(void *eMatrix, int *pRow, int *pCol) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; + MatrixPtr Matrix = (MatrixPtr)eMatrix; -/* Begin `spWhereSingular'. */ - ASSERT( IS_SPARSE( Matrix ) ); + /* Begin `spWhereSingular'. */ + assert( IS_SPARSE( Matrix ) ); - if (Matrix->Error == spSINGULAR OR Matrix->Error == spZERO_DIAG) - { *pRow = Matrix->SingularRow; + if (Matrix->Error == spSINGULAR || Matrix->Error == spZERO_DIAG) + { + *pRow = Matrix->SingularRow; *pCol = Matrix->SingularCol; } else *pRow = *pCol = 0; @@ -816,9 +773,9 @@ MatrixPtr Matrix = (MatrixPtr)eMatrix; * the matrix is returned. * * >>> Arguments: - * eMatrix (char *) + * eMatrix (void *) * Pointer to matrix. - * External (BOOLEAN) + * External (int) * 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 @@ -826,15 +783,12 @@ MatrixPtr Matrix = (MatrixPtr)eMatrix; */ int -spGetSize( eMatrix, External ) - -char *eMatrix; -BOOLEAN External; +spGetSize(void *eMatrix, int External) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; + MatrixPtr Matrix = (MatrixPtr)eMatrix; -/* Begin `spGetSize'. */ - ASSERT( IS_SPARSE( Matrix ) ); + /* Begin `spGetSize'. */ + assert( IS_SPARSE( Matrix ) ); #if TRANSLATE if (External) @@ -859,31 +813,27 @@ MatrixPtr Matrix = (MatrixPtr)eMatrix; * Forces matrix to be either real or complex. * * >>> Arguments: - * eMatrix (char *) + * eMatrix (void *) * Pointer to matrix. */ void -spSetReal( eMatrix ) - -char *eMatrix; +spSetReal(void *eMatrix) { -/* Begin `spSetReal'. */ + /* Begin `spSetReal'. */ - ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) AND REAL); + assert( IS_SPARSE( (MatrixPtr)eMatrix )); ((MatrixPtr)eMatrix)->Complex = NO; return; } void -spSetComplex( eMatrix ) - -char *eMatrix; +spSetComplex(void *eMatrix) { -/* Begin `spSetComplex'. */ + /* Begin `spSetComplex'. */ - ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) AND spCOMPLEX); + assert( IS_SPARSE( (MatrixPtr)eMatrix )); ((MatrixPtr)eMatrix)->Complex = YES; return; } @@ -904,40 +854,34 @@ char *eMatrix; * of original elements can be returned. * * >>> Arguments: - * eMatrix (char *) + * eMatrix (void *) * Pointer to matrix. */ int -spFillinCount( eMatrix ) - -char *eMatrix; +spFillinCount(void *eMatrix) { -/* Begin `spFillinCount'. */ + /* Begin `spFillinCount'. */ - ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) ); + assert( IS_SPARSE( (MatrixPtr)eMatrix ) ); return ((MatrixPtr)eMatrix)->Fillins; } int -spElementCount( eMatrix ) - -char *eMatrix; +spElementCount(void *eMatrix) { -/* Begin `spElementCount'. */ + /* Begin `spElementCount'. */ - ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) ); + assert( IS_SPARSE( (MatrixPtr)eMatrix ) ); return ((MatrixPtr)eMatrix)->Elements; } int -spOriginalCount( eMatrix ) - -char *eMatrix; +spOriginalCount(void *eMatrix) { -/* Begin `spOriginalCount'. */ + /* Begin `spOriginalCount'. */ - ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) ); + assert( IS_SPARSE( (MatrixPtr)eMatrix ) ); return ((MatrixPtr)eMatrix)->Originals; } diff --git a/src/maths/sparse/spbuild.c b/src/maths/sparse/spbuild.c index 3fd64db41..07d3e9c69 100644 --- a/src/maths/sparse/spbuild.c +++ b/src/maths/sparse/spbuild.c @@ -56,6 +56,7 @@ * spDefs.h * Matrix type and macro definitions for the sparse matrix routines. */ +#include #define spINSIDE_SPARSE #include "spconfig.h" @@ -69,17 +70,9 @@ /* * Function declarations */ - -#ifdef __STDC__ static void Translate( MatrixPtr, int*, int* ); static void EnlargeMatrix( MatrixPtr, int ); static void ExpandTranslationArrays( MatrixPtr, int ); -#else /* __STDC__ */ -static void Translate(); -static void EnlargeMatrix(); -static void ExpandTranslationArrays(); -#endif /* __STDC__ */ - @@ -100,45 +93,45 @@ static void ExpandTranslationArrays(); */ void -spClear( eMatrix ) - -char *eMatrix; +spClear(void *eMatrix ) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register ElementPtr pElement; -register int I; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + ElementPtr pElement; + int I; -/* Begin `spClear'. */ - ASSERT( IS_SPARSE( Matrix ) ); + /* 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; - } - } + /* Clear matrix. */ + if (Matrix->PreviousMatrixWasComplex || 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; - } - } + { + for (I = Matrix->Size; I > 0; I--) + { + pElement = Matrix->FirstInCol[I]; + while (pElement != NULL) + { + pElement->Real = 0.0; + pElement = pElement->NextInCol; + } + } } -/* Empty the trash. */ + /* Empty the trash. */ Matrix->TrashCan.Real = 0.0; -#if spCOMPLEX Matrix->TrashCan.Imag = 0.0; -#endif Matrix->Error = spOKAY; Matrix->Factored = NO; @@ -194,24 +187,21 @@ register int I; */ RealNumber * -spGetElement( eMatrix, Row, Col ) - -char *eMatrix; -int Row, Col; +spGetElement(void *eMatrix, int Row, int Col) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -RealNumber *pElement; -ElementPtr spcFindElementInCol(); -void Translate(); + MatrixPtr Matrix = (MatrixPtr)eMatrix; + RealNumber *pElement; + ElementPtr spcFindElementInCol(); + void Translate(); -/* Begin `spGetElement'. */ - ASSERT( IS_SPARSE( Matrix ) AND Row >= 0 AND Col >= 0 ); + /* Begin `spGetElement'. */ + assert( IS_SPARSE( Matrix ) && Row >= 0 && Col >= 0 ); - if ((Row == 0) OR (Col == 0)) - return &Matrix->TrashCan.Real; + if ((Row == 0) || (Col == 0)) + return &Matrix->TrashCan.Real; -#if NOT TRANSLATE - ASSERT(Matrix->NeedsOrdering); +#if !TRANSLATE + assert(Matrix->NeedsOrdering); #endif #if TRANSLATE @@ -219,40 +209,38 @@ void Translate(); if (Matrix->Error == spNO_MEMORY) return NULL; #endif -#if NOT TRANSLATE -#if NOT EXPANDABLE - ASSERT(Row <= Matrix->Size AND Col <= Matrix->Size); +#if !TRANSLATE +#if !EXPANDABLE + assert(Row <= Matrix->Size && Col <= Matrix->Size); #endif #if EXPANDABLE -/* Re-size Matrix if necessary. */ - if ((Row > Matrix->Size) OR (Col > Matrix->Size)) - EnlargeMatrix( Matrix, MAX(Row, Col) ); + /* Re-size Matrix if necessary. */ + if ((Row > Matrix->Size) || (Col > Matrix->Size)) + EnlargeMatrix( Matrix, MAX(Row, Col) ); if (Matrix->Error == spNO_MEMORY) return NULL; #endif #endif -/* - * The condition part of the following if statement tests to see if the - * element resides along the diagonal, if it does then it tests to see - * if the element has been created yet (Diag pointer not NULL). The - * pointer to the element is then assigned to Element after it is cast - * into a pointer to a RealNumber. This casting makes the pointer into - * a pointer to Real. This statement depends on the fact that Real - * is the first record in the MatrixElement structure. - */ + /* The condition part of the following if statement tests to see + * if the element resides along the diagonal, if it does then it + * tests to see if the element has been created yet (Diag pointer + * not NULL). The pointer to the element is then assigned to + * Element after it is cast into a pointer to a RealNumber. This + * casting makes the pointer into a pointer to Real. This + * statement depends on the fact that Real is the first record in + * the MatrixElement structure. */ - if ((Row != Col) OR ((pElement = (RealNumber *)Matrix->Diag[Row]) == NULL)) + if ((Row != Col) || ((pElement = (RealNumber *)Matrix->Diag[Row]) == NULL)) { -/* - * Element does not exist or does not reside along diagonal. Search - * column for element. As in the if statement above, the pointer to the - * element which is returned by spcFindElementInCol is cast into a - * pointer to Real, a RealNumber. - */ - pElement = (RealNumber*)spcFindElementInCol( Matrix, - &(Matrix->FirstInCol[Col]), - Row, Col, YES ); + /* Element does not exist or does not reside along diagonal. + * Search column for element. As in the if statement above, + * the pointer to the element which is returned by + * spcFindElementInCol is cast into a pointer to Real, a + * RealNumber. */ + pElement = (RealNumber*)spcFindElementInCol( Matrix, + &(Matrix->FirstInCol[Col]), + Row, Col, YES ); } return pElement; } @@ -290,7 +278,7 @@ void Translate(); * Row being searched for. * Col (int) * Column being searched. - * CreateIfMissing (BOOLEAN) + * CreateIfMissing (int) * Indicates what to do if element is not found, create one or return a * NULL pointer. * @@ -300,40 +288,37 @@ void Translate(); */ ElementPtr -spcFindElementInCol( Matrix, LastAddr, Row, Col, CreateIfMissing ) - -MatrixPtr Matrix; -register ElementPtr *LastAddr; -register int Row; -int Col; -BOOLEAN CreateIfMissing; +spcFindElementInCol(MatrixPtr Matrix, ElementPtr *LastAddr, + int Row, int Col, int CreateIfMissing) { -register ElementPtr pElement; -ElementPtr spcCreateElement(); + ElementPtr pElement; + ElementPtr spcCreateElement(); -/* Begin `spcFindElementInCol'. */ + /* Begin `spcFindElementInCol'. */ pElement = *LastAddr; -/* Search for element. */ + /* Search for element. */ while (pElement != NULL) - { if (pElement->Row < Row) + { + if (pElement->Row < Row) { -/* Have not reached element yet. */ - LastAddr = &(pElement->NextInCol); - pElement = pElement->NextInCol; + /* Have not reached element yet. */ + LastAddr = &(pElement->NextInCol); + pElement = pElement->NextInCol; } - else if (pElement->Row == Row) + else if (pElement->Row == Row) { -/* Reached element. */ - return pElement; + /* Reached element. */ + return pElement; } - else break; /* while loop */ + else break; /* while loop */ } -/* Element does not exist and must be created. */ + /* Element does not exist and must be created. */ if (CreateIfMissing) - return spcCreateElement( Matrix, Row, Col, LastAddr, NO ); - else return NULL; + return spcCreateElement( Matrix, Row, Col, LastAddr, NO ); + else + return NULL; } @@ -378,41 +363,39 @@ ElementPtr spcCreateElement(); */ static void -Translate( Matrix, Row, Col ) - -MatrixPtr Matrix; -int *Row, *Col; +Translate(MatrixPtr Matrix, int *Row, int *Col) { -register int IntRow, IntCol, ExtRow, ExtCol; + int IntRow, IntCol, ExtRow, ExtCol; -/* Begin `Translate'. */ + /* Begin `Translate'. */ ExtRow = *Row; ExtCol = *Col; -/* Expand translation arrays if necessary. */ - if ((ExtRow > Matrix->AllocatedExtSize) OR + /* Expand translation arrays if necessary. */ + if ((ExtRow > Matrix->AllocatedExtSize) || (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)) + /* Set ExtSize if necessary. */ + if ((ExtRow > Matrix->ExtSize) || (ExtCol > Matrix->ExtSize)) Matrix->ExtSize = MAX(ExtRow, ExtCol); -/* Translate external row or node number to internal row or node number. */ + /* Translate external row or node number to internal row or node number. */ if ((IntRow = Matrix->ExtToIntRowMap[ExtRow]) == -1) - { Matrix->ExtToIntRowMap[ExtRow] = ++Matrix->CurrentSize; + { + Matrix->ExtToIntRowMap[ExtRow] = ++Matrix->CurrentSize; Matrix->ExtToIntColMap[ExtRow] = Matrix->CurrentSize; IntRow = Matrix->CurrentSize; -#if NOT EXPANDABLE - ASSERT(IntRow <= Matrix->Size); +#if !EXPANDABLE + assert(IntRow <= Matrix->Size); #endif #if EXPANDABLE -/* Re-size Matrix if necessary. */ + /* Re-size Matrix if necessary. */ if (IntRow > Matrix->Size) EnlargeMatrix( Matrix, IntRow ); if (Matrix->Error == spNO_MEMORY) return; @@ -422,18 +405,19 @@ register int IntRow, IntCol, ExtRow, ExtCol; Matrix->IntToExtColMap[IntRow] = ExtRow; } -/* Translate external column or node number to internal column or node number.*/ + /* Translate external column or node number to internal column or node number.*/ if ((IntCol = Matrix->ExtToIntColMap[ExtCol]) == -1) - { Matrix->ExtToIntRowMap[ExtCol] = ++Matrix->CurrentSize; + { + Matrix->ExtToIntRowMap[ExtCol] = ++Matrix->CurrentSize; Matrix->ExtToIntColMap[ExtCol] = Matrix->CurrentSize; IntCol = Matrix->CurrentSize; -#if NOT EXPANDABLE - ASSERT(IntCol <= Matrix->Size); +#if !EXPANDABLE + assert(IntCol <= Matrix->Size); #endif #if EXPANDABLE -/* Re-size Matrix if necessary. */ + /* Re-size Matrix if necessary. */ if (IntCol > Matrix->Size) EnlargeMatrix( Matrix, IntCol ); if (Matrix->Error == spNO_MEMORY) return; @@ -490,24 +474,19 @@ register int IntRow, IntCol, ExtRow, ExtCol; */ int -spGetAdmittance( Matrix, Node1, Node2, Template ) - -char *Matrix; -int Node1, Node2; -struct spTemplate *Template; +spGetAdmittance(void *Matrix, int Node1, int Node2, + struct spTemplate *Template) { - -/* Begin `spGetAdmittance'. */ + /* Begin `spGetAdmittance'. */ Template->Element1 = spGetElement(Matrix, Node1, Node1 ); Template->Element2 = spGetElement(Matrix, Node2, Node2 ); Template->Element3Negated = spGetElement( Matrix, Node2, Node1 ); Template->Element4Negated = spGetElement( Matrix, Node1, Node2 ); - if - ( (Template->Element1 == NULL) - OR (Template->Element2 == NULL) - OR (Template->Element3Negated == NULL) - OR (Template->Element4Negated == NULL) - ) return spNO_MEMORY; + if ((Template->Element1 == NULL) || + (Template->Element2 == NULL) || + (Template->Element3Negated == NULL) || + (Template->Element4Negated == NULL)) + return spNO_MEMORY; if (Node1 == 0) SWAP( RealNumber*, Template->Element1, Template->Element2 ); @@ -578,23 +557,19 @@ struct spTemplate *Template; */ int -spGetQuad( Matrix, Row1, Row2, Col1, Col2, Template ) - -char *Matrix; -int Row1, Row2, Col1, Col2; -struct spTemplate *Template; +spGetQuad(void *Matrix, int Row1, int Row2, int Col1, int Col2, + struct spTemplate *Template) { -/* Begin `spGetQuad'. */ + /* Begin `spGetQuad'. */ Template->Element1 = spGetElement( Matrix, Row1, Col1); Template->Element2 = spGetElement( Matrix, Row2, Col2 ); Template->Element3Negated = spGetElement( Matrix, Row2, Col1 ); Template->Element4Negated = spGetElement( Matrix, Row1, Col2 ); - if - ( (Template->Element1 == NULL) - OR (Template->Element2 == NULL) - OR (Template->Element3Negated == NULL) - OR (Template->Element4Negated == NULL) - ) return spNO_MEMORY; + if ((Template->Element1 == NULL) || + (Template->Element2 == NULL) || + (Template->Element3Negated == NULL) || + (Template->Element4Negated == NULL)) + return spNO_MEMORY; if (Template->Element1 == &((MatrixPtr)Matrix)->TrashCan.Real) SWAP( RealNumber *, Template->Element1, Template->Element2 ); @@ -653,23 +628,19 @@ struct spTemplate *Template; */ int -spGetOnes(Matrix, Pos, Neg, Eqn, Template) - -char *Matrix; -int Pos, Neg, Eqn; -struct spTemplate *Template; +spGetOnes(void *Matrix, int Pos, int Neg, int Eqn, + struct spTemplate *Template) { -/* Begin `spGetOnes'. */ + /* Begin `spGetOnes'. */ Template->Element4Negated = spGetElement( Matrix, Neg, Eqn ); Template->Element3Negated = spGetElement( Matrix, Eqn, Neg ); Template->Element2 = spGetElement( Matrix, Pos, Eqn ); Template->Element1 = spGetElement( Matrix, Eqn, Pos ); - if - ( (Template->Element1 == NULL) - OR (Template->Element2 == NULL) - OR (Template->Element3Negated == NULL) - OR (Template->Element4Negated == NULL) - ) return spNO_MEMORY; + if ((Template->Element1 == NULL) || + (Template->Element2 == NULL) || + (Template->Element3Negated == NULL) || + (Template->Element4Negated == NULL)) + return spNO_MEMORY; spADD_REAL_QUAD( *Template, 1.0 ); return spOKAY; @@ -703,7 +674,7 @@ struct spTemplate *Template; * This contains the address of the pointer to the element just above the * one being created. It is used to speed the search and it is updated with * address of the created element. - * Fillin (BOOLEAN) + * Fillin (int) * Flag that indicates if created element is to be a fill-in. * * >>> Local variables: @@ -723,78 +694,73 @@ struct spTemplate *Template; */ ElementPtr -spcCreateElement( Matrix, Row, Col, LastAddr, Fillin ) - -MatrixPtr Matrix; -int Row; -register int Col; -register ElementPtr *LastAddr; -BOOLEAN Fillin; +spcCreateElement(MatrixPtr Matrix, int Row, int Col, + ElementPtr *LastAddr, int Fillin) { -register ElementPtr pElement, pLastElement; -ElementPtr pCreatedElement, spcGetElement(), spcGetFillin(); + ElementPtr pElement, pLastElement; + ElementPtr pCreatedElement, spcGetElement(), spcGetFillin(); -/* Begin `spcCreateElement'. */ + /* Begin `spcCreateElement'. */ if (Matrix->RowsLinked) { -/* Row pointers cannot be ignored. */ + /* Row pointers cannot be ignored. */ if (Fillin) - { pElement = spcGetFillin( Matrix ); + { + pElement = spcGetFillin( Matrix ); Matrix->Fillins++; } else - { pElement = spcGetElement( Matrix ); + { + pElement = spcGetElement( Matrix ); Matrix->Originals++; Matrix->NeedsOrdering = YES; } if (pElement == NULL) return NULL; -/* If element is on diagonal, store pointer in Diag. */ + /* If element is on diagonal, store pointer in Diag. */ if (Row == Col) Matrix->Diag[Row] = pElement; -/* Initialize Element. */ + /* Initialize Element. */ pCreatedElement = pElement; pElement->Row = Row; pElement->Col = Col; pElement->Real = 0.0; -#if spCOMPLEX pElement->Imag = 0.0; -#endif #if INITIALIZE pElement->pInitInfo = NULL; #endif -/* Splice element into column. */ + /* Splice element into column. */ pElement->NextInCol = *LastAddr; *LastAddr = pElement; - /* Search row for proper element position. */ + /* Search row for proper element position. */ pElement = Matrix->FirstInRow[Row]; pLastElement = NULL; while (pElement != NULL) { -/* Search for element row position. */ + /* Search for element row position. */ if (pElement->Col < Col) { -/* Have not reached desired element. */ + /* Have not reached desired element. */ pLastElement = pElement; pElement = pElement->NextInRow; } else pElement = NULL; } -/* Splice element into row. */ + /* Splice element into row. */ pElement = pCreatedElement; if (pLastElement == NULL) { -/* Element is first in row. */ + /* Element is first in row. */ pElement->NextInRow = Matrix->FirstInRow[Row]; Matrix->FirstInRow[Row] = pElement; } else -/* Element is not first in row. */ { + /* Element is not first in row. */ pElement->NextInRow = pLastElement->NextInRow; pLastElement->NextInRow = pElement; } @@ -802,34 +768,30 @@ ElementPtr pCreatedElement, spcGetElement(), spcGetFillin(); } else { -/* - * Matrix has not been factored yet. Thus get element rather than fill-in. - * Also, row pointers can be ignored. - */ + /* Matrix has not been factored yet. Thus get element rather + * than fill-in. Also, row pointers can be ignored. */ -/* Allocate memory for Element. */ + /* Allocate memory for Element. */ pElement = spcGetElement( Matrix ); Matrix->Originals++; if (pElement == NULL) return NULL; -/* If element is on diagonal, store pointer in Diag. */ + /* If element is on diagonal, store pointer in Diag. */ if (Row == Col) Matrix->Diag[Row] = pElement; -/* Initialize Element. */ + /* Initialize Element. */ pCreatedElement = pElement; pElement->Row = Row; #if DEBUG pElement->Col = Col; #endif pElement->Real = 0.0; -#if spCOMPLEX pElement->Imag = 0.0; -#endif #if INITIALIZE pElement->pInitInfo = NULL; #endif -/* Splice element into column. */ + /* Splice element into column. */ pElement->NextInCol = *LastAddr; *LastAddr = pElement; } @@ -865,30 +827,29 @@ ElementPtr pCreatedElement, spcGetElement(), spcGetFillin(); * currently being operated upon. * FirstInRowArray (ArrayOfElementPtrs) * A pointer to the FirstInRow array. Same as Matrix->FirstInRow but - * resides in a register and requires less indirection so is faster to + * resides in a and requires less indirection so is faster to * use. * Col (int) * Column currently being operated upon. */ void -spcLinkRows( Matrix ) - -MatrixPtr Matrix; +spcLinkRows(MatrixPtr Matrix) { -register ElementPtr pElement, *FirstInRowEntry; -register ArrayOfElementPtrs FirstInRowArray; -register int Col; + ElementPtr pElement, *FirstInRowEntry; + ArrayOfElementPtrs FirstInRowArray; + int Col; -/* Begin `spcLinkRows'. */ + /* Begin `spcLinkRows'. */ FirstInRowArray = Matrix->FirstInRow; for (Col = Matrix->Size; Col >= 1; Col--) { -/* Generate row links for the elements in the Col'th column. */ + /* Generate row links for the elements in the Col'th column. */ pElement = Matrix->FirstInCol[Col]; while (pElement != NULL) - { pElement->Col = Col; + { + pElement->Col = Col; FirstInRowEntry = &FirstInRowArray[pElement->Row]; pElement->NextInRow = *FirstInRowEntry; *FirstInRowEntry = pElement; @@ -923,48 +884,48 @@ register int Col; */ static void -EnlargeMatrix( Matrix, NewSize ) - -MatrixPtr Matrix; -register int NewSize; +EnlargeMatrix(MatrixPtr Matrix, int NewSize) { -register int I, OldAllocatedSize = Matrix->AllocatedSize; + int I, OldAllocatedSize = Matrix->AllocatedSize; -/* Begin `EnlargeMatrix'. */ + /* Begin `EnlargeMatrix'. */ Matrix->Size = NewSize; if (NewSize <= OldAllocatedSize) return; -/* Expand the matrix frame. */ + /* Expand the matrix frame. */ NewSize = MAX( NewSize, EXPANSION_FACTOR * OldAllocatedSize ); Matrix->AllocatedSize = NewSize; if (( REALLOC(Matrix->IntToExtColMap, int, NewSize+1)) == NULL) - { Matrix->Error = spNO_MEMORY; + { + Matrix->Error = spNO_MEMORY; return; } if (( REALLOC(Matrix->IntToExtRowMap, int, NewSize+1)) == NULL) - { Matrix->Error = spNO_MEMORY; + { + Matrix->Error = spNO_MEMORY; return; } if (( REALLOC(Matrix->Diag, ElementPtr, NewSize+1)) == NULL) - { Matrix->Error = spNO_MEMORY; + { + Matrix->Error = spNO_MEMORY; return; } if (( REALLOC(Matrix->FirstInCol, ElementPtr, NewSize+1)) == NULL) - { Matrix->Error = spNO_MEMORY; + { + Matrix->Error = spNO_MEMORY; return; } if (( REALLOC(Matrix->FirstInRow, ElementPtr, NewSize+1)) == NULL) - { Matrix->Error = spNO_MEMORY; + { + Matrix->Error = spNO_MEMORY; return; } -/* - * Destroy the Markowitz and Intermediate vectors, they will be recreated - * in spOrderAndFactor(). - */ + /* Destroy the Markowitz and Intermediate vectors, they will be + * recreated in spOrderAndFactor(). */ FREE( Matrix->MarkowitzRow ); FREE( Matrix->MarkowitzCol ); FREE( Matrix->MarkowitzProd ); @@ -973,9 +934,10 @@ register int I, OldAllocatedSize = Matrix->AllocatedSize; FREE( Matrix->Intermediate ); Matrix->InternalVectorsAllocated = NO; -/* Initialize the new portion of the vectors. */ + /* Initialize the new portion of the vectors. */ for (I = OldAllocatedSize+1; I <= NewSize; I++) - { Matrix->IntToExtColMap[I] = I; + { + Matrix->IntToExtColMap[I] = I; Matrix->IntToExtRowMap[I] = I; Matrix->Diag[I] = NULL; Matrix->FirstInRow[I] = NULL; @@ -1012,35 +974,35 @@ register int I, OldAllocatedSize = Matrix->AllocatedSize; */ static void -ExpandTranslationArrays( Matrix, NewSize ) - -MatrixPtr Matrix; -register int NewSize; +ExpandTranslationArrays(MatrixPtr Matrix, int NewSize) { -register int I, OldAllocatedSize = Matrix->AllocatedExtSize; + int I, OldAllocatedSize = Matrix->AllocatedExtSize; -/* Begin `ExpandTranslationArrays'. */ + /* Begin `ExpandTranslationArrays'. */ Matrix->ExtSize = NewSize; if (NewSize <= OldAllocatedSize) return; -/* Expand the translation arrays ExtToIntRowMap and ExtToIntColMap. */ + /* Expand the translation arrays ExtToIntRowMap and ExtToIntColMap. */ NewSize = MAX( NewSize, EXPANSION_FACTOR * OldAllocatedSize ); Matrix->AllocatedExtSize = NewSize; if (( REALLOC(Matrix->ExtToIntRowMap, int, NewSize+1)) == NULL) - { Matrix->Error = spNO_MEMORY; + { + Matrix->Error = spNO_MEMORY; return; } if (( REALLOC(Matrix->ExtToIntColMap, int, NewSize+1)) == NULL) - { Matrix->Error = spNO_MEMORY; + { + Matrix->Error = spNO_MEMORY; return; } -/* Initialize the new portion of the vectors. */ + /* Initialize the new portion of the vectors. */ for (I = OldAllocatedSize+1; I <= NewSize; I++) - { Matrix->ExtToIntRowMap[I] = -1; + { + Matrix->ExtToIntRowMap[I] = -1; Matrix->ExtToIntColMap[I] = -1; } @@ -1090,72 +1052,68 @@ register int I, OldAllocatedSize = Matrix->AllocatedExtSize; */ void -spInstallInitInfo( pElement, pInitInfo ) - -RealNumber *pElement; -char *pInitInfo; +spInstallInitInfo(RealNumber *pElement, char *pInitInfo) { -/* Begin `spInstallInitInfo'. */ - ASSERT(pElement != NULL); + /* Begin `spInstallInitInfo'. */ + assert(pElement != NULL); ((ElementPtr)pElement)->pInitInfo = pInitInfo; } char * -spGetInitInfo( pElement ) - -RealNumber *pElement; +spGetInitInfo(RealNumber *pElement) { -/* Begin `spGetInitInfo'. */ - ASSERT(pElement != NULL); + /* Begin `spGetInitInfo'. */ + assert(pElement != NULL); return (char *)((ElementPtr)pElement)->pInitInfo; } int -spInitialize( eMatrix, pInit ) - -char *eMatrix; -int (*pInit)(); +spInitialize(void *eMatrix, int (*pInit)()) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register ElementPtr pElement; -int J, Error, Col; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + ElementPtr pElement; + int J, Error, Col; -/* Begin `spInitialize'. */ - ASSERT( IS_SPARSE( Matrix ) ); + /* Begin `spInitialize'. */ + assert( IS_SPARSE( Matrix ) ); -#if spCOMPLEX -/* Clear imaginary part of matrix if matrix is real but was complex. */ - if (Matrix->PreviousMatrixWasComplex AND NOT Matrix->Complex) - { for (J = Matrix->Size; J > 0; J--) - { pElement = Matrix->FirstInCol[J]; + /* Clear imaginary part of matrix if matrix is real but was complex. */ + if (Matrix->PreviousMatrixWasComplex && !Matrix->Complex) + { + for (J = Matrix->Size; J > 0; J--) + { + pElement = Matrix->FirstInCol[J]; while (pElement != NULL) - { pElement->Imag = 0.0; + { + pElement->Imag = 0.0; pElement = pElement->NextInCol; } } } -#endif /* spCOMPLEX */ -/* Initialize the matrix. */ + /* Initialize the matrix. */ for (J = Matrix->Size; J > 0; J--) - { pElement = Matrix->FirstInCol[J]; + { + pElement = Matrix->FirstInCol[J]; Col = Matrix->IntToExtColMap[J]; while (pElement != NULL) - { if (pElement->pInitInfo == NULL) - { pElement->Real = 0.0; -# if spCOMPLEX - pElement->Imag = 0.0; -# endif + { + if (pElement->pInitInfo == NULL) + { + pElement->Real = 0.0; + pElement->Imag = 0.0; } else - { Error = (*pInit)((RealNumber *)pElement, pElement->pInitInfo, + { + Error = (*pInit)((RealNumber *)pElement, pElement->pInitInfo, Matrix->IntToExtRowMap[pElement->Row], Col); if (Error) - { Matrix->Error = spFATAL; + { + Matrix->Error = spFATAL; return Error; } @@ -1164,11 +1122,9 @@ int J, Error, Col; } } -/* Empty the trash. */ + /* Empty the trash. */ Matrix->TrashCan.Real = 0.0; -#if spCOMPLEX Matrix->TrashCan.Imag = 0.0; -#endif Matrix->Error = spOKAY; Matrix->Factored = NO; diff --git a/src/maths/sparse/spcombin.c b/src/maths/sparse/spcombin.c index fe70cec65..9679509ee 100644 --- a/src/maths/sparse/spcombin.c +++ b/src/maths/sparse/spcombin.c @@ -58,20 +58,10 @@ #if COMBINE -#ifdef __STDC__ -#if spSEPARATED_COMPLEX_VECTORS static void CombineComplexMatrix( MatrixPtr, RealVector, RealVector, RealVector, RealVector ); -#else -static void CombineComplexMatrix( MatrixPtr, RealVector, RealVector ); -#endif static void ClearBuffer( MatrixPtr, int, int, ElementPtr ); static void ClearComplexBuffer( MatrixPtr, int, int, ElementPtr ); -#else /* __STDC__ */ -static void CombineComplexMatrix(); -static void ClearBuffer(); -static void ClearComplexBuffer(); -#endif /* __STDC__ */ /* * COMBINE MATRICES ON A MULTIPROCESSOR @@ -91,10 +81,10 @@ static void ClearComplexBuffer(); static double Buffer[SPBSIZE]; void -spCombine( eMatrix, RHS, Spare IMAG_VECTORS ) +spCombine( eMatrix, RHS, Spare, iRHS, iSolution ) char *eMatrix; -RealVector RHS, Spare IMAG_VECTORS; +RealVector RHS, Spare, iRHS, iSolution; { MatrixPtr Matrix = (MatrixPtr)eMatrix; register ElementPtr pElement; @@ -105,18 +95,15 @@ struct ElementListNodeStruct *pListNode; long type = MT_COMBINE, length = Matrix->Size + 1; /* Begin `spCombine'. */ - ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored ); - if (NOT Matrix->InternalVectorsAllocated) + assert( IS_VALID(Matrix) && !Matrix->Factored ); + if (!Matrix->InternalVectorsAllocated) spcCreateInternalVectors( Matrix ); -#if spCOMPLEX if (Matrix->Complex) { - CombineComplexMatrix( Matrix, RHS, Spare IMAG_VECTORS ); + CombineComplexMatrix( Matrix, RHS, Spare, iRHS, iSolution ); return; } -#endif -#if REAL Size = Matrix->Size; /* Mark original non-zeroes. */ @@ -167,15 +154,14 @@ long type = MT_COMBINE, length = Matrix->Size + 1; DGOP_( &type, RHS, &length, "+" ); return; -#endif /* REAL */ } -#if spCOMPLEX + static void -CombineComplexMatrix( Matrix, RHS, Spare IMAG_VECTORS ) +CombineComplexMatrix( Matrix, RHS, Spare, iRHS, iSolution ) MatrixPtr Matrix; -RealVector RHS, Spare IMAG_VECTORS; +RealVector RHS, Spare, iRHS, iSolution; { register ElementPtr pElement; register int I, Size; @@ -185,7 +171,7 @@ struct ElementListNodeStruct *pListNode; long type = MT_COMBINE, length = Matrix->Size + 1; /* Begin `CombineComplexMatrix'. */ - ASSERT(Matrix->Complex); + assert(Matrix->Complex); Size = Matrix->Size; /* Mark original non-zeroes. */ @@ -235,19 +221,13 @@ long type = MT_COMBINE, length = Matrix->Size + 1; } /* Sum all RHS's together */ -#if spSEPARATED_COMPLEX_VECTORS DGOP_( &type, RHS, &length, "+" ); DGOP_( &type, iRHS, &length, "+" ); -#else - length *= 2; - DGOP_( &type, RHS, &length, "+" ); -#endif return; } -#endif /* spCOMPLEX */ -#if REAL + static void ClearBuffer( Matrix, NumElems, StartCol, StartElement ) @@ -281,9 +261,8 @@ ElementPtr StartElement; pElement = pElement->NextInCol; } } -#endif REAL -#if spCOMPLEX + static void ClearComplexBuffer( Matrix, DataCount, StartCol, StartElement ) @@ -318,5 +297,4 @@ ElementPtr StartElement; pElement = pElement->NextInCol; } } -#endif /* spCOMPLEX */ #endif /* COMBINE */ diff --git a/src/include/spconfig.h b/src/maths/sparse/spconfig.h similarity index 86% rename from src/include/spconfig.h rename to src/maths/sparse/spconfig.h index 17ecf3dea..71470e5ba 100644 --- a/src/include/spconfig.h +++ b/src/maths/sparse/spconfig.h @@ -62,15 +62,6 @@ * less storage is required. There is often a noticeable speed * penalty when using single precision. Sparse internally refers * to a spREAL as a RealNumber. - * REAL - * This specifies that the routines are expected to handle real - * systems of equations. The routines can be compiled to handle - * both real and complex systems at the same time, but there is a - * slight speed and memory advantage if the routines are complied - * to handle only real systems of equations. - * spCOMPLEX - * This specifies that the routines will be complied to handle - * complex systems of equations. * EXPANDABLE * Setting this compiler flag true (1) makes the matrix * expandable before it has been factored. If the matrix is @@ -129,28 +120,6 @@ * threshold. When DIAGONAL_PIVOTING is turned off, the following * options and constants are not used: MODIFIED_MARKOWITZ, * MAX_MARKOWITZ_TIES, and TIES_MULTIPLIER. - * ARRAY_OFFSET - * This determines whether arrays start at an index of zero or one. - * This option is necessitated by the fact that standard C - * convention dictates that arrays begin with an index of zero but - * the standard mathematic convention states that arrays begin with - * an index of one. So if you prefer to start your arrays with - * zero, or your calling Sparse from FORTRAN, set ARRAY_OFFSET to - * NO or 0. Otherwise, set ARRAY_OFFSET to YES or 1. Note that if - * you use an offset of one, the arrays that you pass to Sparse - * must have an allocated length of one plus the size of the - * matrix. ARRAY_OFFSET must be either 0 or 1, no other offsets - * are valid. - * spSEPARATED_COMPLEX_VECTORS - * This specifies the format for complex vectors. If this is set - * false then a complex vector is made up of one double sized - * array of RealNumber's in which the real and imaginary numbers - * are placed in the alternately array in the array. In other - * words, the first entry would be Complex[1].Real, then comes - * Complex[1].Imag, then Complex[1].Real, etc. If - * spSEPARATED_COMPLEX_VECTORS is set true, then each complex - * vector is represented by two arrays of RealNumbers, one with - * the real terms, the other with the imaginary. [NO] * MODIFIED_MARKOWITZ * This specifies that the modified Markowitz method of pivot * selection is to be used. The modified Markowitz method differs @@ -229,33 +198,19 @@ * MULTIPLICATION * This specifies that the routines to multiply the unfactored * matrix by a vector should be compiled. - * FORTRAN - * This specifies that the FORTRAN interface routines should be - * compiled. When interfacing to FORTRAN programs, the ARRAY_OFFSET - * options should be set to NO. * DEBUG * This specifies that additional error checking will be compiled. * The type of error checked are those that are common when the * matrix routines are first integrated into a user's program. Once * the routines have been integrated in and are running smoothly, this * option should be turned off. - * spCOMPATIBILITY - * This specifies that Sparse1.3 should be configured to be upward - * compatible from Sparse1.2. This option is not suggested for use - * in new software. Sparse1.3, when configured to be compatible with - * Sparse1.2, is not completely compatible. In particular, if - * recompiling the calling program, it is necessary to change the - * of the Sparse include files. This option will go away in future - * versions of Sparse. [0] */ /* Begin options. */ -#define REAL YES #define EXPANDABLE YES #define TRANSLATE YES #define INITIALIZE NO #define DIAGONAL_PIVOTING YES -#define ARRAY_OFFSET YES #define MODIFIED_MARKOWITZ NO #define DELETE NO #define STRIP NO @@ -270,7 +225,6 @@ #define STABILITY NO #define CONDITION NO #define PSEUDOCONDITION NO -#define FORTRAN NO #ifdef HAS_MINDATA # define DEBUG NO #else @@ -283,14 +237,6 @@ * constants YES an NO are not defined in spMatrix.h to avoid conflicts * with user code, so use 0 for NO and 1 for YES. */ -#endif /* spINSIDE_SPARSE */ -#define spCOMPLEX 1 -#define spSEPARATED_COMPLEX_VECTORS 1 -#define spCOMPATIBILITY 0 -#ifdef spINSIDE_SPARSE - - - diff --git a/src/maths/sparse/spdefs.h b/src/maths/sparse/spdefs.h index 333f32277..bb865d728 100644 --- a/src/maths/sparse/spdefs.h +++ b/src/maths/sparse/spdefs.h @@ -1,3 +1,5 @@ +#ifndef _SPDEFS_H +#define _SPDEFS_H /* * DATA STRUCTURE AND MACRO DEFINITIONS for Sparse. * @@ -36,96 +38,29 @@ * IMPORTS */ +#include #include -#include "ngspice.h" #undef ABORT #undef MALLOC #undef FREE #undef REALLOC -/* - * If running lint, change some of the compiler options to get a more - * complete inspection. - */ - -#ifdef lint -#undef REAL -#undef spCOMPLEX -#undef EXPANDABLE -#undef TRANSLATE -#undef INITIALIZE -#undef DELETE -#undef STRIP -#undef MODIFIED_NODAL -#undef QUAD_ELEMENT -#undef TRANSPOSE -#undef SCALING -#undef DOCUMENTATION -#undef MULTIPLICATION -#undef DETERMINANT -#undef CONDITION -#undef PSEUDOCONDITION -#undef FORTRAN -#undef DEBUG -#undef spCOMPATIBILITY - - - -#define REAL YES -#define spCOMPLEX YES -#define EXPANDABLE YES -#define TRANSLATE YES -#define INITIALIZE YES -#define DELETE YES -#define STRIP YES -#define MODIFIED_NODAL YES -#define QUAD_ELEMENT YES -#define TRANSPOSE YES -#define SCALING YES -#define DOCUMENTATION YES -#define MULTIPLICATION YES -#define DETERMINANT YES -#define CONDITION YES -#define PSEUDOCONDITION YES -#define FORTRAN YES -#define DEBUG YES -#define spCOMPATIBILITY YES - -#define LINT YES -#else /* not lint */ -#define LINT NO -#endif /* not lint */ - - - - - /* * MACRO DEFINITIONS * * Macros are distinguished by using solely capital letters in their - * identifiers. This contrasts with C defined identifiers which are strictly - * lower case, and program variable and procedure names which use both upper - * and lower case. - */ + * identifiers. This contrasts with C defined identifiers which are + * strictly lower case, and program variable and procedure names + * which use both upper and lower case. */ /* Begin macros. */ -/* Boolean data type */ -#define BOOLEAN int #define NO 0 #define YES 1 -#define NOT ! -#define AND && -#define OR || -/* NULL pointer */ -#ifndef NULL -#define NULL 0 -#endif #define SPARSE_ID 0x772773 /* Arbitrary (is Sparse on phone). */ #define IS_SPARSE(matrix) ((matrix) != NULL && \ @@ -151,11 +86,7 @@ #define SWAP(type, a, b) {type swapx; swapx = a; a = b; b = swapx;} /* Macro function that returns the approx absolute value of a complex number. */ -#if spCOMPLEX #define ELEMENT_MAG(ptr) (ABS((ptr)->Real) + ABS((ptr)->Imag)) -#else -#define ELEMENT_MAG(ptr) ((ptr)->Real < 0.0 ? -(ptr)->Real : (ptr)->Real) -#endif /* Complex assignment statements. */ #define CMPLX_ASSIGN(to,from) \ @@ -317,114 +248,6 @@ (from_a).Imag * (from_b).Real; \ } -/* - * Macro functions that provide complex division. - */ - -/* Complex division: to = num / den */ -#define CMPLX_DIV(to,num,den) \ -{ RealNumber r_, s_; \ - if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \ - ((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \ - { r_ = (den).Imag / (den).Real; \ - s_ = (den).Real + r_*(den).Imag; \ - (to).Real = ((num).Real + r_*(num).Imag)/s_; \ - (to).Imag = ((num).Imag - r_*(num).Real)/s_; \ - } \ - else \ - { r_ = (den).Real / (den).Imag; \ - s_ = (den).Imag + r_*(den).Real; \ - (to).Real = (r_*(num).Real + (num).Imag)/s_; \ - (to).Imag = (r_*(num).Imag - (num).Real)/s_; \ - } \ -} - -/* Complex division and assignment: num /= den */ -#define CMPLX_DIV_ASSIGN(num,den) \ -{ RealNumber r_, s_, t_; \ - if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \ - ((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \ - { r_ = (den).Imag / (den).Real; \ - s_ = (den).Real + r_*(den).Imag; \ - t_ = ((num).Real + r_*(num).Imag)/s_; \ - (num).Imag = ((num).Imag - r_*(num).Real)/s_; \ - (num).Real = t_; \ - } \ - else \ - { r_ = (den).Real / (den).Imag; \ - s_ = (den).Imag + r_*(den).Real; \ - t_ = (r_*(num).Real + (num).Imag)/s_; \ - (num).Imag = (r_*(num).Imag - (num).Real)/s_; \ - (num).Real = t_; \ - } \ -} - -/* Complex reciprocation: to = 1.0 / den */ -#define CMPLX_RECIPROCAL(to,den) \ -{ RealNumber r_; \ - if (((den).Real >= (den).Imag AND (den).Real > -(den).Imag) OR \ - ((den).Real < (den).Imag AND (den).Real <= -(den).Imag)) \ - { r_ = (den).Imag / (den).Real; \ - (to).Imag = -r_*((to).Real = 1.0/((den).Real + r_*(den).Imag)); \ - } \ - else \ - { r_ = (den).Real / (den).Imag; \ - (to).Real = -r_*((to).Imag = -1.0/((den).Imag + r_*(den).Real));\ - } \ -} - - - - - - -/* - * ASSERT and ABORT - * - * Macro used to assert that if the code is working correctly, then - * a condition must be true. If not, then execution is terminated - * and an error message is issued stating that there is an internal - * error and giving the file and line number. These assertions are - * not evaluated unless the DEBUG flag is true. - */ - -#if DEBUG -#define ASSERT(condition) if (NOT(condition)) ABORT() -#else -#define ASSERT(condition) -#endif - -#if DEBUG -#define ABORT() \ -{ (void)fflush(stdout); \ - (void)fprintf(stderr, "sparse: panic in file `%s' at line %d.\n", \ - __FILE__, __LINE__); \ - (void)fflush(stderr); \ - abort(); \ -} -#else -#define ABORT() -#endif - - - - - -/* - * IMAGINARY VECTORS - * - * The imaginary vectors iRHS and iSolution are only needed when the - * options spCOMPLEX and spSEPARATED_COMPLEX_VECTORS are set. The following - * macro makes it easy to include or exclude these vectors as needed. - */ - -#if spCOMPLEX AND spSEPARATED_COMPLEX_VECTORS -#define IMAG_VECTORS , iRHS, iSolution -#define IMAG_RHS , iRHS -#else -#define IMAG_VECTORS -#define IMAG_RHS -#endif #define ALLOC(type,number) ((type *)tmalloc((unsigned)(sizeof(type)*(number)))) #define REALLOC(ptr,type,number) \ @@ -491,10 +314,9 @@ /* Begin `MatrixElement'. */ struct MatrixElement -{ RealNumber Real; -#if spCOMPLEX +{ + RealNumber Real; RealNumber Imag; -#endif int Row; int Col; struct MatrixElement *NextInRow; @@ -532,7 +354,8 @@ typedef ElementPtr *ArrayOfElementPtrs; /* Begin `AllocationRecord'. */ struct AllocationRecord -{ char *AllocatedPtr; +{ + char *AllocatedPtr; struct AllocationRecord *NextRecord; }; @@ -567,7 +390,8 @@ typedef struct AllocationRecord *AllocationListPtr; /* Begin `FillinListNodeStruct'. */ struct FillinListNodeStruct -{ ElementPtr pFillinList; +{ + ElementPtr pFillinList; int NumberOfFillinsInList; struct FillinListNodeStruct *Next; }; @@ -575,7 +399,8 @@ struct FillinListNodeStruct /* Similar to above, but keeps track of the original Elements */ /* Begin `ElementListNodeStruct'. */ struct ElementListNodeStruct -{ ElementPtr pElementList; +{ + ElementPtr pElementList; int NumberOfElementsInList; struct ElementListNodeStruct *Next; }; @@ -613,7 +438,7 @@ struct ElementListNodeStruct * grow to when EXPANDABLE is set true and AllocatedSize is the largest * the matrix can get without requiring that the matrix frame be * reallocated. - * Complex (BOOLEAN) + * Complex (int) * The flag which indicates whether the matrix is complex (true) or * real. * CurrentSize (int) @@ -622,12 +447,12 @@ struct ElementListNodeStruct * rows and columns that have elements in them. * Diag (ArrayOfElementPtrs) * Array of pointers that points to the diagonal elements. - * DoCmplxDirect (BOOLEAN *) + * DoCmplxDirect (int *) * Array of flags, one for each column in matrix. If a flag is true * then corresponding column in a complex matrix should be eliminated * in spFactor() using direct addressing (rather than indirect * addressing). - * DoRealDirect (BOOLEAN *) + * DoRealDirect (int *) * Array of flags, one for each column in matrix. If a flag is true * then corresponding column in a real matrix should be eliminated * in spFactor() using direct addressing (rather than indirect @@ -644,7 +469,7 @@ struct ElementListNodeStruct * ExtToIntRowMap (int []) * An array that is used to convert external row numbers to internal * external row numbers. Present only if TRANSLATE option is set true. - * Factored (BOOLEAN) + * Factored (int) * Indicates if matrix has been factored. This flag is set true in * spFactor() and spOrderAndFactor() and set false in spCreate() * and spClear(). @@ -665,7 +490,7 @@ struct ElementListNodeStruct * array used during forward and backward substitution. It is * commonly called y when the forward and backward substitution process is * denoted Ax = b => Ly = b and Ux = y. - * InternalVectorsAllocated (BOOLEAN) + * InternalVectorsAllocated (int) * A flag that indicates whether the Markowitz vectors and the * Intermediate vector have been created. * These vectors are created in spcCreateInternalVectors(). @@ -689,19 +514,19 @@ struct ElementListNodeStruct * The maximum number of off-diagonal element in the rows of L, the * lower triangular matrix. This quantity is used when computing an * estimate of the roundoff error in the matrix. - * NeedsOrdering (BOOLEAN) + * NeedsOrdering (int) * This is a flag that signifies that the matrix needs to be ordered * or reordered. NeedsOrdering is set true in spCreate() and * spGetElement() or spGetAdmittance() if new elements are added to the * matrix after it has been previously factored. It is set false in * spOrderAndFactor(). - * NumberOfInterchangesIsOdd (BOOLEAN) + * NumberOfInterchangesIsOdd (int) * Flag that indicates the sum of row and column interchange counts * is an odd number. Used when determining the sign of the determinant. * Originals (int) * The number of original elements (total elements minus fill ins) * present in matrix. - * Partitioned (BOOLEAN) + * Partitioned (int) * This flag indicates that the columns of the matrix have been * partitioned into two groups. Those that will be addressed directly * and those that will be addressed indirectly in spFactor(). @@ -711,7 +536,7 @@ struct ElementListNodeStruct * Row pivot was chosen from. * PivotSelectionMethod (char) * Character that indicates which pivot search method was successful. - * PreviousMatrixWasComplex (BOOLEAN) + * PreviousMatrixWasComplex (int) * This flag in needed to determine how to clear the matrix. When * dealing with real matrices, it is important that the imaginary terms * in the matrix elements be zero. Thus, if the previous matrix was @@ -720,11 +545,11 @@ struct ElementListNodeStruct * RelThreshold (RealNumber) * The magnitude an element must have relative to others in its row * to be considered as a pivot candidate, except as a last resort. - * Reordered (BOOLEAN) + * Reordered (int) * This flag signifies that the matrix has been reordered. It * is cleared in spCreate(), set in spMNA_Preorder() and * spOrderAndFactor() and is used in spPrint(). - * RowsLinked (BOOLEAN) + * RowsLinked (int) * A flag that indicates whether the row pointers exist. The AddByIndex * routines do not generate the row pointers, which are needed by some * of the other routines, such as spOrderAndFactor() and spScale(). @@ -781,43 +606,44 @@ struct ElementListNodeStruct /* Begin `MatrixFrame'. */ struct MatrixFrame -{ RealNumber AbsThreshold; +{ + RealNumber AbsThreshold; int AllocatedSize; int AllocatedExtSize; - BOOLEAN Complex; + int Complex; int CurrentSize; ArrayOfElementPtrs Diag; - BOOLEAN *DoCmplxDirect; - BOOLEAN *DoRealDirect; + int *DoCmplxDirect; + int *DoRealDirect; int Elements; int Error; int ExtSize; int *ExtToIntColMap; int *ExtToIntRowMap; - BOOLEAN Factored; + int Factored; int Fillins; ArrayOfElementPtrs FirstInCol; ArrayOfElementPtrs FirstInRow; unsigned long ID; RealVector Intermediate; - BOOLEAN InternalVectorsAllocated; + int InternalVectorsAllocated; int *IntToExtColMap; int *IntToExtRowMap; int *MarkowitzRow; int *MarkowitzCol; long *MarkowitzProd; int MaxRowCountInLowerTri; - BOOLEAN NeedsOrdering; - BOOLEAN NumberOfInterchangesIsOdd; + int NeedsOrdering; + int NumberOfInterchangesIsOdd; int Originals; - BOOLEAN Partitioned; + int Partitioned; int PivotsOriginalCol; int PivotsOriginalRow; char PivotSelectionMethod; - BOOLEAN PreviousMatrixWasComplex; + int PreviousMatrixWasComplex; RealNumber RelThreshold; - BOOLEAN Reordered; - BOOLEAN RowsLinked; + int Reordered; + int RowsLinked; int SingularCol; int SingularRow; int Singletons; @@ -844,7 +670,6 @@ typedef struct MatrixFrame *MatrixPtr; * Function declarations */ -#ifdef __STDC__ extern ElementPtr spcGetElement( MatrixPtr ); extern ElementPtr spcGetFillin( MatrixPtr ); extern ElementPtr spcFindElementInCol( MatrixPtr, ElementPtr*, int, int, int ); @@ -853,13 +678,6 @@ extern void spcCreateInternalVectors( MatrixPtr ); extern void spcLinkRows( MatrixPtr ); extern void spcColExchange( MatrixPtr, int, int ); extern void spcRowExchange( MatrixPtr, int, int ); -#else /* __STDC__ */ -extern ElementPtr spcGetElement(); -extern ElementPtr spcGetFillin(); -extern ElementPtr spcFindElementInCol(); -extern ElementPtr spcCreateElement(); -extern void spcCreateInternalVectors(); -extern void spcLinkRows(); -extern void spcColExchange(); -extern void spcRowExchange(); -#endif /* __STDC__ */ + +#endif + diff --git a/src/maths/sparse/spfactor.c b/src/maths/sparse/spfactor.c index e12b1526e..0bcf6a7ad 100644 --- a/src/maths/sparse/spfactor.c +++ b/src/maths/sparse/spfactor.c @@ -55,6 +55,7 @@ * spDefs.h * Matrix type and macro definitions for the sparse matrix routines. */ +#include #define spINSIDE_SPARSE #include "spconfig.h" @@ -155,7 +156,7 @@ static int ZeroPivot( MatrixPtr, int ); * pivot an element that has suffered heavy cancellation and as a * result mainly consists of roundoff error. Once a valid * threshold is given, it becomes the new default. - * DiagPivoting (BOOLEAN) + * DiagPivoting (int) * A flag indicating that pivot selection should be confined to the * diagonal if possible. If DiagPivoting is nonzero and if * DIAGONAL_PIVOTING is enabled pivots will be chosen only from @@ -177,7 +178,7 @@ static int ZeroPivot( MatrixPtr, int ); * >>> Local variables: * pPivot (ElementPtr) * Pointer to the element being used as a pivot. - * ReorderingRequired (BOOLEAN) + * ReorderingRequired (int) * Flag that indicates whether reordering is required. * * >>> Possible errors: @@ -188,86 +189,86 @@ static int ZeroPivot( MatrixPtr, int ); */ int -spOrderAndFactor( eMatrix, RHS, RelThreshold, AbsThreshold, DiagPivoting ) - -char *eMatrix; -RealNumber RHS[], RelThreshold, AbsThreshold; -BOOLEAN DiagPivoting; +spOrderAndFactor(void *eMatrix, RealNumber RHS[], RealNumber RelThreshold, + RealNumber AbsThreshold, int DiagPivoting) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -ElementPtr pPivot; -int Step, Size, ReorderingRequired; -ElementPtr SearchForPivot(); -RealNumber LargestInCol, FindLargestInCol(); + MatrixPtr Matrix = (MatrixPtr)eMatrix; + ElementPtr pPivot; + int Step, Size, ReorderingRequired; + ElementPtr SearchForPivot(); + RealNumber LargestInCol, FindLargestInCol(); -/* Begin `spOrderAndFactor'. */ - ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored); + /* Begin `spOrderAndFactor'. */ + assert( IS_VALID(Matrix) && !Matrix->Factored); Matrix->Error = spOKAY; Size = Matrix->Size; - if (RelThreshold <= 0.0) RelThreshold = Matrix->RelThreshold; - if (RelThreshold > 1.0) RelThreshold = Matrix->RelThreshold; + if (RelThreshold <= 0.0) + RelThreshold = Matrix->RelThreshold; + if (RelThreshold > 1.0) + RelThreshold = Matrix->RelThreshold; Matrix->RelThreshold = RelThreshold; - if (AbsThreshold < 0.0) AbsThreshold = Matrix->AbsThreshold; + if (AbsThreshold < 0.0) + AbsThreshold = Matrix->AbsThreshold; Matrix->AbsThreshold = AbsThreshold; ReorderingRequired = NO; - if (NOT Matrix->NeedsOrdering) + if (!Matrix->NeedsOrdering) { -/* Matrix has been factored before and reordering is not required. */ + /* Matrix has been factored before and reordering is not required. */ for (Step = 1; Step <= Size; Step++) - { pPivot = Matrix->Diag[Step]; + { + pPivot = Matrix->Diag[Step]; LargestInCol = FindLargestInCol(pPivot->NextInCol); if ((LargestInCol * RelThreshold < ELEMENT_MAG(pPivot))) - { if (Matrix->Complex) + { + if (Matrix->Complex) ComplexRowColElimination( Matrix, pPivot ); else RealRowColElimination( Matrix, pPivot ); } else - { ReorderingRequired = YES; + { + ReorderingRequired = YES; break; /* for loop */ } } - if (NOT ReorderingRequired) + if (!ReorderingRequired) goto Done; else { -/* - * A pivot was not large enough to maintain accuracy, - * so a partial reordering is required. - */ + /* A pivot was not large enough to maintain accuracy, so a + * partial reordering is required. */ #if (ANNOTATE >= ON_STRANGE_BEHAVIOR) printf("Reordering, Step = %1d\n", Step); #endif } - } /* End of if(NOT Matrix->NeedsOrdering) */ + } /* End of if(!Matrix->NeedsOrdering) */ else { -/* - * This is the first time the matrix has been factored. These few statements - * indicate to the rest of the code that a full reodering is required rather - * than a partial reordering, which occurs during a failure of a fast - * factorization. - */ + /* This is the first time the matrix has been factored. These + * few statements indicate to the rest of the code that a full + * reodering is required rather than a partial reordering, + * which occurs during a failure of a fast factorization. */ Step = 1; - if (NOT Matrix->RowsLinked) + if (!Matrix->RowsLinked) spcLinkRows( Matrix ); - if (NOT Matrix->InternalVectorsAllocated) + if (!Matrix->InternalVectorsAllocated) spcCreateInternalVectors( Matrix ); if (Matrix->Error >= spFATAL) return Matrix->Error; } -/* Form initial Markowitz products. */ + /* Form initial Markowitz products. */ CountMarkowitz( Matrix, RHS, Step ); MarkowitzProducts( Matrix, Step ); Matrix->MaxRowCountInLowerTri = -1; -/* Perform reordering and factorization. */ + /* Perform reordering and factorization. */ for (; Step <= Size; Step++) - { pPivot = SearchForPivot( Matrix, Step, DiagPivoting ); + { + pPivot = SearchForPivot( Matrix, Step, DiagPivoting ); if (pPivot == NULL) return MatrixIsSingular( Matrix, Step ); ExchangeRowsAndCols( Matrix, pPivot, Step ); @@ -284,7 +285,7 @@ RealNumber LargestInCol, FindLargestInCol(); #endif } -Done: + Done: Matrix->NeedsOrdering = NO; Matrix->Reordered = YES; Matrix->Factored = YES; @@ -304,15 +305,15 @@ Done: * This routine is the companion routine to spOrderAndFactor(). * Unlike spOrderAndFactor(), spFactor() cannot change the ordering. * It is also faster than spOrderAndFactor(). The standard way of - * using these two routines is to first use spOrderAndFactor() for the - * initial factorization. For subsequent factorizations, spFactor() - * is used if there is some assurance that little growth will occur - * (say for example, that the matrix is diagonally dominant). If - * spFactor() is called for the initial factorization of the matrix, - * then spOrderAndFactor() is automatically called with the default - * threshold. This routine uses "row at a time" LU factorization. - * Pivots are associated with the lower triangular matrix and the - * diagonals of the upper triangular matrix are ones. + * using these two routines is to first use spOrderAndFactor() for + * the initial factorization. For subsequent factorizations, + * spFactor() is used if there is some assurance that little growth + * will occur (say for example, that the matrix is diagonally + * dominant). If spFactor() is called for the initial factorization + * of the matrix, then spOrderAndFactor() is automatically called + * with the default threshold. This routine uses "row at a time" LU + * factorization. Pivots are associated with the lower triangular + * matrix and the diagonals of the upper triangular matrix are ones. * * >>> Returned: * The error code is returned. Possible errors are listed below. @@ -326,94 +327,98 @@ Done: * spSINGULAR * spZERO_DIAG * spSMALL_PIVOT - * Error is cleared in this function. - */ + * Error is cleared in this function. */ int -spFactor( eMatrix ) - -char *eMatrix; +spFactor(void *eMatrix) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register ElementPtr pElement; -register ElementPtr pColumn; -register int Step, Size; -RealNumber Mult; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + ElementPtr pElement; + ElementPtr pColumn; + int Step, Size; + RealNumber Mult; -/* Begin `spFactor'. */ - ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored); + /* Begin `spFactor'. */ + assert( IS_VALID(Matrix) && !Matrix->Factored); if (Matrix->NeedsOrdering) - { return spOrderAndFactor( eMatrix, (RealVector)NULL, + { + return spOrderAndFactor( eMatrix, (RealVector)NULL, 0.0, 0.0, DIAG_PIVOTING_AS_DEFAULT ); } - if (NOT Matrix->Partitioned) spPartition( eMatrix, spDEFAULT_PARTITION ); -#if spCOMPLEX - if (Matrix->Complex) return FactorComplexMatrix( Matrix ); -#endif + if (!Matrix->Partitioned) spPartition( eMatrix, spDEFAULT_PARTITION ); + if (Matrix->Complex) + return FactorComplexMatrix( Matrix ); -#if REAL Size = Matrix->Size; if (Matrix->Diag[1]->Real == 0.0) return ZeroPivot( Matrix, 1 ); Matrix->Diag[1]->Real = 1.0 / Matrix->Diag[1]->Real; -/* Start factorization. */ + /* Start factorization. */ for (Step = 2; Step <= Size; Step++) - { if (Matrix->DoRealDirect[Step]) - { /* Update column using direct addressing scatter-gather. */ - register RealNumber *Dest = (RealNumber *)Matrix->Intermediate; + { + if (Matrix->DoRealDirect[Step]) + { + /* Update column using direct addressing scatter-gather. */ + RealNumber *Dest = (RealNumber *)Matrix->Intermediate; -/* Scatter. */ + /* Scatter. */ pElement = Matrix->FirstInCol[Step]; while (pElement != NULL) - { Dest[pElement->Row] = pElement->Real; + { + Dest[pElement->Row] = pElement->Real; pElement = pElement->NextInCol; } -/* Update column. */ + /* Update column. */ pColumn = Matrix->FirstInCol[Step]; while (pColumn->Row < Step) - { pElement = Matrix->Diag[pColumn->Row]; + { + pElement = Matrix->Diag[pColumn->Row]; pColumn->Real = Dest[pColumn->Row] * pElement->Real; while ((pElement = pElement->NextInCol) != NULL) Dest[pElement->Row] -= pColumn->Real * pElement->Real; pColumn = pColumn->NextInCol; } -/* Gather. */ + /* Gather. */ pElement = Matrix->Diag[Step]->NextInCol; while (pElement != NULL) - { pElement->Real = Dest[pElement->Row]; + { + pElement->Real = Dest[pElement->Row]; pElement = pElement->NextInCol; } -/* Check for singular matrix. */ + /* Check for singular matrix. */ if (Dest[Step] == 0.0) return ZeroPivot( Matrix, Step ); Matrix->Diag[Step]->Real = 1.0 / Dest[Step]; } else - { /* Update column using indirect addressing scatter-gather. */ - register RealNumber **pDest = (RealNumber **)Matrix->Intermediate; + { + /* Update column using indirect addressing scatter-gather. */ + RealNumber **pDest = (RealNumber **)Matrix->Intermediate; -/* Scatter. */ + /* Scatter. */ pElement = Matrix->FirstInCol[Step]; while (pElement != NULL) - { pDest[pElement->Row] = &pElement->Real; + { + pDest[pElement->Row] = &pElement->Real; pElement = pElement->NextInCol; } -/* Update column. */ + /* Update column. */ pColumn = Matrix->FirstInCol[Step]; while (pColumn->Row < Step) - { pElement = Matrix->Diag[pColumn->Row]; + { + pElement = Matrix->Diag[pColumn->Row]; Mult = (*pDest[pColumn->Row] *= pElement->Real); while ((pElement = pElement->NextInCol) != NULL) *pDest[pElement->Row] -= Mult * pElement->Real; pColumn = pColumn->NextInCol; } -/* Check for singular matrix. */ + /* Check for singular matrix. */ if (Matrix->Diag[Step]->Real == 0.0) return ZeroPivot( Matrix, Step ); Matrix->Diag[Step]->Real = 1.0 / Matrix->Diag[Step]->Real; @@ -422,7 +427,6 @@ RealNumber Mult; Matrix->Factored = YES; return (Matrix->Error = spOKAY); -#endif /* REAL */ } @@ -430,7 +434,6 @@ RealNumber Mult; -#if spCOMPLEX /* * FACTOR COMPLEX MATRIX * @@ -454,87 +457,97 @@ FactorComplexMatrix( Matrix ) MatrixPtr Matrix; { -register ElementPtr pElement; -register ElementPtr pColumn; -register int Step, Size; -ComplexNumber Mult, Pivot; + ElementPtr pElement; + ElementPtr pColumn; + int Step, Size; + ComplexNumber Mult, Pivot; -/* Begin `FactorComplexMatrix'. */ - ASSERT(Matrix->Complex); + /* Begin `FactorComplexMatrix'. */ + assert(Matrix->Complex); Size = Matrix->Size; pElement = Matrix->Diag[1]; if (ELEMENT_MAG(pElement) == 0.0) return ZeroPivot( Matrix, 1 ); -/* Cmplx expr: *pPivot = 1.0 / *pPivot. */ + /* Cmplx expr: *pPivot = 1.0 / *pPivot. */ CMPLX_RECIPROCAL( *pElement, *pElement ); -/* Start factorization. */ + /* Start factorization. */ for (Step = 2; Step <= Size; Step++) - { if (Matrix->DoCmplxDirect[Step]) - { /* Update column using direct addressing scatter-gather. */ - register ComplexNumber *Dest; + { + if (Matrix->DoCmplxDirect[Step]) + { + /* Update column using direct addressing scatter-gather. */ + ComplexNumber *Dest; Dest = (ComplexNumber *)Matrix->Intermediate; -/* Scatter. */ + /* Scatter. */ pElement = Matrix->FirstInCol[Step]; while (pElement != NULL) - { Dest[pElement->Row] = *(ComplexNumber *)pElement; + { + Dest[pElement->Row] = *(ComplexNumber *)pElement; pElement = pElement->NextInCol; } -/* Update column. */ + /* Update column. */ pColumn = Matrix->FirstInCol[Step]; while (pColumn->Row < Step) - { pElement = Matrix->Diag[pColumn->Row]; + { + pElement = Matrix->Diag[pColumn->Row]; /* Cmplx expr: Mult = Dest[pColumn->Row] * (1.0 / *pPivot). */ CMPLX_MULT(Mult, Dest[pColumn->Row], *pElement); CMPLX_ASSIGN(*pColumn, Mult); while ((pElement = pElement->NextInCol) != NULL) - { /* Cmplx expr: Dest[pElement->Row] -= Mult * pElement */ + { + /* Cmplx expr: Dest[pElement->Row] -= Mult * pElement */ CMPLX_MULT_SUBT_ASSIGN(Dest[pElement->Row],Mult,*pElement); } pColumn = pColumn->NextInCol; } -/* Gather. */ + /* Gather. */ pElement = Matrix->Diag[Step]->NextInCol; while (pElement != NULL) - { *(ComplexNumber *)pElement = Dest[pElement->Row]; + { + *(ComplexNumber *)pElement = Dest[pElement->Row]; pElement = pElement->NextInCol; } -/* Check for singular matrix. */ + /* Check for singular matrix. */ Pivot = Dest[Step]; if (CMPLX_1_NORM(Pivot) == 0.0) return ZeroPivot( Matrix, Step ); CMPLX_RECIPROCAL( *Matrix->Diag[Step], Pivot ); } else - { /* Update column using direct addressing scatter-gather. */ - register ComplexNumber **pDest; + { + /* Update column using direct addressing scatter-gather. */ + ComplexNumber **pDest; pDest = (ComplexNumber **)Matrix->Intermediate; -/* Scatter. */ + /* Scatter. */ pElement = Matrix->FirstInCol[Step]; while (pElement != NULL) - { pDest[pElement->Row] = (ComplexNumber *)pElement; + { + pDest[pElement->Row] = (ComplexNumber *)pElement; pElement = pElement->NextInCol; } -/* Update column. */ + /* Update column. */ pColumn = Matrix->FirstInCol[Step]; while (pColumn->Row < Step) - { pElement = Matrix->Diag[pColumn->Row]; + { + pElement = Matrix->Diag[pColumn->Row]; /* Cmplx expr: Mult = *pDest[pColumn->Row] * (1.0 / *pPivot). */ CMPLX_MULT(Mult, *pDest[pColumn->Row], *pElement); CMPLX_ASSIGN(*pDest[pColumn->Row], Mult); while ((pElement = pElement->NextInCol) != NULL) - { /* Cmplx expr: *pDest[pElement->Row] -= Mult * pElement */ - CMPLX_MULT_SUBT_ASSIGN(*pDest[pElement->Row],Mult,*pElement); + { + /* Cmplx expr: *pDest[pElement->Row] -= Mult * pElement */ + CMPLX_MULT_SUBT_ASSIGN(*pDest[pElement->Row],Mult,*pElement); } pColumn = pColumn->NextInCol; } -/* Check for singular matrix. */ + /* Check for singular matrix. */ pElement = Matrix->Diag[Step]; if (ELEMENT_MAG(pElement) == 0.0) return ZeroPivot( Matrix, Step ); CMPLX_RECIPROCAL( *pElement, *pElement ); @@ -544,7 +557,6 @@ ComplexNumber Mult, Pivot; Matrix->Factored = YES; return (Matrix->Error = spOKAY); } -#endif /* spCOMPLEX */ @@ -559,97 +571,95 @@ ComplexNumber Mult, Pivot; * which addressing mode is fastest. This information is used in * spFactor() to speed the factorization. * - * When factoring a previously ordered matrix using spFactor(), Sparse - * operates on a row-at-a-time basis. For speed, on each step, the - * row being updated is copied into a full vector and the operations - * are performed on that vector. This can be done one of two ways, - * either using direct addressing or indirect addressing. Direct - * addressing is fastest when the matrix is relatively dense and - * indirect addressing is best when the matrix is quite sparse. The - * user selects the type of partition used with Mode. If Mode is set - * to spDIRECT_PARTITION, then the all rows are placed in the direct - * addressing partition. Similarly, if Mode is set to + * When factoring a previously ordered matrix using spFactor(), + * Sparse operates on a row-at-a-time basis. For speed, on each + * step, the row being updated is copied into a full vector and the + * operations are performed on that vector. This can be done one of + * two ways, either using direct addressing or indirect addressing. + * Direct addressing is fastest when the matrix is relatively dense + * and indirect addressing is best when the matrix is quite sparse. + * The user selects the type of partition used with Mode. If Mode is + * set to spDIRECT_PARTITION, then the all rows are placed in the + * direct addressing partition. Similarly, if Mode is set to * spINDIRECT_PARTITION, then the all rows are placed in the indirect * addressing partition. By setting Mode to spAUTO_PARTITION, the * user allows Sparse to select the partition for each row * individually. spFactor() generally runs faster if Sparse is * allowed to choose its own partitioning, however choosing a - * partition is expensive. The time required to choose a partition is - * of the same order of the cost to factor the matrix. If you plan to - * factor a large number of matrices with the same structure, it is - * best to let Sparse choose the partition. Otherwise, you should - * choose the partition based on the predicted density of the matrix. + * partition is expensive. The time required to choose a partition + * is of the same order of the cost to factor the matrix. If you + * plan to factor a large number of matrices with the same structure, + * it is best to let Sparse choose the partition. Otherwise, you + * should choose the partition based on the predicted density of the + * matrix. * * >>> Arguments: * Matrix (char *) * Pointer to matrix. * Mode (int) * Mode must be one of three special codes: spDIRECT_PARTITION, - * spINDIRECT_PARTITION, or spAUTO_PARTITION. - */ + * spINDIRECT_PARTITION, or spAUTO_PARTITION. */ void -spPartition( eMatrix, Mode ) - -char *eMatrix; -int Mode; +spPartition(void *eMatrix, int Mode) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register ElementPtr pElement, pColumn; -register int Step, Size; -register int *Nc, *No, *Nm; -BOOLEAN *DoRealDirect, *DoCmplxDirect; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + ElementPtr pElement, pColumn; + int Step, Size; + int *Nc, *No, *Nm; + int *DoRealDirect, *DoCmplxDirect; -/* Begin `spPartition'. */ - ASSERT( IS_SPARSE( Matrix ) ); + /* Begin `spPartition'. */ + assert( IS_SPARSE( Matrix ) ); if (Matrix->Partitioned) return; Size = Matrix->Size; DoRealDirect = Matrix->DoRealDirect; DoCmplxDirect = Matrix->DoCmplxDirect; Matrix->Partitioned = YES; -/* If partition is specified by the user, this is easy. */ + /* If partition is specified by the user, this is easy. */ if (Mode == spDEFAULT_PARTITION) Mode = DEFAULT_PARTITION; if (Mode == spDIRECT_PARTITION) - { for (Step = 1; Step <= Size; Step++) -#if REAL + { + for (Step = 1; Step <= Size; Step++) { DoRealDirect[Step] = YES; -#endif -#if spCOMPLEX - DoCmplxDirect[Step] = YES; -#endif + DoCmplxDirect[Step] = YES; + } return; } else if (Mode == spINDIRECT_PARTITION) - { for (Step = 1; Step <= Size; Step++) -#if REAL + { + for (Step = 1; Step <= Size; Step++) + { DoRealDirect[Step] = NO; -#endif -#if spCOMPLEX - DoCmplxDirect[Step] = NO; -#endif + DoCmplxDirect[Step] = NO; + } return; } - else ASSERT( Mode == spAUTO_PARTITION ); + else + assert( Mode == spAUTO_PARTITION ); -/* Otherwise, count all operations needed in when factoring matrix. */ + /* Otherwise, count all operations needed in when factoring matrix. */ Nc = (int *)Matrix->MarkowitzRow; No = (int *)Matrix->MarkowitzCol; Nm = (int *)Matrix->MarkowitzProd; -/* Start mock-factorization. */ + /* Start mock-factorization. */ for (Step = 1; Step <= Size; Step++) - { Nc[Step] = No[Step] = Nm[Step] = 0; + { + Nc[Step] = No[Step] = Nm[Step] = 0; pElement = Matrix->FirstInCol[Step]; while (pElement != NULL) - { Nc[Step]++; + { + Nc[Step]++; pElement = pElement->NextInCol; } pColumn = Matrix->FirstInCol[Step]; while (pColumn->Row < Step) - { pElement = Matrix->Diag[pColumn->Row]; + { + pElement = Matrix->Diag[pColumn->Row]; Nm[Step]++; while ((pElement = pElement->NextInCol) != NULL) No[Step]++; @@ -659,53 +669,42 @@ BOOLEAN *DoRealDirect, *DoCmplxDirect; for (Step = 1; Step <= Size; Step++) { -/* - * The following are just estimates based on a count on the number of - * machine instructions used on each machine to perform the various - * tasks. It was assumed that each machine instruction required the - * same amount of time (I don't believe this is TRUE for the VAX, and - * have no idea if this is TRUE for the 68000 family). For optimum - * performance, these numbers should be tuned to the machine. - * Nc is the number of nonzero elements in the column. - * Nm is the number of multipliers in the column. - * No is the number of operations in the inner loop. - */ + /* The following are just estimates based on a count on the + * number of machine instructions used on each machine to + * perform the various tasks. It was assumed that each + * machine instruction required the same amount of time (I + * don't believe this is TRUE for the VAX, and have no idea if + * this is TRUE for the 68000 family). For optimum + * performance, these numbers should be tuned to the machine. + * + * Nc is the number of nonzero elements in the column. + * Nm is the number of multipliers in the column. + * No is the number of operations in the inner loop. */ #define generic #ifdef hp9000s300 -#if REAL DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]); -#endif -#if spCOMPLEX /* On the hp350, it is never profitable to use direct for complex. */ DoCmplxDirect[Step] = NO; -#endif #undef generic #endif #ifdef vax -#if REAL DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]); -#endif -#if spCOMPLEX DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7*Nc[Step] - 4*Nm[Step]); -#endif #undef generic #endif #ifdef generic -#if REAL DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]); -#endif -#if spCOMPLEX DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7*Nc[Step] - 4*Nm[Step]); -#endif #undef generic #endif } #if (ANNOTATE == FULL) - { int Ops = 0; + { + int Ops = 0; for (Step = 1; Step <= Size; Step++) Ops += No[Step]; printf("Operation count for inner loop of factorization = %d.\n", Ops); @@ -742,54 +741,46 @@ spcCreateInternalVectors( Matrix ) MatrixPtr Matrix; { -int Size; + int Size; -/* Begin `spcCreateInternalVectors'. */ -/* Create Markowitz arrays. */ + /* Begin `spcCreateInternalVectors'. */ + /* Create Markowitz arrays. */ Size= Matrix->Size; if (Matrix->MarkowitzRow == NULL) - { if (( Matrix->MarkowitzRow = ALLOC(int, Size+1)) == NULL) - Matrix->Error = spNO_MEMORY; + { + if (( Matrix->MarkowitzRow = ALLOC(int, Size+1)) == NULL) + Matrix->Error = spNO_MEMORY; } if (Matrix->MarkowitzCol == NULL) - { if (( Matrix->MarkowitzCol = ALLOC(int, Size+1)) == NULL) - Matrix->Error = spNO_MEMORY; + { + if (( Matrix->MarkowitzCol = ALLOC(int, Size+1)) == NULL) + Matrix->Error = spNO_MEMORY; } if (Matrix->MarkowitzProd == NULL) - { if (( Matrix->MarkowitzProd = ALLOC(long, Size+2)) == NULL) - Matrix->Error = spNO_MEMORY; + { + if (( Matrix->MarkowitzProd = ALLOC(long, Size+2)) == NULL) + Matrix->Error = spNO_MEMORY; } -/* Create DoDirect vectors for use in spFactor(). */ -#if REAL - if (Matrix->DoRealDirect == NULL) - { if (( Matrix->DoRealDirect = ALLOC(BOOLEAN, Size+1)) == NULL) - Matrix->Error = spNO_MEMORY; + /* Create DoDirect vectors for use in spFactor(). */ + if (Matrix->DoRealDirect == NULL) { + if (( Matrix->DoRealDirect = ALLOC(int, Size+1)) == NULL) + Matrix->Error = spNO_MEMORY; } -#endif -#if spCOMPLEX - if (Matrix->DoCmplxDirect == NULL) - { if (( Matrix->DoCmplxDirect = ALLOC(BOOLEAN, Size+1)) == NULL) - Matrix->Error = spNO_MEMORY; + if (Matrix->DoCmplxDirect == NULL) { + if (( Matrix->DoCmplxDirect = ALLOC(int, Size+1)) == NULL) + Matrix->Error = spNO_MEMORY; } -#endif -/* Create Intermediate vectors for use in MatrixSolve. */ -#if spCOMPLEX - if (Matrix->Intermediate == NULL) - { if ((Matrix->Intermediate = ALLOC(RealNumber,2*(Size+1))) == NULL) - Matrix->Error = spNO_MEMORY; + /* Create Intermediate vectors for use in MatrixSolve. */ + if (Matrix->Intermediate == NULL) { + if ((Matrix->Intermediate = ALLOC(RealNumber,2*(Size+1))) == NULL) + Matrix->Error = spNO_MEMORY; } -#else - if (Matrix->Intermediate == NULL) - { if ((Matrix->Intermediate = ALLOC(RealNumber, Size+1)) == NULL) - Matrix->Error = spNO_MEMORY; - } -#endif if (Matrix->Error != spNO_MEMORY) - Matrix->InternalVectorsAllocated = YES; + Matrix->InternalVectorsAllocated = YES; return; } @@ -827,71 +818,44 @@ int Size; */ static void -CountMarkowitz( Matrix, RHS, Step ) - -MatrixPtr Matrix; -register RealVector RHS; -int Step; +CountMarkowitz(MatrixPtr Matrix, RealVector RHS, int Step) { -register int Count, I, Size = Matrix->Size; -register ElementPtr pElement; -int ExtRow; + int Count, I, Size = Matrix->Size; + ElementPtr pElement; + int ExtRow; -/* Begin `CountMarkowitz'. */ + /* Begin `CountMarkowitz'. */ -/* Correct array pointer for ARRAY_OFFSET. */ -#if NOT ARRAY_OFFSET -#if spSEPARATED_COMPLEX_VECTORS OR NOT spCOMPLEX - if (RHS != NULL) --RHS; -#else - if (RHS != NULL) - { if (Matrix->Complex) RHS -= 2; - else --RHS; - } -#endif -#endif - -/* Generate MarkowitzRow Count for each row. */ - for (I = Step; I <= Size; I++) - { -/* Set Count to -1 initially to remove count due to pivot element. */ + /* Generate MarkowitzRow Count for each row. */ + for (I = Step; I <= Size; I++) { + /* Set Count to -1 initially to remove count due to pivot element. */ Count = -1; pElement = Matrix->FirstInRow[I]; - while (pElement != NULL AND pElement->Col < Step) + while (pElement != NULL && pElement->Col < Step) pElement = pElement->NextInRow; - while (pElement != NULL) - { Count++; + while (pElement != NULL) { + Count++; pElement = pElement->NextInRow; } -/* Include nonzero elements in the RHS vector. */ + /* Include nonzero elements in the RHS vector. */ ExtRow = Matrix->IntToExtRowMap[I]; -#if spSEPARATED_COMPLEX_VECTORS OR NOT spCOMPLEX if (RHS != NULL) - if (RHS[ExtRow] != 0.0) Count++; -#else - if (RHS != NULL) - { if (Matrix->Complex) - { if ((RHS[2*ExtRow] != 0.0) OR (RHS[2*ExtRow+1] != 0.0)) - Count++; - } - else if (RHS[I] != 0.0) Count++; - } -#endif + if (RHS[ExtRow] != 0.0) + Count++; Matrix->MarkowitzRow[I] = Count; } -/* Generate the MarkowitzCol count for each column. */ - for (I = Step; I <= Size; I++) - { -/* Set Count to -1 initially to remove count due to pivot element. */ + /* Generate the MarkowitzCol count for each column. */ + for (I = Step; I <= Size; I++) { + /* Set Count to -1 initially to remove count due to pivot element. */ Count = -1; pElement = Matrix->FirstInCol[I]; - while (pElement != NULL AND pElement->Row < Step) + while (pElement != NULL && pElement->Row < Step) pElement = pElement->NextInCol; - while (pElement != NULL) - { Count++; + while (pElement != NULL) { + Count++; pElement = pElement->NextInCol; } Matrix->MarkowitzCol[I] = Count; @@ -937,36 +901,31 @@ int ExtRow; */ static void -MarkowitzProducts( Matrix, Step ) - -MatrixPtr Matrix; -int Step; +MarkowitzProducts(MatrixPtr Matrix, int Step) { -register int I, *pMarkowitzRow, *pMarkowitzCol; -register long Product, *pMarkowitzProduct; -register int Size = Matrix->Size; -double fProduct; + int I, *pMarkowitzRow, *pMarkowitzCol; + long Product, *pMarkowitzProduct; + int Size = Matrix->Size; + double fProduct; -/* Begin `MarkowitzProducts'. */ + /* Begin `MarkowitzProducts'. */ Matrix->Singletons = 0; pMarkowitzProduct = &(Matrix->MarkowitzProd[Step]); pMarkowitzRow = &(Matrix->MarkowitzRow[Step]); pMarkowitzCol = &(Matrix->MarkowitzCol[Step]); - for (I = Step; I <= Size; I++) - { -/* If chance of overflow, use real numbers. */ - if ((*pMarkowitzRow > LARGEST_SHORT_INTEGER AND *pMarkowitzCol != 0) OR - (*pMarkowitzCol > LARGEST_SHORT_INTEGER AND *pMarkowitzRow != 0)) - { fProduct = (double)(*pMarkowitzRow++) * (double)(*pMarkowitzCol++); + for (I = Step; I <= Size; I++) { + /* If chance of overflow, use real numbers. */ + if ((*pMarkowitzRow > LARGEST_SHORT_INTEGER && *pMarkowitzCol != 0) || + (*pMarkowitzCol > LARGEST_SHORT_INTEGER && *pMarkowitzRow != 0)) { + fProduct = (double)(*pMarkowitzRow++) * (double)(*pMarkowitzCol++); if (fProduct >= LARGEST_LONG_INTEGER) *pMarkowitzProduct++ = LARGEST_LONG_INTEGER; else *pMarkowitzProduct++ = fProduct; - } - else - { Product = *pMarkowitzRow++ * *pMarkowitzCol++; + } else { + Product = *pMarkowitzRow++ * *pMarkowitzCol++; if ((*pMarkowitzProduct++ = Product) == 0) Matrix->Singletons++; } @@ -1029,19 +988,21 @@ SearchForPivot( Matrix, Step, DiagPivoting ) MatrixPtr Matrix; int Step, DiagPivoting; { -register ElementPtr ChosenPivot; +ElementPtr ChosenPivot; ElementPtr SearchForSingleton(); ElementPtr QuicklySearchDiagonal(); ElementPtr SearchDiagonal(); ElementPtr SearchEntireMatrix(); -/* Begin `SearchForPivot'. */ + /* Begin `SearchForPivot'. */ -/* If singletons exist, look for an acceptable one to use as pivot. */ + /* If singletons exist, look for an acceptable one to use as pivot. */ if (Matrix->Singletons) - { ChosenPivot = SearchForSingleton( Matrix, Step ); + { + ChosenPivot = SearchForSingleton( Matrix, Step ); if (ChosenPivot != NULL) - { Matrix->PivotSelectionMethod = 's'; + { + Matrix->PivotSelectionMethod = 's'; return ChosenPivot; } } @@ -1049,7 +1010,7 @@ ElementPtr SearchEntireMatrix(); #if DIAGONAL_PIVOTING if (DiagPivoting) { -/* + /* * Either no singletons exist or they weren't acceptable. Take quick first * pass at searching diagonal. First search for element on diagonal of * remaining submatrix with smallest Markowitz product, then check to see @@ -1057,23 +1018,25 @@ ElementPtr SearchEntireMatrix(); */ ChosenPivot = QuicklySearchDiagonal( Matrix, Step ); if (ChosenPivot != NULL) - { Matrix->PivotSelectionMethod = 'q'; + { + Matrix->PivotSelectionMethod = 'q'; return ChosenPivot; } -/* + /* * Quick search of diagonal failed, carefully search diagonal and check each * pivot candidate numerically before even tentatively accepting it. */ ChosenPivot = SearchDiagonal( Matrix, Step ); if (ChosenPivot != NULL) - { Matrix->PivotSelectionMethod = 'd'; + { + Matrix->PivotSelectionMethod = 'd'; return ChosenPivot; } } #endif /* DIAGONAL_PIVOTING */ -/* No acceptable pivot found yet, search entire matrix. */ + /* No acceptable pivot found yet, search entire matrix. */ ChosenPivot = SearchEntireMatrix( Matrix, Step ); Matrix->PivotSelectionMethod = 'e'; @@ -1130,22 +1093,22 @@ SearchForSingleton( Matrix, Step ) MatrixPtr Matrix; int Step; { -register ElementPtr ChosenPivot; -register int I; -register long *pMarkowitzProduct; + ElementPtr ChosenPivot; + int I; + long *pMarkowitzProduct; int Singletons; RealNumber PivotMag, FindBiggestInColExclude(); -/* Begin `SearchForSingleton'. */ -/* Initialize pointer that is to scan through MarkowitzProduct vector. */ + /* Begin `SearchForSingleton'. */ + /* Initialize pointer that is to scan through MarkowitzProduct vector. */ pMarkowitzProduct = &(Matrix->MarkowitzProd[Matrix->Size+1]); Matrix->MarkowitzProd[Matrix->Size+1] = Matrix->MarkowitzProd[Step]; -/* Decrement the count of available singletons, on the assumption that an + /* Decrement the count of available singletons, on the assumption that an * acceptable one will be found. */ Singletons = Matrix->Singletons--; -/* + /* * Assure that following while loop will always terminate, this is just * preventive medicine, if things are working right this should never * be needed. @@ -1154,9 +1117,9 @@ RealNumber PivotMag, FindBiggestInColExclude(); while (Singletons-- > 0) { -/* Singletons exist, find them. */ + /* Singletons exist, find them. */ -/* + /* * This is tricky. Am using a pointer to sequentially step through the * MarkowitzProduct array. Search terminates when singleton (Product = 0) * is found. Note that the conditional in the while statement @@ -1176,7 +1139,8 @@ RealNumber PivotMag, FindBiggestInColExclude(); */ while ( *pMarkowitzProduct-- ) - { /* + { + /* * N bottles of beer on the wall; * N bottles of beer. * you take one down and pass it around; @@ -1185,51 +1149,56 @@ RealNumber PivotMag, FindBiggestInColExclude(); } I = pMarkowitzProduct - Matrix->MarkowitzProd + 1; -/* Assure that I is valid. */ + /* Assure that I is valid. */ if (I < Step) break; /* while (Singletons-- > 0) */ if (I > Matrix->Size) I = Step; -/* Singleton has been found in either/both row or/and column I. */ + /* Singleton has been found in either/both row or/and column I. */ if ((ChosenPivot = Matrix->Diag[I]) != NULL) { -/* Singleton lies on the diagonal. */ + /* Singleton lies on the diagonal. */ PivotMag = ELEMENT_MAG(ChosenPivot); if - ( PivotMag > Matrix->AbsThreshold AND + ( PivotMag > Matrix->AbsThreshold && PivotMag > Matrix->RelThreshold * FindBiggestInColExclude( Matrix, ChosenPivot, Step ) ) return ChosenPivot; } else { -/* Singleton does not lie on diagonal, find it. */ + /* Singleton does not lie on diagonal, find it. */ if (Matrix->MarkowitzCol[I] == 0) - { ChosenPivot = Matrix->FirstInCol[I]; - while ((ChosenPivot != NULL) AND (ChosenPivot->Row < Step)) + { + ChosenPivot = Matrix->FirstInCol[I]; + while ((ChosenPivot != NULL) && (ChosenPivot->Row < Step)) ChosenPivot = ChosenPivot->NextInCol; if (ChosenPivot != NULL) - { /* Reduced column has no elements, matrix is singular. */ + { + /* Reduced column has no elements, matrix is singular. */ break; } PivotMag = ELEMENT_MAG(ChosenPivot); if - ( PivotMag > Matrix->AbsThreshold AND + ( PivotMag > Matrix->AbsThreshold && PivotMag > Matrix->RelThreshold * FindBiggestInColExclude( Matrix, ChosenPivot, Step ) ) return ChosenPivot; else - { if (Matrix->MarkowitzRow[I] == 0) - { ChosenPivot = Matrix->FirstInRow[I]; - while((ChosenPivot != NULL) AND (ChosenPivot->ColMarkowitzRow[I] == 0) + { + ChosenPivot = Matrix->FirstInRow[I]; + while((ChosenPivot != NULL) && (ChosenPivot->ColNextInRow; if (ChosenPivot != NULL) - { /* Reduced row has no elements, matrix is singular. */ + { + /* Reduced row has no elements, matrix is singular. */ break; } PivotMag = ELEMENT_MAG(ChosenPivot); if - ( PivotMag > Matrix->AbsThreshold AND + ( PivotMag > Matrix->AbsThreshold && PivotMag > Matrix->RelThreshold * FindBiggestInColExclude( Matrix, ChosenPivot, @@ -1239,26 +1208,28 @@ RealNumber PivotMag, FindBiggestInColExclude(); } } else - { ChosenPivot = Matrix->FirstInRow[I]; - while ((ChosenPivot != NULL) AND (ChosenPivot->Col < Step)) + { + ChosenPivot = Matrix->FirstInRow[I]; + while ((ChosenPivot != NULL) && (ChosenPivot->Col < Step)) ChosenPivot = ChosenPivot->NextInRow; if (ChosenPivot != NULL) - { /* Reduced row has no elements, matrix is singular. */ + { + /* Reduced row has no elements, matrix is singular. */ break; } PivotMag = ELEMENT_MAG(ChosenPivot); if - ( PivotMag > Matrix->AbsThreshold AND + ( PivotMag > Matrix->AbsThreshold && PivotMag > Matrix->RelThreshold * FindBiggestInColExclude( Matrix, ChosenPivot, Step ) ) return ChosenPivot; } } -/* Singleton not acceptable (too small), try another. */ + /* Singleton not acceptable (too small), try another. */ } /* end of while(lSingletons>0) */ -/* + /* * All singletons were unacceptable. Restore Matrix->Singletons count. * Initial assumption that an acceptable singleton would be found was wrong. */ @@ -1352,24 +1323,24 @@ QuicklySearchDiagonal( Matrix, Step ) MatrixPtr Matrix; int Step; { -register long MinMarkowitzProduct, *pMarkowitzProduct; -register ElementPtr pDiag, pOtherInRow, pOtherInCol; +long MinMarkowitzProduct, *pMarkowitzProduct; + ElementPtr pDiag, pOtherInRow, pOtherInCol; int I, NumberOfTies; ElementPtr ChosenPivot, TiedElements[MAX_MARKOWITZ_TIES + 1]; RealNumber Magnitude, LargestInCol, Ratio, MaxRatio; RealNumber LargestOffDiagonal; RealNumber FindBiggestInColExclude(); -/* Begin `QuicklySearchDiagonal'. */ + /* Begin `QuicklySearchDiagonal'. */ NumberOfTies = -1; MinMarkowitzProduct = LARGEST_LONG_INTEGER; pMarkowitzProduct = &(Matrix->MarkowitzProd[Matrix->Size+2]); Matrix->MarkowitzProd[Matrix->Size+1] = Matrix->MarkowitzProd[Step]; -/* Assure that following while loop will always terminate. */ + /* Assure that following while loop will always terminate. */ Matrix->MarkowitzProd[Step-1] = -1; -/* + /* * This is tricky. Am using a pointer in the inner while loop to * sequentially step through the MarkowitzProduct array. Search * terminates when the Markowitz product of zero placed at location @@ -1389,8 +1360,10 @@ RealNumber FindBiggestInColExclude(); */ for(;;) /* Endless for loop. */ - { while (MinMarkowitzProduct < *(--pMarkowitzProduct)) - { /* + { + while (MinMarkowitzProduct < *(--pMarkowitzProduct)) + { + /* * N bottles of beer on the wall; * N bottles of beer. * You take one down and pass it around; @@ -1400,7 +1373,7 @@ RealNumber FindBiggestInColExclude(); I = pMarkowitzProduct - Matrix->MarkowitzProd; -/* Assure that I is valid; if I < Step, terminate search. */ + /* Assure that I is valid; if I < Step, terminate search. */ if (I < Step) break; /* Endless for loop */ if (I > Matrix->Size) I = Step; @@ -1411,35 +1384,40 @@ RealNumber FindBiggestInColExclude(); if (*pMarkowitzProduct == 1) { -/* Case where only one element exists in row and column other than diagonal. */ + /* Case where only one element exists in row and column other than diagonal. */ -/* Find off diagonal elements. */ + /* Find off diagonal elements. */ pOtherInRow = pDiag->NextInRow; pOtherInCol = pDiag->NextInCol; - if (pOtherInRow == NULL AND pOtherInCol == NULL) - { pOtherInRow = Matrix->FirstInRow[I]; + if (pOtherInRow == NULL && pOtherInCol == NULL) + { + pOtherInRow = Matrix->FirstInRow[I]; while(pOtherInRow != NULL) - { if (pOtherInRow->Col >= Step AND pOtherInRow->Col != I) + { + if (pOtherInRow->Col >= Step && pOtherInRow->Col != I) break; pOtherInRow = pOtherInRow->NextInRow; } pOtherInCol = Matrix->FirstInCol[I]; while(pOtherInCol != NULL) - { if (pOtherInCol->Row >= Step AND pOtherInCol->Row != I) + { + if (pOtherInCol->Row >= Step && pOtherInCol->Row != I) break; pOtherInCol = pOtherInCol->NextInCol; } } -/* Accept diagonal as pivot if diagonal is larger than off diagonals and the + /* Accept diagonal as pivot if diagonal is larger than off diagonals and the * off diagonals are placed symmetricly. */ - if (pOtherInRow != NULL AND pOtherInCol != NULL) - { if (pOtherInRow->Col == pOtherInCol->Row) - { LargestOffDiagonal = MAX(ELEMENT_MAG(pOtherInRow), + if (pOtherInRow != NULL && pOtherInCol != NULL) + { + if (pOtherInRow->Col == pOtherInCol->Row) + { + LargestOffDiagonal = MAX(ELEMENT_MAG(pOtherInRow), ELEMENT_MAG(pOtherInCol)); if (Magnitude >= LargestOffDiagonal) { -/* Accept pivot, it is unlikely to contribute excess error. */ + /* Accept pivot, it is unlikely to contribute excess error. */ return pDiag; } } @@ -1448,37 +1426,40 @@ RealNumber FindBiggestInColExclude(); if (*pMarkowitzProduct < MinMarkowitzProduct) { -/* Notice strict inequality in test. This is a new smallest MarkowitzProduct. */ + /* Notice strict inequality in test. This is a new smallest MarkowitzProduct. */ TiedElements[0] = pDiag; MinMarkowitzProduct = *pMarkowitzProduct; NumberOfTies = 0; } else { -/* This case handles Markowitz ties. */ + /* This case handles Markowitz ties. */ if (NumberOfTies < MAX_MARKOWITZ_TIES) - { TiedElements[++NumberOfTies] = pDiag; + { + TiedElements[++NumberOfTies] = pDiag; if (NumberOfTies >= MinMarkowitzProduct * TIES_MULTIPLIER) break; /* Endless for loop */ } } } /* End of endless for loop. */ -/* Test to see if any element was chosen as a pivot candidate. */ + /* Test to see if any element was chosen as a pivot candidate. */ if (NumberOfTies < 0) return NULL; -/* Determine which of tied elements is best numerically. */ + /* Determine which of tied elements is best numerically. */ ChosenPivot = NULL; MaxRatio = 1.0 / Matrix->RelThreshold; for (I = 0; I <= NumberOfTies; I++) - { pDiag = TiedElements[I]; + { + pDiag = TiedElements[I]; Magnitude = ELEMENT_MAG(pDiag); LargestInCol = FindBiggestInColExclude( Matrix, pDiag, Step ); Ratio = LargestInCol / Magnitude; if (Ratio < MaxRatio) - { ChosenPivot = pDiag; + { + ChosenPivot = pDiag; MaxRatio = Ratio; } } @@ -1551,23 +1532,23 @@ QuicklySearchDiagonal( Matrix, Step ) MatrixPtr Matrix; int Step; { -register long MinMarkowitzProduct, *pMarkowitzProduct; -register ElementPtr pDiag; +long MinMarkowitzProduct, *pMarkowitzProduct; + ElementPtr pDiag; int I; ElementPtr ChosenPivot, pOtherInRow, pOtherInCol; RealNumber Magnitude, LargestInCol, LargestOffDiagonal; RealNumber FindBiggestInColExclude(); -/* Begin `QuicklySearchDiagonal'. */ + /* Begin `QuicklySearchDiagonal'. */ ChosenPivot = NULL; MinMarkowitzProduct = LARGEST_LONG_INTEGER; pMarkowitzProduct = &(Matrix->MarkowitzProd[Matrix->Size+2]); Matrix->MarkowitzProd[Matrix->Size+1] = Matrix->MarkowitzProd[Step]; -/* Assure that following while loop will always terminate. */ + /* Assure that following while loop will always terminate. */ Matrix->MarkowitzProd[Step-1] = -1; -/* + /* * This is tricky. Am using a pointer in the inner while loop to * sequentially step through the MarkowitzProduct array. Search * terminates when the Markowitz product of zero placed at location @@ -1587,13 +1568,15 @@ RealNumber FindBiggestInColExclude(); */ for (;;) /* Endless for loop. */ - { while (*(--pMarkowitzProduct) >= MinMarkowitzProduct) - { /* Just passing through. */ + { + while (*(--pMarkowitzProduct) >= MinMarkowitzProduct) + { + /* Just passing through. */ } I = pMarkowitzProduct - Matrix->MarkowitzProd; -/* Assure that I is valid; if I < Step, terminate search. */ + /* Assure that I is valid; if I < Step, terminate search. */ if (I < Step) break; /* Endless for loop */ if (I > Matrix->Size) I = Step; @@ -1604,35 +1587,40 @@ RealNumber FindBiggestInColExclude(); if (*pMarkowitzProduct == 1) { -/* Case where only one element exists in row and column other than diagonal. */ + /* Case where only one element exists in row and column other than diagonal. */ -/* Find off-diagonal elements. */ + /* Find off-diagonal elements. */ pOtherInRow = pDiag->NextInRow; pOtherInCol = pDiag->NextInCol; - if (pOtherInRow == NULL AND pOtherInCol == NULL) - { pOtherInRow = Matrix->FirstInRow[I]; + if (pOtherInRow == NULL && pOtherInCol == NULL) + { + pOtherInRow = Matrix->FirstInRow[I]; while(pOtherInRow != NULL) - { if (pOtherInRow->Col >= Step AND pOtherInRow->Col != I) + { + if (pOtherInRow->Col >= Step && pOtherInRow->Col != I) break; pOtherInRow = pOtherInRow->NextInRow; } pOtherInCol = Matrix->FirstInCol[I]; while(pOtherInCol != NULL) - { if (pOtherInCol->Row >= Step AND pOtherInCol->Row != I) + { + if (pOtherInCol->Row >= Step && pOtherInCol->Row != I) break; pOtherInCol = pOtherInCol->NextInCol; } } -/* Accept diagonal as pivot if diagonal is larger than off-diagonals and the + /* Accept diagonal as pivot if diagonal is larger than off-diagonals and the * off-diagonals are placed symmetricly. */ - if (pOtherInRow != NULL AND pOtherInCol != NULL) - { if (pOtherInRow->Col == pOtherInCol->Row) - { LargestOffDiagonal = MAX(ELEMENT_MAG(pOtherInRow), + if (pOtherInRow != NULL && pOtherInCol != NULL) + { + if (pOtherInRow->Col == pOtherInCol->Row) + { + LargestOffDiagonal = MAX(ELEMENT_MAG(pOtherInRow), ELEMENT_MAG(pOtherInCol)); if (Magnitude >= LargestOffDiagonal) { -/* Accept pivot, it is unlikely to contribute excess error. */ + /* Accept pivot, it is unlikely to contribute excess error. */ return pDiag; } } @@ -1644,7 +1632,8 @@ RealNumber FindBiggestInColExclude(); } /* End of endless for loop. */ if (ChosenPivot != NULL) - { LargestInCol = FindBiggestInColExclude( Matrix, ChosenPivot, Step ); + { + LargestInCol = FindBiggestInColExclude( Matrix, ChosenPivot, Step ); if( ELEMENT_MAG(ChosenPivot) <= Matrix->RelThreshold * LargestInCol ) ChosenPivot = NULL; } @@ -1688,7 +1677,7 @@ RealNumber FindBiggestInColExclude(); * ChosenPivot (ElementPtr) * Pointer to the element that has been chosen to be the pivot. * Size (int) - * Local version of size which is placed in a register to increase speed. + * Local version of size which is placed in a to increase speed. * Magnitude (RealNumber) * Absolute value of diagonal element. * MinMarkowitzProduct (long) @@ -1714,24 +1703,24 @@ static ElementPtr SearchDiagonal( Matrix, Step ) MatrixPtr Matrix; -register int Step; +int Step; { -register int J; -register long MinMarkowitzProduct, *pMarkowitzProduct; -register int I; -register ElementPtr pDiag; -int NumberOfTies, Size = Matrix->Size; -ElementPtr ChosenPivot; -RealNumber Magnitude, Ratio, RatioOfAccepted, LargestInCol; -RealNumber FindBiggestInColExclude(); + int J; + long MinMarkowitzProduct, *pMarkowitzProduct; + int I; + ElementPtr pDiag; + int NumberOfTies, Size = Matrix->Size; + ElementPtr ChosenPivot; + RealNumber Magnitude, Ratio, RatioOfAccepted, LargestInCol; + RealNumber FindBiggestInColExclude(); -/* Begin `SearchDiagonal'. */ + /* Begin `SearchDiagonal'. */ ChosenPivot = NULL; MinMarkowitzProduct = LARGEST_LONG_INTEGER; pMarkowitzProduct = &(Matrix->MarkowitzProd[Size+2]); Matrix->MarkowitzProd[Size+1] = Matrix->MarkowitzProd[Step]; -/* Start search of diagonal. */ + /* Start search of diagonal. */ for (J = Size+1; J > Step; J--) { if (*(--pMarkowitzProduct) > MinMarkowitzProduct) @@ -1745,14 +1734,15 @@ RealNumber FindBiggestInColExclude(); if ((Magnitude = ELEMENT_MAG(pDiag)) <= Matrix->AbsThreshold) continue; /* for loop */ -/* Test to see if diagonal's magnitude is acceptable. */ + /* Test to see if diagonal's magnitude is acceptable. */ LargestInCol = FindBiggestInColExclude( Matrix, pDiag, Step ); if (Magnitude <= Matrix->RelThreshold * LargestInCol) continue; /* for loop */ if (*pMarkowitzProduct < MinMarkowitzProduct) { -/* Notice strict inequality in test. This is a new smallest MarkowitzProduct. */ + /* Notice strict inequality in test. This is a new + smallest MarkowitzProduct. */ ChosenPivot = pDiag; MinMarkowitzProduct = *pMarkowitzProduct; RatioOfAccepted = LargestInCol / Magnitude; @@ -1760,11 +1750,12 @@ RealNumber FindBiggestInColExclude(); } else { -/* This case handles Markowitz ties. */ + /* This case handles Markowitz ties. */ NumberOfTies++; Ratio = LargestInCol / Magnitude; if (Ratio < RatioOfAccepted) - { ChosenPivot = pDiag; + { + ChosenPivot = pDiag; RatioOfAccepted = Ratio; } if (NumberOfTies >= MinMarkowitzProduct * TIES_MULTIPLIER) @@ -1812,7 +1803,7 @@ RealNumber FindBiggestInColExclude(); * LargestElementMag (RealNumber) * Magnitude of the largest element yet found in the reduced submatrix. * Size (int) - * Local version of Size; placed in a register for speed. + * Local version of Size; placed in a for speed. * Magnitude (RealNumber) * Absolute value of diagonal element. * MinMarkowitzProduct (long) @@ -1845,24 +1836,25 @@ SearchEntireMatrix( Matrix, Step ) MatrixPtr Matrix; int Step; { -register int I, Size = Matrix->Size; -register ElementPtr pElement; -int NumberOfTies; -long Product, MinMarkowitzProduct; -ElementPtr ChosenPivot, pLargestElement; -RealNumber Magnitude, LargestElementMag, Ratio, RatioOfAccepted, LargestInCol; -RealNumber FindLargestInCol(); + int I, Size = Matrix->Size; + ElementPtr pElement; + int NumberOfTies; + long Product, MinMarkowitzProduct; + ElementPtr ChosenPivot, pLargestElement; + RealNumber Magnitude, LargestElementMag, Ratio, RatioOfAccepted, LargestInCol; + RealNumber FindLargestInCol(); -/* Begin `SearchEntireMatrix'. */ + /* Begin `SearchEntireMatrix'. */ ChosenPivot = NULL; LargestElementMag = 0.0; MinMarkowitzProduct = LARGEST_LONG_INTEGER; -/* Start search of matrix on column by column basis. */ + /* Start search of matrix on column by column basis. */ for (I = Step; I <= Size; I++) - { pElement = Matrix->FirstInCol[I]; + { + pElement = Matrix->FirstInCol[I]; - while (pElement != NULL AND pElement->Row < Step) + while (pElement != NULL && pElement->Row < Step) pElement = pElement->NextInCol; if((LargestInCol = FindLargestInCol(pElement)) == 0.0) @@ -1870,26 +1862,30 @@ RealNumber FindLargestInCol(); while (pElement != NULL) { -/* Check to see if element is the largest encountered so far. If so, record - its magnitude and address. */ + /* Check to see if element is the largest encountered so + far. If so, record its magnitude and address. */ if ((Magnitude = ELEMENT_MAG(pElement)) > LargestElementMag) - { LargestElementMag = Magnitude; + { + LargestElementMag = Magnitude; pLargestElement = pElement; } -/* Calculate element's MarkowitzProduct. */ + /* Calculate element's MarkowitzProduct. */ Product = Matrix->MarkowitzRow[pElement->Row] * - Matrix->MarkowitzCol[pElement->Col]; + Matrix->MarkowitzCol[pElement->Col]; -/* Test to see if element is acceptable as a pivot candidate. */ - if ((Product <= MinMarkowitzProduct) AND - (Magnitude > Matrix->RelThreshold * LargestInCol) AND + /* Test to see if element is acceptable as a pivot + candidate. */ + if ((Product <= MinMarkowitzProduct) && + (Magnitude > Matrix->RelThreshold * LargestInCol) && (Magnitude > Matrix->AbsThreshold)) { -/* Test to see if element has lowest MarkowitzProduct yet found, or whether it - is tied with an element found earlier. */ + /* Test to see if element has lowest MarkowitzProduct + yet found, or whether it is tied with an element + found earlier. */ if (Product < MinMarkowitzProduct) { -/* Notice strict inequality in test. This is a new smallest MarkowitzProduct. */ + /* Notice strict inequality in test. This is a new + smallest MarkowitzProduct. */ ChosenPivot = pElement; MinMarkowitzProduct = Product; RatioOfAccepted = LargestInCol / Magnitude; @@ -1897,11 +1893,12 @@ RealNumber FindLargestInCol(); } else { -/* This case handles Markowitz ties. */ + /* This case handles Markowitz ties. */ NumberOfTies++; Ratio = LargestInCol / Magnitude; if (Ratio < RatioOfAccepted) - { ChosenPivot = pElement; + { + ChosenPivot = pElement; RatioOfAccepted = Ratio; } if (NumberOfTies >= MinMarkowitzProduct * TIES_MULTIPLIER) @@ -1915,7 +1912,8 @@ RealNumber FindLargestInCol(); if (ChosenPivot != NULL) return ChosenPivot; if (LargestElementMag == 0.0) - { Matrix->Error = spSINGULAR; + { + Matrix->Error = spSINGULAR; return NULL; } @@ -1937,14 +1935,14 @@ RealNumber FindLargestInCol(); /* * DETERMINE THE MAGNITUDE OF THE LARGEST ELEMENT IN A COLUMN * - * This routine searches a column and returns the magnitude of the largest - * element. This routine begins the search at the element pointed to by - * pElement, the parameter. + * This routine searches a column and returns the magnitude of the + * largest element. This routine begins the search at the element + * pointed to by pElement, the parameter. * - * The search is conducted by starting at the element specified by a pointer, - * which should be one below the diagonal, and moving down the column. On - * the way down the column, the magnitudes of the elements are tested to see - * if they are the largest yet found. + * The search is conducted by starting at the element specified by a + * pointer, which should be one below the diagonal, and moving down + * the column. On the way down the column, the magnitudes of the + * elements are tested to see if they are the largest yet found. * * >>> Returned: * The magnitude of the largest element in the column below and including @@ -1959,20 +1957,20 @@ RealNumber FindLargestInCol(); * Largest (RealNumber) * The magnitude of the largest element. * Magnitude (RealNumber) - * The magnitude of the currently active element. - */ + * The magnitude of the currently active element. */ static RealNumber FindLargestInCol( pElement ) -register ElementPtr pElement; + ElementPtr pElement; { -RealNumber Magnitude, Largest = 0.0; + RealNumber Magnitude, Largest = 0.0; -/* Begin `FindLargestInCol'. */ -/* Search column for largest element beginning at Element. */ + /* Begin `FindLargestInCol'. */ + /* Search column for largest element beginning at Element. */ while (pElement != NULL) - { if ((Magnitude = ELEMENT_MAG(pElement)) > Largest) + { + if ((Magnitude = ELEMENT_MAG(pElement)) > Largest) Largest = Magnitude; pElement = pElement->NextInCol; } @@ -2032,32 +2030,34 @@ static RealNumber FindBiggestInColExclude( Matrix, pElement, Step ) MatrixPtr Matrix; -register ElementPtr pElement; -register int Step; + ElementPtr pElement; + int Step; { -register int Row; + int Row; int Col; RealNumber Largest, Magnitude; -/* Begin `FindBiggestInColExclude'. */ + /* Begin `FindBiggestInColExclude'. */ Row = pElement->Row; Col = pElement->Col; pElement = Matrix->FirstInCol[Col]; -/* Travel down column until reduced submatrix is entered. */ - while ((pElement != NULL) AND (pElement->Row < Step)) + /* Travel down column until reduced submatrix is entered. */ + while ((pElement != NULL) && (pElement->Row < Step)) pElement = pElement->NextInCol; -/* Initialize the variable Largest. */ + /* Initialize the variable Largest. */ if (pElement->Row != Row) Largest = ELEMENT_MAG(pElement); else Largest = 0.0; -/* Search rest of column for largest element, avoiding excluded element. */ + /* Search rest of column for largest element, avoiding excluded element. */ while ((pElement = pElement->NextInCol) != NULL) - { if ((Magnitude = ELEMENT_MAG(pElement)) > Largest) - { if (pElement->Row != Row) + { + if ((Magnitude = ELEMENT_MAG(pElement)) > Largest) + { + if (pElement->Row != Row) Largest = Magnitude; } } @@ -2110,23 +2110,24 @@ ExchangeRowsAndCols( Matrix, pPivot, Step ) MatrixPtr Matrix; ElementPtr pPivot; -register int Step; +int Step; { -register int Row, Col; + int Row, Col; long OldMarkowitzProd_Step, OldMarkowitzProd_Row, OldMarkowitzProd_Col; ElementPtr spcFindElementInCol(); -/* Begin `ExchangeRowsAndCols'. */ + /* Begin `ExchangeRowsAndCols'. */ Row = pPivot->Row; Col = pPivot->Col; Matrix->PivotsOriginalRow = Row; Matrix->PivotsOriginalCol = Col; - if ((Row == Step) AND (Col == Step)) return; + if ((Row == Step) && (Col == Step)) return; -/* Exchange rows and columns. */ + /* Exchange rows and columns. */ if (Row == Col) - { spcRowExchange( Matrix, Step, Row ); + { + spcRowExchange( Matrix, Step, Row ); spcColExchange( Matrix, Step, Col ); SWAP( long, Matrix->MarkowitzProd[Step], Matrix->MarkowitzProd[Row] ); SWAP( ElementPtr, Matrix->Diag[Row], Matrix->Diag[Step] ); @@ -2134,39 +2135,43 @@ ElementPtr spcFindElementInCol(); else { -/* Initialize variables that hold old Markowitz products. */ + /* Initialize variables that hold old Markowitz products. */ OldMarkowitzProd_Step = Matrix->MarkowitzProd[Step]; OldMarkowitzProd_Row = Matrix->MarkowitzProd[Row]; OldMarkowitzProd_Col = Matrix->MarkowitzProd[Col]; -/* Exchange rows. */ + /* Exchange rows. */ if (Row != Step) - { spcRowExchange( Matrix, Step, Row ); + { + spcRowExchange( Matrix, Step, Row ); Matrix->NumberOfInterchangesIsOdd = - NOT Matrix->NumberOfInterchangesIsOdd; + !Matrix->NumberOfInterchangesIsOdd; Matrix->MarkowitzProd[Row] = Matrix->MarkowitzRow[Row] * Matrix->MarkowitzCol[Row]; -/* Update singleton count. */ + /* Update singleton count. */ if ((Matrix->MarkowitzProd[Row]==0) != (OldMarkowitzProd_Row==0)) - { if (OldMarkowitzProd_Row == 0) + { + if (OldMarkowitzProd_Row == 0) Matrix->Singletons--; else Matrix->Singletons++; } } -/* Exchange columns. */ + /* Exchange columns. */ if (Col != Step) - { spcColExchange( Matrix, Step, Col ); + { + spcColExchange( Matrix, Step, Col ); Matrix->NumberOfInterchangesIsOdd = - NOT Matrix->NumberOfInterchangesIsOdd; + !Matrix->NumberOfInterchangesIsOdd; Matrix->MarkowitzProd[Col] = Matrix->MarkowitzCol[Col] * Matrix->MarkowitzRow[Col]; -/* Update singleton count. */ + /* Update singleton count. */ if ((Matrix->MarkowitzProd[Col]==0) != (OldMarkowitzProd_Col==0)) - { if (OldMarkowitzProd_Col == 0) + { + if (OldMarkowitzProd_Col == 0) Matrix->Singletons--; else Matrix->Singletons++; @@ -2177,7 +2182,8 @@ ElementPtr spcFindElementInCol(); Col, Col, NO ); } if (Row != Step) - { Matrix->Diag[Row] = spcFindElementInCol( Matrix, + { + Matrix->Diag[Row] = spcFindElementInCol( Matrix, Matrix->FirstInCol+Row, Row, Row, NO ); } @@ -2185,11 +2191,12 @@ ElementPtr spcFindElementInCol(); Matrix->FirstInCol+Step, Step, Step, NO ); -/* Update singleton count. */ + /* Update singleton count. */ Matrix->MarkowitzProd[Step] = Matrix->MarkowitzCol[Step] * Matrix->MarkowitzRow[Step]; if ((Matrix->MarkowitzProd[Step]==0) != (OldMarkowitzProd_Step==0)) - { if (OldMarkowitzProd_Step == 0) + { + if (OldMarkowitzProd_Step == 0) Matrix->Singletons--; else Matrix->Singletons++; @@ -2241,44 +2248,49 @@ spcRowExchange( Matrix, Row1, Row2 ) MatrixPtr Matrix; int Row1, Row2; { -register ElementPtr Row1Ptr, Row2Ptr; + ElementPtr Row1Ptr, Row2Ptr; int Column; ElementPtr Element1, Element2; -/* Begin `spcRowExchange'. */ + /* Begin `spcRowExchange'. */ if (Row1 > Row2) SWAP(int, Row1, Row2); Row1Ptr = Matrix->FirstInRow[Row1]; Row2Ptr = Matrix->FirstInRow[Row2]; - while (Row1Ptr != NULL OR Row2Ptr != NULL) + while (Row1Ptr != NULL || Row2Ptr != NULL) { -/* Exchange elements in rows while traveling from left to right. */ + /* Exchange elements in rows while traveling from left to right. */ if (Row1Ptr == NULL) - { Column = Row2Ptr->Col; + { + Column = Row2Ptr->Col; Element1 = NULL; Element2 = Row2Ptr; Row2Ptr = Row2Ptr->NextInRow; } else if (Row2Ptr == NULL) - { Column = Row1Ptr->Col; + { + Column = Row1Ptr->Col; Element1 = Row1Ptr; Element2 = NULL; Row1Ptr = Row1Ptr->NextInRow; } else if (Row1Ptr->Col < Row2Ptr->Col) - { Column = Row1Ptr->Col; + { + Column = Row1Ptr->Col; Element1 = Row1Ptr; Element2 = NULL; Row1Ptr = Row1Ptr->NextInRow; } else if (Row1Ptr->Col > Row2Ptr->Col) - { Column = Row2Ptr->Col; + { + Column = Row2Ptr->Col; Element1 = NULL; Element2 = Row2Ptr; Row2Ptr = Row2Ptr->NextInRow; } else /* Row1Ptr->Col == Row2Ptr->Col */ - { Column = Row1Ptr->Col; + { + Column = Row1Ptr->Col; Element1 = Row1Ptr; Element2 = Row2Ptr; Row1Ptr = Row1Ptr->NextInRow; @@ -2286,7 +2298,7 @@ ElementPtr Element1, Element2; } ExchangeColElements( Matrix, Row1, Element1, Row2, Element2, Column); - } /* end of while(Row1Ptr != NULL OR Row2Ptr != NULL) */ + } /* end of while(Row1Ptr != NULL || Row2Ptr != NULL) */ if (Matrix->InternalVectorsAllocated) SWAP( int, Matrix->MarkowitzRow[Row1], Matrix->MarkowitzRow[Row2]); @@ -2343,44 +2355,49 @@ spcColExchange( Matrix, Col1, Col2 ) MatrixPtr Matrix; int Col1, Col2; { -register ElementPtr Col1Ptr, Col2Ptr; + ElementPtr Col1Ptr, Col2Ptr; int Row; ElementPtr Element1, Element2; -/* Begin `spcColExchange'. */ + /* Begin `spcColExchange'. */ if (Col1 > Col2) SWAP(int, Col1, Col2); Col1Ptr = Matrix->FirstInCol[Col1]; Col2Ptr = Matrix->FirstInCol[Col2]; - while (Col1Ptr != NULL OR Col2Ptr != NULL) + while (Col1Ptr != NULL || Col2Ptr != NULL) { -/* Exchange elements in rows while traveling from top to bottom. */ + /* Exchange elements in rows while traveling from top to bottom. */ if (Col1Ptr == NULL) - { Row = Col2Ptr->Row; + { + Row = Col2Ptr->Row; Element1 = NULL; Element2 = Col2Ptr; Col2Ptr = Col2Ptr->NextInCol; } else if (Col2Ptr == NULL) - { Row = Col1Ptr->Row; + { + Row = Col1Ptr->Row; Element1 = Col1Ptr; Element2 = NULL; Col1Ptr = Col1Ptr->NextInCol; } else if (Col1Ptr->Row < Col2Ptr->Row) - { Row = Col1Ptr->Row; + { + Row = Col1Ptr->Row; Element1 = Col1Ptr; Element2 = NULL; Col1Ptr = Col1Ptr->NextInCol; } else if (Col1Ptr->Row > Col2Ptr->Row) - { Row = Col2Ptr->Row; + { + Row = Col2Ptr->Row; Element1 = NULL; Element2 = Col2Ptr; Col2Ptr = Col2Ptr->NextInCol; } else /* Col1Ptr->Row == Col2Ptr->Row */ - { Row = Col1Ptr->Row; + { + Row = Col1Ptr->Row; Element1 = Col1Ptr; Element2 = Col2Ptr; Col1Ptr = Col1Ptr->NextInCol; @@ -2388,7 +2405,7 @@ ElementPtr Element1, Element2; } ExchangeRowElements( Matrix, Col1, Element1, Col2, Element2, Row); - } /* end of while(Col1Ptr != NULL OR Col2Ptr != NULL) */ + } /* end of while(Col1Ptr != NULL || Col2Ptr != NULL) */ if (Matrix->InternalVectorsAllocated) SWAP( int, Matrix->MarkowitzCol[Col1], Matrix->MarkowitzCol[Col2]); @@ -2448,39 +2465,42 @@ static void ExchangeColElements( Matrix, Row1, Element1, Row2, Element2, Column ) MatrixPtr Matrix; -register ElementPtr Element1, Element2; + ElementPtr Element1, Element2; int Row1, Row2, Column; { ElementPtr *ElementAboveRow1, *ElementAboveRow2; ElementPtr ElementBelowRow1, ElementBelowRow2; -register ElementPtr pElement; + ElementPtr pElement; -/* Begin `ExchangeColElements'. */ -/* Search to find the ElementAboveRow1. */ + /* Begin `ExchangeColElements'. */ + /* Search to find the ElementAboveRow1. */ ElementAboveRow1 = &(Matrix->FirstInCol[Column]); pElement = *ElementAboveRow1; while (pElement->Row < Row1) - { ElementAboveRow1 = &(pElement->NextInCol); + { + ElementAboveRow1 = &(pElement->NextInCol); pElement = *ElementAboveRow1; } if (Element1 != NULL) - { ElementBelowRow1 = Element1->NextInCol; + { + ElementBelowRow1 = Element1->NextInCol; if (Element2 == NULL) { -/* Element2 does not exist, move Element1 down to Row2. */ - if ( ElementBelowRow1 != NULL AND ElementBelowRow1->Row < Row2 ) + /* Element2 does not exist, move Element1 down to Row2. */ + if ( ElementBelowRow1 != NULL && ElementBelowRow1->Row < Row2 ) { -/* Element1 must be removed from linked list and moved. */ + /* Element1 must be removed from linked list and moved. */ *ElementAboveRow1 = ElementBelowRow1; -/* Search column for Row2. */ + /* Search column for Row2. */ pElement = ElementBelowRow1; do - { ElementAboveRow2 = &(pElement->NextInCol); + { + ElementAboveRow2 = &(pElement->NextInCol); pElement = *ElementAboveRow2; - } while (pElement != NULL AND pElement->Row < Row2); + } while (pElement != NULL && pElement->Row < Row2); -/* Place Element1 in Row2. */ + /* Place Element1 in Row2. */ *ElementAboveRow2 = Element1; Element1->NextInCol = pElement; *ElementAboveRow1 =ElementBelowRow1; @@ -2489,26 +2509,27 @@ register ElementPtr pElement; } else { -/* Element2 does exist, and the two elements must be exchanged. */ + /* Element2 does exist, and the two elements must be exchanged. */ if ( ElementBelowRow1->Row == Row2) { -/* Element2 is just below Element1, exchange them. */ + /* Element2 is just below Element1, exchange them. */ Element1->NextInCol = Element2->NextInCol; Element2->NextInCol = Element1; *ElementAboveRow1 = Element2; } else { -/* Element2 is not just below Element1 and must be searched for. */ + /* Element2 is not just below Element1 and must be searched for. */ pElement = ElementBelowRow1; do - { ElementAboveRow2 = &(pElement->NextInCol); + { + ElementAboveRow2 = &(pElement->NextInCol); pElement = *ElementAboveRow2; } while (pElement->Row < Row2); ElementBelowRow2 = Element2->NextInCol; -/* Switch Element1 and Element2. */ + /* Switch Element1 and Element2. */ *ElementAboveRow1 = Element2; Element2->NextInCol = ElementBelowRow1; *ElementAboveRow2 = Element1; @@ -2520,19 +2541,21 @@ register ElementPtr pElement; } else { -/* Element1 does not exist. */ + /* Element1 does not exist. */ ElementBelowRow1 = pElement; -/* Find Element2. */ + /* Find Element2. */ if (ElementBelowRow1->Row != Row2) - { do - { ElementAboveRow2 = &(pElement->NextInCol); + { + do + { + ElementAboveRow2 = &(pElement->NextInCol); pElement = *ElementAboveRow2; } while (pElement->Row < Row2); ElementBelowRow2 = Element2->NextInCol; -/* Move Element2 to Row1. */ + /* Move Element2 to Row1. */ *ElementAboveRow2 = Element2->NextInCol; *ElementAboveRow1 = Element2; Element2->NextInCol = ElementBelowRow1; @@ -2591,38 +2614,41 @@ ExchangeRowElements( Matrix, Col1, Element1, Col2, Element2, Row ) MatrixPtr Matrix; int Col1, Col2, Row; -register ElementPtr Element1, Element2; +ElementPtr Element1, Element2; { ElementPtr *ElementLeftOfCol1, *ElementLeftOfCol2; ElementPtr ElementRightOfCol1, ElementRightOfCol2; -register ElementPtr pElement; + ElementPtr pElement; -/* Begin `ExchangeRowElements'. */ -/* Search to find the ElementLeftOfCol1. */ + /* Begin `ExchangeRowElements'. */ + /* Search to find the ElementLeftOfCol1. */ ElementLeftOfCol1 = &(Matrix->FirstInRow[Row]); pElement = *ElementLeftOfCol1; while (pElement->Col < Col1) - { ElementLeftOfCol1 = &(pElement->NextInRow); + { + ElementLeftOfCol1 = &(pElement->NextInRow); pElement = *ElementLeftOfCol1; } if (Element1 != NULL) - { ElementRightOfCol1 = Element1->NextInRow; + { + ElementRightOfCol1 = Element1->NextInRow; if (Element2 == NULL) { -/* Element2 does not exist, move Element1 to right to Col2. */ - if ( ElementRightOfCol1 != NULL AND ElementRightOfCol1->Col < Col2 ) + /* Element2 does not exist, move Element1 to right to Col2. */ + if ( ElementRightOfCol1 != NULL && ElementRightOfCol1->Col < Col2 ) { -/* Element1 must be removed from linked list and moved. */ + /* Element1 must be removed from linked list and moved. */ *ElementLeftOfCol1 = ElementRightOfCol1; -/* Search Row for Col2. */ + /* Search Row for Col2. */ pElement = ElementRightOfCol1; do - { ElementLeftOfCol2 = &(pElement->NextInRow); + { + ElementLeftOfCol2 = &(pElement->NextInRow); pElement = *ElementLeftOfCol2; - } while (pElement != NULL AND pElement->Col < Col2); + } while (pElement != NULL && pElement->Col < Col2); -/* Place Element1 in Col2. */ + /* Place Element1 in Col2. */ *ElementLeftOfCol2 = Element1; Element1->NextInRow = pElement; *ElementLeftOfCol1 =ElementRightOfCol1; @@ -2631,26 +2657,27 @@ register ElementPtr pElement; } else { -/* Element2 does exist, and the two elements must be exchanged. */ + /* Element2 does exist, and the two elements must be exchanged. */ if ( ElementRightOfCol1->Col == Col2) { -/* Element2 is just right of Element1, exchange them. */ + /* Element2 is just right of Element1, exchange them. */ Element1->NextInRow = Element2->NextInRow; Element2->NextInRow = Element1; *ElementLeftOfCol1 = Element2; } else { -/* Element2 is not just right of Element1 and must be searched for. */ + /* Element2 is not just right of Element1 and must be searched for. */ pElement = ElementRightOfCol1; do - { ElementLeftOfCol2 = &(pElement->NextInRow); + { + ElementLeftOfCol2 = &(pElement->NextInRow); pElement = *ElementLeftOfCol2; } while (pElement->Col < Col2); ElementRightOfCol2 = Element2->NextInRow; -/* Switch Element1 and Element2. */ + /* Switch Element1 and Element2. */ *ElementLeftOfCol1 = Element2; Element2->NextInRow = ElementRightOfCol1; *ElementLeftOfCol2 = Element1; @@ -2662,19 +2689,21 @@ register ElementPtr pElement; } else { -/* Element1 does not exist. */ + /* Element1 does not exist. */ ElementRightOfCol1 = pElement; -/* Find Element2. */ + /* Find Element2. */ if (ElementRightOfCol1->Col != Col2) - { do - { ElementLeftOfCol2 = &(pElement->NextInRow); + { + do + { + ElementLeftOfCol2 = &(pElement->NextInRow); pElement = *ElementLeftOfCol2; } while (pElement->Col < Col2); ElementRightOfCol2 = Element2->NextInRow; -/* Move Element2 to Col1. */ + /* Move Element2 to Col1. */ *ElementLeftOfCol2 = Element2->NextInRow; *ElementLeftOfCol1 = Element2; Element2->NextInRow = ElementRightOfCol1; @@ -2725,19 +2754,19 @@ static void RealRowColElimination( Matrix, pPivot ) MatrixPtr Matrix; -register ElementPtr pPivot; + ElementPtr pPivot; { -#if REAL -register ElementPtr pSub; -register int Row; -register ElementPtr pLower, pUpper; -extern ElementPtr CreateFillin(); + ElementPtr pSub; + int Row; + ElementPtr pLower, pUpper; + extern ElementPtr CreateFillin(); -/* Begin `RealRowColElimination'. */ + /* Begin `RealRowColElimination'. */ -/* Test for zero pivot. */ + /* Test for zero pivot. */ if (ABS(pPivot->Real) == 0.0) - { (void)MatrixIsSingular( Matrix, pPivot->Row ); + { + (void)MatrixIsSingular( Matrix, pPivot->Row ); return; } pPivot->Real = 1.0 / pPivot->Real; @@ -2745,23 +2774,26 @@ extern ElementPtr CreateFillin(); pUpper = pPivot->NextInRow; while (pUpper != NULL) { -/* Calculate upper triangular element. */ + /* Calculate upper triangular element. */ pUpper->Real *= pPivot->Real; pSub = pUpper->NextInCol; pLower = pPivot->NextInCol; while (pLower != NULL) - { Row = pLower->Row; + { + Row = pLower->Row; -/* Find element in row that lines up with current lower triangular element. */ - while (pSub != NULL AND pSub->Row < Row) + /* Find element in row that lines up with current lower triangular element. */ + while (pSub != NULL && pSub->Row < Row) pSub = pSub->NextInCol; -/* Test to see if desired element was not found, if not, create fill-in. */ - if (pSub == NULL OR pSub->Row > Row) - { pSub = CreateFillin( Matrix, Row, pUpper->Col ); + /* Test to see if desired element was not found, if not, create fill-in. */ + if (pSub == NULL || pSub->Row > Row) + { + pSub = CreateFillin( Matrix, Row, pUpper->Col ); if (pSub == NULL) - { Matrix->Error = spNO_MEMORY; + { + Matrix->Error = spNO_MEMORY; return; } } @@ -2772,7 +2804,6 @@ extern ElementPtr CreateFillin(); pUpper = pUpper->NextInRow; } return; -#endif /* REAL */ } @@ -2814,19 +2845,19 @@ static void ComplexRowColElimination( Matrix, pPivot ) MatrixPtr Matrix; -register ElementPtr pPivot; + ElementPtr pPivot; { -#if spCOMPLEX -register ElementPtr pSub; -register int Row; -register ElementPtr pLower, pUpper; -ElementPtr CreateFillin(); + ElementPtr pSub; + int Row; + ElementPtr pLower, pUpper; + ElementPtr CreateFillin(); -/* Begin `ComplexRowColElimination'. */ + /* Begin `ComplexRowColElimination'. */ -/* Test for zero pivot. */ + /* Test for zero pivot. */ if (ELEMENT_MAG(pPivot) == 0.0) - { (void)MatrixIsSingular( Matrix, pPivot->Row ); + { + (void)MatrixIsSingular( Matrix, pPivot->Row ); return; } CMPLX_RECIPROCAL(*pPivot, *pPivot); @@ -2834,29 +2865,32 @@ ElementPtr CreateFillin(); pUpper = pPivot->NextInRow; while (pUpper != NULL) { -/* Calculate upper triangular element. */ -/* Cmplx expr: *pUpper = *pUpper * (1.0 / *pPivot). */ + /* Calculate upper triangular element. */ + /* Cmplx expr: *pUpper = *pUpper * (1.0 / *pPivot). */ CMPLX_MULT_ASSIGN(*pUpper, *pPivot); pSub = pUpper->NextInCol; pLower = pPivot->NextInCol; while (pLower != NULL) - { Row = pLower->Row; + { + Row = pLower->Row; -/* Find element in row that lines up with current lower triangular element. */ - while (pSub != NULL AND pSub->Row < Row) + /* Find element in row that lines up with current lower triangular element. */ + while (pSub != NULL && pSub->Row < Row) pSub = pSub->NextInCol; -/* Test to see if desired element was not found, if not, create fill-in. */ - if (pSub == NULL OR pSub->Row > Row) - { pSub = CreateFillin( Matrix, Row, pUpper->Col ); + /* Test to see if desired element was not found, if not, create fill-in. */ + if (pSub == NULL || pSub->Row > Row) + { + pSub = CreateFillin( Matrix, Row, pUpper->Col ); if (pSub == NULL) - { Matrix->Error = spNO_MEMORY; + { + Matrix->Error = spNO_MEMORY; return; } } -/* Cmplx expr: pElement -= *pUpper * pLower. */ + /* Cmplx expr: pElement -= *pUpper * pLower. */ CMPLX_MULT_SUBT_ASSIGN(*pSub, *pUpper, *pLower); pSub = pSub->NextInCol; pLower = pLower->NextInCol; @@ -2864,7 +2898,6 @@ ElementPtr CreateFillin(); pUpper = pUpper->NextInRow; } return; -#endif /* spCOMPLEX */ } @@ -2900,22 +2933,25 @@ UpdateMarkowitzNumbers( Matrix, pPivot ) MatrixPtr Matrix; ElementPtr pPivot; { -register int Row, Col; -register ElementPtr ColPtr, RowPtr; -register int *MarkoRow = Matrix->MarkowitzRow, *MarkoCol = Matrix->MarkowitzCol; -double Product; + int Row, Col; + ElementPtr ColPtr, RowPtr; + int *MarkoRow = Matrix->MarkowitzRow; + int *MarkoCol = Matrix->MarkowitzCol; + double Product; -/* Begin `UpdateMarkowitzNumbers'. */ + /* Begin `UpdateMarkowitzNumbers'. */ -/* Update Markowitz numbers. */ + /* Update Markowitz numbers. */ for (ColPtr = pPivot->NextInCol; ColPtr != NULL; ColPtr = ColPtr->NextInCol) - { Row = ColPtr->Row; + { + Row = ColPtr->Row; --MarkoRow[Row]; -/* Form Markowitz product while being cautious of overflows. */ - if ((MarkoRow[Row] > LARGEST_SHORT_INTEGER AND MarkoCol[Row] != 0) OR - (MarkoCol[Row] > LARGEST_SHORT_INTEGER AND MarkoRow[Row] != 0)) - { Product = MarkoCol[Row] * MarkoRow[Row]; + /* Form Markowitz product while being cautious of overflows. */ + if ((MarkoRow[Row] > LARGEST_SHORT_INTEGER && MarkoCol[Row] != 0) || + (MarkoCol[Row] > LARGEST_SHORT_INTEGER && MarkoRow[Row] != 0)) + { + Product = MarkoCol[Row] * MarkoRow[Row]; if (Product >= LARGEST_LONG_INTEGER) Matrix->MarkowitzProd[Row] = LARGEST_LONG_INTEGER; else @@ -2926,21 +2962,23 @@ double Product; Matrix->Singletons++; } - for (RowPtr = pPivot->NextInRow; RowPtr != NULL; RowPtr = RowPtr->NextInRow) - { Col = RowPtr->Col; + for (RowPtr = pPivot->NextInRow; + RowPtr != NULL; + RowPtr = RowPtr->NextInRow) { + Col = RowPtr->Col; --MarkoCol[Col]; -/* Form Markowitz product while being cautious of overflows. */ - if ((MarkoRow[Col] > LARGEST_SHORT_INTEGER AND MarkoCol[Col] != 0) OR - (MarkoCol[Col] > LARGEST_SHORT_INTEGER AND MarkoRow[Col] != 0)) - { Product = MarkoCol[Col] * MarkoRow[Col]; + /* Form Markowitz product while being cautious of overflows. */ + if ((MarkoRow[Col] > LARGEST_SHORT_INTEGER && MarkoCol[Col] != 0) || + (MarkoCol[Col] > LARGEST_SHORT_INTEGER && MarkoRow[Col] != 0)) { + Product = MarkoCol[Col] * MarkoRow[Col]; if (Product >= LARGEST_LONG_INTEGER) Matrix->MarkowitzProd[Col] = LARGEST_LONG_INTEGER; else Matrix->MarkowitzProd[Col] = Product; } else Matrix->MarkowitzProd[Col] = MarkoRow[Col] * MarkoCol[Col]; - if ((MarkoCol[Col] == 0) AND (MarkoRow[Col] != 0)) + if ((MarkoCol[Col] == 0) && (MarkoRow[Col] != 0)) Matrix->Singletons++; } return; @@ -2986,36 +3024,38 @@ static ElementPtr CreateFillin( Matrix, Row, Col ) MatrixPtr Matrix; -register int Row; +int Row; int Col; { -register ElementPtr pElement, *ppElementAbove; + ElementPtr pElement, *ppElementAbove; ElementPtr spcCreateElement(); -/* Begin `CreateFillin'. */ + /* Begin `CreateFillin'. */ -/* Find Element above fill-in. */ + /* Find Element above fill-in. */ ppElementAbove = &Matrix->FirstInCol[Col]; pElement = *ppElementAbove; while (pElement != NULL) - { if (pElement->Row < Row) - { ppElementAbove = &pElement->NextInCol; + { + if (pElement->Row < Row) + { + ppElementAbove = &pElement->NextInCol; pElement = *ppElementAbove; } else break; /* while loop */ } -/* End of search, create the element. */ + /* End of search, create the element. */ pElement = spcCreateElement( Matrix, Row, Col, ppElementAbove, YES ); -/* Update Markowitz counts and products. */ + /* Update Markowitz counts and products. */ Matrix->MarkowitzProd[Row] = ++Matrix->MarkowitzRow[Row] * Matrix->MarkowitzCol[Row]; - if ((Matrix->MarkowitzRow[Row] == 1) AND (Matrix->MarkowitzCol[Row] != 0)) + if ((Matrix->MarkowitzRow[Row] == 1) && (Matrix->MarkowitzCol[Row] != 0)) Matrix->Singletons--; Matrix->MarkowitzProd[Col] = ++Matrix->MarkowitzCol[Col] * Matrix->MarkowitzRow[Col]; - if ((Matrix->MarkowitzRow[Col] != 0) AND (Matrix->MarkowitzCol[Col] == 1)) + if ((Matrix->MarkowitzRow[Col] != 0) && (Matrix->MarkowitzCol[Col] == 1)) Matrix->Singletons--; return pElement; @@ -3050,7 +3090,7 @@ MatrixIsSingular( Matrix, Step ) MatrixPtr Matrix; int Step; { -/* Begin `MatrixIsSingular'. */ + /* Begin `MatrixIsSingular'. */ Matrix->SingularRow = Matrix->IntToExtRowMap[ Step ]; Matrix->SingularCol = Matrix->IntToExtColMap[ Step ]; @@ -3064,7 +3104,7 @@ ZeroPivot( Matrix, Step ) MatrixPtr Matrix; int Step; { -/* Begin `ZeroPivot'. */ + /* Begin `ZeroPivot'. */ Matrix->SingularRow = Matrix->IntToExtRowMap[ Step ]; Matrix->SingularCol = Matrix->IntToExtColMap[ Step ]; @@ -3093,13 +3133,14 @@ int Step; { int I; -/* Begin `WriteStatus'. */ + /* Begin `WriteStatus'. */ printf("Step = %1d ", Step); printf("Pivot found at %1d,%1d using ", Matrix->PivotsOriginalRow, Matrix->PivotsOriginalCol); switch(Matrix->PivotSelectionMethod) - { case 's': printf("SearchForSingleton\n"); break; + { + case 's': printf("SearchForSingleton\n"); break; case 'q': printf("QuicklySearchDiagonal\n"); break; case 'd': printf("SearchDiagonal\n"); break; case 'e': printf("SearchEntireMatrix\n"); break; @@ -3142,8 +3183,6 @@ int I; printf("%2d ", Matrix->ExtToIntColMap[I]); printf("\n\n"); -/* spPrint((char *)Matrix, NO, YES); */ - return; } diff --git a/src/maths/sparse/spoutput.c b/src/maths/sparse/spoutput.c index b42b58104..9511c72fd 100644 --- a/src/maths/sparse/spoutput.c +++ b/src/maths/sparse/spoutput.c @@ -45,6 +45,7 @@ * spDefs.h * Matrix type and macro definitions for the sparse matrix routines. */ +#include #define spINSIDE_SPARSE #include "spconfig.h" @@ -130,30 +131,30 @@ int Printer_Width = PRINTER_WIDTH; */ void -spPrint( eMatrix, PrintReordered, Data, Header ) - -char *eMatrix; -int PrintReordered, Data, Header; +spPrint(void *eMatrix, int PrintReordered, int Data, int Header) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register int J = 0; -int I, Row, Col, Size, Top, StartCol = 1, StopCol, Columns, ElementCount = 0; -double Magnitude, SmallestDiag, SmallestElement; -double LargestElement = 0.0, LargestDiag = 0.0; -ElementPtr pElement, *pImagElements; -int *PrintOrdToIntRowMap, *PrintOrdToIntColMap; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + int J = 0; + int I, Row, Col, Size, Top; + int StartCol = 1, StopCol, Columns, ElementCount = 0; + double Magnitude, SmallestDiag, SmallestElement; + double LargestElement = 0.0, LargestDiag = 0.0; + ElementPtr pElement, *pImagElements; + int *PrintOrdToIntRowMap, *PrintOrdToIntColMap; -/* Begin `spPrint'. */ - ASSERT( IS_SPARSE( Matrix ) ); + /* Begin `spPrint'. */ + assert( IS_SPARSE( Matrix ) ); Size = Matrix->Size; CALLOC(pImagElements, ElementPtr, Printer_Width / 10 + 1); if ( pImagElements == NULL) - { Matrix->Error = spNO_MEMORY; + { + Matrix->Error = spNO_MEMORY; FREE(pImagElements); return; } -/* Create a packed external to internal row and column translation array. */ + /* Create a packed external to internal row and column translation + array. */ # if TRANSLATE Top = Matrix->AllocatedExtSize; #else @@ -161,37 +162,43 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap; #endif CALLOC( PrintOrdToIntRowMap, int, Top + 1 ); if ( PrintOrdToIntRowMap == NULL) - { Matrix->Error = spNO_MEMORY; + { + Matrix->Error = spNO_MEMORY; FREE(pImagElements); return; } CALLOC( PrintOrdToIntColMap, int, Top + 1 ); if (PrintOrdToIntColMap == NULL) - { Matrix->Error = spNO_MEMORY; + { + Matrix->Error = spNO_MEMORY; FREE(pImagElements); FREE(PrintOrdToIntRowMap); return; } for (I = 1; I <= Size; I++) - { PrintOrdToIntRowMap[ Matrix->IntToExtRowMap[I] ] = I; + { + PrintOrdToIntRowMap[ Matrix->IntToExtRowMap[I] ] = I; PrintOrdToIntColMap[ Matrix->IntToExtColMap[I] ] = I; } -/* Pack the arrays. */ + /* Pack the arrays. */ for (J = 1, I = 1; I <= Top; I++) - { if (PrintOrdToIntRowMap[I] != 0) + { + if (PrintOrdToIntRowMap[I] != 0) PrintOrdToIntRowMap[ J++ ] = PrintOrdToIntRowMap[ I ]; } for (J = 1, I = 1; I <= Top; I++) - { if (PrintOrdToIntColMap[I] != 0) + { + if (PrintOrdToIntColMap[I] != 0) PrintOrdToIntColMap[ J++ ] = PrintOrdToIntColMap[ I ]; } -/* Print header. */ + /* Print header. */ if (Header) - { printf("MATRIX SUMMARY\n\n"); + { + printf("MATRIX SUMMARY\n\n"); printf("Size of matrix = %1d x %1d.\n", Size, Size); - if ( Matrix->Reordered AND PrintReordered ) + if ( Matrix->Reordered && PrintReordered ) printf("Matrix has been reordered.\n"); putchar('\n'); @@ -204,29 +211,30 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap; SmallestDiag = SmallestElement; } -/* Determine how many columns to use. */ + /* Determine how many columns to use. */ Columns = Printer_Width; if (Header) Columns -= 5; if (Data) Columns = (Columns+1) / 10; -/* - * Print matrix by printing groups of complete columns until all the columns - * are printed. - */ + /* Print matrix by printing groups of complete columns until all + * the columns are printed. */ J = 0; while ( J <= Size ) - -/* Calculatestrchr of last column to printed in this group. */ - { StopCol = StartCol + Columns - 1; + { + /* Calculatestrchr of last column to printed in this group. */ + StopCol = StartCol + Columns - 1; if (StopCol > Size) StopCol = Size; -/* Label the columns. */ + /* Label the columns. */ if (Header) - { if (Data) - { printf(" "); + { + if (Data) + { + printf(" "); for (I = StartCol; I <= StopCol; I++) - { if (PrintReordered) + { + if (PrintReordered) Col = I; else Col = PrintOrdToIntColMap[I]; @@ -235,64 +243,70 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap; printf("\n\n"); } else - { if (PrintReordered) + { + if (PrintReordered) printf("Columns %1d to %1d.\n",StartCol,StopCol); else - { printf("Columns %1d to %1d.\n", - Matrix->IntToExtColMap[ PrintOrdToIntColMap[StartCol] ], - Matrix->IntToExtColMap[ PrintOrdToIntColMap[StopCol] ]); + { + printf("Columns %1d to %1d.\n", + Matrix->IntToExtColMap[ PrintOrdToIntColMap[StartCol] ], + Matrix->IntToExtColMap[ PrintOrdToIntColMap[StopCol] ]); } } } -/* Print every row ... */ + /* Print every row ... */ for (I = 1; I <= Size; I++) - { if (PrintReordered) + { + if (PrintReordered) Row = I; else Row = PrintOrdToIntRowMap[I]; if (Header) - { if (PrintReordered AND NOT Data) + { + if (PrintReordered && !Data) printf("%4d", I); else printf("%4d", Matrix->IntToExtRowMap[ Row ]); - if (NOT Data) putchar(' '); + if (!Data) putchar(' '); } -/* ... in each column of the group. */ + /* ... in each column of the group. */ for (J = StartCol; J <= StopCol; J++) - { if (PrintReordered) + { + if (PrintReordered) Col = J; else Col = PrintOrdToIntColMap[J]; pElement = Matrix->FirstInCol[Col]; - while(pElement != NULL AND pElement->Row != Row) + while(pElement != NULL && pElement->Row != Row) pElement = pElement->NextInCol; if (Data) pImagElements[J - StartCol] = pElement; if (pElement != NULL) - -/* Case where element exists */ - { if (Data) + { + /* Case where element exists */ + if (Data) printf(" %9.3g", (double)pElement->Real); else putchar('x'); -/* Update status variables */ + /* Update status variables */ if ( (Magnitude = ELEMENT_MAG(pElement)) > LargestElement ) LargestElement = Magnitude; - if ((Magnitude < SmallestElement) AND (Magnitude != 0.0)) + if ((Magnitude < SmallestElement) && (Magnitude != 0.0)) SmallestElement = Magnitude; ElementCount++; } -/* Case where element is structurally zero */ + /* Case where element is structurally zero */ else - { if (Data) + { + if (Data) printf(" ..."); else putchar('.'); @@ -300,55 +314,61 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap; } putchar('\n'); -#if spCOMPLEX - if (Matrix->Complex AND Data) - { printf(" "); + if (Matrix->Complex && Data) + { + printf(" "); for (J = StartCol; J <= StopCol; J++) - { if (pImagElements[J - StartCol] != NULL) - { printf(" %8.2gj", + { + if (pImagElements[J - StartCol] != NULL) + { + printf(" %8.2gj", (double)pImagElements[J-StartCol]->Imag); } else printf(" "); } putchar('\n'); } -#endif /* spCOMPLEX */ } -/* Calculatestrchr of first column in next group. */ + /* Calculatestrchr of first column in next group. */ StartCol = StopCol; StartCol++; putchar('\n'); } if (Header) - { printf("\nLargest element in matrix = %-1.4g.\n", LargestElement); + { + printf("\nLargest element in matrix = %-1.4g.\n", LargestElement); printf("Smallest element in matrix = %-1.4g.\n", SmallestElement); -/* Search for largest and smallest diagonal values */ + /* Search for largest and smallest diagonal values */ for (I = 1; I <= Size; I++) - { if (Matrix->Diag[I] != NULL) - { Magnitude = ELEMENT_MAG( Matrix->Diag[I] ); + { + if (Matrix->Diag[I] != NULL) + { + Magnitude = ELEMENT_MAG( Matrix->Diag[I] ); if ( Magnitude > LargestDiag ) LargestDiag = Magnitude; if ( Magnitude < SmallestDiag ) SmallestDiag = Magnitude; } } - /* Print the largest and smallest diagonal values */ + /* Print the largest and smallest diagonal values */ if ( Matrix->Factored ) - { printf("\nLargest diagonal element = %-1.4g.\n", LargestDiag); + { + printf("\nLargest diagonal element = %-1.4g.\n", LargestDiag); printf("Smallest diagonal element = %-1.4g.\n", SmallestDiag); } else - { printf("\nLargest pivot element = %-1.4g.\n", LargestDiag); + { + printf("\nLargest pivot element = %-1.4g.\n", LargestDiag); printf("Smallest pivot element = %-1.4g.\n", SmallestDiag); } - /* Calculate and print sparsity and number of fill-ins created. */ + /* Calculate and print sparsity and number of fill-ins created. */ printf("\nDensity = %2.2f%%.\n", ((double)(ElementCount * 100)) / - ((double)(Size * Size))); + ((double)(Size * Size))); printf("Number of originals = %1d.\n", Matrix->Originals); - if (NOT Matrix->NeedsOrdering) + if (!Matrix->NeedsOrdering) printf("Number of fill-ins = %1d.\n", Matrix->Fillins); } putchar('\n'); @@ -387,14 +407,14 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap; * Name of file into which matrix is to be written. * Label (char *) * String that is transferred to file and is used as a label. - * Reordered (BOOLEAN) + * Reordered (int) * Specifies whether matrix should be output in reordered form, * or in original order. - * Data (BOOLEAN) + * Data (int) * Indicates that the element values should be output along with * the indices for each element. This parameter must be TRUE if * matrix is to be read by the sparse test program. - * Header (BOOLEAN) + * Header (int) * Indicates that header is desired. This parameter must be TRUE if * matrix is to be read by the sparse test program. * @@ -412,113 +432,124 @@ int *PrintOrdToIntRowMap, *PrintOrdToIntColMap; */ int -spFileMatrix( eMatrix, File, Label, Reordered, Data, Header ) - -char *eMatrix, *Label, *File; -int Reordered, Data, Header; +spFileMatrix(void *eMatrix, char *File, char *Label, int Reordered, + int Data, int Header) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register int I, Size; -register ElementPtr pElement; -int Row, Col, Err; -FILE *pMatrixFile, *fopen(); + MatrixPtr Matrix = (MatrixPtr)eMatrix; + int I, Size; + ElementPtr pElement; + int Row, Col, Err; + FILE *pMatrixFile; -/* Begin `spFileMatrix'. */ - ASSERT( IS_SPARSE( Matrix ) ); + /* Begin `spFileMatrix'. */ + assert( IS_SPARSE( Matrix ) ); -/* Open file matrix file in write mode. */ + /* Open file matrix file in write mode. */ if ((pMatrixFile = fopen(File, "w")) == NULL) return 0; -/* Output header. */ + /* Output header. */ Size = Matrix->Size; if (Header) - { if (Matrix->Factored AND Data) - { Err = fprintf - ( pMatrixFile, - "Warning : The following matrix is factored in to LU form.\n" - ); - if (Err < 0) return 0; + { + if (Matrix->Factored && Data) + { + Err = fprintf(pMatrixFile, + "Warning : The following matrix is " + "factored in to LU form.\n"); + if (Err < 0) + return 0; } - if (fprintf(pMatrixFile, "%s\n", Label) < 0) return 0; + if (fprintf(pMatrixFile, "%s\n", Label) < 0) + return 0; Err = fprintf( pMatrixFile, "%d\t%s\n", Size, - (Matrix->Complex ? "complex" : "real")); - if (Err < 0) return 0; + (Matrix->Complex ? "complex" : "real")); + if (Err < 0) + return 0; } -/* Output matrix. */ - if (NOT Data) - { for (I = 1; I <= Size; I++) - { pElement = Matrix->FirstInCol[I]; + /* Output matrix. */ + if (!Data) + { + for (I = 1; I <= Size; I++) + { + pElement = Matrix->FirstInCol[I]; while (pElement != NULL) - { if (Reordered) - { Row = pElement->Row; + { + if (Reordered) + { + Row = pElement->Row; Col = I; } else - { Row = Matrix->IntToExtRowMap[pElement->Row]; + { + Row = Matrix->IntToExtRowMap[pElement->Row]; Col = Matrix->IntToExtColMap[I]; } pElement = pElement->NextInCol; if (fprintf(pMatrixFile, "%d\t%d\n", Row, Col) < 0) return 0; } } -/* Output terminator, a line of zeros. */ + /* Output terminator, a line of zeros. */ if (Header) if (fprintf(pMatrixFile, "0\t0\n") < 0) return 0; } -#if spCOMPLEX - if (Data AND Matrix->Complex) - { for (I = 1; I <= Size; I++) - { pElement = Matrix->FirstInCol[I]; + if (Data && Matrix->Complex) + { + for (I = 1; I <= Size; I++) + { + pElement = Matrix->FirstInCol[I]; while (pElement != NULL) - { if (Reordered) - { Row = pElement->Row; + { + if (Reordered) + { + Row = pElement->Row; Col = I; } else - { Row = Matrix->IntToExtRowMap[pElement->Row]; + { + Row = Matrix->IntToExtRowMap[pElement->Row]; Col = Matrix->IntToExtColMap[I]; } Err = fprintf - ( pMatrixFile,"%d\t%d\t%-.15g\t%-.15g\n", - Row, Col, (double)pElement->Real, (double)pElement->Imag - ); + ( pMatrixFile,"%d\t%d\t%-.15g\t%-.15g\n", + Row, Col, (double)pElement->Real, (double)pElement->Imag + ); if (Err < 0) return 0; pElement = pElement->NextInCol; } } -/* Output terminator, a line of zeros. */ + /* Output terminator, a line of zeros. */ if (Header) if (fprintf(pMatrixFile,"0\t0\t0.0\t0.0\n") < 0) return 0; } -#endif /* spCOMPLEX */ -#if REAL - if (Data AND NOT Matrix->Complex) - { for (I = 1; I <= Size; I++) - { pElement = Matrix->FirstInCol[I]; + if (Data && !Matrix->Complex) + { + for (I = 1; I <= Size; I++) + { + pElement = Matrix->FirstInCol[I]; while (pElement != NULL) - { Row = Matrix->IntToExtRowMap[pElement->Row]; + { + Row = Matrix->IntToExtRowMap[pElement->Row]; Col = Matrix->IntToExtColMap[I]; Err = fprintf - ( pMatrixFile,"%d\t%d\t%-.15g\n", - Row, Col, (double)pElement->Real - ); + ( pMatrixFile,"%d\t%d\t%-.15g\n", + Row, Col, (double)pElement->Real + ); if (Err < 0) return 0; pElement = pElement->NextInCol; } } -/* Output terminator, a line of zeros. */ + /* Output terminator, a line of zeros. */ if (Header) if (fprintf(pMatrixFile,"0\t0\t0.0\n") < 0) return 0; } -#endif /* REAL */ -/* Close file. */ + /* Close file. */ if (fclose(pMatrixFile) < 0) return 0; return 1; } @@ -547,11 +578,9 @@ FILE *pMatrixFile, *fopen(); * File (char *) * Name of file into which matrix is to be written. * RHS (RealNumber []) - * Right-hand side vector. This is only the real portion if - * spSEPARATED_COMPLEX_VECTORS is TRUE. + * Right-hand side vector, real portion * iRHS (RealNumber []) - * Right-hand side vector, imaginary portion. Not necessary if matrix - * is real or if spSEPARATED_COMPLEX_VECTORS is set FALSE. + * Right-hand side vector, imaginary portion. * * >>> Local variables: * pMatrixFile (FILE *) @@ -559,86 +588,47 @@ FILE *pMatrixFile, *fopen(); * Size (int) * The size of the matrix. * - * >>> Obscure Macros - * IMAG_RHS - * Replaces itself with `, iRHS' if the options spCOMPLEX and - * spSEPARATED_COMPLEX_VECTORS are set, otherwise it disappears - * without a trace. */ int -spFileVector( eMatrix, File, RHS IMAG_RHS ) - -char *eMatrix, *File; -RealVector RHS IMAG_RHS; +spFileVector(void *eMatrix, char *File, RealVector RHS, RealVector iRHS) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register int I, Size, Err; -FILE *pMatrixFile; -FILE *fopen(); + MatrixPtr Matrix = (MatrixPtr)eMatrix; + int I, Size, Err; + FILE *pMatrixFile; + FILE *fopen(); -/* Begin `spFileVector'. */ - ASSERT( IS_SPARSE( Matrix ) AND RHS != NULL) + /* Begin `spFileVector'. */ + assert( IS_SPARSE( Matrix ) && RHS != NULL); -/* Open File in append mode. */ - if ((pMatrixFile = fopen(File,"a")) == NULL) - return 0; + /* Open File in append mode. */ + pMatrixFile = fopen(File,"a"); + if (pMatrixFile == NULL) + return 0; -/* Correct array pointers for ARRAY_OFFSET. */ -#if NOT ARRAY_OFFSET -#if spCOMPLEX - if (Matrix->Complex) - { -#if spSEPARATED_COMPLEX_VECTORS - ASSERT(iRHS != NULL) - --RHS; - --iRHS; -#else - RHS -= 2; -#endif - } - else -#endif /* spCOMPLEX */ - --RHS; -#endif /* NOT ARRAY_OFFSET */ - - -/* Output vector. */ + /* Output vector. */ Size = Matrix->Size; -#if spCOMPLEX if (Matrix->Complex) { -#if spSEPARATED_COMPLEX_VECTORS for (I = 1; I <= Size; I++) - { Err = fprintf - ( pMatrixFile, "%-.15g\t%-.15g\n", - (double)RHS[I], (double)iRHS[I] - ); + { + Err = fprintf + ( pMatrixFile, "%-.15g\t%-.15g\n", + (double)RHS[I], (double)iRHS[I] + ); if (Err < 0) return 0; } -#else - for (I = 1; I <= Size; I++) - { Err = fprintf - ( pMatrixFile, "%-.15g\t%-.15g\n", - (double)RHS[2*I], (double)RHS[2*I+1] - ); - if (Err < 0) return 0; - } -#endif } -#endif /* spCOMPLEX */ -#if REAL AND spCOMPLEX else -#endif -#if REAL - { for (I = 1; I <= Size; I++) - { if (fprintf(pMatrixFile, "%-.15g\n", (double)RHS[I]) < 0) + { + for (I = 1; I <= Size; I++) + { + if (fprintf(pMatrixFile, "%-.15g\n", (double)RHS[I]) < 0) return 0; } } -#endif /* REAL */ -/* Close file. */ + /* Close file. */ if (fclose(pMatrixFile) < 0) return 0; return 1; } @@ -688,27 +678,25 @@ FILE *fopen(); */ int -spFileStats( eMatrix, File, Label ) - -char *eMatrix, *File, *Label; +spFileStats(void *eMatrix, char *File, char *Label) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register int Size, I; -register ElementPtr pElement; -int NumberOfElements; -RealNumber Data, LargestElement, SmallestElement; -FILE *pStatsFile, *fopen(); + MatrixPtr Matrix = (MatrixPtr)eMatrix; + int Size, I; + ElementPtr pElement; + int NumberOfElements; + RealNumber Data, LargestElement, SmallestElement; + FILE *pStatsFile, *fopen(); -/* Begin `spFileStats'. */ - ASSERT( IS_SPARSE( Matrix ) ); + /* Begin `spFileStats'. */ + assert( IS_SPARSE( Matrix ) ); -/* Open File in append mode. */ + /* Open File in append mode. */ if ((pStatsFile = fopen(File, "a")) == NULL) return 0; -/* Output statistics. */ + /* Output statistics. */ Size = Matrix->Size; - if (NOT Matrix->Factored) + if (!Matrix->Factored) fprintf(pStatsFile, "Matrix has not been factored.\n"); fprintf(pStatsFile, "||| Starting new matrix |||\n"); fprintf(pStatsFile, "%s\n", Label); @@ -718,19 +706,21 @@ FILE *pStatsFile, *fopen(); fprintf(pStatsFile, "Matrix is real.\n"); fprintf(pStatsFile," Size = %d\n",Size); -/* Search matrix. */ + /* Search matrix. */ NumberOfElements = 0; LargestElement = 0.0; SmallestElement = LARGEST_REAL; for (I = 1; I <= Size; I++) - { pElement = Matrix->FirstInCol[I]; + { + pElement = Matrix->FirstInCol[I]; while (pElement != NULL) - { NumberOfElements++; + { + NumberOfElements++; Data = ELEMENT_MAG(pElement); if (Data > LargestElement) LargestElement = Data; - if (Data < SmallestElement AND Data != 0.0) + if (Data < SmallestElement && Data != 0.0) SmallestElement = Data; pElement = pElement->NextInCol; } @@ -738,7 +728,7 @@ FILE *pStatsFile, *fopen(); SmallestElement = MIN( SmallestElement, LargestElement ); -/* Output remaining statistics. */ + /* Output remaining statistics. */ fprintf(pStatsFile, " Initial number of elements = %d\n", NumberOfElements - Matrix->Fillins); fprintf(pStatsFile, @@ -758,7 +748,7 @@ FILE *pStatsFile, *fopen(); fprintf(pStatsFile," Largest Element = %e\n", LargestElement); fprintf(pStatsFile," Smallest Element = %e\n\n\n", SmallestElement); -/* Close file. */ + /* Close file. */ (void)fclose(pStatsFile); return 1; } diff --git a/src/maths/sparse/spsmp.c b/src/maths/sparse/spsmp.c index 972a752a4..eab788435 100644 --- a/src/maths/sparse/spsmp.c +++ b/src/maths/sparse/spsmp.c @@ -47,12 +47,10 @@ * To be compatible with SPICE, the following Sparse compiler options * (in spConfig.h) should be set as shown below: * - * REAL YES * EXPANDABLE YES * TRANSLATE NO * INITIALIZE NO or YES, YES for use with test prog. * DIAGONAL_PIVOTING YES - * ARRAY_OFFSET YES * MODIFIED_MARKOWITZ NO * DELETE NO * STRIP NO @@ -66,11 +64,7 @@ * STABILITY NO * CONDITION NO * PSEUDOCONDITION NO - * FORTRAN NO * DEBUG YES - * spCOMPLEX 1 - * spSEPARATED_COMPLEX_VECTORS 1 - * spCOMPATIBILITY 0 * * spREAL double */ @@ -100,63 +94,56 @@ * Spice3's matrix macro definitions. */ -#include "ngspice.h" +#include +#include + +#include #include -#include "spmatrix.h" -#include "smpdefs.h" + +#include +#include + #include "spdefs.h" -/* #define NO 0 */ -/* #define YES 1 */ -#ifdef __STDC__ -static void LoadGmin( char * /*eMatrix*/, double /*Gmin*/ ); -#else -static void LoadGmin( ); -#endif +static void LoadGmin(SMPmatrix *eMatrix, double Gmin); + /* * SMPaddElt() */ int -SMPaddElt( Matrix, Row, Col, Value ) -SMPmatrix *Matrix; -int Row, Col; -double Value; +SMPaddElt(SMPmatrix *Matrix, int Row, int Col, double Value) { - *spGetElement( (char *)Matrix, Row, Col ) = Value; - return spError( (char *)Matrix ); + *spGetElement( (void *)Matrix, Row, Col ) = Value; + return spError( (void *)Matrix ); } /* * SMPmakeElt() */ double * -SMPmakeElt( Matrix, Row, Col ) -SMPmatrix *Matrix; -int Row, Col; +SMPmakeElt(SMPmatrix *Matrix, int Row, int Col) { - return spGetElement( (char *)Matrix, Row, Col ); + return spGetElement( (void *)Matrix, Row, Col ); } /* * SMPcClear() */ void -SMPcClear( Matrix ) -SMPmatrix *Matrix; +SMPcClear(SMPmatrix *Matrix) { - spClear( (char *)Matrix ); + spClear( (void *)Matrix ); } /* * SMPclear() */ void -SMPclear( Matrix ) -SMPmatrix *Matrix; +SMPclear(SMPmatrix *Matrix) { - spClear( (char *)Matrix ); + spClear( (void *)Matrix ); } /* @@ -164,12 +151,10 @@ SMPmatrix *Matrix; */ /*ARGSUSED*/ int -SMPcLUfac( Matrix, PivTol ) -SMPmatrix *Matrix; -double PivTol; +SMPcLUfac(SMPmatrix *Matrix, double PivTol) { - spSetComplex( (char *)Matrix ); - return spFactor( (char *)Matrix ); + spSetComplex( (void *)Matrix ); + return spFactor( (void *)Matrix ); } /* @@ -177,27 +162,23 @@ double PivTol; */ /*ARGSUSED*/ int -SMPluFac( Matrix, PivTol, Gmin ) -SMPmatrix *Matrix; -double PivTol, Gmin; +SMPluFac(SMPmatrix *Matrix, double PivTol, double Gmin) { - spSetReal( (char *)Matrix ); - LoadGmin( (char *)Matrix, Gmin ); - return spFactor( (char *)Matrix ); + spSetReal( (void *)Matrix ); + LoadGmin( (void *)Matrix, Gmin ); + return spFactor( (void *)Matrix ); } /* * SMPcReorder() */ int -SMPcReorder( Matrix, PivTol, PivRel, NumSwaps ) -SMPmatrix *Matrix; -double PivTol, PivRel; -int *NumSwaps; +SMPcReorder(SMPmatrix *Matrix, double PivTol, double PivRel, + int *NumSwaps) { *NumSwaps = 1; - spSetComplex( (char *)Matrix ); - return spOrderAndFactor( (char *)Matrix, (spREAL*)NULL, + spSetComplex( (void *)Matrix ); + return spOrderAndFactor( (void *)Matrix, (spREAL*)NULL, (spREAL)PivRel, (spREAL)PivTol, YES ); } @@ -205,13 +186,11 @@ int *NumSwaps; * SMPreorder() */ int -SMPreorder( Matrix, PivTol, PivRel, Gmin ) -SMPmatrix *Matrix; -double PivTol, PivRel, Gmin; +SMPreorder(SMPmatrix *Matrix, double PivTol, double PivRel, double Gmin) { - spSetReal( (char *)Matrix ); - LoadGmin( (char *)Matrix, Gmin ); - return spOrderAndFactor( (char *)Matrix, (spREAL*)NULL, + spSetReal( (void *)Matrix ); + LoadGmin( (void *)Matrix, Gmin ); + return spOrderAndFactor( (void *)Matrix, (spREAL*)NULL, (spREAL)PivRel, (spREAL)PivTol, YES ); } @@ -219,53 +198,47 @@ double PivTol, PivRel, Gmin; * SMPcaSolve() */ void -SMPcaSolve( Matrix, RHS, iRHS, Spare, iSpare) -SMPmatrix *Matrix; -double RHS[], iRHS[], Spare[], iSpare[]; +SMPcaSolve(SMPmatrix *Matrix, double RHS[], double iRHS[], + double Spare[], double iSpare[]) { - spSolveTransposed( (char *)Matrix, RHS, RHS, iRHS, iRHS ); + spSolveTransposed( (void *)Matrix, RHS, RHS, iRHS, iRHS ); } /* * SMPcSolve() */ void -SMPcSolve( Matrix, RHS, iRHS, Spare, iSpare) -SMPmatrix *Matrix; -double RHS[], iRHS[], Spare[], iSpare[]; +SMPcSolve(SMPmatrix *Matrix, double RHS[], double iRHS[], + double Spare[], double iSpare[]) { - spSolve( (char *)Matrix, RHS, RHS, iRHS, iRHS ); + spSolve( (void *)Matrix, RHS, RHS, iRHS, iRHS ); } /* * SMPsolve() */ void -SMPsolve( Matrix, RHS, Spare ) -SMPmatrix *Matrix; -double RHS[], Spare[]; +SMPsolve(SMPmatrix *Matrix, double RHS[], double Spare[]) { - spSolve( (char *)Matrix, RHS, RHS, (spREAL*)NULL, (spREAL*)NULL ); + spSolve( (void *)Matrix, RHS, RHS, (spREAL*)NULL, (spREAL*)NULL ); } /* * SMPmatSize() */ int -SMPmatSize( Matrix ) -SMPmatrix *Matrix; +SMPmatSize(SMPmatrix *Matrix) { - return spGetSize( (char *)Matrix, 1 ); + return spGetSize( (void *)Matrix, 1 ); } /* * SMPnewMatrix() */ int -SMPnewMatrix( pMatrix ) -SMPmatrix **pMatrix; +SMPnewMatrix(SMPmatrix **pMatrix) { -int Error; + int Error; *pMatrix = (SMPmatrix *)spCreate( 0, 1, &Error ); return Error; } @@ -274,21 +247,19 @@ int Error; * SMPdestroy() */ void -SMPdestroy( Matrix ) -SMPmatrix *Matrix; +SMPdestroy(SMPmatrix *Matrix) { - spDestroy( (char *)Matrix ); + spDestroy( (void *)Matrix ); } /* * SMPpreOrder() */ int -SMPpreOrder( Matrix ) -SMPmatrix *Matrix; +SMPpreOrder(SMPmatrix *Matrix) { - spMNA_Preorder( (char *)Matrix ); - return spError( (char *)Matrix ); + spMNA_Preorder( (void *)Matrix ); + return spError( (void *)Matrix ); } /* @@ -296,22 +267,18 @@ SMPmatrix *Matrix; */ /*ARGSUSED*/ void -SMPprint( Matrix, File ) -SMPmatrix *Matrix; -FILE *File; +SMPprint(SMPmatrix *Matrix, FILE *File) { - spPrint( (char *)Matrix, 0, 1, 1 ); + spPrint( (void *)Matrix, 0, 1, 1 ); } /* * SMPgetError() */ void -SMPgetError( Matrix, Col, Row) -SMPmatrix *Matrix; -int *Row, *Col; +SMPgetError(SMPmatrix *Matrix, int *Col, int *Row) { - spWhereSingular( (char *)Matrix, Row, Col ); + spWhereSingular( (void *)Matrix, Row, Col ); } /* @@ -319,29 +286,23 @@ int *Row, *Col; * note: obsolete for Spice3d2 and later */ int -SMPcProdDiag( Matrix, pMantissa, pExponent) -SMPmatrix *Matrix; -SPcomplex *pMantissa; -int *pExponent; +SMPcProdDiag(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent) { - spDeterminant( (char *)Matrix, pExponent, &(pMantissa->real), + spDeterminant( (void *)Matrix, pExponent, &(pMantissa->real), &(pMantissa->imag) ); - return spError( (char *)Matrix ); + return spError( (void *)Matrix ); } /* * SMPcDProd() */ int -SMPcDProd( Matrix, pMantissa, pExponent) -SMPmatrix *Matrix; -SPcomplex *pMantissa; -int *pExponent; +SMPcDProd(SMPmatrix *Matrix, SPcomplex *pMantissa, int *pExponent) { double re, im, x, y, z; int p; - spDeterminant( (char *)Matrix, &p, &re, &im); + spDeterminant( (void *)Matrix, &p, &re, &im); #ifndef M_LN2 #define M_LN2 0.69314718055994530942 @@ -411,7 +372,7 @@ int *pExponent; printf("Determinant 10->2: (%20g,%20g)^%d\n", pMantissa->real, pMantissa->imag, *pExponent); #endif - return spError( (char *)Matrix ); + return spError( (void *)Matrix ); } @@ -433,17 +394,15 @@ int *pExponent; */ static void -LoadGmin( eMatrix, Gmin ) -char *eMatrix; -register double Gmin; +LoadGmin(SMPmatrix *eMatrix, double Gmin) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register int I; -register ArrayOfElementPtrs Diag; -register ElementPtr diag; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + int I; + ArrayOfElementPtrs Diag; + ElementPtr diag; -/* Begin `LoadGmin'. */ - ASSERT( IS_SPARSE( Matrix ) ); + /* Begin `LoadGmin'. */ + assert( IS_SPARSE( Matrix ) ); if (Gmin != 0.0) { Diag = Matrix->Diag; @@ -468,17 +427,13 @@ register ElementPtr diag; */ SMPelement * -SMPfindElt( eMatrix, Row, Col, CreateIfMissing ) - -SMPmatrix *eMatrix; -int Row, Col; -int CreateIfMissing; +SMPfindElt(SMPmatrix *eMatrix, int Row, int Col, int CreateIfMissing) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -ElementPtr Element; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + ElementPtr Element; -/* Begin `SMPfindElt'. */ - ASSERT( IS_SPARSE( Matrix ) ); + /* Begin `SMPfindElt'. */ + assert( IS_SPARSE( Matrix ) ); Row = Matrix->ExtToIntRowMap[Row]; Col = Matrix->ExtToIntColMap[Col]; Element = Matrix->FirstInCol[Col]; @@ -492,10 +447,9 @@ ElementPtr Element; * SMPcZeroCol() */ int -SMPcZeroCol( Matrix, Col ) -MatrixPtr Matrix; -int Col; +SMPcZeroCol(SMPmatrix *eMatrix, int Col) { + MatrixPtr Matrix = (MatrixPtr)eMatrix; ElementPtr Element; Col = Matrix->ExtToIntColMap[Col]; @@ -508,17 +462,16 @@ int Col; Element->Imag = 0.0; } - return spError( (char *)Matrix ); + return spError( (void *)Matrix ); } /* * SMPcAddCol() */ int -SMPcAddCol( Matrix, Accum_Col, Addend_Col ) -MatrixPtr Matrix; -int Accum_Col, Addend_Col; +SMPcAddCol(SMPmatrix *eMatrix, int Accum_Col, int Addend_Col) { + MatrixPtr Matrix = (MatrixPtr)eMatrix; ElementPtr Accum, Addend, *Prev; Accum_Col = Matrix->ExtToIntColMap[Accum_Col]; @@ -541,17 +494,16 @@ int Accum_Col, Addend_Col; Addend = Addend->NextInCol; } - return spError( (char *)Matrix ); + return spError( (void *)Matrix ); } /* * SMPzeroRow() */ int -SMPzeroRow( Matrix, Row ) -MatrixPtr Matrix; -int Row; +SMPzeroRow(SMPmatrix *eMatrix, int Row) { + MatrixPtr Matrix = (MatrixPtr)eMatrix; ElementPtr Element; Row = Matrix->ExtToIntColMap[Row]; @@ -559,8 +511,7 @@ int Row; if (Matrix->RowsLinked == NO) spcLinkRows(Matrix); -#if spCOMPLEX - if (Matrix->PreviousMatrixWasComplex OR Matrix->Complex) { + if (Matrix->PreviousMatrixWasComplex || Matrix->Complex) { for (Element = Matrix->FirstInRow[Row]; Element != NULL; Element = Element->NextInRow) @@ -568,9 +519,7 @@ int Row; Element->Real = 0.0; Element->Imag = 0.0; } - } else -#endif - { + } else { for (Element = Matrix->FirstInRow[Row]; Element != NULL; Element = Element->NextInRow) @@ -579,7 +528,7 @@ int Row; } } - return spError( (char *)Matrix ); + return spError( (void *)Matrix ); } #ifdef PARALLEL_ARCH @@ -587,24 +536,20 @@ int Row; * SMPcombine() */ void -SMPcombine( Matrix, RHS, Spare ) -SMPmatrix *Matrix; -double RHS[], Spare[]; +SMPcombine(SMPmatrix *Matrix, double RHS[], double Spare[]) { - spSetReal( (char *)Matrix ); - spCombine( (char *)Matrix, RHS, Spare, (spREAL*)NULL, (spREAL*)NULL ); + spSetReal( (void *)Matrix ); + spCombine( (void *)Matrix, RHS, Spare, (spREAL*)NULL, (spREAL*)NULL ); } /* * SMPcCombine() */ void -SMPcCombine( Matrix, RHS, Spare, iRHS, iSpare ) -SMPmatrix *Matrix; -double RHS[], Spare[]; -double iRHS[], iSpare[]; +SMPcCombine(SMPmatrix *Matrix, double RHS[], double Spare[], + double iRHS[], double iSpare[]) { - spSetComplex( (char *)Matrix ); - spCombine( (char *)Matrix, RHS, Spare, iRHS, iSpare ); + spSetComplex( (void *)Matrix ); + spCombine( (void *)Matrix, RHS, Spare, iRHS, iSpare ); } #endif /* PARALLEL_ARCH */ diff --git a/src/maths/sparse/spsolve.c b/src/maths/sparse/spsolve.c index 9fa83605f..7db3c85be 100644 --- a/src/maths/sparse/spsolve.c +++ b/src/maths/sparse/spsolve.c @@ -45,6 +45,7 @@ * spDefs.h * Matrix type and macro definitions for the sparse matrix routines. */ +#include #define spINSIDE_SPARSE #include "spconfig.h" @@ -58,20 +59,10 @@ * Function declarations */ -#ifdef __STDC__ -#if spSEPARATED_COMPLEX_VECTORS static void SolveComplexMatrix( MatrixPtr, RealVector, RealVector, RealVector, RealVector ); static void SolveComplexTransposedMatrix( MatrixPtr, RealVector, RealVector, RealVector, RealVector ); -#else -static void SolveComplexMatrix( MatrixPtr, RealVector, RealVector ); -static void SolveComplexTransposedMatrix( MatrixPtr, RealVector, RealVector ); -#endif -#else /* __STDC__ */ -static void SolveComplexMatrix(); -static void SolveComplexTransposedMatrix(); -#endif /* __STDC__ */ @@ -102,13 +93,10 @@ static void SolveComplexTransposedMatrix(); * iRHS (RealVector) * iRHS is the imaginary portion of the input data array, the right * hand side. This data is undisturbed and may be reused for other solves. - * This argument is only necessary if matrix is complex and if - * spSEPARATED_COMPLEX_VECTOR is set TRUE. * iSolution (RealVector) * iSolution is the imaginary portion of the output data array. This * routine is constructed such that iRHS and iSolution can be - * the same array. This argument is only necessary if matrix is complex - * and if spSEPARATED_COMPLEX_VECTOR is set TRUE. + * the same array. * * >>> Local variables: * Intermediate (RealVector) @@ -131,89 +119,77 @@ static void SolveComplexTransposedMatrix(); * Size of matrix. Made local to reduce indirection. * Temp (RealNumber) * Temporary storage for entries in arrays. - * - * >>> 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. */ /*VARARGS3*/ void -spSolve( eMatrix, RHS, Solution IMAG_VECTORS ) - -char *eMatrix; -RealVector RHS, Solution IMAG_VECTORS; +spSolve(void *eMatrix, RealVector RHS, RealVector Solution, + RealVector iRHS, RealVector iSolution) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register ElementPtr pElement; -register RealVector Intermediate; -register RealNumber Temp; -register int I, *pExtOrder, Size; -ElementPtr pPivot; -void SolveComplexMatrix(); + MatrixPtr Matrix = (MatrixPtr)eMatrix; + ElementPtr pElement; + RealVector Intermediate; + RealNumber Temp; + int I, *pExtOrder, Size; + ElementPtr pPivot; + void SolveComplexMatrix(); -/* Begin `spSolve'. */ - ASSERT( IS_VALID(Matrix) AND IS_FACTORED(Matrix) ); + /* Begin `spSolve'. */ + assert( IS_VALID(Matrix) && IS_FACTORED(Matrix) ); -#if spCOMPLEX if (Matrix->Complex) - { SolveComplexMatrix( Matrix, RHS, Solution IMAG_VECTORS ); + { + SolveComplexMatrix( Matrix, RHS, Solution, iRHS, iSolution ); return; } -#endif -#if REAL Intermediate = Matrix->Intermediate; Size = Matrix->Size; -/* Correct array pointers for ARRAY_OFFSET. */ -#if NOT ARRAY_OFFSET - --RHS; - --Solution; -#endif - -/* Initialize Intermediate vector. */ + /* Initialize Intermediate vector. */ pExtOrder = &Matrix->IntToExtRowMap[Size]; for (I = Size; I > 0; I--) Intermediate[I] = RHS[*(pExtOrder--)]; -/* Forward elimination. Solves Lc = b.*/ + /* Forward elimination. Solves Lc = b.*/ for (I = 1; I <= Size; I++) - { -/* This step of the elimination is skipped if Temp equals zero. */ + { + + /* This step of the elimination is skipped if Temp equals zero. */ if ((Temp = Intermediate[I]) != 0.0) - { pPivot = Matrix->Diag[I]; + { + pPivot = Matrix->Diag[I]; Intermediate[I] = (Temp *= pPivot->Real); pElement = pPivot->NextInCol; while (pElement != NULL) - { Intermediate[pElement->Row] -= Temp * pElement->Real; + { + Intermediate[pElement->Row] -= Temp * pElement->Real; pElement = pElement->NextInCol; } } } -/* Backward Substitution. Solves Ux = c.*/ + /* Backward Substitution. Solves Ux = c.*/ for (I = Size; I > 0; I--) - { Temp = Intermediate[I]; + { + Temp = Intermediate[I]; pElement = Matrix->Diag[I]->NextInRow; while (pElement != NULL) - { Temp -= pElement->Real * Intermediate[pElement->Col]; + { + Temp -= pElement->Real * Intermediate[pElement->Col]; pElement = pElement->NextInRow; } Intermediate[I] = Temp; } -/* Unscramble Intermediate vector while placing data in to Solution vector. */ + /* Unscramble Intermediate vector while placing data in to Solution vector. */ pExtOrder = &Matrix->IntToExtColMap[Size]; for (I = Size; I > 0; I--) Solution[*(pExtOrder--)] = Intermediate[I]; return; -#endif /* REAL */ } @@ -226,7 +202,6 @@ void SolveComplexMatrix(); -#if spCOMPLEX /* * SOLVE COMPLEX MATRIX EQUATION * @@ -280,69 +255,50 @@ void SolveComplexMatrix(); * Size of matrix. Made local to reduce indirection. * Temp (ComplexNumber) * Temporary storage for entries in arrays. - * - * >>> 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 -SolveComplexMatrix( Matrix, RHS, Solution IMAG_VECTORS ) +SolveComplexMatrix( Matrix, RHS, Solution , iRHS, iSolution ) MatrixPtr Matrix; -RealVector RHS, Solution IMAG_VECTORS; +RealVector RHS, Solution , iRHS, iSolution; { -register ElementPtr pElement; -register ComplexVector Intermediate; -register int I, *pExtOrder, Size; -ElementPtr pPivot; -ComplexNumber Temp; + ElementPtr pElement; + ComplexVector Intermediate; + int I, *pExtOrder, Size; + ElementPtr pPivot; + ComplexNumber Temp; -/* Begin `SolveComplexMatrix'. */ + /* Begin `SolveComplexMatrix'. */ Size = Matrix->Size; Intermediate = (ComplexVector)Matrix->Intermediate; -/* 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. */ + /* Initialize Intermediate vector. */ pExtOrder = &Matrix->IntToExtRowMap[Size]; -#if spSEPARATED_COMPLEX_VECTORS for (I = Size; I > 0; I--) - { Intermediate[I].Real = RHS[*(pExtOrder)]; + { + Intermediate[I].Real = RHS[*(pExtOrder)]; Intermediate[I].Imag = iRHS[*(pExtOrder--)]; } -#else - ExtVector = (ComplexVector)RHS; - for (I = Size; I > 0; I--) - Intermediate[I] = ExtVector[*(pExtOrder--)]; -#endif -/* Forward substitution. Solves Lc = b.*/ + /* Forward substitution. Solves Lc = b.*/ for (I = 1; I <= Size; I++) - { Temp = Intermediate[I]; + { + Temp = Intermediate[I]; -/* This step of the substitution is skipped if Temp equals zero. */ - if ((Temp.Real != 0.0) OR (Temp.Imag != 0.0)) - { pPivot = Matrix->Diag[I]; -/* Cmplx expr: Temp *= (1.0 / Pivot). */ + /* This step of the substitution is skipped if Temp equals zero. */ + if ((Temp.Real != 0.0) || (Temp.Imag != 0.0)) + { + pPivot = Matrix->Diag[I]; + /* Cmplx expr: Temp *= (1.0 / Pivot). */ CMPLX_MULT_ASSIGN(Temp, *pPivot); Intermediate[I] = Temp; pElement = pPivot->NextInCol; while (pElement != NULL) { -/* Cmplx expr: Intermediate[Element->Row] -= Temp * *Element. */ + /* Cmplx expr: Intermediate[Element->Row] -= Temp * *Element. */ CMPLX_MULT_SUBT_ASSIGN(Intermediate[pElement->Row], Temp, *pElement); pElement = pElement->NextInCol; @@ -350,37 +306,32 @@ ComplexNumber Temp; } } -/* Backward Substitution. Solves Ux = c.*/ + /* Backward Substitution. Solves Ux = c.*/ for (I = Size; I > 0; I--) - { Temp = Intermediate[I]; + { + Temp = Intermediate[I]; pElement = Matrix->Diag[I]->NextInRow; while (pElement != NULL) { -/* Cmplx expr: Temp -= *Element * Intermediate[Element->Col]. */ + /* Cmplx expr: Temp -= *Element * Intermediate[Element->Col]. */ CMPLX_MULT_SUBT_ASSIGN(Temp, *pElement,Intermediate[pElement->Col]); pElement = pElement->NextInRow; } Intermediate[I] = Temp; } -/* Unscramble Intermediate vector while placing data in to Solution vector. */ + /* Unscramble Intermediate vector while placing data in to Solution vector. */ pExtOrder = &Matrix->IntToExtColMap[Size]; -#if spSEPARATED_COMPLEX_VECTORS for (I = Size; I > 0; I--) - { Solution[*(pExtOrder)] = Intermediate[I].Real; + { + Solution[*(pExtOrder)] = Intermediate[I].Real; iSolution[*(pExtOrder--)] = Intermediate[I].Imag; } -#else - ExtVector = (ComplexVector)Solution; - for (I = Size; I > 0; I--) - ExtVector[*(pExtOrder--)] = Intermediate[I]; -#endif return; } -#endif /* spCOMPLEX */ @@ -448,88 +399,77 @@ ComplexNumber Temp; * Size of matrix. Made local to reduce indirection. * Temp (RealNumber) * Temporary storage for entries in arrays. - * - * >>> 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. */ /*VARARGS3*/ void -spSolveTransposed( eMatrix, RHS, Solution IMAG_VECTORS ) - -char *eMatrix; -RealVector RHS, Solution IMAG_VECTORS; +spSolveTransposed(void *eMatrix, RealVector RHS, RealVector Solution, + RealVector iRHS, RealVector iSolution) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register ElementPtr pElement; -register RealVector Intermediate; -register int I, *pExtOrder, Size; -ElementPtr pPivot; -RealNumber Temp; -void SolveComplexTransposedMatrix(); + MatrixPtr Matrix = (MatrixPtr)eMatrix; + ElementPtr pElement; + RealVector Intermediate; + int I, *pExtOrder, Size; + ElementPtr pPivot; + RealNumber Temp; + void SolveComplexTransposedMatrix(); -/* Begin `spSolveTransposed'. */ - ASSERT( IS_VALID(Matrix) AND IS_FACTORED(Matrix) ); + /* Begin `spSolveTransposed'. */ + assert( IS_VALID(Matrix) && IS_FACTORED(Matrix) ); -#if spCOMPLEX if (Matrix->Complex) - { SolveComplexTransposedMatrix( Matrix, RHS, Solution IMAG_VECTORS ); + { + SolveComplexTransposedMatrix( Matrix, RHS, Solution , iRHS, iSolution ); return; } -#endif -#if REAL Size = Matrix->Size; Intermediate = Matrix->Intermediate; -/* Correct array pointers for ARRAY_OFFSET. */ -#if NOT ARRAY_OFFSET - --RHS; - --Solution; -#endif - -/* Initialize Intermediate vector. */ + /* Initialize Intermediate vector. */ pExtOrder = &Matrix->IntToExtColMap[Size]; for (I = Size; I > 0; I--) Intermediate[I] = RHS[*(pExtOrder--)]; -/* Forward elimination. */ + /* Forward elimination. */ for (I = 1; I <= Size; I++) - { -/* This step of the elimination is skipped if Temp equals zero. */ + { + + /* This step of the elimination is skipped if Temp equals zero. */ if ((Temp = Intermediate[I]) != 0.0) - { pElement = Matrix->Diag[I]->NextInRow; + { + pElement = Matrix->Diag[I]->NextInRow; while (pElement != NULL) - { Intermediate[pElement->Col] -= Temp * pElement->Real; + { + Intermediate[pElement->Col] -= Temp * pElement->Real; pElement = pElement->NextInRow; } } } -/* Backward Substitution. */ + /* Backward Substitution. */ for (I = Size; I > 0; I--) - { pPivot = Matrix->Diag[I]; + { + pPivot = Matrix->Diag[I]; Temp = Intermediate[I]; pElement = pPivot->NextInCol; while (pElement != NULL) - { Temp -= pElement->Real * Intermediate[pElement->Row]; + { + Temp -= pElement->Real * Intermediate[pElement->Row]; pElement = pElement->NextInCol; } Intermediate[I] = Temp * pPivot->Real; } -/* Unscramble Intermediate vector while placing data in to Solution vector. */ + /* Unscramble Intermediate vector while placing data in to + Solution vector. */ pExtOrder = &Matrix->IntToExtRowMap[Size]; for (I = Size; I > 0; I--) Solution[*(pExtOrder--)] = Intermediate[I]; return; -#endif /* REAL */ } #endif /* TRANSPOSE */ @@ -542,7 +482,7 @@ void SolveComplexTransposedMatrix(); -#if TRANSPOSE AND spCOMPLEX +#if TRANSPOSE /* * SOLVE COMPLEX TRANSPOSED MATRIX EQUATION * @@ -560,23 +500,18 @@ void SolveComplexTransposedMatrix(); * RHS (RealVector) * RHS is the input data array, the right hand * side. This data is undisturbed and may be reused for other solves. - * This vector is only the real portion if the matrix is complex and - * spSEPARATED_COMPLEX_VECTORS is set TRUE. + * This vector is only the real portion if the matrix is complex. * Solution (RealVector) * Solution is the real portion of the output data array. This routine * is constructed such that RHS and Solution can be the same array. - * This vector is only the real portion if the matrix is complex and - * spSEPARATED_COMPLEX_VECTORS is set TRUE. + * This vector is only the real portion if the matrix is complex. * iRHS (RealVector) * iRHS is the imaginary portion of the input data array, the right * hand side. This data is undisturbed and may be reused for other solves. - * If either spCOMPLEX or spSEPARATED_COMPLEX_VECTOR is set FALSE, there - * is no need to supply this array. * iSolution (RealVector) * iSolution is the imaginary portion of the output data array. This * routine is constructed such that iRHS and iSolution can be - * the same array. If spCOMPLEX or spSEPARATED_COMPLEX_VECTOR is set - * FALSE, there is no need to supply this array. + * the same array. * * >>> Local variables: * Intermediate (ComplexVector) @@ -599,65 +534,46 @@ void SolveComplexTransposedMatrix(); * Size of matrix. Made local to reduce indirection. * Temp (ComplexNumber) * Temporary storage for entries in arrays. - * - * >>> 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 -SolveComplexTransposedMatrix(Matrix, RHS, Solution IMAG_VECTORS ) +SolveComplexTransposedMatrix(Matrix, RHS, Solution , iRHS, iSolution ) MatrixPtr Matrix; -RealVector RHS, Solution IMAG_VECTORS; +RealVector RHS, Solution , iRHS, iSolution; { -register ElementPtr pElement; -register ComplexVector Intermediate; -register int I, *pExtOrder, Size; -ElementPtr pPivot; -ComplexNumber Temp; + ElementPtr pElement; + ComplexVector Intermediate; + int I, *pExtOrder, Size; + ElementPtr pPivot; + ComplexNumber Temp; -/* Begin `SolveComplexTransposedMatrix'. */ + /* Begin `SolveComplexTransposedMatrix'. */ Size = Matrix->Size; Intermediate = (ComplexVector)Matrix->Intermediate; -/* 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. */ + /* Initialize Intermediate vector. */ pExtOrder = &Matrix->IntToExtColMap[Size]; -#if spSEPARATED_COMPLEX_VECTORS for (I = Size; I > 0; I--) - { Intermediate[I].Real = RHS[*(pExtOrder)]; + { + Intermediate[I].Real = RHS[*(pExtOrder)]; Intermediate[I].Imag = iRHS[*(pExtOrder--)]; } -#else - ExtVector = (ComplexVector)RHS; - for (I = Size; I > 0; I--) - Intermediate[I] = ExtVector[*(pExtOrder--)]; -#endif -/* Forward elimination. */ + /* Forward elimination. */ for (I = 1; I <= Size; I++) - { Temp = Intermediate[I]; + { + Temp = Intermediate[I]; -/* This step of the elimination is skipped if Temp equals zero. */ - if ((Temp.Real != 0.0) OR (Temp.Imag != 0.0)) - { pElement = Matrix->Diag[I]->NextInRow; + /* This step of the elimination is skipped if Temp equals zero. */ + if ((Temp.Real != 0.0) || (Temp.Imag != 0.0)) + { + pElement = Matrix->Diag[I]->NextInRow; while (pElement != NULL) { -/* Cmplx expr: Intermediate[Element->Col] -= Temp * *Element. */ + /* Cmplx expr: Intermediate[Element->Col] -= Temp * *Element. */ CMPLX_MULT_SUBT_ASSIGN( Intermediate[pElement->Col], Temp, *pElement); pElement = pElement->NextInRow; @@ -665,37 +581,34 @@ ComplexNumber Temp; } } -/* Backward Substitution. */ + /* Backward Substitution. */ for (I = Size; I > 0; I--) - { pPivot = Matrix->Diag[I]; + { + pPivot = Matrix->Diag[I]; Temp = Intermediate[I]; pElement = pPivot->NextInCol; while (pElement != NULL) { -/* Cmplx expr: Temp -= Intermediate[Element->Row] * *Element. */ + /* Cmplx expr: Temp -= Intermediate[Element->Row] * *Element. */ CMPLX_MULT_SUBT_ASSIGN(Temp,Intermediate[pElement->Row],*pElement); pElement = pElement->NextInCol; } -/* Cmplx expr: Intermediate = Temp * (1.0 / *pPivot). */ + /* Cmplx expr: Intermediate = Temp * (1.0 / *pPivot). */ CMPLX_MULT(Intermediate[I], Temp, *pPivot); } -/* Unscramble Intermediate vector while placing data in to Solution vector. */ + /* Unscramble Intermediate vector while placing data in to + Solution vector. */ pExtOrder = &Matrix->IntToExtRowMap[Size]; -#if spSEPARATED_COMPLEX_VECTORS for (I = Size; I > 0; I--) - { Solution[*(pExtOrder)] = Intermediate[I].Real; + { + Solution[*(pExtOrder)] = Intermediate[I].Real; iSolution[*(pExtOrder--)] = Intermediate[I].Imag; } -#else - ExtVector = (ComplexVector)Solution; - for (I = Size; I > 0; I--) - ExtVector[*(pExtOrder--)] = Intermediate[I]; -#endif return; } -#endif /* TRANSPOSE AND spCOMPLEX */ +#endif /* TRANSPOSE */ diff --git a/src/maths/sparse/sputils.c b/src/maths/sparse/sputils.c index 1821d9ebf..7aa20c329 100644 --- a/src/maths/sparse/sputils.c +++ b/src/maths/sparse/sputils.c @@ -58,6 +58,7 @@ * spDefs.h * Matrix type and macro definitions for the sparse matrix routines. */ +#include #define spINSIDE_SPARSE #include "spconfig.h" @@ -68,21 +69,14 @@ -/* - * Function declarations - */ +/* Function declarations */ static int CountTwins( MatrixPtr, int, ElementPtr*, ElementPtr* ); static void SwapCols( MatrixPtr, ElementPtr, ElementPtr ); -#if spSEPARATED_COMPLEX_VECTORS static void ComplexMatrixMultiply( MatrixPtr, RealVector, RealVector, RealVector, RealVector ); static void ComplexTransposedMatrixMultiply( MatrixPtr, RealVector, RealVector, RealVector, RealVector); -#else -static void ComplexMatrixMultiply( MatrixPtr, RealVector, RealVector ); -static void ComplexTransposedMatrixMultiply(MatrixPtr, RealVector, RealVector); -#endif @@ -167,59 +161,65 @@ static void ComplexTransposedMatrixMultiply(MatrixPtr, RealVector, RealVector); * pTwin2 (ElementPtr) * Pointer to the twin found in the row belonging to the zero diagonal. * belonging to the zero diagonal. - * AnotherPassNeeded (BOOLEAN) + * AnotherPassNeeded (int) * Flag indicating that at least one zero diagonal with symmetric twins * remain. * StartAt (int) * Column number of first zero diagonal with symmetric twins. - * Swapped (BOOLEAN) + * Swapped (int) * Flag indicating that columns were swapped on this pass. * Twins (int) * Number of symmetric twins corresponding to current zero diagonal. */ void -spMNA_Preorder( eMatrix ) - -char *eMatrix; +spMNA_Preorder(void *eMatrix) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register int J, Size; -ElementPtr pTwin1, pTwin2; -int Twins, StartAt = 1; -BOOLEAN Swapped, AnotherPassNeeded; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + int J, Size; + ElementPtr pTwin1, pTwin2; + int Twins, StartAt = 1; + int Swapped, AnotherPassNeeded; -/* Begin `spMNA_Preorder'. */ - ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored ); + /* Begin `spMNA_Preorder'. */ + assert( IS_VALID(Matrix) && !Matrix->Factored ); if (Matrix->RowsLinked) return; Size = Matrix->Size; Matrix->Reordered = YES; do - { AnotherPassNeeded = Swapped = NO; + { + AnotherPassNeeded = Swapped = NO; -/* Search for zero diagonals with lone twins. */ + /* Search for zero diagonals with lone twins. */ for (J = StartAt; J <= Size; J++) - { if (Matrix->Diag[J] == NULL) - { Twins = CountTwins( Matrix, J, &pTwin1, &pTwin2 ); + { + if (Matrix->Diag[J] == NULL) + { + Twins = CountTwins( Matrix, J, &pTwin1, &pTwin2 ); if (Twins == 1) - { /* Lone twins found, swap rows. */ + { + /* Lone twins found, swap rows. */ SwapCols( Matrix, pTwin1, pTwin2 ); Swapped = YES; } - else if ((Twins > 1) AND NOT AnotherPassNeeded) - { AnotherPassNeeded = YES; + else if ((Twins > 1) && !AnotherPassNeeded) + { + AnotherPassNeeded = YES; StartAt = J; } } } -/* All lone twins are gone, look for zero diagonals with multiple twins. */ + /* All lone twins are gone, look for zero diagonals with multiple twins. */ if (AnotherPassNeeded) - { for (J = StartAt; NOT Swapped AND (J <= Size); J++) - { if (Matrix->Diag[J] == NULL) - { Twins = CountTwins( Matrix, J, &pTwin1, &pTwin2 ); + { + for (J = StartAt; !Swapped && (J <= Size); J++) + { + if (Matrix->Diag[J] == NULL) + { + Twins = CountTwins( Matrix, J, &pTwin1, &pTwin2 ); SwapCols( Matrix, pTwin1, pTwin2 ); Swapped = YES; } @@ -247,20 +247,23 @@ MatrixPtr Matrix; int Col; ElementPtr *ppTwin1, *ppTwin2; { -int Row, Twins = 0; -ElementPtr pTwin1, pTwin2; + int Row, Twins = 0; + ElementPtr pTwin1, pTwin2; -/* Begin `CountTwins'. */ + /* Begin `CountTwins'. */ pTwin1 = Matrix->FirstInCol[Col]; while (pTwin1 != NULL) - { if (ABS(pTwin1->Real) == 1.0) - { Row = pTwin1->Row; + { + if (ABS(pTwin1->Real) == 1.0) + { + Row = pTwin1->Row; pTwin2 = Matrix->FirstInCol[Row]; - while ((pTwin2 != NULL) AND (pTwin2->Row != Col)) + while ((pTwin2 != NULL) && (pTwin2->Row != Col)) pTwin2 = pTwin2->NextInCol; - if ((pTwin2 != NULL) AND (ABS(pTwin2->Real) == 1.0)) - { /* Found symmetric twins. */ + if ((pTwin2 != NULL) && (ABS(pTwin2->Real) == 1.0)) + { + /* Found symmetric twins. */ if (++Twins >= 2) return Twins; (*ppTwin1 = pTwin1)->Col = Col; (*ppTwin2 = pTwin2)->Col = Row; @@ -287,9 +290,9 @@ SwapCols( Matrix, pTwin1, pTwin2 ) MatrixPtr Matrix; ElementPtr pTwin1, pTwin2; { -int Col1 = pTwin1->Col, Col2 = pTwin2->Col; + int Col1 = pTwin1->Col, Col2 = pTwin2->Col; -/* Begin `SwapCols'. */ + /* Begin `SwapCols'. */ SWAP (ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]); SWAP (int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]); @@ -300,7 +303,7 @@ int Col1 = pTwin1->Col, Col2 = pTwin2->Col; Matrix->Diag[Col1] = pTwin2; Matrix->Diag[Col2] = pTwin1; - Matrix->NumberOfInterchangesIsOdd = NOT Matrix->NumberOfInterchangesIsOdd; + Matrix->NumberOfInterchangesIsOdd = !Matrix->NumberOfInterchangesIsOdd; return; } #endif /* MODIFIED_NODAL */ @@ -374,62 +377,59 @@ void spScale( eMatrix, RHS_ScaleFactors, SolutionScaleFactors ) char *eMatrix; -register RealVector RHS_ScaleFactors, SolutionScaleFactors; + RealVector RHS_ScaleFactors, SolutionScaleFactors; { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register ElementPtr pElement; -register int I, lSize, *pExtOrder; -RealNumber ScaleFactor; -void ScaleComplexMatrix(); + MatrixPtr Matrix = (MatrixPtr)eMatrix; + ElementPtr pElement; + int I, lSize, *pExtOrder; + RealNumber ScaleFactor; + void ScaleComplexMatrix(); -/* Begin `spScale'. */ - ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored ); - if (NOT Matrix->RowsLinked) spcLinkRows( Matrix ); + /* Begin `spScale'. */ + assert( IS_VALID(Matrix) && !Matrix->Factored ); + if (!Matrix->RowsLinked) spcLinkRows( Matrix ); -#if spCOMPLEX if (Matrix->Complex) - { ScaleComplexMatrix( Matrix, RHS_ScaleFactors, SolutionScaleFactors ); + { + ScaleComplexMatrix( Matrix, RHS_ScaleFactors, SolutionScaleFactors ); return; } -#endif -#if REAL lSize = Matrix->Size; -/* Correct pointers to arrays for ARRAY_OFFSET */ -#if NOT ARRAY_OFFSET - --RHS_ScaleFactors; - --SolutionScaleFactors; -#endif - -/* Scale Rows */ + /* Scale Rows */ pExtOrder = &Matrix->IntToExtRowMap[1]; for (I = 1; I <= lSize; I++) - { if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0) - { pElement = Matrix->FirstInRow[I]; + { + if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0) + { + pElement = Matrix->FirstInRow[I]; while (pElement != NULL) - { pElement->Real *= ScaleFactor; + { + pElement->Real *= ScaleFactor; pElement = pElement->NextInRow; } } } -/* Scale Columns */ + /* Scale Columns */ pExtOrder = &Matrix->IntToExtColMap[1]; for (I = 1; I <= lSize; I++) - { if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0) - { pElement = Matrix->FirstInCol[I]; + { + if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0) + { + pElement = Matrix->FirstInCol[I]; while (pElement != NULL) - { pElement->Real *= ScaleFactor; + { + pElement->Real *= ScaleFactor; pElement = pElement->NextInCol; } } } return; -#endif /* REAL */ } #endif /* SCALING */ @@ -441,7 +441,7 @@ void ScaleComplexMatrix(); -#if spCOMPLEX AND SCALING +#if SCALING /* * SCALE COMPLEX MATRIX * @@ -502,43 +502,43 @@ static void ScaleComplexMatrix( Matrix, RHS_ScaleFactors, SolutionScaleFactors ) MatrixPtr Matrix; -register RealVector RHS_ScaleFactors, SolutionScaleFactors; + RealVector RHS_ScaleFactors, SolutionScaleFactors; { -register ElementPtr pElement; -register int I, lSize, *pExtOrder; -RealNumber ScaleFactor; + ElementPtr pElement; + int I, lSize, *pExtOrder; + RealNumber ScaleFactor; -/* Begin `ScaleComplexMatrix'. */ + /* Begin `ScaleComplexMatrix'. */ lSize = Matrix->Size; -/* Correct pointers to arrays for ARRAY_OFFSET */ -#if NOT ARRAY_OFFSET - --RHS_ScaleFactors; - --SolutionScaleFactors; -#endif - -/* Scale Rows */ + /* Scale Rows */ pExtOrder = &Matrix->IntToExtRowMap[1]; for (I = 1; I <= lSize; I++) - { if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0) - { pElement = Matrix->FirstInRow[I]; + { + if ((ScaleFactor = RHS_ScaleFactors[*(pExtOrder++)]) != 1.0) + { + pElement = Matrix->FirstInRow[I]; while (pElement != NULL) - { pElement->Real *= ScaleFactor; + { + pElement->Real *= ScaleFactor; pElement->Imag *= ScaleFactor; pElement = pElement->NextInRow; } } } -/* Scale Columns */ + /* Scale Columns */ pExtOrder = &Matrix->IntToExtColMap[1]; for (I = 1; I <= lSize; I++) - { if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0) - { pElement = Matrix->FirstInCol[I]; + { + if ((ScaleFactor = SolutionScaleFactors[*(pExtOrder++)]) != 1.0) + { + pElement = Matrix->FirstInCol[I]; while (pElement != NULL) - { pElement->Real *= ScaleFactor; + { + pElement->Real *= ScaleFactor; pElement->Imag *= ScaleFactor; pElement = pElement->NextInCol; } @@ -546,7 +546,7 @@ RealNumber ScaleFactor; } return; } -#endif /* SCALING AND spCOMPLEX */ +#endif /* SCALING */ @@ -573,55 +573,37 @@ RealNumber ScaleFactor; * Solution is the vector being multiplied by the matrix. * iRHS (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. + * what is being solved for. * iSolution (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. + * by the matrix. */ void -spMultiply( eMatrix, RHS, Solution IMAG_VECTORS ) - -char *eMatrix; -RealVector RHS, Solution IMAG_VECTORS; +spMultiply(void *eMatrix, RealVector RHS, RealVector Solution, + RealVector iRHS, RealVector iSolution) { -register ElementPtr pElement; -register RealVector Vector; -register RealNumber Sum; -register int I, *pExtOrder; -MatrixPtr Matrix = (MatrixPtr)eMatrix; -extern void ComplexMatrixMultiply(); + ElementPtr pElement; + RealVector Vector; + RealNumber Sum; + int I, *pExtOrder; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + extern void ComplexMatrixMultiply(); -/* Begin `spMultiply'. */ - ASSERT( IS_SPARSE( Matrix ) AND NOT Matrix->Factored ); - if (NOT Matrix->RowsLinked) + /* Begin `spMultiply'. */ + assert( IS_SPARSE( Matrix ) && !Matrix->Factored ); + if (!Matrix->RowsLinked) spcLinkRows(Matrix); - if (NOT Matrix->InternalVectorsAllocated) + if (!Matrix->InternalVectorsAllocated) spcCreateInternalVectors( Matrix ); -#if spCOMPLEX if (Matrix->Complex) - { ComplexMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS ); + { + ComplexMatrixMultiply( Matrix, RHS, Solution , iRHS, iSolution ); return; } -#endif -#if REAL -#if NOT ARRAY_OFFSET -/* Correct array pointers for ARRAY_OFFSET. */ - --RHS; - --Solution; -#endif - -/* Initialize Intermediate vector with reordered Solution vector. */ + /* Initialize Intermediate vector with reordered Solution vector. */ Vector = Matrix->Intermediate; pExtOrder = &Matrix->IntToExtColMap[Matrix->Size]; for (I = Matrix->Size; I > 0; I--) @@ -629,17 +611,18 @@ extern void ComplexMatrixMultiply(); pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size]; for (I = Matrix->Size; I > 0; I--) - { pElement = Matrix->FirstInRow[I]; + { + pElement = Matrix->FirstInRow[I]; Sum = 0.0; while (pElement != NULL) - { Sum += pElement->Real * Vector[pElement->Col]; + { + Sum += pElement->Real * Vector[pElement->Col]; pElement = pElement->NextInRow; } RHS[*pExtOrder--] = Sum; } return; -#endif /* REAL */ } #endif /* MULTIPLICATION */ @@ -649,7 +632,7 @@ extern void ComplexMatrixMultiply(); -#if spCOMPLEX AND MULTIPLICATION +#if MULTIPLICATION /* * COMPLEX MATRIX MULTIPLICATION * @@ -662,86 +645,58 @@ extern void ComplexMatrixMultiply(); * Pointer to the matrix. * RHS (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 (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. + * Solution is the vector being multiplied by the matrix. * iRHS (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. + * what is being solved for. * iSolution (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. + * by the matrix. */ static void -ComplexMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS ) +ComplexMatrixMultiply( Matrix, RHS, Solution , iRHS, iSolution ) MatrixPtr Matrix; -RealVector RHS, Solution IMAG_VECTORS; +RealVector RHS, Solution , iRHS, iSolution; { -register ElementPtr pElement; -register ComplexVector Vector; -ComplexNumber Sum; -register int I, *pExtOrder; + ElementPtr pElement; + ComplexVector Vector; + ComplexNumber Sum; + int I, *pExtOrder; -/* Begin `ComplexMatrixMultiply'. */ + /* Begin `ComplexMatrixMultiply'. */ -/* 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. */ + /* Initialize Intermediate vector with reordered Solution vector. */ Vector = (ComplexVector)Matrix->Intermediate; pExtOrder = &Matrix->IntToExtColMap[Matrix->Size]; -#if spSEPARATED_COMPLEX_VECTORS for (I = Matrix->Size; I > 0; I--) - { Vector[I].Real = Solution[*pExtOrder]; + { + 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->IntToExtRowMap[Matrix->Size]; for (I = Matrix->Size; I > 0; I--) - { pElement = Matrix->FirstInRow[I]; + { + pElement = Matrix->FirstInRow[I]; Sum.Real = Sum.Imag = 0.0; while (pElement != NULL) - { /* Cmplx expression : Sum += Element * Vector[Col] */ + { + /* Cmplx expression : Sum += Element * Vector[Col] */ CMPLX_MULT_ADD_ASSIGN( Sum, *pElement, Vector[pElement->Col] ); pElement = pElement->NextInRow; } -#if spSEPARATED_COMPLEX_VECTORS RHS[*pExtOrder] = Sum.Real; iRHS[*pExtOrder--] = Sum.Imag; -#else - ((ComplexVector)RHS)[*pExtOrder--] = Sum; -#endif } return; } -#endif /* spCOMPLEX AND MULTIPLICATION */ +#endif /* MULTIPLICATION */ @@ -750,7 +705,7 @@ register int I, *pExtOrder; -#if MULTIPLICATION AND TRANSPOSE +#if MULTIPLICATION && TRANSPOSE /* * TRANSPOSED MATRIX MULTIPLICATION * @@ -768,53 +723,35 @@ register int I, *pExtOrder; * Solution is the vector being multiplied by the matrix. * iRHS (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. + * what is being solved for. * iSolution (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. + * by the matrix. */ void -spMultTransposed( eMatrix, RHS, Solution IMAG_VECTORS ) - -char *eMatrix; -RealVector RHS, Solution IMAG_VECTORS; +spMultTransposed(void *eMatrix, RealVector RHS, RealVector Solution, + RealVector iRHS, RealVector iSolution) { -register ElementPtr pElement; -register RealVector Vector; -register RealNumber Sum; -register int I, *pExtOrder; -MatrixPtr Matrix = (MatrixPtr)eMatrix; -extern void ComplexTransposedMatrixMultiply(); + ElementPtr pElement; + RealVector Vector; + RealNumber Sum; + int I, *pExtOrder; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + extern void ComplexTransposedMatrixMultiply(); -/* Begin `spMultTransposed'. */ - ASSERT( IS_SPARSE( Matrix ) AND NOT Matrix->Factored ); - if (NOT Matrix->InternalVectorsAllocated) + /* Begin `spMultTransposed'. */ + assert( IS_SPARSE( Matrix ) && !Matrix->Factored ); + if (!Matrix->InternalVectorsAllocated) spcCreateInternalVectors( Matrix ); -#if spCOMPLEX if (Matrix->Complex) - { ComplexTransposedMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS ); + { + ComplexTransposedMatrixMultiply( Matrix, RHS, Solution , iRHS, iSolution ); return; } -#endif -#if REAL -#if NOT ARRAY_OFFSET -/* Correct array pointers for ARRAY_OFFSET. */ - --RHS; - --Solution; -#endif - -/* Initialize Intermediate vector with reordered Solution vector. */ + /* Initialize Intermediate vector with reordered Solution vector. */ Vector = Matrix->Intermediate; pExtOrder = &Matrix->IntToExtRowMap[Matrix->Size]; for (I = Matrix->Size; I > 0; I--) @@ -822,19 +759,20 @@ extern void ComplexTransposedMatrixMultiply(); pExtOrder = &Matrix->IntToExtColMap[Matrix->Size]; for (I = Matrix->Size; I > 0; I--) - { pElement = Matrix->FirstInCol[I]; + { + pElement = Matrix->FirstInCol[I]; Sum = 0.0; while (pElement != NULL) - { Sum += pElement->Real * Vector[pElement->Row]; + { + Sum += pElement->Real * Vector[pElement->Row]; pElement = pElement->NextInCol; } RHS[*pExtOrder--] = Sum; } return; -#endif /* REAL */ } -#endif /* MULTIPLICATION AND TRANSPOSE */ +#endif /* MULTIPLICATION && TRANSPOSE */ @@ -842,7 +780,7 @@ extern void ComplexTransposedMatrixMultiply(); -#if spCOMPLEX AND MULTIPLICATION AND TRANSPOSE +#if MULTIPLICATION && TRANSPOSE /* * COMPLEX TRANSPOSED MATRIX MULTIPLICATION * @@ -855,86 +793,58 @@ extern void ComplexTransposedMatrixMultiply(); * Pointer to the matrix. * RHS (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 (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. + * Solution is the vector being multiplied by the matrix. * iRHS (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. + * what is being solved for. * iSolution (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. + * by the matrix. */ static void -ComplexTransposedMatrixMultiply( Matrix, RHS, Solution IMAG_VECTORS ) +ComplexTransposedMatrixMultiply( Matrix, RHS, Solution , iRHS, iSolution ) MatrixPtr Matrix; -RealVector RHS, Solution IMAG_VECTORS; +RealVector RHS, Solution , iRHS, iSolution; { -register ElementPtr pElement; -register ComplexVector Vector; -ComplexNumber Sum; -register int I, *pExtOrder; + ElementPtr pElement; + ComplexVector Vector; + ComplexNumber Sum; + int I, *pExtOrder; -/* Begin `ComplexMatrixMultiply'. */ + /* Begin `ComplexMatrixMultiply'. */ -/* 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. */ + /* 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].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]; + { + pElement = Matrix->FirstInCol[I]; Sum.Real = Sum.Imag = 0.0; while (pElement != NULL) - { /* Cmplx expression : Sum += Element * Vector[Row] */ + { + /* 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 */ +#endif /* MULTIPLICATION && TRANSPOSE */ @@ -979,66 +889,60 @@ register int I, *pExtOrder; * Norm (RealNumber) * L-infinity norm of a complex number. * Size (int) - * Local storage for Matrix->Size. Placed in a register for speed. + * Local storage for Matrix->Size. Placed in a for speed. * Temp (RealNumber) * Temporary storage for real portion of determinant. */ -#if spCOMPLEX void -spDeterminant( eMatrix, pExponent, pDeterminant, piDeterminant ) -RealNumber *piDeterminant; -#else -void -spDeterminant( eMatrix, pExponent, pDeterminant ) -#endif - -char *eMatrix; -register RealNumber *pDeterminant; -int *pExponent; +spDeterminant(void *eMatrix, int *pExponent, RealNumber *pDeterminant, + RealNumber *piDeterminant) { -register MatrixPtr Matrix = (MatrixPtr)eMatrix; -register int I, Size; -RealNumber Norm, nr, ni; -ComplexNumber Pivot, cDeterminant; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + int I, Size; + RealNumber Norm, nr, ni; + ComplexNumber Pivot, cDeterminant; #define NORM(a) (nr = ABS((a).Real), ni = ABS((a).Imag), MAX (nr,ni)) -/* Begin `spDeterminant'. */ - ASSERT( IS_SPARSE( Matrix ) AND IS_FACTORED(Matrix) ); + /* Begin `spDeterminant'. */ + assert( IS_SPARSE( Matrix ) && IS_FACTORED(Matrix) ); *pExponent = 0; if (Matrix->Error == spSINGULAR) - { *pDeterminant = 0.0; -#if spCOMPLEX + { + *pDeterminant = 0.0; 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.Real = 1.0; cDeterminant.Imag = 0.0; while (++I <= Size) - { CMPLX_RECIPROCAL( Pivot, *Matrix->Diag[I] ); + { + CMPLX_RECIPROCAL( Pivot, *Matrix->Diag[I] ); CMPLX_MULT_ASSIGN( cDeterminant, Pivot ); -/* Scale Determinant. */ + /* Scale Determinant. */ Norm = NORM( cDeterminant ); if (Norm != 0.0) - { while (Norm >= 1.0e12) - { cDeterminant.Real *= 1.0e-12; + { + 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.Real *= 1.0e12; cDeterminant.Imag *= 1.0e12; *pExponent -= 12; Norm = NORM( cDeterminant ); @@ -1046,17 +950,20 @@ ComplexNumber Pivot, cDeterminant; } } -/* Scale Determinant again, this time to be between 1.0 <= x < 10.0. */ + /* 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; + { + 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.Real *= 10.0; cDeterminant.Imag *= 10.0; (*pExponent)--; Norm = NORM( cDeterminant ); @@ -1068,45 +975,49 @@ ComplexNumber Pivot, cDeterminant; *pDeterminant = cDeterminant.Real; *piDeterminant = cDeterminant.Imag; } -#endif /* spCOMPLEX */ -#if REAL AND spCOMPLEX else -#endif -#if REAL - { /* Real Case. */ + { + /* Real Case. */ *pDeterminant = 1.0; while (++I <= Size) - { *pDeterminant /= Matrix->Diag[I]->Real; + { + *pDeterminant /= Matrix->Diag[I]->Real; -/* Scale Determinant. */ + /* Scale Determinant. */ if (*pDeterminant != 0.0) - { while (ABS(*pDeterminant) >= 1.0e12) - { *pDeterminant *= 1.0e-12; + { + while (ABS(*pDeterminant) >= 1.0e12) + { + *pDeterminant *= 1.0e-12; *pExponent += 12; } while (ABS(*pDeterminant) < 1.0e-12) - { *pDeterminant *= 1.0e12; + { + *pDeterminant *= 1.0e12; *pExponent -= 12; } } } -/* Scale Determinant again, this time to be between 1.0 <= x < 10.0. */ + /* 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; + { + while (ABS(*pDeterminant) >= 10.0) + { + *pDeterminant *= 0.1; (*pExponent)++; } while (ABS(*pDeterminant) < 1.0) - { *pDeterminant *= 10.0; + { + *pDeterminant *= 10.0; (*pExponent)--; } } if (Matrix->NumberOfInterchangesIsOdd) *pDeterminant = -*pDeterminant; } -#endif /* REAL */ } #endif /* DETERMINANT */ @@ -1147,25 +1058,27 @@ spStripFills( eMatrix ) char *eMatrix; { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -struct FillinListNodeStruct *pListNode; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + struct FillinListNodeStruct *pListNode; -/* Begin `spStripFills'. */ - ASSERT( IS_SPARSE( Matrix ) ); + /* 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; + /* Mark the fill-ins. */ + { + ElementPtr pFillin, pLastFillin; pListNode = Matrix->LastFillinListNode = Matrix->FirstFillinListNode; Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList; Matrix->NextAvailFillin = pListNode->pFillinList; while (pListNode != NULL) - { pFillin = pListNode->pFillinList; + { + pFillin = pListNode->pFillinList; pLastFillin = &(pFillin[ pListNode->NumberOfFillinsInList - 1 ]); while (pFillin <= pLastFillin) (pFillin++)->Row = 0; @@ -1173,16 +1086,20 @@ struct FillinListNodeStruct *pListNode; } } -/* Unlink fill-ins by searching for elements marked with Row = 0. */ - { register ElementPtr pElement, *ppElement; - register int I, Size = Matrix->Size; + /* Unlink fill-ins by searching for elements marked with Row = 0. */ + { + ElementPtr pElement, *ppElement; + int I, Size = Matrix->Size; -/* Unlink fill-ins in all columns. */ + /* Unlink fill-ins in all columns. */ for (I = 1; I <= Size; I++) - { ppElement = &(Matrix->FirstInCol[I]); + { + ppElement = &(Matrix->FirstInCol[I]); while ((pElement = *ppElement) != NULL) - { if (pElement->Row == 0) - { *ppElement = pElement->NextInCol; /* Unlink fill-in. */ + { + if (pElement->Row == 0) + { + *ppElement = pElement->NextInCol; /* Unlink fill-in. */ if (Matrix->Diag[pElement->Col] == pElement) Matrix->Diag[pElement->Col] = NULL; } @@ -1191,11 +1108,13 @@ struct FillinListNodeStruct *pListNode; } } -/* Unlink fill-ins in all rows. */ + /* Unlink fill-ins in all rows. */ for (I = 1; I <= Size; I++) - { ppElement = &(Matrix->FirstInRow[I]); + { + ppElement = &(Matrix->FirstInRow[I]); while ((pElement = *ppElement) != NULL) - { if (pElement->Row == 0) + { + if (pElement->Row == 0) *ppElement = pElement->NextInRow; /* Unlink fill-in. */ else ppElement = &pElement->NextInRow; /* Skip element. */ @@ -1205,18 +1124,18 @@ struct FillinListNodeStruct *pListNode; return; } -/* Same as above, but strips entire matrix without destroying the frame. - * This assumes that the matrix will be replaced with one of the same size. - */ +/* Same as above, but strips entire matrix without destroying the + * frame. This assumes that the matrix will be replaced with one of + * the same size. */ void spStripMatrix( eMatrix ) char *eMatrix; { -MatrixPtr Matrix = (MatrixPtr)eMatrix; + MatrixPtr Matrix = (MatrixPtr)eMatrix; -/* Begin `spStripMatrix'. */ - ASSERT( IS_SPARSE( Matrix ) ); + /* Begin `spStripMatrix'. */ + assert( IS_SPARSE( Matrix ) ); if (Matrix->Elements == 0) return; Matrix->RowsLinked = NO; Matrix->NeedsOrdering = YES; @@ -1224,8 +1143,9 @@ MatrixPtr Matrix = (MatrixPtr)eMatrix; Matrix->Originals = 0; Matrix->Fillins = 0; -/* Reset the element lists. */ - { register ElementPtr pElement; + /* Reset the element lists. */ + { + ElementPtr pElement; struct ElementListNodeStruct *pListNode; pListNode = Matrix->LastElementListNode = Matrix->FirstElementListNode; @@ -1233,8 +1153,9 @@ MatrixPtr Matrix = (MatrixPtr)eMatrix; Matrix->NextAvailElement = pListNode->pElementList; } -/* Reset the fill-in lists. */ - { register ElementPtr pFillin; + /* Reset the fill-in lists. */ + { + ElementPtr pFillin; struct FillinListNodeStruct *pListNode; pListNode = Matrix->LastFillinListNode = Matrix->FirstFillinListNode; @@ -1242,10 +1163,12 @@ MatrixPtr Matrix = (MatrixPtr)eMatrix; Matrix->NextAvailFillin = pListNode->pFillinList; } -/* Reset the Row, Column and Diag pointers */ - { register int I, Size = Matrix->Size; + /* Reset the Row, Column and Diag pointers */ + { + int I, Size = Matrix->Size; for (I = 1; I <= Size; I++) - { Matrix->FirstInRow[I] = NULL; + { + Matrix->FirstInRow[I] = NULL; Matrix->FirstInCol[I] = NULL; Matrix->Diag[I] = NULL; } @@ -1259,7 +1182,7 @@ MatrixPtr Matrix = (MatrixPtr)eMatrix; -#if TRANSLATE AND DELETE +#if TRANSLATE && DELETE /* * DELETE A ROW AND COLUMN FROM THE MATRIX * @@ -1295,55 +1218,53 @@ MatrixPtr Matrix = (MatrixPtr)eMatrix; */ void -spDeleteRowAndCol( eMatrix, Row, Col ) - -char *eMatrix; -int Row, Col; +spDeleteRowAndCol(char *eMatrix, int Row, int Col ) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register ElementPtr pElement, *ppElement, pLastElement; -int Size, ExtRow, ExtCol; -ElementPtr spcFindElementInCol(); + MatrixPtr Matrix = (MatrixPtr)eMatrix; + ElementPtr pElement, *ppElement, pLastElement; + int Size, ExtRow, ExtCol; + ElementPtr spcFindElementInCol(); -/* Begin `spDeleteRowAndCol'. */ + /* Begin `spDeleteRowAndCol'. */ - ASSERT( IS_SPARSE(Matrix) AND Row > 0 AND Col > 0 ); - ASSERT( Row <= Matrix->ExtSize AND Col <= Matrix->ExtSize ); + assert( IS_SPARSE(Matrix) && Row > 0 && Col > 0 ); + assert( Row <= Matrix->ExtSize && Col <= Matrix->ExtSize ); Size = Matrix->Size; ExtRow = Row; ExtCol = Col; - if (NOT Matrix->RowsLinked) spcLinkRows( Matrix ); + if (!Matrix->RowsLinked) + spcLinkRows( Matrix ); Row = Matrix->ExtToIntRowMap[Row]; Col = Matrix->ExtToIntColMap[Col]; - ASSERT( Row > 0 AND Col > 0 ); + assert( Row > 0 && Col > 0 ); -/* Move Row so that it is the last row in the matrix. */ + /* 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. */ + /* Move Col so that it is the last column in the matrix. */ if (Col != Size) spcColExchange( Matrix, Col, Size ); -/* Correct Diag pointers. */ + /* Correct Diag pointers. */ if (Row == Col) - SWAP( ElementPtr, Matrix->Diag[Row], Matrix->Diag[Size] ) - else - { Matrix->Diag[Row] = spcFindElementInCol( Matrix, Matrix->FirstInCol+Row, - Row, Row, NO ); - Matrix->Diag[Col] = spcFindElementInCol( Matrix, Matrix->FirstInCol+Col, - Col, Col, NO ); + SWAP( ElementPtr, Matrix->Diag[Row], Matrix->Diag[Size] ); + else { + Matrix->Diag[Row] = spcFindElementInCol(Matrix, + Matrix->FirstInCol+Row, + Row, Row, NO ); + Matrix->Diag[Col] = spcFindElementInCol(Matrix, + Matrix->FirstInCol+Col, + Col, Col, NO ); } -/* - * Delete last row and column of the matrix. - */ -/* Break the column links to every element in the last row. */ + /* 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) + 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. */ @@ -1351,20 +1272,20 @@ ElementPtr spcFindElementInCol(); pLastElement = pLastElement->NextInRow; } -/* Break the row links to every element in the last column. */ + /* 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. */ + 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. */ + ppElement = &pElement->NextInRow; /* Skip element. */ } pLastElement = pLastElement->NextInCol; } -/* Clean up some details. */ + /* Clean up some details. */ Matrix->Size = Size - 1; Matrix->Diag[Size] = NULL; Matrix->FirstInRow[Size] = NULL; @@ -1393,12 +1314,12 @@ ElementPtr spcFindElementInCol(); * 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 + * 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. * * >>> Returns: @@ -1407,35 +1328,33 @@ ElementPtr spcFindElementInCol(); * * >>> Arguments: * eMatrix (char *) - * Pointer to the matrix. - */ + * Pointer to the matrix. */ RealNumber -spPseudoCondition( eMatrix ) - -char *eMatrix; +spPseudoCondition(char *eMatrix) { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register int I; -register ArrayOfElementPtrs Diag; -RealNumber MaxPivot, MinPivot, Mag; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + int I; + ArrayOfElementPtrs Diag; + RealNumber MaxPivot, MinPivot, Mag; -/* Begin `spPseudoCondition'. */ + /* Begin `spPseudoCondition'. */ - ASSERT( IS_SPARSE(Matrix) AND IS_FACTORED(Matrix) ); - if (Matrix->Error == spSINGULAR OR Matrix->Error == spZERO_DIAG) + assert( IS_SPARSE(Matrix) && IS_FACTORED(Matrix) ); + if (Matrix->Error == spSINGULAR || 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; + { + Mag = ELEMENT_MAG( Diag[I] ); + if (Mag > MaxPivot) + MaxPivot = Mag; + else if (Mag < MinPivot) + MinPivot = Mag; } - ASSERT( MaxPivot > 0.0); + assert( MaxPivot > 0.0); return MaxPivot / MinPivot; } #endif @@ -1452,20 +1371,22 @@ RealNumber MaxPivot, MinPivot, Mag; * ESTIMATE CONDITION NUMBER * * 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(). + * 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. + * rather than the lower. This difference should be of no + * importance. * * References: * A.K. Cline, C.B. Moler, G.W. Stewart, J.H. Wilkinson. An estimate @@ -1498,8 +1419,7 @@ RealNumber MaxPivot, MinPivot, Mag; * * >>> Possible errors: * spSINGULAR - * spNO_MEMORY - */ + * spNO_MEMORY */ RealNumber spCondition( eMatrix, NormOfMatrix, pError ) @@ -1508,60 +1428,53 @@ char *eMatrix; RealNumber 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(); + MatrixPtr Matrix = (MatrixPtr)eMatrix; + ElementPtr pElement; + RealVector T, Tm; + 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'. */ + /* Begin `spCondition'. */ - ASSERT( IS_SPARSE(Matrix) AND IS_FACTORED(Matrix) ); + assert( IS_SPARSE(Matrix) && IS_FACTORED(Matrix) ); *pError = Matrix->Error; if (Matrix->Error >= spFATAL) return 0.0; if (NormOfMatrix == 0.0) - { *pError = spSINGULAR; + { + *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. - */ + /* + * 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. */ + /* Forward elimination. Solves Lw = e while choosing e. */ E = 1.0; for (I = 1; I <= Size; I++) - { pPivot = Matrix->Diag[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) ); + { + ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(Wm) ); for (K = Size; K > 0; K--) T[K] *= ScaleFactor; E *= ScaleFactor; Em *= ScaleFactor; @@ -1571,10 +1484,12 @@ RealNumber Linpack, OLeary, InvNormOfInverse, ComplexCondition(); ASp = ABS(T[I] - Em); ASm = ABS(Em + T[I]); -/* Update T for both values of W, minus value is placed in Tm. */ + /* Update T for both values of W, minus value is placed in + Tm. */ pElement = pPivot->NextInCol; while (pElement != NULL) - { Row = pElement->Row; + { + Row = pElement->Row; Tm[Row] = T[Row] - (Wm * pElement->Real); T[Row] -= (Wp * pElement->Real); ASp += ABS(T[Row]); @@ -1582,115 +1497,127 @@ RealNumber Linpack, OLeary, InvNormOfInverse, ComplexCondition(); pElement = pElement->NextInCol; } -/* If minus value causes more growth, overwrite T with its values. */ + /* If minus value causes more growth, overwrite T with its + values. */ if (ASm > ASp) - { T[I] = Wm; + { + T[I] = Wm; pElement = pPivot->NextInCol; while (pElement != NULL) - { T[pElement->Row] = Tm[pElement->Row]; + { + 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. */ + /* 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; + { + for (I = Size; I > 0; I--) T[I] *= ScaleFactor; E *= ScaleFactor; } -/* Backward Substitution. Solves Uy = w.*/ + /* Backward Substitution. Solves Uy = w.*/ for (I = Size; I >= 1; I--) - { pElement = Matrix->Diag[I]->NextInRow; + { + pElement = Matrix->Diag[I]->NextInRow; while (pElement != NULL) - { T[I] -= pElement->Real * T[pElement->Col]; + { + T[I] -= pElement->Real * T[pElement->Col]; pElement = pElement->NextInRow; } if (ABS(T[I]) > SLACK) - { ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) ); + { + 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. */ + /* 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; + { + 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. */ + /* 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*. - */ + /* + * Part 2. + * + * A* z = y where the * represents the transpose. Recall that A = + * LU implies A* = U* L*. */ -/* Forward elimination, U* v = y. */ + /* Forward elimination, U* v = y. */ for (I = 1; I <= Size; I++) - { pElement = Matrix->Diag[I]->NextInRow; + { + pElement = Matrix->Diag[I]->NextInRow; while (pElement != NULL) - { T[pElement->Col] -= T[I] * pElement->Real; + { + T[pElement->Col] -= T[I] * pElement->Real; pElement = pElement->NextInRow; } if (ABS(T[I]) > SLACK) - { ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) ); + { + 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. */ + /* 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; + { + for (I = Size; I > 0; I--) T[I] *= ScaleFactor; ASy *= ScaleFactor; } -/* Backward Substitution, L* z = v. */ + /* Backward Substitution, L* z = v. */ for (I = Size; I >= 1; I--) - { pPivot = Matrix->Diag[I]; + { + pPivot = Matrix->Diag[I]; pElement = pPivot->NextInCol; while (pElement != NULL) - { T[I] -= pElement->Real * T[pElement->Row]; + { + T[I] -= pElement->Real * T[pElement->Row]; pElement = pElement->NextInCol; } T[I] *= pPivot->Real; if (ABS(T[I]) > SLACK) - { ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) ); + { + 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 z. */ + /* Compute 1-norm of T, which now contains z. */ for (ASz = 0.0, I = Size; I > 0; I--) ASz += ABS(T[I]); -#if NOT spCOMPLEX - FREE( Tm ); -#endif - Linpack = ASy / ASz; OLeary = E / MaxY; InvNormOfInverse = MIN( Linpack, OLeary ); return InvNormOfInverse / NormOfMatrix; -#endif /* REAL */ } -#if spCOMPLEX /* * ESTIMATE CONDITION NUMBER * @@ -1719,44 +1646,48 @@ MatrixPtr Matrix; RealNumber NormOfMatrix; int *pError; { -register ElementPtr pElement; -register ComplexVector T, Tm; -register int I, K, Row; -ElementPtr pPivot; -int Size; -RealNumber E, Em, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor; -RealNumber Linpack, OLeary, InvNormOfInverse; -ComplexNumber Wp, Wm; + ElementPtr pElement; + ComplexVector T, Tm; + int I, K, Row; + ElementPtr pPivot; + int Size; + RealNumber E, Em, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor; + RealNumber Linpack, OLeary, InvNormOfInverse; + ComplexNumber Wp, Wm; -/* Begin `ComplexCondition'. */ + /* Begin `ComplexCondition'. */ Size = Matrix->Size; T = (ComplexVector)Matrix->Intermediate; Tm = ALLOC( ComplexNumber, Size+1 ); if (Tm == NULL) - { *pError = spNO_MEMORY; + { + *pError = spNO_MEMORY; return 0.0; } for (I = Size; I > 0; I--) T[I].Real = T[I].Imag = 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. - */ + /* + * 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. */ + /* Forward elimination. Solves Lw = e while choosing e. */ E = 1.0; for (I = 1; I <= Size; I++) - { pPivot = Matrix->Diag[I]; + { + pPivot = Matrix->Diag[I]; if (T[I].Real < 0.0) Em = -E; else Em = E; Wm = T[I]; Wm.Real += Em; ASm = CMPLX_1_NORM( Wm ); CMPLX_MULT_ASSIGN( Wm, *pPivot ); if (CMPLX_1_NORM(Wm) > SLACK) - { ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(Wm) ); + { + ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(Wm) ); for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor ); E *= ScaleFactor; Em *= ScaleFactor; @@ -1768,10 +1699,11 @@ ComplexNumber Wp, Wm; ASp = CMPLX_1_NORM( Wp ); CMPLX_MULT_ASSIGN( Wp, *pPivot ); -/* Update T for both values of W, minus value is placed in Tm. */ + /* Update T for both values of W, minus value is placed in Tm. */ pElement = pPivot->NextInCol; while (pElement != NULL) - { Row = pElement->Row; + { + Row = pElement->Row; /* Cmplx expr: Tm[Row] = T[Row] - (Wp * *pElement). */ CMPLX_MULT_SUBT( Tm[Row], Wm, *pElement, T[Row] ); /* Cmplx expr: T[Row] -= Wp * *pElement. */ @@ -1781,100 +1713,116 @@ ComplexNumber Wp, Wm; pElement = pElement->NextInCol; } -/* If minus value causes more growth, overwrite T with its values. */ + /* If minus value causes more growth, overwrite T with its + values. */ if (ASm > ASp) - { T[I] = Wm; + { + T[I] = Wm; pElement = pPivot->NextInCol; while (pElement != NULL) - { T[pElement->Row] = Tm[pElement->Row]; + { + 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. */ + /* 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 += CMPLX_1_NORM(T[I]); ScaleFactor = 1.0 / (SLACK * ASw); if (ScaleFactor < 0.5) - { for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor ); + { + for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor ); E *= ScaleFactor; } -/* Backward Substitution. Solves Uy = w.*/ + /* Backward Substitution. Solves Uy = w.*/ for (I = Size; I >= 1; I--) - { pElement = Matrix->Diag[I]->NextInRow; + { + pElement = Matrix->Diag[I]->NextInRow; while (pElement != NULL) - { /* Cmplx expr: T[I] -= T[pElement->Col] * *pElement. */ + { + /* Cmplx expr: T[I] -= T[pElement->Col] * *pElement. */ CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Col], *pElement ); pElement = pElement->NextInRow; } if (CMPLX_1_NORM(T[I]) > SLACK) - { ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) ); + { + ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) ); for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor ); E *= ScaleFactor; } } -/* Compute 1-norm of T, which now contains y, and scale ||T|| to 1/SLACK. */ + /* 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 += CMPLX_1_NORM(T[I]); ScaleFactor = 1.0 / (SLACK * ASy); if (ScaleFactor < 0.5) - { for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor ); + { + for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor ); ASy = 1.0 / SLACK; E *= ScaleFactor; } -/* Compute infinity-norm of T for O'Leary's estimate. */ + /* Compute infinity-norm of T for O'Leary's estimate. */ for (MaxY = 0.0, I = Size; I > 0; I--) if (MaxY < CMPLX_1_NORM(T[I])) MaxY = CMPLX_1_NORM(T[I]); -/* - * Part 2. A* z = y where the * represents the transpose. - * Recall that A = LU implies A* = U* L*. - */ + /* Part 2. A* z = y where the * represents the transpose. Recall + * that A = LU implies A* = U* L*. */ -/* Forward elimination, U* v = y. */ + /* Forward elimination, U* v = y. */ for (I = 1; I <= Size; I++) - { pElement = Matrix->Diag[I]->NextInRow; + { + pElement = Matrix->Diag[I]->NextInRow; while (pElement != NULL) - { /* Cmplx expr: T[pElement->Col] -= T[I] * *pElement. */ + { + /* Cmplx expr: T[pElement->Col] -= T[I] * *pElement. */ CMPLX_MULT_SUBT_ASSIGN( T[pElement->Col], T[I], *pElement ); pElement = pElement->NextInRow; } if (CMPLX_1_NORM(T[I]) > SLACK) - { ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) ); + { + ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) ); for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor ); ASy *= ScaleFactor; } } -/* Compute 1-norm of T, which now contains v, and scale ||T|| to 1/SLACK. */ + /* 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 += CMPLX_1_NORM(T[I]); ScaleFactor = 1.0 / (SLACK * ASv); if (ScaleFactor < 0.5) - { for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor ); + { + for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor ); ASy *= ScaleFactor; } -/* Backward Substitution, L* z = v. */ + /* Backward Substitution, L* z = v. */ for (I = Size; I >= 1; I--) - { pPivot = Matrix->Diag[I]; + { + pPivot = Matrix->Diag[I]; pElement = pPivot->NextInCol; while (pElement != NULL) - { /* Cmplx expr: T[I] -= T[pElement->Row] * *pElement. */ + { + /* Cmplx expr: T[I] -= T[pElement->Row] * *pElement. */ CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Row], *pElement ); pElement = pElement->NextInCol; } CMPLX_MULT_ASSIGN( T[I], *pPivot ); if (CMPLX_1_NORM(T[I]) > SLACK) - { ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) ); + { + ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) ); for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor ); ASy *= ScaleFactor; } } -/* Compute 1-norm of T, which now contains z. */ + /* Compute 1-norm of T, which now contains z. */ for (ASz = 0.0, I = Size; I > 0; I--) ASz += CMPLX_1_NORM(T[I]); FREE( Tm ); @@ -1884,7 +1832,6 @@ ComplexNumber Wp, Wm; InvNormOfInverse = MIN( Linpack, OLeary ); return InvNormOfInverse / NormOfMatrix; } -#endif /* spCOMPLEX */ @@ -1911,42 +1858,44 @@ spNorm( eMatrix ) char *eMatrix; { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register ElementPtr pElement; -register int I; -RealNumber Max = 0.0, AbsRowSum; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + ElementPtr pElement; + int I; + RealNumber Max = 0.0, AbsRowSum; -/* Begin `spNorm'. */ - ASSERT( IS_SPARSE(Matrix) AND NOT IS_FACTORED(Matrix) ); - if (NOT Matrix->RowsLinked) spcLinkRows( Matrix ); + /* Begin `spNorm'. */ + assert( IS_SPARSE(Matrix) && !IS_FACTORED(Matrix) ); + if (!Matrix->RowsLinked) spcLinkRows( Matrix ); -/* Compute row sums. */ -#if REAL - if (NOT Matrix->Complex) - { for (I = Matrix->Size; I > 0; I--) - { pElement = Matrix->FirstInRow[I]; + /* Compute row sums. */ + if (!Matrix->Complex) + { + for (I = Matrix->Size; I > 0; I--) + { + pElement = Matrix->FirstInRow[I]; AbsRowSum = 0.0; while (pElement != NULL) - { AbsRowSum += ABS( pElement->Real ); + { + AbsRowSum += ABS( pElement->Real ); pElement = pElement->NextInRow; } if (Max < AbsRowSum) Max = AbsRowSum; } } -#endif -#if spCOMPLEX - if (Matrix->Complex) - { for (I = Matrix->Size; I > 0; I--) - { pElement = Matrix->FirstInRow[I]; + else + { + for (I = Matrix->Size; I > 0; I--) + { + pElement = Matrix->FirstInRow[I]; AbsRowSum = 0.0; while (pElement != NULL) - { AbsRowSum += CMPLX_1_NORM( *pElement ); + { + AbsRowSum += CMPLX_1_NORM( *pElement ); pElement = pElement->NextInRow; } if (Max < AbsRowSum) Max = AbsRowSum; } } -#endif return Max; } #endif /* CONDITION */ @@ -1961,20 +1910,21 @@ RealNumber Max = 0.0, AbsRowSum; * STABILITY OF FACTORIZATION * * The following routines are used to gauge the stability of a - * factorization. If the factorization is determined to be too unstable, - * then the matrix should be reordered. The routines compute quantities - * that are needed in the computation of a bound on the error attributed - * to any one element in the matrix during the factorization. In other - * words, there is a matrix E = [e_ij] of error terms such that A+E = LU. - * This routine finds a bound on |e_ij|. Erisman & Reid [1] showed that - * |e_ij| < 3.01 u rho m_ij, where u is the machine rounding unit, - * rho = max a_ij where the max is taken over every row i, column j, and - * step k, and m_ij is the number of multiplications required in the - * computation of l_ij if i > j or u_ij otherwise. Barlow [2] showed that - * rho < max_i || l_i ||_p max_j || u_j ||_q where 1/p + 1/q = 1. + * factorization. If the factorization is determined to be too + * unstable, then the matrix should be reordered. The routines + * compute quantities that are needed in the computation of a bound + * on the error attributed to any one element in the matrix during + * the factorization. In other words, there is a matrix E = [e_ij] + * of error terms such that A+E = LU. This routine finds a bound on + * |e_ij|. Erisman & Reid [1] showed that |e_ij| < 3.01 u rho m_ij, + * where u is the machine rounding unit, rho = max a_ij where the max + * is taken over every row i, column j, and step k, and m_ij is the + * number of multiplications required in the computation of l_ij if i + * > j or u_ij otherwise. Barlow [2] showed that rho < max_i || l_i + * ||_p max_j || u_j ||_q where 1/p + 1/q = 1. * - * The first routine finds the magnitude on the largest element in the - * matrix. If the matrix has not yet been factored, the largest + * The first routine finds the magnitude on the largest element in + * the matrix. If the matrix has not yet been factored, the largest * element is found by direct search. If the matrix is factored, a * bound on the largest element in any of the reduced submatrices is * computed using Barlow with p = oo and q = 1. The ratio of these @@ -1989,14 +1939,15 @@ RealNumber Max = 0.0, AbsRowSum; * * Using only the size of the matrix as an upper bound on m_ij and * Barlow's bound, the user can estimate the size of the matrix error - * terms e_ij using the bound of Erisman and Reid. The second routine - * computes a tighter bound (with more work) based on work by Gear - * [3], |e_ij| < 1.01 u rho (t c^3 + (1 + t)c^2) where t is the + * terms e_ij using the bound of Erisman and Reid. The second + * routine computes a tighter bound (with more work) based on work by + * Gear [3], |e_ij| < 1.01 u rho (t c^3 + (1 + t)c^2) where t is the * threshold and c is the maximum number of off-diagonal elements in * any row of L. The expensive part of computing this bound is - * determining the maximum number of off-diagonals in L, which changes - * only when the order of the matrix changes. This number is computed - * and saved, and only recomputed if the matrix is reordered. + * determining the maximum number of off-diagonals in L, which + * changes only when the order of the matrix changes. This number is + * computed and saved, and only recomputed if the matrix is + * reordered. * * [1] A. M. Erisman, J. K. Reid. Monitoring the stability of the * triangular factorization of a sparse matrix. Numerische @@ -2007,8 +1958,7 @@ RealNumber Max = 0.0, AbsRowSum; * and Statistical Computing." Vol. 7, No. 1, January 1986, pp 166-168. * * [3] I. S. Duff, A. M. Erisman, J. K. Reid. "Direct Methods for Sparse - * Matrices." Oxford 1986. pp 99. - */ + * Matrices." Oxford 1986. pp 99. */ /* * LARGEST ELEMENT IN MATRIX @@ -2028,98 +1978,110 @@ spLargestElement( eMatrix ) char *eMatrix; { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register int I; -RealNumber Mag, AbsColSum, Max = 0.0, MaxRow = 0.0, MaxCol = 0.0; -RealNumber Pivot; -ComplexNumber cPivot; -register ElementPtr pElement, pDiag; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + int I; + RealNumber Mag, AbsColSum, Max = 0.0, MaxRow = 0.0, MaxCol = 0.0; + RealNumber Pivot; + ComplexNumber cPivot; + ElementPtr pElement, pDiag; -/* Begin `spLargestElement'. */ - ASSERT( IS_SPARSE(Matrix) ); + /* Begin `spLargestElement'. */ + assert( IS_SPARSE(Matrix) ); -#if REAL - if (Matrix->Factored AND NOT Matrix->Complex) - { if (Matrix->Error == spSINGULAR) return 0.0; + if (Matrix->Factored && !Matrix->Complex) + { + if (Matrix->Error == spSINGULAR) return 0.0; -/* Find the bound on the size of the largest element over all factorization. */ + /* Find the bound on the size of the largest element over all + factorization. */ for (I = 1; I <= Matrix->Size; I++) - { pDiag = Matrix->Diag[I]; + { + pDiag = Matrix->Diag[I]; -/* Lower triangular matrix. */ + /* Lower triangular matrix. */ Pivot = 1.0 / pDiag->Real; Mag = ABS( Pivot ); if (Mag > MaxRow) MaxRow = Mag; pElement = Matrix->FirstInRow[I]; while (pElement != pDiag) - { Mag = ABS( pElement->Real ); + { + Mag = ABS( pElement->Real ); if (Mag > MaxRow) MaxRow = Mag; pElement = pElement->NextInRow; } -/* Upper triangular matrix. */ + /* Upper triangular matrix. */ pElement = Matrix->FirstInCol[I]; AbsColSum = 1.0; /* Diagonal of U is unity. */ while (pElement != pDiag) - { AbsColSum += ABS( pElement->Real ); + { + AbsColSum += ABS( pElement->Real ); pElement = pElement->NextInCol; } if (AbsColSum > MaxCol) MaxCol = AbsColSum; } } - else if (NOT Matrix->Complex) - { for (I = 1; I <= Matrix->Size; I++) - { pElement = Matrix->FirstInCol[I]; + else if (!Matrix->Complex) + { + for (I = 1; I <= Matrix->Size; I++) + { + pElement = Matrix->FirstInCol[I]; while (pElement != NULL) - { Mag = ABS( pElement->Real ); + { + Mag = ABS( pElement->Real ); if (Mag > Max) Max = Mag; pElement = pElement->NextInCol; } } return Max; } -#endif -#if spCOMPLEX - if (Matrix->Factored AND Matrix->Complex) - { if (Matrix->Error == spSINGULAR) return 0.0; + if (Matrix->Factored && Matrix->Complex) + { + if (Matrix->Error == spSINGULAR) return 0.0; -/* Find the bound on the size of the largest element over all factorization. */ + /* Find the bound on the size of the largest element over all + factorization. */ for (I = 1; I <= Matrix->Size; I++) - { pDiag = Matrix->Diag[I]; + { + pDiag = Matrix->Diag[I]; -/* Lower triangular matrix. */ + /* Lower triangular matrix. */ CMPLX_RECIPROCAL( cPivot, *pDiag ); Mag = CMPLX_1_NORM( cPivot ); if (Mag > MaxRow) MaxRow = Mag; pElement = Matrix->FirstInRow[I]; while (pElement != pDiag) - { Mag = CMPLX_1_NORM( *pElement ); + { + Mag = CMPLX_1_NORM( *pElement ); if (Mag > MaxRow) MaxRow = Mag; pElement = pElement->NextInRow; } -/* Upper triangular matrix. */ + /* Upper triangular matrix. */ pElement = Matrix->FirstInCol[I]; AbsColSum = 1.0; /* Diagonal of U is unity. */ while (pElement != pDiag) - { AbsColSum += CMPLX_1_NORM( *pElement ); + { + AbsColSum += CMPLX_1_NORM( *pElement ); pElement = pElement->NextInCol; } if (AbsColSum > MaxCol) MaxCol = AbsColSum; } } else if (Matrix->Complex) - { for (I = 1; I <= Matrix->Size; I++) - { pElement = Matrix->FirstInCol[I]; + { + for (I = 1; I <= Matrix->Size; I++) + { + pElement = Matrix->FirstInCol[I]; while (pElement != NULL) - { Mag = CMPLX_1_NORM( *pElement ); + { + Mag = CMPLX_1_NORM( *pElement ); if (Mag > Max) Max = Mag; pElement = pElement->NextInCol; } } return Max; } -#endif return MaxRow * MaxCol; } @@ -2148,24 +2110,27 @@ spRoundoff( eMatrix, Rho ) char *eMatrix; RealNumber Rho; { -MatrixPtr Matrix = (MatrixPtr)eMatrix; -register ElementPtr pElement; -register int Count, I, MaxCount = 0; -RealNumber Reid, Gear; + MatrixPtr Matrix = (MatrixPtr)eMatrix; + ElementPtr pElement; + int Count, I, MaxCount = 0; + RealNumber Reid, Gear; -/* Begin `spRoundoff'. */ - ASSERT( IS_SPARSE(Matrix) AND IS_FACTORED(Matrix) ); + /* Begin `spRoundoff'. */ + assert( IS_SPARSE(Matrix) && IS_FACTORED(Matrix) ); -/* Compute Barlow's bound if it is not given. */ + /* Compute Barlow's bound if it is not given. */ if (Rho < 0.0) Rho = spLargestElement( eMatrix ); -/* Find the maximum number of off-diagonals in L if not previously computed. */ + /* Find the maximum number of off-diagonals in L if not previously computed. */ if (Matrix->MaxRowCountInLowerTri < 0) - { for (I = Matrix->Size; I > 0; I--) - { pElement = Matrix->FirstInRow[I]; + { + for (I = Matrix->Size; I > 0; I--) + { + pElement = Matrix->FirstInRow[I]; Count = 0; while (pElement->Col < I) - { Count++; + { + Count++; pElement = pElement->NextInRow; } if (Count > MaxCount) MaxCount = Count; @@ -2174,7 +2139,7 @@ RealNumber Reid, Gear; } else MaxCount = Matrix->MaxRowCountInLowerTri; -/* Compute error bound. */ + /* Compute error bound. */ Gear = 1.01*((MaxCount + 1) * Matrix->RelThreshold + 1.0) * SQR(MaxCount); Reid = 3.01 * Matrix->Size; @@ -2215,13 +2180,14 @@ spErrorMessage( eMatrix, Stream, Originator ) char *eMatrix, *Originator; FILE *Stream; { -int Row, Col, Error; + int Row, Col, Error; -/* Begin `spErrorMessage'. */ + /* Begin `spErrorMessage'. */ if (eMatrix == NULL) Error = spNO_MEMORY; else - { ASSERT(((MatrixPtr)eMatrix)->ID == SPARSE_ID); + { + assert(((MatrixPtr)eMatrix)->ID == SPARSE_ID); Error = ((MatrixPtr)eMatrix)->Error; } @@ -2232,29 +2198,34 @@ int Row, Col, Error; fprintf( Stream, "fatal error, "); else fprintf( Stream, "warning, "); -/* - * Print particular error message. - * Do not use switch statement because error codes may not be unique. - */ + + /* Print particular error message. Do not use switch statement + * because error codes may not be unique. */ if (Error == spPANIC) fprintf( Stream, "Sparse called improperly.\n"); else if (Error == spNO_MEMORY) fprintf( Stream, "insufficient memory available.\n"); else if (Error == spSINGULAR) - { spWhereSingular( eMatrix, &Row, &Col ); + { + spWhereSingular( eMatrix, &Row, &Col ); fprintf( Stream, "singular matrix detected at row %d and column %d.\n", Row, Col); } else if (Error == spZERO_DIAG) - { spWhereSingular( eMatrix, &Row, &Col ); + { + spWhereSingular( eMatrix, &Row, &Col ); fprintf( Stream, "zero diagonal detected at row %d and column %d.\n", Row, Col); } else if (Error == spSMALL_PIVOT) - { fprintf( Stream, - "unable to find a pivot that is larger than absolute threshold.\n"); + { + fprintf( Stream, + "unable to find a pivot that is larger than absolute threshold.\n"); + } + else + { + abort(); } - else ABORT(); return; } #endif /* DOCUMENTATION */