Index: trunk/psLib/src/math/psMinimize.c
===================================================================
--- trunk/psLib/src/math/psMinimize.c	(revision 1798)
+++ trunk/psLib/src/math/psMinimize.c	(revision 1819)
@@ -9,6 +9,6 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-11 03:03:36 $
+ *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-16 19:30:24 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -865,13 +865,15 @@
                   psMinimizeLMFunc func)
 {
-    psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32);
-    psVector *perm = psVectorAlloc(params->n, PS_TYPE_F32);
-    psVector *newParams = psVectorAlloc(params->n, PS_TYPE_F32);
+    psVector *beta = psVectorAlloc(params->n, PS_TYPE_F64);
+    psVector *perm = psVectorAlloc(params->n, PS_TYPE_F64);
+    psVector *newParamsF64 = psVectorAlloc(params->n, PS_TYPE_F64);
+    psVector *newParamsF32 = psVectorAlloc(params->n, PS_TYPE_F32);
     psVector *origParams = psVectorAlloc(params->n, PS_TYPE_F32);
     psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32);
-    psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F32);
-    psImage *aOut = psImageAlloc(params->n, params->n, PS_TYPE_F32);
+    psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F64);
+    psImage *aOut = psImageAlloc(params->n, params->n, PS_TYPE_F64);
     psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
     psVector *newDeriv = psVectorAlloc(params->n, PS_TYPE_F32);
+    psVector *NRparams = psVectorAlloc(params->n, PS_TYPE_F32);
     int i;
     int j;
@@ -879,5 +881,12 @@
     float newValue;
     float oldValue;
-    float lamda = 1.0;
+    float lamda = 0.00005;
+
+    psTrace(".psLib.dataManip.psMinimizeLM", 4,
+            "---- psMinimizeLM() begin ----\n");
+    psTrace(".psLib.dataManip.psMinimizeLM", 6,
+            "min->maxIter is %d\n", min->maxIter);
+    psTrace(".psLib.dataManip.psMinimizeLM", 6,
+            "min->tol is %f\n", min->tol);
 
     for (i=0;i<params->n;i++) {
@@ -885,50 +894,97 @@
     }
 
-    min->lastDelta = -HUGE;
+    min->lastDelta = HUGE;
+    min->lastDelta = 12345.0;
     min->iter = 0;
-    min->value = func(deriv, params, coords);
-    for (i=0;i<params->n;i++) {
-        if ((paramMask != NULL) &&
-                (paramMask->data.U8[i] != 0)) {
-            deriv->data.F32[i] = 0.0;
-        }
-    }
 
     while ((min->lastDelta > min->tol) &&
             (min->iter < min->maxIter)) {
+        psTrace(".psLib.dataManip.psMinimizeLM", 4,
+                "------------------------------------------------------\n");
+        psTrace(".psLib.dataManip.psMinimizeLM", 4,
+                "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
+
+        min->value = func(deriv, params, coords);
 
         for (i=0;i<params->n;i++) {
-            beta->data.F32[i] = deriv->data.F32[i];
-        }
-
+            psTrace(".psLib.dataManip.psMinimizeLM", 6,
+                    "params->data.F32[%d] is %f.\n", i, params->data.F32[i]);
+        }
+        for (i=0;i<params->n;i++) {
+            psTrace(".psLib.dataManip.psMinimizeLM", 6,
+                    "deriv->data.F32[%d] is %f.\n", i, deriv->data.F32[i]);
+        }
+        psTrace(".psLib.dataManip.psMinimizeLM", 6,
+                "min->value is (%f)\n", min->value);
+
+        for (i=0;i<params->n;i++) {
+            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
+                deriv->data.F32[i] = 0.0;
+            }
+        }
+
+        // Calculate the BETA vector.
+        for (i=0;i<params->n;i++) {
+            beta->data.F64[i] = (float) deriv->data.F32[i];
+            psTrace(".psLib.dataManip.psMinimizeLM", 6,
+                    "beta->data.F64[%d] is %f.\n", i, beta->data.F64[i]);
+        }
+
+        // Calculate the ALPHA matrix.
         for (j=0;j<params->n;j++) {
             for (k=0;k<params->n;k++) {
-                alpha->data.F32[j][k] = deriv->data.F32[j] * deriv->data.F32[j];
-            }
-        }
+                alpha->data.F32[j][k] = deriv->data.F32[k] *
+                                        deriv->data.F32[j];
+            }
+        }
+
+        // Calculate the matrix A.
         for (j=0;j<params->n;j++) {
             for (k=0;k<params->n;k++) {
-                if (j != k) {
-                    A->data.F32[j][k] = alpha->data.F32[j][k];
+                if (j == k) {
+                    A->data.F64[j][k] =
+                        (double) ((1.0 + lamda) * alpha->data.F32[j][k]);
                 } else {
-                    A->data.F32[j][k] = (1.0 + lamda) * alpha->data.F32[j][k];
+                    A->data.F64[j][k] = (double) alpha->data.F32[j][k];
                 }
             }
         }
+        for (j=0;j<params->n;j++) {
+            psTrace(".psLib.dataManip.psMinimizeLM", 6, "Matrix A[][]:\n");
+            for (k=0;k<params->n;k++) {
+                psTrace(".psLib.dataManip.psMinimizeLM", 6, "%f ", A->data.F64[j][k]);
+            }
+            psTrace(".psLib.dataManip.psMinimizeLM", 6, "Matrix A[][]:\n");
+        }
+
+        // Solve A * alpha = Beta
         aOut = psMatrixLUD(aOut, perm, A);
-        newParams = psMatrixLUSolve(newParams, aOut, beta, perm);
-
+        newParamsF64 = psMatrixLUSolve(newParamsF64, aOut, beta, perm);
+
+        // Mask any masked parameters.
         for (i=0;i<params->n;i++) {
-            if ((paramMask != NULL) &&
-                    (paramMask->data.U8[i] != 0)) {
-                newParams->data.F32[i] = origParams->data.F32[i];
-            }
+            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
+                newParamsF64->data.F64[i] = (double) origParams->data.F32[i];
+            }
+            newParamsF32->data.F32[i] = (float) newParamsF64->data.F64[i];
+
+            psTrace(".psLib.dataManip.psMinimizeLM", 6,
+                    "newParamsF32->data.F32[%d] is %f.\n", i, newParamsF32->data.F32[i]);
+            NRparams->data.F32[i] = params->data.F32[i] - newParamsF32->data.F32[i];
+        }
+
+        psTrace(".psLib.dataManip.psMinimizeLM", 6,
+                "Calling func() with new parameters:\n");
+        for (i=0;i<params->n;i++) {
+            psTrace(".psLib.dataManip.psMinimizeLM", 6,
+                    "NRparams->data.F32[%d] is %f.\n", i, NRparams->data.F32[i]);
         }
 
         oldValue = min->value;
-        newValue = func(newDeriv, newParams, coords);
+        newValue = func(newDeriv, NRparams, coords);
+        psTrace(".psLib.dataManip.psMinimizeLM", 4,
+                "old/new values are (%f, %f)\n", oldValue, newValue);
         for (i=0;i<params->n;i++) {
-            if ((paramMask != NULL) &&
-                    (paramMask->data.U8[i] != 0)) {
+            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
                 deriv->data.F32[i] = 0.0;
             }
@@ -942,6 +998,6 @@
             // before the last function evaluation.
             for (i=0;i<params->n;i++) {
-                params->data.F32[i] = newParams->data.F32[i];
-                deriv->data.F32[i] = newDeriv->data.F32[i];
+                //                params->data.F32[i] = (float) newParamsF64->data.F64[i];
+                params->data.F32[i] = (float) NRparams->data.F32[i];
             }
             min->value = func(deriv, params, coords);
@@ -951,9 +1007,12 @@
             lamda*= 10.0;
         }
+        psTrace(".psLib.dataManip.psMinimizeLM", 4,
+                "lamda is %f\n", lamda);
         min->iter++;
     }
     psFree(beta);
     psFree(perm);
-    psFree(newParams);
+    psFree(newParamsF64);
+    psFree(newParamsF32);
     psFree(origParams);
     psFree(alpha);
@@ -968,4 +1027,6 @@
     }
 
+    psTrace(".psLib.dataManip.psMinimizeLM", 4,
+            "---- psMinimizeLM() end (false) ----\n");
     return(false);
 }
@@ -1299,5 +1360,5 @@
         psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
                 "psDetermineBracket(): iteration %d\n", iter);
-        if ((b < a) && (b < c)) {
+        if ((fb < fa) && (fb < fc)) {
             bracket->data.F32[0] = a;
             bracket->data.F32[1] = b;
@@ -1738,5 +1799,7 @@
 
                 psTrace(".psLib.dataManip.psMinimizePowell", 6,
-                        "WOW (%f %f %f)\n", u->data.F32[i], Q->data.F32[i], params->data.F32[i]);
+                        "u=Q-P (%f = %f - %f)\n", u->data.F32[i],
+                        Q->data.F32[i],
+                        params->data.F32[i]);
 
             } else {
@@ -1747,6 +1810,8 @@
         // 5: Move Q only in the direction u, and minimize the function.
         psTrace(".psLib.dataManip.psMinimizePowell", 6, "HERE 01\n");
-        psTrace(".psLib.dataManip.psMinimizePowell", 6,
-                "u[0] is %f\n", u->data.F32[0]);
+        for (i=0;i<numDims;i++) {
+            psTrace(".psLib.dataManip.psMinimizePowell", 6,
+                    "u[i] is %f\n", u->data.F32[i]);
+        }
 
         mul = psLineMin(min, params, u, paramMask, coords, func);