Commit 506f3a3f authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files

Merge branch 'master' of git+ssh://scm.gforge.inria.fr//gitroot/scalfmm/scalfmm

# By Berenger Bramas (5) and Pierre Blanchard (2)
# Via Pierre Blanchard (2) and Berenger Bramas (1)
* 'master' of git+ssh://scm.gforge.inria.fr//gitroot/scalfmm/scalfmm:
  HMat: fix ChebSymKernel...
  HMat: Factorize ACA routines; Provided routine for recompression (ChebSym keeps its own recompression scheme, namely with storage in K); getRank can be found in namespace FSvd.
  test hmat clustering and rank
  use cell class in the antidiagonal at leaf level
  use inverse permutations
  new a plus equal
  add a custom well separated clustering
parents 90e5b470 9cd38b6a
// ===================================================================================
// 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 FACABLOCK_HPP
#define FACABLOCK_HPP
#include "Utils/FBlas.hpp"
#include "Utils/FAca.hpp"
/*! Choose either \a FULLY_PIVOTED_ACA or \a PARTIALLY_PIVOTED_ACA */
#define FULLY_PIVOTED_ACA
//#define PARTIALLY_PIVOTED_ACA
/*! Choose \a RECOMPRESSED_ACA */
//#define RECOMPRESSED_ACA
template <class FReal, int ORDER = 14>
class FACABlock{
protected:
// members
FReal* block;
FReal* U;
FReal* VT;
int nbRows;
int nbCols;
int level;
unsigned int rank;
FReal accuracy;
FACABlock(const FACABlock&) = delete;
FACABlock& operator=(const FACABlock&) = delete;
public:
FACABlock()
: block(nullptr), U(nullptr), VT(nullptr), nbRows(0), nbCols(0), level(0), rank(0), accuracy(FMath::pow(10.0,static_cast<FReal>(-ORDER))) {
}
// ctor
template <class ViewerClass>
void fill(const ViewerClass& viewer, const int inLevel){
clear();
// Allocate memory
level = inLevel;
nbRows = viewer.getNbRows();
nbCols = viewer.getNbCols();
#if (defined FULLY_PIVOTED_ACA)
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);
}
}
#endif
// SVD specific (col major)
rank = std::min(nbRows,nbCols);
//S = new FReal[rank];
FReal* _U ;// = new FReal[nbRows*nbCols]; // Call to computeSVD() copies block into U
FReal* _V;// = new FReal[rank*nbCols];
#if (defined FULLY_PIVOTED_ACA)
// Perform fully pivoted ACA
FAca::fACA<FReal>(block, nbRows, nbCols, accuracy, _U, _V, rank);
#elif (defined PARTIAL_PIVOTED_ACA)
// TODO
// Perform partially pivoted ACA
//FAca::fACA<FReal>(viewer, nbRows, nbCols, accuracy, _U, _V, rank);
#endif
// display rank
std::cout << "rank after ACA=" << rank << std::endl;
#if (defined RECOMPRESSED_ACA)
// Recompression by QR+SVD
FAca::recompress(_U,_V,nbRows,nbCols,accuracy,U,VT,rank);
// display rank after recompression
std::cout << " rank after QR+SVD=" << rank << std::endl;
#else
// Resize U and VT
U = new FReal[nbRows*rank];
for(int idxRow = 0 ; idxRow < nbRows ; ++idxRow)
for(int idxCol = 0 ; idxCol < rank ; ++idxCol)
U[idxCol*nbRows+idxRow] = _U[idxCol*nbRows+idxRow];
VT = new FReal[rank*nbCols];
for(int idxRow = 0 ; idxRow < rank ; ++idxRow)
for(int idxCol = 0 ; idxCol < nbCols ; ++idxCol)
VT[idxCol*rank+idxRow] = _V[idxRow*nbCols+idxCol]; // transposition
#endif
// Free memory
delete [] _U;
delete [] _V;
};
// dtor
~FACABlock(){
// Free memory
clear();
};
void clear(){
nbRows = 0;
nbCols = 0;
level = 0;
rank = 0;
delete[] block;
block = 0;
delete[] U;
U = 0;
delete[] VT;
VT = 0;
}
void gemv(FReal res[], const FReal vec[], const FReal scale = FReal(1.)) const {
//// Apply (dense) block
//FReal* res_dense = new FReal[nbRows];
//FBlas::copy(nbRows,res,res_dense);
//FBlas::gemva(nbRows, nbCols, scale, const_cast<FReal*>(block), const_cast<FReal*>(vec), res_dense);
// Apply low-rank block
FReal* VTvec = new FReal[rank];
FBlas::setzero(rank,VTvec);
// Apply VT
FBlas::gemv(rank, nbCols, scale, const_cast<FReal*>(VT), const_cast<FReal*>(vec), VTvec);
// Apply U
FBlas::gemva(nbRows, rank, scale, const_cast<FReal*>(U), const_cast<FReal*>(VTvec), res);
}
void gemm(FReal res[], const FReal mat[], const int nbRhs, const FReal scale = FReal(1.)) const {
//// Apply (dense) block
//FBlas::gemma(nbRows, nbCols, nbRhs, scale, const_cast<FReal*>(block), nbRows, const_cast<FReal*>(mat), nbCols, res, nbRows);
// Apply low-rank block
FReal* VTmat = new FReal[nbCols*nbRhs];
// Apply VT
FBlas::gemm(rank, nbCols, nbRhs, scale, const_cast<FReal*>(VT), rank, const_cast<FReal*>(mat), nbCols, VTmat, rank);
// Apply U
FBlas::gemma(nbRows, rank, nbRhs, scale, const_cast<FReal*>(U), nbRows, const_cast<FReal*>(VTmat), rank, res, nbRows);
}
};
#endif // FACABLOCK_HPP
// [--END--]
......@@ -222,7 +222,9 @@ public:
}
int getRank() const{
return rank;
}
};
#endif // FSVDBLOCK_HPP
......
......@@ -118,7 +118,7 @@ public:
free(croot);
}
void fillClusterTree(FClusterTree<FReal>* ctree){
void fillClusterTree(FClusterTree<FReal>* ctree) const {
int* permsOrigToNew = new int[dim];
int* permsNewToOrig = new int[dim];
......
......@@ -29,6 +29,19 @@ protected:
std::pair<int, int>* gclusters;
public:
static FReal GetDefaultRadius(const int inDim, const FReal inDistMat[]){
FReal avg = 0;
FReal min = std::numeric_limits<FReal>::max();
const FReal dimSquare = FReal(inDim*inDim);
for(int idxCol = 0 ; idxCol < inDim ; ++idxCol){
for(int idxRow = 0 ; idxRow < inDim ; ++idxRow){
avg += (inDistMat[idxCol*inDim + idxRow]/dimSquare);
min = (min < inDistMat[idxCol*inDim + idxRow]? min:inDistMat[idxCol*inDim + idxRow]);
}
}
return (avg+min) / FReal(2);
}
FGraphThreshold(const int inDim, const FReal inDistMat[], const FReal inTreshold)
: dim(inDim), treshold(inTreshold),
permutations(new int[dim]),
......@@ -211,7 +224,7 @@ public:
delete[] gclusters;
}
void fillClusterTree(FClusterTree<FReal>* ctree){
void fillClusterTree(FClusterTree<FReal>* ctree) const {
int* permsOrigToNew = new int[dim];
int* permsNewToOrig = new int[dim];
......
#ifndef FMAXDISTCUT_HPP
#define FMAXDISTCUT_HPP
// @SCALFMM_PRIVATE
#include "./Utils/FGlobal.hpp"
#include "./Utils/FAssert.hpp"
#include "./Utils/FMath.hpp"
#include "./Containers/FBoolArray.hpp"
#include "FClusterTree.hpp"
#include <stack>
#include <vector>
#include <functional>
#include <queue>
#include <limits>
#include <memory>
template <class FReal>
class FMaxDistCut {
protected:
const int dim;
int* permutations;
std::pair<int, int>* gclusters;
public:
FMaxDistCut(const int inDim, const FReal inDistMat[])
: dim(inDim),
permutations(new int[dim]),
gclusters (new std::pair<int, int>[dim-1]){
for(int idx = 0 ; idx < dim ; ++idx){
permutations[idx] = idx;
}
std::queue<std::pair<int, int>> intervals;
intervals.push(std::pair<int,int>(0, dim));
int idxCluster = dim-2;
int idxChildCluster = dim-3;
int countToUnknowns = 0;
while(intervals.size()){
FAssertLF(idxCluster >= 0);
const std::pair<int,int> currentCluster = intervals.front();
intervals.pop();
const int sizeInterval = currentCluster.second-currentCluster.first;
FAssertLF(sizeInterval != 1);
int firstIdxMax = permutations[0 + currentCluster.first];
int secondIdxMax = permutations[1 + currentCluster.first];
for(int idxCol = 0 ; idxCol < sizeInterval ; ++idxCol){
const int idxColReal = permutations[idxCol+currentCluster.first];
for(int idxRow = idxCol+1 ; idxRow < sizeInterval ; ++idxRow){
const int idxRowReal = permutations[idxRow+currentCluster.first];
if(inDistMat[idxColReal*dim+ idxRowReal]
> inDistMat[firstIdxMax*dim + secondIdxMax]){
firstIdxMax = idxColReal;
secondIdxMax = idxRowReal;
}
}
}
{
int idxLeft = 0;
int idxRight = sizeInterval-1;
while(idxLeft <= idxRight){
const int idxLeftReal = permutations[idxLeft + currentCluster.first];
const int idxRightReal = permutations[idxRight + currentCluster.first];
if(inDistMat[idxLeftReal*dim+ firstIdxMax]
< inDistMat[idxLeftReal*dim + secondIdxMax]){
idxLeft += 1;
}
else if(inDistMat[idxRightReal*dim+ firstIdxMax]
>= inDistMat[idxRightReal*dim + secondIdxMax]){
idxRight -= 1;
}
else{
const int unkPerm = permutations[idxLeft+currentCluster.first];
permutations[idxLeft+currentCluster.first] = permutations[idxRight+currentCluster.first];
permutations[idxRight+currentCluster.first] = unkPerm;
idxLeft += 1;
idxRight -= 1;
}
}
// idxLeft is on the first 1
if(idxLeft == 1){
gclusters[idxCluster].first = permutations[currentCluster.first];
countToUnknowns += 1;
FAssertLF(countToUnknowns <= dim);
}
else if(idxLeft > 1){
FAssertLF(idxChildCluster >= 0);
gclusters[idxCluster].first = (-idxChildCluster)-1;
idxChildCluster -= 1;
intervals.push(std::pair<int,int>(currentCluster.first, currentCluster.first+idxLeft));
}
if(idxRight == sizeInterval-2){
gclusters[idxCluster].second = permutations[sizeInterval-1+currentCluster.first];
countToUnknowns += 1;
FAssertLF(countToUnknowns <= dim);
}
else if(idxRight < sizeInterval-2){
FAssertLF(idxChildCluster >= 0);
gclusters[idxCluster].second = (-idxChildCluster)-1;
idxChildCluster -= 1;
intervals.push(std::pair<int,int>(currentCluster.first+idxLeft, currentCluster.first+sizeInterval));
}
}
idxCluster -= 1;
}
FAssertLF(idxCluster == -1);
FAssertLF(idxChildCluster == -1);
FAssertLF(countToUnknowns == dim);
}
~FMaxDistCut(){
delete[] permutations;
delete[] gclusters;
}
void fillClusterTree(FClusterTree<FReal>* ctree) const {
int* permsOrigToNew = new int[dim];
int* permsNewToOrig = new int[dim];
{
std::stack<int> depthFirst;
depthFirst.push(gclusters[dim-2].second);
depthFirst.push(gclusters[dim-2].first);
int idxPerm = 0;
while(depthFirst.size()){
const int current = depthFirst.top();
depthFirst.pop();
if(0 <= current){
permsOrigToNew[current] = idxPerm;
permsNewToOrig[idxPerm] = current;
FAssertLF(permsNewToOrig[idxPerm] == permutations[idxPerm]);
idxPerm += 1;
}
else{
depthFirst.push(gclusters[-1-current].second);
depthFirst.push(gclusters[-1-current].first);
}
}
}
typename FClusterTree<FReal>::Leaf* leaves = new typename FClusterTree<FReal>::Leaf[dim];
{
for(int idxUnk = 0 ; idxUnk < dim ; ++idxUnk){
leaves[idxUnk].id = idxUnk;
leaves[idxUnk].offset= idxUnk;
leaves[idxUnk].size = 1;
}
}
typename FClusterTree<FReal>::Node* clusters = new typename FClusterTree<FReal>::Node[dim-1];
for(int idxCluster = 0 ; idxCluster < dim-1 ; ++idxCluster){
typename FClusterTree<FReal>::Node& currentNd = clusters[idxCluster];
const std::pair<int,int>& srcNd = gclusters[idxCluster];
currentNd.id = (-idxCluster)-1;
currentNd.size = 0;
currentNd.left = srcNd.first;
if(0 <= srcNd.first){
currentNd.size += 1;
leaves[srcNd.first].parent = currentNd.id;
}
else{
currentNd.size += clusters[(-srcNd.first)-1].size;
clusters[(-srcNd.first)-1].parent = currentNd.id;
}
currentNd.right = srcNd.second;
if(0 <= srcNd.second){
currentNd.size += 1;
leaves[srcNd.second].parent = currentNd.id;
}
else{
currentNd.size += clusters[(-srcNd.second)-1].size;
clusters[(-srcNd.second)-1].parent = currentNd.id;
}
currentNd.score = 0;
}
clusters[dim-2].parent = 0;
ctree->setData(dim, dim-1, clusters, dim, leaves, permsOrigToNew, permsNewToOrig);
// We do not deallocate, ctree is in charge of this
}
};
#endif // FMAXDISTCUT_HPP
......@@ -63,10 +63,10 @@ public:
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);
cells = new CellNode*[height];
FSetToZeros(cells, height);
nbCells = new int[height];
FSetToZeros(nbCells, height);
for(int idxLevel = 1 ; idxLevel < height-1 ; ++idxLevel){
const int nbCellsAtLevel = FMath::pow2(idxLevel);
......@@ -87,22 +87,42 @@ public:
cells[idxLevel][idxCell].infos.level = idxLevel;
}
}
nbLeaves = FMath::pow2(height-1);
{
const int idxLevel = height-1;
const int nbCellsAtLevel = nbLeaves;
cells[idxLevel] = new CellNode[nbCellsAtLevel];
nbCells[idxLevel] = nbCellsAtLevel;
totalNbBlocks += nbCellsAtLevel;
const int nbLeavesInDirection = nbLeaves;
for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){
const int rowLeafNumber = (idxLeaf^1);
const int colLeafNumber = idxLeaf;
cells[idxLevel][idxLeaf].infos.row = ((rowLeafNumber*dim)/nbLeavesInDirection);
cells[idxLevel][idxLeaf].infos.col = ((colLeafNumber*dim)/nbLeavesInDirection);
cells[idxLevel][idxLeaf].infos.nbRows = (((rowLeafNumber+1)*dim)/nbLeavesInDirection)
- cells[idxLevel][idxLeaf].infos.row;
cells[idxLevel][idxLeaf].infos.nbCols = (((colLeafNumber+1)*dim)/nbLeavesInDirection)
- cells[idxLevel][idxLeaf].infos.col;
cells[idxLevel][idxLeaf].infos.level = height-1;
}
}
nbLeaves = FMath::pow2(height);
leaves = new LeafNode[nbLeaves];
totalNbBlocks += nbLeaves;
{
const int nbLeavesInDirection = (nbLeaves/2);
const int nbLeavesInDirection = nbLeaves;
for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){
const int rowLeafNumber = ((idxLeaf/4)*2) + (idxLeaf&1?1:0);
const int colLeafNumber = ((idxLeaf/4)*2) + (idxLeaf&2?1:0);
const int rowLeafNumber = idxLeaf;
const int colLeafNumber = idxLeaf;
leaves[idxLeaf].infos.row = ((rowLeafNumber*dim)/nbLeavesInDirection);
leaves[idxLeaf].infos.col = ((colLeafNumber*dim)/nbLeavesInDirection);
leaves[idxLeaf].infos.nbRows = (((rowLeafNumber+1)*dim)/nbLeavesInDirection)
- leaves[idxLeaf].infos.row;
leaves[idxLeaf].infos.nbCols = (((colLeafNumber+1)*dim)/nbLeavesInDirection)
- leaves[idxLeaf].infos.col;
leaves[idxLeaf].infos.level = height-1;
leaves[idxLeaf].infos.level = height;
}
}
}
......@@ -124,26 +144,44 @@ public:
}
FAssertLF(offsetUnknowns[nbPartitions] == dim);
nbLeaves = FMath::pow2(height);
nbLeaves = FMath::pow2(height-1);
leaves = new LeafNode[nbLeaves];
totalNbBlocks += nbLeaves;
{
for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){
const int rowLeafNumber = ((idxLeaf/4)*2) + (idxLeaf&1?1:0);
const int colLeafNumber = ((idxLeaf/4)*2) + (idxLeaf&2?1:0);
const int rowLeafNumber = idxLeaf;
const int colLeafNumber = idxLeaf;
leaves[idxLeaf].infos.row = offsetUnknowns[rowLeafNumber];
leaves[idxLeaf].infos.col = offsetUnknowns[colLeafNumber];
leaves[idxLeaf].infos.nbRows = offsetUnknowns[rowLeafNumber+1] - offsetUnknowns[rowLeafNumber];
leaves[idxLeaf].infos.nbCols = offsetUnknowns[colLeafNumber+1] - offsetUnknowns[colLeafNumber];
leaves[idxLeaf].infos.level = height-1;
leaves[idxLeaf].infos.level = height;
}
}
cells = new CellNode*[height-1];
FSetToZeros(cells, height-1);
nbCells = new int[height-1];
FSetToZeros(nbCells, height-1);
cells = new CellNode*[height];
FSetToZeros(cells, height);
nbCells = new int[height];
FSetToZeros(nbCells, height);
{
const int idxLevel = height-1;
const int nbCellsAtLevel = nbLeaves;
cells[idxLevel] = new CellNode[nbCellsAtLevel];
nbCells[idxLevel] = nbCellsAtLevel;
totalNbBlocks += nbCellsAtLevel;
for(int idxLeaf = 0 ; idxLeaf < nbLeaves ; ++idxLeaf){
const int rowLeafNumber = (idxLeaf^1);
const int colLeafNumber = idxLeaf;
cells[idxLevel][idxLeaf].infos.row = offsetUnknowns[rowLeafNumber];
cells[idxLevel][idxLeaf].infos.col = offsetUnknowns[colLeafNumber];
cells[idxLevel][idxLeaf].infos.nbRows = offsetUnknowns[rowLeafNumber+1] - offsetUnknowns[rowLeafNumber];
cells[idxLevel][idxLeaf].infos.nbCols = offsetUnknowns[colLeafNumber+1] - offsetUnknowns[colLeafNumber];
cells[idxLevel][idxLeaf].infos.level = height-1;
}
}
for(int idxLevel = height-2 ; idxLevel >= 1 ; --idxLevel){
const int nbCellsAtLevel = FMath::pow2(idxLevel);
......@@ -173,7 +211,7 @@ public:
}
~FStaticDiagonalBisection(){
for(int idxLevel = 0 ; idxLevel < height-1 ; ++idxLevel){
for(int idxLevel = 0 ; idxLevel < height ; ++idxLevel){
delete[] cells[idxLevel];
}
delete[] cells;
......@@ -186,7 +224,7 @@ public:
}
void forAllBlocksDescriptor(std::function<void(const FBlockDescriptor&)> callback){
for(int idxLevel = 1 ; idxLevel < height-1 ; ++idxLevel){
for(int idxLevel = 1 ; idxLevel < height ; ++idxLevel){
for(int idxCell = 0 ; idxCell < nbCells[idxLevel] ; ++idxCell){
callback(cells[idxLevel][idxCell].infos);
}
......@@ -197,7 +235,7 @@ public:
}
void forAllCellBlocks(std::function<void(const FBlockDescriptor&, CellClass&)> callback){
for(int idxLevel = 1 ; idxLevel < height-1 ; ++idxLevel){
for(int idxLevel = 1 ; idxLevel < height ; ++idxLevel){
for(int idxCell = 0 ; idxCell < nbCells[idxLevel] ; ++idxCell){
callback(cells[idxLevel][idxCell].infos,
cells[idxLevel][idxCell].cell);
......@@ -214,7 +252,7 @@ public:
template <class MatrixClass>
void fillBlocks(MatrixClass& matrix){
for(int idxLevel = 1 ; idxLevel < height-1 ; ++idxLevel){
for(int idxLevel = 1 ; idxLevel < height ; ++idxLevel){
for(int idxCell = 0 ; idxCell < nbCells[idxLevel] ; ++idxCell){
cells[idxLevel