/[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 3549 by jfenwick, Wed Jul 27 06:19:10 2011 UTC revision 3550 by jfenwick, Fri Aug 19 01:51:35 2011 UTC
# Line 32  Line 32 
32    
33  #ifdef MKLRANDOM  #ifdef MKLRANDOM
34  #include <mkl_vsl.h>  #include <mkl_vsl.h>
35    #else
36    #include <boost/random/mersenne_twister.hpp>
37  #endif  #endif
38    
39  using namespace std;  using namespace std;
# Line 846  DataExpanded::getVectorRO() const Line 847  DataExpanded::getVectorRO() const
847  }  }
848    
849    
850    #ifndef MKLRANDOM
851    namespace {
852        
853    boost::mt19937 base;        // used to seed all the other generators  
854    vector<boost::mt19937> gens;
855    
856    
857    void seedGens(long seed)
858    {
859    #ifdef _OPENMP
860        int numthreads=omp_get_max_threads();
861    #else
862        int numthreads=1;
863    #endif
864        if (gens.size()==0)     // we haven't instantiated the generators yet  
865        {
866            gens.resize(numthreads);    // yes this means all the generators will be owned by one thread in a NUMA sense      
867        }                   // I don't think it is worth a more complicated solution at this point
868        if (seed!=0)
869        {
870           base.seed(seed);        
871           for (int i=0;i<numthreads;++i)
872           {
873                gens[i].seed(base());   // initialise each generator with successive random values      
874           }      
875        }
876    }
877      
878      
879    }
880    #endif
881    
882  // Idea here is to create an array of seeds by feeding the original seed into the random generator  // Idea here is to create an array of seeds by feeding the original seed into the random generator
883  // The code at the beginning of the function to compute the seed if one is given is  // The code at the beginning of the function to compute the seed if one is given is
884  // just supposed to introduce some variety (and ensure that multiple ranks don't get the same seed).  // just supposed to introduce some variety (and ensure that multiple ranks don't get the same seed).
# Line 873  void DataExpanded::randomFill(long seed) Line 906  void DataExpanded::randomFill(long seed)
906      // now we need to consider MPI since we don't want each rank to start with the same seed      // now we need to consider MPI since we don't want each rank to start with the same seed
907      seed+=getFunctionSpace().getDomain()->getMPIRank()*getFunctionSpace().getDomain()->getMPISize()*3;      seed+=getFunctionSpace().getDomain()->getMPIRank()*getFunctionSpace().getDomain()->getMPISize()*3;
908      prevseed=seed;      prevseed=seed;
909    
910    #ifdef MKLRANDOM
911    
912  #ifdef _OPENMP  #ifdef _OPENMP
913      int numthreads=omp_get_max_threads();      int numthreads=omp_get_max_threads();
914  #else  #else
915      int numthreads=1;      int numthreads=1;
916  #endif  #endif
   
 #ifdef MKLRANDOM  
917      double* seeds=new double[numthreads];      double* seeds=new double[numthreads];
918      VSLStreamStatePtr sstream;      VSLStreamStatePtr sstream;
919    
# Line 906  void DataExpanded::randomFill(long seed) Line 940  void DataExpanded::randomFill(long seed)
940      }      }
941      delete[] seeds;      delete[] seeds;
942  #else  #else
943      srand(seed);      boost::mt19937::result_type RMAX=base.max();
944      unsigned* seeds=new unsigned[numthreads];      seedGens(seed);
     for (int i=0;i<numthreads;++i)  
     {  
     seeds[i]=rand();  
     }  
945      DataVector&  dv=getVectorRW();      DataVector&  dv=getVectorRW();
946      long i;      long i;
947      const size_t dvsize=dv.size();      const size_t dvsize=dv.size();
948        
949      #pragma omp parallel private(i)      #pragma omp parallel private(i)
950      {      {
951      int tnum=0;      int tnum=0;
952      #ifdef _OPENMP      #ifdef _OPENMP
953      tnum=omp_get_thread_num();      tnum=omp_get_thread_num();
954      #endif      #endif
955      unsigned info=seeds[tnum];      boost::mt19937& generator=gens[tnum];
956            
957          #pragma omp for schedule(static)          #pragma omp for schedule(static)
958          for (i=0;i<dvsize;++i)          for (i=0;i<dvsize;++i)
959          {          {
960          
961          
962        
963          
964          
965  #ifdef _WIN32  #ifdef _WIN32
966          dv[i]=((double)rand())/RAND_MAX;          dv[i]=((double)generator())/RMAX;
967  #else  #else
968          dv[i]=((double)rand_r(&info))/RAND_MAX;          dv[i]=((double)generator())/RMAX;
969  #endif  #endif
970          }          }
971      }      }
     delete[] seeds;  
972  #endif  #endif
973  }  }
974    

Legend:
Removed from v.3549  
changed lines
  Added in v.3550

  ViewVC Help
Powered by ViewVC 1.1.26