Index: /trunk/psLib/src/math/psSparse.c
===================================================================
--- /trunk/psLib/src/math/psSparse.c	(revision 9990)
+++ /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 */
Index: /trunk/psLib/src/math/psSparse.h
===================================================================
--- /trunk/psLib/src/math/psSparse.h	(revision 9990)
+++ /trunk/psLib/src/math/psSparse.h	(revision 9991)
@@ -65,3 +65,34 @@
                        );
 
+# if (0)
+    // The border elements of a sparse matrix equation:
+    // 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)
+    typedef struct
+    {
+        psImage *Bij;                      // Aij contains the populated elements of the matrix
+        psImage *Qii;                      // Bfj contains the elements of the vector Bf
+        int Nelem;                         // Number of elements (long dimension of Bij, 0-j)
+        int Nrows;                         // Number of rows (size of Qii)
+    }
+psSparseBorder;
+
+
+// allocate a sparse matrix structure
+psSparse *psSparseBorderAlloc(int Nrows, int Nelem);
+
+// add a new matrix element
+// user should only add elements above the diagonal
+bool psSparseBorderMatrixElement(psSparse *sparse, // Matrix to which to add
+                                 int i, int j, // Matrix indices at which to add
+                                 float value  // Value to add
+                                );
+
+// define a new sparse matrix equation vector element
+void psSparseBorderVectorElement(psSparse *sparse, // Matrix to which to add
+                                 int i, int j, // Index to add
+                                 float value // Value to add
+                                );
+# endif /* gene's dev work */
+
 #endif /* PS_SPARSE_H */