easy-iso

commit 0fd21eba6944988e88b4fc8e9a6801b6e7ddf768
parent 6ced3abfff5224d46f8359a6fdd293eeff936b22
Author: Fabian Wermelinger <fabianw@mavt.ethz.ch>
Date:   Wed, 30 Nov 2016 20:39:46 +0100

added MPI wavelet reader

Diffstat:
M
src/Interpolator.h
 | 
33
+++++++++++++++++----------------
M
src/InterpolatorMPI.h
 | 
57
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++

2 files changed, 74 insertions(+), 16 deletions(-)
diff --git a/src/Interpolator.h b/src/Interpolator.h
@@ -217,8 +217,7 @@ void Interpolator<TKernel>::_load_wavelet(ArgumentParser& p)
         abort();
     }
 
-    const string fname(filename);
-    Reader_WaveletCompression myreader(fname, byte_swap, wavelet_type);
+    Reader_WaveletCompression myreader(filename, byte_swap, wavelet_type);
     myreader.load_file();
     const int NBX = myreader.xblocks();
     const int NBY = myreader.yblocks();
@@ -228,22 +227,24 @@ void Interpolator<TKernel>::_load_wavelet(ArgumentParser& p)
     Real* const data = new Real[num_elem];
 
     Real blockdata[_BLOCKSIZE_][_BLOCKSIZE_][_BLOCKSIZE_];
-    for (int iz=0; iz < NBZ; ++iz)
-        for (int iy=0; iy < NBY; ++iy)
-            for (int ix=0; ix < NBX; ++ix)
+    const int nBlocks = NBX * NBY * NBZ;
+    for (int i = 0; i < nBlocks; ++i)
+    {
+        const int ix = i%NBX;
+        const int iy = (i/NBX)%NBY;
+        const int iz = (i/(NBX*NBY))%NBZ;
+        const double zratio = myreader.load_block2(ix, iy, iz, blockdata);
+
+        for (int z=0; z < _BLOCKSIZE_; ++z)
+            for (int y=0; y < _BLOCKSIZE_; ++y)
             {
-                const double zratio = myreader.load_block2(ix, iy, iz, blockdata);
-
-                for (int z=0; z < _BLOCKSIZE_; ++z)
-                    for (int y=0; y < _BLOCKSIZE_; ++y)
-                    {
-                        assert(iy*_BLOCKSIZE_+y < maxDim[1]);
-                        assert(iz*_BLOCKSIZE_+z < maxDim[2]);
-                        const size_t offset = _BLOCKSIZE_*(static_cast<size_t>(ix) + NBX*(y+iy*_BLOCKSIZE_ + (z+static_cast<size_t>(iz)*_BLOCKSIZE_)*NBY*_BLOCKSIZE_));
-                        Real* const dst = data + offset;
-                        memcpy(dst, &blockdata[z][y][0], _BLOCKSIZE_*sizeof(Real));
-                    }
+                assert(iy*_BLOCKSIZE_+y < maxDim[1]);
+                assert(iz*_BLOCKSIZE_+z < maxDim[2]);
+                const size_t offset = _BLOCKSIZE_*(static_cast<size_t>(ix) + NBX*(y+iy*_BLOCKSIZE_ + (z+static_cast<size_t>(iz)*_BLOCKSIZE_)*NBY*_BLOCKSIZE_));
+                Real* const dst = data + offset;
+                memcpy(dst, &blockdata[z][y][0], _BLOCKSIZE_*sizeof(Real));
             }
+    }
 
     this->m_data = std::move(Matrix_t(maxDim[0], maxDim[1], maxDim[2], data));
     delete [] data;
diff --git a/src/InterpolatorMPI.h b/src/InterpolatorMPI.h
@@ -324,6 +324,63 @@ void InterpolatorMPI<Tinterp>::_load_h5_MPI(ArgumentParser& parser)
 template <typename Tinterp>
 void InterpolatorMPI<Tinterp>::_load_wavelet_MPI(ArgumentParser& p)
 {
+#ifdef _USE_CUBISMZ_
+    const bool byte_swap   = p("-swap").asBool(false);
+    const int wavelet_type = p("-wtype").asInt(1);
+    const std::string filename = p("file").asString("");
+    if (filename == "")
+    {
+        if (m_isroot)
+            std::cerr << "ERROR: InterpolatorMPI: -file is not specified. No input file given" << std::endl;
+        abort();
+    }
+
+    // Reader_WaveletCompressionMPI myreader(m_comm, filename, byte_swap, wavelet_type);
+    Reader_WaveletCompressionMPI myreader(MPI::COMM_WORLD, filename, byte_swap, wavelet_type);
+    myreader.load_file();
+    const int NBX = myreader.xblocks();
+    const int NBY = myreader.yblocks();
+    const int NBZ = myreader.zblocks();
+    assert(NBX % m_PESize[0] == 0);
+    assert(NBY % m_PESize[1] == 0);
+    assert(NBZ % m_PESize[2] == 0);
+    const int myNBX = NBX/m_PESize[0];
+    const int myNBY = NBY/m_PESize[1];
+    const int myNBZ = NBZ/m_PESize[2];
+
+    const int maxDim[3] = {myNBX*_BLOCKSIZE_, myNBY*_BLOCKSIZE_, myNBZ*_BLOCKSIZE_};
+    const size_t num_elem = static_cast<size_t>(maxDim[0]) * maxDim[1] * maxDim[2];
+    Real* const data = new Real[num_elem];
+
+    Real blockdata[_BLOCKSIZE_][_BLOCKSIZE_][_BLOCKSIZE_];
+    const int nBlocks = myNBX * myNBY * myNBZ;
+    for (int i = 0; i < nBlocks; ++i)
+    {
+        const int ix = i%myNBX;
+        const int iy = (i/myNBX)%myNBY;
+        const int iz = (i/(myNBX*myNBY))%myNBZ;
+        const double zratio = myreader.load_block2(
+                ix+myNBX*m_myPEIndex[0],
+                iy+myNBY*m_myPEIndex[1],
+                iz+myNBZ*m_myPEIndex[2],
+                blockdata);
+
+        for (int z=0; z < _BLOCKSIZE_; ++z)
+            for (int y=0; y < _BLOCKSIZE_; ++y)
+            {
+                assert(iy*_BLOCKSIZE_+y < maxDim[1]);
+                assert(iz*_BLOCKSIZE_+z < maxDim[2]);
+                const size_t offset = _BLOCKSIZE_*(static_cast<size_t>(ix) + myNBX*(y+iy*_BLOCKSIZE_ + (z+static_cast<size_t>(iz)*_BLOCKSIZE_)*myNBY*_BLOCKSIZE_));
+                Real* const dst = data + offset;
+                memcpy(dst, &blockdata[z][y][0], _BLOCKSIZE_*sizeof(Real));
+            }
+    }
+
+    this->m_data = std::move(Matrix_t(maxDim[0], maxDim[1], maxDim[2], data));
+    delete [] data;
+#else
+    fprintf(stderr, "WARNING: Executable was compiled without wavelet compressor support...\n");
+#endif /* _USE_CUBISMZ_ */
 }