Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14709)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14710)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.50 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-08-30 21:45:26 $
+ *  @version $Revision: 1.51 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-30 21:52:56 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -621,9 +621,10 @@
 
         // Generate least-squares vector and matrix
+        bool bad = false;           // Are there bad values?
         for (int j = 0; j < numKernels; j++) {
             psKernel *jConv = convolutions->data[j]; // Convolution for j-th element
 
-            // Generate upper diagonals of matrix
-            for (int k = 0; k < numKernels; k++) {
+            // Generate matrix
+            for (int k = j; k < numKernels; k++) {
                 psKernel *kConv = convolutions->data[k]; // Convolution for k-th element
                 double sumCC = 0.0; // Sum of the convolution products
@@ -633,14 +634,19 @@
                     }
                 }
+
+                if (!isfinite(sumCC) {
+                    bad = true;
+                }
                 // Generate the pseudo-convolutions from the spatial polynomial terms
                 for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) {
                     for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder;
                          jxOrder++, jIndex += numKernels) {
-                        for (int kyOrder = 0, kIndex = k; kyOrder <= spatialOrder; kyOrder++) {
-                            for (int kxOrder = 0; kxOrder <= spatialOrder - kyOrder;
+                        for (int kyOrder = jyOrder, kIndex = k; kyOrder <= spatialOrder; kyOrder++) {
+                            for (int kxOrder = jxOrder; kxOrder <= spatialOrder - kyOrder;
                                  kxOrder++, kIndex += numKernels) {
-                                stampMatrix->data.F64[jIndex][kIndex] = sumCC *
-                                    polyValues->data.F64[jyOrder][jxOrder] *
-                                    polyValues->data.F64[kyOrder][kxOrder];
+                                double convPoly = sumCC * polyValues->data.F64[jyOrder][jxOrder] *
+                                    polyValues->data.F64[kyOrder][kxOrder]; // Temporary value
+                                stampMatrix->data.F64[jIndex][kIndex] = convPoly;
+                                stampMatrix->data.F64[kIndex][jIndex] = convPoly;
                             }
                         }
@@ -659,10 +665,15 @@
                 }
             }
+            if (!isfinite(sumC) || !isfinite(sumIC)) {
+                bad = true;
+            }
             for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) {
                 for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder;
                      jxOrder++, jIndex += numKernels) {
                     stampVector->data.F64[jIndex] = sumIC * polyValues->data.F64[jyOrder][jxOrder];
-                    stampMatrix->data.F64[jIndex][bgIndex] = sumC * polyValues->data.F64[jyOrder][jxOrder];
-                    stampMatrix->data.F64[bgIndex][jIndex] = sumC * polyValues->data.F64[jyOrder][jxOrder];
+
+                    double convPoly = sumC * polyValues->data.F64[jyOrder][jxOrder]; // Value
+                    stampMatrix->data.F64[jIndex][bgIndex] = convPoly;
+                    stampMatrix->data.F64[bgIndex][jIndex] = convPoly;
                 }
             }
@@ -681,36 +692,9 @@
             }
         }
+        if (!isfinite(sum1) || !isfinite(sumI)) {
+            bad = true;
+        }
         stampMatrix->data.F64[bgIndex][bgIndex] = sum1;
         stampVector->data.F64[bgIndex] = sumI;
-
-        // Fill in lower diagonals of symmetric matrix, while checking for bad values
-        // Note, there are two symmetries going on here --- the matrix is C_i C_j P_i P_j
-        bool bad = false;           // Are there bad values?
-#if 0
-        for (int j = 0; j < numKernels; j++) {
-            for (int jSpatial = 0, jIndex = j; jSpatial < numSpatial; jSpatial++, jIndex += numKernels) {
-                for (int k = 0; k < j; k++) {
-                    for (int kSpatial = 0, kIndex = k; kSpatial < numSpatial;
-                         kSpatial++, kIndex += numKernels) {
-                        double value = stampMatrix->data.F64[kIndex][jIndex]; // Value of matrix
-                        stampMatrix->data.F64[jIndex][kIndex] = value;
-                        if (!isfinite(value)) {
-                            bad = true;
-                        }
-                    }
-                }
-
-                stampMatrix->data.F64[bgIndex][jIndex] = stampMatrix->data.F64[jIndex][bgIndex];
-                if (!isfinite(stampMatrix->data.F64[jIndex][bgIndex]) ||
-                    !isfinite(stampMatrix->data.F64[jIndex][jIndex]) ||
-                    !isfinite(stampVector->data.F64[jIndex])) {
-                    bad = true;
-                }
-            }
-        }
-#endif
-        if (!isfinite(stampVector->data.F64[bgIndex])) {
-            bad = true;
-        }
 
         if (bad) {