Merge branch 'master' of git+ssh://scm.gforge.inria.fr//gitroot//scalfmm/scalfmm

Addons/HMat/CMakeLists.txt

@@ -12,6 +12,8 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${SCALFMM_CXX_FLAGS}")
 ENABLE_LANGUAGE(CXX C)
 MESSAGE(STATUS " CXX  ${CMAKE_CXX_COMPILER_ID}" )
 
+# Options
+OPTION( SCALFMM_ADDON_HMAT "Set to ON to build ScalFMM HMat interface" OFF )
 
 # if ask to build addon
 if(SCALFMM_ADDON_HMAT)

Addons/HMat/Src/Utils/FMatrixIO.hpp

+// ===================================================================================
+// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
+// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
+// This software is a computer program whose purpose is to compute the FMM.
+//
+// This software is governed by the CeCILL-C and LGPL licenses and
+// abiding by the rules of distribution of free software.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public and CeCILL-C Licenses for more details.
+// "http://www.cecill.info".
+// "http://www.gnu.org/licenses".
+// ===================================================================================
+// 
+// @SCALFMM_PRIVATE
+// 
+#ifndef FMATRIXIO_HPP
+#define FMATRIXIO_HPP
+
+// std includes
+#include <numeric>
+#include <stdexcept>
+#include <string>
+#include <sstream>
+#include <fstream>
+#include <typeinfo>
+
+// ScalFMM includes
+#include "Utils/FGlobal.hpp"
+#include "Utils/FParameters.hpp"
+#include "Utils/FMath.hpp"
+
+/**
+ * @brief The FMatrixIO class
+ * Beware! Only handles sizes (nrows/ncols) representable by int32. 
+ */
+class FMatrixIO {
+
+public:
+
+    // Write matrix to binary file.
+    template<class FReal>
+    static void write(const FSize nrows, const FSize ncols, FReal* matrix, const std::string filename){
+
+        // store into binary file
+        std::ofstream stream(filename.c_str(),
+                             std::ios::out | std::ios::binary | std::ios::trunc);
+        if (stream.good()) {
+            stream.seekp(0);
+            // 1) write number of rows (int)
+            //is_int(nrows);
+            int _nrows = int(nrows);
+            stream.write(reinterpret_cast<char*>(&_nrows), sizeof(int));
+            // 2) write number of cols (int)
+            //is_int(ncols);
+            int _ncols = int(ncols);
+            stream.write(reinterpret_cast<char*>(&_ncols), sizeof(int));
+            // 1) write matrix
+            stream.write(reinterpret_cast<char*>(matrix), sizeof(FReal)*nrows*ncols);
+        }   
+        else throw std::runtime_error("File could not be opened to write");
+        stream.close();
+
+    }
+
+    // Read matrix from binary file. 
+    // Also allocates memory! Please set matrix to nullptr before calling this function. 
+    // nrows is just here to control the values that is read from file
+    // ncols  is updated with the value read in the file (because it is not always known!)
+    template<class FReal>
+    static void read(const FSize nrows, FSize &ncols, FReal* &matrix, const std::string filename, const bool readSize = true){
+
+        // start reading process
+        if (matrix) throw std::runtime_error("Matrix is already set!");
+
+        std::ifstream stream(filename.c_str(),
+                             std::ios::in | std::ios::binary | std::ios::ate);
+        const std::ifstream::pos_type size = stream.tellg();
+        if (size<=0) {
+            std::cout << "Info: The requested binary file " << filename
+                      << " does not yet exist. Compute it now ... " << std::endl;
+            return;
+        } 
+        if (stream.good()) {
+            stream.seekg(0);
+            if(readSize){
+                // 1) read number of rows (int)
+                int _nrows;
+                stream.read(reinterpret_cast<char*>(&_nrows), sizeof(int));
+                //is_int(nrows);
+                if (_nrows!=int(nrows)) throw std::runtime_error("Read nrows and input nrows do not correspond");
+                // 2) read number of cols (int)
+                int _ncols;
+                stream.read(reinterpret_cast<char*>(&_ncols), sizeof(int));
+                //is_int(ncols);
+                ncols=_ncols;
+            }
+            // 3) read matrix
+            matrix = new FReal[nrows*ncols];
+            stream.read(reinterpret_cast<char*>(matrix), sizeof(FReal)*nrows*ncols);
+        }   else throw std::runtime_error("File could not be opened to read");
+        stream.close();
+
+    }
+
+};
+
+
+#endif // FMATRIXIO_HPP

Addons/HMat/Tests/testGenerateKernelMatrix.cpp

+
+// @SCALFMM_PRIVATE
+
+#include <iostream>
+#include <fstream>
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <random> // for normal distribution generator
+#include <algorithm> // for sort
+
+
+// ScalFMM includes
+#include "Files/FFmaGenericLoader.hpp"
+#include "Utils/FGlobal.hpp"
+#include "Utils/FTic.hpp"
+#include "Utils/FParameters.hpp"
+#include "Utils/FMemUtils.hpp"
+#include "Utils/FBlas.hpp" // for FBlas::potrf (and QR,SVD...)
+#include "Utils/FMath.hpp"
+#include "Utils/FParameterNames.hpp"
+
+#include "Kernels/Interpolation/FInterpMatrixKernel.hpp" // for kernel matrices
+
+// not mandatory but useful to define some flags
+#include "Core/FFmmAlgorithm.hpp"
+#include "Core/FFmmAlgorithmThread.hpp"
+
+// ScalFMM/HMat includes
+#include "../Src/Utils/FMatrixIO.hpp"
+
+/**
+* In this file we explicitely build a kernel matrix and then store it in a binary file. 
+*
+* Author: Pierre Blanchard (pierre.blanchard@inria.fr)
+* Date created: November 26th, 2015 
+*/
+
+
+int main(int argc, char* argv[])
+{
+    std::cout << "Addons/HMat: Build a kernel matrix and then store it in a binary file." << std::endl;
+
+    // Example for unitCube5000
+    // ./Addons/HMat/testGenerateKernelMatrix -binin -N 5000 -d unitCube5000
+
+    ////////////////////////////////////////////////////////////////////
+    /// Timers 
+    FTic time; 
+    time.tic();
+
+    ////////////////////////////////////////////////////////////////////
+    /// Parameters
+    // Verbose (print matrices)
+    const int verbose = FParameters::getValue(argc, argv, "-v", 0);
+    std::cout<< "Verbose level: " << verbose <<std::endl;
+    // size of the grid
+    const FSize matrixSize = FParameters::getValue(argc,argv,"-N", 1000); 
+    std::cout<< "Full matrix size: " << matrixSize <<std::endl;
+    // Precision
+    typedef double FReal;
+    if(sizeof(FReal)==4)
+        std::cout<< "Precision: Single Float" <<std::endl;
+    if(sizeof(FReal)==8)
+        std::cout<< "Precision: Double Float" <<std::endl;
+
+    // Memory
+    std::cout<< "Memory requirement: " << sizeof(FReal) * matrixSize * matrixSize *1e-6 << " MBytes" <<std::endl;
+
+    // Data path
+    const std::string ioPath = FParameters::getStr(argc,argv,"-path", std::string("../Addons/HMat/Data/").c_str());
+
+    // Read geometry
+    const std::string distributionName(FParameters::getStr(argc,argv,"-d",   "unitCube1000"));
+    std::string distributionFileName = ioPath + distributionName;
+    if(  FParameters::existParameter(argc, argv, FParameterDefinitions::InputBinFormat.options)){
+        distributionFileName += ".bfma";
+    }
+    else {
+        distributionFileName += ".fma";
+    }
+
+    // open particle file
+    FFmaGenericLoader<FReal> loader(distributionFileName.c_str());
+    if(!loader.isOpen()) throw std::runtime_error("Particle distribution file couldn't be opened!");
+
+    ////////////////////////////////////////////////////////////////////
+    /// Load 3D grid from distribution file
+    FPoint<FReal>* grid = new FPoint<FReal>[matrixSize];
+    FReal* pV = new FReal[matrixSize];
+
+    for(FSize idxPart = 0 ; idxPart < matrixSize ; ++idxPart){
+        FPoint<FReal> position;
+        FReal physicalValue = 0.0;
+        loader.fillParticle(&position,&physicalValue);
+        grid[idxPart]=position;
+        pV[idxPart]=physicalValue;
+    }
+
+
+
+    ////////////////////////////////////////////////////////////////////
+    /// Build kernel matrix K
+    
+    // Interaction kernel evaluator
+    typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
+    const MatrixKernelClass MatrixKernel;
+    const std::string MatrixKernelID = MatrixKernelClass::getID();//.c_str();
+
+    // Allocate memory
+    FReal* K = new FReal[matrixSize*matrixSize];
+
+    // Build kernel matrix
+    FTic timeAssK;
+    for(FSize idxRow = 0 ; idxRow < matrixSize  ; ++idxRow)
+        for(FSize idxCol = 0 ; idxCol < matrixSize  ; ++idxCol)
+            K[idxRow*matrixSize+idxCol] = MatrixKernel.evaluate(grid[idxRow],
+                                                                grid[idxCol]);
+    double tAssK = timeAssK.tacAndElapsed();
+    std::cout << "... took @tAssK = "<< tAssK <<"\n";
+
+    ////////////////////////////////////////////////////////////////////
+    /// Write kernel matrix in binary file
+
+    // Output file name
+    const std::string matrixName((distributionName + "_" + MatrixKernelID).c_str()); // if nothing specified then use file associated with n=50
+    const std::string fileName = ioPath + matrixName + ".bin";
+
+    // Write 
+    std::cout<< "Write matrix in binary file: " << fileName << "\n";
+    FTic timeWriteMat;
+    double tWriteMat;
+    timeWriteMat.tic();
+
+    FMatrixIO::write<FReal>(matrixSize,matrixSize,K,fileName);
+
+    tWriteMat = timeWriteMat.tacAndElapsed();
+    std::cout << "... took @tWriteMat = "<< tWriteMat <<"\n";
+
+    /// Free memory
+    delete[] K;
+
+    return 0;
+}