Index: trunk/psLib/src/math/psSparse.c
===================================================================
--- trunk/psLib/src/math/psSparse.c	(revision 9730)
+++ trunk/psLib/src/math/psSparse.c	(revision 9991)
@@ -56,9 +56,9 @@
     PS_ASSERT_INT_NONNEGATIVE(j, false);
 
+    // EAM : this is now a fatal error
     if (i < j) {
-        psLogMsg(__func__, PS_LOG_WARN, "i=%d, j=%d refers to a sub-diagonal element; values switched.\n", i, j);
-        int temp = i;
-        i = j;
-        j = temp;
+        // psError(PS_ERR_UNKNOWN, true, "i=%d, j=%d refers to a sub-diagonal element. not allowed!");
+        psAbort (__func__, "i=%d, j=%d refers to a sub-diagonal element. not allowed!");
+        return false;
     }
 
@@ -209,2 +209,65 @@
     return true;
 }
+
+# if (0)
+
+    /*** psSparseBorder Functions : these are used to solve a matrix equation of the form:
+      A x = f  where A is partitioned into:
+      A = |S B| where Q is a low-rank square matrix (N<20) 
+          |B Q| and B is a rectangular band (technically this is B and B^T)
+      and S is a sparse matrix.  
+    */
+
+    static void psSparseBorderFree(psSparse *sparse)
+{
+    if (!sparse) {
+        return;
+    }
+    psFree(sparse->Bij);
+    psFree(sparse->Qii);
+    return;
+}
+
+// allocate a sparse matrix container for Nrows, with Nelem slots allocated
+psSparseBorder *psSparseBorderAlloc (int Nrows, int Nelem)
+{
+    psSparseBorder *sparse = (psSparseBorder *)psAlloc(sizeof(psSparseBorder));
+    psMemSetDeallocator(sparse, (psFreeFunc) psSparseBorderFree);
+
+    sparse->Nrows = Nrows;
+    sparse->Nelem = Nelem;
+
+    sparse->Bij = psImageAlloc(Nelem, Nrows, PS_DATA_F32);
+    psInitImage (sparse->Bij, 0.0);
+
+    sparse->Qii = psImageAlloc(Nrows, Nrows, PS_DATA_F32);
+    psInitImage (sparse->Qii, 0.0);
+
+    return sparse;
+}
+
+// add elements to sparse->Qii
+bool psSparseBorderMatrixElement(psSparseBorder *sparse, int i, int j, float value)
+{
+    PS_ASSERT_PTR_NON_NULL(sparse, false);
+    PS_ASSERT_PTR_NON_NULL(sparse->Qii, false);
+    PS_ASSERT_INT_NONNEGATIVE(i, false);
+    PS_ASSERT_INT_NONNEGATIVE(j, false);
+
+    // check i,j against sparse->Qii->nX,nY
+    sparse->Qii->data.F32[i][j] = value;
+    return true;
+}
+
+// add elements to sparse->Bij
+bool psSparseBorderVectorElement(psSparseBorder *sparse, int i, int j, float value)
+{
+    PS_ASSERT_PTR_NON_NULL(sparse, false);
+    PS_ASSERT_PTR_NON_NULL(sparse->Qii, false);
+    PS_ASSERT_INT_NONNEGATIVE(i, false);
+    PS_ASSERT_INT_NONNEGATIVE(j, false);
+
+    sparse->Bij->data.F32[i][j] = value;
+    return true;
+}
+# endif /* gene's dev work */