/[escript]/trunk/escript/src/DataExpanded.cpp
ViewVC logotype

Diff of /trunk/escript/src/DataExpanded.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3468 by jfenwick, Tue Feb 22 06:38:57 2011 UTC revision 3506 by jfenwick, Wed May 11 01:59:45 2011 UTC
# Line 17  Line 17 
17  #include "DataException.h"  #include "DataException.h"
18  #include "DataConstant.h"  #include "DataConstant.h"
19  #include "DataTagged.h"  #include "DataTagged.h"
20    #include <limits>
21    
22  #include "esysUtils/Esys_MPI.h"  #include "esysUtils/Esys_MPI.h"
23    
# Line 27  Line 28 
28  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
29  #include "DataMaths.h"  #include "DataMaths.h"
30    
31    //#define MKLRANDOM
32    
33    #ifdef MKLRANDOM
34    #include <mkl_vsl.h>
35    
36    #endif
37    
38  using namespace std;  using namespace std;
39  using namespace boost::python;  using namespace boost::python;
40  using namespace boost;  using namespace boost;
# Line 837  DataExpanded::getVectorRO() const Line 845  DataExpanded::getVectorRO() const
845      return m_data.getData();      return m_data.getData();
846  }  }
847    
848  void DataExpanded::randomFill(double seed)  
849    // Idea here is to create an array of seeds by feeding the original seed into the random generator
850    // The code at the beginning of the function to compute the seed if one is given is
851    // just supposed to introduce some variety (and ensure that multiple ranks don't get the same seed).
852    // I make no claim about how well these initial seeds are distributed
853    void DataExpanded::randomFill(long seed)
854  {  {
855      CHECK_FOR_EX_WRITE      CHECK_FOR_EX_WRITE
856      if (seed==0)      static unsigned prevseed=0; // So if we create a bunch of objects we don't get the same start seed
857        if (seed==0)        // for each one
858        {
859        if (prevseed==0)
860        {
861            time_t s=time(0);
862            seed=s;
863        }
864        else
865        {
866            seed=prevseed+419;  // these numbers are arbitrary
867            if (seed>3040101)       // I want to avoid overflow on 32bit systems
868            {
869            seed=((int)(seed)%0xABCD)+1;
870            }
871        }
872        }
873        // now we need to consider MPI since we don't want each rank to start with the same seed
874        seed+=getFunctionSpace().getDomain()->getMPIRank()*getFunctionSpace().getDomain()->getMPISize()*3;
875        prevseed=seed;
876    #ifdef _OPENMP
877        int numthreads=omp_get_max_threads();
878    #else
879        int numthreads=1;
880    #endif
881    
882    #ifdef MKLRANDOM
883        double* seeds=new double[numthreads];
884        VSLStreamStatePtr sstream;
885    
886        int status=vslNewStream(&sstream, VSL_BRNG_MT19937, seed);  // use a Mersenne Twister
887        numeric_limits<double> dlim;
888        vdRngUniform(VSL_METHOD_DUNIFORM_STD, sstream , numthreads, seeds, -1, 1);
889        vslDeleteStream(&sstream);
890        DataVector& dv=getVectorRW();
891        size_t dvsize=dv.size();
892        #pragma omp parallel
893      {      {
894      time_t s=time(0);      int tnum=0;
895      seed=s;      #ifdef _OPENMP
896        tnum=omp_get_thread_num();
897        #endif
898        VSLStreamStatePtr stream;
899        // the 12345 is a hack to give us a better chance of getting different integer seeds.
900            int status=vslNewStream(&stream, VSL_BRNG_MT19937, seeds[tnum]*12345);  // use a Mersenne Twister
901        int bigchunk=(dvsize/numthreads+1);
902        int smallchunk=dvsize-bigchunk*(numthreads-1);
903        int chunksize=(tnum<(numthreads-1))?bigchunk:smallchunk;
904            vdRngUniform(VSL_METHOD_DUNIFORM_STD, stream, chunksize, &(dv[bigchunk*tnum]), 0,1);
905            vslDeleteStream(&stream);
906      }      }
907        delete[] seeds;
908    #else
909      srand(seed);      srand(seed);
910        unsigned* seeds=new unsigned[numthreads];
911        for (int i=0;i<numthreads;++i)
912        {
913        seeds[i]=rand();
914        }
915      DataVector&  dv=getVectorRW();      DataVector&  dv=getVectorRW();
916      for (long i=0;i<dv.size();++i)      long i;
917        const size_t dvsize=dv.size();
918        #pragma omp parallel private(i)
919      {      {
920      dv[i]=(double)rand()/RAND_MAX;      int tnum=0;
921        #ifdef _OPENMP
922        tnum=omp_get_thread_num();
923        #endif
924        unsigned info=seeds[tnum];
925        
926            #pragma omp for schedule(static)
927            for (i=0;i<dvsize;++i)
928            {
929            dv[i]=(double)rand_r(&info)/RAND_MAX;
930            }
931      }      }
932        delete[] seeds;
933    #endif
934  }  }
935    
936  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.3468  
changed lines
  Added in v.3506

  ViewVC Help
Powered by ViewVC 1.1.26