Changeset 6891 for trunk/pois/src/poisSolveEquation.c

Timestamp:

Apr 18, 2006, 12:42:12 PM (20 years ago)

Author:

Paul Price

Message:

Updating to psLib rel11. Most changes are setting the 'n' element of the vector or array upon allocation. Needed to change the way the image subsets are done: not sure why we switched back to the old method of trim regions, but it seems we did.

File:

• trunk/pois/src/poisSolveEquation.c (modified) (4 diffs)

trunk/pois/src/poisSolveEquation.c

@@ -4 +4 @@
 #include "pois.h"
 
-psVector *poisSolveEquation(psVector *solution, // Solution vector, or NULL
-                            const psArray *stamps, // Array of stamps
-                            const poisConfig *config
+psVector *poisSolveEquation(psVector *solution, // Solution vector, or NULL
+                            const psArray *stamps, // Array of stamps
+                            const poisConfig *config
     )
 {
@@ -14 +14 @@
 
     int size = (2 * config->xKernel + 1) * (2 * config->yKernel + 1) * (config->spatialOrder + 1) *
-        (config->spatialOrder + 2) / 2 + 1; // Size of matrix and vector
+        (config->spatialOrder + 2) / 2 + 1; // Size of matrix and vector
 
     assert(!solution || solution->n == size);
@@ -20 +20 @@
     psImage *matrix = psImageAlloc(size, size, PS_TYPE_F64);
     psVector *vector = psVectorAlloc(size, PS_TYPE_F64);
+    vector->n = size;
     for (int j = 0; j < size; j++) {
-        for (int i = 0; i < size; i++) {
-            matrix->data.F64[j][i] = 0.0;
-        }
-        vector->data.F64[j] = 0.0;
+        for (int i = 0; i < size; i++) {
+            matrix->data.F64[j][i] = 0.0;
+        }
+        vector->data.F64[j] = 0.0;
     }
 
     for (int s = 0; s < stamps->n; s++) {
-        poisStamp *stamp = stamps->data[s]; // Stamp of interest
-        psImage *stampMatrix = stamp->matrix; // Corresponding matrix
-        psVector *stampVector = stamp->vector; // Corresponding vector
+        poisStamp *stamp = stamps->data[s]; // Stamp of interest
+        psImage *stampMatrix = stamp->matrix; // Corresponding matrix
+        psVector *stampVector = stamp->vector; // Corresponding vector
 
-        // Check inputs
-        assert(matrix->type.type == PS_TYPE_F64);
-        assert(vector->type.type == PS_TYPE_F64);
-        assert(matrix->numCols == size && matrix->numRows == size);
-        assert(vector->n == size);
+        // Check inputs
+        assert(matrix->type.type == PS_TYPE_F64);
+        assert(vector->type.type == PS_TYPE_F64);
+        assert(matrix->numCols == size && matrix->numRows == size);
+        assert(vector->n == size);
 
-        if (stamp->status == POIS_STAMP_USED) {
-            (void)psBinaryOp(matrix, matrix, "+", stampMatrix);
-            (void)psBinaryOp(vector, vector, "+", stampVector);
-        }
+        if (stamp->status == POIS_STAMP_USED) {
+            (void)psBinaryOp(matrix, matrix, "+", stampMatrix);
+            (void)psBinaryOp(vector, vector, "+", stampVector);
+        }
     }
 
@@ -53 +54 @@
     printf("Matrix:\n");
     for (int i = 0; i < matrix->numCols; i++) {
-        for (int j = 0; j < matrix->numRows; j++) {
-            printf("%f ", matrix->data.F64[i][j]);
-        }
-        printf("\n");
+        for (int j = 0; j < matrix->numRows; j++) {
+            printf("%f ", matrix->data.F64[i][j]);
+        }
+        printf("\n");
     }
 #endif