Commit dd9742df authored by BRAMAS Berenger's avatar BRAMAS Berenger
Browse files

make the hmat demo working

parent 18666b75
......@@ -22,45 +22,60 @@
#include "Utils/FBlas.hpp"
template <class FReal>
class FDenseMatrix
{
class FDenseMatrix{
protected:
// members
FReal* block;
const int height;
const int width;
const int leading_dimension;
const int level;
int nbRows;
int nbCols;
int level;
FDenseMatrix(const FDenseMatrix&) = delete;
FDenseMatrix& operator=(const FDenseMatrix&) = delete;
public:
FDenseMatrix()
: block(nullptr), nbRows(0), nbCols(0), level(0) {
}
// ctor
explicit FDenseMatrix(const FReal* in_block, const int in_height, const int in_width, const int in_leading_dimension /*h<N*/, const int in_level)
: height(in_height), width(in_width), leading_dimension(in_leading_dimension), level(in_level)
{
template <class ViewerClass>
void fill(const ViewerClass& viewer, const int inLevel){
clear();
// Allocate memory
block = new FReal[height*width];
// Copy
FBlas::copy(height*width,in_block,block);
level = inLevel;
nbRows = viewer.getNbRows();
nbCols = viewer.getNbCols();
block = new FReal[nbRows*nbCols];
for(int idxRow = 0 ; idxRow < nbRows ; ++idxRow){
for(int idxCol = 0 ; idxCol < nbCols ; ++idxCol){
block[idxCol*nbRows+idxRow] = viewer.getValue(idxRow,idxCol);
}
}
};
// dtor
~FDenseMatrix(){
// Free memory
delete[] block;
clear();
};
void MV(FReal res[], const FReal vec[], const FReal scale = FReal(1.)) {
FBlas::gemv(height, width, scale, block, vec, res);
void clear(){
nbRows = 0;
nbCols = 0;
nbCols = 0;
delete[] block;
block = 0;
}
void gemv(FReal res[], const FReal vec[], const FReal scale = FReal(1.)) const {
FBlas::gemva(nbRows, nbCols, scale, const_cast<FReal*>(block), const_cast<FReal*>(vec), res);
}
void gemm(FReal res[], const FReal vec[], const int nbrhs, const FReal scale = FReal(1.)) const {
FBlas::gemva(nbRows, nbCols, scale, const_cast<FReal*>(block), const_cast<FReal*>(vec), res);
}
};
#endif // FDENSEMATRIX_HPP
......
......@@ -136,21 +136,23 @@ public:
void fillBlocks(MatrixClass& matrix){
for(int idxLevel = 1 ; idxLevel < height-1 ; ++idxLevel){
for(int idxCell = 0 ; idxCell < nbCells[idxLevel] ; ++idxCell){
cells[idxLevel][idxCell].fill(matrix.getBlock(
cells[idxLevel][idxCell].cell.fill(matrix.getBlock(
cells[idxLevel][idxCell].infos.row,
cells[idxLevel][idxCell].infos.col,
cells[idxLevel][idxCell].infos.nbRows,
cells[idxLevel][idxCell].infos.nbCols
));
),
cells[idxLevel][idxCell].infos.level);
}
}
for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){
leaves[idxLeaf].fill(matrix.getBlock(
leaves[idxLeaf].leaf.fill(matrix.getBlock(
leaves[idxLeaf].infos.row,
leaves[idxLeaf].infos.col,
leaves[idxLeaf].infos.nbRows,
leaves[idxLeaf].infos.nbCols
));
),
leaves[idxLeaf].infos.level);
}
}
......@@ -158,47 +160,49 @@ public:
void fillBlocks(const MatrixClass& matrix){
for(int idxLevel = 1 ; idxLevel < height-1 ; ++idxLevel){
for(int idxCell = 0 ; idxCell < nbCells[idxLevel] ; ++idxCell){
cells[idxLevel][idxCell].fill(matrix.getBlock(
cells[idxLevel][idxCell].cell.fill(matrix.getBlock(
cells[idxLevel][idxCell].infos.row,
cells[idxLevel][idxCell].infos.col,
cells[idxLevel][idxCell].infos.nbRows,
cells[idxLevel][idxCell].infos.nbCols
));
),
cells[idxLevel][idxCell].infos.level);
}
}
for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){
leaves[idxLeaf].fill(matrix.getBlock(
leaves[idxLeaf].leaf.fill(matrix.getBlock(
leaves[idxLeaf].infos.row,
leaves[idxLeaf].infos.col,
leaves[idxLeaf].infos.nbRows,
leaves[idxLeaf].infos.nbCols
));
),
leaves[idxLeaf].infos.level);
}
}
void gemv(FReal res[], const FReal vec[]) const {
for(int idxLevel = 1 ; idxLevel < height-1 ; ++idxLevel){
for(int idxCell = 0 ; idxCell < nbCells[idxLevel] ; ++idxCell){
cells[idxLevel][idxCell].gemv(&res[cells[idxLevel][idxCell].infos.col],
&vec[cells[idxLevel][idxCell].infos.row]);
cells[idxLevel][idxCell].cell.gemv(&res[cells[idxLevel][idxCell].infos.row],
&vec[cells[idxLevel][idxCell].infos.col]);
}
}
for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){
leaves[idxLeaf].gemv(&res[leaves[idxLeaf].infos.col],
&vec[leaves[idxLeaf].infos.row]);
leaves[idxLeaf].leaf.gemv(&res[leaves[idxLeaf].infos.row],
&vec[leaves[idxLeaf].infos.col]);
}
}
void gemm(FReal res[], const FReal mat[], const int nbRhs) const {
for(int idxLevel = 1 ; idxLevel < height-1 ; ++idxLevel){
for(int idxCell = 0 ; idxCell < nbCells[idxLevel] ; ++idxCell){
cells[idxLevel][idxCell].gemm(&res[cells[idxLevel][idxCell].infos.col],
cells[idxLevel][idxCell].cell.gemm(&res[cells[idxLevel][idxCell].infos.col],
&mat[cells[idxLevel][idxCell].infos.row],
nbRhs, dim);
}
}
for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){
leaves[idxLeaf].gemm(&res[leaves[idxLeaf].infos.col],
leaves[idxLeaf].leaf.gemm(&res[leaves[idxLeaf].infos.col],
&mat[leaves[idxLeaf].infos.row],
nbRhs, dim);
}
......
#ifndef FMATGRID_HPP
#define FMATGRID_HPP
#ifndef FMATDENSE_HPP
#define FMATDENSE_HPP
// @SCALFMM_PRIVATE
......@@ -7,40 +7,58 @@
#include "Utils/FAssert.hpp"
#include "FDenseBlockWrapper.hpp"
#include "../Utils/FHUtils.hpp"
template <class FReal>
class FMatGrid {
class FMatDense {
protected:
int matDim;
FReal* values;
FMatGrid(const FMatGrid&) = delete;
FMatGrid& operator=(const FMatGrid&) = delete;
FMatDense(const FMatDense&) = delete;
FMatDense& operator=(const FMatDense&) = delete;
public:
using BlockDescriptor = FDenseBlockWrapper<FReal>;
FMatGrid(const int inDim, const FReal* inValues)
explicit FMatDense(const int inDim, const FReal* inValues = nullptr)
: matDim(inDim), values(nullptr){
values = new FReal[matDim*matDim];
for(int idxVal = 0 ; idxVal < matDim*matDim ; ++idxVal){
values[idxVal] = inValues[idxVal];
if(inValues){
for(int idxVal = 0 ; idxVal < matDim*matDim ; ++idxVal){
values[idxVal] = inValues[idxVal];
}
}
else{
FSetToZeros(values, matDim*matDim);
}
}
~FMatGrid(){
~FMatDense(){
delete[] values;
}
const FReal& getVal(const int idxRow , const int idxCol) const{
return values[idxCol*matDim+idxRow];
}
FReal& getVal(const int idxRow , const int idxCol) {
return values[idxCol*matDim+idxRow];
}
void setVal(const int idxRow , const int idxCol, const FReal& val) {
values[idxCol*matDim+idxRow] = val;
}
FDenseBlockWrapper<FReal> getBlock(const int rowIdx, const int colIdx, const int nbRows, const int nbCols) const {
// static_assert(std::is_move_constructible<BlockClass>::value, "The block class must be movable");
// static_assert(std::is_move_assignable<BlockClass>::value, "The block class must be movable");
FAssertLF(0 < rowIdx);
FAssertLF(0 < colIdx);
FAssertLF(0 < nbRows);
FAssertLF(0 < nbCols);
FAssertLF(rowIdx + nbRows <= matDim);
FAssertLF(colIdx + nbRows <= matDim);
return FDenseBlockWrapper<FReal>(&values[colIdx*matDim+rowIdx], nbRows, nbCols, matDim);
}
};
#endif // FMATGRID_HPP
#endif // FMATDENSE_HPP
......@@ -2,8 +2,9 @@
// @SCALFMM_PRIVATE
#include "../Src/Containers/FDiv2Bissection.hpp"
#include "../Src/Viewers/FMatGrid.hpp"
#include "../Src/Utils/FSvgRect.hpp"
#include "../Src/Viewers/FDenseBlockWrapper.hpp"
#include "../Src/Blocks/FDenseMatrix.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FParameterNames.hpp"
......@@ -35,12 +36,11 @@ int main(int argc, char** argv){
std::cout << "Config : outputfile = " << outputfile << "\n";
typedef double FReal;
typedef int LeafClass;
typedef int CellClass;
typedef FDenseMatrix<FReal> LeafClass;
typedef FDenseMatrix<FReal> CellClass;
typedef FDiv2Bissection<FReal, LeafClass, CellClass> GridClass;
GridClass bissection(dim, height);
{
FSvgRect output(outputfile, dim);
......
......@@ -2,12 +2,14 @@
// @SCALFMM_PRIVATE
#include "../Src/Containers/FDiv2Bissection.hpp"
#include "../Src/Viewers/FMatGrid.hpp"
#include "../Src/Viewers/FMatDense.hpp"
#include "../Src/Blocks/FDenseMatrix.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FParameterNames.hpp"
#include <memory>
int main(int argc, char** argv){
static const FParameterNames DimParam = {
{"-N", "-nb", "-dim"} ,
......@@ -18,9 +20,7 @@ int main(int argc, char** argv){
"Number of dissection (+1)."
};
FHelpDescribeAndExit(argc, argv,
"Test the bisection.",SvgOutParam,
DimParam,HeightParam);
FHelpDescribeAndExit(argc, argv, "Test the bisection.", DimParam,HeightParam);
const int dim = FParameters::getValue(argc, argv, DimParam.options, 100);
const int height = FParameters::getValue(argc, argv, HeightParam.options, 4);
......@@ -28,15 +28,48 @@ int main(int argc, char** argv){
std::cout << "Config : dim = " << dim << "\n";
std::cout << "Config : height = " << height << "\n";
typedef double FReal;
typedef FMatDense<FReal> MatrixClass;
MatrixClass matrix(dim);
for(int idxRow = 0; idxRow < dim ; ++idxRow){
for(int idxCol = 0; idxCol < dim ; ++idxCol){
matrix.setVal(idxRow, idxCol, 1./(FMath::Abs(FReal(idxRow-idxCol))+1.));
}
}
std::unique_ptr<FReal[]> vec(new FReal[dim]);
for(int idxVal = 0 ; idxVal < dim ; ++idxVal){
vec[idxVal] = 1.0;
}
std::unique_ptr<FReal[]> resTest(new FReal[dim]);
FSetToZeros(resTest.get(), dim);
{
for(int idxRow = 0; idxRow < dim ; ++idxRow){
for(int idxCol = 0; idxCol < dim ; ++idxCol){
resTest[idxRow] += vec[idxCol] * matrix.getVal(idxRow, idxCol);
}
}
}
{
typedef double FReal;
typedef FMatGrid<FReal> MatrixClass;
typedef FDenseMatrix<FReal> LeafClass;
typedef FDenseMatrix<FReal> CellClass;
typedef FDiv2Bissection<FReal, LeafClass, CellClass> GridClass;
GridClass bissection(dim, height);
bissection.fillBlocks(matrix);
GridClass grid(dim, height);
grid.fillBlocks(matrix);
std::unique_ptr<FReal[]> resDense(new FReal[dim]);
FSetToZeros(resDense.get(), dim);
grid.gemv(resDense.get(), vec.get());
FMath::FAccurater<FReal> testDense(resTest.get(), resDense.get(), dim);
std::cout << "Test Dense, Error = " << testDense << "\n";
}
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