/[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 2792 by jfenwick, Tue Dec 1 05:02:18 2009 UTC revision 3981 by jfenwick, Fri Sep 21 02:47:54 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2009 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16    
17  #include "Data.h"  #include "Data.h"
# Line 17  Line 19 
19  #include "DataException.h"  #include "DataException.h"
20  #include "DataConstant.h"  #include "DataConstant.h"
21  #include "DataTagged.h"  #include "DataTagged.h"
22    #include <limits>
23    
24    #include "esysUtils/Esys_MPI.h"
25    
26  #ifdef USE_NETCDF  #ifdef USE_NETCDF
27  #include <netcdfcpp.h>  #include <netcdfcpp.h>
28  #endif  #endif
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #endif  
29    
30  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
31  #include "DataMaths.h"  #include "DataMaths.h"
32    
33    //#define MKLRANDOM
34    
35    #ifdef MKLRANDOM
36    #include <mkl_vsl.h>
37    #else
38    #include <boost/random/mersenne_twister.hpp>
39    #endif
40    
41  using namespace std;  using namespace std;
42  using namespace boost::python;  using namespace boost::python;
43  using namespace boost;  using namespace boost;
# Line 161  DataExpanded::DataExpanded(const Functio Line 172  DataExpanded::DataExpanded(const Functio
172    
173  }  }
174    
175    DataExpanded::DataExpanded(const FunctionSpace& what,
176                               const DataTypes::ShapeType &shape,
177                               const double v)
178      : parent(what,shape)
179    {
180         ValueType& vec=m_data.getData();
181         //
182         // create the view of the data
183         initialise(what.getNumSamples(),what.getNumDPPSample());
184         // now we copy this value to all elements
185         const int L=getLength();
186         int i;
187         #pragma omp parallel for schedule(static) private(i)
188         for (i=0;i<L;++i)
189         {
190            vec[i]=v;
191         }
192    }
193    
194    
195    
196  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
197  {  {
198  }  }
# Line 663  DataExpanded::dump(const std::string fil Line 695  DataExpanded::dump(const std::string fil
695     const DataTypes::ShapeType& shape = getShape();     const DataTypes::ShapeType& shape = getShape();
696     int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();     int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
697     int mpi_num=getFunctionSpace().getDomain()->getMPISize();     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
698  #ifdef PASO_MPI  #ifdef ESYS_MPI
699     MPI_Status status;     MPI_Status status;
700  #endif  #endif
701    
702  #ifdef PASO_MPI  #ifdef ESYS_MPI
703     /* Serialize NetCDF I/O */     /* Serialize NetCDF I/O */
704     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);
705  #endif  #endif
# Line 727  DataExpanded::dump(const std::string fil Line 759  DataExpanded::dump(const std::string fil
759       if (! (var->put(d_ptr,dims)) )       if (! (var->put(d_ptr,dims)) )
760          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
761     }     }
762  #ifdef PASO_MPI  #ifdef ESYS_MPI
763     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);
764  #endif  #endif
765     #else     #else
# Line 817  DataExpanded::getVectorRO() const Line 849  DataExpanded::getVectorRO() const
849  }  }
850    
851    
852    #ifndef MKLRANDOM
853    namespace {
854        
855    boost::mt19937 base;        // used to seed all the other generators  
856    vector<boost::mt19937> gens;
857    
858    
859    void seedGens(long seed)
860    {
861    #ifdef _OPENMP
862        int numthreads=omp_get_max_threads();
863    #else
864        int numthreads=1;
865    #endif
866        if (gens.size()==0)     // we haven't instantiated the generators yet  
867        {
868            gens.resize(numthreads);    // yes this means all the generators will be owned by one thread in a NUMA sense      
869        }                   // I don't think it is worth a more complicated solution at this point
870        if (seed!=0)
871        {
872           base.seed((uint32_t)seed);   // without this cast, icc gets confused    
873           for (int i=0;i<numthreads;++i)
874           {
875            uint32_t b=base();
876                gens[i].seed(b);    // initialise each generator with successive random values      
877           }      
878        }
879    }
880      
881      
882    }
883    #endif
884    
885    // Idea here is to create an array of seeds by feeding the original seed into the random generator
886    // The code at the beginning of the function to compute the seed if one is given is
887    // just supposed to introduce some variety (and ensure that multiple ranks don't get the same seed).
888    // I make no claim about how well these initial seeds are distributed
889    void DataExpanded::randomFill(long seed)
890    {
891        CHECK_FOR_EX_WRITE
892        static unsigned prevseed=0; // So if we create a bunch of objects we don't get the same start seed
893        if (seed==0)        // for each one
894        {
895        if (prevseed==0)
896        {
897            time_t s=time(0);
898            seed=s;
899        }
900        else
901        {
902            seed=prevseed+419;  // these numbers are arbitrary
903            if (seed>3040101)       // I want to avoid overflow on 32bit systems
904            {
905            seed=((int)(seed)%0xABCD)+1;
906            }
907        }
908        }
909        // now we need to consider MPI since we don't want each rank to start with the same seed
910        seed+=getFunctionSpace().getDomain()->getMPIRank()*getFunctionSpace().getDomain()->getMPISize()*3;
911        prevseed=seed;
912    
913    #ifdef MKLRANDOM
914    
915    #ifdef _OPENMP
916        int numthreads=omp_get_max_threads();
917    #else
918        int numthreads=1;
919    #endif
920        double* seeds=new double[numthreads];
921        VSLStreamStatePtr sstream;
922    
923        int status=vslNewStream(&sstream, VSL_BRNG_MT19937, seed);  // use a Mersenne Twister
924        numeric_limits<double> dlim;
925        vdRngUniform(VSL_METHOD_DUNIFORM_STD, sstream , numthreads, seeds, -1, 1);
926        vslDeleteStream(&sstream);
927        DataVector& dv=getVectorRW();
928        size_t dvsize=dv.size();
929        #pragma omp parallel
930        {
931        int tnum=0;
932        #ifdef _OPENMP
933        tnum=omp_get_thread_num();
934        #endif
935        VSLStreamStatePtr stream;
936        // the 12345 is a hack to give us a better chance of getting different integer seeds.
937            int status=vslNewStream(&stream, VSL_BRNG_MT19937, seeds[tnum]*12345);  // use a Mersenne Twister
938        int bigchunk=(dvsize/numthreads+1);
939        int smallchunk=dvsize-bigchunk*(numthreads-1);
940        int chunksize=(tnum<(numthreads-1))?bigchunk:smallchunk;
941            vdRngUniform(VSL_METHOD_DUNIFORM_STD, stream, chunksize, &(dv[bigchunk*tnum]), 0,1);
942            vslDeleteStream(&stream);
943        }
944        delete[] seeds;
945    #else
946        boost::mt19937::result_type RMAX=base.max();
947        seedGens(seed);
948        DataVector&  dv=getVectorRW();
949        long i;
950        const size_t dvsize=dv.size();
951        
952        #pragma omp parallel private(i)
953        {
954        int tnum=0;
955        #ifdef _OPENMP
956        tnum=omp_get_thread_num();
957        #endif
958        boost::mt19937& generator=gens[tnum];
959        
960            #pragma omp for schedule(static)
961            for (i=0;i<dvsize;++i)
962            {
963          
964          
965        
966          
967          
968    #ifdef _WIN32
969            dv[i]=((double)generator())/RMAX;
970    #else
971            dv[i]=((double)generator())/RMAX;
972    #endif
973            }
974        }
975    #endif
976    }
977    
978  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.2792  
changed lines
  Added in v.3981

  ViewVC Help
Powered by ViewVC 1.1.26