Index: /trunk/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtractionEquation.c	(revision 24619)
+++ /trunk/psModules/src/imcombine/pmSubtractionEquation.c	(revision 24620)
@@ -15,5 +15,5 @@
 #include "pmSubtractionVisual.h"
 
-//#define TESTING                         // TESTING output for debugging; may not work with threads!
+// #define TESTING                         // TESTING output for debugging; may not work with threads!
 
 #define USE_VARIANCE                    // Include variance in equation?
@@ -789,5 +789,23 @@
         for (int i = 0; i < stamps->num; i++) {
             pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+
             if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
+
+#ifdef TESTING
+	      // XXX double-check for NAN in data:
+		for (int iy = 0; iy < stamp->matrix1->numRows; iy++) {
+		    for (int ix = 0; ix < stamp->matrix1->numCols; ix++) {
+			if (!isfinite(stamp->matrix1->data.F64[iy][ix])) {
+			    fprintf (stderr, "WARNING: NAN in matrix1\n");
+			}
+		    }
+		}
+		for (int ix = 0; ix < stamp->vector1->n; ix++) {
+		    if (!isfinite(stamp->vector1->data.F64[ix])) {
+			fprintf (stderr, "WARNING: NAN in vector1\n");
+		    }
+		}
+#endif
+
                 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1);
                 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1);
@@ -798,7 +816,21 @@
             }
         }
+
+#ifdef TESTING
+	for (int ix = 0; ix < sumVector->n; ix++) {
+	    if (!isfinite(sumVector->data.F64[ix])) {
+		fprintf (stderr, "WARNING: NAN in vector1\n");
+	    }
+	}
+#endif
+
         calculatePenalty(sumVector, kernels);
 
 #ifdef TESTING
+	for (int ix = 0; ix < sumVector->n; ix++) {
+	    if (!isfinite(sumVector->data.F64[ix])) {
+		fprintf (stderr, "WARNING: NAN in vector1\n");
+	    }
+	}
         {
             psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
@@ -832,4 +864,12 @@
         }
         kernels->solution1 = psMatrixLUSolution(kernels->solution1, luMatrix, sumVector, permutation);
+
+	// XXX double-check for NAN in data:
+	for (int ix = 0; ix < kernels->solution1->n; ix++) {
+	    if (!isfinite(kernels->solution1->data.F64[ix])) {
+		fprintf (stderr, "WARNING: NAN in vector1\n");
+	    }
+	}
+
 
         psFree(sumVector);