more adaption to existing solution

This commit is contained in:
dwarning 2023-06-16 15:01:19 +02:00
parent feb1cef086
commit 4a3b629439
7 changed files with 201 additions and 196 deletions

View File

@ -133,9 +133,9 @@ spCreate(
spError *pError spError *pError
) )
{ {
register unsigned SizePlusOne; unsigned SizePlusOne;
register MatrixPtr Matrix; MatrixPtr Matrix;
register int I; int I;
int AllocatedSize; int AllocatedSize;
/* Begin `spCreate'. */ /* Begin `spCreate'. */
@ -201,7 +201,7 @@ int AllocatedSize;
Matrix->ElementsRemaining = 0; Matrix->ElementsRemaining = 0;
Matrix->FillinsRemaining = 0; Matrix->FillinsRemaining = 0;
RecordAllocation( Matrix, (void *)Matrix ); RecordAllocation( Matrix, Matrix );
if (Matrix->Error == spNO_MEMORY) goto MemoryError; /* FIXME: Use of memory after free */ if (Matrix->Error == spNO_MEMORY) goto MemoryError; /* FIXME: Use of memory after free */
/* Take out the trash. */ /* Take out the trash. */
@ -322,7 +322,7 @@ ElementPtr pElement;
/* Allocate block of MatrixElements if necessary. */ /* Allocate block of MatrixElements if necessary. */
if (Matrix->ElementsRemaining == 0) if (Matrix->ElementsRemaining == 0)
{ pElement = SP_MALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION); { pElement = SP_MALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, (void *)pElement ); RecordAllocation( Matrix, pElement );
if (Matrix->Error == spNO_MEMORY) return NULL; if (Matrix->Error == spNO_MEMORY) return NULL;
Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION; Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION;
Matrix->NextAvailElement = pElement; Matrix->NextAvailElement = pElement;
@ -383,21 +383,21 @@ ElementPtr pElement;
/* Allocate block of MatrixElements for elements. */ /* Allocate block of MatrixElements for elements. */
pElement = SP_MALLOC(struct MatrixElement, InitialNumberOfElements); pElement = SP_MALLOC(struct MatrixElement, InitialNumberOfElements);
RecordAllocation( Matrix, (void *)pElement ); RecordAllocation( Matrix, pElement );
if (Matrix->Error == spNO_MEMORY) return; if (Matrix->Error == spNO_MEMORY) return;
Matrix->ElementsRemaining = InitialNumberOfElements; Matrix->ElementsRemaining = InitialNumberOfElements;
Matrix->NextAvailElement = pElement; Matrix->NextAvailElement = pElement;
/* Allocate block of MatrixElements for fill-ins. */ /* Allocate block of MatrixElements for fill-ins. */
pElement = SP_MALLOC(struct MatrixElement, NumberOfFillinsExpected); pElement = SP_MALLOC(struct MatrixElement, NumberOfFillinsExpected);
RecordAllocation( Matrix, (void *)pElement ); RecordAllocation( Matrix, pElement );
if (Matrix->Error == spNO_MEMORY) return; if (Matrix->Error == spNO_MEMORY) return;
Matrix->FillinsRemaining = NumberOfFillinsExpected; Matrix->FillinsRemaining = NumberOfFillinsExpected;
Matrix->NextAvailFillin = pElement; Matrix->NextAvailFillin = pElement;
/* Allocate a fill-in list structure. */ /* Allocate a fill-in list structure. */
Matrix->FirstFillinListNode = SP_MALLOC(struct FillinListNodeStruct,1); Matrix->FirstFillinListNode = SP_MALLOC(struct FillinListNodeStruct,1);
RecordAllocation( Matrix, (void *)Matrix->FirstFillinListNode ); RecordAllocation( Matrix, Matrix->FirstFillinListNode );
if (Matrix->Error == spNO_MEMORY) return; if (Matrix->Error == spNO_MEMORY) return;
Matrix->LastFillinListNode = Matrix->FirstFillinListNode; Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
@ -465,14 +465,14 @@ ElementPtr pFillins;
{ {
/* Allocate block of fill-ins. */ /* Allocate block of fill-ins. */
pFillins = SP_MALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION); pFillins = SP_MALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION);
RecordAllocation( Matrix, (void *)pFillins ); RecordAllocation( Matrix, pFillins );
if (Matrix->Error == spNO_MEMORY) return NULL; if (Matrix->Error == spNO_MEMORY) return NULL;
Matrix->FillinsRemaining = ELEMENTS_PER_ALLOCATION; Matrix->FillinsRemaining = ELEMENTS_PER_ALLOCATION;
Matrix->NextAvailFillin = pFillins; Matrix->NextAvailFillin = pFillins;
/* Allocate a fill-in list structure. */ /* Allocate a fill-in list structure. */
pListNode->Next = SP_MALLOC(struct FillinListNodeStruct,1); pListNode->Next = SP_MALLOC(struct FillinListNodeStruct,1);
RecordAllocation( Matrix, (void *)pListNode->Next ); RecordAllocation( Matrix, pListNode->Next );
if (Matrix->Error == spNO_MEMORY) return NULL; if (Matrix->Error == spNO_MEMORY) return NULL;
Matrix->LastFillinListNode = pListNode = pListNode->Next; Matrix->LastFillinListNode = pListNode = pListNode->Next;
@ -505,7 +505,7 @@ ElementPtr pFillins;
* >>> Arguments: * >>> Arguments:
* Matrix <input> (MatrixPtr) * Matrix <input> (MatrixPtr)
* Pointer to the matrix. * Pointer to the matrix.
* AllocatedPtr <input> (void *) * AllocatedPtr <input>
* The pointer returned by malloc or calloc. These pointers are saved in * The pointer returned by malloc or calloc. These pointers are saved in
* a list so that they can be easily freed. * a list so that they can be easily freed.
* *
@ -573,8 +573,8 @@ RecordAllocation(
static void static void
AllocateBlockOfAllocationList( MatrixPtr Matrix ) AllocateBlockOfAllocationList( MatrixPtr Matrix )
{ {
register int I; int I;
register AllocationListPtr ListPtr; AllocationListPtr ListPtr;
/* Begin `AllocateBlockOfAllocationList'. */ /* Begin `AllocateBlockOfAllocationList'. */
/* Allocate block of records for allocation list. */ /* Allocate block of records for allocation list. */
@ -596,7 +596,7 @@ register AllocationListPtr ListPtr;
} }
/* Record allocation of space for allocation list on allocation list. */ /* Record allocation of space for allocation list on allocation list. */
Matrix->TopOfAllocationList->AllocatedPtr = (void *)ListPtr; Matrix->TopOfAllocationList->AllocatedPtr = ListPtr;
Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION; Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION;
return; return;
@ -629,7 +629,7 @@ register AllocationListPtr ListPtr;
void void
spDestroy( MatrixPtr Matrix ) spDestroy( MatrixPtr Matrix )
{ {
register AllocationListPtr ListPtr, NextListPtr; AllocationListPtr ListPtr, NextListPtr;
/* Begin `spDestroy'. */ /* Begin `spDestroy'. */
ASSERT_IS_SPARSE( Matrix ); ASSERT_IS_SPARSE( Matrix );

View File

@ -104,8 +104,8 @@ static void EnlargeMatrix( MatrixPtr, int );
void void
spClear( MatrixPtr Matrix ) spClear( MatrixPtr Matrix )
{ {
register ElementPtr pElement; ElementPtr pElement;
register int I; int I;
/* Begin `spClear'. */ /* Begin `spClear'. */
ASSERT_IS_SPARSE( Matrix ); ASSERT_IS_SPARSE( Matrix );
@ -184,7 +184,7 @@ spFindElement(
int Col int Col
) )
{ {
register ElementPtr pElement; ElementPtr pElement;
int StartAt = 0; int StartAt = 0;
long int Min = LARGEST_LONG_INTEGER; long int Min = LARGEST_LONG_INTEGER;
#define BorderRight 0 /* Start at left border, move right. */ #define BorderRight 0 /* Start at left border, move right. */
@ -449,7 +449,7 @@ Translate(
int *Col int *Col
) )
{ {
register int IntRow, IntCol, ExtRow, ExtCol; int IntRow, IntCol, ExtRow, ExtCol;
/* Begin `Translate'. */ /* Begin `Translate'. */
ExtRow = *Row; ExtRow = *Row;
@ -649,7 +649,7 @@ spGetQuad(
OR (Template->Element4Negated == NULL) OR (Template->Element4Negated == NULL)
) return spNO_MEMORY; ) return spNO_MEMORY;
if (Template->Element1 == &((MatrixPtr)Matrix)->TrashCan.Real) if (Template->Element1 == &(Matrix)->TrashCan.Real)
SWAP( RealNumber *, Template->Element1, Template->Element2 ); SWAP( RealNumber *, Template->Element1, Template->Element2 );
return spOKAY; return spOKAY;
@ -756,10 +756,10 @@ spGetOnes(
ElementPtr ElementPtr
spcFindDiag( spcFindDiag(
MatrixPtr Matrix, MatrixPtr Matrix,
register int Index int Index
) )
{ {
register ElementPtr pElement; ElementPtr pElement;
/* Begin `spcFindDiag'. */ /* Begin `spcFindDiag'. */
pElement = Matrix->FirstInCol[Index]; pElement = Matrix->FirstInCol[Index];
@ -823,18 +823,18 @@ ElementPtr
spcCreateElement( spcCreateElement(
MatrixPtr Matrix, MatrixPtr Matrix,
int Row, int Row,
register int Col, int Col,
register ElementPtr *ppToLeft, ElementPtr *ppToLeft,
register ElementPtr *ppAbove, ElementPtr *ppAbove,
BOOLEAN Fillin BOOLEAN Fillin
) )
{ {
register ElementPtr pElement, pCreatedElement; ElementPtr pElement, pCreatedElement;
/* Begin `spcCreateElement'. */ /* Begin `spcCreateElement'. */
/* Find element immediately above the desired element. */ /* Find element immediately above the desired element. */
pElement = *ppAbove; pElement = *ppAbove; /* FIXME: Dereference of null pointer */
while ((pElement != NULL) AND (pElement->Row < Row)) while ((pElement != NULL) AND (pElement->Row < Row))
{ ppAbove = &pElement->NextInCol; { ppAbove = &pElement->NextInCol;
pElement = *ppAbove; pElement = *ppAbove;
@ -937,7 +937,7 @@ register ElementPtr pElement, pCreatedElement;
* currently being operated upon. * currently being operated upon.
* FirstInRowArray (ArrayOfElementPtrs) * FirstInRowArray (ArrayOfElementPtrs)
* A pointer to the FirstInRow array. Same as Matrix->FirstInRow but * 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. * use.
* Col (int) * Col (int)
* Column currently being operated upon. * Column currently being operated upon.
@ -946,9 +946,9 @@ register ElementPtr pElement, pCreatedElement;
void void
spcLinkRows( MatrixPtr Matrix ) spcLinkRows( MatrixPtr Matrix )
{ {
register ElementPtr pElement, *FirstInRowEntry; ElementPtr pElement, *FirstInRowEntry;
register ArrayOfElementPtrs FirstInRowArray; ArrayOfElementPtrs FirstInRowArray;
register int Col; int Col;
/* Begin `spcLinkRows'. */ /* Begin `spcLinkRows'. */
FirstInRowArray = Matrix->FirstInRow; FirstInRowArray = Matrix->FirstInRow;
@ -998,10 +998,10 @@ register int Col;
static void static void
EnlargeMatrix( EnlargeMatrix(
MatrixPtr Matrix, MatrixPtr Matrix,
register int NewSize int NewSize
) )
{ {
register int I, OldAllocatedSize = Matrix->AllocatedSize; int I, OldAllocatedSize = Matrix->AllocatedSize;
/* Begin `EnlargeMatrix'. */ /* Begin `EnlargeMatrix'. */
Matrix->Size = NewSize; Matrix->Size = NewSize;
@ -1087,10 +1087,10 @@ register int I, OldAllocatedSize = Matrix->AllocatedSize;
static void static void
ExpandTranslationArrays( ExpandTranslationArrays(
MatrixPtr Matrix, MatrixPtr Matrix,
register int NewSize int NewSize
) )
{ {
register int I, OldAllocatedSize = Matrix->AllocatedExtSize; int I, OldAllocatedSize = Matrix->AllocatedExtSize;
/* Begin `ExpandTranslationArrays'. */ /* Begin `ExpandTranslationArrays'. */
Matrix->ExtSize = NewSize; Matrix->ExtSize = NewSize;
@ -1176,7 +1176,7 @@ spInitialize(
) )
) )
{ {
register ElementPtr pElement; ElementPtr pElement;
int J, Error, Col; int J, Error, Col;
/* Begin `spInitialize'. */ /* Begin `spInitialize'. */

View File

@ -228,6 +228,11 @@ RealNumber LargestInCol, FindLargestInCol();
/* 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++) for (Step = 1; Step <= Size; Step++)
{ pPivot = Matrix->Diag[Step]; { pPivot = Matrix->Diag[Step];
if (!pPivot) {
fprintf(stderr, "Warning: spfactor.c, 230, Pivot for step = %d not found\n", Step);
Matrix->NeedsOrdering = YES;
break; /* for loop */
}
LargestInCol = FindLargestInCol(pPivot->NextInCol); LargestInCol = FindLargestInCol(pPivot->NextInCol);
if ((LargestInCol * RelThreshold < ELEMENT_MAG(pPivot))) if ((LargestInCol * RelThreshold < ELEMENT_MAG(pPivot)))
{ if (Matrix->Complex) { if (Matrix->Complex)
@ -337,9 +342,9 @@ spError
spFactor( MatrixPtr Matrix ) spFactor( MatrixPtr Matrix )
{ {
#if REAL #if REAL
register ElementPtr pElement; ElementPtr pElement;
register ElementPtr pColumn; ElementPtr pColumn;
register int Step, Size; int Step, Size;
RealNumber Mult; RealNumber Mult;
#endif #endif
/* Begin `spFactor'. */ /* Begin `spFactor'. */
@ -348,7 +353,7 @@ RealNumber Mult;
ASSERT_IS_NOT_FACTORED( Matrix ); ASSERT_IS_NOT_FACTORED( Matrix );
if (Matrix->NeedsOrdering) if (Matrix->NeedsOrdering)
{ return spOrderAndFactor( Matrix, (RealVector)NULL, { return spOrderAndFactor( Matrix, NULL,
0.0, 0.0, DIAG_PIVOTING_AS_DEFAULT ); 0.0, 0.0, DIAG_PIVOTING_AS_DEFAULT );
} }
if (NOT Matrix->Partitioned) spPartition( Matrix, spDEFAULT_PARTITION ); if (NOT Matrix->Partitioned) spPartition( Matrix, spDEFAULT_PARTITION );
@ -359,6 +364,11 @@ RealNumber Mult;
#if REAL #if REAL
Size = Matrix->Size; Size = Matrix->Size;
if (Size == 0) {
Matrix->Factored = YES;
return (Matrix->Error = spOKAY);
}
if (Matrix->Diag[1]->Real == 0.0) return ZeroPivot( Matrix, 1 ); if (Matrix->Diag[1]->Real == 0.0) return ZeroPivot( Matrix, 1 );
Matrix->Diag[1]->Real = 1.0 / Matrix->Diag[1]->Real; Matrix->Diag[1]->Real = 1.0 / Matrix->Diag[1]->Real;
@ -366,7 +376,7 @@ RealNumber Mult;
for (Step = 2; Step <= Size; Step++) for (Step = 2; Step <= Size; Step++)
{ if (Matrix->DoRealDirect[Step]) { if (Matrix->DoRealDirect[Step])
{ /* Update column using direct addressing scatter-gather. */ { /* Update column using direct addressing scatter-gather. */
register RealNumber *Dest = (RealNumber *)Matrix->Intermediate; RealNumber *Dest = (RealNumber *)Matrix->Intermediate;
/* Scatter. */ /* Scatter. */
pElement = Matrix->FirstInCol[Step]; pElement = Matrix->FirstInCol[Step];
@ -398,7 +408,7 @@ RealNumber Mult;
} }
else else
{ /* Update column using indirect addressing scatter-gather. */ { /* Update column using indirect addressing scatter-gather. */
register RealNumber **pDest = (RealNumber **)Matrix->Intermediate; RealNumber **pDest = (RealNumber **)Matrix->Intermediate;
/* Scatter. */ /* Scatter. */
pElement = Matrix->FirstInCol[Step]; pElement = Matrix->FirstInCol[Step];
@ -456,15 +466,21 @@ RealNumber Mult;
static int static int
FactorComplexMatrix( MatrixPtr Matrix ) FactorComplexMatrix( MatrixPtr Matrix )
{ {
register ElementPtr pElement; ElementPtr pElement;
register ElementPtr pColumn; ElementPtr pColumn;
register int Step, Size; int Step, Size;
ComplexNumber Mult, Pivot; ComplexNumber Mult, Pivot;
/* Begin `FactorComplexMatrix'. */ /* Begin `FactorComplexMatrix'. */
ASSERT(Matrix->Complex); ASSERT(Matrix->Complex);
Size = Matrix->Size; Size = Matrix->Size;
if (Size == 0) {
Matrix->Factored = YES;
return (Matrix->Error = spOKAY);
}
pElement = Matrix->Diag[1]; pElement = Matrix->Diag[1];
if (ELEMENT_MAG(pElement) == 0.0) return ZeroPivot( Matrix, 1 ); if (ELEMENT_MAG(pElement) == 0.0) return ZeroPivot( Matrix, 1 );
/* Cmplx expr: *pPivot = 1.0 / *pPivot. */ /* Cmplx expr: *pPivot = 1.0 / *pPivot. */
@ -474,7 +490,7 @@ ComplexNumber Mult, Pivot;
for (Step = 2; Step <= Size; Step++) for (Step = 2; Step <= Size; Step++)
{ if (Matrix->DoCmplxDirect[Step]) { if (Matrix->DoCmplxDirect[Step])
{ /* Update column using direct addressing scatter-gather. */ { /* Update column using direct addressing scatter-gather. */
register ComplexNumber *Dest; ComplexNumber *Dest;
Dest = (ComplexNumber *)Matrix->Intermediate; Dest = (ComplexNumber *)Matrix->Intermediate;
/* Scatter. */ /* Scatter. */
@ -512,7 +528,7 @@ ComplexNumber Mult, Pivot;
} }
else else
{ /* Update column using direct addressing scatter-gather. */ { /* Update column using direct addressing scatter-gather. */
register ComplexNumber **pDest; ComplexNumber **pDest;
pDest = (ComplexNumber **)Matrix->Intermediate; pDest = (ComplexNumber **)Matrix->Intermediate;
/* Scatter. */ /* Scatter. */
@ -593,10 +609,10 @@ spPartition(
int Mode int Mode
) )
{ {
register ElementPtr pElement, pColumn; ElementPtr pElement, pColumn;
register int Step, Size; int Step, Size;
register int *Nc, *No; int *Nc, *No;
register long *Nm; long *Nm;
#if spCOMPLEX #if spCOMPLEX
BOOLEAN *DoCmplxDirect; BOOLEAN *DoCmplxDirect;
#endif #endif
@ -831,12 +847,12 @@ int Size;
static void static void
CountMarkowitz( CountMarkowitz(
MatrixPtr Matrix, MatrixPtr Matrix,
register RealVector RHS, RealVector RHS,
int Step int Step
) )
{ {
register int Count, I, Size = Matrix->Size; int Count, I, Size = Matrix->Size;
register ElementPtr pElement; ElementPtr pElement;
int ExtRow; int ExtRow;
/* Begin `CountMarkowitz'. */ /* Begin `CountMarkowitz'. */
@ -944,9 +960,9 @@ MarkowitzProducts(
int Step int Step
) )
{ {
register int I, *pMarkowitzRow, *pMarkowitzCol; int I, *pMarkowitzRow, *pMarkowitzCol;
register long Product, *pMarkowitzProduct; long Product, *pMarkowitzProduct;
register int Size = Matrix->Size; int Size = Matrix->Size;
double fProduct; double fProduct;
/* Begin `MarkowitzProducts'. */ /* Begin `MarkowitzProducts'. */
@ -1032,11 +1048,7 @@ SearchForPivot(
BOOLEAN DiagPivoting BOOLEAN DiagPivoting
) )
{ {
register ElementPtr ChosenPivot; ElementPtr ChosenPivot;
ElementPtr SearchForSingleton();
ElementPtr QuicklySearchDiagonal();
ElementPtr SearchDiagonal();
ElementPtr SearchEntireMatrix();
/* Begin `SearchForPivot'. */ /* Begin `SearchForPivot'. */
@ -1133,11 +1145,11 @@ SearchForSingleton(
int Step int Step
) )
{ {
register ElementPtr ChosenPivot; ElementPtr ChosenPivot;
register int I; int I;
register long *pMarkowitzProduct; long *pMarkowitzProduct;
int Singletons; int Singletons;
RealNumber PivotMag, FindBiggestInColExclude(); RealNumber PivotMag;
/* Begin `SearchForSingleton'. */ /* Begin `SearchForSingleton'. */
/* Initialize pointer that is to scan through MarkowitzProduct vector. */ /* Initialize pointer that is to scan through MarkowitzProduct vector. */
@ -1355,8 +1367,8 @@ QuicklySearchDiagonal(
int Step int Step
) )
{ {
register long MinMarkowitzProduct, *pMarkowitzProduct; long MinMarkowitzProduct, *pMarkowitzProduct;
register ElementPtr pDiag, pOtherInRow, pOtherInCol; ElementPtr pDiag, pOtherInRow, pOtherInCol;
int I, NumberOfTies; int I, NumberOfTies;
ElementPtr ChosenPivot, TiedElements[MAX_MARKOWITZ_TIES + 1]; ElementPtr ChosenPivot, TiedElements[MAX_MARKOWITZ_TIES + 1];
RealNumber Magnitude, LargestInCol, Ratio, MaxRatio; RealNumber Magnitude, LargestInCol, Ratio, MaxRatio;
@ -1554,12 +1566,11 @@ QuicklySearchDiagonal(
int Step int Step
) )
{ {
register long MinMarkowitzProduct, *pMarkowitzProduct; long MinMarkowitzProduct, *pMarkowitzProduct;
register ElementPtr pDiag; ElementPtr pDiag;
int I; int I;
ElementPtr ChosenPivot, pOtherInRow, pOtherInCol; ElementPtr ChosenPivot, pOtherInRow, pOtherInCol;
RealNumber Magnitude, LargestInCol, LargestOffDiagonal; RealNumber Magnitude, LargestInCol, LargestOffDiagonal;
RealNumber FindBiggestInColExclude();
/* Begin `QuicklySearchDiagonal'. */ /* Begin `QuicklySearchDiagonal'. */
ChosenPivot = NULL; ChosenPivot = NULL;
@ -1691,7 +1702,7 @@ RealNumber FindBiggestInColExclude();
* ChosenPivot (ElementPtr) * ChosenPivot (ElementPtr)
* Pointer to the element that has been chosen to be the pivot. * Pointer to the element that has been chosen to be the pivot.
* Size (int) * 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) * Magnitude (RealNumber)
* Absolute value of diagonal element. * Absolute value of diagonal element.
* MinMarkowitzProduct (long) * MinMarkowitzProduct (long)
@ -1716,17 +1727,16 @@ RealNumber FindBiggestInColExclude();
static ElementPtr static ElementPtr
SearchDiagonal( SearchDiagonal(
MatrixPtr Matrix, MatrixPtr Matrix,
register int Step int Step
) )
{ {
register int J; int J;
register long MinMarkowitzProduct, *pMarkowitzProduct; long MinMarkowitzProduct, *pMarkowitzProduct;
register int I; int I;
register ElementPtr pDiag; ElementPtr pDiag;
int NumberOfTies = 0, Size = Matrix->Size; int NumberOfTies = 0, Size = Matrix->Size;
ElementPtr ChosenPivot; ElementPtr ChosenPivot;
RealNumber Magnitude, Ratio, RatioOfAccepted = 0, LargestInCol; RealNumber Magnitude, Ratio, RatioOfAccepted = 0, LargestInCol;
RealNumber FindBiggestInColExclude();
/* Begin `SearchDiagonal'. */ /* Begin `SearchDiagonal'. */
ChosenPivot = NULL; ChosenPivot = NULL;
@ -1815,7 +1825,7 @@ RealNumber FindBiggestInColExclude();
* LargestElementMag (RealNumber) * LargestElementMag (RealNumber)
* Magnitude of the largest element yet found in the reduced submatrix. * Magnitude of the largest element yet found in the reduced submatrix.
* Size (int) * Size (int)
* Local version of Size; placed in a register for speed. * Local version of Size; placed in a for speed.
* Magnitude (RealNumber) * Magnitude (RealNumber)
* Absolute value of diagonal element. * Absolute value of diagonal element.
* MinMarkowitzProduct (long) * MinMarkowitzProduct (long)
@ -1848,8 +1858,8 @@ SearchEntireMatrix(
int Step int Step
) )
{ {
register int I, Size = Matrix->Size; int I, Size = Matrix->Size;
register ElementPtr pElement; ElementPtr pElement;
int NumberOfTies = 0; int NumberOfTies = 0;
long Product, MinMarkowitzProduct; long Product, MinMarkowitzProduct;
ElementPtr ChosenPivot, pLargestElement = NULL; ElementPtr ChosenPivot, pLargestElement = NULL;
@ -1966,7 +1976,7 @@ RealNumber FindLargestInCol();
*/ */
static RealNumber static RealNumber
FindLargestInCol( register ElementPtr pElement ) FindLargestInCol( ElementPtr pElement )
{ {
RealNumber Magnitude, Largest = 0.0; RealNumber Magnitude, Largest = 0.0;
@ -2032,11 +2042,11 @@ RealNumber Magnitude, Largest = 0.0;
static RealNumber static RealNumber
FindBiggestInColExclude( FindBiggestInColExclude(
MatrixPtr Matrix, MatrixPtr Matrix,
register ElementPtr pElement, ElementPtr pElement,
register int Step int Step
) )
{ {
register int Row; int Row;
int Col; int Col;
RealNumber Largest, Magnitude; RealNumber Largest, Magnitude;
@ -2110,10 +2120,10 @@ static void
ExchangeRowsAndCols( ExchangeRowsAndCols(
MatrixPtr Matrix, MatrixPtr Matrix,
ElementPtr pPivot, ElementPtr pPivot,
register int Step int Step
) )
{ {
register int Row, Col; int Row, Col;
long OldMarkowitzProd_Step, OldMarkowitzProd_Row, OldMarkowitzProd_Col; long OldMarkowitzProd_Step, OldMarkowitzProd_Row, OldMarkowitzProd_Col;
/* Begin `ExchangeRowsAndCols'. */ /* Begin `ExchangeRowsAndCols'. */
@ -2237,7 +2247,7 @@ spcRowExchange(
int Row2 int Row2
) )
{ {
register ElementPtr Row1Ptr, Row2Ptr; ElementPtr Row1Ptr, Row2Ptr;
int Column; int Column;
ElementPtr Element1, Element2; ElementPtr Element1, Element2;
@ -2340,7 +2350,7 @@ spcColExchange(
int Col2 int Col2
) )
{ {
register ElementPtr Col1Ptr, Col2Ptr; ElementPtr Col1Ptr, Col2Ptr;
int Row; int Row;
ElementPtr Element1, Element2; ElementPtr Element1, Element2;
@ -2445,15 +2455,15 @@ static void
ExchangeColElements( ExchangeColElements(
MatrixPtr Matrix, MatrixPtr Matrix,
int Row1, int Row1,
register ElementPtr Element1, ElementPtr Element1,
int Row2, int Row2,
register ElementPtr Element2, ElementPtr Element2,
int Column int Column
) )
{ {
ElementPtr *ElementAboveRow1, *ElementAboveRow2; ElementPtr *ElementAboveRow1, *ElementAboveRow2;
ElementPtr ElementBelowRow1, ElementBelowRow2; ElementPtr ElementBelowRow1, ElementBelowRow2;
register ElementPtr pElement; ElementPtr pElement;
/* Begin `ExchangeColElements'. */ /* Begin `ExchangeColElements'. */
/* Search to find the ElementAboveRow1. */ /* Search to find the ElementAboveRow1. */
@ -2590,15 +2600,15 @@ static void
ExchangeRowElements( ExchangeRowElements(
MatrixPtr Matrix, MatrixPtr Matrix,
int Col1, int Col1,
register ElementPtr Element1, ElementPtr Element1,
int Col2, int Col2,
register ElementPtr Element2, ElementPtr Element2,
int Row int Row
) )
{ {
ElementPtr *ElementLeftOfCol1, *ElementLeftOfCol2; ElementPtr *ElementLeftOfCol1, *ElementLeftOfCol2;
ElementPtr ElementRightOfCol1, ElementRightOfCol2; ElementPtr ElementRightOfCol1, ElementRightOfCol2;
register ElementPtr pElement; ElementPtr pElement;
/* Begin `ExchangeRowElements'. */ /* Begin `ExchangeRowElements'. */
/* Search to find the ElementLeftOfCol1. */ /* Search to find the ElementLeftOfCol1. */
@ -2727,13 +2737,13 @@ register ElementPtr pElement;
static void static void
RealRowColElimination( RealRowColElimination(
MatrixPtr Matrix, MatrixPtr Matrix,
register ElementPtr pPivot ElementPtr pPivot
) )
{ {
#if REAL #if REAL
register ElementPtr pSub, *ppAbove; ElementPtr pSub, *ppAbove;
register int Row; int Row;
register ElementPtr pLower, pUpper; ElementPtr pLower, pUpper;
/* Begin `RealRowColElimination'. */ /* Begin `RealRowColElimination'. */
@ -2819,13 +2829,13 @@ register ElementPtr pLower, pUpper;
static void static void
ComplexRowColElimination( ComplexRowColElimination(
MatrixPtr Matrix, MatrixPtr Matrix,
register ElementPtr pPivot ElementPtr pPivot
) )
{ {
#if spCOMPLEX #if spCOMPLEX
register ElementPtr pSub, *ppAbove; ElementPtr pSub, *ppAbove;
register int Row; int Row;
register ElementPtr pLower, pUpper; ElementPtr pLower, pUpper;
/* Begin `ComplexRowColElimination'. */ /* Begin `ComplexRowColElimination'. */
@ -2845,7 +2855,7 @@ register ElementPtr pLower, pUpper;
pSub = pUpper->NextInCol; pSub = pUpper->NextInCol;
pLower = pPivot->NextInCol; pLower = pPivot->NextInCol;
ppAbove = &pUpper->NextInCol; ppAbove = &pUpper->NextInCol;
while (pLower != NULL) while (pLower != NULL)
{ Row = pLower->Row; { Row = pLower->Row;
@ -2909,9 +2919,9 @@ UpdateMarkowitzNumbers(
ElementPtr pPivot ElementPtr pPivot
) )
{ {
register int Row, Col; int Row, Col;
register ElementPtr ColPtr, RowPtr; ElementPtr ColPtr, RowPtr;
register int *MarkoRow = Matrix->MarkowitzRow, *MarkoCol = Matrix->MarkowitzCol; int *MarkoRow = Matrix->MarkowitzRow, *MarkoCol = Matrix->MarkowitzCol;
double Product; double Product;
/* Begin `UpdateMarkowitzNumbers'. */ /* Begin `UpdateMarkowitzNumbers'. */
@ -3074,8 +3084,6 @@ int I;
printf("%2d ", Matrix->ExtToIntColMap[I]); printf("%2d ", Matrix->ExtToIntColMap[I]);
printf("\n\n"); printf("\n\n");
/* spPrint((char *)Matrix, NO, YES); */
return; return;
} }

View File

@ -142,7 +142,7 @@ spPrint(
int Header int Header
) )
{ {
register int J = 0; int J = 0;
int I, Row, Col, Size, Top, StartCol = 1, StopCol, Columns, ElementCount = 0; int I, Row, Col, Size, Top, StartCol = 1, StopCol, Columns, ElementCount = 0;
double Magnitude, SmallestDiag = 0.0, SmallestElement = 0.0; double Magnitude, SmallestDiag = 0.0, SmallestElement = 0.0;
double LargestElement = 0.0, LargestDiag = 0.0; double LargestElement = 0.0, LargestDiag = 0.0;
@ -418,8 +418,8 @@ spFileMatrix(
int Header int Header
) )
{ {
register int I, Size; int I, Size;
register ElementPtr pElement; ElementPtr pElement;
int Row, Col, Err; int Row, Col, Err;
FILE *pMatrixFile; FILE *pMatrixFile;
@ -571,9 +571,9 @@ spFileVector(
#endif #endif
) )
{ {
register int I, Size; int I, Size;
#if spCOMPLEX #if spCOMPLEX
register int Err; int Err;
#endif #endif
FILE *pMatrixFile; FILE *pMatrixFile;
@ -581,9 +581,12 @@ FILE *pMatrixFile;
ASSERT_IS_SPARSE( Matrix ); ASSERT_IS_SPARSE( Matrix );
vASSERT( RHS != NULL, "Vector missing" ); vASSERT( RHS != NULL, "Vector missing" );
/* Open File in append mode. */ if (File) {
if ((pMatrixFile = fopen(File,"a")) == NULL) /* Open File in write mode. */
return 0; pMatrixFile = fopen(File,"w");
if (pMatrixFile == NULL)
return 0;
}
/* Correct array pointers for ARRAY_OFFSET. */ /* Correct array pointers for ARRAY_OFFSET. */
#if NOT ARRAY_OFFSET #if NOT ARRAY_OFFSET
@ -694,8 +697,8 @@ spFileStats(
char *Label char *Label
) )
{ {
register int Size, I; int Size, I;
register ElementPtr pElement; ElementPtr pElement;
int NumberOfElements; int NumberOfElements;
RealNumber Data, LargestElement, SmallestElement; RealNumber Data, LargestElement, SmallestElement;
FILE *pStatsFile; FILE *pStatsFile;

View File

@ -285,12 +285,11 @@ SMPmatrix *Matrix)
int int
SMPnewMatrix( SMPnewMatrix(
SMPmatrix **pMatrix, SMPmatrix **pMatrix,
int dummy) int size)
{ {
int Error; int Error;
NG_IGNORE(dummy);
*pMatrix = (SMPmatrix *)spCreate( 0, 1, &Error ); *pMatrix = (SMPmatrix *)spCreate( size, 1, &Error );
return Error; return Error;
} }
@ -562,18 +561,18 @@ SMPcAddCol(SMPmatrix *Matrix, int Accum_Col, int Addend_Col)
Accum = *Prev; Accum = *Prev;
while (Addend != NULL) { while (Addend != NULL) {
while (Accum && Accum->Row < Addend->Row) { while (Accum && Accum->Row < Addend->Row) {
Prev = &Accum->NextInCol; Prev = &Accum->NextInCol;
Accum = *Prev; Accum = *Prev;
} }
if (!Accum || Accum->Row > Addend->Row) { if (!Accum || Accum->Row > Addend->Row) {
Accum = spcCreateElement(Matrix, Addend->Row, Accum_Col, Prev, 0, 0); Accum = spcCreateElement(Matrix, Addend->Row, Accum_Col, Prev, 0, NO);
} }
Accum->Real += Addend->Real; Accum->Real += Addend->Real;
#if spCOMPLEX #if spCOMPLEX
Accum->Imag += Addend->Imag; Accum->Imag += Addend->Imag;
#endif #endif
Addend = Addend->NextInCol; Addend = Addend->NextInCol;
} }
return spErrorState( Matrix ); return spErrorState( Matrix );

View File

@ -156,13 +156,12 @@ spSolve(
) )
{ {
#if REAL #if REAL
register ElementPtr pElement; ElementPtr pElement;
register RealVector Intermediate; RealVector Intermediate;
register RealNumber Temp; RealNumber Temp;
register int I, *pExtOrder, Size; int I, *pExtOrder, Size;
ElementPtr pPivot; ElementPtr pPivot;
#endif #endif
void SolveComplexMatrix();
/* Begin `spSolve'. */ /* Begin `spSolve'. */
ASSERT_IS_SPARSE( Matrix ); ASSERT_IS_SPARSE( Matrix );
@ -301,12 +300,12 @@ SolveComplexMatrix(
# endif # endif
) )
{ {
register ElementPtr pElement; ElementPtr pElement;
register ComplexVector Intermediate; ComplexVector Intermediate;
register int I, *pExtOrder, Size; int I, *pExtOrder, Size;
ElementPtr pPivot; ElementPtr pPivot;
#if NOT spSEPARATED_COMPLEX_VECTORS #if NOT spSEPARATED_COMPLEX_VECTORS
register ComplexVector ExtVector; ComplexVector ExtVector;
#endif #endif
ComplexNumber Temp; ComplexNumber Temp;
@ -471,13 +470,12 @@ spSolveTransposed(
) )
{ {
#if REAL #if REAL
register ElementPtr pElement; ElementPtr pElement;
register RealVector Intermediate; RealVector Intermediate;
register int I, *pExtOrder, Size; int I, *pExtOrder, Size;
ElementPtr pPivot; ElementPtr pPivot;
RealNumber Temp; RealNumber Temp;
#endif #endif
void SolveComplexTransposedMatrix();
/* Begin `spSolveTransposed'. */ /* Begin `spSolveTransposed'. */
ASSERT_IS_SPARSE( Matrix ); ASSERT_IS_SPARSE( Matrix );
@ -618,11 +616,11 @@ SolveComplexTransposedMatrix(
# endif # endif
) )
{ {
register ElementPtr pElement; ElementPtr pElement;
register ComplexVector Intermediate; ComplexVector Intermediate;
register int I, *pExtOrder, Size; int I, *pExtOrder, Size;
#if NOT spSEPARATED_COMPLEX_VECTORS #if NOT spSEPARATED_COMPLEX_VECTORS
register ComplexVector ExtVector; ComplexVector ExtVector;
#endif #endif
ElementPtr pPivot; ElementPtr pPivot;
ComplexNumber Temp; ComplexNumber Temp;

View File

@ -203,7 +203,7 @@ static RealNumber ComplexCondition( MatrixPtr, RealNumber, int* );
void void
spMNA_Preorder( MatrixPtr Matrix ) spMNA_Preorder( MatrixPtr Matrix )
{ {
register int J, Size; int J, Size;
ElementPtr pTwin1, pTwin2; ElementPtr pTwin1, pTwin2;
int Twins, StartAt = 1; int Twins, StartAt = 1;
BOOLEAN Swapped, AnotherPassNeeded; BOOLEAN Swapped, AnotherPassNeeded;
@ -397,10 +397,9 @@ spScale(
spREAL SolutionScaleFactors[] spREAL SolutionScaleFactors[]
) )
{ {
register ElementPtr pElement; ElementPtr pElement;
register int I, lSize, *pExtOrder; int I, lSize, *pExtOrder;
RealNumber ScaleFactor; RealNumber ScaleFactor;
//void ScaleComplexMatrix();
/* Begin `spScale'. */ /* Begin `spScale'. */
ASSERT_IS_SPARSE( Matrix ); ASSERT_IS_SPARSE( Matrix );
@ -523,12 +522,12 @@ RealNumber ScaleFactor;
static void static void
ScaleComplexMatrix( ScaleComplexMatrix(
MatrixPtr Matrix, MatrixPtr Matrix,
register RealVector RHS_ScaleFactors, RealVector RHS_ScaleFactors,
register RealVector SolutionScaleFactors RealVector SolutionScaleFactors
) )
{ {
register ElementPtr pElement; ElementPtr pElement;
register int I, lSize, *pExtOrder; int I, lSize, *pExtOrder;
RealNumber ScaleFactor; RealNumber ScaleFactor;
/* Begin `ScaleComplexMatrix'. */ /* Begin `ScaleComplexMatrix'. */
@ -613,12 +612,11 @@ spMultiply(
) )
{ {
#if REAL #if REAL
register ElementPtr pElement; ElementPtr pElement;
register RealVector Vector; RealVector Vector;
register RealNumber Sum; RealNumber Sum;
register int I, *pExtOrder; int I, *pExtOrder;
#endif #endif
extern void ComplexMatrixMultiply();
/* Begin `spMultiply'. */ /* Begin `spMultiply'. */
ASSERT_IS_SPARSE( Matrix ); ASSERT_IS_SPARSE( Matrix );
@ -710,10 +708,10 @@ ComplexMatrixMultiply(
#endif #endif
) )
{ {
register ElementPtr pElement; ElementPtr pElement;
register ComplexVector Vector; ComplexVector Vector;
ComplexNumber Sum; ComplexNumber Sum;
register int I, *pExtOrder; int I, *pExtOrder;
/* Begin `ComplexMatrixMultiply'. */ /* Begin `ComplexMatrixMultiply'. */
@ -805,11 +803,10 @@ spMultTransposed(
) )
{ {
#if REAL #if REAL
register ElementPtr pElement; ElementPtr pElement;
register RealVector Vector; RealVector Vector;
register RealNumber Sum; RealNumber Sum;
register int I, *pExtOrder; int I, *pExtOrder;
extern void ComplexTransposedMatrixMultiply();
#endif #endif
/* Begin `spMultTransposed'. */ /* Begin `spMultTransposed'. */
@ -906,10 +903,10 @@ ComplexTransposedMatrixMultiply(
#endif #endif
) )
{ {
register ElementPtr pElement; ElementPtr pElement;
register ComplexVector Vector; ComplexVector Vector;
ComplexNumber Sum; ComplexNumber Sum;
register int I, *pExtOrder; int I, *pExtOrder;
/* Begin `ComplexTransposedMatrixMultiply'. */ /* Begin `ComplexTransposedMatrixMultiply'. */
@ -999,7 +996,7 @@ register int I, *pExtOrder;
* Norm (RealNumber) * Norm (RealNumber)
* L-infinity norm of a complex number. * L-infinity norm of a complex number.
* Size (int) * 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) * Temp (RealNumber)
* Temporary storage for real portion of determinant. * Temporary storage for real portion of determinant.
*/ */
@ -1014,7 +1011,7 @@ spDeterminant(
#endif #endif
) )
{ {
register int I, Size; int I, Size;
RealNumber Norm, nr, ni; RealNumber Norm, nr, ni;
ComplexNumber Pivot, cDeterminant; ComplexNumber Pivot, cDeterminant;
@ -1170,7 +1167,7 @@ struct FillinListNodeStruct *pListNode;
Matrix->Fillins = 0; Matrix->Fillins = 0;
/* Mark the fill-ins. */ /* Mark the fill-ins. */
{ register ElementPtr pFillin, pLastFillin; { ElementPtr pFillin, pLastFillin;
pListNode = Matrix->LastFillinListNode = Matrix->FirstFillinListNode; pListNode = Matrix->LastFillinListNode = Matrix->FirstFillinListNode;
Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList; Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList;
@ -1186,8 +1183,8 @@ struct FillinListNodeStruct *pListNode;
} }
/* Unlink fill-ins by searching for elements marked with Row = 0. */ /* Unlink fill-ins by searching for elements marked with Row = 0. */
{ register ElementPtr pElement, *ppElement; { ElementPtr pElement, *ppElement;
register int I, Size = Matrix->Size; int I, Size = Matrix->Size;
/* Unlink fill-ins in all columns. */ /* Unlink fill-ins in all columns. */
for (I = 1; I <= Size; I++) for (I = 1; I <= Size; I++)
@ -1263,7 +1260,7 @@ spDeleteRowAndCol(
int Col int Col
) )
{ {
register ElementPtr pElement, *ppElement, pLastElement; ElementPtr pElement, *ppElement, pLastElement;
int Size, ExtRow, ExtCol; int Size, ExtRow, ExtCol;
/* Begin `spDeleteRowAndCol'. */ /* Begin `spDeleteRowAndCol'. */
@ -1370,8 +1367,8 @@ int Size, ExtRow, ExtCol;
spREAL spREAL
spPseudoCondition( MatrixPtr Matrix ) spPseudoCondition( MatrixPtr Matrix )
{ {
register int I; int I;
register ArrayOfElementPtrs Diag; ArrayOfElementPtrs Diag;
RealNumber MaxPivot, MinPivot, Mag; RealNumber MaxPivot, MinPivot, Mag;
/* Begin `spPseudoCondition'. */ /* Begin `spPseudoCondition'. */
@ -1458,9 +1455,9 @@ spCondition(
int *pError int *pError
) )
{ {
register ElementPtr pElement; ElementPtr pElement;
register RealVector T, Tm; RealVector T, Tm;
register int I, K, Row; int I, K, Row;
ElementPtr pPivot; ElementPtr pPivot;
int Size; int Size;
RealNumber E, Em, Wp, Wm, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor; RealNumber E, Em, Wp, Wm, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor;
@ -1626,7 +1623,7 @@ RealNumber Linpack, OLeary, InvNormOfInverse, ComplexCondition();
for (ASz = 0.0, I = Size; I > 0; I--) ASz += ABS(T[I]); for (ASz = 0.0, I = Size; I > 0; I--) ASz += ABS(T[I]);
#if NOT spCOMPLEX #if NOT spCOMPLEX
FREE( Tm ); SP_FREE( Tm );
#endif #endif
Linpack = ASy / ASz; Linpack = ASy / ASz;
@ -1669,9 +1666,9 @@ ComplexCondition(
int *pError int *pError
) )
{ {
register ElementPtr pElement; ElementPtr pElement;
register ComplexVector T, Tm; ComplexVector T, Tm;
register int I, K, Row; int I, K, Row;
ElementPtr pPivot; ElementPtr pPivot;
int Size; int Size;
RealNumber E, Em, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor; RealNumber E, Em, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor;
@ -1827,7 +1824,7 @@ ComplexNumber Wp, Wm;
/* 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]); for (ASz = 0.0, I = Size; I > 0; I--) ASz += CMPLX_1_NORM(T[I]);
FREE( Tm ); SP_FREE( Tm );
Linpack = ASy / ASz; Linpack = ASy / ASz;
OLeary = E / MaxY; OLeary = E / MaxY;
@ -1854,8 +1851,8 @@ ComplexNumber Wp, Wm;
spREAL spREAL
spNorm( MatrixPtr Matrix ) spNorm( MatrixPtr Matrix )
{ {
register ElementPtr pElement; ElementPtr pElement;
register int I; int I;
RealNumber Max = 0.0, AbsRowSum; RealNumber Max = 0.0, AbsRowSum;
/* Begin `spNorm'. */ /* Begin `spNorm'. */
@ -1968,11 +1965,11 @@ RealNumber Max = 0.0, AbsRowSum;
spREAL spREAL
spLargestElement( MatrixPtr Matrix ) spLargestElement( MatrixPtr Matrix )
{ {
register int I; int I;
RealNumber Mag, AbsColSum, Max = 0.0, MaxRow = 0.0, MaxCol = 0.0; RealNumber Mag, AbsColSum, Max = 0.0, MaxRow = 0.0, MaxCol = 0.0;
RealNumber Pivot; RealNumber Pivot;
ComplexNumber cPivot; ComplexNumber cPivot;
register ElementPtr pElement, pDiag; ElementPtr pElement, pDiag;
/* Begin `spLargestElement'. */ /* Begin `spLargestElement'. */
ASSERT_IS_SPARSE( Matrix ); ASSERT_IS_SPARSE( Matrix );
@ -2089,8 +2086,8 @@ spRoundoff(
spREAL Rho spREAL Rho
) )
{ {
register ElementPtr pElement; ElementPtr pElement;
register int Count, I, MaxCount = 0; int Count, I, MaxCount = 0;
RealNumber Reid, Gear; RealNumber Reid, Gear;
/* Begin `spRoundoff'. */ /* Begin `spRoundoff'. */