/[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 2881 by jfenwick, Thu Jan 28 02:03:15 2010 UTC revision 3552 by jfenwick, Mon Aug 22 00:05:17 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"
23    
24  #ifdef USE_NETCDF  #ifdef USE_NETCDF
25  #include <netcdfcpp.h>  #include <netcdfcpp.h>
26  #endif  #endif
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #endif  
27    
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    #else
36    #include <boost/random/mersenne_twister.hpp>
37    #endif
38    
39  using namespace std;  using namespace std;
40  using namespace boost::python;  using namespace boost::python;
41  using namespace boost;  using namespace boost;
# Line 161  DataExpanded::DataExpanded(const Functio Line 170  DataExpanded::DataExpanded(const Functio
170    
171  }  }
172    
173    DataExpanded::DataExpanded(const FunctionSpace& what,
174                               const DataTypes::ShapeType &shape,
175                               const double v)
176      : parent(what,shape)
177    {
178         ValueType& vec=m_data.getData();
179         //
180         // create the view of the data
181         initialise(what.getNumSamples(),what.getNumDPPSample());
182         // now we copy this value to all elements
183         const int L=getLength();
184         int i;
185         #pragma omp parallel for schedule(static) private(i)
186         for (i=0;i<L;++i)
187         {
188            vec[i]=v;
189         }
190    }
191    
192    
193    
194  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
195  {  {
196  }  }
# Line 663  DataExpanded::dump(const std::string fil Line 693  DataExpanded::dump(const std::string fil
693     const DataTypes::ShapeType& shape = getShape();     const DataTypes::ShapeType& shape = getShape();
694     int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();     int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
695     int mpi_num=getFunctionSpace().getDomain()->getMPISize();     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
696  #ifdef PASO_MPI  #ifdef ESYS_MPI
697     MPI_Status status;     MPI_Status status;
698  #endif  #endif
699    
700  #ifdef PASO_MPI  #ifdef ESYS_MPI
701     /* Serialize NetCDF I/O */     /* Serialize NetCDF I/O */
702     if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);     if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
703  #endif  #endif
# Line 727  DataExpanded::dump(const std::string fil Line 757  DataExpanded::dump(const std::string fil
757       if (! (var->put(d_ptr,dims)) )       if (! (var->put(d_ptr,dims)) )
758          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
759     }     }
760  #ifdef PASO_MPI  #ifdef ESYS_MPI
761     if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);     if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
762  #endif  #endif
763     #else     #else
# Line 817  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
884    // 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).
886    // I make no claim about how well these initial seeds are distributed
887    void DataExpanded::randomFill(long seed)
888    {
889        CHECK_FOR_EX_WRITE
890        static unsigned prevseed=0; // So if we create a bunch of objects we don't get the same start seed
891        if (seed==0)        // for each one
892        {
893        if (prevseed==0)
894        {
895            time_t s=time(0);
896            seed=s;
897        }
898        else
899        {
900            seed=prevseed+419;  // these numbers are arbitrary
901            if (seed>3040101)       // I want to avoid overflow on 32bit systems
902            {
903            seed=((int)(seed)%0xABCD)+1;
904            }
905        }
906        }
907        // 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;
909        prevseed=seed;
910    
911    #ifdef MKLRANDOM
912    
913    #ifdef _OPENMP
914        int numthreads=omp_get_max_threads();
915    #else
916        int numthreads=1;
917    #endif
918        double* seeds=new double[numthreads];
919        VSLStreamStatePtr sstream;
920    
921        int status=vslNewStream(&sstream, VSL_BRNG_MT19937, seed);  // use a Mersenne Twister
922        numeric_limits<double> dlim;
923        vdRngUniform(VSL_METHOD_DUNIFORM_STD, sstream , numthreads, seeds, -1, 1);
924        vslDeleteStream(&sstream);
925        DataVector& dv=getVectorRW();
926        size_t dvsize=dv.size();
927        #pragma omp parallel
928        {
929        int tnum=0;
930        #ifdef _OPENMP
931        tnum=omp_get_thread_num();
932        #endif
933        VSLStreamStatePtr stream;
934        // the 12345 is a hack to give us a better chance of getting different integer seeds.
935            int status=vslNewStream(&stream, VSL_BRNG_MT19937, seeds[tnum]*12345);  // use a Mersenne Twister
936        int bigchunk=(dvsize/numthreads+1);
937        int smallchunk=dvsize-bigchunk*(numthreads-1);
938        int chunksize=(tnum<(numthreads-1))?bigchunk:smallchunk;
939            vdRngUniform(VSL_METHOD_DUNIFORM_STD, stream, chunksize, &(dv[bigchunk*tnum]), 0,1);
940            vslDeleteStream(&stream);
941        }
942        delete[] seeds;
943    #else
944        boost::mt19937::result_type RMAX=base.max();
945        seedGens(seed);
946        DataVector&  dv=getVectorRW();
947        long i;
948        const size_t dvsize=dv.size();
949        
950        #pragma omp parallel private(i)
951        {
952        int tnum=0;
953        #ifdef _OPENMP
954        tnum=omp_get_thread_num();
955        #endif
956        boost::mt19937& generator=gens[tnum];
957        
958            #pragma omp for schedule(static)
959            for (i=0;i<dvsize;++i)
960            {
961          
962          
963        
964          
965          
966    #ifdef _WIN32
967            dv[i]=((double)generator())/RMAX;
968    #else
969            dv[i]=((double)generator())/RMAX;
970    #endif
971            }
972        }
973    #endif
974    }
975    
976  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.2881  
changed lines
  Added in v.3552

  ViewVC Help
Powered by ViewVC 1.1.26