Changeset 10017 for trunk/psLib/src/math/psSparse.c

Timestamp:

Nov 16, 2006, 8:09:44 AM (19 years ago)

Author:

magnier

Message:

finished psSparseBorder, runs tests ok

File:

• trunk/psLib/src/math/psSparse.c (modified) (7 diffs)

trunk/psLib/src/math/psSparse.c

@@ -266 +266 @@
 
     // check i,j against border->Tjj->nX,nY
-    border->Tjj->data.F32[i][j] = value;
+    border->Tjj->data.F32[j][i] = value;
     return true;
 }
@@ -278 +278 @@
     PS_ASSERT_INT_NONNEGATIVE(j, false);
 
-    border->Bij->data.F32[i][j] = value;
+    border->Bij->data.F32[j][i] = value;
     return true;
 }
@@ -307 +307 @@
         double value = 0;
         for (int i = 0; i < Nrows; i++) {
-            value += border->Bij->data.F32[i][j] * xVec->data.F32[i];
+            value += border->Bij->data.F32[j][i] * xVec->data.F32[i];
         }
         dG->data.F32[j] = value;
@@ -329 +329 @@
         double value = 0;
         for (int j = 0; j < Nborder; j++) {
-            value += border->Bij->data.F32[i][j] * yVec->data.F32[j];
+            value += border->Bij->data.F32[j][i] * yVec->data.F32[j];
         }
         dF->data.F32[i] = value;
@@ -350 +350 @@
         double value = 0;
         for (int j = 0; j < Nborder; j++) {
-            value += border->Tjj->data.F32[j][i] * yVec->data.F32[j];
+            value += border->Tjj->data.F32[i][j] * yVec->data.F32[j];
         }
         dG->data.F32[i] = value;
@@ -372 +372 @@
 
 // perform the operation dF = B^T*y
-bool psSparseBorderLowerDelta (psSparseBorder *border, psVector *dG)
-{
-
+psVector *psSparseBorderLowerDelta (psVector *Go, psSparseBorder *border, psVector *dG)
+{
     int Nborder = border->Nborder;
 
+    Go = psVectorRecycle (Go, Nborder, PS_TYPE_F64);
+    psVectorInit (Go, 0.0);
+
     for (int i = 0; i < Nborder; i++) {
-        border->Gj->data.F32[i] -= dG->data.F32[i];
-    }
-
-    return true;
+        Go->data.F64[i] = border->Gj->data.F32[i] - dG->data.F32[i];
+    }
+    return Go;
 }
 
@@ -417 +418 @@
     psVector *yVec = *yFit;
 
+    // save the original value of Bfj, alloc other temp vectors
+    // XXX be careful about TYPE requirements of support functions...
+    psVector *Fo = psVectorCopy (NULL, border->sparse->Bfj, PS_TYPE_F32);
+    psVector *dF = psVectorAlloc (border->Nrows, PS_TYPE_F32);
+    psVector *Go = psVectorAlloc (border->Nborder, PS_TYPE_F64);
+    psVector *dG = psVectorAlloc (border->Nborder, PS_TYPE_F32);
+    psImage *square = psImageAlloc (border->Nborder, border->Nborder, PS_TYPE_F64);
+
     for (int j = 0; j < Niter; j++) {
 
         // solve Sx = f
-        xVec = psSparseSolve(xVec, constraint, border->sparse, Niter);
-
-        psVector *dG = psSparseBorderLowerProduct (NULL, border, xVec);
-        psSparseBorderLowerDelta (border, dG);
-
-        // XXX not defined
-        psImage *square = psImageCopy (NULL, border->Tjj, PS_TYPE_F32);
-        psVectorCopy (yVec, border->Gj, PS_TYPE_F32);
-        psMatrixGJSolve (square, yVec);
-        psFree (square);
-
-        psVector *dF = psSparseBorderUpperProduct (NULL, border, yVec);
+        xVec = psSparseSolve(xVec, constraint, border->sparse, 4);
+
+        dG = psSparseBorderLowerProduct (dG, border, xVec);
+
+        // XXX i have an error here: i'm not
+        // I need to reset to the original value of g before subtracting
+        // XXX this function returns an psVector:F64 so we can us it in GJ
+        Go = psSparseBorderLowerDelta (Go, border, dG);
+
+        // use gauss-jordan to solve the lower square
+        // square is modified by GJSolve, so re-copy each time
+        square = psImageCopy (square, border->Tjj, PS_TYPE_F64);
+        psMatrixGJSolve (square, Go);
+        yVec = psVectorCopy (yVec, Go, PS_TYPE_F32);
+
+        // calculate the delta relative to the original Bfj:
+        border->sparse->Bfj = psVectorCopy (border->sparse->Bfj, Fo, PS_TYPE_F32);
+        dF = psSparseBorderUpperProduct (dF, border, yVec);
         psSparseBorderUpperDelta (border, dF);
-
-        psFree (dF);
-        psFree (dG);
-    }
+    }
+
+    psFree (dG);
+    psFree (Go);
+    psFree (dF);
+    psFree (Fo);
+    psFree (square);
 
     *xFit = xVec;