Changeset 23984 for trunk/psLib/src/math/psMinimizeLMM.c

Timestamp:

Apr 28, 2009, 10:59:54 AM (17 years ago)

Author:

eugene

Message:

added some test code to verify GJ elimination: move to a real test

File:

• trunk/psLib/src/math/psMinimizeLMM.c (modified) (2 diffs)

trunk/psLib/src/math/psMinimizeLMM.c

@@ -80 +80 @@
     assert (alpha->numCols == alpha->numRows);
 
+# define TESTGJ 0
+# if (TESTGJ)
+    lambda = 0.0;
+#endif
+
     // set new guess values, applying (1+lambda) scaling to pivots
     Beta = psVectorCopy(Beta, beta, PS_TYPE_F32);
@@ -92 +97 @@
         return(false);
     }
+
+// XXX check that the GJ solver works:
+# if (TESTGJ) 
+    psImage *out = psImageAlloc (alpha->numRows, alpha->numCols, PS_TYPE_F32);
+    for (int oy = 0; oy < out->numRows; oy++) {
+        for (int ox = 0; ox < out->numCols; ox++) {
+            float value = 0;
+            for (int i = 0; i < alpha->numCols; i++) {
+                value += alpha->data.F32[i][ox]*Alpha->data.F32[oy][i];
+            }
+            out->data.F32[oy][ox] = value;
+        }
+    }
+
+    psVector *vect = psVectorAlloc (beta->n, PS_TYPE_F32);
+    for (int oy = 0; oy < vect->n; oy++) {
+        float value = 0;
+        for (int i = 0; i < alpha->numCols; i++) {
+            value += alpha->data.F32[oy][i]*Beta->data.F32[i];
+        }
+        vect->data.F32[oy] = value;
+    }
+    
+    psFree (out);
+    psFree (vect);
+
+# endif
 
     // check for non-finite entries in result