From 4a3b6294391b58b0e5f7bd6d38c8c75ca7efecb1 Mon Sep 17 00:00:00 2001 From: dwarning Date: Fri, 16 Jun 2023 15:01:19 +0200 Subject: [PATCH] more adaption to existing solution --- src/maths/sparse/spAllocate.c | 30 +++---- src/maths/sparse/spBuild.c | 42 +++++----- src/maths/sparse/spFactor.c | 154 ++++++++++++++++++---------------- src/maths/sparse/spOutput.c | 23 ++--- src/maths/sparse/spSMP.c | 27 +++--- src/maths/sparse/spSolve.c | 32 ++++--- src/maths/sparse/spUtils.c | 89 ++++++++++---------- 7 files changed, 201 insertions(+), 196 deletions(-) diff --git a/src/maths/sparse/spAllocate.c b/src/maths/sparse/spAllocate.c index b8afd886f..88d7a7b75 100644 --- a/src/maths/sparse/spAllocate.c +++ b/src/maths/sparse/spAllocate.c @@ -133,9 +133,9 @@ spCreate( spError *pError ) { -register unsigned SizePlusOne; -register MatrixPtr Matrix; -register int I; +unsigned SizePlusOne; +MatrixPtr Matrix; +int I; int AllocatedSize; /* Begin `spCreate'. */ @@ -201,7 +201,7 @@ int AllocatedSize; Matrix->ElementsRemaining = 0; Matrix->FillinsRemaining = 0; - RecordAllocation( Matrix, (void *)Matrix ); + RecordAllocation( Matrix, Matrix ); if (Matrix->Error == spNO_MEMORY) goto MemoryError; /* FIXME: Use of memory after free */ /* Take out the trash. */ @@ -322,7 +322,7 @@ ElementPtr pElement; /* Allocate block of MatrixElements if necessary. */ if (Matrix->ElementsRemaining == 0) { pElement = SP_MALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION); - RecordAllocation( Matrix, (void *)pElement ); + RecordAllocation( Matrix, pElement ); if (Matrix->Error == spNO_MEMORY) return NULL; Matrix->ElementsRemaining = ELEMENTS_PER_ALLOCATION; Matrix->NextAvailElement = pElement; @@ -383,21 +383,21 @@ ElementPtr pElement; /* Allocate block of MatrixElements for elements. */ pElement = SP_MALLOC(struct MatrixElement, InitialNumberOfElements); - RecordAllocation( Matrix, (void *)pElement ); + RecordAllocation( Matrix, pElement ); if (Matrix->Error == spNO_MEMORY) return; Matrix->ElementsRemaining = InitialNumberOfElements; Matrix->NextAvailElement = pElement; /* Allocate block of MatrixElements for fill-ins. */ pElement = SP_MALLOC(struct MatrixElement, NumberOfFillinsExpected); - RecordAllocation( Matrix, (void *)pElement ); + RecordAllocation( Matrix, pElement ); if (Matrix->Error == spNO_MEMORY) return; Matrix->FillinsRemaining = NumberOfFillinsExpected; Matrix->NextAvailFillin = pElement; /* Allocate a fill-in list structure. */ Matrix->FirstFillinListNode = SP_MALLOC(struct FillinListNodeStruct,1); - RecordAllocation( Matrix, (void *)Matrix->FirstFillinListNode ); + RecordAllocation( Matrix, Matrix->FirstFillinListNode ); if (Matrix->Error == spNO_MEMORY) return; Matrix->LastFillinListNode = Matrix->FirstFillinListNode; @@ -465,14 +465,14 @@ ElementPtr pFillins; { /* Allocate block of fill-ins. */ pFillins = SP_MALLOC(struct MatrixElement, ELEMENTS_PER_ALLOCATION); - RecordAllocation( Matrix, (void *)pFillins ); + RecordAllocation( Matrix, pFillins ); if (Matrix->Error == spNO_MEMORY) return NULL; Matrix->FillinsRemaining = ELEMENTS_PER_ALLOCATION; Matrix->NextAvailFillin = pFillins; /* Allocate a fill-in list structure. */ pListNode->Next = SP_MALLOC(struct FillinListNodeStruct,1); - RecordAllocation( Matrix, (void *)pListNode->Next ); + RecordAllocation( Matrix, pListNode->Next ); if (Matrix->Error == spNO_MEMORY) return NULL; Matrix->LastFillinListNode = pListNode = pListNode->Next; @@ -505,7 +505,7 @@ ElementPtr pFillins; * >>> Arguments: * Matrix (MatrixPtr) * Pointer to the matrix. - * AllocatedPtr (void *) + * AllocatedPtr * The pointer returned by malloc or calloc. These pointers are saved in * a list so that they can be easily freed. * @@ -573,8 +573,8 @@ RecordAllocation( static void AllocateBlockOfAllocationList( MatrixPtr Matrix ) { -register int I; -register AllocationListPtr ListPtr; +int I; +AllocationListPtr ListPtr; /* Begin `AllocateBlockOfAllocationList'. */ /* Allocate block of records for allocation list. */ @@ -596,7 +596,7 @@ register AllocationListPtr ListPtr; } /* Record allocation of space for allocation list on allocation list. */ - Matrix->TopOfAllocationList->AllocatedPtr = (void *)ListPtr; + Matrix->TopOfAllocationList->AllocatedPtr = ListPtr; Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION; return; @@ -629,7 +629,7 @@ register AllocationListPtr ListPtr; void spDestroy( MatrixPtr Matrix ) { -register AllocationListPtr ListPtr, NextListPtr; +AllocationListPtr ListPtr, NextListPtr; /* Begin `spDestroy'. */ ASSERT_IS_SPARSE( Matrix ); diff --git a/src/maths/sparse/spBuild.c b/src/maths/sparse/spBuild.c index 2ce9471be..ff767d0f2 100644 --- a/src/maths/sparse/spBuild.c +++ b/src/maths/sparse/spBuild.c @@ -104,8 +104,8 @@ static void EnlargeMatrix( MatrixPtr, int ); void spClear( MatrixPtr Matrix ) { -register ElementPtr pElement; -register int I; +ElementPtr pElement; +int I; /* Begin `spClear'. */ ASSERT_IS_SPARSE( Matrix ); @@ -184,7 +184,7 @@ spFindElement( int Col ) { -register ElementPtr pElement; +ElementPtr pElement; int StartAt = 0; long int Min = LARGEST_LONG_INTEGER; #define BorderRight 0 /* Start at left border, move right. */ @@ -449,7 +449,7 @@ Translate( int *Col ) { -register int IntRow, IntCol, ExtRow, ExtCol; +int IntRow, IntCol, ExtRow, ExtCol; /* Begin `Translate'. */ ExtRow = *Row; @@ -649,7 +649,7 @@ spGetQuad( OR (Template->Element4Negated == NULL) ) return spNO_MEMORY; - if (Template->Element1 == &((MatrixPtr)Matrix)->TrashCan.Real) + if (Template->Element1 == &(Matrix)->TrashCan.Real) SWAP( RealNumber *, Template->Element1, Template->Element2 ); return spOKAY; @@ -756,10 +756,10 @@ spGetOnes( ElementPtr spcFindDiag( MatrixPtr Matrix, - register int Index + int Index ) { -register ElementPtr pElement; +ElementPtr pElement; /* Begin `spcFindDiag'. */ pElement = Matrix->FirstInCol[Index]; @@ -823,18 +823,18 @@ ElementPtr spcCreateElement( MatrixPtr Matrix, int Row, - register int Col, - register ElementPtr *ppToLeft, - register ElementPtr *ppAbove, + int Col, + ElementPtr *ppToLeft, + ElementPtr *ppAbove, BOOLEAN Fillin ) { -register ElementPtr pElement, pCreatedElement; +ElementPtr pElement, pCreatedElement; /* Begin `spcCreateElement'. */ /* Find element immediately above the desired element. */ - pElement = *ppAbove; + pElement = *ppAbove; /* FIXME: Dereference of null pointer */ while ((pElement != NULL) AND (pElement->Row < Row)) { ppAbove = &pElement->NextInCol; pElement = *ppAbove; @@ -937,7 +937,7 @@ register ElementPtr pElement, pCreatedElement; * 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. @@ -946,9 +946,9 @@ register ElementPtr pElement, pCreatedElement; void spcLinkRows( MatrixPtr Matrix ) { -register ElementPtr pElement, *FirstInRowEntry; -register ArrayOfElementPtrs FirstInRowArray; -register int Col; +ElementPtr pElement, *FirstInRowEntry; +ArrayOfElementPtrs FirstInRowArray; +int Col; /* Begin `spcLinkRows'. */ FirstInRowArray = Matrix->FirstInRow; @@ -998,10 +998,10 @@ register int Col; static void EnlargeMatrix( MatrixPtr Matrix, - register int NewSize + int NewSize ) { -register int I, OldAllocatedSize = Matrix->AllocatedSize; +int I, OldAllocatedSize = Matrix->AllocatedSize; /* Begin `EnlargeMatrix'. */ Matrix->Size = NewSize; @@ -1087,10 +1087,10 @@ register int I, OldAllocatedSize = Matrix->AllocatedSize; static void ExpandTranslationArrays( MatrixPtr Matrix, - register int NewSize + int NewSize ) { -register int I, OldAllocatedSize = Matrix->AllocatedExtSize; +int I, OldAllocatedSize = Matrix->AllocatedExtSize; /* Begin `ExpandTranslationArrays'. */ Matrix->ExtSize = NewSize; @@ -1176,7 +1176,7 @@ spInitialize( ) ) { -register ElementPtr pElement; +ElementPtr pElement; int J, Error, Col; /* Begin `spInitialize'. */ diff --git a/src/maths/sparse/spFactor.c b/src/maths/sparse/spFactor.c index 6ed0add0f..a91de02c8 100644 --- a/src/maths/sparse/spFactor.c +++ b/src/maths/sparse/spFactor.c @@ -228,6 +228,11 @@ RealNumber LargestInCol, FindLargestInCol(); /* Matrix has been factored before and reordering is not required. */ for (Step = 1; Step <= Size; 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); if ((LargestInCol * RelThreshold < ELEMENT_MAG(pPivot))) { if (Matrix->Complex) @@ -337,9 +342,9 @@ spError spFactor( MatrixPtr Matrix ) { #if REAL -register ElementPtr pElement; -register ElementPtr pColumn; -register int Step, Size; +ElementPtr pElement; +ElementPtr pColumn; +int Step, Size; RealNumber Mult; #endif /* Begin `spFactor'. */ @@ -348,7 +353,7 @@ RealNumber Mult; ASSERT_IS_NOT_FACTORED( Matrix ); if (Matrix->NeedsOrdering) - { return spOrderAndFactor( Matrix, (RealVector)NULL, + { return spOrderAndFactor( Matrix, NULL, 0.0, 0.0, DIAG_PIVOTING_AS_DEFAULT ); } if (NOT Matrix->Partitioned) spPartition( Matrix, spDEFAULT_PARTITION ); @@ -359,6 +364,11 @@ RealNumber Mult; #if REAL Size = Matrix->Size; + if (Size == 0) { + Matrix->Factored = YES; + return (Matrix->Error = spOKAY); + } + if (Matrix->Diag[1]->Real == 0.0) return ZeroPivot( Matrix, 1 ); Matrix->Diag[1]->Real = 1.0 / Matrix->Diag[1]->Real; @@ -366,7 +376,7 @@ RealNumber Mult; for (Step = 2; Step <= Size; Step++) { if (Matrix->DoRealDirect[Step]) { /* Update column using direct addressing scatter-gather. */ - register RealNumber *Dest = (RealNumber *)Matrix->Intermediate; + RealNumber *Dest = (RealNumber *)Matrix->Intermediate; /* Scatter. */ pElement = Matrix->FirstInCol[Step]; @@ -398,7 +408,7 @@ RealNumber Mult; } else { /* Update column using indirect addressing scatter-gather. */ - register RealNumber **pDest = (RealNumber **)Matrix->Intermediate; + RealNumber **pDest = (RealNumber **)Matrix->Intermediate; /* Scatter. */ pElement = Matrix->FirstInCol[Step]; @@ -456,15 +466,21 @@ RealNumber Mult; static int FactorComplexMatrix( MatrixPtr Matrix ) { -register ElementPtr pElement; -register ElementPtr pColumn; -register int Step, Size; +ElementPtr pElement; +ElementPtr pColumn; +int Step, Size; ComplexNumber Mult, Pivot; /* Begin `FactorComplexMatrix'. */ ASSERT(Matrix->Complex); Size = Matrix->Size; + + if (Size == 0) { + Matrix->Factored = YES; + return (Matrix->Error = spOKAY); + } + pElement = Matrix->Diag[1]; if (ELEMENT_MAG(pElement) == 0.0) return ZeroPivot( Matrix, 1 ); /* Cmplx expr: *pPivot = 1.0 / *pPivot. */ @@ -474,7 +490,7 @@ ComplexNumber Mult, Pivot; for (Step = 2; Step <= Size; Step++) { if (Matrix->DoCmplxDirect[Step]) { /* Update column using direct addressing scatter-gather. */ - register ComplexNumber *Dest; + ComplexNumber *Dest; Dest = (ComplexNumber *)Matrix->Intermediate; /* Scatter. */ @@ -512,7 +528,7 @@ ComplexNumber Mult, Pivot; } else { /* Update column using direct addressing scatter-gather. */ - register ComplexNumber **pDest; + ComplexNumber **pDest; pDest = (ComplexNumber **)Matrix->Intermediate; /* Scatter. */ @@ -593,10 +609,10 @@ spPartition( int Mode ) { -register ElementPtr pElement, pColumn; -register int Step, Size; -register int *Nc, *No; -register long *Nm; +ElementPtr pElement, pColumn; +int Step, Size; +int *Nc, *No; +long *Nm; #if spCOMPLEX BOOLEAN *DoCmplxDirect; #endif @@ -831,12 +847,12 @@ int Size; static void CountMarkowitz( MatrixPtr Matrix, - register RealVector RHS, + RealVector RHS, int Step ) { -register int Count, I, Size = Matrix->Size; -register ElementPtr pElement; +int Count, I, Size = Matrix->Size; +ElementPtr pElement; int ExtRow; /* Begin `CountMarkowitz'. */ @@ -944,9 +960,9 @@ MarkowitzProducts( int Step ) { -register int I, *pMarkowitzRow, *pMarkowitzCol; -register long Product, *pMarkowitzProduct; -register int Size = Matrix->Size; +int I, *pMarkowitzRow, *pMarkowitzCol; +long Product, *pMarkowitzProduct; +int Size = Matrix->Size; double fProduct; /* Begin `MarkowitzProducts'. */ @@ -1032,11 +1048,7 @@ SearchForPivot( BOOLEAN DiagPivoting ) { -register ElementPtr ChosenPivot; -ElementPtr SearchForSingleton(); -ElementPtr QuicklySearchDiagonal(); -ElementPtr SearchDiagonal(); -ElementPtr SearchEntireMatrix(); +ElementPtr ChosenPivot; /* Begin `SearchForPivot'. */ @@ -1133,11 +1145,11 @@ SearchForSingleton( int Step ) { -register ElementPtr ChosenPivot; -register int I; -register long *pMarkowitzProduct; +ElementPtr ChosenPivot; +int I; +long *pMarkowitzProduct; int Singletons; -RealNumber PivotMag, FindBiggestInColExclude(); +RealNumber PivotMag; /* Begin `SearchForSingleton'. */ /* Initialize pointer that is to scan through MarkowitzProduct vector. */ @@ -1355,8 +1367,8 @@ QuicklySearchDiagonal( 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; @@ -1554,12 +1566,11 @@ QuicklySearchDiagonal( 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'. */ ChosenPivot = NULL; @@ -1691,7 +1702,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) @@ -1716,17 +1727,16 @@ RealNumber FindBiggestInColExclude(); static ElementPtr SearchDiagonal( MatrixPtr Matrix, - register int Step + int Step ) { -register int J; -register long MinMarkowitzProduct, *pMarkowitzProduct; -register int I; -register ElementPtr pDiag; +int J; +long MinMarkowitzProduct, *pMarkowitzProduct; +int I; +ElementPtr pDiag; int NumberOfTies = 0, Size = Matrix->Size; ElementPtr ChosenPivot; RealNumber Magnitude, Ratio, RatioOfAccepted = 0, LargestInCol; -RealNumber FindBiggestInColExclude(); /* Begin `SearchDiagonal'. */ ChosenPivot = NULL; @@ -1815,7 +1825,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) @@ -1848,8 +1858,8 @@ SearchEntireMatrix( int Step ) { -register int I, Size = Matrix->Size; -register ElementPtr pElement; +int I, Size = Matrix->Size; +ElementPtr pElement; int NumberOfTies = 0; long Product, MinMarkowitzProduct; ElementPtr ChosenPivot, pLargestElement = NULL; @@ -1966,7 +1976,7 @@ RealNumber FindLargestInCol(); */ static RealNumber -FindLargestInCol( register ElementPtr pElement ) +FindLargestInCol( ElementPtr pElement ) { RealNumber Magnitude, Largest = 0.0; @@ -2032,11 +2042,11 @@ RealNumber Magnitude, Largest = 0.0; static RealNumber FindBiggestInColExclude( MatrixPtr Matrix, - register ElementPtr pElement, - register int Step + ElementPtr pElement, + int Step ) { -register int Row; +int Row; int Col; RealNumber Largest, Magnitude; @@ -2110,10 +2120,10 @@ static void ExchangeRowsAndCols( MatrixPtr Matrix, ElementPtr pPivot, - register int Step + int Step ) { -register int Row, Col; +int Row, Col; long OldMarkowitzProd_Step, OldMarkowitzProd_Row, OldMarkowitzProd_Col; /* Begin `ExchangeRowsAndCols'. */ @@ -2237,7 +2247,7 @@ spcRowExchange( int Row2 ) { -register ElementPtr Row1Ptr, Row2Ptr; +ElementPtr Row1Ptr, Row2Ptr; int Column; ElementPtr Element1, Element2; @@ -2340,7 +2350,7 @@ spcColExchange( int Col2 ) { -register ElementPtr Col1Ptr, Col2Ptr; +ElementPtr Col1Ptr, Col2Ptr; int Row; ElementPtr Element1, Element2; @@ -2445,15 +2455,15 @@ static void ExchangeColElements( MatrixPtr Matrix, int Row1, - register ElementPtr Element1, + ElementPtr Element1, int Row2, - register ElementPtr Element2, + ElementPtr Element2, int Column ) { ElementPtr *ElementAboveRow1, *ElementAboveRow2; ElementPtr ElementBelowRow1, ElementBelowRow2; -register ElementPtr pElement; +ElementPtr pElement; /* Begin `ExchangeColElements'. */ /* Search to find the ElementAboveRow1. */ @@ -2590,15 +2600,15 @@ static void ExchangeRowElements( MatrixPtr Matrix, int Col1, - register ElementPtr Element1, + ElementPtr Element1, int Col2, - register ElementPtr Element2, + ElementPtr Element2, int Row ) { ElementPtr *ElementLeftOfCol1, *ElementLeftOfCol2; ElementPtr ElementRightOfCol1, ElementRightOfCol2; -register ElementPtr pElement; + ElementPtr pElement; /* Begin `ExchangeRowElements'. */ /* Search to find the ElementLeftOfCol1. */ @@ -2727,13 +2737,13 @@ register ElementPtr pElement; static void RealRowColElimination( MatrixPtr Matrix, - register ElementPtr pPivot + ElementPtr pPivot ) { #if REAL -register ElementPtr pSub, *ppAbove; -register int Row; -register ElementPtr pLower, pUpper; +ElementPtr pSub, *ppAbove; +int Row; +ElementPtr pLower, pUpper; /* Begin `RealRowColElimination'. */ @@ -2819,13 +2829,13 @@ register ElementPtr pLower, pUpper; static void ComplexRowColElimination( MatrixPtr Matrix, - register ElementPtr pPivot + ElementPtr pPivot ) { #if spCOMPLEX -register ElementPtr pSub, *ppAbove; -register int Row; -register ElementPtr pLower, pUpper; +ElementPtr pSub, *ppAbove; +int Row; +ElementPtr pLower, pUpper; /* Begin `ComplexRowColElimination'. */ @@ -2845,7 +2855,7 @@ register ElementPtr pLower, pUpper; pSub = pUpper->NextInCol; pLower = pPivot->NextInCol; - ppAbove = &pUpper->NextInCol; + ppAbove = &pUpper->NextInCol; while (pLower != NULL) { Row = pLower->Row; @@ -2909,9 +2919,9 @@ UpdateMarkowitzNumbers( ElementPtr pPivot ) { -register int Row, Col; -register ElementPtr ColPtr, RowPtr; -register int *MarkoRow = Matrix->MarkowitzRow, *MarkoCol = Matrix->MarkowitzCol; +int Row, Col; +ElementPtr ColPtr, RowPtr; +int *MarkoRow = Matrix->MarkowitzRow, *MarkoCol = Matrix->MarkowitzCol; double Product; /* Begin `UpdateMarkowitzNumbers'. */ @@ -3074,8 +3084,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 a9c4567be..7928fbe43 100644 --- a/src/maths/sparse/spOutput.c +++ b/src/maths/sparse/spOutput.c @@ -142,7 +142,7 @@ spPrint( int Header ) { -register int J = 0; +int J = 0; int I, Row, Col, Size, Top, StartCol = 1, StopCol, Columns, ElementCount = 0; double Magnitude, SmallestDiag = 0.0, SmallestElement = 0.0; double LargestElement = 0.0, LargestDiag = 0.0; @@ -418,8 +418,8 @@ spFileMatrix( int Header ) { -register int I, Size; -register ElementPtr pElement; +int I, Size; +ElementPtr pElement; int Row, Col, Err; FILE *pMatrixFile; @@ -571,9 +571,9 @@ spFileVector( #endif ) { -register int I, Size; +int I, Size; #if spCOMPLEX -register int Err; +int Err; #endif FILE *pMatrixFile; @@ -581,9 +581,12 @@ FILE *pMatrixFile; ASSERT_IS_SPARSE( Matrix ); vASSERT( RHS != NULL, "Vector missing" ); -/* Open File in append mode. */ - if ((pMatrixFile = fopen(File,"a")) == NULL) - return 0; + if (File) { + /* Open File in write mode. */ + pMatrixFile = fopen(File,"w"); + if (pMatrixFile == NULL) + return 0; + } /* Correct array pointers for ARRAY_OFFSET. */ #if NOT ARRAY_OFFSET @@ -694,8 +697,8 @@ spFileStats( char *Label ) { -register int Size, I; -register ElementPtr pElement; +int Size, I; +ElementPtr pElement; int NumberOfElements; RealNumber Data, LargestElement, SmallestElement; FILE *pStatsFile; diff --git a/src/maths/sparse/spSMP.c b/src/maths/sparse/spSMP.c index d1cbbf8e4..d16f76966 100644 --- a/src/maths/sparse/spSMP.c +++ b/src/maths/sparse/spSMP.c @@ -285,12 +285,11 @@ SMPmatrix *Matrix) int SMPnewMatrix( SMPmatrix **pMatrix, -int dummy) +int size) { -int Error; - NG_IGNORE(dummy); + int Error; - *pMatrix = (SMPmatrix *)spCreate( 0, 1, &Error ); + *pMatrix = (SMPmatrix *)spCreate( size, 1, &Error ); return Error; } @@ -562,18 +561,18 @@ SMPcAddCol(SMPmatrix *Matrix, int Accum_Col, int Addend_Col) Accum = *Prev; while (Addend != NULL) { - while (Accum && Accum->Row < Addend->Row) { - Prev = &Accum->NextInCol; - Accum = *Prev; - } - if (!Accum || Accum->Row > Addend->Row) { - Accum = spcCreateElement(Matrix, Addend->Row, Accum_Col, Prev, 0, 0); - } - Accum->Real += Addend->Real; + while (Accum && Accum->Row < Addend->Row) { + Prev = &Accum->NextInCol; + Accum = *Prev; + } + if (!Accum || Accum->Row > Addend->Row) { + Accum = spcCreateElement(Matrix, Addend->Row, Accum_Col, Prev, 0, NO); + } + Accum->Real += Addend->Real; #if spCOMPLEX - Accum->Imag += Addend->Imag; + Accum->Imag += Addend->Imag; #endif - Addend = Addend->NextInCol; + Addend = Addend->NextInCol; } return spErrorState( Matrix ); diff --git a/src/maths/sparse/spSolve.c b/src/maths/sparse/spSolve.c index 047776df3..bdf412294 100644 --- a/src/maths/sparse/spSolve.c +++ b/src/maths/sparse/spSolve.c @@ -156,13 +156,12 @@ spSolve( ) { #if REAL -register ElementPtr pElement; -register RealVector Intermediate; -register RealNumber Temp; -register int I, *pExtOrder, Size; +ElementPtr pElement; +RealVector Intermediate; +RealNumber Temp; +int I, *pExtOrder, Size; ElementPtr pPivot; #endif -void SolveComplexMatrix(); /* Begin `spSolve'. */ ASSERT_IS_SPARSE( Matrix ); @@ -301,12 +300,12 @@ SolveComplexMatrix( # endif ) { -register ElementPtr pElement; -register ComplexVector Intermediate; -register int I, *pExtOrder, Size; +ElementPtr pElement; +ComplexVector Intermediate; +int I, *pExtOrder, Size; ElementPtr pPivot; #if NOT spSEPARATED_COMPLEX_VECTORS -register ComplexVector ExtVector; +ComplexVector ExtVector; #endif ComplexNumber Temp; @@ -471,13 +470,12 @@ spSolveTransposed( ) { #if REAL -register ElementPtr pElement; -register RealVector Intermediate; -register int I, *pExtOrder, Size; +ElementPtr pElement; +RealVector Intermediate; +int I, *pExtOrder, Size; ElementPtr pPivot; RealNumber Temp; #endif -void SolveComplexTransposedMatrix(); /* Begin `spSolveTransposed'. */ ASSERT_IS_SPARSE( Matrix ); @@ -618,11 +616,11 @@ SolveComplexTransposedMatrix( # endif ) { -register ElementPtr pElement; -register ComplexVector Intermediate; -register int I, *pExtOrder, Size; +ElementPtr pElement; +ComplexVector Intermediate; +int I, *pExtOrder, Size; #if NOT spSEPARATED_COMPLEX_VECTORS -register ComplexVector ExtVector; +ComplexVector ExtVector; #endif ElementPtr pPivot; ComplexNumber Temp; diff --git a/src/maths/sparse/spUtils.c b/src/maths/sparse/spUtils.c index a932e086f..0ff8bb673 100644 --- a/src/maths/sparse/spUtils.c +++ b/src/maths/sparse/spUtils.c @@ -203,7 +203,7 @@ static RealNumber ComplexCondition( MatrixPtr, RealNumber, int* ); void spMNA_Preorder( MatrixPtr Matrix ) { -register int J, Size; +int J, Size; ElementPtr pTwin1, pTwin2; int Twins, StartAt = 1; BOOLEAN Swapped, AnotherPassNeeded; @@ -397,10 +397,9 @@ spScale( spREAL SolutionScaleFactors[] ) { -register ElementPtr pElement; -register int I, lSize, *pExtOrder; +ElementPtr pElement; +int I, lSize, *pExtOrder; RealNumber ScaleFactor; -//void ScaleComplexMatrix(); /* Begin `spScale'. */ ASSERT_IS_SPARSE( Matrix ); @@ -523,12 +522,12 @@ RealNumber ScaleFactor; static void ScaleComplexMatrix( MatrixPtr Matrix, - register RealVector RHS_ScaleFactors, - register RealVector SolutionScaleFactors + RealVector RHS_ScaleFactors, + RealVector SolutionScaleFactors ) { -register ElementPtr pElement; -register int I, lSize, *pExtOrder; +ElementPtr pElement; +int I, lSize, *pExtOrder; RealNumber ScaleFactor; /* Begin `ScaleComplexMatrix'. */ @@ -613,12 +612,11 @@ spMultiply( ) { #if REAL -register ElementPtr pElement; -register RealVector Vector; -register RealNumber Sum; -register int I, *pExtOrder; +ElementPtr pElement; +RealVector Vector; +RealNumber Sum; +int I, *pExtOrder; #endif -extern void ComplexMatrixMultiply(); /* Begin `spMultiply'. */ ASSERT_IS_SPARSE( Matrix ); @@ -710,10 +708,10 @@ ComplexMatrixMultiply( #endif ) { -register ElementPtr pElement; -register ComplexVector Vector; +ElementPtr pElement; +ComplexVector Vector; ComplexNumber Sum; -register int I, *pExtOrder; +int I, *pExtOrder; /* Begin `ComplexMatrixMultiply'. */ @@ -805,11 +803,10 @@ spMultTransposed( ) { #if REAL -register ElementPtr pElement; -register RealVector Vector; -register RealNumber Sum; -register int I, *pExtOrder; -extern void ComplexTransposedMatrixMultiply(); +ElementPtr pElement; +RealVector Vector; +RealNumber Sum; +int I, *pExtOrder; #endif /* Begin `spMultTransposed'. */ @@ -906,10 +903,10 @@ ComplexTransposedMatrixMultiply( #endif ) { -register ElementPtr pElement; -register ComplexVector Vector; +ElementPtr pElement; +ComplexVector Vector; ComplexNumber Sum; -register int I, *pExtOrder; +int I, *pExtOrder; /* Begin `ComplexTransposedMatrixMultiply'. */ @@ -999,7 +996,7 @@ 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. */ @@ -1014,7 +1011,7 @@ spDeterminant( #endif ) { -register int I, Size; +int I, Size; RealNumber Norm, nr, ni; ComplexNumber Pivot, cDeterminant; @@ -1170,7 +1167,7 @@ struct FillinListNodeStruct *pListNode; Matrix->Fillins = 0; /* Mark the fill-ins. */ - { register ElementPtr pFillin, pLastFillin; + { ElementPtr pFillin, pLastFillin; pListNode = Matrix->LastFillinListNode = Matrix->FirstFillinListNode; Matrix->FillinsRemaining = pListNode->NumberOfFillinsInList; @@ -1186,8 +1183,8 @@ struct FillinListNodeStruct *pListNode; } /* Unlink fill-ins by searching for elements marked with Row = 0. */ - { register ElementPtr pElement, *ppElement; - register int I, Size = Matrix->Size; + { ElementPtr pElement, *ppElement; + int I, Size = Matrix->Size; /* Unlink fill-ins in all columns. */ for (I = 1; I <= Size; I++) @@ -1263,7 +1260,7 @@ spDeleteRowAndCol( int Col ) { -register ElementPtr pElement, *ppElement, pLastElement; +ElementPtr pElement, *ppElement, pLastElement; int Size, ExtRow, ExtCol; /* Begin `spDeleteRowAndCol'. */ @@ -1370,8 +1367,8 @@ int Size, ExtRow, ExtCol; spREAL spPseudoCondition( MatrixPtr Matrix ) { - register int I; - register ArrayOfElementPtrs Diag; + int I; + ArrayOfElementPtrs Diag; RealNumber MaxPivot, MinPivot, Mag; /* Begin `spPseudoCondition'. */ @@ -1458,9 +1455,9 @@ spCondition( int *pError ) { -register ElementPtr pElement; -register RealVector T, Tm; -register int I, K, Row; +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; @@ -1626,7 +1623,7 @@ RealNumber Linpack, OLeary, InvNormOfInverse, ComplexCondition(); for (ASz = 0.0, I = Size; I > 0; I--) ASz += ABS(T[I]); #if NOT spCOMPLEX - FREE( Tm ); + SP_FREE( Tm ); #endif Linpack = ASy / ASz; @@ -1669,9 +1666,9 @@ ComplexCondition( int *pError ) { -register ElementPtr pElement; -register ComplexVector T, Tm; -register int I, K, Row; +ElementPtr pElement; +ComplexVector T, Tm; +int I, K, Row; ElementPtr pPivot; int Size; 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. */ for (ASz = 0.0, I = Size; I > 0; I--) ASz += CMPLX_1_NORM(T[I]); - FREE( Tm ); + SP_FREE( Tm ); Linpack = ASy / ASz; OLeary = E / MaxY; @@ -1854,8 +1851,8 @@ ComplexNumber Wp, Wm; spREAL spNorm( MatrixPtr Matrix ) { -register ElementPtr pElement; -register int I; +ElementPtr pElement; +int I; RealNumber Max = 0.0, AbsRowSum; /* Begin `spNorm'. */ @@ -1968,11 +1965,11 @@ RealNumber Max = 0.0, AbsRowSum; spREAL spLargestElement( MatrixPtr Matrix ) { -register int I; +int I; RealNumber Mag, AbsColSum, Max = 0.0, MaxRow = 0.0, MaxCol = 0.0; RealNumber Pivot; ComplexNumber cPivot; -register ElementPtr pElement, pDiag; +ElementPtr pElement, pDiag; /* Begin `spLargestElement'. */ ASSERT_IS_SPARSE( Matrix ); @@ -2089,8 +2086,8 @@ spRoundoff( spREAL Rho ) { -register ElementPtr pElement; -register int Count, I, MaxCount = 0; +ElementPtr pElement; +int Count, I, MaxCount = 0; RealNumber Reid, Gear; /* Begin `spRoundoff'. */