Index: trunk/psModules/src/objects/pmSourceFitSet.c
===================================================================
--- trunk/psModules/src/objects/pmSourceFitSet.c	(revision 20937)
+++ trunk/psModules/src/objects/pmSourceFitSet.c	(revision 21163)
@@ -6,6 +6,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-12-08 02:51:14 $
+ *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-01-24 20:52:26 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -45,5 +45,30 @@
 /********************* Source Model Set Functions ***************************/
 
-static pmSourceFitSetData *thisSet = NULL;
+// these functions currently need to use this static variable because of the way psMinimizeLMM
+// is implemented.  We could re-work that structure, but for now it is probably easier to make
+// this thread safe by pre-allocating separate static variables for each thread.
+
+static psArray *fitSets = NULL;
+static pthread_mutex_t fitSetInitMutex = PTHREAD_MUTEX_INITIALIZER;
+
+// call this before launching the threads
+bool pmSourceFitSetInit (int nThreads) {
+
+    if (!fitSets) {
+	fitSets = psArrayAlloc (PS_MAX (1, nThreads));
+    }
+
+    // the allocated elements should be NULL on psArrayAlloc, 
+    // and a previously allocated array of fitSets should have been cleared 
+    // before pmSourceFitSetInit is called
+    for (int i = 0; i < fitSets->n; i++) {
+	psAssert (fitSets->data[i] == NULL, "failure to init or clear fitSets?");
+    }
+    return true;
+}
+
+void pmSourceFitSetDone () {
+    psFree (fitSets);
+}
 
 static void pmSourceFitSetDataFree (pmSourceFitSetData *set) {
@@ -62,8 +87,9 @@
     psMemSetDeallocator(set, (psFreeFunc) pmSourceFitSetDataFree);
 
-    set->modelSet = psMemIncrRefCounter (modelSet);
-    set->paramSet = psArrayAlloc (modelSet->n);
-    set->derivSet = psArrayAlloc (modelSet->n);
+    set->modelSet  = psMemIncrRefCounter (modelSet);
+    set->paramSet  = psArrayAlloc (modelSet->n);
+    set->derivSet  = psArrayAlloc (modelSet->n);
     set->nParamSet = 0;
+    set->thread    = pthread_self();
 
     for (int i = 0; i < modelSet->n; i++) {
@@ -88,4 +114,88 @@
 }
 
+pmSourceFitSetData *pmSourceFitSetDataSet (psArray *modelSet) {
+
+    pmSourceFitSetData *thisSet = NULL;
+
+    psAssert (fitSets, "pmSourceFitSetInit not called");
+
+    // find the fitSet used by this thread
+    pthread_t id = pthread_self();
+
+    // is our ID is already on the stack, abort.
+    // we do not need to lock to do this....
+    for (int i = 0; i < fitSets->n; i++) {
+	thisSet = fitSets->data[i];
+	if (!thisSet) continue;
+	if (thisSet->thread == id) {
+	    break;
+	}
+	thisSet = NULL;
+    }
+    psAssert (thisSet == NULL, "pmSourceFitSetDataSet() called but previous entry not cleared");
+
+    thisSet = pmSourceFitSetDataAlloc(modelSet);
+
+    pthread_mutex_lock (&fitSetInitMutex);
+	  
+    // find an entry that is NULL and place it there
+    for (int i = 0; i < fitSets->n; i++) {
+	if (fitSets->data[i]) continue;
+	fitSets->data[i] = thisSet;
+	pthread_mutex_unlock (&fitSetInitMutex);
+	return thisSet;
+    }
+    pthread_mutex_unlock (&fitSetInitMutex);
+    psAbort ("no empty slot for new pmSourceFitSetData");
+    return NULL;
+}
+
+pmSourceFitSetData *pmSourceFitSetDataGet () {
+
+    psAssert (fitSets, "pmSourceFitSetInit not called");
+
+    // find the fitSet used by this thread
+    pthread_t id = pthread_self();
+
+    // can we find our fitSet?  we do not need to lock to do this....
+    pmSourceFitSetData *thisSet = NULL;
+    for (int i = 0; i < fitSets->n; i++) {
+	thisSet = fitSets->data[i];
+	if (!thisSet) continue;
+	if (thisSet->thread == id) {
+	    break;
+	}
+	thisSet = NULL;
+    }
+    psAssert (thisSet != NULL, "pmSourceFitSetDataGet() called, but no entry found");
+
+    return thisSet;
+}
+
+void pmSourceFitSetDataClear () {
+
+    int i;
+
+    psAssert (fitSets, "pmSourceFitSetInit not called");
+
+    // find the fitSet used by this thread
+    pthread_t id = pthread_self();
+
+    // can we find our fitSet?  we do not need to lock to do this....
+    pmSourceFitSetData *thisSet = NULL;
+    for (i = 0; i < fitSets->n; i++) {
+	thisSet = fitSets->data[i];
+	if (!thisSet) continue;
+	if (thisSet->thread == id) {
+	    break;
+	}
+	thisSet = NULL;
+    }
+    psAssert (thisSet != NULL, "pmSourceFitSetDataClear() called, but no entry found");
+
+    psFree (thisSet);
+    fitSets->data[i] = NULL;
+    return;
+}
 
 // this function is called with the full set of parameters and the beta values in a single vector
@@ -93,5 +203,6 @@
                                 float *betas)
 {
-    PS_ASSERT_PTR_NON_NULL(thisSet, false);
+    PS_ASSERT_PTR_NON_NULL(fitSets, false);
+    pmSourceFitSetData *thisSet = pmSourceFitSetDataGet();
 
     // nParam is the parameter in the full sequence.  determine which single model this comes from
@@ -190,4 +301,5 @@
 }
 
+// set the model parameters for this fit set
 bool pmSourceFitSetValues (pmSourceFitSetData *set, const psVector *dparam, 
                            const psVector *param, pmSource *source,
@@ -242,5 +354,6 @@
 psF32 pmSourceFitSetFunction(psVector *deriv, const psVector *param, const psVector *x)
 {
-    PS_ASSERT_PTR_NON_NULL(thisSet, NAN);
+    pmSourceFitSetData *thisSet = pmSourceFitSetDataGet();
+
     float chisqSum = 0.0;
     float chisqOne = 0.0;
@@ -391,6 +504,6 @@
     yErr->n = nPix;
 
-    // create the FitSet and set the initial parameter guesses
-    thisSet = pmSourceFitSetDataAlloc (modelSet);
+    // create the FitSet for this thread and set the initial parameter guesses
+    pmSourceFitSetData *thisSet = pmSourceFitSetDataSet(modelSet);
 
     // define param and deriv vectors for complete set of parameters
@@ -431,6 +544,5 @@
         psFree (params);
         psFree(constraint);
-        psFree (thisSet);
-        thisSet = NULL;
+	pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets
         return(false);
     }
@@ -488,7 +600,5 @@
     psFree(params);
     psFree(dparams);
-    psFree(thisSet);
-
-    thisSet = NULL;
+    pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets
 
     bool rc = (onPic && fitStatus);