Index: /trunk/psLib/src/math/psSparse.c
===================================================================
--- /trunk/psLib/src/math/psSparse.c	(revision 10000)
+++ /trunk/psLib/src/math/psSparse.c	(revision 10001)
@@ -229,44 +229,144 @@
 }
 
-// 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);
+// allocate a sparse matrix border container for a system with Nrows + Nborder elements
+// the supplied psSparse has a size of Nrows
+psSparseBorder *psSparseBorderAlloc(psSparse *sparse, int Nborder)
+{
+    psSparseBorder *border = (psSparseBorder *)psAlloc(sizeof(psSparseBorder));
+    psMemSetDeallocator(border, (psFreeFunc) psSparseBorderFree);
+
+    int Nrows = sparse->Nrows;
+    border->Nrows = Nrows;
+    border->Nborder = Nborder;
+
+    border->Bij = psImageAlloc(Nrows, Nborder, PS_DATA_F32);
+    psImageInit (border->Bij, 0.0);
+
+    border->Qii = psImageAlloc(Nborder, Nborder, PS_DATA_F32);
+    psImageInit (border->Qii, 0.0);
+
+    border->Yi = psVectorAlloc(Nborder, PS_DATA_F32);
+    psVectorInit (border->Yi, 0.0);
+
+    border->Gi = psVectorAlloc(Nborder, PS_DATA_F32);
+    psVectorInit (border->Gi, 0.0);
+
+    return border;
+}
+
+// add elements to border->Qii
+bool psSparseBorderElementQ(psSparseBorder *border, int i, int j, float value)
+{
+    PS_ASSERT_PTR_NON_NULL(border, false);
+    PS_ASSERT_PTR_NON_NULL(border->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);
+    // check i,j against border->Qii->nX,nY
+    border->Qii->data.F32[i][j] = value;
+    return true;
+}
+
+// add elements to border->Bij
+bool psSparseBorderElementB(psSparseBorder *border, int i, int j, float value)
+{
+    PS_ASSERT_PTR_NON_NULL(border, false);
+    PS_ASSERT_PTR_NON_NULL(border->Bij, false);
     PS_ASSERT_INT_NONNEGATIVE(i, false);
     PS_ASSERT_INT_NONNEGATIVE(j, false);
 
-    sparse->Bij->data.F32[i][j] = value;
+    border->Bij->data.F32[i][j] = value;
+    return true;
+}
+
+// add elements to border->Yi
+bool psSparseBorderElementY(psSparseBorder *border, int i, float value)
+{
+    PS_ASSERT_PTR_NON_NULL(border, false);
+    PS_ASSERT_PTR_NON_NULL(border->Yi, false);
+    PS_ASSERT_INT_NONNEGATIVE(i, false);
+
+    border->Yi->data.F32[i] = value;
+    return true;
+}
+
+// add elements to border->Gi
+bool psSparseBorderElementY(psSparseBorder *border, int i, float value)
+{
+    PS_ASSERT_PTR_NON_NULL(border, false);
+    PS_ASSERT_PTR_NON_NULL(border->Gi, false);
+    PS_ASSERT_INT_NONNEGATIVE(i, false);
+
+    border->Gi->data.F32[i] = value;
+    return true;
+}
+
+// perform the operation dG = B*x
+psVector *psSparseBorderUpperProduct (psSparse *border)
+{
+
+    int Nborder = border->Nborder;
+    int Nrow = border->Nrows;
+
+    psVector *dG = psVectorAlloc (Nborder, PS_DATA_F32);
+    psVectorInit (dG, 0.0);
+
+    for (int j = 0; j < Nborder; j++) {
+        double value = 0;
+        for (int i = 0; i < Nrows; i++) {
+            value += border->Bij->data.F32[i][j] * x->data.F32[i];
+        }
+        dG->data.F32[j] = value;
+    }
+
+    return dG;
+}
+
+// perform the operation dF = B^T*y
+psVector *psSparseBorderLowerProduct (psSparse *border)
+{
+
+    int Nborder = border->Nborder;
+    int Nrow = border->Nrows;
+
+    psVector *dF = psVectorAlloc (Nrows, PS_DATA_F32);
+    psVectorInit (dF, 0.0);
+
+    for (int i = 0; i < Nrows; i++) {
+        double value = 0;
+        for (int j = 0; j < Nborder; j++) {
+            value += border->Bij->data.F32[i][j] * border->Yi->data.F32[j];
+        }
+        dF->data.F32[i] = value;
+    }
+
+    return dF;
+}
+
+// perform the operation dF = B^T*y
+psVector *psSparseBorderUpperDelta (psSparse *border, psVector *dF)
+{
+
+    int Nborder = border->Nborder;
+    int Nrow = border->Nrows;
+
+    for (int i = 0; i < Nrows; i++) {
+        border->sparse->Bfj->data.F32[i] -= dF->data.F32[i];
+    }
+
+    return true;
+}
+
+// perform the operation dF = B^T*y
+psVector *psSparseBorderLowerDelta (psSparse *border, psVector *dG)
+{
+
+    int Nborder = border->Nborder;
+    int Nrow = border->Nrows;
+
+    for (int i = 0; i < Nborder; i++) {
+        border->Gi->data.F32[i] -= dG->data.F32[i];
+    }
+
     return true;
 }
Index: /trunk/psLib/src/math/psSparse.h
===================================================================
--- /trunk/psLib/src/math/psSparse.h	(revision 10000)
+++ /trunk/psLib/src/math/psSparse.h	(revision 10001)
@@ -71,8 +71,11 @@
     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)
+        psSparse *sparse;  // corresponding sparse matrix equation
+        psImage *Bij;   // Aij contains the populated elements of the matrix
+        psImage *Qii;   // Bfj contains the elements of the vector Bf
+        psVector *Yi;   // lower independent var
+        psVector *Gi;   // lower dependent var
+        int Nrows;   // Number of rows (long dimension of Bij, 0-j)
+        int Nborder;   // Number of border elements (size of Qii)
     }
 psSparseBorder;
@@ -80,5 +83,5 @@
 
 // allocate a sparse matrix structure
-psSparse *psSparseBorderAlloc(int Nrows, int Nelem);
+psSparseBorder *psSparseBorderAlloc(psSparse *sparse, int Nborder);
 
 // add a new matrix element