update to build block matrix

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);
-        }
     }
 };
 

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{

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);
     }
 };
 

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]);

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]);
         {

Addons/HMat/Tests/testPartitionsMappingGemvBlock.cpp

+// ===================================================================================
+// 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;
+}
+
+
+