Index: /trunk/psLib/src/dataManip/Makefile
===================================================================
--- /trunk/psLib/src/dataManip/Makefile	(revision 1184)
+++ /trunk/psLib/src/dataManip/Makefile	(revision 1185)
@@ -13,6 +13,6 @@
            psMatrixVectorArithmetic.o \
            psFFT.o \
-           psImageIO.o
-#           psMinimize.o \
+           psImageIO.o \
+           psMinimize.o
 
 all: $(TARGET_STATIC)
Index: /trunk/psLib/src/dataManip/psMinimize.c
===================================================================
--- /trunk/psLib/src/dataManip/psMinimize.c	(revision 1184)
+++ /trunk/psLib/src/dataManip/psMinimize.c	(revision 1185)
@@ -7,4 +7,5 @@
 
 #include <gsl/gsl_multifit_nlin.h>
+#include <gsl/gsl_multimin.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
@@ -38,5 +39,4 @@
 }
 psModelData;
-
 // The first argument to evalModel() and
 // d_evalModel() specifies the data point.
@@ -47,17 +47,241 @@
 
 
+typedef struct
+{
+    int paramCount;   // Number of non-masked parameters.
+    psVector *restrict       initialGuess;
+    const psVector *restrict  coord;
+    const psVector *restrict paramMask;
+    float (*evalModel) (const psVector *, const psVector *);
+    float (*d_evalModel) (const psVector *, const psVector *, int);
+}
+psModelData2;
+
 
 /******************************************************************************
-    This routine must minimize an arbitrary function.
+p_psMinFunc(*v, *params): This routine calls the user supplied function to
+be minimized.  The "v" argument of this function corresponds to the
+parameters of the function to be minimized.
+ 
+ 
+ *****************************************************************************/
+double p_psMinFunc(const gsl_vector *v,
+                   void *params)
+{
+    int i;    // Loop index variable.
+    int j;    // Loop index variable.
+    float tmpf;    // Temporary floating point variable.
+    const psVector *restrict coord   = ((psModelData2 *) params)->coord;
+    const psVector *restrict mask     = ((psModelData2 *) params)->paramMask;
+    psVector *restrict initialGuess = ((psModelData2 *)params)->initialGuess;
+    float (*evalModel)(const psVector *, const psVector *) =
+        ((psModelData2 *) params)->evalModel;
+    psVector *inputParameterList = NULL;
+
+    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
+
+    if (mask != NULL) {
+        j = 0;
+        for (i=0;i<mask->n;i++) {
+            if (mask->data.U8[i] != 0) {
+                inputParameterList->data.F32[i] = initialGuess->data.F32[i];
+            } else {
+                inputParameterList->data.F32[i] = gsl_vector_get(v, j++);
+            }
+        }
+    } else {
+        for (i=0;i<initialGuess->n;i++) {
+            inputParameterList->data.F32[i] = gsl_vector_get(v, i);
+        }
+    }
+    //printf("HMMM: (%.1f %.1f %.1f %.1f)\n",
+    //       inputParameterList->data.F32[0],
+    //       inputParameterList->data.F32[1],
+    //       coord->data.F32[0],
+    //       coord->data.F32[1]);
+
+    tmpf = evalModel(inputParameterList, coord);
+    psFree(inputParameterList);
+    //printf("Called p_psMinFunc()\n");
+    return(tmpf);
+}
+
+void p_psMinFuncDeriv(const gsl_vector *v,
+                      void *params,
+                      gsl_vector *df)
+{
+    int i;    // Loop index variable.
+    int j;    // Loop index variable.
+    float tmpf;    // Temporary floating point variable.
+    const psVector *restrict coord   = ((psModelData2 *) params)->coord;
+    const psVector *restrict mask     = ((psModelData2 *) params)->paramMask;
+    psVector *restrict initialGuess = ((psModelData2 *)params)->initialGuess;
+    float (*d_evalModel)(const psVector *, const psVector *, int) =
+        ((psModelData2 *) params)->d_evalModel;
+    psVector *inputParameterList = NULL;
+
+    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
+
+    if (mask != NULL) {
+        j = 0;
+        for (i=0;i<mask->n;i++) {
+            if (mask->data.U8[i] != 0) {
+                inputParameterList->data.F32[i] = initialGuess->data.F32[i];
+            } else {
+                inputParameterList->data.F32[i] = gsl_vector_get(v, j++);
+            }
+        }
+    } else {
+        for (i=0;i<initialGuess->n;i++) {
+            inputParameterList->data.F32[i] = gsl_vector_get(v, i);
+        }
+    }
+
+    for (i=0;i<initialGuess->n;i++) {
+        tmpf = d_evalModel(inputParameterList, coord, i);
+        gsl_vector_set(df, i, tmpf);
+    }
+    psFree(inputParameterList);
+}
+
+/* Compute both f and df together. */
+void p_psMinFuncFuncDeriv(const gsl_vector *x,
+                          void *params,
+                          double *f,
+                          gsl_vector *df)
+{
+    *f = p_psMinFunc(x, params);
+    p_psMinFuncDeriv(x, params, df);
+}
+
+
+/******************************************************************************
+psMinimize(initialGuess, myFunction, myFunctionDeriv, coord, paramMask):
+ 
+This routine must minimize an arbitrary function; it must determine the set
+of parameters of that function such that the 
+ 
  *****************************************************************************/
 psVector *
-psMinimize(float (*myFunction)(const psVector *restrict),
-           psVector *restrict initialGuess,
-           psVector *restrict paramMask)
-{
-    return(NULL);
-}
-
-
+psMinimize(psVector *restrict initialGuess,
+           float (*myFunction)(const psVector *restrict, const psVector *restrict),
+           float (*myFunctionDeriv)(const psVector *restrict, const psVector *restrict, int),
+           const psVector *restrict coord,
+           const psVector *restrict paramMask)
+{
+    int status;
+    int i = 0;
+    int j = 0;
+    int iter = 0;
+    gsl_multimin_function_fdf f;
+    double *xInit = NULL;
+    const gsl_multimin_fdfminimizer_type *T;
+    gsl_multimin_fdfminimizer *s;
+    psModelData2 inputData;
+    gsl_vector *x;
+
+    inputData.initialGuess = initialGuess;
+    inputData.coord = coord;
+    inputData.paramMask = paramMask;
+    inputData.evalModel = myFunction;
+    inputData.d_evalModel = myFunctionDeriv;
+    inputData.paramCount = 0;
+
+    printf("psMinimize(): initialGuess->n is %d\n", initialGuess->n);
+    printf("psMinimize(): coord->n is %d\n", coord->n);
+    printf("psMinimize(): coord is (%.1f, %.1f)\n", coord->data.F32[0],
+           coord->data.F32[1]);
+
+    if (paramMask != NULL) {
+        for (i=0;i<paramMask->n;i++) {
+            if (paramMask->data.U8[i] != 0) {
+                inputData.paramCount++;
+            }
+        }
+    } else {
+        inputData.paramCount= initialGuess->n;
+    }
+
+    // The initial guess at the parameters for the function are written into
+    // the vector inputParameterList.  If the paramMask is not NULL, then those
+    // parameters are masked out.
+
+    xInit = (double *) psAlloc(inputData.paramCount * sizeof(double));
+    if (paramMask != NULL) {
+        j = 0;
+        for (i=0;i<initialGuess->n;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                xInit[j++] = initialGuess->data.F32[i];
+            }
+        }
+    } else {
+        for (i=0;i<initialGuess->n;i++) {
+            xInit[i] = initialGuess->data.F32[i];
+        }
+    }
+    f.f = &p_psMinFunc;
+    f.df = &p_psMinFuncDeriv;
+    f.fdf = &p_psMinFuncFuncDeriv;
+    f.n = inputData.paramCount;
+    // GUS: is this correct?
+    //f.params = xInit;
+    f.params = &inputData;
+
+
+    printf("inputData.paramCount is %d\n", inputData.paramCount);
+    x = gsl_vector_alloc(inputData.paramCount);
+    for (i=0;i<inputData.paramCount;i++) {
+        gsl_vector_set(x, i, xInit[i]);
+    }
+    T = gsl_multimin_fdfminimizer_conjugate_fr;
+    s = gsl_multimin_fdfminimizer_alloc(T, inputData.paramCount);
+
+    gsl_multimin_fdfminimizer_set(s, &f, x, 0.01, 1e-4);
+
+    do {
+        iter++;
+        status = gsl_multimin_fdfminimizer_iterate(s);
+
+        if (status)
+            break;
+
+        status = gsl_multimin_test_gradient(s->gradient, 1e-3);
+
+        if (status == GSL_SUCCESS)
+            printf ("Minimum found at:\n");
+
+    } while (status == GSL_CONTINUE && iter < 100);
+
+    // In the above steps we had removed the masked elements from the
+    // the solver.  This next code blocks puts those masked elements
+    // into the solution.
+    if (paramMask != NULL) {
+        j = 0;
+        for (i=0;i<initialGuess->n;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                initialGuess->data.F32[i] = gsl_vector_get(s->x, j++);
+            } else {
+                initialGuess->data.F32[i] = initialGuess->data.F32[i];
+            }
+        }
+    } else {
+        for (i=0;i<initialGuess->n;i++) {
+            initialGuess->data.F32[i] = gsl_vector_get(s->x, i);
+        }
+    }
+    psFree(xInit);
+    return(initialGuess);
+}
+
+
+
+
+
+// The first argument to evalModel() and
+// d_evalModel() specifies the data point.
+// It must have the same size as the second
+// dimension of *domain.
+// The second argument must have the same size
+// as *initialGuess and *paramMask.
 
 /******************************************************************************
Index: /trunk/psLib/src/dataManip/psMinimize.h
===================================================================
--- /trunk/psLib/src/dataManip/psMinimize.h	(revision 1184)
+++ /trunk/psLib/src/dataManip/psMinimize.h	(revision 1185)
@@ -14,8 +14,9 @@
 /** This routine must minimize a non-linear function */
 psVector *
-psMinimize(float (*myFunction)(const psVector *restrict), ///< Function to minimize
-           psVector *restrict initialGuess, ///< Initial guess
-           psVector *restrict paramMask //!< 1 = fit for parameter, 0 = hold parameter constant
-          );
+psMinimize(psVector *restrict initialGuess,
+           float (*myFunction)(const psVector *restrict, const psVector *restrict),
+           float (*myFunctionDeriv)(const psVector *restrict, const psVector *restrict, int),
+           const psVector *restrict coord,
+           const psVector *restrict paramMask);
 
 
Index: /trunk/psLib/src/math/psMinimize.c
===================================================================
--- /trunk/psLib/src/math/psMinimize.c	(revision 1184)
+++ /trunk/psLib/src/math/psMinimize.c	(revision 1185)
@@ -7,4 +7,5 @@
 
 #include <gsl/gsl_multifit_nlin.h>
+#include <gsl/gsl_multimin.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
@@ -38,5 +39,4 @@
 }
 psModelData;
-
 // The first argument to evalModel() and
 // d_evalModel() specifies the data point.
@@ -47,17 +47,241 @@
 
 
+typedef struct
+{
+    int paramCount;   // Number of non-masked parameters.
+    psVector *restrict       initialGuess;
+    const psVector *restrict  coord;
+    const psVector *restrict paramMask;
+    float (*evalModel) (const psVector *, const psVector *);
+    float (*d_evalModel) (const psVector *, const psVector *, int);
+}
+psModelData2;
+
 
 /******************************************************************************
-    This routine must minimize an arbitrary function.
+p_psMinFunc(*v, *params): This routine calls the user supplied function to
+be minimized.  The "v" argument of this function corresponds to the
+parameters of the function to be minimized.
+ 
+ 
+ *****************************************************************************/
+double p_psMinFunc(const gsl_vector *v,
+                   void *params)
+{
+    int i;    // Loop index variable.
+    int j;    // Loop index variable.
+    float tmpf;    // Temporary floating point variable.
+    const psVector *restrict coord   = ((psModelData2 *) params)->coord;
+    const psVector *restrict mask     = ((psModelData2 *) params)->paramMask;
+    psVector *restrict initialGuess = ((psModelData2 *)params)->initialGuess;
+    float (*evalModel)(const psVector *, const psVector *) =
+        ((psModelData2 *) params)->evalModel;
+    psVector *inputParameterList = NULL;
+
+    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
+
+    if (mask != NULL) {
+        j = 0;
+        for (i=0;i<mask->n;i++) {
+            if (mask->data.U8[i] != 0) {
+                inputParameterList->data.F32[i] = initialGuess->data.F32[i];
+            } else {
+                inputParameterList->data.F32[i] = gsl_vector_get(v, j++);
+            }
+        }
+    } else {
+        for (i=0;i<initialGuess->n;i++) {
+            inputParameterList->data.F32[i] = gsl_vector_get(v, i);
+        }
+    }
+    //printf("HMMM: (%.1f %.1f %.1f %.1f)\n",
+    //       inputParameterList->data.F32[0],
+    //       inputParameterList->data.F32[1],
+    //       coord->data.F32[0],
+    //       coord->data.F32[1]);
+
+    tmpf = evalModel(inputParameterList, coord);
+    psFree(inputParameterList);
+    //printf("Called p_psMinFunc()\n");
+    return(tmpf);
+}
+
+void p_psMinFuncDeriv(const gsl_vector *v,
+                      void *params,
+                      gsl_vector *df)
+{
+    int i;    // Loop index variable.
+    int j;    // Loop index variable.
+    float tmpf;    // Temporary floating point variable.
+    const psVector *restrict coord   = ((psModelData2 *) params)->coord;
+    const psVector *restrict mask     = ((psModelData2 *) params)->paramMask;
+    psVector *restrict initialGuess = ((psModelData2 *)params)->initialGuess;
+    float (*d_evalModel)(const psVector *, const psVector *, int) =
+        ((psModelData2 *) params)->d_evalModel;
+    psVector *inputParameterList = NULL;
+
+    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
+
+    if (mask != NULL) {
+        j = 0;
+        for (i=0;i<mask->n;i++) {
+            if (mask->data.U8[i] != 0) {
+                inputParameterList->data.F32[i] = initialGuess->data.F32[i];
+            } else {
+                inputParameterList->data.F32[i] = gsl_vector_get(v, j++);
+            }
+        }
+    } else {
+        for (i=0;i<initialGuess->n;i++) {
+            inputParameterList->data.F32[i] = gsl_vector_get(v, i);
+        }
+    }
+
+    for (i=0;i<initialGuess->n;i++) {
+        tmpf = d_evalModel(inputParameterList, coord, i);
+        gsl_vector_set(df, i, tmpf);
+    }
+    psFree(inputParameterList);
+}
+
+/* Compute both f and df together. */
+void p_psMinFuncFuncDeriv(const gsl_vector *x,
+                          void *params,
+                          double *f,
+                          gsl_vector *df)
+{
+    *f = p_psMinFunc(x, params);
+    p_psMinFuncDeriv(x, params, df);
+}
+
+
+/******************************************************************************
+psMinimize(initialGuess, myFunction, myFunctionDeriv, coord, paramMask):
+ 
+This routine must minimize an arbitrary function; it must determine the set
+of parameters of that function such that the 
+ 
  *****************************************************************************/
 psVector *
-psMinimize(float (*myFunction)(const psVector *restrict),
-           psVector *restrict initialGuess,
-           psVector *restrict paramMask)
-{
-    return(NULL);
-}
-
-
+psMinimize(psVector *restrict initialGuess,
+           float (*myFunction)(const psVector *restrict, const psVector *restrict),
+           float (*myFunctionDeriv)(const psVector *restrict, const psVector *restrict, int),
+           const psVector *restrict coord,
+           const psVector *restrict paramMask)
+{
+    int status;
+    int i = 0;
+    int j = 0;
+    int iter = 0;
+    gsl_multimin_function_fdf f;
+    double *xInit = NULL;
+    const gsl_multimin_fdfminimizer_type *T;
+    gsl_multimin_fdfminimizer *s;
+    psModelData2 inputData;
+    gsl_vector *x;
+
+    inputData.initialGuess = initialGuess;
+    inputData.coord = coord;
+    inputData.paramMask = paramMask;
+    inputData.evalModel = myFunction;
+    inputData.d_evalModel = myFunctionDeriv;
+    inputData.paramCount = 0;
+
+    printf("psMinimize(): initialGuess->n is %d\n", initialGuess->n);
+    printf("psMinimize(): coord->n is %d\n", coord->n);
+    printf("psMinimize(): coord is (%.1f, %.1f)\n", coord->data.F32[0],
+           coord->data.F32[1]);
+
+    if (paramMask != NULL) {
+        for (i=0;i<paramMask->n;i++) {
+            if (paramMask->data.U8[i] != 0) {
+                inputData.paramCount++;
+            }
+        }
+    } else {
+        inputData.paramCount= initialGuess->n;
+    }
+
+    // The initial guess at the parameters for the function are written into
+    // the vector inputParameterList.  If the paramMask is not NULL, then those
+    // parameters are masked out.
+
+    xInit = (double *) psAlloc(inputData.paramCount * sizeof(double));
+    if (paramMask != NULL) {
+        j = 0;
+        for (i=0;i<initialGuess->n;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                xInit[j++] = initialGuess->data.F32[i];
+            }
+        }
+    } else {
+        for (i=0;i<initialGuess->n;i++) {
+            xInit[i] = initialGuess->data.F32[i];
+        }
+    }
+    f.f = &p_psMinFunc;
+    f.df = &p_psMinFuncDeriv;
+    f.fdf = &p_psMinFuncFuncDeriv;
+    f.n = inputData.paramCount;
+    // GUS: is this correct?
+    //f.params = xInit;
+    f.params = &inputData;
+
+
+    printf("inputData.paramCount is %d\n", inputData.paramCount);
+    x = gsl_vector_alloc(inputData.paramCount);
+    for (i=0;i<inputData.paramCount;i++) {
+        gsl_vector_set(x, i, xInit[i]);
+    }
+    T = gsl_multimin_fdfminimizer_conjugate_fr;
+    s = gsl_multimin_fdfminimizer_alloc(T, inputData.paramCount);
+
+    gsl_multimin_fdfminimizer_set(s, &f, x, 0.01, 1e-4);
+
+    do {
+        iter++;
+        status = gsl_multimin_fdfminimizer_iterate(s);
+
+        if (status)
+            break;
+
+        status = gsl_multimin_test_gradient(s->gradient, 1e-3);
+
+        if (status == GSL_SUCCESS)
+            printf ("Minimum found at:\n");
+
+    } while (status == GSL_CONTINUE && iter < 100);
+
+    // In the above steps we had removed the masked elements from the
+    // the solver.  This next code blocks puts those masked elements
+    // into the solution.
+    if (paramMask != NULL) {
+        j = 0;
+        for (i=0;i<initialGuess->n;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                initialGuess->data.F32[i] = gsl_vector_get(s->x, j++);
+            } else {
+                initialGuess->data.F32[i] = initialGuess->data.F32[i];
+            }
+        }
+    } else {
+        for (i=0;i<initialGuess->n;i++) {
+            initialGuess->data.F32[i] = gsl_vector_get(s->x, i);
+        }
+    }
+    psFree(xInit);
+    return(initialGuess);
+}
+
+
+
+
+
+// The first argument to evalModel() and
+// d_evalModel() specifies the data point.
+// It must have the same size as the second
+// dimension of *domain.
+// The second argument must have the same size
+// as *initialGuess and *paramMask.
 
 /******************************************************************************
Index: /trunk/psLib/src/math/psMinimize.h
===================================================================
--- /trunk/psLib/src/math/psMinimize.h	(revision 1184)
+++ /trunk/psLib/src/math/psMinimize.h	(revision 1185)
@@ -14,8 +14,9 @@
 /** This routine must minimize a non-linear function */
 psVector *
-psMinimize(float (*myFunction)(const psVector *restrict), ///< Function to minimize
-           psVector *restrict initialGuess, ///< Initial guess
-           psVector *restrict paramMask //!< 1 = fit for parameter, 0 = hold parameter constant
-          );
+psMinimize(psVector *restrict initialGuess,
+           float (*myFunction)(const psVector *restrict, const psVector *restrict),
+           float (*myFunctionDeriv)(const psVector *restrict, const psVector *restrict, int),
+           const psVector *restrict coord,
+           const psVector *restrict paramMask);
 
 
Index: /trunk/psLib/test/dataManip/Makefile
===================================================================
--- /trunk/psLib/test/dataManip/Makefile	(revision 1184)
+++ /trunk/psLib/test/dataManip/Makefile	(revision 1185)
@@ -3,6 +3,6 @@
 ##  Makefile:   test/sysUtils
 ##
-##  $Revision: 1.29 $  $Name: not supported by cvs2svn $
-##  $Date: 2004-07-02 01:08:54 $
+##  $Revision: 1.30 $  $Name: not supported by cvs2svn $
+##  $Date: 2004-07-07 02:38:32 $
 ##
 ##  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -40,6 +40,8 @@
  tst_psFunc00 \
  tst_psFunc01 \
+ tst_psMinimize00 \
  tst_psMinimize01 \
- tst_psMinimize00 \
+ tst_psMinimize02 \
+ tst_psMinimize03 \
  tst_psImageStats00 \
  tst_psImageStats01 \
@@ -77,4 +79,6 @@
 tst_psMinimize00:	tst_psMinimize00.o
 tst_psMinimize01:	tst_psMinimize01.o
+tst_psMinimize02:	tst_psMinimize02.o
+tst_psMinimize03:	tst_psMinimize03.o
 tst_psImageStats00:	tst_psImageStats00.o
 tst_psImageStats01:	tst_psImageStats01.o
Index: /trunk/psLib/test/dataManip/tst_psMinimize02.c
===================================================================
--- /trunk/psLib/test/dataManip/tst_psMinimize02.c	(revision 1185)
+++ /trunk/psLib/test/dataManip/tst_psMinimize02.c	(revision 1185)
@@ -0,0 +1,164 @@
+/*****************************************************************************
+    This routine must ensure that the psMinimize function correctly
+    minimizes an arbitrary function.
+ *****************************************************************************/
+#include <stdio.h>
+#include "pslib.h"
+#include "psTest.h"
+#include "psMemory.h"
+#include "psVector.h"
+#include "psImage.h"
+#include "psFunctions.h"
+#include "psMinimize.h"
+#include <math.h>
+#define DATA_WIDTH 2
+#define NUM_PARAMS 2
+
+float myFunc0(const psVector *restrict myParams,
+              const psVector *restrict myCoords)
+{
+    float P0 = myCoords->data.F32[0];
+    float P1 = myCoords->data.F32[1];
+    float x = myParams->data.F32[0];
+    float y = myParams->data.F32[1];
+    float tmp = 0.0;
+
+    tmp = 10.0 * (x - P0) * (x - P0) + 20.0 * (y - P1) * (y - P1) + 30.0;
+    //    printf("--------- myFunc((%.1f %.1f) %.1f %.1f) is %.1f ---------\n",
+    //    printf("--------- myFunc() params: (%.1f %.1f) coords: (%.1f %.1f) is %.1f ---------\n",
+    //           x, y, P0, P1, tmp);
+
+    return(tmp);
+}
+
+float myFuncDeriv0(const psVector *restrict myParams,
+                   const psVector *restrict myCoords,
+                   int whichParamDeriv)
+{
+    float P0 = myCoords->data.F32[0];
+    float P1 = myCoords->data.F32[1];
+    float x = myParams->data.F32[0];
+    float y = myParams->data.F32[1];
+    float tmp = 0.0;
+
+    if (whichParamDeriv == 0) {
+        tmp = 20.0 * (x - P0);
+    } else if (whichParamDeriv == 1) {
+        tmp = 40.0 * (y - P1);
+    }
+
+    //    printf("--------- myFuncDeriv((%.1f %.1f) %.1f %.1f (%d)) is %.1f ---------\n",
+    //            x, y, myParams->data.F32[0], myParams->data.F32[1], whichParamDeriv, tmp);
+    return(tmp);
+}
+
+
+float myFunc(const psVector *restrict myParams,
+             const psVector *restrict myCoords)
+{
+    float x = myCoords->data.F32[0];
+    float y = myCoords->data.F32[1];
+    float A = myParams->data.F32[0];
+    float B = myParams->data.F32[1];
+    float C = myParams->data.F32[2];
+    float D = myParams->data.F32[3];
+    float tmp = 0.0;
+
+    tmp = (A * x) + (B*x*x) + (C*y*y) + (D*x*x*y);
+    printf("--------- myFunc((%.1f %.1f) %.1f %.1f %.1f %.1f) is %.1f ---------\n",
+           x, y, A, B, C, D, tmp);
+
+    return(tmp);
+}
+
+float myFuncDeriv(const psVector *restrict myParams,
+                  const psVector *restrict myCoords,
+                  int whichParamDeriv)
+{
+    float x = myCoords->data.F32[0];
+    float y = myCoords->data.F32[1];
+    float tmp = 0.0;
+
+    if (whichParamDeriv == 0) {
+        tmp = 1.0;
+        tmp = x;
+    } else if (whichParamDeriv == 1) {
+        tmp = x;
+        tmp = x*x;
+    } else if (whichParamDeriv == 2) {
+        tmp = y;
+        tmp = y*y;
+    } else if (whichParamDeriv == 3) {
+        tmp = x * x * y;
+    }
+
+    printf("--------- myFuncDeriv((%.1f %.1f) %.1f %.1f %.1f %.1f (%d)) is %.1f ---------\n",
+           x, y, myParams->data.F32[0], myParams->data.F32[1],
+           myParams->data.F32[2], myParams->data.F32[3], whichParamDeriv, tmp);
+    return(tmp);
+}
+
+
+int main()
+{
+    psVector *initialGuess = NULL;
+    psVector *paramMask = NULL;
+    psVector *coord = NULL;
+    int i = 0;
+    int currentId = psMemGetId();
+    int testStatus = true;
+    int memLeaks = 0;
+    psVector *theParams = NULL;
+
+    initialGuess = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);
+    paramMask = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);
+    coord = psVectorAlloc(DATA_WIDTH, PS_TYPE_F32);
+
+
+    //--------- myFunc((0.0 1.0) -8.6 -7.6 14.6 5.0) is 6.0 ---------
+
+    // Build the data.
+    initialGuess->data.F32[0] = 5.0;
+    initialGuess->data.F32[1] = 7.0;
+
+    for (i=0;i<DATA_WIDTH;i++) {
+        coord->data.F32[i] = (float) i;
+    }
+    coord->data.F32[0] = 1.0;
+    coord->data.F32[1] = 2.0;
+
+    printPositiveTestHeader(stdout,
+                            "psMinimize functions",
+                            "psMinimize(): no masks");
+
+    for (i=0;i<NUM_PARAMS;i++)
+        printf("Initial parameter %d is: %.1f\n", i, initialGuess->data.F32[i]);
+
+    theParams = psMinimize(initialGuess,
+                           myFunc0,
+                           myFuncDeriv0,
+                           coord,
+                           NULL);
+
+    for (i=0;i<NUM_PARAMS;i++)
+        printf("Parameter %d is: %.1f\n", i, theParams->data.F32[i]);
+
+    psMemCheckCorruption(1);
+    printFooter(stdout,
+                "psMinimize functions",
+                "psMinimize(): no masks",
+                testStatus);
+
+    psMemCheckCorruption(1);
+    psFree(initialGuess);
+    psFree(paramMask);
+    psFree(coord);
+
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,NULL);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+
+    return (!testStatus);
+}
Index: /trunk/psLib/test/dataManip/tst_psMinimize03.c
===================================================================
--- /trunk/psLib/test/dataManip/tst_psMinimize03.c	(revision 1185)
+++ /trunk/psLib/test/dataManip/tst_psMinimize03.c	(revision 1185)
@@ -0,0 +1,168 @@
+/*****************************************************************************
+    This routine must ensure that the psMinimize function correctly
+    minimizes an arbitrary function.
+ *****************************************************************************/
+#include <stdio.h>
+#include "pslib.h"
+#include "psTest.h"
+#include "psMemory.h"
+#include "psVector.h"
+#include "psImage.h"
+#include "psFunctions.h"
+#include "psMinimize.h"
+#include <math.h>
+#define DATA_WIDTH 2
+#define NUM_PARAMS 4
+
+float myFunc0(const psVector *restrict myParams,
+              const psVector *restrict myCoords)
+{
+    float P0 = myCoords->data.F32[0];
+    float P1 = myCoords->data.F32[1];
+    float x = myParams->data.F32[0];
+    float y = myParams->data.F32[1];
+    float tmp = 0.0;
+
+    tmp = 10.0 * (x - P0) * (x - P0) + 20.0 * (y - P1) * (y - P1) + 30.0;
+    //    printf("--------- myFunc((%.1f %.1f) %.1f %.1f) is %.1f ---------\n",
+    //    printf("--------- myFunc() params: (%.1f %.1f) coords: (%.1f %.1f) is %.1f ---------\n",
+    //           x, y, P0, P1, tmp);
+
+    return(tmp);
+}
+
+float myFuncDeriv0(const psVector *restrict myParams,
+                   const psVector *restrict myCoords,
+                   int whichParamDeriv)
+{
+    float P0 = myCoords->data.F32[0];
+    float P1 = myCoords->data.F32[1];
+    float x = myParams->data.F32[0];
+    float y = myParams->data.F32[1];
+    float tmp = 0.0;
+
+    if (whichParamDeriv == 0) {
+        tmp = 20.0 * (x - P0);
+    } else if (whichParamDeriv == 1) {
+        tmp = 40.0 * (y - P1);
+    }
+
+    //    printf("--------- myFuncDeriv((%.1f %.1f) %.1f %.1f (%d)) is %.1f ---------\n",
+    //            x, y, myParams->data.F32[0], myParams->data.F32[1], whichParamDeriv, tmp);
+    return(tmp);
+}
+
+
+float myFunc(const psVector *restrict myParams,
+             const psVector *restrict myCoords)
+{
+    float x = myCoords->data.F32[0];
+    float y = myCoords->data.F32[1];
+    float A = myParams->data.F32[0];
+    float B = myParams->data.F32[1];
+    float C = myParams->data.F32[2];
+    float D = myParams->data.F32[3];
+    float tmp = 0.0;
+
+    tmp = (A * x) + (B*x*x) + (C*y*y) + (D*x*x*y);
+    printf("-- myFunc() params: (%.1f %.1f %.1f %.1f) coords: (%.1f %.1f) is %.1f ---------\n",
+           A, B, C, D, x, y, tmp);
+
+    return(tmp);
+}
+
+float myFuncDeriv(const psVector *restrict myParams,
+                  const psVector *restrict myCoords,
+                  int whichParamDeriv)
+{
+    float x = myCoords->data.F32[0];
+    float y = myCoords->data.F32[1];
+    float tmp = 0.0;
+
+    if (whichParamDeriv == 0) {
+        tmp = 1.0;
+        tmp = x;
+    } else if (whichParamDeriv == 1) {
+        tmp = x;
+        tmp = x*x;
+    } else if (whichParamDeriv == 2) {
+        tmp = y;
+        tmp = y*y;
+    } else if (whichParamDeriv == 3) {
+        tmp = x * x * y;
+    }
+
+    printf("--------- myFuncDeriv((%.1f %.1f) %.1f %.1f %.1f %.1f (%d)) is %.1f ---------\n",
+           x, y, myParams->data.F32[0], myParams->data.F32[1],
+           myParams->data.F32[2], myParams->data.F32[3], whichParamDeriv, tmp);
+    return(tmp);
+}
+
+
+int main()
+{
+    psVector *initialGuess = NULL;
+    psVector *paramMask = NULL;
+    psVector *coord = NULL;
+    int i = 0;
+    int currentId = psMemGetId();
+    int testStatus = true;
+    int memLeaks = 0;
+    psVector *theParams = NULL;
+
+    initialGuess = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);
+    paramMask = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);
+    coord = psVectorAlloc(DATA_WIDTH, PS_TYPE_F32);
+
+
+    //--------- myFunc((0.0 1.0) -8.6 -7.6 14.6 5.0) is 6.0 ---------
+
+    // Build the data.
+    initialGuess->data.F32[0] = 2.0;
+    initialGuess->data.F32[1] = 3.0;
+    initialGuess->data.F32[2] = 4.0;
+    initialGuess->data.F32[3] = 5.0;
+
+    for (i=0;i<DATA_WIDTH;i++) {
+        coord->data.F32[i] = (float) i;
+    }
+    coord->data.F32[0] = 1.0;
+    coord->data.F32[1] = 2.0;
+
+    printPositiveTestHeader(stdout,
+                            "psMinimize functions",
+                            "psMinimize(): no masks");
+
+    for (i=0;i<NUM_PARAMS;i++)
+        printf("Initial parameter %d is: %.1f\n", i, initialGuess->data.F32[i]);
+
+    theParams = psMinimize(initialGuess,
+                           myFunc,
+                           myFuncDeriv,
+                           coord,
+                           NULL);
+
+    for (i=0;i<NUM_PARAMS;i++)
+        printf("Parameter %d is: %.1f\n", i, theParams->data.F32[i]);
+
+    printf("Function is %.1f\n", myFunc(theParams, coord));
+
+    psMemCheckCorruption(1);
+    printFooter(stdout,
+                "psMinimize functions",
+                "psMinimize(): no masks",
+                testStatus);
+
+    psMemCheckCorruption(1);
+    psFree(initialGuess);
+    psFree(paramMask);
+    psFree(coord);
+
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,NULL);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+
+    return (!testStatus);
+}