/[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 3317 by caltinay, Thu Oct 28 00:50:41 2010 UTC revision 3544 by jfenwick, Wed Jul 27 06:19:10 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 161  DataExpanded::DataExpanded(const Functio Line 169  DataExpanded::DataExpanded(const Functio
169    
170  }  }
171    
172    DataExpanded::DataExpanded(const FunctionSpace& what,
173                               const DataTypes::ShapeType &shape,
174                               const double v)
175      : parent(what,shape)
176    {
177         ValueType& vec=m_data.getData();
178         //
179         // create the view of the data
180         initialise(what.getNumSamples(),what.getNumDPPSample());
181         // now we copy this value to all elements
182         const int L=getLength();
183         int i;
184         #pragma omp parallel for schedule(static) private(i)
185         for (i=0;i<L;++i)
186         {
187            vec[i]=v;
188         }
189    }
190    
191    
192    
193  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
194  {  {
195  }  }
# Line 817  DataExpanded::getVectorRO() const Line 846  DataExpanded::getVectorRO() const
846  }  }
847    
848    
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
856        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        int tnum=0;
895        #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);
910        unsigned* seeds=new unsigned[numthreads];
911        for (int i=0;i<numthreads;++i)
912        {
913        seeds[i]=rand();
914        }
915        DataVector&  dv=getVectorRW();
916        long i;
917        const size_t dvsize=dv.size();
918        #pragma omp parallel private(i)
919        {
920        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    #ifdef _WIN32
930            dv[i]=((double)rand())/RAND_MAX;
931    #else
932            dv[i]=((double)rand_r(&info))/RAND_MAX;
933    #endif
934            }
935        }
936        delete[] seeds;
937    #endif
938    }
939    
940  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.3317  
changed lines
  Added in v.3544

  ViewVC Help
Powered by ViewVC 1.1.26