start working on hmat demo

Addons/HMat/CMakeLists.txt

@@ -43,7 +43,7 @@ if(SCALFMM_ADDON_HMAT)
         file( GLOB hpp_in_dir Src/*.hpp Src/*.h)
         INSTALL( FILES ${hpp_in_dir} DESTINATION include/ScalFmm/HMat )
 
-        file( GLOB_RECURSE source_tests_files Tests/*.c )
+        file( GLOB_RECURSE source_tests_files Tests/*.cpp )
         INCLUDE_DIRECTORIES( ${CMAKE_BINARY_DIR}/Src )
 
         # Then build test files

Addons/HMat/Src/Containers/FMatGrid.hpp

+#ifndef FMATGRID_HPP
+#define FMATGRID_HPP
+
+// @SCALFMM_PRIVATE
+
+#include "Utils/FGlobal.hpp"
+#include "Utils/FAssert.hpp"
+
+
+template <class FReal>
+class FDenseBlockWrapper{
+protected:
+    const FReal* values;
+    const int nbRows;
+    const int nbCols;
+    const int leadingDim;
+
+public:
+    FDenseBlockWrapper(const FReal* inValues, const int inNbRows, const int inNbCols, const int inLeading)
+        : values(inValues), nbRows(inNbRows), nbCols(inNbCols), leadingDim(inLeading){
+    }
+
+    int getNbRows() const {
+        return nbRows;
+    }
+
+    int getNbCols() const {
+        return nbCols;
+    }
+
+    FReal getValue(const int rowIdx, const int colIdx) const{
+        FAssertLF(rowIdx < nbRows);
+        FAssertLF(colIdx < nbCols);
+        return values[colIdx*leadingDim + rowIdx];
+    }
+
+    constexpr bool existsForReal() const{
+        return true;
+    }
+};
+
+
+template <class FReal>
+class FMatGrid {
+protected:
+    int matDim;
+    FReal* values;
+
+    FMatGrid(const FMatGrid&) = delete;
+    FMatGrid& operator=(const FMatGrid&) = delete;
+
+public:
+    using BlockDescriptor = FDenseBlockWrapper<FReal>;
+
+    FMatGrid(const int inDim, const FReal* inValues)
+        : matDim(inDim), values(nullptr){
+        values = new FReal[matDim*matDim];
+        for(int idxVal = 0 ; idxVal < matDim*matDim ; ++idxVal){
+            values[idxVal] = inValues[idxVal];
+        }
+    }
+
+    ~FMatGrid(){
+        delete[] values;
+    }
+
+    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(rowIdx + nbRows <= matDim);
+        FAssertLF(colIdx + nbRows <= matDim);
+        return FDenseBlockWrapper<FReal>(&values[colIdx*matDim+rowIdx], nbRows, nbCols, matDim);
+    }
+};
+
+#endif // FMATGRID_HPP
+

Addons/HMat/Src/Core/FDiv2Bissection.hpp

+#ifndef FDIV2BISSECTION_HPP
+#define FDIV2BISSECTION_HPP
+
+// @SCALFMM_PRIVATE
+
+#include "Utils/FGlobal.hpp"
+#include "Utils/FMath.hpp"
+#include "Utils/FAssert.hpp"
+
+#include "../Utils/FHUtils.hpp"
+
+#include <functional>
+
+struct FBlockDescriptor {
+    int x, y, height, width, level;
+};
+
+template <class FReal, class LeafClass, class CellClass >
+class FDiv2Bissection {
+protected:
+    struct LeafNode {
+        FBlockDescriptor infos;
+        LeafClass leaf;
+    };
+
+    struct CellNode {
+        FBlockDescriptor infos;
+        CellClass cell;
+    };
+
+    const int dim;
+    const int height;
+
+    int* nbCells;
+    CellNode** cells;
+
+    int nbLeaves;
+    LeafNode* leaves;
+
+    int totalNbBlocks;
+
+public:
+    explicit FDiv2Bissection(const int inDim, const int inHeight)
+        : dim(inDim), height(inHeight),
+          nbCells(0), cells(nullptr),
+          nbLeaves(0), leaves(nullptr),
+          totalNbBlocks(0){
+        FAssertLF(FMath::pow2(height) <= inDim);
+
+        cells   = new CellNode*[height-1];
+        FSetToZeros(cells, height-1);
+        nbCells = new int[height-1];
+        FSetToZeros(nbCells, height-1);
+
+        int previousWidth = dim;
+        for(int idxLevel = 1 ; idxLevel < height-1 ; ++idxLevel){
+            const int nbCellsAtLevel = FMath::pow2(idxLevel);
+            cells[idxLevel]   = new CellNode[nbCellsAtLevel];
+            nbCells[idxLevel] = nbCellsAtLevel;
+            totalNbBlocks += nbCellsAtLevel;
+
+            const int blockDim = previousWidth/2;
+            for(int idxCell = 0 ; idxCell < nbCellsAtLevel ; ++idxCell){
+                const int xOffset = blockDim * (idxCell&1? idxCell-1 : idxCell+1);
+                const int yOffset = idxCell * blockDim;
+                cells[idxLevel][idxCell].infos.x = xOffset;
+                cells[idxLevel][idxCell].infos.y = yOffset;
+                cells[idxLevel][idxCell].infos.width = blockDim;
+                cells[idxLevel][idxCell].infos.height= blockDim;
+                cells[idxLevel][idxCell].infos.level = idxLevel;
+            }
+
+            previousWidth = blockDim;
+        }
+
+        nbLeaves = FMath::pow2(height);
+        leaves   = new LeafNode[nbLeaves];
+        totalNbBlocks += nbLeaves;
+        {
+            const int blockDim = previousWidth/2;
+            for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){
+                const int corner = (idxLeaf/4)*previousWidth;
+                int xOffset = corner + blockDim * (idxLeaf&1?1:0);
+                int yOffset = corner + blockDim * (idxLeaf&2?1:0);
+                leaves[idxLeaf].infos.x = xOffset;
+                leaves[idxLeaf].infos.y = yOffset;
+                leaves[idxLeaf].infos.width = blockDim;
+                leaves[idxLeaf].infos.height= blockDim;
+                leaves[idxLeaf].infos.level = height-1;
+            }
+        }
+    }
+
+    ~FDiv2Bissection(){
+        for(int idxLevel = 0 ; idxLevel < height-1 ; ++idxLevel){
+            delete[] cells[idxLevel];
+        }
+        delete[] cells;
+        delete[] nbCells;
+        delete[] leaves;
+    }
+
+    int getNbBlocks() const {
+        return totalNbBlocks;
+    }
+
+    void forAllBlocksDescriptor(std::function<void(const FBlockDescriptor&)> callback){
+        for(int idxLevel = 1 ; idxLevel < height-1 ; ++idxLevel){
+            for(int idxCell = 0 ; idxCell < nbCells[idxLevel] ; ++idxCell){
+                callback(cells[idxLevel][idxCell].infos);
+            }
+        }
+        for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){
+            callback(leaves[idxLeaf].infos);
+        }
+    }
+
+    void forAllCellBlocks(std::function<void(const FBlockDescriptor&, CellClass&)> callback){
+        for(int idxLevel = 1 ; idxLevel < height-1 ; ++idxLevel){
+            for(int idxCell = 0 ; idxCell < nbCells[idxLevel] ; ++idxCell){
+                callback(cells[idxLevel][idxCell].infos,
+                         cells[idxLevel][idxCell].cell);
+            }
+        }
+    }
+
+    void forAllLeafBlocks(std::function<void(const FBlockDescriptor&, LeafClass&)> callback){
+        for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){
+            callback(leaves[idxLeaf].infos,
+                     leaves[idxLeaf].leaf);
+        }
+    }
+};
+
+
+#endif // FDIV2BISSECTION_HPP
+

Addons/HMat/Src/Utils/FHUtils.hpp

+#ifndef FHUTILS_HPP
+#define FHUTILS_HPP
+
+#include "Utils/FGlobal.hpp"
+
+#include <cstring>
+
+template <class Type>
+void FSetToZeros(Type array[], const int length){
+    memset(array, 0, length*sizeof(Type));
+}
+
+#endif // FHUTILS_HPP
+

Addons/HMat/Src/Utils/FSvgRect.hpp

+#ifndef FSVGRECT_HPP
+#define FSVGRECT_HPP
+
+#include "Utils/FGlobal.hpp"
+#include "Utils/FAssert.hpp"
+
+#include <vector>
+#include <algorithm>
+
+
+class FSvgRect {
+protected:
+    int dim;
+    int svgSize;
+    int margin;
+
+    int rectCounter;
+    double pixelperunknown;
+    int descriptiony;
+
+    FILE* fsvg;
+
+    std::vector<std::tuple<int, int, int>> colors;
+
+public:
+    FSvgRect(const char inFilename[], const int inDim,
+             const int inSvgSize = 2048, const int inMargin = 50)
+        : dim(inDim), svgSize(inSvgSize), margin(inMargin), rectCounter(0),
+          pixelperunknown(0), descriptiony(0), fsvg(NULL){
+        fsvg = fopen(inFilename, "w");
+        FAssertLF(fsvg);
+
+        pixelperunknown = double(svgSize-margin-margin*2)/double(dim);
+
+        fprintf(fsvg, "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n");
+        fprintf(fsvg, "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.0//EN\" \"http://www.w3.org/TR/2001/REC-SVG-20010904/DTD/svg10.dtd\">\n"
+                "<svg xmlns=\"http://www.w3.org/2000/svg\"  xmlns:xlink=\"http://www.w3.org/1999/xlink\"  width=\"%d\" height=\"%d\">\n",
+                svgSize, svgSize);
+        fprintf(fsvg, "<rect x=\"%d\" y=\"%d\" width=\"%d\" height=\"%d\" style=\"fill:rgb(%d,%d,%d);stroke-width:1;stroke:rgb(0,0,0)\" />\n",
+                            0, 0,svgSize, svgSize, 255, 255, 255);
+
+        descriptiony = int(svgSize - margin);
+        fprintf(fsvg, "<text x=\"%d\" y=\"%d\" fill=\"red\" id=\"res\" font-size=\"55\">Description: </text>\n", 0, descriptiony);
+
+    }
+
+    ~FSvgRect(){
+        fprintf(fsvg, "</svg>");
+        fclose(fsvg);
+    }
+
+    void addRect(const int inX, const int inY, const int inWidth, const int inHeight, const int inLevel = 0){
+        if(int(colors.size()) <= inLevel){
+            for(int idxColor = int(colors.size()) ; idxColor <= inLevel ; ++idxColor){
+                colors.push_back(std::tuple<int, int, int>(255*drand48(),255*drand48(),255*drand48()));
+            }
+        }
+        fprintf(fsvg, "<rect x=\"%d\" y=\"%d\" width=\"%d\" height=\"%d\" style=\"fill:rgb(%d,%d,%d);stroke-width:1;stroke:rgb(0,0,0)\" />\n",
+                            int(margin + inX*pixelperunknown),
+                            int(margin + inY*pixelperunknown),
+                            int(inWidth*pixelperunknown),
+                            int(inHeight*pixelperunknown),
+                            std::get<0>(colors[inLevel]),std::get<1>(colors[inLevel]),std::get<2>(colors[inLevel]));
+
+        rectCounter += 1;
+    }
+
+    template <class... T>
+    void addRectWithLegend(const int inX, const int inY, const int inWidth, const int inHeight, int inLevel = -1){
+        if(int(colors.size()) <= inLevel){
+            for(int idxColor = int(colors.size()) ; idxColor <= inLevel ; ++idxColor){
+                colors.push_back(std::tuple<int, int, int>(255*drand48(),255*drand48(),255*drand48()));
+            }
+        }
+
+        fprintf(fsvg, "<text x=\"%d\" y=\"%d\" fill=\"red\" id=\"res%d\" font-size=\"55\">Description: X = %d, Y = %d, W = %d, H = %d, L = % d </text>\n",
+                0, descriptiony + svgSize, rectCounter,
+                inX, inY, inWidth, inHeight, inLevel);
+
+        fprintf(fsvg, "<rect x=\"%d\" y=\"%d\" width=\"%d\" height=\"%d\" style=\"fill:rgb(%d,%d,%d);stroke-width:1;stroke:rgb(0,0,0)\" ",
+                int(margin + inX*pixelperunknown),
+                int(margin + inY*pixelperunknown),
+                int(inWidth*pixelperunknown),
+                int(inHeight*pixelperunknown),
+                std::get<0>(colors[inLevel]),std::get<1>(colors[inLevel]),std::get<2>(colors[inLevel]));
+
+        fprintf(fsvg, " onmousedown=\"evt.target.parentNode.getElementById('res%d').setAttribute('y', '%d');\" ",
+                rectCounter,descriptiony);
+        fprintf(fsvg, " onmouseup=\"evt.target.parentNode.getElementById('res%d').setAttribute('y', '%d');\" ",
+                rectCounter,descriptiony + svgSize);
+
+        fprintf(fsvg, "/>\n");
+
+        rectCounter += 1;
+    }
+};
+
+#endif // FSVGRECT_HPP
+

Addons/HMat/Tests/testDiv2Bissection.cpp

+
+// @SCALFMM_PRIVATE
+
+#include "../Src/Core/FDiv2Bissection.hpp"
+#include "../Src/Containers/FMatGrid.hpp"
+#include "../Src/Utils/FSvgRect.hpp"
+
+#include "Utils/FParameters.hpp"
+#include "Utils/FParameterNames.hpp"
+
+int main(int argc, char** argv){
+    static const FParameterNames SvgOutParam = {
+        {"-fout", "--out"} ,
+         "Svg output filename."
+    };
+    static const FParameterNames DimParam = {
+        {"-N", "-nb", "-dim"} ,
+         "Dim of the matrix."
+    };
+    static const FParameterNames HeightParam = {
+        {"-h", "-height"} ,
+         "Number of dissection (+1)."
+    };
+
+    FHelpDescribeAndExit(argc, argv,
+            "Test the bisection.",SvgOutParam,
+            DimParam,HeightParam);
+
+    const int dim = FParameters::getValue(argc, argv, DimParam.options, 100);
+    const int height = FParameters::getValue(argc, argv, HeightParam.options, 4);
+    const char* outputfile = FParameters::getStr(argc, argv, SvgOutParam.options, "/tmp/example.svg");
+
+    std::cout << "Config : dim = " << dim << "\n";
+    std::cout << "Config : height = " << height << "\n";
+    std::cout << "Config : outputfile = " << outputfile << "\n";
+
+    FDiv2Bissection<double, int, int> bissection(dim, height);
+
+    {
+        FSvgRect output(outputfile, dim);
+
+        bissection.forAllBlocksDescriptor([&](const FBlockDescriptor& info){
+            output.addRectWithLegend(info.x, info.y, info.width, info.height, info.level);
+        });
+    }
+
+//    const double array[] = {1, 2, 3, 4,
+//                             1, 2, 3, 4,
+//                             1, 2, 3, 4,
+//                             1, 2, 3, 4};
+//    FMatGrid<double> denseMat(4 , array);
+
+//    FMatGrid<double>::BlockDescriptor blockDescriptor = denseMat.getBlock(1,1,1,1);
+
+
+    return 0;
+}