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  <input>    (MatrixPtr)
  *      Pointer to the matrix.
- *  AllocatedPtr  <input>  (void *)
+ *  AllocatedPtr  <input>  
  *      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'. */