File fasthenry-3.0.patch of Package fasthenry

Overview Repositories Revisions Requests Users Attributes Meta

File fasthenry-3.0.patch of Package fasthenry (Revision 1fde13fe4b3f7976b0aba7cc6ec2a760)

Currently displaying revision 1fde13fe4b3f7976b0aba7cc6ec2a760 , Show latest

diff -uNr fasthenry-3.0.orig/src/fasthenry/sparse/spAllocate.c fasthenry-3.0/src/fasthenry/sparse/spAllocate.c
--- fasthenry-3.0.orig/src/fasthenry/sparse/spAllocate.c	2004-08-03 17:17:17.000000000 -0500
+++ fasthenry-3.0/src/fasthenry/sparse/spAllocate.c	2009-08-05 13:31:48.000000000 -0500
@@ -108,206 +108,110 @@
 
 
 
-
-
 
 /*
- *  MATRIX ALLOCATION
- *
- *  Allocates and initializes the data structures associated with a matrix.
+ *  ADD A BLOCK OF SLOTS TO ALLOCATION LIST     
  *
- *  >>> Returned:
- *  A pointer to the matrix is returned cast into the form of a pointer to
- *  a character.  This pointer is then passed and used by the other matrix
- *  routines to refer to a particular matrix.  If an error occurs, the NULL
- *  pointer is returned.
+ *  This routine increases the size of the allocation list.
  *
  *  >>> Arguments:
- *  Size  <input>  (int)
- *      Size of matrix or estimate of size of matrix if matrix is EXPANDABLE.
- *  Complex  <input>  (int)
- *      Type of matrix.  If Complex is 0 then the matrix is real, otherwise
- *      the matrix will be complex.  Note that if the routines are not set up
- *      to handle the type of matrix requested, then a spPANIC error will occur.
- *      Further note that if a matrix will be both real and complex, it must
- *      be specified here as being complex.
- *  pError  <output>  (int *)
- *      Returns error flag, needed because function spError() will not work
- *      correctly if spCreate() returns NULL.
+ *  Matrix  <input>    (MatrixPtr)
+ *      Pointer to the matrix.
  *
  *  >>> Local variables:
- *  AllocatedSize  (int)
- *      The size of the matrix being allocated.
- *  Matrix  (MatrixPtr)
- *      A pointer to the matrix frame being created.
+ *  ListPtr  (AllocationListPtr)
+ *      Pointer to the list that contains the pointers to segments of memory
+ *      that were allocated by the operating system for the current matrix.
  *
  *  >>> Possible errors:
  *  spNO_MEMORY
- *  spPANIC
- *  Error is cleared in this routine.
  */
 
-char *
-spCreate( Size, Complex, pError )
+static
+AllocateBlockOfAllocationList( Matrix )
 
-int  Size, *pError;
-BOOLEAN  Complex;
+MatrixPtr Matrix;
 {
-register  unsigned  SizePlusOne;
-register  MatrixPtr  Matrix;
 register  int  I;
-int  AllocatedSize;
-
-/* Begin `spCreate'. */
-/* Clear error flag. */
-    *pError = spOKAY;
-
-/* Test for valid size. */
-    if ((Size < 0) OR (Size == 0 AND NOT EXPANDABLE))
-    {   *pError = spPANIC;
-        return NULL;
-    }
+register  AllocationListPtr  ListPtr;
 
-/* Test for valid type. */
-#if NOT spCOMPLEX
-    if (Complex)
-    {   *pError = spPANIC;
-        return NULL;
-    }
-#endif
-#if NOT REAL
-    if (NOT Complex)
-    {   *pError = spPANIC;
-        return NULL;
+/* Begin `AllocateBlockOfAllocationList'. */
+/* Allocate block of records for allocation list. */
+    ListPtr = ALLOC(struct AllocationRecord, (ELEMENTS_PER_ALLOCATION+1));
+    if (ListPtr == NULL)
+    {   Matrix->Error = spNO_MEMORY;
+        return;
     }
-#endif
 
-/* Create Matrix. */
-    AllocatedSize = MAX( Size, MINIMUM_ALLOCATED_SIZE );
-    SizePlusOne = (unsigned)(AllocatedSize + 1);
+/* String entries of allocation list into singly linked list.  List is linked
+   such that any record points to the one before it. */
 
-    if ((Matrix = ALLOC(struct MatrixFrame, 1)) == NULL)
-    {   *pError = spNO_MEMORY;
-        return NULL;
+    ListPtr->NextRecord = Matrix->TopOfAllocationList;
+    Matrix->TopOfAllocationList = ListPtr;
+    ListPtr += ELEMENTS_PER_ALLOCATION;
+    for (I = ELEMENTS_PER_ALLOCATION; I > 0; I--)
+    {    ListPtr->NextRecord = ListPtr - 1;
+         ListPtr--;
     }
 
-/* Initialize matrix */
-    Matrix->ID = SPARSE_ID;
-    Matrix->Complex = Complex;
-    Matrix->PreviousMatrixWasComplex = Complex;
-    Matrix->Factored = NO;
-    Matrix->Elements = 0;
-    Matrix->Error = *pError;
-    Matrix->Fillins = 0;
-    Matrix->Reordered = NO;
-    Matrix->NeedsOrdering = YES;
-    Matrix->NumberOfInterchangesIsOdd = NO;
-    Matrix->Partitioned = NO;
-    Matrix->RowsLinked = NO;
-    Matrix->InternalVectorsAllocated = NO;
-    Matrix->SingularCol = 0;
-    Matrix->SingularRow = 0;
-    Matrix->Size = Size;
-    Matrix->AllocatedSize = AllocatedSize;
-    Matrix->ExtSize = Size;
-    Matrix->AllocatedExtSize = AllocatedSize;
-    Matrix->CurrentSize = 0;
-    Matrix->ExtToIntColMap = NULL;
-    Matrix->ExtToIntRowMap = NULL;
-    Matrix->IntToExtColMap = NULL;
-    Matrix->IntToExtRowMap = NULL;
-    Matrix->MarkowitzRow = NULL;
-    Matrix->MarkowitzCol = NULL;
-    Matrix->MarkowitzProd = NULL;
-    Matrix->DoCmplxDirect = NULL;
-    Matrix->DoRealDirect = NULL;
-    Matrix->Intermediate = NULL;
-    Matrix->RelThreshold = DEFAULT_THRESHOLD;
-    Matrix->AbsThreshold = 0.0;
-
-    Matrix->TopOfAllocationList = NULL;
-    Matrix->RecordsRemaining = 0;
-    Matrix->ElementsRemaining = 0;
-    Matrix->FillinsRemaining = 0;
-
-    RecordAllocation( Matrix, (char *)Matrix );
-    if (Matrix->Error == spNO_MEMORY) goto MemoryError;
+/* Record allocation of space for allocation list on allocation list. */
+    Matrix->TopOfAllocationList->AllocatedPtr = (char *)ListPtr;
+    Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION;
 
-/* Take out the trash. */
-    Matrix->TrashCan.Real = 0.0;
-#if spCOMPLEX
-    Matrix->TrashCan.Imag = 0.0;
-#endif
-    Matrix->TrashCan.Row = 0;
-    Matrix->TrashCan.Col = 0;
-    Matrix->TrashCan.NextInRow = NULL;
-    Matrix->TrashCan.NextInCol = NULL;
-#if INITIALIZE
-    Matrix->TrashCan.pInitInfo = NULL;
-#endif
+    return;
+}
 
-/* Allocate space in memory for Diag pointer vector. */
-    CALLOC( Matrix->Diag, ElementPtr, SizePlusOne);
-    if (Matrix->Diag == NULL)
-        goto MemoryError;
 
-/* Allocate space in memory for FirstInCol pointer vector. */
-    CALLOC( Matrix->FirstInCol, ElementPtr, SizePlusOne);
-    if (Matrix->FirstInCol == NULL)
-        goto MemoryError;
 
-/* Allocate space in memory for FirstInRow pointer vector. */
-    CALLOC( Matrix->FirstInRow, ElementPtr, SizePlusOne);
-    if (Matrix->FirstInRow == NULL)
-        goto MemoryError;
 
-/* Allocate space in memory for IntToExtColMap vector. */
-    if (( Matrix->IntToExtColMap = ALLOC(int, SizePlusOne)) == NULL)
-        goto MemoryError;
+
+/*
+ *  RECORD A MEMORY ALLOCATION
+ *
+ *  This routine is used to record all memory allocations so that the memory
+ *  can be freed later.
+ *
+ *  >>> Arguments:
+ *  Matrix  <input>    (MatrixPtr)
+ *      Pointer to the matrix.
+ *  AllocatedPtr  <input>  (char *)
+ *      The pointer returned by malloc or calloc.  These pointers are saved in
+ *      a list so that they can be easily freed.
+ *
+ *  >>> Possible errors:
+ *  spNO_MEMORY
+ */
 
-/* Allocate space in memory for IntToExtRowMap vector. */
-    if (( Matrix->IntToExtRowMap = ALLOC(int, SizePlusOne)) == NULL)
-        goto MemoryError;
+static
+RecordAllocation( Matrix, AllocatedPtr )
 
-/* Initialize MapIntToExt vectors. */
-    for (I = 1; I <= AllocatedSize; I++)
-    {   Matrix->IntToExtRowMap[I] = I;
-        Matrix->IntToExtColMap[I] = I;
+MatrixPtr Matrix;
+char  *AllocatedPtr;
+{
+/* Begin `RecordAllocation'. */
+/*
+ * If Allocated pointer is NULL, assume that malloc returned a NULL pointer,
+ * which indicates a spNO_MEMORY error.
+ */
+    if (AllocatedPtr == NULL)
+    {   Matrix->Error = spNO_MEMORY;
+        return;
     }
 
-#if TRANSLATE
-/* Allocate space in memory for ExtToIntColMap vector. */
-    if (( Matrix->ExtToIntColMap = ALLOC(int, SizePlusOne)) == NULL)
-        goto MemoryError;
-
-/* Allocate space in memory for ExtToIntRowMap vector. */
-    if (( Matrix->ExtToIntRowMap = ALLOC(int, SizePlusOne)) == NULL)
-        goto MemoryError;
-
-/* Initialize MapExtToInt vectors. */
-    for (I = 1; I <= AllocatedSize; I++)
-    {   Matrix->ExtToIntColMap[I] = -1;
-        Matrix->ExtToIntRowMap[I] = -1;
+/* Allocate block of MatrixElements if necessary. */
+    if (Matrix->RecordsRemaining == 0)
+    {   AllocateBlockOfAllocationList( Matrix );
+        if (Matrix->Error == spNO_MEMORY)
+        {   FREE(AllocatedPtr);
+            return;
+        }
     }
-    Matrix->ExtToIntColMap[0] = 0;
-    Matrix->ExtToIntRowMap[0] = 0;
-#endif
 
-/* Allocate space for fill-ins and initial set of elements. */
-    InitializeElementBlocks( Matrix, SPACE_FOR_ELEMENTS*AllocatedSize,
-                                     SPACE_FOR_FILL_INS*AllocatedSize );
-    if (Matrix->Error == spNO_MEMORY)
-        goto MemoryError;
-
-    return (char *)Matrix;
-
-MemoryError:
+/* Add Allocated pointer to Allocation List. */
+    (++Matrix->TopOfAllocationList)->AllocatedPtr = AllocatedPtr;
+    Matrix->RecordsRemaining--;
+    return;
 
-/* Deallocate matrix and return no pointer to matrix if there is not enough
-   memory. */
-    *pError = spNO_MEMORY;
-    spDestroy( (char *)Matrix);
-    return NULL;
 }
 
 
@@ -371,6 +275,10 @@
 
 
 
+
+
+
+
 
 /*
  *  ELEMENT ALLOCATION INITIALIZATION
@@ -515,66 +423,9 @@
     }
 #endif
 
-/* Update Fill-in counter and return pointer to Fill-in. */
-    Matrix->FillinsRemaining--;
-    return Matrix->NextAvailFillin++;
-}
-
-
-
-
-
-
-
-
-
-/*
- *  RECORD A MEMORY ALLOCATION
- *
- *  This routine is used to record all memory allocations so that the memory
- *  can be freed later.
- *
- *  >>> Arguments:
- *  Matrix  <input>    (MatrixPtr)
- *      Pointer to the matrix.
- *  AllocatedPtr  <input>  (char *)
- *      The pointer returned by malloc or calloc.  These pointers are saved in
- *      a list so that they can be easily freed.
- *
- *  >>> Possible errors:
- *  spNO_MEMORY
- */
-
-static
-RecordAllocation( Matrix, AllocatedPtr )
-
-MatrixPtr Matrix;
-char  *AllocatedPtr;
-{
-/* Begin `RecordAllocation'. */
-/*
- * If Allocated pointer is NULL, assume that malloc returned a NULL pointer,
- * which indicates a spNO_MEMORY error.
- */
-    if (AllocatedPtr == NULL)
-    {   Matrix->Error = spNO_MEMORY;
-        return;
-    }
-
-/* Allocate block of MatrixElements if necessary. */
-    if (Matrix->RecordsRemaining == 0)
-    {   AllocateBlockOfAllocationList( Matrix );
-        if (Matrix->Error == spNO_MEMORY)
-        {   FREE(AllocatedPtr);
-            return;
-        }
-    }
-
-/* Add Allocated pointer to Allocation List. */
-    (++Matrix->TopOfAllocationList)->AllocatedPtr = AllocatedPtr;
-    Matrix->RecordsRemaining--;
-    return;
-
+/* Update Fill-in counter and return pointer to Fill-in. */
+    Matrix->FillinsRemaining--;
+    return Matrix->NextAvailFillin++;
 }
 
 
@@ -583,58 +434,8 @@
 
 
 
-
-/*
- *  ADD A BLOCK OF SLOTS TO ALLOCATION LIST     
- *
- *  This routine increases the size of the allocation list.
- *
- *  >>> Arguments:
- *  Matrix  <input>    (MatrixPtr)
- *      Pointer to the matrix.
- *
- *  >>> Local variables:
- *  ListPtr  (AllocationListPtr)
- *      Pointer to the list that contains the pointers to segments of memory
- *      that were allocated by the operating system for the current matrix.
- *
- *  >>> Possible errors:
- *  spNO_MEMORY
- */
-
-static
-AllocateBlockOfAllocationList( Matrix )
-
-MatrixPtr Matrix;
-{
-register  int  I;
-register  AllocationListPtr  ListPtr;
-
-/* Begin `AllocateBlockOfAllocationList'. */
-/* Allocate block of records for allocation list. */
-    ListPtr = ALLOC(struct AllocationRecord, (ELEMENTS_PER_ALLOCATION+1));
-    if (ListPtr == NULL)
-    {   Matrix->Error = spNO_MEMORY;
-        return;
-    }
-
-/* String entries of allocation list into singly linked list.  List is linked
-   such that any record points to the one before it. */
-
-    ListPtr->NextRecord = Matrix->TopOfAllocationList;
-    Matrix->TopOfAllocationList = ListPtr;
-    ListPtr += ELEMENTS_PER_ALLOCATION;
-    for (I = ELEMENTS_PER_ALLOCATION; I > 0; I--)
-    {    ListPtr->NextRecord = ListPtr - 1;
-         ListPtr--;
-    }
 
-/* Record allocation of space for allocation list on allocation list. */
-    Matrix->TopOfAllocationList->AllocatedPtr = (char *)ListPtr;
-    Matrix->RecordsRemaining = ELEMENTS_PER_ALLOCATION;
 
-    return;
-}
 
 
 
@@ -902,3 +703,206 @@
     ASSERT( IS_SPARSE( (MatrixPtr)eMatrix ) );
     return ((MatrixPtr)eMatrix)->Elements;
 }
+
+
+
+
+/*
+ *  MATRIX ALLOCATION
+ *
+ *  Allocates and initializes the data structures associated with a matrix.
+ *
+ *  >>> Returned:
+ *  A pointer to the matrix is returned cast into the form of a pointer to
+ *  a character.  This pointer is then passed and used by the other matrix
+ *  routines to refer to a particular matrix.  If an error occurs, the NULL
+ *  pointer is returned.
+ *
+ *  >>> Arguments:
+ *  Size  <input>  (int)
+ *      Size of matrix or estimate of size of matrix if matrix is EXPANDABLE.
+ *  Complex  <input>  (int)
+ *      Type of matrix.  If Complex is 0 then the matrix is real, otherwise
+ *      the matrix will be complex.  Note that if the routines are not set up
+ *      to handle the type of matrix requested, then a spPANIC error will occur.
+ *      Further note that if a matrix will be both real and complex, it must
+ *      be specified here as being complex.
+ *  pError  <output>  (int *)
+ *      Returns error flag, needed because function spError() will not work
+ *      correctly if spCreate() returns NULL.
+ *
+ *  >>> Local variables:
+ *  AllocatedSize  (int)
+ *      The size of the matrix being allocated.
+ *  Matrix  (MatrixPtr)
+ *      A pointer to the matrix frame being created.
+ *
+ *  >>> Possible errors:
+ *  spNO_MEMORY
+ *  spPANIC
+ *  Error is cleared in this routine.
+ */
+
+char *
+spCreate( Size, Complex, pError )
+
+int  Size, *pError;
+BOOLEAN  Complex;
+{
+register  unsigned  SizePlusOne;
+register  MatrixPtr  Matrix;
+register  int  I;
+int  AllocatedSize;
+
+/* Begin `spCreate'. */
+/* Clear error flag. */
+    *pError = spOKAY;
+
+/* Test for valid size. */
+    if ((Size < 0) OR (Size == 0 AND NOT EXPANDABLE))
+    {   *pError = spPANIC;
+        return NULL;
+    }
+
+/* Test for valid type. */
+#if NOT spCOMPLEX
+    if (Complex)
+    {   *pError = spPANIC;
+        return NULL;
+    }
+#endif
+#if NOT REAL
+    if (NOT Complex)
+    {   *pError = spPANIC;
+        return NULL;
+    }
+#endif
+
+/* Create Matrix. */
+    AllocatedSize = MAX( Size, MINIMUM_ALLOCATED_SIZE );
+    SizePlusOne = (unsigned)(AllocatedSize + 1);
+
+    if ((Matrix = ALLOC(struct MatrixFrame, 1)) == NULL)
+    {   *pError = spNO_MEMORY;
+        return NULL;
+    }
+
+/* Initialize matrix */
+    Matrix->ID = SPARSE_ID;
+    Matrix->Complex = Complex;
+    Matrix->PreviousMatrixWasComplex = Complex;
+    Matrix->Factored = NO;
+    Matrix->Elements = 0;
+    Matrix->Error = *pError;
+    Matrix->Fillins = 0;
+    Matrix->Reordered = NO;
+    Matrix->NeedsOrdering = YES;
+    Matrix->NumberOfInterchangesIsOdd = NO;
+    Matrix->Partitioned = NO;
+    Matrix->RowsLinked = NO;
+    Matrix->InternalVectorsAllocated = NO;
+    Matrix->SingularCol = 0;
+    Matrix->SingularRow = 0;
+    Matrix->Size = Size;
+    Matrix->AllocatedSize = AllocatedSize;
+    Matrix->ExtSize = Size;
+    Matrix->AllocatedExtSize = AllocatedSize;
+    Matrix->CurrentSize = 0;
+    Matrix->ExtToIntColMap = NULL;
+    Matrix->ExtToIntRowMap = NULL;
+    Matrix->IntToExtColMap = NULL;
+    Matrix->IntToExtRowMap = NULL;
+    Matrix->MarkowitzRow = NULL;
+    Matrix->MarkowitzCol = NULL;
+    Matrix->MarkowitzProd = NULL;
+    Matrix->DoCmplxDirect = NULL;
+    Matrix->DoRealDirect = NULL;
+    Matrix->Intermediate = NULL;
+    Matrix->RelThreshold = DEFAULT_THRESHOLD;
+    Matrix->AbsThreshold = 0.0;
+
+    Matrix->TopOfAllocationList = NULL;
+    Matrix->RecordsRemaining = 0;
+    Matrix->ElementsRemaining = 0;
+    Matrix->FillinsRemaining = 0;
+
+    RecordAllocation( Matrix, (char *)Matrix );
+    if (Matrix->Error == spNO_MEMORY) goto MemoryError;
+
+/* Take out the trash. */
+    Matrix->TrashCan.Real = 0.0;
+#if spCOMPLEX
+    Matrix->TrashCan.Imag = 0.0;
+#endif
+    Matrix->TrashCan.Row = 0;
+    Matrix->TrashCan.Col = 0;
+    Matrix->TrashCan.NextInRow = NULL;
+    Matrix->TrashCan.NextInCol = NULL;
+#if INITIALIZE
+    Matrix->TrashCan.pInitInfo = NULL;
+#endif
+
+/* Allocate space in memory for Diag pointer vector. */
+    CALLOC( Matrix->Diag, ElementPtr, SizePlusOne);
+    if (Matrix->Diag == NULL)
+        goto MemoryError;
+
+/* Allocate space in memory for FirstInCol pointer vector. */
+    CALLOC( Matrix->FirstInCol, ElementPtr, SizePlusOne);
+    if (Matrix->FirstInCol == NULL)
+        goto MemoryError;
+
+/* Allocate space in memory for FirstInRow pointer vector. */
+    CALLOC( Matrix->FirstInRow, ElementPtr, SizePlusOne);
+    if (Matrix->FirstInRow == NULL)
+        goto MemoryError;
+
+/* Allocate space in memory for IntToExtColMap vector. */
+    if (( Matrix->IntToExtColMap = ALLOC(int, SizePlusOne)) == NULL)
+        goto MemoryError;
+
+/* Allocate space in memory for IntToExtRowMap vector. */
+    if (( Matrix->IntToExtRowMap = ALLOC(int, SizePlusOne)) == NULL)
+        goto MemoryError;
+
+/* Initialize MapIntToExt vectors. */
+    for (I = 1; I <= AllocatedSize; I++)
+    {   Matrix->IntToExtRowMap[I] = I;
+        Matrix->IntToExtColMap[I] = I;
+    }
+
+#if TRANSLATE
+/* Allocate space in memory for ExtToIntColMap vector. */
+    if (( Matrix->ExtToIntColMap = ALLOC(int, SizePlusOne)) == NULL)
+        goto MemoryError;
+
+/* Allocate space in memory for ExtToIntRowMap vector. */
+    if (( Matrix->ExtToIntRowMap = ALLOC(int, SizePlusOne)) == NULL)
+        goto MemoryError;
+
+/* Initialize MapExtToInt vectors. */
+    for (I = 1; I <= AllocatedSize; I++)
+    {   Matrix->ExtToIntColMap[I] = -1;
+        Matrix->ExtToIntRowMap[I] = -1;
+    }
+    Matrix->ExtToIntColMap[0] = 0;
+    Matrix->ExtToIntRowMap[0] = 0;
+#endif
+
+/* Allocate space for fill-ins and initial set of elements. */
+    InitializeElementBlocks( Matrix, SPACE_FOR_ELEMENTS*AllocatedSize,
+                                     SPACE_FOR_FILL_INS*AllocatedSize );
+    if (Matrix->Error == spNO_MEMORY)
+        goto MemoryError;
+
+    return (char *)Matrix;
+
+MemoryError:
+
+/* Deallocate matrix and return no pointer to matrix if there is not enough
+   memory. */
+    *pError = spNO_MEMORY;
+    spDestroy( (char *)Matrix);
+    return NULL;
+}
+
diff -uNr fasthenry-3.0.orig/src/fasthenry/sparse/spBuild.c fasthenry-3.0/src/fasthenry/sparse/spBuild.c
--- fasthenry-3.0.orig/src/fasthenry/sparse/spBuild.c	2004-08-03 17:17:18.000000000 -0500
+++ fasthenry-3.0/src/fasthenry/sparse/spBuild.c	2009-08-05 13:36:03.000000000 -0500
@@ -176,6 +176,170 @@
 
 
 
+#if TRANSLATE
+
+/*
+ *  EXPAND TRANSLATION ARRAYS
+ *
+ *  Increases the size arrays that are used to translate external to internal
+ *  row and column numbers.
+ *
+ *  >>> Arguments:
+ *  Matrix  <input>    (MatrixPtr)
+ *      Pointer to the matrix.
+ *  NewSize  <input>  (int)
+ *     The new size of the translation arrays.
+ *
+ *  >>> Local variables:
+ *  OldAllocatedSize  (int)
+ *     The allocated size of the translation arrays before being expanded.
+ */
+
+static
+ExpandTranslationArrays( Matrix, NewSize )
+
+MatrixPtr Matrix;
+register int  NewSize;
+{
+register int I, OldAllocatedSize = Matrix->AllocatedExtSize;
+
+/* Begin `ExpandTranslationArrays'. */
+    Matrix->ExtSize = NewSize;
+
+    if (NewSize <= OldAllocatedSize)
+        return;
+
+/* Expand the translation arrays ExtToIntRowMap and ExtToIntColMap. */
+    NewSize = MAX( NewSize, EXPANSION_FACTOR * OldAllocatedSize );
+    Matrix->AllocatedExtSize = NewSize;
+
+    if (( REALLOC(Matrix->ExtToIntRowMap, int, NewSize+1)) == NULL)
+    {   Matrix->Error = spNO_MEMORY;
+        return;
+    }
+    if (( REALLOC(Matrix->ExtToIntColMap, int, NewSize+1)) == NULL)
+    {   Matrix->Error = spNO_MEMORY;
+        return;
+    }
+
+/* Initialize the new portion of the vectors. */
+    for (I = OldAllocatedSize+1; I <= NewSize; I++)
+    {   Matrix->ExtToIntRowMap[I] = -1;
+        Matrix->ExtToIntColMap[I] = -1;
+    }
+
+    return;
+}
+#endif
+
+
+#if TRANSLATE
+
+/*
+ *  TRANSLATE EXTERNAL INDICES TO INTERNAL
+ *
+ *  Convert internal row and column numbers to internal row and column numbers.
+ *  Also updates Ext/Int maps.
+ *
+ *
+ *  >>> Arguments:
+ *  Matrix  <input>    (MatrixPtr)
+ *      Pointer to the matrix.
+ *  Row  <input/output>  (int *)
+ *     Upon entry Row is either a external row number of an external node
+ *     number.  Upon entry, the internal equivalent is supplied.
+ *  Col  <input/output>  (int *)
+ *     Upon entry Column is either a external column number of an external node
+ *     number.  Upon entry, the internal equivalent is supplied.
+ *
+ *  >>> Local variables:
+ *  ExtCol  (int)
+ *     Temporary variable used to hold the external column or node number
+ *     during the external to internal column number translation.
+ *  ExtRow  (int)
+ *     Temporary variable used to hold the external row or node number during
+ *     the external to internal row number translation.
+ *  IntCol  (int)
+ *     Temporary variable used to hold the internal column or node number
+ *     during the external to internal column number translation.
+ *  IntRow  (int)
+ *     Temporary variable used to hold the internal row or node number during
+ *     the external to internal row number translation.
+ */
+
+static void
+Translate( Matrix, Row, Col )
+
+MatrixPtr Matrix;
+int  *Row, *Col;
+{
+register int IntRow, IntCol, ExtRow, ExtCol;
+
+/* Begin `Translate'. */
+    ExtRow = *Row;
+    ExtCol = *Col;
+
+/* Expand translation arrays if necessary. */
+    if ((ExtRow > Matrix->AllocatedExtSize) OR
+        (ExtCol > Matrix->AllocatedExtSize))
+    {
+        ExpandTranslationArrays( Matrix, MAX(ExtRow, ExtCol) );
+        if (Matrix->Error == spNO_MEMORY) return;
+    }
+
+/* Set ExtSize if necessary. */
+    if ((ExtRow > Matrix->ExtSize) OR (ExtCol > Matrix->ExtSize))
+        Matrix->ExtSize = MAX(ExtRow, ExtCol);
+
+/* Translate external row or node number to internal row or node number. */
+    if ((IntRow = Matrix->ExtToIntRowMap[ExtRow]) == -1)
+    {   Matrix->ExtToIntRowMap[ExtRow] = ++Matrix->CurrentSize;
+        Matrix->ExtToIntColMap[ExtRow] = Matrix->CurrentSize;
+        IntRow = Matrix->CurrentSize;
+
+#if NOT EXPANDABLE
+        ASSERT(IntRow <= Matrix->Size);
+#endif
+
+#if EXPANDABLE
+/* Re-size Matrix if necessary. */
+        if (IntRow > Matrix->Size)
+            EnlargeMatrix( Matrix, IntRow );
+        if (Matrix->Error == spNO_MEMORY) return;
+#endif
+
+        Matrix->IntToExtRowMap[IntRow] = ExtRow;
+        Matrix->IntToExtColMap[IntRow] = ExtRow;
+    }
+
+/* Translate external column or node number to internal column or node number.*/
+    if ((IntCol = Matrix->ExtToIntColMap[ExtCol]) == -1)
+    {   Matrix->ExtToIntRowMap[ExtCol] = ++Matrix->CurrentSize;
+        Matrix->ExtToIntColMap[ExtCol] = Matrix->CurrentSize;
+        IntCol = Matrix->CurrentSize;
+
+#if NOT EXPANDABLE
+        ASSERT(IntCol <= Matrix->Size);
+#endif
+
+#if EXPANDABLE
+/* Re-size Matrix if necessary. */
+        if (IntCol > Matrix->Size)
+            EnlargeMatrix( Matrix, IntCol );
+        if (Matrix->Error == spNO_MEMORY) return;
+#endif
+
+        Matrix->IntToExtRowMap[IntCol] = ExtCol;
+        Matrix->IntToExtColMap[IntCol] = ExtCol;
+    }
+
+    *Row = IntRow;
+    *Col = IntCol;
+    return;
+}
+#endif
+
+
 
 
 
@@ -368,112 +532,6 @@
 
 
 
-#if TRANSLATE
-
-/*
- *  TRANSLATE EXTERNAL INDICES TO INTERNAL
- *
- *  Convert internal row and column numbers to internal row and column numbers.
- *  Also updates Ext/Int maps.
- *
- *
- *  >>> Arguments:
- *  Matrix  <input>    (MatrixPtr)
- *      Pointer to the matrix.
- *  Row  <input/output>  (int *)
- *     Upon entry Row is either a external row number of an external node
- *     number.  Upon entry, the internal equivalent is supplied.
- *  Col  <input/output>  (int *)
- *     Upon entry Column is either a external column number of an external node
- *     number.  Upon entry, the internal equivalent is supplied.
- *
- *  >>> Local variables:
- *  ExtCol  (int)
- *     Temporary variable used to hold the external column or node number
- *     during the external to internal column number translation.
- *  ExtRow  (int)
- *     Temporary variable used to hold the external row or node number during
- *     the external to internal row number translation.
- *  IntCol  (int)
- *     Temporary variable used to hold the internal column or node number
- *     during the external to internal column number translation.
- *  IntRow  (int)
- *     Temporary variable used to hold the internal row or node number during
- *     the external to internal row number translation.
- */
-
-static void
-Translate( Matrix, Row, Col )
-
-MatrixPtr Matrix;
-int  *Row, *Col;
-{
-register int IntRow, IntCol, ExtRow, ExtCol;
-
-/* Begin `Translate'. */
-    ExtRow = *Row;
-    ExtCol = *Col;
-
-/* Expand translation arrays if necessary. */
-    if ((ExtRow > Matrix->AllocatedExtSize) OR
-        (ExtCol > Matrix->AllocatedExtSize))
-    {
-        ExpandTranslationArrays( Matrix, MAX(ExtRow, ExtCol) );
-        if (Matrix->Error == spNO_MEMORY) return;
-    }
-
-/* Set ExtSize if necessary. */
-    if ((ExtRow > Matrix->ExtSize) OR (ExtCol > Matrix->ExtSize))
-        Matrix->ExtSize = MAX(ExtRow, ExtCol);
-
-/* Translate external row or node number to internal row or node number. */
-    if ((IntRow = Matrix->ExtToIntRowMap[ExtRow]) == -1)
-    {   Matrix->ExtToIntRowMap[ExtRow] = ++Matrix->CurrentSize;
-        Matrix->ExtToIntColMap[ExtRow] = Matrix->CurrentSize;
-        IntRow = Matrix->CurrentSize;
-
-#if NOT EXPANDABLE
-        ASSERT(IntRow <= Matrix->Size);
-#endif
-
-#if EXPANDABLE
-/* Re-size Matrix if necessary. */
-        if (IntRow > Matrix->Size)
-            EnlargeMatrix( Matrix, IntRow );
-        if (Matrix->Error == spNO_MEMORY) return;
-#endif
-
-        Matrix->IntToExtRowMap[IntRow] = ExtRow;
-        Matrix->IntToExtColMap[IntRow] = ExtRow;
-    }
-
-/* Translate external column or node number to internal column or node number.*/
-    if ((IntCol = Matrix->ExtToIntColMap[ExtCol]) == -1)
-    {   Matrix->ExtToIntRowMap[ExtCol] = ++Matrix->CurrentSize;
-        Matrix->ExtToIntColMap[ExtCol] = Matrix->CurrentSize;
-        IntCol = Matrix->CurrentSize;
-
-#if NOT EXPANDABLE
-        ASSERT(IntCol <= Matrix->Size);
-#endif
-
-#if EXPANDABLE
-/* Re-size Matrix if necessary. */
-        if (IntCol > Matrix->Size)
-            EnlargeMatrix( Matrix, IntCol );
-        if (Matrix->Error == spNO_MEMORY) return;
-#endif
-
-        Matrix->IntToExtRowMap[IntCol] = ExtCol;
-        Matrix->IntToExtColMap[IntCol] = ExtCol;
-    }
-
-    *Row = IntRow;
-    *Col = IntCol;
-    return;
-}
-#endif
-
 
 
 
@@ -1014,62 +1072,6 @@
 
 
 
-#if TRANSLATE
-
-/*
- *  EXPAND TRANSLATION ARRAYS
- *
- *  Increases the size arrays that are used to translate external to internal
- *  row and column numbers.
- *
- *  >>> Arguments:
- *  Matrix  <input>    (MatrixPtr)
- *      Pointer to the matrix.
- *  NewSize  <input>  (int)
- *     The new size of the translation arrays.
- *
- *  >>> Local variables:
- *  OldAllocatedSize  (int)
- *     The allocated size of the translation arrays before being expanded.
- */
-
-static
-ExpandTranslationArrays( Matrix, NewSize )
-
-MatrixPtr Matrix;
-register int  NewSize;
-{
-register int I, OldAllocatedSize = Matrix->AllocatedExtSize;
-
-/* Begin `ExpandTranslationArrays'. */
-    Matrix->ExtSize = NewSize;
-
-    if (NewSize <= OldAllocatedSize)
-        return;
-
-/* Expand the translation arrays ExtToIntRowMap and ExtToIntColMap. */
-    NewSize = MAX( NewSize, EXPANSION_FACTOR * OldAllocatedSize );
-    Matrix->AllocatedExtSize = NewSize;
-
-    if (( REALLOC(Matrix->ExtToIntRowMap, int, NewSize+1)) == NULL)
-    {   Matrix->Error = spNO_MEMORY;
-        return;
-    }
-    if (( REALLOC(Matrix->ExtToIntColMap, int, NewSize+1)) == NULL)
-    {   Matrix->Error = spNO_MEMORY;
-        return;
-    }
-
-/* Initialize the new portion of the vectors. */
-    for (I = OldAllocatedSize+1; I <= NewSize; I++)
-    {   Matrix->ExtToIntRowMap[I] = -1;
-        Matrix->ExtToIntColMap[I] = -1;
-    }
-
-    return;
-}
-#endif
-
 
 
 
diff -uNr fasthenry-3.0.orig/src/fasthenry/sparse/spFactor.c fasthenry-3.0/src/fasthenry/sparse/spFactor.c
--- fasthenry-3.0.orig/src/fasthenry/sparse/spFactor.c	2004-08-03 17:17:18.000000000 -0500
+++ fasthenry-3.0/src/fasthenry/sparse/spFactor.c	2009-08-05 14:06:10.000000000 -0500
@@ -106,229 +106,59 @@
 
 
 
-
-
 
 /*
- *  ORDER AND FACTOR MATRIX
+ *  ZERO PIVOT ENCOUNTERED
  *
- *  This routine chooses a pivot order for the matrix and factors it
- *  into LU form.  It handles both the initial factorization and subsequent
- *  factorizations when a reordering is desired.  This is handled in a manner
- *  that is transparent to the user.  The routine uses a variation of
- *  Gauss's method where the pivots are associated with L and the
- *  diagonal terms of U are one.
+ *  This routine is called when a singular matrix is found.  It then
+ *  records the current row and column and exits.
  *
  *  >>> Returned:
- *  The error code is returned.  Possible errors are listed below.
+ *  The error code spSINGULAR or spZERO_DIAG is returned.
  *
  *  >>> Arguments:
- *  Matrix  <input>  (char *)
+ *  Matrix  <input>  (MatrixPtr)
  *      Pointer to matrix.
- *  RHS  <input>  (RealVector)
- *      Representative right-hand side vector that is used to determine
- *      pivoting order when the right hand side vector is sparse.  If
- *      RHS is a NULL pointer then the RHS vector is assumed to
- *      be full and it is not used when determining the pivoting
- *      order.
- *  RelThreshold  <input>  (RealNumber)
- *      This number determines what the pivot relative threshold will
- *      be.  It should be between zero and one.  If it is one then the
- *      pivoting method becomes complete pivoting, which is very slow
- *      and tends to fill up the matrix.  If it is set close to zero
- *      the pivoting method becomes strict Markowitz with no
- *      threshold.  The pivot threshold is used to eliminate pivot
- *      candidates that would cause excessive element growth if they
- *      were used.  Element growth is the cause of roundoff error.
- *      Element growth occurs even in well-conditioned matrices.
- *      Setting the RelThreshold large will reduce element growth and
- *      roundoff error, but setting it too large will cause execution
- *      time to be excessive and will result in a large number of
- *      fill-ins.  If this occurs, accuracy can actually be degraded
- *      because of the large number of operations required on the
- *      matrix due to the large number of fill-ins.  A good value seems
- *      to be 0.001.  The default is chosen by giving a value larger
- *      than one or less than or equal to zero.  This value should be
- *      increased and the matrix resolved if growth is found to be
- *      excessive.  Changing the pivot threshold does not improve
- *      performance on matrices where growth is low, as is often the
- *      case with ill-conditioned matrices.  Once a valid threshold is
- *      given, it becomes the new default.  The default value of
- *      RelThreshold was choosen for use with nearly diagonally
- *      dominant matrices such as node- and modified-node admittance
- *      matrices.  For these matrices it is usually best to use
- *      diagonal pivoting.  For matrices without a strong diagonal, it
- *      is usually best to use a larger threshold, such as 0.01 or
- *      0.1.
- *  AbsThreshold  <input>  (RealNumber)
- *      The absolute magnitude an element must have to be considered
- *      as a pivot candidate, except as a last resort.  This number
- *      should be set significantly smaller than the smallest diagonal
- *      element that is is expected to be placed in the matrix.  If
- *      there is no reasonable prediction for the lower bound on these
- *      elements, then AbsThreshold should be set to zero.
- *      AbsThreshold is used to reduce the possibility of choosing as a
- *      pivot an element that has suffered heavy cancellation and as a
- *      result mainly consists of roundoff error.  Once a valid
- *      threshold is given, it becomes the new default.
- *  DiagPivoting  <input>  (BOOLEAN)
- *      A flag indicating that pivot selection should be confined to the
- *      diagonal if possible.  If DiagPivoting is nonzero and if
- *      DIAGONAL_PIVOTING is enabled pivots will be chosen only from
- *      the diagonal unless there are no diagonal elements that satisfy
- *      the threshold criteria.  Otherwise, the entire reduced
- *      submatrix is searched when looking for a pivot.  The diagonal
- *      pivoting in Sparse is efficient and well refined, while the
- *      off-diagonal pivoting is not.  For symmetric and near symmetric
- *      matrices, it is best to use diagonal pivoting because it
- *      results in the best performance when reordering the matrix and
- *      when factoring the matrix without ordering.  If there is a
- *      considerable amount of nonsymmetry in the matrix, then
- *      off-diagonal pivoting may result in a better equation ordering
- *      simply because there are more pivot candidates to choose from.
- *      A better ordering results in faster subsequent factorizations.
- *      However, the initial pivot selection process takes considerably
- *      longer for off-diagonal pivoting.
- *
- *  >>> Local variables:
- *  pPivot  (ElementPtr)
- *      Pointer to the element being used as a pivot.
- *  ReorderingRequired  (BOOLEAN)
- *      Flag that indicates whether reordering is required.
- *
- *  >>> Possible errors:
- *  spNO_MEMORY
- *  spSINGULAR
- *  spSMALL_PIVOT
- *  Error is cleared in this function.
+ *  Step  <input>  (int)
+ *      Index of diagonal that is zero.
  */
 
-int
-spOrderAndFactor( eMatrix, RHS, RelThreshold, AbsThreshold, DiagPivoting )
+static int
+MatrixIsSingular( Matrix, Step )
 
-char *eMatrix;
-RealNumber  RHS[], RelThreshold, AbsThreshold;
-BOOLEAN DiagPivoting;
+MatrixPtr  Matrix;
+int  Step;
 {
-MatrixPtr  Matrix = (MatrixPtr)eMatrix;
-ElementPtr  pPivot;
-int  Step, Size, ReorderingRequired;
-ElementPtr SearchForPivot();
-RealNumber LargestInCol, FindLargestInCol();
-
-/* Begin `spOrderAndFactor'. */
-    ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored);
-
-    Matrix->Error = spOKAY;
-    Size = Matrix->Size;
-    if (RelThreshold <= 0.0) RelThreshold = Matrix->RelThreshold;
-    if (RelThreshold > 1.0) RelThreshold = Matrix->RelThreshold;
-    Matrix->RelThreshold = RelThreshold;
-    if (AbsThreshold < 0.0) AbsThreshold = Matrix->AbsThreshold;
-    Matrix->AbsThreshold = AbsThreshold;
-    ReorderingRequired = NO;
-
-    if (NOT Matrix->NeedsOrdering)
-    {
-/* Matrix has been factored before and reordering is not required. */
-        for (Step = 1; Step <= Size; Step++)
-        {   pPivot = Matrix->Diag[Step];
-            LargestInCol = FindLargestInCol(pPivot->NextInCol);
-            if ((LargestInCol * RelThreshold < ELEMENT_MAG(pPivot)))
-            {   if (Matrix->Complex)
-                    ComplexRowColElimination( Matrix, pPivot );
-                else
-                    RealRowColElimination( Matrix, pPivot );
-            }
-            else
-            {   ReorderingRequired = YES;
-                break; /* for loop */
-            }
-        }
-        if (NOT ReorderingRequired)
-            goto Done;
-        else
-        {
-/*
- * A pivot was not large enough to maintain accuracy,
- * so a partial reordering is required.
- */
-
-#if ANNOTATE >= ON_STRANGE_BEHAVIOR
-            printf("Reordering,  Step = %1d\n", Step);
-#endif
-        }
-    } /* End of if(NOT Matrix->NeedsOrdering) */
-    else
-    {
-/*
- * This is the first time the matrix has been factored.  These few statements
- * indicate to the rest of the code that a full reodering is required rather
- * than a partial reordering, which occurs during a failure of a fast
- * factorization.
- */
-        Step = 1;
-        if (NOT Matrix->RowsLinked)
-            spcLinkRows( Matrix );
-        if (NOT Matrix->InternalVectorsAllocated)
-            CreateInternalVectors( Matrix );
-        if (Matrix->Error >= spFATAL)
-            return Matrix->Error;
-    }
-
-/* Form initial Markowitz products. */
-    CountMarkowitz( Matrix, RHS, Step );
-    MarkowitzProducts( Matrix, Step );
-    Matrix->MaxRowCountInLowerTri = -1;
-
-/* Perform reordering and factorization. */
-    for (; Step <= Size; Step++)
-    {   pPivot = SearchForPivot( Matrix, Step, DiagPivoting );
-        if (pPivot == NULL) return MatrixIsSingular( Matrix, Step );
-        ExchangeRowsAndCols( Matrix, pPivot, Step );
+/* Begin `MatrixIsSingular'. */
 
-        if (Matrix->Complex)
-            ComplexRowColElimination( Matrix, pPivot );
-        else
-            RealRowColElimination( Matrix, pPivot );
+    Matrix->SingularRow = Matrix->IntToExtRowMap[ Step ];
+    Matrix->SingularCol = Matrix->IntToExtColMap[ Step ];
+    return (Matrix->Error = spSINGULAR);
+}
 
-        if (Matrix->Error >= spFATAL) return Matrix->Error;
-        UpdateMarkowitzNumbers( Matrix, pPivot );
 
-#if ANNOTATE == FULL
-        WriteStatus( Matrix, Step );
-#endif
-    }
+static int
+ZeroPivot( Matrix, Step )
 
-Done:
-    Matrix->NeedsOrdering = NO;
-    Matrix->Reordered = YES;
-    Matrix->Factored = YES;
+MatrixPtr  Matrix;
+int  Step;
+{
+/* Begin `ZeroPivot'. */
 
-    return Matrix->Error;
+    Matrix->SingularRow = Matrix->IntToExtRowMap[ Step ];
+    Matrix->SingularCol = Matrix->IntToExtColMap[ Step ];
+    return (Matrix->Error = spZERO_DIAG);
 }
 
 
 
-
-
-
 
+#if spCOMPLEX
 /*
- *  FACTOR MATRIX
+ *  FACTOR COMPLEX MATRIX
  *
- *  This routine is the companion routine to spOrderAndFactor().
- *  Unlike spOrderAndFactor(), spFactor() cannot change the ordering.
- *  It is also faster than spOrderAndFactor().  The standard way of
- *  using these two routines is to first use spOrderAndFactor() for the
- *  initial factorization.  For subsequent factorizations, spFactor()
- *  is used if there is some assurance that little growth will occur
- *  (say for example, that the matrix is diagonally dominant).  If
- *  spFactor() is called for the initial factorization of the matrix,
- *  then spOrderAndFactor() is automatically called with the default
- *  threshold.  This routine uses "row at a time" LU factorization.
- *  Pivots are associated with the lower triangular matrix and the
- *  diagonals of the upper triangular matrix are ones.
+ *  This routine is the companion routine to spFactor(), it
+ *  handles complex matrices.  It is otherwise identical.
  *
  *  >>> Returned:
  *  The error code is returned.  Possible errors are listed below.
@@ -338,145 +168,22 @@
  *      Pointer to matrix.
  *
  *  >>> Possible errors:
- *  spNO_MEMORY
  *  spSINGULAR
- *  spZERO_DIAG
- *  spSMALL_PIVOT
  *  Error is cleared in this function.
  */
 
-int
-spFactor( eMatrix )
+static int
+FactorComplexMatrix( Matrix )
 
-char *eMatrix;
+MatrixPtr  Matrix;
 {
-MatrixPtr  Matrix = (MatrixPtr)eMatrix;
 register  ElementPtr  pElement;
 register  ElementPtr  pColumn;
 register  int  Step, Size;
-RealNumber Mult;
-
-/* Begin `spFactor'. */
-    ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored);
+ComplexNumber Mult, Pivot;
 
-    if (Matrix->NeedsOrdering)
-    {   return spOrderAndFactor( eMatrix, (RealVector)NULL,
-                                 0.0, 0.0, DIAG_PIVOTING_AS_DEFAULT );
-    }
-    if (NOT Matrix->Partitioned) spPartition( eMatrix, spDEFAULT_PARTITION );
-#if spCOMPLEX
-    if (Matrix->Complex) return FactorComplexMatrix( Matrix );
-#endif
-
-#if REAL
-    Size = Matrix->Size;
-
-    if (Matrix->Diag[1]->Real == 0.0) return ZeroPivot( Matrix, 1 );
-    Matrix->Diag[1]->Real = 1.0 / Matrix->Diag[1]->Real;
-
-/* Start factorization. */
-    for (Step = 2; Step <= Size; Step++)
-    {   if (Matrix->DoRealDirect[Step])
-        {   /* Update column using direct addressing scatter-gather. */
-            register RealNumber *Dest = (RealNumber *)Matrix->Intermediate;
-
-/* Scatter. */
-            pElement = Matrix->FirstInCol[Step];
-            while (pElement != NULL)
-            {   Dest[pElement->Row] = pElement->Real;
-                pElement = pElement->NextInCol;
-            }
-
-/* Update column. */
-            pColumn = Matrix->FirstInCol[Step];
-            while (pColumn->Row < Step)
-            {   pElement = Matrix->Diag[pColumn->Row];
-                pColumn->Real = Dest[pColumn->Row] * pElement->Real;
-                while ((pElement = pElement->NextInCol) != NULL)
-                    Dest[pElement->Row] -= pColumn->Real * pElement->Real;
-                pColumn = pColumn->NextInCol;
-            }
-
-/* Gather. */
-            pElement = Matrix->Diag[Step]->NextInCol;
-            while (pElement != NULL)
-            {   pElement->Real = Dest[pElement->Row];
-                pElement = pElement->NextInCol;
-            }
-
-/* Check for singular matrix. */
-            if (Dest[Step] == 0.0) return ZeroPivot( Matrix, Step );
-            Matrix->Diag[Step]->Real = 1.0 / Dest[Step];
-        }
-        else
-        {   /* Update column using indirect addressing scatter-gather. */
-            register RealNumber **pDest = (RealNumber **)Matrix->Intermediate;
-
-/* Scatter. */
-            pElement = Matrix->FirstInCol[Step];
-            while (pElement != NULL)
-            {   pDest[pElement->Row] = &pElement->Real;
-                pElement = pElement->NextInCol;
-            }
-
-/* Update column. */
-            pColumn = Matrix->FirstInCol[Step];
-            while (pColumn->Row < Step)
-            {   pElement = Matrix->Diag[pColumn->Row];
-                Mult = (*pDest[pColumn->Row] *= pElement->Real);
-                while ((pElement = pElement->NextInCol) != NULL)
-                    *pDest[pElement->Row] -= Mult * pElement->Real;
-                pColumn = pColumn->NextInCol;
-            }
-
-/* Check for singular matrix. */
-            if (Matrix->Diag[Step]->Real == 0.0)
-                return ZeroPivot( Matrix, Step );
-            Matrix->Diag[Step]->Real = 1.0 / Matrix->Diag[Step]->Real;
-        }
-    }
-
-    Matrix->Factored = YES;
-    return (Matrix->Error = spOKAY);
-#endif /* REAL */
-}
-
-
-
-
-
-
-#if spCOMPLEX
-/*
- *  FACTOR COMPLEX MATRIX
- *
- *  This routine is the companion routine to spFactor(), it
- *  handles complex matrices.  It is otherwise identical.
- *
- *  >>> Returned:
- *  The error code is returned.  Possible errors are listed below.
- *
- *  >>> Arguments:
- *  Matrix  <input>  (char *)
- *      Pointer to matrix.
- *
- *  >>> Possible errors:
- *  spSINGULAR
- *  Error is cleared in this function.
- */
-
-static int
-FactorComplexMatrix( Matrix )
-
-MatrixPtr  Matrix;
-{
-register  ElementPtr  pElement;
-register  ElementPtr  pColumn;
-register  int  Step, Size;
-ComplexNumber Mult, Pivot;
-
-/* Begin `FactorComplexMatrix'. */
-    ASSERT(Matrix->Complex);
+/* Begin `FactorComplexMatrix'. */
+    ASSERT(Matrix->Complex);
 
     Size = Matrix->Size;
     pElement = Matrix->Diag[1];
@@ -565,256 +272,83 @@
 
 
 
-
 
 /*
- *  PARTITION MATRIX
- *
- *  This routine determines the cost to factor each row using both
- *  direct and indirect addressing and decides, on a row-by-row basis,
- *  which addressing mode is fastest.  This information is used in
- *  spFactor() to speed the factorization.
+ *  CREATE INTERNAL VECTORS
  *
- *  When factoring a previously ordered matrix using spFactor(), Sparse
- *  operates on a row-at-a-time basis.  For speed, on each step, the
- *  row being updated is copied into a full vector and the operations
- *  are performed on that vector.  This can be done one of two ways,
- *  either using direct addressing or indirect addressing.  Direct
- *  addressing is fastest when the matrix is relatively dense and
- *  indirect addressing is best when the matrix is quite sparse.  The
- *  user selects the type of partition used with Mode.  If Mode is set
- *  to spDIRECT_PARTITION, then the all rows are placed in the direct
- *  addressing partition.  Similarly, if Mode is set to
- *  spINDIRECT_PARTITION, then the all rows are placed in the indirect
- *  addressing partition.  By setting Mode to spAUTO_PARTITION, the
- *  user allows Sparse to select the partition for each row
- *  individually.  spFactor() generally runs faster if Sparse is
- *  allowed to choose its own partitioning, however choosing a
- *  partition is expensive.  The time required to choose a partition is
- *  of the same order of the cost to factor the matrix.  If you plan to
- *  factor a large number of matrices with the same structure, it is
- *  best to let Sparse choose the partition.  Otherwise, you should
- *  choose the partition based on the predicted density of the matrix.
+ *  Creates the Markowitz and Intermediate vectors.
  *
  *  >>> Arguments:
- *  Matrix  <input>  (char *)
+ *  Matrix  <input>  (MatrixPtr)
  *      Pointer to matrix.
- *  Mode  <input>  (int)
- *      Mode must be one of three special codes: spDIRECT_PARTITION,
- *      spINDIRECT_PARTITION, or spAUTO_PARTITION.
+ *
+ *  >>> Local variables:
+ *  SizePlusOne  (unsigned)
+ *      Size of the arrays to be allocated.
+ *
+ *  >>> Possible errors:
+ *  spNO_MEMORY
  */
 
-void
-spPartition( eMatrix, Mode )
+static
+CreateInternalVectors( Matrix )
 
-char *eMatrix;
-int Mode;
+MatrixPtr  Matrix;
 {
-MatrixPtr  Matrix = (MatrixPtr)eMatrix;
-register  ElementPtr  pElement, pColumn;
-register  int  Step, Size;
-register  int  *Nc, *No, *Nm;
-BOOLEAN *DoRealDirect, *DoCmplxDirect;
+int  Size;
 
-/* Begin `spPartition'. */
-    ASSERT( IS_SPARSE( Matrix ) );
-    if (Matrix->Partitioned) return;
-    Size = Matrix->Size;
-    DoRealDirect = Matrix->DoRealDirect;
-    DoCmplxDirect = Matrix->DoCmplxDirect;
-    Matrix->Partitioned = YES;
+/* Begin `CreateInternalVectors'. */
+/* Create Markowitz arrays. */
+    Size= Matrix->Size;
 
-/* If partition is specified by the user, this is easy. */
-    if (Mode == spDEFAULT_PARTITION) Mode = DEFAULT_PARTITION;
-    if (Mode == spDIRECT_PARTITION)
-    {   for (Step = 1; Step <= Size; Step++)
+    if (Matrix->MarkowitzRow == NULL)
+    {   if (( Matrix->MarkowitzRow = ALLOC(int, Size+1)) == NULL)
+            Matrix->Error = spNO_MEMORY;
+    }
+    if (Matrix->MarkowitzCol == NULL)
+    {   if (( Matrix->MarkowitzCol = ALLOC(int, Size+1)) == NULL)
+            Matrix->Error = spNO_MEMORY;
+    }
+    if (Matrix->MarkowitzProd == NULL)
+    {   if (( Matrix->MarkowitzProd = ALLOC(long, Size+2)) == NULL)
+            Matrix->Error = spNO_MEMORY;
+    }
+
+/* Create DoDirect vectors for use in spFactor(). */
 #if REAL
-            DoRealDirect[Step] = YES;
+    if (Matrix->DoRealDirect == NULL)
+    {   if (( Matrix->DoRealDirect = ALLOC(BOOLEAN, Size+1)) == NULL)
+            Matrix->Error = spNO_MEMORY;
+    }
 #endif
 #if spCOMPLEX
-            DoCmplxDirect[Step] = YES;
-#endif
-        return;
+    if (Matrix->DoCmplxDirect == NULL)
+    {   if (( Matrix->DoCmplxDirect = ALLOC(BOOLEAN, Size+1)) == NULL)
+            Matrix->Error = spNO_MEMORY;
     }
-    else if (Mode == spINDIRECT_PARTITION)
-    {   for (Step = 1; Step <= Size; Step++)
-#if REAL
-            DoRealDirect[Step] = NO;
 #endif
+
+/* Create Intermediate vectors for use in MatrixSolve. */
 #if spCOMPLEX
-            DoCmplxDirect[Step] = NO;
-#endif
-        return;
+    if (Matrix->Intermediate == NULL)
+    {   if ((Matrix->Intermediate = ALLOC(RealNumber,2*(Size+1))) == NULL)
+            Matrix->Error = spNO_MEMORY;
     }
-    else ASSERT( Mode == spAUTO_PARTITION );
-
-/* Otherwise, count all operations needed in when factoring matrix. */
-    Nc = (int *)Matrix->MarkowitzRow;
-    No = (int *)Matrix->MarkowitzCol;
-    Nm = (int *)Matrix->MarkowitzProd;
+#else
+    if (Matrix->Intermediate == NULL)
+    {   if ((Matrix->Intermediate = ALLOC(RealNumber, Size+1)) == NULL)
+            Matrix->Error = spNO_MEMORY;
+    }
+#endif
 
-/* Start mock-factorization. */
-    for (Step = 1; Step <= Size; Step++)
-    {   Nc[Step] = No[Step] = Nm[Step] = 0;
+    if (Matrix->Error != spNO_MEMORY)
+        Matrix->InternalVectorsAllocated = YES;
+    return;
+}
 
-        pElement = Matrix->FirstInCol[Step];
-        while (pElement != NULL)
-        {   Nc[Step]++;
-            pElement = pElement->NextInCol;
-        }
 
-        pColumn = Matrix->FirstInCol[Step];
-        while (pColumn->Row < Step)
-        {   pElement = Matrix->Diag[pColumn->Row];
-            Nm[Step]++;
-            while ((pElement = pElement->NextInCol) != NULL)
-                No[Step]++;
-            pColumn = pColumn->NextInCol;
-        }
-    }
 
-    for (Step = 1; Step <= Size; Step++)
-    {
-/*
- * The following are just estimates based on a count on the number of
- * machine instructions used on each machine to perform the various
- * tasks.  It was assumed that each machine instruction required the
- * same amount of time (I don't believe this is true for the VAX, and
- * have no idea if this is true for the 68000 family).  For optimum
- * performance, these numbers should be tuned to the machine.
- *   Nc is the number of nonzero elements in the column.
- *   Nm is the number of multipliers in the column.
- *   No is the number of operations in the inner loop.
- */
-
-#define generic
-#ifdef hp9000s300
-#if REAL
-        DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
-#endif
-#if spCOMPLEX
-        /* On the hp350, it is never profitable to use direct for complex. */
-        DoCmplxDirect[Step] = NO;
-#endif
-#undef generic
-#endif
-
-#ifdef vax
-#if REAL
-        DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
-#endif
-#if spCOMPLEX
-        DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7*Nc[Step] - 4*Nm[Step]);
-#endif
-#undef generic
-#endif
-
-#ifdef generic
-#if REAL
-        DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
-#endif
-#if spCOMPLEX
-        DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7*Nc[Step] - 4*Nm[Step]);
-#endif
-#undef generic
-#endif
-    }
-
-#if (ANNOTATE == FULL)
-    {   int Ops = 0;
-        for (Step = 1; Step <= Size; Step++)
-            Ops += No[Step];
-        printf("Operation count for inner loop of factorization = %d.\n", Ops);
-    }
-#endif
-    return;
-}
-
-
-
-
-
-
-
-/*
- *  CREATE INTERNAL VECTORS
- *
- *  Creates the Markowitz and Intermediate vectors.
- *
- *  >>> Arguments:
- *  Matrix  <input>  (MatrixPtr)
- *      Pointer to matrix.
- *
- *  >>> Local variables:
- *  SizePlusOne  (unsigned)
- *      Size of the arrays to be allocated.
- *
- *  >>> Possible errors:
- *  spNO_MEMORY
- */
-
-static
-CreateInternalVectors( Matrix )
-
-MatrixPtr  Matrix;
-{
-int  Size;
-
-/* Begin `CreateInternalVectors'. */
-/* Create Markowitz arrays. */
-    Size= Matrix->Size;
-
-    if (Matrix->MarkowitzRow == NULL)
-    {   if (( Matrix->MarkowitzRow = ALLOC(int, Size+1)) == NULL)
-            Matrix->Error = spNO_MEMORY;
-    }
-    if (Matrix->MarkowitzCol == NULL)
-    {   if (( Matrix->MarkowitzCol = ALLOC(int, Size+1)) == NULL)
-            Matrix->Error = spNO_MEMORY;
-    }
-    if (Matrix->MarkowitzProd == NULL)
-    {   if (( Matrix->MarkowitzProd = ALLOC(long, Size+2)) == NULL)
-            Matrix->Error = spNO_MEMORY;
-    }
-
-/* Create DoDirect vectors for use in spFactor(). */
-#if REAL
-    if (Matrix->DoRealDirect == NULL)
-    {   if (( Matrix->DoRealDirect = ALLOC(BOOLEAN, Size+1)) == NULL)
-            Matrix->Error = spNO_MEMORY;
-    }
-#endif
-#if spCOMPLEX
-    if (Matrix->DoCmplxDirect == NULL)
-    {   if (( Matrix->DoCmplxDirect = ALLOC(BOOLEAN, Size+1)) == NULL)
-            Matrix->Error = spNO_MEMORY;
-    }
-#endif
-
-/* Create Intermediate vectors for use in MatrixSolve. */
-#if spCOMPLEX
-    if (Matrix->Intermediate == NULL)
-    {   if ((Matrix->Intermediate = ALLOC(RealNumber,2*(Size+1))) == NULL)
-            Matrix->Error = spNO_MEMORY;
-    }
-#else
-    if (Matrix->Intermediate == NULL)
-    {   if ((Matrix->Intermediate = ALLOC(RealNumber, Size+1)) == NULL)
-            Matrix->Error = spNO_MEMORY;
-    }
-#endif
-
-    if (Matrix->Error != spNO_MEMORY)
-        Matrix->InternalVectorsAllocated = YES;
-    return;
-}
-
-
-
-
-
-
-
+
 /*
  *  COUNT MARKOWITZ
  *
@@ -918,11 +452,6 @@
 
 
 
-
-
-
-
-
 
 /*
  *  MARKOWITZ PRODUCTS
@@ -993,108 +522,83 @@
 
 
 
-
-
-
-
-
-
 
 /*
- *  SEARCH FOR BEST PIVOT
+ *  DETERMINE THE MAGNITUDE OF THE LARGEST ELEMENT IN A COLUMN
+ *  EXCLUDING AN ELEMENT
  *
- *  Performs a search to determine the element with the lowest Markowitz
- *  Product that is also acceptable.  An acceptable element is one that is
- *  larger than the AbsThreshold and at least as large as RelThreshold times
- *  the largest element in the same column.  The first step is to look for
- *  singletons if any exist.  If none are found, then all the diagonals are
- *  searched. The diagonal is searched once quickly using the assumption that
- *  elements on the diagonal are large compared to other elements in their
- *  column, and so the pivot can be chosen only on the basis of the Markowitz
- *  criterion.  After a element has been chosen to be pivot on the basis of
- *  its Markowitz product, it is checked to see if it is large enough.
- *  Waiting to the end of the Markowitz search to check the size of a pivot
- *  candidate saves considerable time, but is not guaranteed to find an
- *  acceptable pivot.  Thus if unsuccessful a second pass of the diagonal is
- *  made.  This second pass checks to see if an element is large enough during
- *  the search, not after it.  If still no acceptable pivot candidate has
- *  been found, the search expands to cover the entire matrix.
+ *  This routine searches a column and returns the magnitude of the largest
+ *  element.  One given element is specifically excluded from the search.
+ *
+ *  The search is conducted by starting at the first element in the column
+ *  and moving down the column until the active part of the matrix is entered,
+ *  i.e. the reduced submatrix.  The rest of the column is then traversed
+ *  looking for the largest element.
  *
  *  >>> Returned:
- *  A pointer to the element chosen to be pivot.  If every element in the
- *  matrix is zero, then NULL is returned.
+ *  The magnitude of the largest element in the active portion of the column,
+ *  excluding the specified element, is returned.
  *
  *  >>> Arguments:
- *  Matrix  <input>  (MatrixPtr)
- *      Pointer to matrix.
+ *  Matrix  <input>    (MatrixPtr)
+ *      Pointer to the matrix.
+ *  pElement  <input>  (ElementPtr)
+ *      The pointer to the element that is to be excluded from search. Column
+ *      to be searched is one that contains this element.  Also used to
+ *      access the elements in the column.
  *  Step  <input>  (int)
- *      The row and column number of the beginning of the reduced submatrix.
+ *      Index of the diagonal currently being eliminated.  Indicates where
+ *      the active part of the matrix begins.
  *
  *  >>> Local variables:
- *  ChosenPivot  (ElementPtr)
- *      Pointer to element that has been chosen to be the pivot.
- *
- *  >>> Possible errors:
- *  spSINGULAR
- *  spSMALL_PIVOT
+ *  Col  (int)
+ *      The number of the column to be searched.  Also the column number of
+ *      the element to be avoided in the search.
+ *  Largest  (RealNumber)
+ *      The magnitude of the largest element.
+ *  Magnitude  (RealNumber)
+ *      The magnitude of the currently active element.
+ *  Row  (int)
+ *      The row number of element to be excluded from the search.
  */
 
-static ElementPtr
-SearchForPivot( Matrix, Step, DiagPivoting )
+static RealNumber
+FindBiggestInColExclude( Matrix, pElement, Step )
 
 MatrixPtr Matrix;
-int Step, DiagPivoting;
+register  ElementPtr  pElement;
+register  int Step;
 {
-register ElementPtr  ChosenPivot;
-ElementPtr  SearchForSingleton();
-ElementPtr  QuicklySearchDiagonal();
-ElementPtr  SearchDiagonal();
-ElementPtr  SearchEntireMatrix();
+register  int  Row;
+int  Col;
+RealNumber  Largest, Magnitude;
 
-/* Begin `SearchForPivot'. */
+/* Begin `FindBiggestInColExclude'. */
+    Row = pElement->Row;
+    Col = pElement->Col;
+    pElement = Matrix->FirstInCol[Col];
 
-/* If singletons exist, look for an acceptable one to use as pivot. */
-    if (Matrix->Singletons)
-    {   ChosenPivot = SearchForSingleton( Matrix, Step );
-        if (ChosenPivot != NULL)
-        {   Matrix->PivotSelectionMethod = 's';
-            return ChosenPivot;
+/* Travel down column until reduced submatrix is entered. */
+    while ((pElement != NULL) AND (pElement->Row < Step))
+        pElement = pElement->NextInCol;
+
+/* Initialize the variable Largest. */
+    if (pElement->Row != Row)
+        Largest = ELEMENT_MAG(pElement);
+    else
+        Largest = 0.0;
+
+/* Search rest of column for largest element, avoiding excluded element. */
+    while ((pElement = pElement->NextInCol) != NULL)
+    {   if ((Magnitude = ELEMENT_MAG(pElement)) > Largest)
+        {   if (pElement->Row != Row)
+                Largest = Magnitude;
         }
     }
 
-#if DIAGONAL_PIVOTING
-    if (DiagPivoting)
-    {
-/*
- * Either no singletons exist or they weren't acceptable.  Take quick first
- * pass at searching diagonal.  First search for element on diagonal of 
- * remaining submatrix with smallest Markowitz product, then check to see
- * if it okay numerically.  If not, QuicklySearchDiagonal fails.
- */
-        ChosenPivot = QuicklySearchDiagonal( Matrix, Step );
-        if (ChosenPivot != NULL)
-        {   Matrix->PivotSelectionMethod = 'q';
-            return ChosenPivot;
-        }
+    return Largest;
+}
 
-/*
- * Quick search of diagonal failed, carefully search diagonal and check each
- * pivot candidate numerically before even tentatively accepting it.
- */
-        ChosenPivot = SearchDiagonal( Matrix, Step );
-        if (ChosenPivot != NULL)
-        {   Matrix->PivotSelectionMethod = 'd';
-            return ChosenPivot;
-        }
-    }
-#endif /* DIAGONAL_PIVOTING */
-
-/* No acceptable pivot found yet, search entire matrix. */
-    ChosenPivot = SearchEntireMatrix( Matrix, Step );
-    Matrix->PivotSelectionMethod = 'e';
-
-    return ChosenPivot;
-}
 
 
 
@@ -1280,6 +784,201 @@
 
 
 
+
+
+/*
+ *  PARTITION MATRIX
+ *
+ *  This routine determines the cost to factor each row using both
+ *  direct and indirect addressing and decides, on a row-by-row basis,
+ *  which addressing mode is fastest.  This information is used in
+ *  spFactor() to speed the factorization.
+ *
+ *  When factoring a previously ordered matrix using spFactor(), Sparse
+ *  operates on a row-at-a-time basis.  For speed, on each step, the
+ *  row being updated is copied into a full vector and the operations
+ *  are performed on that vector.  This can be done one of two ways,
+ *  either using direct addressing or indirect addressing.  Direct
+ *  addressing is fastest when the matrix is relatively dense and
+ *  indirect addressing is best when the matrix is quite sparse.  The
+ *  user selects the type of partition used with Mode.  If Mode is set
+ *  to spDIRECT_PARTITION, then the all rows are placed in the direct
+ *  addressing partition.  Similarly, if Mode is set to
+ *  spINDIRECT_PARTITION, then the all rows are placed in the indirect
+ *  addressing partition.  By setting Mode to spAUTO_PARTITION, the
+ *  user allows Sparse to select the partition for each row
+ *  individually.  spFactor() generally runs faster if Sparse is
+ *  allowed to choose its own partitioning, however choosing a
+ *  partition is expensive.  The time required to choose a partition is
+ *  of the same order of the cost to factor the matrix.  If you plan to
+ *  factor a large number of matrices with the same structure, it is
+ *  best to let Sparse choose the partition.  Otherwise, you should
+ *  choose the partition based on the predicted density of the matrix.
+ *
+ *  >>> Arguments:
+ *  Matrix  <input>  (char *)
+ *      Pointer to matrix.
+ *  Mode  <input>  (int)
+ *      Mode must be one of three special codes: spDIRECT_PARTITION,
+ *      spINDIRECT_PARTITION, or spAUTO_PARTITION.
+ */
+
+void
+spPartition( eMatrix, Mode )
+
+char *eMatrix;
+int Mode;
+{
+MatrixPtr  Matrix = (MatrixPtr)eMatrix;
+register  ElementPtr  pElement, pColumn;
+register  int  Step, Size;
+register  int  *Nc, *No, *Nm;
+BOOLEAN *DoRealDirect, *DoCmplxDirect;
+
+/* Begin `spPartition'. */
+    ASSERT( IS_SPARSE( Matrix ) );
+    if (Matrix->Partitioned) return;
+    Size = Matrix->Size;
+    DoRealDirect = Matrix->DoRealDirect;
+    DoCmplxDirect = Matrix->DoCmplxDirect;
+    Matrix->Partitioned = YES;
+
+/* If partition is specified by the user, this is easy. */
+    if (Mode == spDEFAULT_PARTITION) Mode = DEFAULT_PARTITION;
+    if (Mode == spDIRECT_PARTITION)
+    {   for (Step = 1; Step <= Size; Step++)
+#if REAL
+            DoRealDirect[Step] = YES;
+#endif
+#if spCOMPLEX
+            DoCmplxDirect[Step] = YES;
+#endif
+        return;
+    }
+    else if (Mode == spINDIRECT_PARTITION)
+    {   for (Step = 1; Step <= Size; Step++)
+#if REAL
+            DoRealDirect[Step] = NO;
+#endif
+#if spCOMPLEX
+            DoCmplxDirect[Step] = NO;
+#endif
+        return;
+    }
+    else ASSERT( Mode == spAUTO_PARTITION );
+
+/* Otherwise, count all operations needed in when factoring matrix. */
+    Nc = (int *)Matrix->MarkowitzRow;
+    No = (int *)Matrix->MarkowitzCol;
+    Nm = (int *)Matrix->MarkowitzProd;
+
+/* Start mock-factorization. */
+    for (Step = 1; Step <= Size; Step++)
+    {   Nc[Step] = No[Step] = Nm[Step] = 0;
+
+        pElement = Matrix->FirstInCol[Step];
+        while (pElement != NULL)
+        {   Nc[Step]++;
+            pElement = pElement->NextInCol;
+        }
+
+        pColumn = Matrix->FirstInCol[Step];
+        while (pColumn->Row < Step)
+        {   pElement = Matrix->Diag[pColumn->Row];
+            Nm[Step]++;
+            while ((pElement = pElement->NextInCol) != NULL)
+                No[Step]++;
+            pColumn = pColumn->NextInCol;
+        }
+    }
+
+    for (Step = 1; Step <= Size; Step++)
+    {
+/*
+ * The following are just estimates based on a count on the number of
+ * machine instructions used on each machine to perform the various
+ * tasks.  It was assumed that each machine instruction required the
+ * same amount of time (I don't believe this is true for the VAX, and
+ * have no idea if this is true for the 68000 family).  For optimum
+ * performance, these numbers should be tuned to the machine.
+ *   Nc is the number of nonzero elements in the column.
+ *   Nm is the number of multipliers in the column.
+ *   No is the number of operations in the inner loop.
+ */
+
+#define generic
+#ifdef hp9000s300
+#if REAL
+        DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
+#endif
+#if spCOMPLEX
+        /* On the hp350, it is never profitable to use direct for complex. */
+        DoCmplxDirect[Step] = NO;
+#endif
+#undef generic
+#endif
+
+#ifdef vax
+#if REAL
+        DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
+#endif
+#if spCOMPLEX
+        DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7*Nc[Step] - 4*Nm[Step]);
+#endif
+#undef generic
+#endif
+
+#ifdef generic
+#if REAL
+        DoRealDirect[Step] = (Nm[Step] + No[Step] > 3*Nc[Step] - 2*Nm[Step]);
+#endif
+#if spCOMPLEX
+        DoCmplxDirect[Step] = (Nm[Step] + No[Step] > 7*Nc[Step] - 4*Nm[Step]);
+#endif
+#undef generic
+#endif
+    }
+
+#if (ANNOTATE == FULL)
+    {   int Ops = 0;
+        for (Step = 1; Step <= Size; Step++)
+            Ops += No[Step];
+        printf("Operation count for inner loop of factorization = %d.\n", Ops);
+    }
+#endif
+    return;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
 #if DIAGONAL_PIVOTING
 #if MODIFIED_MARKOWITZ
@@ -1658,11 +1357,6 @@
 
 
 
-
-
-
-
-
 
 /*
  *  SEARCH DIAGONAL FOR PIVOT WITH MODIFIED MARKOWITZ CRITERION
@@ -1781,39 +1475,81 @@
 
 
 
-
-
-
-
-
-
 
 /*
- *  SEARCH ENTIRE MATRIX FOR BEST PIVOT
+ *  DETERMINE THE MAGNITUDE OF THE LARGEST ELEMENT IN A COLUMN
  *
- *  Performs a search over the entire matrix looking for the acceptable
- *  element with the lowest MarkowitzProduct.  If there are several that
- *  are tied for the smallest MarkowitzProduct, the tie is broken by using
- *  the ratio of the magnitude of the element being considered to the largest
- *  element in the same column.  If no element is acceptable then the largest
- *  element in the reduced submatrix is used as the pivot and the
- *  matrix is declared to be spSMALL_PIVOT.  If the largest element is
- *  zero, the matrix is declared to be spSINGULAR.
+ *  This routine searches a column and returns the magnitude of the largest
+ *  element.  This routine begins the search at the element pointed to by
+ *  pElement, the parameter.
+ *
+ *  The search is conducted by starting at the element specified by a pointer,
+ *  which should be one below the diagonal, and moving down the column.  On
+ *  the way down the column, the magnitudes of the elements are tested to see
+ *  if they are the largest yet found.
  *
  *  >>> Returned:
- *  A pointer to the diagonal element chosen to be pivot.  If no element is
- *  found, then NULL is returned and the matrix is spSINGULAR.
+ *  The magnitude of the largest element in the column below and including
+ *  the one pointed to by the input parameter.
  *
  *  >>> Arguments:
- *  Matrix  <input>  (MatrixPtr)
- *      Pointer to matrix.
- *  Step  <input>  (int)
- *      Index of the diagonal currently being eliminated.
+ *  pElement  <input>  (ElementPtr)
+ *      The pointer to the first element to be tested.  Also, used by the
+ *      routine to access all lower elements in the column.
  *
  *  >>> Local variables:
- *  ChosenPivot  (ElementPtr)
- *      Pointer to the element that has been chosen to be the pivot.
- *  LargestElementMag  (RealNumber)
+ *  Largest  (RealNumber)
+ *      The magnitude of the largest element.
+ *  Magnitude  (RealNumber)
+ *      The magnitude of the currently active element.
+ */
+
+static RealNumber
+FindLargestInCol( pElement )
+
+register  ElementPtr  pElement;
+{
+RealNumber  Magnitude, Largest = 0.0;
+
+/* Begin `FindLargestInCol'. */
+/* Search column for largest element beginning at Element. */
+    while (pElement != NULL)
+    {   if ((Magnitude = ELEMENT_MAG(pElement)) > Largest)
+            Largest = Magnitude;
+        pElement = pElement->NextInCol;
+    }
+
+    return Largest;
+}
+
+
+
+/*
+ *  SEARCH ENTIRE MATRIX FOR BEST PIVOT
+ *
+ *  Performs a search over the entire matrix looking for the acceptable
+ *  element with the lowest MarkowitzProduct.  If there are several that
+ *  are tied for the smallest MarkowitzProduct, the tie is broken by using
+ *  the ratio of the magnitude of the element being considered to the largest
+ *  element in the same column.  If no element is acceptable then the largest
+ *  element in the reduced submatrix is used as the pivot and the
+ *  matrix is declared to be spSMALL_PIVOT.  If the largest element is
+ *  zero, the matrix is declared to be spSINGULAR.
+ *
+ *  >>> Returned:
+ *  A pointer to the diagonal element chosen to be pivot.  If no element is
+ *  found, then NULL is returned and the matrix is spSINGULAR.
+ *
+ *  >>> Arguments:
+ *  Matrix  <input>  (MatrixPtr)
+ *      Pointer to matrix.
+ *  Step  <input>  (int)
+ *      Index of the diagonal currently being eliminated.
+ *
+ *  >>> Local variables:
+ *  ChosenPivot  (ElementPtr)
+ *      Pointer to the element that has been chosen to be the pivot.
+ *  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.
@@ -1931,150 +1667,102 @@
 
 
 
-
-
-
-
-
-
-
-/*
- *  DETERMINE THE MAGNITUDE OF THE LARGEST ELEMENT IN A COLUMN
- *
- *  This routine searches a column and returns the magnitude of the largest
- *  element.  This routine begins the search at the element pointed to by
- *  pElement, the parameter.
- *
- *  The search is conducted by starting at the element specified by a pointer,
- *  which should be one below the diagonal, and moving down the column.  On
- *  the way down the column, the magnitudes of the elements are tested to see
- *  if they are the largest yet found.
- *
- *  >>> Returned:
- *  The magnitude of the largest element in the column below and including
- *  the one pointed to by the input parameter.
- *
- *  >>> Arguments:
- *  pElement  <input>  (ElementPtr)
- *      The pointer to the first element to be tested.  Also, used by the
- *      routine to access all lower elements in the column.
- *
- *  >>> Local variables:
- *  Largest  (RealNumber)
- *      The magnitude of the largest element.
- *  Magnitude  (RealNumber)
- *      The magnitude of the currently active element.
- */
-
-static RealNumber
-FindLargestInCol( pElement )
-
-register  ElementPtr  pElement;
-{
-RealNumber  Magnitude, Largest = 0.0;
-
-/* Begin `FindLargestInCol'. */
-/* Search column for largest element beginning at Element. */
-    while (pElement != NULL)
-    {   if ((Magnitude = ELEMENT_MAG(pElement)) > Largest)
-            Largest = Magnitude;
-        pElement = pElement->NextInCol;
-    }
-
-    return Largest;
-}
-
-
-
-
-
-
-
-
-
 
 /*
- *  DETERMINE THE MAGNITUDE OF THE LARGEST ELEMENT IN A COLUMN
- *  EXCLUDING AN ELEMENT
- *
- *  This routine searches a column and returns the magnitude of the largest
- *  element.  One given element is specifically excluded from the search.
+ *  SEARCH FOR BEST PIVOT
  *
- *  The search is conducted by starting at the first element in the column
- *  and moving down the column until the active part of the matrix is entered,
- *  i.e. the reduced submatrix.  The rest of the column is then traversed
- *  looking for the largest element.
+ *  Performs a search to determine the element with the lowest Markowitz
+ *  Product that is also acceptable.  An acceptable element is one that is
+ *  larger than the AbsThreshold and at least as large as RelThreshold times
+ *  the largest element in the same column.  The first step is to look for
+ *  singletons if any exist.  If none are found, then all the diagonals are
+ *  searched. The diagonal is searched once quickly using the assumption that
+ *  elements on the diagonal are large compared to other elements in their
+ *  column, and so the pivot can be chosen only on the basis of the Markowitz
+ *  criterion.  After a element has been chosen to be pivot on the basis of
+ *  its Markowitz product, it is checked to see if it is large enough.
+ *  Waiting to the end of the Markowitz search to check the size of a pivot
+ *  candidate saves considerable time, but is not guaranteed to find an
+ *  acceptable pivot.  Thus if unsuccessful a second pass of the diagonal is
+ *  made.  This second pass checks to see if an element is large enough during
+ *  the search, not after it.  If still no acceptable pivot candidate has
+ *  been found, the search expands to cover the entire matrix.
  *
  *  >>> Returned:
- *  The magnitude of the largest element in the active portion of the column,
- *  excluding the specified element, is returned.
+ *  A pointer to the element chosen to be pivot.  If every element in the
+ *  matrix is zero, then NULL is returned.
  *
  *  >>> Arguments:
- *  Matrix  <input>    (MatrixPtr)
- *      Pointer to the matrix.
- *  pElement  <input>  (ElementPtr)
- *      The pointer to the element that is to be excluded from search. Column
- *      to be searched is one that contains this element.  Also used to
- *      access the elements in the column.
+ *  Matrix  <input>  (MatrixPtr)
+ *      Pointer to matrix.
  *  Step  <input>  (int)
- *      Index of the diagonal currently being eliminated.  Indicates where
- *      the active part of the matrix begins.
+ *      The row and column number of the beginning of the reduced submatrix.
  *
  *  >>> Local variables:
- *  Col  (int)
- *      The number of the column to be searched.  Also the column number of
- *      the element to be avoided in the search.
- *  Largest  (RealNumber)
- *      The magnitude of the largest element.
- *  Magnitude  (RealNumber)
- *      The magnitude of the currently active element.
- *  Row  (int)
- *      The row number of element to be excluded from the search.
+ *  ChosenPivot  (ElementPtr)
+ *      Pointer to element that has been chosen to be the pivot.
+ *
+ *  >>> Possible errors:
+ *  spSINGULAR
+ *  spSMALL_PIVOT
  */
 
-static RealNumber
-FindBiggestInColExclude( Matrix, pElement, Step )
+static ElementPtr
+SearchForPivot( Matrix, Step, DiagPivoting )
 
 MatrixPtr Matrix;
-register  ElementPtr  pElement;
-register  int Step;
+int Step, DiagPivoting;
 {
-register  int  Row;
-int  Col;
-RealNumber  Largest, Magnitude;
-
-/* Begin `FindBiggestInColExclude'. */
-    Row = pElement->Row;
-    Col = pElement->Col;
-    pElement = Matrix->FirstInCol[Col];
-
-/* Travel down column until reduced submatrix is entered. */
-    while ((pElement != NULL) AND (pElement->Row < Step))
-        pElement = pElement->NextInCol;
+register ElementPtr  ChosenPivot;
+ElementPtr  SearchForSingleton();
+ElementPtr  QuicklySearchDiagonal();
+ElementPtr  SearchDiagonal();
+ElementPtr  SearchEntireMatrix();
 
-/* Initialize the variable Largest. */
-    if (pElement->Row != Row)
-        Largest = ELEMENT_MAG(pElement);
-    else
-        Largest = 0.0;
+/* Begin `SearchForPivot'. */
 
-/* Search rest of column for largest element, avoiding excluded element. */
-    while ((pElement = pElement->NextInCol) != NULL)
-    {   if ((Magnitude = ELEMENT_MAG(pElement)) > Largest)
-        {   if (pElement->Row != Row)
-                Largest = Magnitude;
+/* If singletons exist, look for an acceptable one to use as pivot. */
+    if (Matrix->Singletons)
+    {   ChosenPivot = SearchForSingleton( Matrix, Step );
+        if (ChosenPivot != NULL)
+        {   Matrix->PivotSelectionMethod = 's';
+            return ChosenPivot;
         }
     }
 
-    return Largest;
-}
-
-
-
-
+#if DIAGONAL_PIVOTING
+    if (DiagPivoting)
+    {
+/*
+ * Either no singletons exist or they weren't acceptable.  Take quick first
+ * pass at searching diagonal.  First search for element on diagonal of 
+ * remaining submatrix with smallest Markowitz product, then check to see
+ * if it okay numerically.  If not, QuicklySearchDiagonal fails.
+ */
+        ChosenPivot = QuicklySearchDiagonal( Matrix, Step );
+        if (ChosenPivot != NULL)
+        {   Matrix->PivotSelectionMethod = 'q';
+            return ChosenPivot;
+        }
 
+/*
+ * Quick search of diagonal failed, carefully search diagonal and check each
+ * pivot candidate numerically before even tentatively accepting it.
+ */
+        ChosenPivot = SearchDiagonal( Matrix, Step );
+        if (ChosenPivot != NULL)
+        {   Matrix->PivotSelectionMethod = 'd';
+            return ChosenPivot;
+        }
+    }
+#endif /* DIAGONAL_PIVOTING */
 
+/* No acceptable pivot found yet, search entire matrix. */
+    ChosenPivot = SearchEntireMatrix( Matrix, Step );
+    Matrix->PivotSelectionMethod = 'e';
 
+    return ChosenPivot;
+}
 
 
 
@@ -2204,103 +1892,73 @@
 
 
 
-
-
-
-
-
 
 /*
- *  EXCHANGE ROWS
+ *  CREATE FILL-IN
  *
- *  Performs all required operations to exchange two rows. Those operations
- *  include: swap FirstInRow pointers, fixing up the NextInCol pointers,
- *  swapping row indexes in MatrixElements, and swapping Markowitz row
- *  counts.
+ *  This routine is used to create fill-ins and splice them into the
+ *  matrix.
+ *
+ *  >>> Returns:
+ *  Pointer to fill-in.
  *
  *  >>> Arguments:
  *  Matrix  <input>  (MatrixPtr)
  *      Pointer to the matrix.
- *  Row1  <input>  (int)
- *      Row index of one of the rows, becomes the smallest index.
- *  Row2  <input>  (int)
- *      Row index of the other row, becomes the largest index.
+ *  Col  <input>  (int)
+ *      Column index for element.
+ *  Row  <input>  (int)
+ *      Row index for element.
  *
- *  Local variables:
- *  Column  (int)
- *      Column in which row elements are currently being exchanged.
- *  Row1Ptr  (ElementPtr)
- *      Pointer to an element in Row1.
- *  Row2Ptr  (ElementPtr)
- *      Pointer to an element in Row2.
- *  Element1  (ElementPtr)
- *      Pointer to the element in Row1 to be exchanged.
- *  Element2  (ElementPtr)
- *      Pointer to the element in Row2 to be exchanged.
+ *  >>> Local variables:
+ *  pElement  (ElementPtr)
+ *      Pointer to an element in the matrix.
+ *  ppElementAbove  (ElementPtr *)
+ *      This contains the address of the pointer to the element just above the
+ *      one being created. It is used to speed the search and it is updated
+ *      with address of the created element.
+ *
+ *  >>> Possible errors:
+ *  spNO_MEMORY
  */
 
-spcRowExchange( Matrix, Row1, Row2 )
+static ElementPtr
+CreateFillin( Matrix, Row, Col )
 
 MatrixPtr Matrix;
-int  Row1, Row2;
+register int  Row;
+int  Col;
 {
-register  ElementPtr  Row1Ptr, Row2Ptr;
-int  Column;
-ElementPtr  Element1, Element2;
+register  ElementPtr  pElement, *ppElementAbove;
+ElementPtr  spcCreateElement();
 
-/* Begin `spcRowExchange'. */
-    if (Row1 > Row2)  SWAP(int, Row1, Row2);
+/* Begin `CreateFillin'. */
 
-    Row1Ptr = Matrix->FirstInRow[Row1];
-    Row2Ptr = Matrix->FirstInRow[Row2];
-    while (Row1Ptr != NULL OR Row2Ptr != NULL)
-    {
-/* Exchange elements in rows while traveling from left to right. */
-        if (Row1Ptr == NULL)
-        {   Column = Row2Ptr->Col;
-            Element1 = NULL;
-            Element2 = Row2Ptr;
-            Row2Ptr = Row2Ptr->NextInRow;
-        }
-        else if (Row2Ptr == NULL)
-        {   Column = Row1Ptr->Col;
-            Element1 = Row1Ptr;
-            Element2 = NULL;
-            Row1Ptr = Row1Ptr->NextInRow;
-        }
-        else if (Row1Ptr->Col < Row2Ptr->Col)
-        {   Column = Row1Ptr->Col;
-            Element1 = Row1Ptr;
-            Element2 = NULL;
-            Row1Ptr = Row1Ptr->NextInRow;
-        }
-        else if (Row1Ptr->Col > Row2Ptr->Col)
-        {   Column = Row2Ptr->Col;
-            Element1 = NULL;
-            Element2 = Row2Ptr;
-            Row2Ptr = Row2Ptr->NextInRow;
-        }
-        else   /* Row1Ptr->Col == Row2Ptr->Col */
-        {   Column = Row1Ptr->Col;
-            Element1 = Row1Ptr;
-            Element2 = Row2Ptr;
-            Row1Ptr = Row1Ptr->NextInRow;
-            Row2Ptr = Row2Ptr->NextInRow;
+/* Find Element above fill-in. */
+    ppElementAbove = &Matrix->FirstInCol[Col];
+    pElement = *ppElementAbove;
+    while (pElement != NULL)
+    {   if (pElement->Row < Row)
+        {   ppElementAbove = &pElement->NextInCol;
+            pElement = *ppElementAbove;
         }
+        else break; /* while loop */
+    }
 
-        ExchangeColElements( Matrix, Row1, Element1, Row2, Element2, Column);
-    }  /* end of while(Row1Ptr != NULL OR Row2Ptr != NULL) */
+/* End of search, create the element. */
+    pElement = spcCreateElement( Matrix, Row, Col, ppElementAbove, YES );
 
-    if (Matrix->InternalVectorsAllocated)
-        SWAP( int, Matrix->MarkowitzRow[Row1], Matrix->MarkowitzRow[Row2]);
-    SWAP( ElementPtr, Matrix->FirstInRow[Row1], Matrix->FirstInRow[Row2]);
-    SWAP( int, Matrix->IntToExtRowMap[Row1], Matrix->IntToExtRowMap[Row2]);
-#if TRANSLATE
-    Matrix->ExtToIntRowMap[ Matrix->IntToExtRowMap[Row1] ] = Row1;
-    Matrix->ExtToIntRowMap[ Matrix->IntToExtRowMap[Row2] ] = Row2;
-#endif
+/* Update Markowitz counts and products. */
+    Matrix->MarkowitzProd[Row] = ++Matrix->MarkowitzRow[Row] *
+                                   Matrix->MarkowitzCol[Row];
+    if ((Matrix->MarkowitzRow[Row] == 1) AND (Matrix->MarkowitzCol[Row] != 0))
+        Matrix->Singletons--;
+    Matrix->MarkowitzProd[Col] = ++Matrix->MarkowitzCol[Col] *
+                                   Matrix->MarkowitzRow[Col];
+    if ((Matrix->MarkowitzRow[Col] != 0) AND (Matrix->MarkowitzCol[Col] == 1))
+        Matrix->Singletons--;
 
-    return;
+    return pElement;
 }
 
 
@@ -2312,96 +1970,84 @@
 
 
 /*
- *  EXCHANGE COLUMNS
+ *  PERFORM ROW AND COLUMN ELIMINATION ON REAL MATRIX
  *
- *  Performs all required operations to exchange two columns. Those operations
- *  include: swap FirstInCol pointers, fixing up the NextInRow pointers,
- *  swapping column indexes in MatrixElements, and swapping Markowitz 
- *  column counts.
+ *  Eliminates a single row and column of the matrix and leaves single row of
+ *  the upper triangular matrix and a single column of the lower triangular
+ *  matrix in its wake.  Uses Gauss's method.
  *
- *  >>> Arguments:
+ *  >>> Argument:
  *  Matrix  <input>  (MatrixPtr)
  *      Pointer to the matrix.
- *  Col1  <input>  (int)
- *      Column index of one of the columns, becomes the smallest index.
- *  Col2  <input>  (int)
- *      Column index of the other column, becomes the largest index
+ *  pPivot  <input>  (ElementPtr)
+ *      Pointer to the current pivot.
  *
- *  Local variables:
+ *  >>> Local variables:
+ *  pLower  (ElementPtr)
+ *      Points to matrix element in lower triangular column.
+ *  pSub (ElementPtr)
+ *      Points to elements in the reduced submatrix.
  *  Row  (int)
- *      Row in which column elements are currently being exchanged.
- *  Col1Ptr  (ElementPtr)
- *      Pointer to an element in Col1.
- *  Col2Ptr  (ElementPtr)
- *      Pointer to an element in Col2.
- *  Element1  (ElementPtr)
- *      Pointer to the element in Col1 to be exchanged.
- *  Element2  (ElementPtr)
- *      Pointer to the element in Col2 to be exchanged.
+ *      Row index.
+ *  pUpper  (ElementPtr)
+ *      Points to matrix element in upper triangular row.
+ *
+ *  >>> Possible errors:
+ *  spNO_MEMORY
  */
 
-spcColExchange( Matrix, Col1, Col2 )
+static
+RealRowColElimination( Matrix, pPivot )
 
 MatrixPtr Matrix;
-int  Col1, Col2;
+register  ElementPtr  pPivot;
 {
-register  ElementPtr  Col1Ptr, Col2Ptr;
-int  Row;
-ElementPtr  Element1, Element2;
+#if REAL
+register  ElementPtr  pSub;
+register  int  Row;
+register  ElementPtr  pLower, pUpper;
+extern ElementPtr  CreateFillin();
 
-/* Begin `spcColExchange'. */
-    if (Col1 > Col2)  SWAP(int, Col1, Col2);
+/* Begin `RealRowColElimination'. */
 
-    Col1Ptr = Matrix->FirstInCol[Col1];
-    Col2Ptr = Matrix->FirstInCol[Col2];
-    while (Col1Ptr != NULL OR Col2Ptr != NULL)
+/* Test for zero pivot. */
+    if (ABS(pPivot->Real) == 0.0)
+    {   (void)MatrixIsSingular( Matrix, pPivot->Row );
+        return;
+    }
+    pPivot->Real = 1.0 / pPivot->Real;
+
+    pUpper = pPivot->NextInRow;
+    while (pUpper != NULL)
     {
-/* Exchange elements in rows while traveling from top to bottom. */
-        if (Col1Ptr == NULL)
-        {   Row = Col2Ptr->Row;
-            Element1 = NULL;
-            Element2 = Col2Ptr;
-            Col2Ptr = Col2Ptr->NextInCol;
-        }
-        else if (Col2Ptr == NULL)
-        {   Row = Col1Ptr->Row;
-            Element1 = Col1Ptr;
-            Element2 = NULL;
-            Col1Ptr = Col1Ptr->NextInCol;
-        }
-        else if (Col1Ptr->Row < Col2Ptr->Row)
-        {   Row = Col1Ptr->Row;
-            Element1 = Col1Ptr;
-            Element2 = NULL;
-            Col1Ptr = Col1Ptr->NextInCol;
-        }
-        else if (Col1Ptr->Row > Col2Ptr->Row)
-        {   Row = Col2Ptr->Row;
-            Element1 = NULL;
-            Element2 = Col2Ptr;
-            Col2Ptr = Col2Ptr->NextInCol;
-        }
-        else   /* Col1Ptr->Row == Col2Ptr->Row */
-        {   Row = Col1Ptr->Row;
-            Element1 = Col1Ptr;
-            Element2 = Col2Ptr;
-            Col1Ptr = Col1Ptr->NextInCol;
-            Col2Ptr = Col2Ptr->NextInCol;
-        }
+/* Calculate upper triangular element. */
+        pUpper->Real *= pPivot->Real;
 
-        ExchangeRowElements( Matrix, Col1, Element1, Col2, Element2, Row);
-    }  /* end of while(Col1Ptr != NULL OR Col2Ptr != NULL) */
+        pSub = pUpper->NextInCol;
+        pLower = pPivot->NextInCol;
+        while (pLower != NULL)
+        {   Row = pLower->Row;
 
-    if (Matrix->InternalVectorsAllocated)
-        SWAP( int, Matrix->MarkowitzCol[Col1], Matrix->MarkowitzCol[Col2]);
-    SWAP( ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]);
-    SWAP( int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]);
-#if TRANSLATE
-    Matrix->ExtToIntColMap[ Matrix->IntToExtColMap[Col1] ] = Col1;
-    Matrix->ExtToIntColMap[ Matrix->IntToExtColMap[Col2] ] = Col2;
-#endif
+/* Find element in row that lines up with current lower triangular element. */
+            while (pSub != NULL AND pSub->Row < Row)
+                pSub = pSub->NextInCol;
 
+/* Test to see if desired element was not found, if not, create fill-in. */
+            if (pSub == NULL OR pSub->Row > Row)
+            {   pSub = CreateFillin( Matrix, Row, pUpper->Col );
+                if (pSub == NULL)
+                {   Matrix->Error = spNO_MEMORY;
+                    return;
+                }
+            }
+            pSub->Real -= pUpper->Real * pLower->Real;
+            pSub = pSub->NextInCol;
+            pLower = pLower->NextInCol;
+        }
+        pUpper = pUpper->NextInRow;
+    }
     return;
+#endif /* REAL */
 }
 
 
@@ -2409,279 +2055,168 @@
 
 
 
+
+
+
 
 /*
- *  EXCHANGE TWO ELEMENTS IN A COLUMN
+ *  PERFORM ROW AND COLUMN ELIMINATION ON COMPLEX MATRIX
  *
- *  Performs all required operations to exchange two elements in a column.
- *  Those operations are: restring NextInCol pointers and swapping row indexes
- *  in the MatrixElements.
+ *  Eliminates a single row and column of the matrix and leaves single row of
+ *  the upper triangular matrix and a single column of the lower triangular
+ *  matrix in its wake.  Uses Gauss's method.
  *
- *  >>> Arguments:
+ *  >>> Argument:
  *  Matrix  <input>  (MatrixPtr)
  *      Pointer to the matrix.
- *  Row1  <input>  (int)
- *      Row of top element to be exchanged.
- *  Element1  <input>  (ElementPtr)
- *      Pointer to top element to be exchanged.
- *  Row2  <input>  (int)
- *      Row of bottom element to be exchanged.
- *  Element2  <input>  (ElementPtr)
- *      Pointer to bottom element to be exchanged.
- *  Column <input>  (int)
- *      Column that exchange is to take place in.
+ *  pPivot  <input>  (ElementPtr)
+ *      Pointer to the current pivot.
  *
  *  >>> Local variables:
- *  ElementAboveRow1  (ElementPtr *)
- *      Location of pointer which points to the element above Element1. This
- *      pointer is modified so that it points to correct element on exit.
- *  ElementAboveRow2  (ElementPtr *)
- *      Location of pointer which points to the element above Element2. This
- *      pointer is modified so that it points to correct element on exit.
- *  ElementBelowRow1  (ElementPtr)
- *      Pointer to element below Element1.
- *  ElementBelowRow2  (ElementPtr)
- *      Pointer to element below Element2.
- *  pElement  (ElementPtr)
- *      Pointer used to traverse the column.
+ *  pLower  (ElementPtr)
+ *      Points to matrix element in lower triangular column.
+ *  pSub (ElementPtr)
+ *      Points to elements in the reduced submatrix.
+ *  Row  (int)
+ *      Row index.
+ *  pUpper  (ElementPtr)
+ *      Points to matrix element in upper triangular row.
+ *
+ *  Possible errors:
+ *  spNO_MEMORY
  */
 
 static
-ExchangeColElements( Matrix, Row1, Element1, Row2, Element2, Column )
+ComplexRowColElimination( Matrix, pPivot )
 
 MatrixPtr Matrix;
-register  ElementPtr  Element1, Element2;
-int  Row1, Row2, Column;
+register  ElementPtr  pPivot;
 {
-ElementPtr  *ElementAboveRow1, *ElementAboveRow2;
-ElementPtr  ElementBelowRow1, ElementBelowRow2;
-register  ElementPtr  pElement;
+#if spCOMPLEX
+register  ElementPtr  pSub;
+register  int  Row;
+register  ElementPtr  pLower, pUpper;
+ElementPtr  CreateFillin();
 
-/* Begin `ExchangeColElements'. */
-/* Search to find the ElementAboveRow1. */
-    ElementAboveRow1 = &(Matrix->FirstInCol[Column]);
-    pElement = *ElementAboveRow1;
-    while (pElement->Row < Row1)
-    {   ElementAboveRow1 = &(pElement->NextInCol);
-        pElement = *ElementAboveRow1;
+/* Begin `ComplexRowColElimination'. */
+
+/* Test for zero pivot. */
+    if (ELEMENT_MAG(pPivot) == 0.0)
+    {   (void)MatrixIsSingular( Matrix, pPivot->Row );
+        return;
     }
-    if (Element1 != NULL)
-    {   ElementBelowRow1 = Element1->NextInCol;
-        if (Element2 == NULL)
-        {
-/* Element2 does not exist, move Element1 down to Row2. */
-            if ( ElementBelowRow1 != NULL AND ElementBelowRow1->Row < Row2 )
-            {
-/* Element1 must be removed from linked list and moved. */
-                *ElementAboveRow1 = ElementBelowRow1;
+    CMPLX_RECIPROCAL(*pPivot, *pPivot);
 
-/* Search column for Row2. */
-                pElement = ElementBelowRow1;
-                do
-                {   ElementAboveRow2 = &(pElement->NextInCol);
-                    pElement = *ElementAboveRow2;
-                }   while (pElement != NULL AND pElement->Row < Row2);
+    pUpper = pPivot->NextInRow;
+    while (pUpper != NULL)
+    {
+/* Calculate upper triangular element. */
+/* Cmplx expr: *pUpper = *pUpper * (1.0 / *pPivot). */
+        CMPLX_MULT_ASSIGN(*pUpper, *pPivot);
 
-/* Place Element1 in Row2. */
-                *ElementAboveRow2 = Element1;
-                Element1->NextInCol = pElement;
-                *ElementAboveRow1 =ElementBelowRow1;
-            }
-            Element1->Row = Row2;
-        }
-        else
-        {
-/* Element2 does exist, and the two elements must be exchanged. */
-            if ( ElementBelowRow1->Row == Row2)
-            {
-/* Element2 is just below Element1, exchange them. */
-                Element1->NextInCol = Element2->NextInCol;
-                Element2->NextInCol = Element1;
-                *ElementAboveRow1 = Element2;
-            }
-            else
-            {
-/* Element2 is not just below Element1 and must be searched for. */
-                pElement = ElementBelowRow1;
-                do
-                {   ElementAboveRow2 = &(pElement->NextInCol);
-                    pElement = *ElementAboveRow2;
-                }   while (pElement->Row < Row2);
+        pSub = pUpper->NextInCol;
+        pLower = pPivot->NextInCol;
+        while (pLower != NULL)
+        {   Row = pLower->Row;
 
-                ElementBelowRow2 = Element2->NextInCol;
+/* Find element in row that lines up with current lower triangular element. */
+            while (pSub != NULL AND pSub->Row < Row)
+                pSub = pSub->NextInCol;
 
-/* Switch Element1 and Element2. */
-                *ElementAboveRow1 = Element2;
-                Element2->NextInCol = ElementBelowRow1;
-                *ElementAboveRow2 = Element1;
-                Element1->NextInCol = ElementBelowRow2;
+/* Test to see if desired element was not found, if not, create fill-in. */
+            if (pSub == NULL OR pSub->Row > Row)
+            {   pSub = CreateFillin( Matrix, Row, pUpper->Col );
+                if (pSub == NULL)
+                {   Matrix->Error = spNO_MEMORY;
+                    return;
+                }
             }
-            Element1->Row = Row2;
-            Element2->Row = Row1;
-        }
-    }
-    else
-    {
-/* Element1 does not exist. */
-        ElementBelowRow1 = pElement;
-
-/* Find Element2. */
-        if (ElementBelowRow1->Row != Row2)
-        {   do
-            {   ElementAboveRow2 = &(pElement->NextInCol);
-                pElement = *ElementAboveRow2;
-            }   while (pElement->Row < Row2);
 
-        ElementBelowRow2 = Element2->NextInCol;
-
-/* Move Element2 to Row1. */
-            *ElementAboveRow2 = Element2->NextInCol;
-            *ElementAboveRow1 = Element2;
-            Element2->NextInCol = ElementBelowRow1;
+/* Cmplx expr: pElement -= *pUpper * pLower. */
+            CMPLX_MULT_SUBT_ASSIGN(*pSub, *pUpper, *pLower);
+            pSub = pSub->NextInCol;
+            pLower = pLower->NextInCol;
         }
-        Element2->Row = Row1;
+        pUpper = pUpper->NextInRow;
     }
     return;
+#endif /* spCOMPLEX */
 }
 
 
 
-
-
-
 
 /*
- *  EXCHANGE TWO ELEMENTS IN A ROW
+ *  UPDATE MARKOWITZ NUMBERS
  *
- *  Performs all required operations to exchange two elements in a row.
- *  Those operations are: restring NextInRow pointers and swapping column
- *  indexes in the MatrixElements.
+ *  Updates the Markowitz numbers after a row and column have been eliminated.
+ *  Also updates singleton count.
  *
- *  >>> Arguments:
+ *  >>> Argument:
  *  Matrix  <input>  (MatrixPtr)
  *      Pointer to the matrix.
- *  Col1  <input>  (int)
- *      Col of left-most element to be exchanged.
- *  Element1  <input>  (ElementPtr)
- *      Pointer to left-most element to be exchanged.
- *  Col2  <input>  (int)
- *      Col of right-most element to be exchanged.
- *  Element2  <input>  (ElementPtr)
- *      Pointer to right-most element to be exchanged.
- *  Row <input>  (int)
- *      Row that exchange is to take place in.
+ *  pPivot  <input>  (ElementPtr)
+ *      Pointer to the current pivot.
  *
  *  >>> Local variables:
- *  ElementLeftOfCol1  (ElementPtr *)
- *      Location of pointer which points to the element to the left of
- *      Element1. This pointer is modified so that it points to correct
- *      element on exit.
- *  ElementLeftOfCol2  (ElementPtr *)
- *      Location of pointer which points to the element to the left of
- *      Element2. This pointer is modified so that it points to correct
- *      element on exit.
- *  ElementRightOfCol1  (ElementPtr)
- *      Pointer to element right of Element1.
- *  ElementRightOfCol2  (ElementPtr)
- *      Pointer to element right of Element2.
- *  pElement  (ElementPtr)
- *      Pointer used to traverse the row.
+ *  Row  (int)
+ *      Row index.
+ *  Col  (int)
+ *      Column index.
+ *  ColPtr  (ElementPtr)
+ *      Points to matrix element in upper triangular column.
+ *  RowPtr  (ElementPtr)
+ *      Points to matrix element in lower triangular row.
  */
 
 static
-ExchangeRowElements( Matrix, Col1, Element1, Col2, Element2, Row )
+UpdateMarkowitzNumbers( Matrix, pPivot )
 
 MatrixPtr Matrix;
-int  Col1, Col2, Row;
-register ElementPtr  Element1, Element2;
+ElementPtr  pPivot;
 {
-ElementPtr  *ElementLeftOfCol1, *ElementLeftOfCol2;
-ElementPtr  ElementRightOfCol1, ElementRightOfCol2;
-register   ElementPtr  pElement;
+register  int  Row, Col;
+register  ElementPtr  ColPtr, RowPtr;
+register  int *MarkoRow = Matrix->MarkowitzRow, *MarkoCol = Matrix->MarkowitzCol;
+double Product;
 
-/* Begin `ExchangeRowElements'. */
-/* Search to find the ElementLeftOfCol1. */
-    ElementLeftOfCol1 = &(Matrix->FirstInRow[Row]);
-    pElement = *ElementLeftOfCol1;
-    while (pElement->Col < Col1)
-    {   ElementLeftOfCol1 = &(pElement->NextInRow);
-        pElement = *ElementLeftOfCol1;
-    }
-    if (Element1 != NULL)
-    {   ElementRightOfCol1 = Element1->NextInRow;
-        if (Element2 == NULL)
-        {
-/* Element2 does not exist, move Element1 to right to Col2. */
-            if ( ElementRightOfCol1 != NULL AND ElementRightOfCol1->Col < Col2 )
-            {
-/* Element1 must be removed from linked list and moved. */
-                *ElementLeftOfCol1 = ElementRightOfCol1;
+/* Begin `UpdateMarkowitzNumbers'. */
 
-/* Search Row for Col2. */
-                pElement = ElementRightOfCol1;
-                do
-                {   ElementLeftOfCol2 = &(pElement->NextInRow);
-                    pElement = *ElementLeftOfCol2;
-                }   while (pElement != NULL AND pElement->Col < Col2);
+/* Update Markowitz numbers. */
+    for (ColPtr = pPivot->NextInCol; ColPtr != NULL; ColPtr = ColPtr->NextInCol)
+    {   Row = ColPtr->Row;
+        --MarkoRow[Row];
 
-/* Place Element1 in Col2. */
-                *ElementLeftOfCol2 = Element1;
-                Element1->NextInRow = pElement;
-                *ElementLeftOfCol1 =ElementRightOfCol1;
-            }
-            Element1->Col = Col2;
-        }
-        else
-        {
-/* Element2 does exist, and the two elements must be exchanged. */
-            if ( ElementRightOfCol1->Col == Col2)
-            {
-/* Element2 is just right of Element1, exchange them. */
-                Element1->NextInRow = Element2->NextInRow;
-                Element2->NextInRow = Element1;
-                *ElementLeftOfCol1 = Element2;
-            }
+/* Form Markowitz product while being cautious of overflows. */
+        if ((MarkoRow[Row] > LARGEST_SHORT_INTEGER AND MarkoCol[Row] != 0) OR
+            (MarkoCol[Row] > LARGEST_SHORT_INTEGER AND MarkoRow[Row] != 0))
+        {   Product = MarkoCol[Row] * MarkoRow[Row];
+            if (Product >= LARGEST_LONG_INTEGER)
+                Matrix->MarkowitzProd[Row] = LARGEST_LONG_INTEGER;
             else
-            {
-/* Element2 is not just right of Element1 and must be searched for. */
-                pElement = ElementRightOfCol1;
-                do
-                {   ElementLeftOfCol2 = &(pElement->NextInRow);
-                    pElement = *ElementLeftOfCol2;
-                }   while (pElement->Col < Col2);
-
-                ElementRightOfCol2 = Element2->NextInRow;
-
-/* Switch Element1 and Element2. */
-                *ElementLeftOfCol1 = Element2;
-                Element2->NextInRow = ElementRightOfCol1;
-                *ElementLeftOfCol2 = Element1;
-                Element1->NextInRow = ElementRightOfCol2;
-            }
-            Element1->Col = Col2;
-            Element2->Col = Col1;
+                Matrix->MarkowitzProd[Row] = Product;
         }
+        else Matrix->MarkowitzProd[Row] = MarkoRow[Row] * MarkoCol[Row];
+        if (MarkoRow[Row] == 0)
+            Matrix->Singletons++;
     }
-    else
-    {
-/* Element1 does not exist. */
-        ElementRightOfCol1 = pElement;
 
-/* Find Element2. */
-        if (ElementRightOfCol1->Col != Col2)
-        {   do
-            {   ElementLeftOfCol2 = &(pElement->NextInRow);
-                pElement = *ElementLeftOfCol2;
-            }   while (pElement->Col < Col2);
-
-            ElementRightOfCol2 = Element2->NextInRow;
+    for (RowPtr = pPivot->NextInRow; RowPtr != NULL; RowPtr = RowPtr->NextInRow)
+    {   Col = RowPtr->Col;
+        --MarkoCol[Col];
 
-/* Move Element2 to Col1. */
-            *ElementLeftOfCol2 = Element2->NextInRow;
-            *ElementLeftOfCol1 = Element2;
-            Element2->NextInRow = ElementRightOfCol1;
+/* Form Markowitz product while being cautious of overflows. */
+        if ((MarkoRow[Col] > LARGEST_SHORT_INTEGER AND MarkoCol[Col] != 0) OR
+            (MarkoCol[Col] > LARGEST_SHORT_INTEGER AND MarkoRow[Col] != 0))
+        {   Product = MarkoCol[Col] * MarkoRow[Col];
+            if (Product >= LARGEST_LONG_INTEGER)
+                Matrix->MarkowitzProd[Col] = LARGEST_LONG_INTEGER;
+            else
+                Matrix->MarkowitzProd[Col] = Product;
         }
-        Element2->Col = Col1;
+        else Matrix->MarkowitzProd[Col] = MarkoRow[Col] * MarkoCol[Col];
+        if ((MarkoCol[Col] == 0) AND (MarkoRow[Col] != 0))
+            Matrix->Singletons++;
     }
     return;
 }
@@ -2697,85 +2232,226 @@
 
 
 /*
- *  PERFORM ROW AND COLUMN ELIMINATION ON REAL MATRIX
+ *  ORDER AND FACTOR MATRIX
  *
- *  Eliminates a single row and column of the matrix and leaves single row of
- *  the upper triangular matrix and a single column of the lower triangular
- *  matrix in its wake.  Uses Gauss's method.
+ *  This routine chooses a pivot order for the matrix and factors it
+ *  into LU form.  It handles both the initial factorization and subsequent
+ *  factorizations when a reordering is desired.  This is handled in a manner
+ *  that is transparent to the user.  The routine uses a variation of
+ *  Gauss's method where the pivots are associated with L and the
+ *  diagonal terms of U are one.
  *
- *  >>> Argument:
- *  Matrix  <input>  (MatrixPtr)
- *      Pointer to the matrix.
- *  pPivot  <input>  (ElementPtr)
- *      Pointer to the current pivot.
+ *  >>> Returned:
+ *  The error code is returned.  Possible errors are listed below.
+ *
+ *  >>> Arguments:
+ *  Matrix  <input>  (char *)
+ *      Pointer to matrix.
+ *  RHS  <input>  (RealVector)
+ *      Representative right-hand side vector that is used to determine
+ *      pivoting order when the right hand side vector is sparse.  If
+ *      RHS is a NULL pointer then the RHS vector is assumed to
+ *      be full and it is not used when determining the pivoting
+ *      order.
+ *  RelThreshold  <input>  (RealNumber)
+ *      This number determines what the pivot relative threshold will
+ *      be.  It should be between zero and one.  If it is one then the
+ *      pivoting method becomes complete pivoting, which is very slow
+ *      and tends to fill up the matrix.  If it is set close to zero
+ *      the pivoting method becomes strict Markowitz with no
+ *      threshold.  The pivot threshold is used to eliminate pivot
+ *      candidates that would cause excessive element growth if they
+ *      were used.  Element growth is the cause of roundoff error.
+ *      Element growth occurs even in well-conditioned matrices.
+ *      Setting the RelThreshold large will reduce element growth and
+ *      roundoff error, but setting it too large will cause execution
+ *      time to be excessive and will result in a large number of
+ *      fill-ins.  If this occurs, accuracy can actually be degraded
+ *      because of the large number of operations required on the
+ *      matrix due to the large number of fill-ins.  A good value seems
+ *      to be 0.001.  The default is chosen by giving a value larger
+ *      than one or less than or equal to zero.  This value should be
+ *      increased and the matrix resolved if growth is found to be
+ *      excessive.  Changing the pivot threshold does not improve
+ *      performance on matrices where growth is low, as is often the
+ *      case with ill-conditioned matrices.  Once a valid threshold is
+ *      given, it becomes the new default.  The default value of
+ *      RelThreshold was choosen for use with nearly diagonally
+ *      dominant matrices such as node- and modified-node admittance
+ *      matrices.  For these matrices it is usually best to use
+ *      diagonal pivoting.  For matrices without a strong diagonal, it
+ *      is usually best to use a larger threshold, such as 0.01 or
+ *      0.1.
+ *  AbsThreshold  <input>  (RealNumber)
+ *      The absolute magnitude an element must have to be considered
+ *      as a pivot candidate, except as a last resort.  This number
+ *      should be set significantly smaller than the smallest diagonal
+ *      element that is is expected to be placed in the matrix.  If
+ *      there is no reasonable prediction for the lower bound on these
+ *      elements, then AbsThreshold should be set to zero.
+ *      AbsThreshold is used to reduce the possibility of choosing as a
+ *      pivot an element that has suffered heavy cancellation and as a
+ *      result mainly consists of roundoff error.  Once a valid
+ *      threshold is given, it becomes the new default.
+ *  DiagPivoting  <input>  (BOOLEAN)
+ *      A flag indicating that pivot selection should be confined to the
+ *      diagonal if possible.  If DiagPivoting is nonzero and if
+ *      DIAGONAL_PIVOTING is enabled pivots will be chosen only from
+ *      the diagonal unless there are no diagonal elements that satisfy
+ *      the threshold criteria.  Otherwise, the entire reduced
+ *      submatrix is searched when looking for a pivot.  The diagonal
+ *      pivoting in Sparse is efficient and well refined, while the
+ *      off-diagonal pivoting is not.  For symmetric and near symmetric
+ *      matrices, it is best to use diagonal pivoting because it
+ *      results in the best performance when reordering the matrix and
+ *      when factoring the matrix without ordering.  If there is a
+ *      considerable amount of nonsymmetry in the matrix, then
+ *      off-diagonal pivoting may result in a better equation ordering
+ *      simply because there are more pivot candidates to choose from.
+ *      A better ordering results in faster subsequent factorizations.
+ *      However, the initial pivot selection process takes considerably
+ *      longer for off-diagonal pivoting.
  *
  *  >>> Local variables:
- *  pLower  (ElementPtr)
- *      Points to matrix element in lower triangular column.
- *  pSub (ElementPtr)
- *      Points to elements in the reduced submatrix.
- *  Row  (int)
- *      Row index.
- *  pUpper  (ElementPtr)
- *      Points to matrix element in upper triangular row.
+ *  pPivot  (ElementPtr)
+ *      Pointer to the element being used as a pivot.
+ *  ReorderingRequired  (BOOLEAN)
+ *      Flag that indicates whether reordering is required.
  *
  *  >>> Possible errors:
  *  spNO_MEMORY
+ *  spSINGULAR
+ *  spSMALL_PIVOT
+ *  Error is cleared in this function.
+ */
+
+int
+spOrderAndFactor( eMatrix, RHS, RelThreshold, AbsThreshold, DiagPivoting )
+
+char *eMatrix;
+RealNumber  RHS[], RelThreshold, AbsThreshold;
+BOOLEAN DiagPivoting;
+{
+MatrixPtr  Matrix = (MatrixPtr)eMatrix;
+ElementPtr  pPivot;
+int  Step, Size, ReorderingRequired;
+ElementPtr SearchForPivot();
+RealNumber LargestInCol, FindLargestInCol();
+
+/* Begin `spOrderAndFactor'. */
+    ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored);
+
+    Matrix->Error = spOKAY;
+    Size = Matrix->Size;
+    if (RelThreshold <= 0.0) RelThreshold = Matrix->RelThreshold;
+    if (RelThreshold > 1.0) RelThreshold = Matrix->RelThreshold;
+    Matrix->RelThreshold = RelThreshold;
+    if (AbsThreshold < 0.0) AbsThreshold = Matrix->AbsThreshold;
+    Matrix->AbsThreshold = AbsThreshold;
+    ReorderingRequired = NO;
+
+    if (NOT Matrix->NeedsOrdering)
+    {
+/* Matrix has been factored before and reordering is not required. */
+        for (Step = 1; Step <= Size; Step++)
+        {   pPivot = Matrix->Diag[Step];
+            LargestInCol = FindLargestInCol(pPivot->NextInCol);
+            if ((LargestInCol * RelThreshold < ELEMENT_MAG(pPivot)))
+            {   if (Matrix->Complex)
+                    ComplexRowColElimination( Matrix, pPivot );
+                else
+                    RealRowColElimination( Matrix, pPivot );
+            }
+            else
+            {   ReorderingRequired = YES;
+                break; /* for loop */
+            }
+        }
+        if (NOT ReorderingRequired)
+            goto Done;
+        else
+        {
+/*
+ * A pivot was not large enough to maintain accuracy,
+ * so a partial reordering is required.
  */
 
-static
-RealRowColElimination( Matrix, pPivot )
+#if ANNOTATE >= ON_STRANGE_BEHAVIOR
+            printf("Reordering,  Step = %1d\n", Step);
+#endif
+        }
+    } /* End of if(NOT Matrix->NeedsOrdering) */
+    else
+    {
+/*
+ * This is the first time the matrix has been factored.  These few statements
+ * indicate to the rest of the code that a full reodering is required rather
+ * than a partial reordering, which occurs during a failure of a fast
+ * factorization.
+ */
+        Step = 1;
+        if (NOT Matrix->RowsLinked)
+            spcLinkRows( Matrix );
+        if (NOT Matrix->InternalVectorsAllocated)
+            CreateInternalVectors( Matrix );
+        if (Matrix->Error >= spFATAL)
+            return Matrix->Error;
+    }
+
+/* Form initial Markowitz products. */
+    CountMarkowitz( Matrix, RHS, Step );
+    MarkowitzProducts( Matrix, Step );
+    Matrix->MaxRowCountInLowerTri = -1;
+
+/* Perform reordering and factorization. */
+    for (; Step <= Size; Step++)
+    {   pPivot = SearchForPivot( Matrix, Step, DiagPivoting );
+        if (pPivot == NULL) return MatrixIsSingular( Matrix, Step );
+        ExchangeRowsAndCols( Matrix, pPivot, Step );
+
+        if (Matrix->Complex)
+            ComplexRowColElimination( Matrix, pPivot );
+        else
+            RealRowColElimination( Matrix, pPivot );
+
+        if (Matrix->Error >= spFATAL) return Matrix->Error;
+        UpdateMarkowitzNumbers( Matrix, pPivot );
+
+#if ANNOTATE == FULL
+        WriteStatus( Matrix, Step );
+#endif
+    }
+
+Done:
+    Matrix->NeedsOrdering = NO;
+    Matrix->Reordered = YES;
+    Matrix->Factored = YES;
+
+    return Matrix->Error;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
-MatrixPtr Matrix;
-register  ElementPtr  pPivot;
-{
-#if REAL
-register  ElementPtr  pSub;
-register  int  Row;
-register  ElementPtr  pLower, pUpper;
-extern ElementPtr  CreateFillin();
 
-/* Begin `RealRowColElimination'. */
 
-/* Test for zero pivot. */
-    if (ABS(pPivot->Real) == 0.0)
-    {   (void)MatrixIsSingular( Matrix, pPivot->Row );
-        return;
-    }
-    pPivot->Real = 1.0 / pPivot->Real;
 
-    pUpper = pPivot->NextInRow;
-    while (pUpper != NULL)
-    {
-/* Calculate upper triangular element. */
-        pUpper->Real *= pPivot->Real;
 
-        pSub = pUpper->NextInCol;
-        pLower = pPivot->NextInCol;
-        while (pLower != NULL)
-        {   Row = pLower->Row;
 
-/* Find element in row that lines up with current lower triangular element. */
-            while (pSub != NULL AND pSub->Row < Row)
-                pSub = pSub->NextInCol;
 
-/* Test to see if desired element was not found, if not, create fill-in. */
-            if (pSub == NULL OR pSub->Row > Row)
-            {   pSub = CreateFillin( Matrix, Row, pUpper->Col );
-                if (pSub == NULL)
-                {   Matrix->Error = spNO_MEMORY;
-                    return;
-                }
-            }
-            pSub->Real -= pUpper->Real * pLower->Real;
-            pSub = pSub->NextInCol;
-            pLower = pLower->NextInCol;
-        }
-        pUpper = pUpper->NextInRow;
-    }
-    return;
-#endif /* REAL */
-}
 
 
 
@@ -2786,165 +2462,240 @@
 
 
 /*
- *  PERFORM ROW AND COLUMN ELIMINATION ON COMPLEX MATRIX
+ *  EXCHANGE TWO ELEMENTS IN A COLUMN
  *
- *  Eliminates a single row and column of the matrix and leaves single row of
- *  the upper triangular matrix and a single column of the lower triangular
- *  matrix in its wake.  Uses Gauss's method.
+ *  Performs all required operations to exchange two elements in a column.
+ *  Those operations are: restring NextInCol pointers and swapping row indexes
+ *  in the MatrixElements.
  *
- *  >>> Argument:
+ *  >>> Arguments:
  *  Matrix  <input>  (MatrixPtr)
  *      Pointer to the matrix.
- *  pPivot  <input>  (ElementPtr)
- *      Pointer to the current pivot.
+ *  Row1  <input>  (int)
+ *      Row of top element to be exchanged.
+ *  Element1  <input>  (ElementPtr)
+ *      Pointer to top element to be exchanged.
+ *  Row2  <input>  (int)
+ *      Row of bottom element to be exchanged.
+ *  Element2  <input>  (ElementPtr)
+ *      Pointer to bottom element to be exchanged.
+ *  Column <input>  (int)
+ *      Column that exchange is to take place in.
  *
  *  >>> Local variables:
- *  pLower  (ElementPtr)
- *      Points to matrix element in lower triangular column.
- *  pSub (ElementPtr)
- *      Points to elements in the reduced submatrix.
- *  Row  (int)
- *      Row index.
- *  pUpper  (ElementPtr)
- *      Points to matrix element in upper triangular row.
- *
- *  Possible errors:
- *  spNO_MEMORY
+ *  ElementAboveRow1  (ElementPtr *)
+ *      Location of pointer which points to the element above Element1. This
+ *      pointer is modified so that it points to correct element on exit.
+ *  ElementAboveRow2  (ElementPtr *)
+ *      Location of pointer which points to the element above Element2. This
+ *      pointer is modified so that it points to correct element on exit.
+ *  ElementBelowRow1  (ElementPtr)
+ *      Pointer to element below Element1.
+ *  ElementBelowRow2  (ElementPtr)
+ *      Pointer to element below Element2.
+ *  pElement  (ElementPtr)
+ *      Pointer used to traverse the column.
  */
 
 static
-ComplexRowColElimination( Matrix, pPivot )
+ExchangeColElements( Matrix, Row1, Element1, Row2, Element2, Column )
 
 MatrixPtr Matrix;
-register  ElementPtr  pPivot;
+register  ElementPtr  Element1, Element2;
+int  Row1, Row2, Column;
 {
-#if spCOMPLEX
-register  ElementPtr  pSub;
-register  int  Row;
-register  ElementPtr  pLower, pUpper;
-ElementPtr  CreateFillin();
-
-/* Begin `ComplexRowColElimination'. */
+ElementPtr  *ElementAboveRow1, *ElementAboveRow2;
+ElementPtr  ElementBelowRow1, ElementBelowRow2;
+register  ElementPtr  pElement;
 
-/* Test for zero pivot. */
-    if (ELEMENT_MAG(pPivot) == 0.0)
-    {   (void)MatrixIsSingular( Matrix, pPivot->Row );
-        return;
+/* Begin `ExchangeColElements'. */
+/* Search to find the ElementAboveRow1. */
+    ElementAboveRow1 = &(Matrix->FirstInCol[Column]);
+    pElement = *ElementAboveRow1;
+    while (pElement->Row < Row1)
+    {   ElementAboveRow1 = &(pElement->NextInCol);
+        pElement = *ElementAboveRow1;
     }
-    CMPLX_RECIPROCAL(*pPivot, *pPivot);
+    if (Element1 != NULL)
+    {   ElementBelowRow1 = Element1->NextInCol;
+        if (Element2 == NULL)
+        {
+/* Element2 does not exist, move Element1 down to Row2. */
+            if ( ElementBelowRow1 != NULL AND ElementBelowRow1->Row < Row2 )
+            {
+/* Element1 must be removed from linked list and moved. */
+                *ElementAboveRow1 = ElementBelowRow1;
 
-    pUpper = pPivot->NextInRow;
-    while (pUpper != NULL)
-    {
-/* Calculate upper triangular element. */
-/* Cmplx expr: *pUpper = *pUpper * (1.0 / *pPivot). */
-        CMPLX_MULT_ASSIGN(*pUpper, *pPivot);
+/* Search column for Row2. */
+                pElement = ElementBelowRow1;
+                do
+                {   ElementAboveRow2 = &(pElement->NextInCol);
+                    pElement = *ElementAboveRow2;
+                }   while (pElement != NULL AND pElement->Row < Row2);
 
-        pSub = pUpper->NextInCol;
-        pLower = pPivot->NextInCol;
-        while (pLower != NULL)
-        {   Row = pLower->Row;
+/* Place Element1 in Row2. */
+                *ElementAboveRow2 = Element1;
+                Element1->NextInCol = pElement;
+                *ElementAboveRow1 =ElementBelowRow1;
+            }
+            Element1->Row = Row2;
+        }
+        else
+        {
+/* Element2 does exist, and the two elements must be exchanged. */
+            if ( ElementBelowRow1->Row == Row2)
+            {
+/* Element2 is just below Element1, exchange them. */
+                Element1->NextInCol = Element2->NextInCol;
+                Element2->NextInCol = Element1;
+                *ElementAboveRow1 = Element2;
+            }
+            else
+            {
+/* Element2 is not just below Element1 and must be searched for. */
+                pElement = ElementBelowRow1;
+                do
+                {   ElementAboveRow2 = &(pElement->NextInCol);
+                    pElement = *ElementAboveRow2;
+                }   while (pElement->Row < Row2);
 
-/* Find element in row that lines up with current lower triangular element. */
-            while (pSub != NULL AND pSub->Row < Row)
-                pSub = pSub->NextInCol;
+                ElementBelowRow2 = Element2->NextInCol;
 
-/* Test to see if desired element was not found, if not, create fill-in. */
-            if (pSub == NULL OR pSub->Row > Row)
-            {   pSub = CreateFillin( Matrix, Row, pUpper->Col );
-                if (pSub == NULL)
-                {   Matrix->Error = spNO_MEMORY;
-                    return;
-                }
+/* Switch Element1 and Element2. */
+                *ElementAboveRow1 = Element2;
+                Element2->NextInCol = ElementBelowRow1;
+                *ElementAboveRow2 = Element1;
+                Element1->NextInCol = ElementBelowRow2;
             }
+            Element1->Row = Row2;
+            Element2->Row = Row1;
+        }
+    }
+    else
+    {
+/* Element1 does not exist. */
+        ElementBelowRow1 = pElement;
 
-/* Cmplx expr: pElement -= *pUpper * pLower. */
-            CMPLX_MULT_SUBT_ASSIGN(*pSub, *pUpper, *pLower);
-            pSub = pSub->NextInCol;
-            pLower = pLower->NextInCol;
+/* Find Element2. */
+        if (ElementBelowRow1->Row != Row2)
+        {   do
+            {   ElementAboveRow2 = &(pElement->NextInCol);
+                pElement = *ElementAboveRow2;
+            }   while (pElement->Row < Row2);
+
+        ElementBelowRow2 = Element2->NextInCol;
+
+/* Move Element2 to Row1. */
+            *ElementAboveRow2 = Element2->NextInCol;
+            *ElementAboveRow1 = Element2;
+            Element2->NextInCol = ElementBelowRow1;
         }
-        pUpper = pUpper->NextInRow;
+        Element2->Row = Row1;
     }
     return;
-#endif /* spCOMPLEX */
 }
 
 
 
 
+
+
+
+
+
+
+
 
 /*
- *  UPDATE MARKOWITZ NUMBERS
+ *  EXCHANGE ROWS
  *
- *  Updates the Markowitz numbers after a row and column have been eliminated.
- *  Also updates singleton count.
+ *  Performs all required operations to exchange two rows. Those operations
+ *  include: swap FirstInRow pointers, fixing up the NextInCol pointers,
+ *  swapping row indexes in MatrixElements, and swapping Markowitz row
+ *  counts.
  *
- *  >>> Argument:
+ *  >>> Arguments:
  *  Matrix  <input>  (MatrixPtr)
  *      Pointer to the matrix.
- *  pPivot  <input>  (ElementPtr)
- *      Pointer to the current pivot.
- *
- *  >>> Local variables:
- *  Row  (int)
- *      Row index.
- *  Col  (int)
- *      Column index.
- *  ColPtr  (ElementPtr)
- *      Points to matrix element in upper triangular column.
- *  RowPtr  (ElementPtr)
- *      Points to matrix element in lower triangular row.
+ *  Row1  <input>  (int)
+ *      Row index of one of the rows, becomes the smallest index.
+ *  Row2  <input>  (int)
+ *      Row index of the other row, becomes the largest index.
+ *
+ *  Local variables:
+ *  Column  (int)
+ *      Column in which row elements are currently being exchanged.
+ *  Row1Ptr  (ElementPtr)
+ *      Pointer to an element in Row1.
+ *  Row2Ptr  (ElementPtr)
+ *      Pointer to an element in Row2.
+ *  Element1  (ElementPtr)
+ *      Pointer to the element in Row1 to be exchanged.
+ *  Element2  (ElementPtr)
+ *      Pointer to the element in Row2 to be exchanged.
  */
 
-static
-UpdateMarkowitzNumbers( Matrix, pPivot )
+spcRowExchange( Matrix, Row1, Row2 )
 
 MatrixPtr Matrix;
-ElementPtr  pPivot;
+int  Row1, Row2;
 {
-register  int  Row, Col;
-register  ElementPtr  ColPtr, RowPtr;
-register  int *MarkoRow = Matrix->MarkowitzRow, *MarkoCol = Matrix->MarkowitzCol;
-double Product;
-
-/* Begin `UpdateMarkowitzNumbers'. */
+register  ElementPtr  Row1Ptr, Row2Ptr;
+int  Column;
+ElementPtr  Element1, Element2;
 
-/* Update Markowitz numbers. */
-    for (ColPtr = pPivot->NextInCol; ColPtr != NULL; ColPtr = ColPtr->NextInCol)
-    {   Row = ColPtr->Row;
-        --MarkoRow[Row];
+/* Begin `spcRowExchange'. */
+    if (Row1 > Row2)  SWAP(int, Row1, Row2);
 
-/* Form Markowitz product while being cautious of overflows. */
-        if ((MarkoRow[Row] > LARGEST_SHORT_INTEGER AND MarkoCol[Row] != 0) OR
-            (MarkoCol[Row] > LARGEST_SHORT_INTEGER AND MarkoRow[Row] != 0))
-        {   Product = MarkoCol[Row] * MarkoRow[Row];
-            if (Product >= LARGEST_LONG_INTEGER)
-                Matrix->MarkowitzProd[Row] = LARGEST_LONG_INTEGER;
-            else
-                Matrix->MarkowitzProd[Row] = Product;
+    Row1Ptr = Matrix->FirstInRow[Row1];
+    Row2Ptr = Matrix->FirstInRow[Row2];
+    while (Row1Ptr != NULL OR Row2Ptr != NULL)
+    {
+/* Exchange elements in rows while traveling from left to right. */
+        if (Row1Ptr == NULL)
+        {   Column = Row2Ptr->Col;
+            Element1 = NULL;
+            Element2 = Row2Ptr;
+            Row2Ptr = Row2Ptr->NextInRow;
+        }
+        else if (Row2Ptr == NULL)
+        {   Column = Row1Ptr->Col;
+            Element1 = Row1Ptr;
+            Element2 = NULL;
+            Row1Ptr = Row1Ptr->NextInRow;
+        }
+        else if (Row1Ptr->Col < Row2Ptr->Col)
+        {   Column = Row1Ptr->Col;
+            Element1 = Row1Ptr;
+            Element2 = NULL;
+            Row1Ptr = Row1Ptr->NextInRow;
+        }
+        else if (Row1Ptr->Col > Row2Ptr->Col)
+        {   Column = Row2Ptr->Col;
+            Element1 = NULL;
+            Element2 = Row2Ptr;
+            Row2Ptr = Row2Ptr->NextInRow;
+        }
+        else   /* Row1Ptr->Col == Row2Ptr->Col */
+        {   Column = Row1Ptr->Col;
+            Element1 = Row1Ptr;
+            Element2 = Row2Ptr;
+            Row1Ptr = Row1Ptr->NextInRow;
+            Row2Ptr = Row2Ptr->NextInRow;
         }
-        else Matrix->MarkowitzProd[Row] = MarkoRow[Row] * MarkoCol[Row];
-        if (MarkoRow[Row] == 0)
-            Matrix->Singletons++;
-    }
 
-    for (RowPtr = pPivot->NextInRow; RowPtr != NULL; RowPtr = RowPtr->NextInRow)
-    {   Col = RowPtr->Col;
-        --MarkoCol[Col];
+        ExchangeColElements( Matrix, Row1, Element1, Row2, Element2, Column);
+    }  /* end of while(Row1Ptr != NULL OR Row2Ptr != NULL) */
+
+    if (Matrix->InternalVectorsAllocated)
+        SWAP( int, Matrix->MarkowitzRow[Row1], Matrix->MarkowitzRow[Row2]);
+    SWAP( ElementPtr, Matrix->FirstInRow[Row1], Matrix->FirstInRow[Row2]);
+    SWAP( int, Matrix->IntToExtRowMap[Row1], Matrix->IntToExtRowMap[Row2]);
+#if TRANSLATE
+    Matrix->ExtToIntRowMap[ Matrix->IntToExtRowMap[Row1] ] = Row1;
+    Matrix->ExtToIntRowMap[ Matrix->IntToExtRowMap[Row2] ] = Row2;
+#endif
 
-/* Form Markowitz product while being cautious of overflows. */
-        if ((MarkoRow[Col] > LARGEST_SHORT_INTEGER AND MarkoCol[Col] != 0) OR
-            (MarkoCol[Col] > LARGEST_SHORT_INTEGER AND MarkoRow[Col] != 0))
-        {   Product = MarkoCol[Col] * MarkoRow[Col];
-            if (Product >= LARGEST_LONG_INTEGER)
-                Matrix->MarkowitzProd[Col] = LARGEST_LONG_INTEGER;
-            else
-                Matrix->MarkowitzProd[Col] = Product;
-        }
-        else Matrix->MarkowitzProd[Col] = MarkoRow[Col] * MarkoCol[Col];
-        if ((MarkoCol[Col] == 0) AND (MarkoRow[Col] != 0))
-            Matrix->Singletons++;
-    }
     return;
 }
 
@@ -2954,75 +2705,144 @@
 
 
 
+
+
 
 /*
- *  CREATE FILL-IN
- *
- *  This routine is used to create fill-ins and splice them into the
- *  matrix.
+ *  EXCHANGE TWO ELEMENTS IN A ROW
  *
- *  >>> Returns:
- *  Pointer to fill-in.
+ *  Performs all required operations to exchange two elements in a row.
+ *  Those operations are: restring NextInRow pointers and swapping column
+ *  indexes in the MatrixElements.
  *
  *  >>> Arguments:
  *  Matrix  <input>  (MatrixPtr)
  *      Pointer to the matrix.
- *  Col  <input>  (int)
- *      Column index for element.
- *  Row  <input>  (int)
- *      Row index for element.
+ *  Col1  <input>  (int)
+ *      Col of left-most element to be exchanged.
+ *  Element1  <input>  (ElementPtr)
+ *      Pointer to left-most element to be exchanged.
+ *  Col2  <input>  (int)
+ *      Col of right-most element to be exchanged.
+ *  Element2  <input>  (ElementPtr)
+ *      Pointer to right-most element to be exchanged.
+ *  Row <input>  (int)
+ *      Row that exchange is to take place in.
  *
  *  >>> Local variables:
+ *  ElementLeftOfCol1  (ElementPtr *)
+ *      Location of pointer which points to the element to the left of
+ *      Element1. This pointer is modified so that it points to correct
+ *      element on exit.
+ *  ElementLeftOfCol2  (ElementPtr *)
+ *      Location of pointer which points to the element to the left of
+ *      Element2. This pointer is modified so that it points to correct
+ *      element on exit.
+ *  ElementRightOfCol1  (ElementPtr)
+ *      Pointer to element right of Element1.
+ *  ElementRightOfCol2  (ElementPtr)
+ *      Pointer to element right of Element2.
  *  pElement  (ElementPtr)
- *      Pointer to an element in the matrix.
- *  ppElementAbove  (ElementPtr *)
- *      This contains the address of the pointer to the element just above the
- *      one being created. It is used to speed the search and it is updated
- *      with address of the created element.
- *
- *  >>> Possible errors:
- *  spNO_MEMORY
+ *      Pointer used to traverse the row.
  */
 
-static ElementPtr
-CreateFillin( Matrix, Row, Col )
+static
+ExchangeRowElements( Matrix, Col1, Element1, Col2, Element2, Row )
 
 MatrixPtr Matrix;
-register int  Row;
-int  Col;
+int  Col1, Col2, Row;
+register ElementPtr  Element1, Element2;
 {
-register  ElementPtr  pElement, *ppElementAbove;
-ElementPtr  spcCreateElement();
+ElementPtr  *ElementLeftOfCol1, *ElementLeftOfCol2;
+ElementPtr  ElementRightOfCol1, ElementRightOfCol2;
+register   ElementPtr  pElement;
 
-/* Begin `CreateFillin'. */
+/* Begin `ExchangeRowElements'. */
+/* Search to find the ElementLeftOfCol1. */
+    ElementLeftOfCol1 = &(Matrix->FirstInRow[Row]);
+    pElement = *ElementLeftOfCol1;
+    while (pElement->Col < Col1)
+    {   ElementLeftOfCol1 = &(pElement->NextInRow);
+        pElement = *ElementLeftOfCol1;
+    }
+    if (Element1 != NULL)
+    {   ElementRightOfCol1 = Element1->NextInRow;
+        if (Element2 == NULL)
+        {
+/* Element2 does not exist, move Element1 to right to Col2. */
+            if ( ElementRightOfCol1 != NULL AND ElementRightOfCol1->Col < Col2 )
+            {
+/* Element1 must be removed from linked list and moved. */
+                *ElementLeftOfCol1 = ElementRightOfCol1;
 
-/* Find Element above fill-in. */
-    ppElementAbove = &Matrix->FirstInCol[Col];
-    pElement = *ppElementAbove;
-    while (pElement != NULL)
-    {   if (pElement->Row < Row)
-        {   ppElementAbove = &pElement->NextInCol;
-            pElement = *ppElementAbove;
+/* Search Row for Col2. */
+                pElement = ElementRightOfCol1;
+                do
+                {   ElementLeftOfCol2 = &(pElement->NextInRow);
+                    pElement = *ElementLeftOfCol2;
+                }   while (pElement != NULL AND pElement->Col < Col2);
+
+/* Place Element1 in Col2. */
+                *ElementLeftOfCol2 = Element1;
+                Element1->NextInRow = pElement;
+                *ElementLeftOfCol1 =ElementRightOfCol1;
+            }
+            Element1->Col = Col2;
         }
-        else break; /* while loop */
-    }
+        else
+        {
+/* Element2 does exist, and the two elements must be exchanged. */
+            if ( ElementRightOfCol1->Col == Col2)
+            {
+/* Element2 is just right of Element1, exchange them. */
+                Element1->NextInRow = Element2->NextInRow;
+                Element2->NextInRow = Element1;
+                *ElementLeftOfCol1 = Element2;
+            }
+            else
+            {
+/* Element2 is not just right of Element1 and must be searched for. */
+                pElement = ElementRightOfCol1;
+                do
+                {   ElementLeftOfCol2 = &(pElement->NextInRow);
+                    pElement = *ElementLeftOfCol2;
+                }   while (pElement->Col < Col2);
 
-/* End of search, create the element. */
-    pElement = spcCreateElement( Matrix, Row, Col, ppElementAbove, YES );
+                ElementRightOfCol2 = Element2->NextInRow;
 
-/* Update Markowitz counts and products. */
-    Matrix->MarkowitzProd[Row] = ++Matrix->MarkowitzRow[Row] *
-                                   Matrix->MarkowitzCol[Row];
-    if ((Matrix->MarkowitzRow[Row] == 1) AND (Matrix->MarkowitzCol[Row] != 0))
-        Matrix->Singletons--;
-    Matrix->MarkowitzProd[Col] = ++Matrix->MarkowitzCol[Col] *
-                                   Matrix->MarkowitzRow[Col];
-    if ((Matrix->MarkowitzRow[Col] != 0) AND (Matrix->MarkowitzCol[Col] == 1))
-        Matrix->Singletons--;
+/* Switch Element1 and Element2. */
+                *ElementLeftOfCol1 = Element2;
+                Element2->NextInRow = ElementRightOfCol1;
+                *ElementLeftOfCol2 = Element1;
+                Element1->NextInRow = ElementRightOfCol2;
+            }
+            Element1->Col = Col2;
+            Element2->Col = Col1;
+        }
+    }
+    else
+    {
+/* Element1 does not exist. */
+        ElementRightOfCol1 = pElement;
+
+/* Find Element2. */
+        if (ElementRightOfCol1->Col != Col2)
+        {   do
+            {   ElementLeftOfCol2 = &(pElement->NextInRow);
+                pElement = *ElementLeftOfCol2;
+            }   while (pElement->Col < Col2);
 
-    return pElement;
-}
+            ElementRightOfCol2 = Element2->NextInRow;
 
+/* Move Element2 to Col1. */
+            *ElementLeftOfCol2 = Element2->NextInRow;
+            *ElementLeftOfCol1 = Element2;
+            Element2->NextInRow = ElementRightOfCol1;
+        }
+        Element2->Col = Col1;
+    }
+    return;
+}
 
 
 
@@ -3031,46 +2851,96 @@
 
 
 /*
- *  ZERO PIVOT ENCOUNTERED
- *
- *  This routine is called when a singular matrix is found.  It then
- *  records the current row and column and exits.
+ *  EXCHANGE COLUMNS
  *
- *  >>> Returned:
- *  The error code spSINGULAR or spZERO_DIAG is returned.
+ *  Performs all required operations to exchange two columns. Those operations
+ *  include: swap FirstInCol pointers, fixing up the NextInRow pointers,
+ *  swapping column indexes in MatrixElements, and swapping Markowitz 
+ *  column counts.
  *
  *  >>> Arguments:
  *  Matrix  <input>  (MatrixPtr)
- *      Pointer to matrix.
- *  Step  <input>  (int)
- *      Index of diagonal that is zero.
+ *      Pointer to the matrix.
+ *  Col1  <input>  (int)
+ *      Column index of one of the columns, becomes the smallest index.
+ *  Col2  <input>  (int)
+ *      Column index of the other column, becomes the largest index
+ *
+ *  Local variables:
+ *  Row  (int)
+ *      Row in which column elements are currently being exchanged.
+ *  Col1Ptr  (ElementPtr)
+ *      Pointer to an element in Col1.
+ *  Col2Ptr  (ElementPtr)
+ *      Pointer to an element in Col2.
+ *  Element1  (ElementPtr)
+ *      Pointer to the element in Col1 to be exchanged.
+ *  Element2  (ElementPtr)
+ *      Pointer to the element in Col2 to be exchanged.
  */
 
-static int
-MatrixIsSingular( Matrix, Step )
+spcColExchange( Matrix, Col1, Col2 )
 
-MatrixPtr  Matrix;
-int  Step;
+MatrixPtr Matrix;
+int  Col1, Col2;
 {
-/* Begin `MatrixIsSingular'. */
+register  ElementPtr  Col1Ptr, Col2Ptr;
+int  Row;
+ElementPtr  Element1, Element2;
 
-    Matrix->SingularRow = Matrix->IntToExtRowMap[ Step ];
-    Matrix->SingularCol = Matrix->IntToExtColMap[ Step ];
-    return (Matrix->Error = spSINGULAR);
-}
+/* Begin `spcColExchange'. */
+    if (Col1 > Col2)  SWAP(int, Col1, Col2);
 
+    Col1Ptr = Matrix->FirstInCol[Col1];
+    Col2Ptr = Matrix->FirstInCol[Col2];
+    while (Col1Ptr != NULL OR Col2Ptr != NULL)
+    {
+/* Exchange elements in rows while traveling from top to bottom. */
+        if (Col1Ptr == NULL)
+        {   Row = Col2Ptr->Row;
+            Element1 = NULL;
+            Element2 = Col2Ptr;
+            Col2Ptr = Col2Ptr->NextInCol;
+        }
+        else if (Col2Ptr == NULL)
+        {   Row = Col1Ptr->Row;
+            Element1 = Col1Ptr;
+            Element2 = NULL;
+            Col1Ptr = Col1Ptr->NextInCol;
+        }
+        else if (Col1Ptr->Row < Col2Ptr->Row)
+        {   Row = Col1Ptr->Row;
+            Element1 = Col1Ptr;
+            Element2 = NULL;
+            Col1Ptr = Col1Ptr->NextInCol;
+        }
+        else if (Col1Ptr->Row > Col2Ptr->Row)
+        {   Row = Col2Ptr->Row;
+            Element1 = NULL;
+            Element2 = Col2Ptr;
+            Col2Ptr = Col2Ptr->NextInCol;
+        }
+        else   /* Col1Ptr->Row == Col2Ptr->Row */
+        {   Row = Col1Ptr->Row;
+            Element1 = Col1Ptr;
+            Element2 = Col2Ptr;
+            Col1Ptr = Col1Ptr->NextInCol;
+            Col2Ptr = Col2Ptr->NextInCol;
+        }
 
-static int
-ZeroPivot( Matrix, Step )
+        ExchangeRowElements( Matrix, Col1, Element1, Col2, Element2, Row);
+    }  /* end of while(Col1Ptr != NULL OR Col2Ptr != NULL) */
 
-MatrixPtr  Matrix;
-int  Step;
-{
-/* Begin `ZeroPivot'. */
+    if (Matrix->InternalVectorsAllocated)
+        SWAP( int, Matrix->MarkowitzCol[Col1], Matrix->MarkowitzCol[Col2]);
+    SWAP( ElementPtr, Matrix->FirstInCol[Col1], Matrix->FirstInCol[Col2]);
+    SWAP( int, Matrix->IntToExtColMap[Col1], Matrix->IntToExtColMap[Col2]);
+#if TRANSLATE
+    Matrix->ExtToIntColMap[ Matrix->IntToExtColMap[Col1] ] = Col1;
+    Matrix->ExtToIntColMap[ Matrix->IntToExtColMap[Col2] ] = Col2;
+#endif
 
-    Matrix->SingularRow = Matrix->IntToExtRowMap[ Step ];
-    Matrix->SingularCol = Matrix->IntToExtColMap[ Step ];
-    return (Matrix->Error = spZERO_DIAG);
+    return;
 }
 
 
@@ -3078,6 +2948,14 @@
 
 
 
+
+
+
+
+
+
+
+
 #if ANNOTATE == FULL
 
 /*
@@ -3150,3 +3028,133 @@
 
 }
 #endif /* ANNOTATE == FULL */
+
+
+
+
+/*
+ *  FACTOR MATRIX
+ *
+ *  This routine is the companion routine to spOrderAndFactor().
+ *  Unlike spOrderAndFactor(), spFactor() cannot change the ordering.
+ *  It is also faster than spOrderAndFactor().  The standard way of
+ *  using these two routines is to first use spOrderAndFactor() for the
+ *  initial factorization.  For subsequent factorizations, spFactor()
+ *  is used if there is some assurance that little growth will occur
+ *  (say for example, that the matrix is diagonally dominant).  If
+ *  spFactor() is called for the initial factorization of the matrix,
+ *  then spOrderAndFactor() is automatically called with the default
+ *  threshold.  This routine uses "row at a time" LU factorization.
+ *  Pivots are associated with the lower triangular matrix and the
+ *  diagonals of the upper triangular matrix are ones.
+ *
+ *  >>> Returned:
+ *  The error code is returned.  Possible errors are listed below.
+ *
+ *  >>> Arguments:
+ *  Matrix  <input>  (char *)
+ *      Pointer to matrix.
+ *
+ *  >>> Possible errors:
+ *  spNO_MEMORY
+ *  spSINGULAR
+ *  spZERO_DIAG
+ *  spSMALL_PIVOT
+ *  Error is cleared in this function.
+ */
+
+int
+spFactor( eMatrix )
+
+char *eMatrix;
+{
+MatrixPtr  Matrix = (MatrixPtr)eMatrix;
+register  ElementPtr  pElement;
+register  ElementPtr  pColumn;
+register  int  Step, Size;
+RealNumber Mult;
+
+/* Begin `spFactor'. */
+    ASSERT( IS_VALID(Matrix) AND NOT Matrix->Factored);
+
+    if (Matrix->NeedsOrdering)
+    {   return spOrderAndFactor( eMatrix, (RealVector)NULL,
+                                 0.0, 0.0, DIAG_PIVOTING_AS_DEFAULT );
+    }
+    if (NOT Matrix->Partitioned) spPartition( eMatrix, spDEFAULT_PARTITION );
+#if spCOMPLEX
+    if (Matrix->Complex) return FactorComplexMatrix( Matrix );
+#endif
+
+#if REAL
+    Size = Matrix->Size;
+
+    if (Matrix->Diag[1]->Real == 0.0) return ZeroPivot( Matrix, 1 );
+    Matrix->Diag[1]->Real = 1.0 / Matrix->Diag[1]->Real;
+
+/* Start factorization. */
+    for (Step = 2; Step <= Size; Step++)
+    {   if (Matrix->DoRealDirect[Step])
+        {   /* Update column using direct addressing scatter-gather. */
+            register RealNumber *Dest = (RealNumber *)Matrix->Intermediate;
+
+/* Scatter. */
+            pElement = Matrix->FirstInCol[Step];
+            while (pElement != NULL)
+            {   Dest[pElement->Row] = pElement->Real;
+                pElement = pElement->NextInCol;
+            }
+
+/* Update column. */
+            pColumn = Matrix->FirstInCol[Step];
+            while (pColumn->Row < Step)
+            {   pElement = Matrix->Diag[pColumn->Row];
+                pColumn->Real = Dest[pColumn->Row] * pElement->Real;
+                while ((pElement = pElement->NextInCol) != NULL)
+                    Dest[pElement->Row] -= pColumn->Real * pElement->Real;
+                pColumn = pColumn->NextInCol;
+            }
+
+/* Gather. */
+            pElement = Matrix->Diag[Step]->NextInCol;
+            while (pElement != NULL)
+            {   pElement->Real = Dest[pElement->Row];
+                pElement = pElement->NextInCol;
+            }
+
+/* Check for singular matrix. */
+            if (Dest[Step] == 0.0) return ZeroPivot( Matrix, Step );
+            Matrix->Diag[Step]->Real = 1.0 / Dest[Step];
+        }
+        else
+        {   /* Update column using indirect addressing scatter-gather. */
+            register RealNumber **pDest = (RealNumber **)Matrix->Intermediate;
+
+/* Scatter. */
+            pElement = Matrix->FirstInCol[Step];
+            while (pElement != NULL)
+            {   pDest[pElement->Row] = &pElement->Real;
+                pElement = pElement->NextInCol;
+            }
+
+/* Update column. */
+            pColumn = Matrix->FirstInCol[Step];
+            while (pColumn->Row < Step)
+            {   pElement = Matrix->Diag[pColumn->Row];
+                Mult = (*pDest[pColumn->Row] *= pElement->Real);
+                while ((pElement = pElement->NextInCol) != NULL)
+                    *pDest[pElement->Row] -= Mult * pElement->Real;
+                pColumn = pColumn->NextInCol;
+            }
+
+/* Check for singular matrix. */
+            if (Matrix->Diag[Step]->Real == 0.0)
+                return ZeroPivot( Matrix, Step );
+            Matrix->Diag[Step]->Real = 1.0 / Matrix->Diag[Step]->Real;
+        }
+    }
+
+    Matrix->Factored = YES;
+    return (Matrix->Error = spOKAY);
+#endif /* REAL */
+}
diff -uNr fasthenry-3.0.orig/src/fasthenry/sparse/spSolve.c fasthenry-3.0/src/fasthenry/sparse/spSolve.c
--- fasthenry-3.0.orig/src/fasthenry/sparse/spSolve.c	2004-08-03 17:17:19.000000000 -0500
+++ fasthenry-3.0/src/fasthenry/sparse/spSolve.c	2009-08-05 14:07:26.000000000 -0500
@@ -95,11 +95,10 @@
 #include "spDefs.h"
 
 
-
-
 
+#if spCOMPLEX
 /*
- *  SOLVE MATRIX EQUATION
+ *  SOLVE COMPLEX MATRIX EQUATION
  *
  *  Performs forward elimination and back substitution to find the
  *  unknown vector from the RHS vector and factored matrix.  This
@@ -113,30 +112,31 @@
  *  Matrix  <input>  (char *)
  *      Pointer to matrix.
  *  RHS  <input>  (RealVector)
- *      RHS is the input data array, the right hand side. This data is
- *      undisturbed and may be reused for other solves.
+ *      RHS is the real portion of the input data array, the right hand
+ *      side. This data is undisturbed and may be reused for other solves.
  *  Solution  <output>  (RealVector)
- *      Solution is the output data array. This routine is constructed such that
- *      RHS and Solution can be the same array.
+ *      Solution is the real portion of the output data array. This routine
+ *      is constructed such that RHS and Solution can be the same
+ *      array.
  *  iRHS  <input>  (RealVector)
  *      iRHS is the imaginary portion of the input data array, the right
  *      hand side. This data is undisturbed and may be reused for other solves.
- *      This argument is only necessary if matrix is complex and if
- *      spSEPARATED_COMPLEX_VECTOR is set true.
+ *      If spSEPARATED_COMPLEX_VECTOR is set false, there is no need to
+ *      supply this array.
  *  iSolution  <output>  (RealVector)
  *      iSolution is the imaginary portion of the output data array. This
  *      routine is constructed such that iRHS and iSolution can be
- *      the same array.  This argument is only necessary if matrix is complex
- *      and if spSEPARATED_COMPLEX_VECTOR is set true.
+ *      the same array.  If spSEPARATED_COMPLEX_VECTOR is set false, there is no
+ *      need to supply this array.
  *
  *  >>> Local variables:
- *  Intermediate  (RealVector)
+ *  Intermediate  (ComplexVector)
  *      Temporary storage for use in forward elimination and backward
  *      substitution.  Commonly referred to as c, when the LU factorization
- *      equations are given as  Ax = b, Lc = b, Ux = c Local version of
- *      Matrix->Intermediate, which was created during the initial
- *      factorization in function CreateInternalVectors() in the matrix
- *      factorization module.
+ *      equations are given as  Ax = b, Lc = b, Ux = c.
+ *      Local version of Matrix->Intermediate, which was created during
+ *      the initial factorization in function CreateInternalVectors() in the
+ *      matrix factorization module.
  *  pElement  (ElementPtr)
  *      Pointer used to address elements in both the lower and upper triangle
  *      matrices.
@@ -148,7 +148,7 @@
  *      Pointer that points to current pivot or diagonal element.
  *  Size  (int)
  *      Size of matrix. Made local to reduce indirection.
- *  Temp  (RealNumber)
+ *  Temp  (ComplexNumber)
  *      Temporary storage for entries in arrays.
  *
  *  >>> Obscure Macros
@@ -158,58 +158,64 @@
  *      without a trace.
  */
 
-/*VARARGS3*/
-
-void
-spSolve( eMatrix, RHS, Solution IMAG_VECTORS )
+static void
+SolveComplexMatrix( Matrix, RHS, Solution IMAG_VECTORS )
 
-char *eMatrix;
+MatrixPtr  Matrix;
 RealVector  RHS, Solution IMAG_VECTORS;
 {
-MatrixPtr  Matrix = (MatrixPtr)eMatrix;
 register  ElementPtr  pElement;
-register  RealVector  Intermediate;
-register  RealNumber  Temp;
+register  ComplexVector  Intermediate;
 register  int  I, *pExtOrder, Size;
 ElementPtr  pPivot;
-void SolveComplexMatrix();
-
-/* Begin `spSolve'. */
-    ASSERT( IS_VALID(Matrix) AND IS_FACTORED(Matrix) );
+register ComplexVector  ExtVector;
+ComplexNumber  Temp;
 
-#if spCOMPLEX
-    if (Matrix->Complex)
-    {   SolveComplexMatrix( Matrix, RHS, Solution IMAG_VECTORS );
-        return;
-    }
-#endif
+/* Begin `SolveComplexMatrix'. */
 
-#if REAL
-    Intermediate = Matrix->Intermediate;
     Size = Matrix->Size;
+    Intermediate = (ComplexVector)Matrix->Intermediate;
 
 /* Correct array pointers for ARRAY_OFFSET. */
 #if NOT ARRAY_OFFSET
-    --RHS;
-    --Solution;
+#if spSEPARATED_COMPLEX_VECTORS
+    --RHS;      --iRHS;
+    --Solution; --iSolution;
+#else
+    RHS -= 2; Solution -= 2;
+#endif
 #endif
 
 /* Initialize Intermediate vector. */
     pExtOrder = &Matrix->IntToExtRowMap[Size];
+
+#if spSEPARATED_COMPLEX_VECTORS
     for (I = Size; I > 0; I--)
-        Intermediate[I] = RHS[*(pExtOrder--)];
+    {   Intermediate[I].Real = RHS[*(pExtOrder)];
+        Intermediate[I].Imag = iRHS[*(pExtOrder--)];
+    }
+#else
+    ExtVector = (ComplexVector)RHS;
+    for (I = Size; I > 0; I--)
+        Intermediate[I] = ExtVector[*(pExtOrder--)];
+#endif
 
-/* Forward elimination. Solves Lc = b.*/
+/* Forward substitution. Solves Lc = b.*/
     for (I = 1; I <= Size; I++)
-    {   
-/* This step of the elimination is skipped if Temp equals zero. */
-        if ((Temp = Intermediate[I]) != 0.0)
-        {   pPivot = Matrix->Diag[I];
-            Intermediate[I] = (Temp *= pPivot->Real);
+    {   Temp = Intermediate[I];
 
+/* This step of the substitution is skipped if Temp equals zero. */
+        if ((Temp.Real != 0.0) OR (Temp.Imag != 0.0))
+        {   pPivot = Matrix->Diag[I];
+/* Cmplx expr: Temp *= (1.0 / Pivot). */
+            CMPLX_MULT_ASSIGN(Temp, *pPivot);
+            Intermediate[I] = Temp;
             pElement = pPivot->NextInCol;
             while (pElement != NULL)
-            {   Intermediate[pElement->Row] -= Temp * pElement->Real;
+            {
+/* Cmplx expr: Intermediate[Element->Row] -= Temp * *Element. */
+                CMPLX_MULT_SUBT_ASSIGN(Intermediate[pElement->Row],
+                                       Temp, *pElement);
                 pElement = pElement->NextInCol;
             }
         }
@@ -219,8 +225,11 @@
     for (I = Size; I > 0; I--)
     {   Temp = Intermediate[I];
         pElement = Matrix->Diag[I]->NextInRow;
+
         while (pElement != NULL)
-        {   Temp -= pElement->Real * Intermediate[pElement->Col];
+        {
+/* Cmplx expr: Temp -= *Element * Intermediate[Element->Col]. */
+            CMPLX_MULT_SUBT_ASSIGN(Temp, *pElement,Intermediate[pElement->Col]);
             pElement = pElement->NextInRow;
         }
         Intermediate[I] = Temp;
@@ -228,14 +237,21 @@
 
 /* Unscramble Intermediate vector while placing data in to Solution vector. */
     pExtOrder = &Matrix->IntToExtColMap[Size];
+
+#if spSEPARATED_COMPLEX_VECTORS
     for (I = Size; I > 0; I--)
-        Solution[*(pExtOrder--)] = Intermediate[I];
+    {   Solution[*(pExtOrder)] = Intermediate[I].Real;
+        iSolution[*(pExtOrder--)] = Intermediate[I].Imag;
+    }
+#else
+    ExtVector = (ComplexVector)Solution;
+    for (I = Size; I > 0; I--)
+        ExtVector[*(pExtOrder--)] = Intermediate[I];
+#endif
 
     return;
-#endif /* REAL */
 }
-
-
+#endif /* spCOMPLEX */
 
 
 
@@ -245,9 +261,8 @@
 
 
 
-#if spCOMPLEX
 /*
- *  SOLVE COMPLEX MATRIX EQUATION
+ *  SOLVE MATRIX EQUATION
  *
  *  Performs forward elimination and back substitution to find the
  *  unknown vector from the RHS vector and factored matrix.  This
@@ -261,31 +276,30 @@
  *  Matrix  <input>  (char *)
  *      Pointer to matrix.
  *  RHS  <input>  (RealVector)
- *      RHS is the real portion of the input data array, the right hand
- *      side. This data is undisturbed and may be reused for other solves.
+ *      RHS is the input data array, the right hand side. This data is
+ *      undisturbed and may be reused for other solves.
  *  Solution  <output>  (RealVector)
- *      Solution is the real portion of the output data array. This routine
- *      is constructed such that RHS and Solution can be the same
- *      array.
+ *      Solution is the output data array. This routine is constructed such that
+ *      RHS and Solution can be the same array.
  *  iRHS  <input>  (RealVector)
  *      iRHS is the imaginary portion of the input data array, the right
  *      hand side. This data is undisturbed and may be reused for other solves.
- *      If spSEPARATED_COMPLEX_VECTOR is set false, there is no need to
- *      supply this array.
+ *      This argument is only necessary if matrix is complex and if
+ *      spSEPARATED_COMPLEX_VECTOR is set true.
  *  iSolution  <output>  (RealVector)
  *      iSolution is the imaginary portion of the output data array. This
  *      routine is constructed such that iRHS and iSolution can be
- *      the same array.  If spSEPARATED_COMPLEX_VECTOR is set false, there is no
- *      need to supply this array.
+ *      the same array.  This argument is only necessary if matrix is complex
+ *      and if spSEPARATED_COMPLEX_VECTOR is set true.
  *
  *  >>> Local variables:
- *  Intermediate  (ComplexVector)
+ *  Intermediate  (RealVector)
  *      Temporary storage for use in forward elimination and backward
  *      substitution.  Commonly referred to as c, when the LU factorization
- *      equations are given as  Ax = b, Lc = b, Ux = c.
- *      Local version of Matrix->Intermediate, which was created during
- *      the initial factorization in function CreateInternalVectors() in the
- *      matrix factorization module.
+ *      equations are given as  Ax = b, Lc = b, Ux = c Local version of
+ *      Matrix->Intermediate, which was created during the initial
+ *      factorization in function CreateInternalVectors() in the matrix
+ *      factorization module.
  *  pElement  (ElementPtr)
  *      Pointer used to address elements in both the lower and upper triangle
  *      matrices.
@@ -297,7 +311,7 @@
  *      Pointer that points to current pivot or diagonal element.
  *  Size  (int)
  *      Size of matrix. Made local to reduce indirection.
- *  Temp  (ComplexNumber)
+ *  Temp  (RealNumber)
  *      Temporary storage for entries in arrays.
  *
  *  >>> Obscure Macros
@@ -307,64 +321,58 @@
  *      without a trace.
  */
 
-static void
-SolveComplexMatrix( Matrix, RHS, Solution IMAG_VECTORS )
+/*VARARGS3*/
 
-MatrixPtr  Matrix;
+void
+spSolve( eMatrix, RHS, Solution IMAG_VECTORS )
+
+char *eMatrix;
 RealVector  RHS, Solution IMAG_VECTORS;
 {
+MatrixPtr  Matrix = (MatrixPtr)eMatrix;
 register  ElementPtr  pElement;
-register  ComplexVector  Intermediate;
+register  RealVector  Intermediate;
+register  RealNumber  Temp;
 register  int  I, *pExtOrder, Size;
 ElementPtr  pPivot;
-register ComplexVector  ExtVector;
-ComplexNumber  Temp;
+void SolveComplexMatrix();
 
-/* Begin `SolveComplexMatrix'. */
+/* Begin `spSolve'. */
+    ASSERT( IS_VALID(Matrix) AND IS_FACTORED(Matrix) );
 
+#if spCOMPLEX
+    if (Matrix->Complex)
+    {   SolveComplexMatrix( Matrix, RHS, Solution IMAG_VECTORS );
+        return;
+    }
+#endif
+
+#if REAL
+    Intermediate = Matrix->Intermediate;
     Size = Matrix->Size;
-    Intermediate = (ComplexVector)Matrix->Intermediate;
 
 /* Correct array pointers for ARRAY_OFFSET. */
 #if NOT ARRAY_OFFSET
-#if spSEPARATED_COMPLEX_VECTORS
-    --RHS;      --iRHS;
-    --Solution; --iSolution;
-#else
-    RHS -= 2; Solution -= 2;
-#endif
+    --RHS;
+    --Solution;
 #endif
 
 /* Initialize Intermediate vector. */
     pExtOrder = &Matrix->IntToExtRowMap[Size];
-
-#if spSEPARATED_COMPLEX_VECTORS
-    for (I = Size; I > 0; I--)
-    {   Intermediate[I].Real = RHS[*(pExtOrder)];
-        Intermediate[I].Imag = iRHS[*(pExtOrder--)];
-    }
-#else
-    ExtVector = (ComplexVector)RHS;
     for (I = Size; I > 0; I--)
-        Intermediate[I] = ExtVector[*(pExtOrder--)];
-#endif
+        Intermediate[I] = RHS[*(pExtOrder--)];
 
-/* Forward substitution. Solves Lc = b.*/
+/* Forward elimination. Solves Lc = b.*/
     for (I = 1; I <= Size; I++)
-    {   Temp = Intermediate[I];
-
-/* This step of the substitution is skipped if Temp equals zero. */
-        if ((Temp.Real != 0.0) OR (Temp.Imag != 0.0))
+    {   
+/* This step of the elimination is skipped if Temp equals zero. */
+        if ((Temp = Intermediate[I]) != 0.0)
         {   pPivot = Matrix->Diag[I];
-/* Cmplx expr: Temp *= (1.0 / Pivot). */
-            CMPLX_MULT_ASSIGN(Temp, *pPivot);
-            Intermediate[I] = Temp;
+            Intermediate[I] = (Temp *= pPivot->Real);
+
             pElement = pPivot->NextInCol;
             while (pElement != NULL)
-            {
-/* Cmplx expr: Intermediate[Element->Row] -= Temp * *Element. */
-                CMPLX_MULT_SUBT_ASSIGN(Intermediate[pElement->Row],
-                                       Temp, *pElement);
+            {   Intermediate[pElement->Row] -= Temp * pElement->Real;
                 pElement = pElement->NextInCol;
             }
         }
@@ -374,11 +382,8 @@
     for (I = Size; I > 0; I--)
     {   Temp = Intermediate[I];
         pElement = Matrix->Diag[I]->NextInRow;
-
         while (pElement != NULL)
-        {
-/* Cmplx expr: Temp -= *Element * Intermediate[Element->Col]. */
-            CMPLX_MULT_SUBT_ASSIGN(Temp, *pElement,Intermediate[pElement->Col]);
+        {   Temp -= pElement->Real * Intermediate[pElement->Col];
             pElement = pElement->NextInRow;
         }
         Intermediate[I] = Temp;
@@ -386,21 +391,16 @@
 
 /* Unscramble Intermediate vector while placing data in to Solution vector. */
     pExtOrder = &Matrix->IntToExtColMap[Size];
-
-#if spSEPARATED_COMPLEX_VECTORS
     for (I = Size; I > 0; I--)
-    {   Solution[*(pExtOrder)] = Intermediate[I].Real;
-        iSolution[*(pExtOrder--)] = Intermediate[I].Imag;
-    }
-#else
-    ExtVector = (ComplexVector)Solution;
-    for (I = Size; I > 0; I--)
-        ExtVector[*(pExtOrder--)] = Intermediate[I];
-#endif
+        Solution[*(pExtOrder--)] = Intermediate[I];
 
     return;
+#endif /* REAL */
 }
-#endif /* spCOMPLEX */
+
+
+
+
 
 
 
diff -uNr fasthenry-3.0.orig/src/fasthenry/sparse/spUtils.c fasthenry-3.0/src/fasthenry/sparse/spUtils.c
--- fasthenry-3.0.orig/src/fasthenry/sparse/spUtils.c	2004-08-03 17:17:20.000000000 -0500
+++ fasthenry-3.0/src/fasthenry/sparse/spUtils.c	2009-08-05 14:11:45.000000000 -0500
@@ -1419,47 +1419,20 @@
 
 
 #if CONDITION
+#define SLACK   1e4
+
+
+#if spCOMPLEX
 /*
  *  ESTIMATE CONDITION NUMBER
  *
- *  Computes an estimate of the condition number using a variation on
- *  the LINPACK condition number estimation algorithm.  This quantity is
- *  an indicator of ill-conditioning in the matrix.  To avoid problems
- *  with overflow, the reciprocal of the condition number is returned.
- *  If this number is small, and if the matrix is scaled such that
- *  uncertainties in the RHS and the matrix entries are equilibrated,
- *  then the matrix is ill-conditioned.  If the this number is near
- *  one, the matrix is well conditioned.  This routine must only be
- *  used after a matrix has been factored by spOrderAndFactor() or
- *  spFactor() and before it is cleared by spClear() or spInitialize().
- *
- *  Unlike the LINPACK condition number estimator, this routines
- *  returns the L infinity condition number.  This is an artifact of
- *  Sparse placing ones on the diagonal of the upper triangular matrix
- *  rather than the lower.  This difference should be of no importance.
- *
- *  References:
- *  A.K. Cline, C.B. Moler, G.W. Stewart, J.H. Wilkinson.  An estimate
- *  for the condition number of a matrix.  SIAM Journal on Numerical
- *  Analysis.  Vol. 16, No. 2, pages 368-375, April 1979.
- *  
- *  J.J. Dongarra, C.B. Moler, J.R. Bunch, G.W. Stewart.  LINPACK
- *  User's Guide.  SIAM, 1979.
- *  
- *  Roger G. Grimes, John G. Lewis.  Condition number estimation for
- *  sparse matrices.  SIAM Journal on Scientific and Statistical
- *  Computing.  Vol. 2, No. 4, pages 384-388, December 1981.
- *  
- *  Dianne Prost O'Leary.  Estimating matrix condition numbers.  SIAM
- *  Journal on Scientific and Statistical Computing.  Vol. 1, No. 2,
- *  pages 205-209, June 1980.
+ *  Complex version of spCondition().
  *
  *  >>> Returns:
- *  The reciprocal of the condition number.  If the matrix was singular,
- *  zero is returned.
+ *  The reciprocal of the condition number.
  *
  *  >>> Arguments:
- *  eMatrix  <input>  (char *)
+ *  Matrix  <input>  (MatrixPtr)
  *      Pointer to the matrix.
  *  NormOfMatrix  <input>  (RealNumber)
  *      The L-infinity norm of the unfactored matrix as computed by
@@ -1468,55 +1441,35 @@
  *      Used to return error code.
  *
  *  >>> Possible errors:
- *  spSINGULAR
  *  spNO_MEMORY
  */
 
-RealNumber
-spCondition( eMatrix, NormOfMatrix, pError )
+static RealNumber
+ComplexCondition( Matrix, NormOfMatrix, pError )
 
-char *eMatrix;
+MatrixPtr Matrix;
 RealNumber NormOfMatrix;
 int *pError;
 {
-MatrixPtr  Matrix = (MatrixPtr)eMatrix;
 register ElementPtr pElement;
-register RealVector T, Tm;
+register ComplexVector T, Tm;
 register int I, K, Row;
 ElementPtr pPivot;
 int Size;
-RealNumber E, Em, Wp, Wm, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor;
-RealNumber Linpack, OLeary, InvNormOfInverse, ComplexCondition();
-#define SLACK   1e4
-
-/* Begin `spCondition'. */
-
-    ASSERT( IS_SPARSE(Matrix) AND IS_FACTORED(Matrix) );
-    *pError = Matrix->Error;
-    if (Matrix->Error >= spFATAL) return 0.0;
-    if (NormOfMatrix == 0.0)
-    {   *pError = spSINGULAR;
-        return 0.0;
-    }
+RealNumber E, Em, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor;
+RealNumber Linpack, OLeary, InvNormOfInverse;
+ComplexNumber Wp, Wm;
 
-#if spCOMPLEX
-    if (Matrix->Complex)
-        return ComplexCondition( Matrix, NormOfMatrix, pError );
-#endif
+/* Begin `ComplexCondition'. */
 
-#if REAL
     Size = Matrix->Size;
-    T = Matrix->Intermediate;
-#if spCOMPLEX
-    Tm = Matrix->Intermediate + Size;
-#else
-    Tm = ALLOC( RealNumber, Size+1 );
+    T = (ComplexVector)Matrix->Intermediate;
+    Tm = ALLOC( ComplexNumber, Size+1 );
     if (Tm == NULL)
     {   *pError = spNO_MEMORY;
         return 0.0;
     }
-#endif
-    for (I = Size; I > 0; I--) T[I] = 0.0;
+    for (I = Size; I > 0; I--) T[I].Real = T[I].Imag = 0.0;
 
 /*
  * Part 1.  Ay = e.
@@ -1529,27 +1482,34 @@
     E = 1.0;
     for (I = 1; I <= Size; I++)
     {   pPivot = Matrix->Diag[I];
-        if (T[I] < 0.0) Em = -E; else Em = E;
-        Wm = (Em + T[I]) * pPivot->Real;
-        if (ABS(Wm) > SLACK)
-        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(Wm) );
-            for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
+        if (T[I].Real < 0.0) Em = -E; else Em = E;
+        Wm = T[I];
+        Wm.Real += Em;
+        ASm = CMPLX_1_NORM( Wm );
+        CMPLX_MULT_ASSIGN( Wm, *pPivot );
+        if (CMPLX_1_NORM(Wm) > SLACK)
+        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(Wm) );
+            for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
             E *= ScaleFactor;
             Em *= ScaleFactor;
-            Wm = (Em + T[I]) * pPivot->Real;
+            ASm *= ScaleFactor;
+            SCLR_MULT_ASSIGN( Wm, ScaleFactor );
         }
-        Wp = (T[I] - Em) * pPivot->Real;
-        ASp = ABS(T[I] - Em);
-        ASm = ABS(Em + T[I]);
+        Wp = T[I];
+        Wp.Real -= Em;
+        ASp = CMPLX_1_NORM( Wp );
+        CMPLX_MULT_ASSIGN( Wp, *pPivot );
 
 /* Update T for both values of W, minus value is placed in Tm. */
         pElement = pPivot->NextInCol;
         while (pElement != NULL)
         {   Row = pElement->Row;
-            Tm[Row] = T[Row] - (Wm * pElement->Real);
-            T[Row] -= (Wp * pElement->Real);
-            ASp += ABS(T[Row]);
-            ASm += ABS(Tm[Row]);
+            /* Cmplx expr: Tm[Row] = T[Row] - (Wp * *pElement). */
+            CMPLX_MULT_SUBT( Tm[Row], Wm, *pElement, T[Row] );
+            /* Cmplx expr: T[Row] -= Wp * *pElement. */
+            CMPLX_MULT_SUBT_ASSIGN( T[Row], Wm, *pElement );
+            ASp += CMPLX_1_NORM(T[Row]);
+            ASm += CMPLX_1_NORM(Tm[Row]);
             pElement = pElement->NextInCol;
         }
 
@@ -1566,10 +1526,10 @@
     }
 
 /* Compute 1-norm of T, which now contains w, and scale ||T|| to 1/SLACK. */
-    for (ASw = 0.0, I = Size; I > 0; I--) ASw += ABS(T[I]);
+    for (ASw = 0.0, I = Size; I > 0; I--) ASw += CMPLX_1_NORM(T[I]);
     ScaleFactor = 1.0 / (SLACK * ASw);
     if (ScaleFactor < 0.5)
-    {   for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
+    {   for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
         E *= ScaleFactor;
     }
 
@@ -1577,28 +1537,29 @@
     for (I = Size; I >= 1; I--)
     {   pElement = Matrix->Diag[I]->NextInRow;
         while (pElement != NULL)
-        {   T[I] -= pElement->Real * T[pElement->Col];
+        {   /* Cmplx expr: T[I] -= T[pElement->Col] * *pElement. */
+            CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Col], *pElement );
             pElement = pElement->NextInRow;
         }
-        if (ABS(T[I]) > SLACK)
-        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) );
-            for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
+        if (CMPLX_1_NORM(T[I]) > SLACK)
+        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
+            for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
             E *= ScaleFactor;
         }
     }
 
 /* Compute 1-norm of T, which now contains y, and scale ||T|| to 1/SLACK. */
-    for (ASy = 0.0, I = Size; I > 0; I--) ASy += ABS(T[I]);
+    for (ASy = 0.0, I = Size; I > 0; I--) ASy += CMPLX_1_NORM(T[I]);
     ScaleFactor = 1.0 / (SLACK * ASy);
     if (ScaleFactor < 0.5)
-    {   for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
+    {   for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
         ASy = 1.0 / SLACK;
         E *= ScaleFactor;
     }
 
 /* Compute infinity-norm of T for O'Leary's estimate. */
     for (MaxY = 0.0, I = Size; I > 0; I--)
-        if (MaxY < ABS(T[I])) MaxY = ABS(T[I]);
+        if (MaxY < CMPLX_1_NORM(T[I])) MaxY = CMPLX_1_NORM(T[I]);
 
 /*
  * Part 2.  A* z = y where the * represents the transpose.
@@ -1609,21 +1570,22 @@
     for (I = 1; I <= Size; I++)
     {   pElement = Matrix->Diag[I]->NextInRow;
         while (pElement != NULL)
-        {   T[pElement->Col] -= T[I] * pElement->Real;
+        {   /* Cmplx expr: T[pElement->Col] -= T[I] * *pElement. */
+            CMPLX_MULT_SUBT_ASSIGN( T[pElement->Col], T[I], *pElement );
             pElement = pElement->NextInRow;
         }
-        if (ABS(T[I]) > SLACK)
-        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) );
-            for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
+        if (CMPLX_1_NORM(T[I]) > SLACK)
+        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
+            for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
             ASy *= ScaleFactor;
         }
     }
 
 /* Compute 1-norm of T, which now contains v, and scale ||T|| to 1/SLACK. */
-    for (ASv = 0.0, I = Size; I > 0; I--) ASv += ABS(T[I]);
+    for (ASv = 0.0, I = Size; I > 0; I--) ASv += CMPLX_1_NORM(T[I]);
     ScaleFactor = 1.0 / (SLACK * ASv);
     if (ScaleFactor < 0.5)
-    {   for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
+    {   for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
         ASy *= ScaleFactor;
     }
 
@@ -1632,46 +1594,73 @@
     {   pPivot = Matrix->Diag[I];
         pElement = pPivot->NextInCol;
         while (pElement != NULL)
-        {   T[I] -= pElement->Real * T[pElement->Row];
+        {   /* Cmplx expr: T[I] -= T[pElement->Row] * *pElement. */
+            CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Row], *pElement );
             pElement = pElement->NextInCol;
         }
-        T[I] *= pPivot->Real;
-        if (ABS(T[I]) > SLACK)
-        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) );
-            for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
+        CMPLX_MULT_ASSIGN( T[I], *pPivot );
+        if (CMPLX_1_NORM(T[I]) > SLACK)
+        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
+            for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
             ASy *= ScaleFactor;
         }
     }
 
 /* Compute 1-norm of T, which now contains z. */
-    for (ASz = 0.0, I = Size; I > 0; I--) ASz += ABS(T[I]);
+    for (ASz = 0.0, I = Size; I > 0; I--) ASz += CMPLX_1_NORM(T[I]);
 
-#if NOT spCOMPLEX
     FREE( Tm );
-#endif
 
     Linpack = ASy / ASz;
     OLeary = E / MaxY;
     InvNormOfInverse = MIN( Linpack, OLeary );
     return InvNormOfInverse / NormOfMatrix;
-#endif /* REAL */
 }
+#endif /* spCOMPLEX */
 
 
 
-
-
-#if spCOMPLEX
 /*
  *  ESTIMATE CONDITION NUMBER
  *
- *  Complex version of spCondition().
+ *  Computes an estimate of the condition number using a variation on
+ *  the LINPACK condition number estimation algorithm.  This quantity is
+ *  an indicator of ill-conditioning in the matrix.  To avoid problems
+ *  with overflow, the reciprocal of the condition number is returned.
+ *  If this number is small, and if the matrix is scaled such that
+ *  uncertainties in the RHS and the matrix entries are equilibrated,
+ *  then the matrix is ill-conditioned.  If the this number is near
+ *  one, the matrix is well conditioned.  This routine must only be
+ *  used after a matrix has been factored by spOrderAndFactor() or
+ *  spFactor() and before it is cleared by spClear() or spInitialize().
+ *
+ *  Unlike the LINPACK condition number estimator, this routines
+ *  returns the L infinity condition number.  This is an artifact of
+ *  Sparse placing ones on the diagonal of the upper triangular matrix
+ *  rather than the lower.  This difference should be of no importance.
+ *
+ *  References:
+ *  A.K. Cline, C.B. Moler, G.W. Stewart, J.H. Wilkinson.  An estimate
+ *  for the condition number of a matrix.  SIAM Journal on Numerical
+ *  Analysis.  Vol. 16, No. 2, pages 368-375, April 1979.
+ *  
+ *  J.J. Dongarra, C.B. Moler, J.R. Bunch, G.W. Stewart.  LINPACK
+ *  User's Guide.  SIAM, 1979.
+ *  
+ *  Roger G. Grimes, John G. Lewis.  Condition number estimation for
+ *  sparse matrices.  SIAM Journal on Scientific and Statistical
+ *  Computing.  Vol. 2, No. 4, pages 384-388, December 1981.
+ *  
+ *  Dianne Prost O'Leary.  Estimating matrix condition numbers.  SIAM
+ *  Journal on Scientific and Statistical Computing.  Vol. 1, No. 2,
+ *  pages 205-209, June 1980.
  *
  *  >>> Returns:
- *  The reciprocal of the condition number.
+ *  The reciprocal of the condition number.  If the matrix was singular,
+ *  zero is returned.
  *
  *  >>> Arguments:
- *  Matrix  <input>  (MatrixPtr)
+ *  eMatrix  <input>  (char *)
  *      Pointer to the matrix.
  *  NormOfMatrix  <input>  (RealNumber)
  *      The L-infinity norm of the unfactored matrix as computed by
@@ -1680,35 +1669,55 @@
  *      Used to return error code.
  *
  *  >>> Possible errors:
+ *  spSINGULAR
  *  spNO_MEMORY
  */
 
-static RealNumber
-ComplexCondition( Matrix, NormOfMatrix, pError )
+RealNumber
+spCondition( eMatrix, NormOfMatrix, pError )
 
-MatrixPtr Matrix;
+char *eMatrix;
 RealNumber NormOfMatrix;
 int *pError;
 {
+MatrixPtr  Matrix = (MatrixPtr)eMatrix;
 register ElementPtr pElement;
-register ComplexVector T, Tm;
+register RealVector T, Tm;
 register int I, K, Row;
 ElementPtr pPivot;
 int Size;
-RealNumber E, Em, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor;
-RealNumber Linpack, OLeary, InvNormOfInverse;
-ComplexNumber Wp, Wm;
+RealNumber E, Em, Wp, Wm, ASp, ASm, ASw, ASy, ASv, ASz, MaxY, ScaleFactor;
+RealNumber Linpack, OLeary, InvNormOfInverse, ComplexCondition();
 
-/* Begin `ComplexCondition'. */
 
+/* Begin `spCondition'. */
+
+    ASSERT( IS_SPARSE(Matrix) AND IS_FACTORED(Matrix) );
+    *pError = Matrix->Error;
+    if (Matrix->Error >= spFATAL) return 0.0;
+    if (NormOfMatrix == 0.0)
+    {   *pError = spSINGULAR;
+        return 0.0;
+    }
+
+#if spCOMPLEX
+    if (Matrix->Complex)
+        return ComplexCondition( Matrix, NormOfMatrix, pError );
+#endif
+
+#if REAL
     Size = Matrix->Size;
-    T = (ComplexVector)Matrix->Intermediate;
-    Tm = ALLOC( ComplexNumber, Size+1 );
+    T = Matrix->Intermediate;
+#if spCOMPLEX
+    Tm = Matrix->Intermediate + Size;
+#else
+    Tm = ALLOC( RealNumber, Size+1 );
     if (Tm == NULL)
     {   *pError = spNO_MEMORY;
         return 0.0;
     }
-    for (I = Size; I > 0; I--) T[I].Real = T[I].Imag = 0.0;
+#endif
+    for (I = Size; I > 0; I--) T[I] = 0.0;
 
 /*
  * Part 1.  Ay = e.
@@ -1721,34 +1730,27 @@
     E = 1.0;
     for (I = 1; I <= Size; I++)
     {   pPivot = Matrix->Diag[I];
-        if (T[I].Real < 0.0) Em = -E; else Em = E;
-        Wm = T[I];
-        Wm.Real += Em;
-        ASm = CMPLX_1_NORM( Wm );
-        CMPLX_MULT_ASSIGN( Wm, *pPivot );
-        if (CMPLX_1_NORM(Wm) > SLACK)
-        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(Wm) );
-            for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
+        if (T[I] < 0.0) Em = -E; else Em = E;
+        Wm = (Em + T[I]) * pPivot->Real;
+        if (ABS(Wm) > SLACK)
+        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(Wm) );
+            for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
             E *= ScaleFactor;
             Em *= ScaleFactor;
-            ASm *= ScaleFactor;
-            SCLR_MULT_ASSIGN( Wm, ScaleFactor );
+            Wm = (Em + T[I]) * pPivot->Real;
         }
-        Wp = T[I];
-        Wp.Real -= Em;
-        ASp = CMPLX_1_NORM( Wp );
-        CMPLX_MULT_ASSIGN( Wp, *pPivot );
+        Wp = (T[I] - Em) * pPivot->Real;
+        ASp = ABS(T[I] - Em);
+        ASm = ABS(Em + T[I]);
 
 /* Update T for both values of W, minus value is placed in Tm. */
         pElement = pPivot->NextInCol;
         while (pElement != NULL)
         {   Row = pElement->Row;
-            /* Cmplx expr: Tm[Row] = T[Row] - (Wp * *pElement). */
-            CMPLX_MULT_SUBT( Tm[Row], Wm, *pElement, T[Row] );
-            /* Cmplx expr: T[Row] -= Wp * *pElement. */
-            CMPLX_MULT_SUBT_ASSIGN( T[Row], Wm, *pElement );
-            ASp += CMPLX_1_NORM(T[Row]);
-            ASm += CMPLX_1_NORM(Tm[Row]);
+            Tm[Row] = T[Row] - (Wm * pElement->Real);
+            T[Row] -= (Wp * pElement->Real);
+            ASp += ABS(T[Row]);
+            ASm += ABS(Tm[Row]);
             pElement = pElement->NextInCol;
         }
 
@@ -1765,10 +1767,10 @@
     }
 
 /* Compute 1-norm of T, which now contains w, and scale ||T|| to 1/SLACK. */
-    for (ASw = 0.0, I = Size; I > 0; I--) ASw += CMPLX_1_NORM(T[I]);
+    for (ASw = 0.0, I = Size; I > 0; I--) ASw += ABS(T[I]);
     ScaleFactor = 1.0 / (SLACK * ASw);
     if (ScaleFactor < 0.5)
-    {   for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
+    {   for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
         E *= ScaleFactor;
     }
 
@@ -1776,29 +1778,28 @@
     for (I = Size; I >= 1; I--)
     {   pElement = Matrix->Diag[I]->NextInRow;
         while (pElement != NULL)
-        {   /* Cmplx expr: T[I] -= T[pElement->Col] * *pElement. */
-            CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Col], *pElement );
+        {   T[I] -= pElement->Real * T[pElement->Col];
             pElement = pElement->NextInRow;
         }
-        if (CMPLX_1_NORM(T[I]) > SLACK)
-        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
-            for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
+        if (ABS(T[I]) > SLACK)
+        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) );
+            for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
             E *= ScaleFactor;
         }
     }
 
 /* Compute 1-norm of T, which now contains y, and scale ||T|| to 1/SLACK. */
-    for (ASy = 0.0, I = Size; I > 0; I--) ASy += CMPLX_1_NORM(T[I]);
+    for (ASy = 0.0, I = Size; I > 0; I--) ASy += ABS(T[I]);
     ScaleFactor = 1.0 / (SLACK * ASy);
     if (ScaleFactor < 0.5)
-    {   for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
+    {   for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
         ASy = 1.0 / SLACK;
         E *= ScaleFactor;
     }
 
 /* Compute infinity-norm of T for O'Leary's estimate. */
     for (MaxY = 0.0, I = Size; I > 0; I--)
-        if (MaxY < CMPLX_1_NORM(T[I])) MaxY = CMPLX_1_NORM(T[I]);
+        if (MaxY < ABS(T[I])) MaxY = ABS(T[I]);
 
 /*
  * Part 2.  A* z = y where the * represents the transpose.
@@ -1809,22 +1810,21 @@
     for (I = 1; I <= Size; I++)
     {   pElement = Matrix->Diag[I]->NextInRow;
         while (pElement != NULL)
-        {   /* Cmplx expr: T[pElement->Col] -= T[I] * *pElement. */
-            CMPLX_MULT_SUBT_ASSIGN( T[pElement->Col], T[I], *pElement );
+        {   T[pElement->Col] -= T[I] * pElement->Real;
             pElement = pElement->NextInRow;
         }
-        if (CMPLX_1_NORM(T[I]) > SLACK)
-        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
-            for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
+        if (ABS(T[I]) > SLACK)
+        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) );
+            for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
             ASy *= ScaleFactor;
         }
     }
 
 /* Compute 1-norm of T, which now contains v, and scale ||T|| to 1/SLACK. */
-    for (ASv = 0.0, I = Size; I > 0; I--) ASv += CMPLX_1_NORM(T[I]);
+    for (ASv = 0.0, I = Size; I > 0; I--) ASv += ABS(T[I]);
     ScaleFactor = 1.0 / (SLACK * ASv);
     if (ScaleFactor < 0.5)
-    {   for (I = Size; I > 0; I--) SCLR_MULT_ASSIGN( T[I], ScaleFactor );
+    {   for (I = Size; I > 0; I--) T[I] *= ScaleFactor;
         ASy *= ScaleFactor;
     }
 
@@ -1833,29 +1833,33 @@
     {   pPivot = Matrix->Diag[I];
         pElement = pPivot->NextInCol;
         while (pElement != NULL)
-        {   /* Cmplx expr: T[I] -= T[pElement->Row] * *pElement. */
-            CMPLX_MULT_SUBT_ASSIGN( T[I], T[pElement->Row], *pElement );
+        {   T[I] -= pElement->Real * T[pElement->Row];
             pElement = pElement->NextInCol;
         }
-        CMPLX_MULT_ASSIGN( T[I], *pPivot );
-        if (CMPLX_1_NORM(T[I]) > SLACK)
-        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), CMPLX_1_NORM(T[I]) );
-            for (K = Size; K > 0; K--) SCLR_MULT_ASSIGN( T[K], ScaleFactor );
+        T[I] *= pPivot->Real;
+        if (ABS(T[I]) > SLACK)
+        {   ScaleFactor = 1.0 / MAX( SQR( SLACK ), ABS(T[I]) );
+            for (K = Size; K > 0; K--) T[K] *= ScaleFactor;
             ASy *= ScaleFactor;
         }
     }
 
 /* Compute 1-norm of T, which now contains z. */
-    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 += ABS(T[I]);
 
+#if NOT spCOMPLEX
     FREE( Tm );
+#endif
 
     Linpack = ASy / ASz;
     OLeary = E / MaxY;
     InvNormOfInverse = MIN( Linpack, OLeary );
     return InvNormOfInverse / NormOfMatrix;
+#endif /* REAL */
 }
-#endif /* spCOMPLEX */
+
+
+

Places

File fasthenry-3.0.patch of Package fasthenry (Revision 1fde13fe4b3f7976b0aba7cc6ec2a760)

Places