add memory handling for DKS/GPU. No fitting possible, yet.

src/classes/PFitterFcnDKS.cpp

@@ -81,6 +81,17 @@ Double_t PFitterFcnDKS::operator()(const std::vector<Double_t>& par) const
 {
   Double_t value = 0.0;
 
+  // write parameter to GPU
+  Int_t ierr = fDKS.writeParams(&par[0], par.size());
+
+  // loop over all data sets
+  for (UInt_t i=0; i<fRunListCollection->GetNoOfSingleHisto(); i++) {
+    // set N0, Nbkg
+    // set fun values
+    // set map values
+    // calc chisq/log-mlh
+  }
+
   return value;
 }
 
@@ -122,6 +133,8 @@ void PFitterFcnDKS::CalcExpectedChiSquare(const std::vector<Double_t> &par, Doub
  */
 void PFitterFcnDKS::InitDKS(const UInt_t dksTag)
 {
+  Int_t ierr = 0;
+
   // if any device was allocated before, free the device resources
   FreeDKS();
 
@@ -142,16 +155,89 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag)
 
   // init chisq buffer on the GPU
 
+  // 1) calculated the maximum size for the data needed.
+  Int_t maxSize = -1, size = -1;
+  Int_t parSize = -1, mapSize = -1, funSize = -1;
+  for (UInt_t i=0; i<fRunListCollection->GetNoOfSingleHisto(); i++) {
+    size = fRunListCollection->GetSingleHisto(i)->GetValue()->size();
+    if (maxSize < size)
+      maxSize = size;
+  }
+  for (UInt_t i=0; i<fRunListCollection->GetNoOfAsymmetry(); i++) {
+    size = fRunListCollection->GetAsymmetry(i)->GetValue()->size();
+    if (maxSize < size)
+      maxSize = size;
+  }
+  if (maxSize == -1) {
+    cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get data size to be fitted." << endl;
+    fValid = false;
+    return;
+  }
+
+  // 2) get number of parameters / functions / maps
+  parSize = fRunListCollection->GetNoOfParameters();
+  if (parSize == -1) {
+    cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get number of fit parameters." << endl;
+    fValid = false;
+    return;
+  }
+  funSize = fRunListCollection->GetNoOfFunctions();
+  if (funSize == -1) {
+    cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get number of functions." << endl;
+    fValid = false;
+    return;
+  }
+  mapSize = fRunListCollection->GetNoOfMaps();
+  if (mapSize == -1) {
+    cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get number of maps." << endl;
+    fValid = false;
+    return;
+  }
+  cout << "debug> parSize=" << parSize << ", funSize=" << funSize << ", mapSize=" << mapSize << endl;
+
+  // now ready to init the chisq buffer on the GPU
+  ierr = fDKS.initChiSquare(maxSize, parSize, funSize, mapSize);
+  if (ierr != 0) {
+    cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to allocate the necessary chisq buffer on the GPU." << endl;
+    fValid = false;
+    return;
+  }
+
   // allocate memory for the data on the GPU/CPU and transfer the data sets
+  fMemData.resize(fRunListCollection->GetNoOfSingleHisto());
+  PRunData *runData=0;
+  for (UInt_t i=0; i<fRunListCollection->GetNoOfSingleHisto(); i++) {
+    runData = fRunListCollection->GetSingleHisto(i);
+    if (runData == 0) {
+      cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to get data set (i=" << i << ") from fRunListCollection." << endl;
+      fValid = false;
+      return;
+    }
+    fMemData[i] = fDKS.allocateMemory<Double_t>(runData->GetValue()->size(), ierr);
+    if (ierr != 0) {
+      cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to allocated data set memory (i=" << i << ") on the GPU" << endl;
+      fValid = false;
+      return;
+    }
+    fDKS.writeData<double>(fMemData[i], runData->GetValue(), runData->GetValue()->size());
+  }
 
   // set the function string and compile the program
-  Int_t ierr = fDKS.callCompileProgram(fTheoStr, !fUseChi2);
+  ierr = fDKS.callCompileProgram(fTheoStr, !fUseChi2);
   if (ierr != 0) {
     cerr << ">> PFitterFcnDKS::InitDKS: **ERROR** failed to compile theory!" << endl;
     fValid = false;
     return;
   }
 
+  // keep the N0 and Nbkg for each data set
+  fNidx.resize(fRunListCollection->GetNoOfSingleHisto());
+  for (UInt_t i=0; i<fRunListCollection->GetNoOfSingleHisto(); i++) {
+    fNidx[i].fN0   = fRunListCollection->GetN0Idx(i);
+    fNidx[i].fNbkg = fRunListCollection->GetNbkgIdx(i);
+    cout << "debug> runNo=" << i << ", N0 idx=" << fNidx[i].fN0 << ", Nbkg idx=" << fNidx[i].fNbkg << endl;
+  }
+
   // checks device properties if openCL
   ierr = fDKS.checkMuSRKernels();
   if (ierr != 0) {
@@ -172,5 +258,10 @@ void PFitterFcnDKS::InitDKS(const UInt_t dksTag)
  */
 void PFitterFcnDKS::FreeDKS()
 {
-
+  PRunData *runData=0;
+  for (UInt_t i=0; i<fMemData.size(); i++) {
+    runData = fRunListCollection->GetSingleHisto(i);
+    fDKS.freeMemory<Double_t>(fMemData[i], runData->GetValue()->size());
+  }
+  fMemData.clear();
 }

src/classes/PMsrHandler.cpp

@@ -69,6 +69,8 @@ PMsrHandler::PMsrHandler(const Char_t *fileName, PStartupOptions *startupOptions
 
   fFuncHandler = 0;
 
+  fNoOfMaps = -1;
+
   // check if the file name given is a path-file-name, and if yes, split it into path and file name.
   if (fFileName.Contains("/")) {
     Int_t idx = -1;
@@ -5664,6 +5666,17 @@ Bool_t PMsrHandler::CheckMaps()
     }
   }
 
+  if (result == true) {
+    PIntVector *map = fRuns[0].GetMap();
+    fNoOfMaps = 0;
+    for (UInt_t i=0; i<map->size(); i++)
+      if (map->at(i) != 0)
+        fNoOfMaps++;
+    if (fNoOfMaps == 0)
+      fNoOfMaps = -1;
+  }
+  cout << "debug> fNoOfMaps=" << fNoOfMaps << endl;
+
   // clean up
   mapVec.clear();
   mapBlock.clear();

src/classes/PRunListCollection.cpp

@@ -1062,3 +1062,44 @@ const Char_t* PRunListCollection::GetYAxisTitle(const TString &runName, const UI
   return result;
 }
 
+//--------------------------------------------------------------------------
+// GetN0Idx (public)
+//--------------------------------------------------------------------------
+/**
+ * \brief PRunListCollection::GetN0Idx
+ * \param idx
+ * \return
+ */
+Int_t PRunListCollection::GetN0Idx(UInt_t idx)
+{
+  Int_t N0idx = -1;
+
+  // make sure idx is within proper bounds
+  if (idx >= fRunSingleHistoList.size())
+    return N0idx;
+
+  N0idx = fMsrInfo->GetMsrRunList()->at(idx).GetNormParamNo();
+
+  return N0idx;
+}
+
+//--------------------------------------------------------------------------
+// GetNbkgIdx (public)
+//--------------------------------------------------------------------------
+/**
+ * \brief PRunListCollection::GetNbkgIdx
+ * \param idx
+ * \return
+ */
+Int_t PRunListCollection::GetNbkgIdx(UInt_t idx)
+{
+  Int_t NbkgIdx = -1;
+
+  // make sure idx is within proper bounds
+  if (idx >= fRunSingleHistoList.size())
+    return NbkgIdx;
+
+  NbkgIdx = fMsrInfo->GetMsrRunList()->at(idx).GetBkgFitParamNo();
+
+  return NbkgIdx;
+}

src/include/PFitterFcnDKS.h

@@ -36,6 +36,11 @@
 #include "DKSBaseMuSR.h"
 #include "PRunListCollection.h"
 
+typedef struct {
+  UInt_t fN0;   ///< N0 parameter index
+  UInt_t fNbkg; ///< Nbkg parameter index
+} PNidx;
+
 /**
  * <p>This is the minuit2 interface function class providing the function to be optimized (chisq or log max-likelihood).
  */
@@ -63,6 +68,9 @@ class PFitterFcnDKS : public ROOT::Minuit2::FCNBase
 
     mutable DKSBaseMuSR fDKS;
 
+    vector<void *> fMemData; ///< vector holding the initial addresses of the data sets on the GPU
+    vector<PNidx>  fNidx;    ///< N0 / Nbkg parameter index vector
+
     virtual void InitDKS(const UInt_t dksTag);
     virtual void FreeDKS();
 };

src/include/PMsrHandler.h

@@ -113,6 +113,7 @@ class PMsrHandler
 
     virtual std::string GetDKSTheoryString();
     virtual UInt_t GetDKSTag();
+    virtual Int_t GetNoOfMaps() { return fNoOfMaps; }
 
   private:
     Bool_t                 fFourierOnly; ///< flag indicating if Fourier transform only is wished. If yes, some part of the msr-file blocks are not needed.
@@ -139,6 +140,8 @@ class PMsrHandler
 
     Bool_t fCopyStatisticsBlock; ///< flag, if true: just copy to old statistics block (musrt0), otherwise write a new one (musrfit)
 
+    Int_t fNoOfMaps;
+
     virtual Bool_t HandleFitParameterEntry(PMsrLines &line);
     virtual Bool_t HandleTheoryEntry(PMsrLines &line);
     virtual Bool_t HandleFunctionsEntry(PMsrLines &line);

src/include/PRunListCollection.h

@@ -100,6 +100,12 @@ class PRunListCollection
     virtual const Char_t* GetXAxisTitle(const TString &runName, const UInt_t idx) const;
     virtual const Char_t* GetYAxisTitle(const TString &runName, const UInt_t idx) const;
 
+    virtual Int_t GetNoOfParameters() { return fMsrInfo->GetNoOfParams(); }
+    virtual Int_t GetNoOfFunctions() { return fMsrInfo->GetNoOfFuncs(); }
+    virtual Int_t GetNoOfMaps() { return fMsrInfo->GetNoOfMaps(); }
+    virtual Int_t GetN0Idx(UInt_t idx);
+    virtual Int_t GetNbkgIdx(UInt_t idx);
+
   private:
     PMsrHandler *fMsrInfo;  ///< pointer to the msr-file handler
     PRunDataHandler *fData; ///< pointer to the run-data handler