From 96c6486d3a44d1e836b88d6edfb85fc533709228 Mon Sep 17 00:00:00 2001 From: Berenger Bramas <Berenger.Bramas@inria.fr> Date: Wed, 16 Dec 2015 08:56:02 +0100 Subject: [PATCH] update to build block matrix --- .../Src/Containers/FPartitionsMapping.hpp | 125 +++++------------ .../Src/Viewers/FDenseBlockPermWrapper.hpp | 17 +-- Addons/HMat/Src/Viewers/FMatDensePerm.hpp | 12 +- Addons/HMat/Tests/testPartitionsMapping.cpp | 3 +- .../HMat/Tests/testPartitionsMappingGemv.cpp | 3 +- .../Tests/testPartitionsMappingGemvBlock.cpp | 127 ++++++++++++++++++ 6 files changed, 175 insertions(+), 112 deletions(-) create mode 100644 Addons/HMat/Tests/testPartitionsMappingGemvBlock.cpp diff --git a/Addons/HMat/Src/Containers/FPartitionsMapping.hpp b/Addons/HMat/Src/Containers/FPartitionsMapping.hpp index 0ad41a1e1..f42a3e12d 100644 --- a/Addons/HMat/Src/Containers/FPartitionsMapping.hpp +++ b/Addons/HMat/Src/Containers/FPartitionsMapping.hpp @@ -29,105 +29,75 @@ #include <memory> -template <class FReal, class LeafClass, class CellClass > +template <class FReal, class CellClass > class FPartitionsMapping { protected: - struct LeafNode { - FBlockDescriptor infos; - LeafClass leaf; - }; - struct CellNode { FBlockDescriptor infos; CellClass cell; }; const int dim; + const int nbPartitions; + const int nbCells; - int nbCells; CellNode* cells; - int nbLeaves; - LeafNode* leaves; - - int totalNbBlocks; - FPartitionsMapping(const FPartitionsMapping&) = delete; FPartitionsMapping& operator=(const FPartitionsMapping&) = delete; public: - explicit FPartitionsMapping(const int inDim, const int partitions[], const int nbPartitions, - const int ratioForLeaf = 0) + explicit FPartitionsMapping(const int inDim, const int partitions[], const int inNbPartitions) : dim(inDim), - nbCells(0), cells(nullptr), - nbLeaves(0), leaves(nullptr), - totalNbBlocks(0){ + nbPartitions(inNbPartitions), + nbCells(inNbPartitions*inNbPartitions), + cells(nullptr){ FAssertLF(nbPartitions <= inDim); FAssertLF(1 <= nbPartitions); - for(int idxPartRow = 0 ; idxPartRow < nbPartitions ; ++idxPartRow){ - for(int idxPartCol = 0 ; idxPartCol < nbPartitions ; ++idxPartCol){ - if(idxPartRow == idxPartCol - || partitions[idxPartRow]*partitions[idxPartCol] < ratioForLeaf){ - nbLeaves += 1; - } - else{ - nbCells += 1; - } - totalNbBlocks += 1; - } + std::unique_ptr<int[]> partitionsOffset(new int[nbPartitions]); + partitionsOffset[0] = 0; + for(int idxPart = 1 ; idxPart < nbPartitions ; ++idxPart){ + partitionsOffset[idxPart] = partitionsOffset[idxPart-1] + partitions[idxPart-1]; } - leaves = new LeafNode[nbLeaves]; cells = new CellNode[nbCells]; - int idxLeaf = 0; - int idxCell = 0; - int offsetRows = 0; - for(int idxPartRow = 0 ; idxPartRow < nbPartitions ; ++idxPartRow){ - int offsetCols = 0; - for(int idxPartCol = 0 ; idxPartCol < nbPartitions ; ++idxPartCol){ - if(idxPartRow == idxPartCol - || partitions[idxPartRow]*partitions[idxPartCol] < ratioForLeaf){ - leaves[idxLeaf].infos.row = offsetRows; - leaves[idxLeaf].infos.col = offsetCols; - leaves[idxLeaf].infos.nbRows = partitions[idxPartRow]; - leaves[idxLeaf].infos.nbCols = partitions[idxPartCol]; - leaves[idxLeaf].infos.level = 0; - idxLeaf += 1; - } - else{ - cells[idxCell].infos.row = offsetRows; - cells[idxCell].infos.col = offsetCols; - cells[idxCell].infos.nbRows = partitions[idxPartRow]; - cells[idxCell].infos.nbCols = partitions[idxPartCol]; - cells[idxCell].infos.level = 1; - idxCell += 1; - } - offsetCols += partitions[idxPartCol]; + for(int idxPartCol = 0 ; idxPartCol < nbPartitions ; ++idxPartCol){ + for(int idxPartRow = 0 ; idxPartRow < nbPartitions ; ++idxPartRow){ + cells[idxPartCol*nbPartitions + idxPartRow].infos.row = partitionsOffset[idxPartRow]; + cells[idxPartCol*nbPartitions + idxPartRow].infos.col = partitionsOffset[idxPartCol]; + cells[idxPartCol*nbPartitions + idxPartRow].infos.nbRows = partitions[idxPartRow]; + cells[idxPartCol*nbPartitions + idxPartRow].infos.nbCols = partitions[idxPartCol]; + cells[idxPartCol*nbPartitions + idxPartRow].infos.level = 0; } - FAssertLF(offsetCols == dim); - offsetRows += partitions[idxPartRow]; } - FAssertLF(offsetRows == dim); } ~FPartitionsMapping(){ delete[] cells; - delete[] leaves; } int getNbBlocks() const { - return totalNbBlocks; + return nbCells; + } + + CellClass& getCell(const int idxRowPart, const int idxColPart){ + return cells[idxColPart*nbPartitions + idxRowPart].cell; + } + + const CellClass& getCell(const int idxRowPart, const int idxColPart) const { + return cells[idxColPart*nbPartitions + idxRowPart].cell; + } + + const FBlockDescriptor& getCellInfo(const int idxRowPart, const int idxColPart) const { + return cells[idxColPart*nbPartitions + idxRowPart].infos; } void forAllBlocksDescriptor(std::function<void(const FBlockDescriptor&)> callback){ for(int idxCell = 0 ; idxCell < nbCells ; ++idxCell){ callback(cells[idxCell].infos); } - for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){ - callback(leaves[idxLeaf].infos); - } } void forAllCellBlocks(std::function<void(const FBlockDescriptor&, CellClass&)> callback){ @@ -137,12 +107,6 @@ public: } } - void forAllLeafBlocks(std::function<void(const FBlockDescriptor&, LeafClass&)> callback){ - for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){ - callback(leaves[idxLeaf].infos, - leaves[idxLeaf].leaf); - } - } template <class MatrixClass> void fillBlocks(MatrixClass& matrix){ @@ -155,15 +119,6 @@ public: ), cells[idxCell].infos.level); } - for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){ - 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); - } } template <class MatrixClass> @@ -177,15 +132,6 @@ public: ), cells[idxCell].infos.level); } - for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){ - 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 { @@ -193,10 +139,6 @@ public: cells[idxCell].cell.gemv(&res[cells[idxCell].infos.row], &vec[cells[idxCell].infos.col]); } - for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){ - 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 { @@ -205,11 +147,6 @@ public: &mat[cells[idxCell].infos.row], nbRhs, dim); } - for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){ - leaves[idxLeaf].leaf.gemm(&res[leaves[idxLeaf].infos.col], - &mat[leaves[idxLeaf].infos.row], - nbRhs, dim); - } } }; diff --git a/Addons/HMat/Src/Viewers/FDenseBlockPermWrapper.hpp b/Addons/HMat/Src/Viewers/FDenseBlockPermWrapper.hpp index 1cd783caf..2604814c8 100644 --- a/Addons/HMat/Src/Viewers/FDenseBlockPermWrapper.hpp +++ b/Addons/HMat/Src/Viewers/FDenseBlockPermWrapper.hpp @@ -3,19 +3,20 @@ // @SCALFMM_PRIVATE -template <class FReal> +template <class FReal, class SrcMatrixClass > class FDenseBlockPermWrapper{ protected: - const FReal* values; + const SrcMatrixClass& matrix; + const int row; + const int col; const int nbRows; const int nbCols; - const int leadingDim; - const int* permuteVector; public: - FDenseBlockPermWrapper(const FReal* inValues, const int inNbRows, const int inNbCols, - const int inLeading, const int* inPermuteVector) - : values(inValues), nbRows(inNbRows), nbCols(inNbCols), leadingDim(inLeading), permuteVector(inPermuteVector){ + FDenseBlockPermWrapper(const SrcMatrixClass& inMatrix, const int inRow, const int inCol, + const int inNbRows, const int inNbCols) + : matrix(inMatrix), row(inRow), col(inCol), + nbRows(inNbRows), nbCols(inNbCols) { } int getNbRows() const { @@ -29,7 +30,7 @@ public: FReal getValue(const int rowIdx, const int colIdx) const{ FAssertLF(rowIdx < nbRows); FAssertLF(colIdx < nbCols); - return values[permuteVector[colIdx]*leadingDim + permuteVector[rowIdx]]; + return matrix.getVal(row+rowIdx, col+colIdx); } constexpr bool existsForReal() const{ diff --git a/Addons/HMat/Src/Viewers/FMatDensePerm.hpp b/Addons/HMat/Src/Viewers/FMatDensePerm.hpp index 184cdf972..0fb3ade0e 100644 --- a/Addons/HMat/Src/Viewers/FMatDensePerm.hpp +++ b/Addons/HMat/Src/Viewers/FMatDensePerm.hpp @@ -21,7 +21,7 @@ protected: FMatDensePerm(const FMatDensePerm&) = delete; FMatDensePerm& operator=(const FMatDensePerm&) = delete; public: - using BlockDescriptor = FDenseBlockPermWrapper<FReal>; + using BlockDescriptor = FDenseBlockPermWrapper<FReal, FMatDensePerm<FReal> >; explicit FMatDensePerm(const int inDim, const FReal* inValues = nullptr, const int* inPermOrigToNew = nullptr) : matDim(inDim), values(nullptr){ @@ -79,25 +79,25 @@ public: } const FReal& getVal(const int idxRow , const int idxCol) const{ - return values[idxCol*matDim+idxRow]; + return values[permOrigToNew[idxCol]*matDim+permOrigToNew[idxRow]]; } FReal& getVal(const int idxRow , const int idxCol) { - return values[idxCol*matDim+idxRow]; + return values[permOrigToNew[idxCol]*matDim+permOrigToNew[idxRow]]; } void setVal(const int idxRow , const int idxCol, const FReal& val) { - values[idxCol*matDim+idxRow] = val; + values[permOrigToNew[idxCol]*matDim+permOrigToNew[idxRow]] = val; } - FDenseBlockPermWrapper<FReal> getBlock(const int rowIdx, const int colIdx, const int nbRows, const int nbCols) const { + FDenseBlockPermWrapper<FReal, FMatDensePerm<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 < nbRows); FAssertLF(0 < nbCols); FAssertLF(rowIdx + nbRows <= matDim); FAssertLF(colIdx + nbCols <= matDim); - return FDenseBlockPermWrapper<FReal>(&values[colIdx*matDim+rowIdx], nbRows, nbCols, matDim, permOrigToNew); + return FDenseBlockPermWrapper<FReal, FMatDensePerm<FReal> >(*this, rowIdx, colIdx, nbRows, nbCols); } }; diff --git a/Addons/HMat/Tests/testPartitionsMapping.cpp b/Addons/HMat/Tests/testPartitionsMapping.cpp index 58376fc83..cede5a9f1 100644 --- a/Addons/HMat/Tests/testPartitionsMapping.cpp +++ b/Addons/HMat/Tests/testPartitionsMapping.cpp @@ -51,9 +51,8 @@ int main(int argc, char** argv){ std::cout << "Config : outputdir = " << outputdir << "\n"; typedef double FReal; - typedef FDenseBlock<FReal> LeafClass; typedef FDenseBlock<FReal> CellClass; - typedef FPartitionsMapping<FReal, LeafClass, CellClass> GridClass; + typedef FPartitionsMapping<FReal, CellClass> GridClass; { std::unique_ptr<int[]> partitions(new int[nbPartitions]); diff --git a/Addons/HMat/Tests/testPartitionsMappingGemv.cpp b/Addons/HMat/Tests/testPartitionsMappingGemv.cpp index 8fd3827e0..7b4d6918e 100644 --- a/Addons/HMat/Tests/testPartitionsMappingGemv.cpp +++ b/Addons/HMat/Tests/testPartitionsMappingGemv.cpp @@ -70,9 +70,8 @@ int main(int argc, char** argv){ } { - typedef FDenseBlock<FReal> LeafClass; typedef FDenseBlock<FReal> CellClass; - typedef FPartitionsMapping<FReal, LeafClass, CellClass> GridClass; + typedef FPartitionsMapping<FReal, CellClass> GridClass; std::unique_ptr<int[]> partitions(new int[nbPartitions]); { diff --git a/Addons/HMat/Tests/testPartitionsMappingGemvBlock.cpp b/Addons/HMat/Tests/testPartitionsMappingGemvBlock.cpp new file mode 100644 index 000000000..fcb9f951b --- /dev/null +++ b/Addons/HMat/Tests/testPartitionsMappingGemvBlock.cpp @@ -0,0 +1,127 @@ +// =================================================================================== +// 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 + +#include "../Src/Containers/FPartitionsMapping.hpp" +#include "../Src/Viewers/FMatDensePerm.hpp" +#include "../Src/Blocks/FDenseBlock.hpp" + +#include "Utils/FParameters.hpp" +#include "Utils/FParameterNames.hpp" + +#include <memory> + +int main(int argc, char** argv){ + static const FParameterNames DimParam = { + {"-N", "-nb", "-dim"} , + "Dim of the matrix." + }; + static const FParameterNames PartitionsParam = { + {"-part", "-parts", "-nbparts"} , + "Number of partitions." + }; + + FHelpDescribeAndExit(argc, argv, "Test the bisection.", DimParam, FParameterDefinitions::OctreeHeight); + + const int dim = FParameters::getValue(argc, argv, DimParam.options, 100); + const int nbPartitions = FParameters::getValue(argc, argv, PartitionsParam.options, 5); + + std::cout << "Config : dim = " << dim << "\n"; + std::cout << "Config : nbPartitions = " << nbPartitions << "\n"; + + + typedef double FReal; + typedef FMatDensePerm<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 FDenseBlock<FReal> CellClass; + typedef FPartitionsMapping<FReal, CellClass> GridClass; + + std::unique_ptr<int[]> partitions(new int[nbPartitions]); + { + int nbValuesLeft = dim; + for(int idxPartition = 0 ; idxPartition < nbPartitions-1 ; ++idxPartition){ + partitions[idxPartition] = FMath::Max(1, int(drand48()*(nbValuesLeft-(nbPartitions-idxPartition)))); + nbValuesLeft -= partitions[idxPartition]; + } + partitions[nbPartitions-1] = nbValuesLeft; + } + + GridClass grid(dim, partitions.get(), nbPartitions); + { // Here we fill the block manually + // We consider a fack permutation + std::unique_ptr<int[]> permutations(new int[dim]); + for(int idx = 0 ; idx < dim ; ++idx){ + permutations[idx] = idx; + } + // Set permutation to matrix + matrix.setPermutOrigToNew(permutations.get()); + + // We iterate on the blocks + for(int idxColBlock = 0 ; idxColBlock < nbPartitions ; ++idxColBlock){ + for(int idxRowBlock = 0 ; idxRowBlock < nbPartitions ; ++idxRowBlock){ + // We get the corresponding class + //>> CellClass& cl = grid.getCell(idxRowBlock, idxColBlock); + const FBlockDescriptor& info = grid.getCellInfo(idxRowBlock, idxColBlock); + + // We iterate on its values + for(int idxColVal = 0 ; idxColVal < info.nbCols ; ++idxColVal){ + for(int idxRowVal = 0 ; idxRowVal < info.nbRows ; ++idxRowVal){ + //>> const FReal srcVal = matrix.getVal(idxRowVal, idxColVal); + } + } + } + } + } + + 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 partitions mapping, Error = " << testDense << "\n"; + } + + return 0; +} + + + -- GitLab