/[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 3544 by jfenwick, Wed Jul 27 06:19:10 2011 UTC revision 3911 by jfenwick, Thu Jun 14 01:01:03 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2010 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# 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((uint32_t)seed);   // without this cast, icc gets confused    
871           for (int i=0;i<numthreads;++i)
872           {
873            uint32_t b=base();
874                gens[i].seed(b);    // initialise each generator with successive random values      
875           }      
876        }
877    }
878      
879      
880    }
881    #endif
882    
883  // 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
884  // 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
885  // 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 907  void DataExpanded::randomFill(long seed)
907      // 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
908      seed+=getFunctionSpace().getDomain()->getMPIRank()*getFunctionSpace().getDomain()->getMPISize()*3;      seed+=getFunctionSpace().getDomain()->getMPIRank()*getFunctionSpace().getDomain()->getMPISize()*3;
909      prevseed=seed;      prevseed=seed;
910    
911    #ifdef MKLRANDOM
912    
913  #ifdef _OPENMP  #ifdef _OPENMP
914      int numthreads=omp_get_max_threads();      int numthreads=omp_get_max_threads();
915  #else  #else
916      int numthreads=1;      int numthreads=1;
917  #endif  #endif
   
 #ifdef MKLRANDOM  
918      double* seeds=new double[numthreads];      double* seeds=new double[numthreads];
919      VSLStreamStatePtr sstream;      VSLStreamStatePtr sstream;
920    
# Line 906  void DataExpanded::randomFill(long seed) Line 941  void DataExpanded::randomFill(long seed)
941      }      }
942      delete[] seeds;      delete[] seeds;
943  #else  #else
944      srand(seed);      boost::mt19937::result_type RMAX=base.max();
945      unsigned* seeds=new unsigned[numthreads];      seedGens(seed);
     for (int i=0;i<numthreads;++i)  
     {  
     seeds[i]=rand();  
     }  
946      DataVector&  dv=getVectorRW();      DataVector&  dv=getVectorRW();
947      long i;      long i;
948      const size_t dvsize=dv.size();      const size_t dvsize=dv.size();
949        
950      #pragma omp parallel private(i)      #pragma omp parallel private(i)
951      {      {
952      int tnum=0;      int tnum=0;
953      #ifdef _OPENMP      #ifdef _OPENMP
954      tnum=omp_get_thread_num();      tnum=omp_get_thread_num();
955      #endif      #endif
956      unsigned info=seeds[tnum];      boost::mt19937& generator=gens[tnum];
957            
958          #pragma omp for schedule(static)          #pragma omp for schedule(static)
959          for (i=0;i<dvsize;++i)          for (i=0;i<dvsize;++i)
960          {          {
961          
962          
963        
964          
965          
966  #ifdef _WIN32  #ifdef _WIN32
967          dv[i]=((double)rand())/RAND_MAX;          dv[i]=((double)generator())/RMAX;
968  #else  #else
969          dv[i]=((double)rand_r(&info))/RAND_MAX;          dv[i]=((double)generator())/RMAX;
970  #endif  #endif
971          }          }
972      }      }
     delete[] seeds;  
973  #endif  #endif
974  }  }
975    

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

  ViewVC Help
Powered by ViewVC 1.1.26