Commit e9b18eae authored by BLANCHARD Pierre's avatar BLANCHARD Pierre
Browse files

HMat: Provided some spatial grids along with programs to generate associated...

HMat: Provided some spatial grids along with programs to generate associated kernel matrices; fixed CMakeLists.
parent f0a38873
......@@ -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)
......
// ===================================================================================
// 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
// @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;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment