Commit cd66ed6a authored by BLANCHARD Pierre's avatar BLANCHARD Pierre
Browse files

HMat: Allow truncation in SVDBlocks (default accuracy is 1e-14); Add timings;...

HMat: Allow truncation in SVDBlocks (default accuracy is 1e-14); Add timings; Fixed minor merge conflicts.
parent d96ccef1
......@@ -58,7 +58,36 @@ static void computeSVD(const FSize nbRows, const FSize nbCols, const FReal* A, F
delete[] WORK;
}
template <class FReal>
/*
* Compute SVD $A=USV'$ and return $V'$, $S$ and $U$.
* \param
*/
template<class FReal>
static void computeNumericalRank(int &rank, const FReal* S, const FReal epsilon){
// verbose
const bool verbose = false;
// init
const FSize maxRank = rank;
FReal sumSigma2 = FReal(0.0);
for(int idxRow = 0 ; idxRow < rank ; ++idxRow)
sumSigma2+=S[idxRow]*S[idxRow];
// set rank to 1
rank = 1;
// increase
FReal sumSigma2r = S[0]*S[0];
FReal eps2 = epsilon*epsilon;
while(sumSigma2r<(1-eps2)*sumSigma2 && rank<maxRank){
sumSigma2r+=S[rank]*S[rank];
rank++;
}
}
template <class FReal, int ORDER = 14>
class FSVDBlock{
protected:
// members
......@@ -70,12 +99,13 @@ protected:
int nbCols;
int level;
int rank;
FReal accuracy;
FSVDBlock(const FSVDBlock&) = delete;
FSVDBlock& operator=(const FSVDBlock&) = delete;
public:
FSVDBlock()
: block(nullptr), U(nullptr), S(nullptr), VT(nullptr), nbRows(0), nbCols(0), level(0), rank(0) {
: block(nullptr), U(nullptr), S(nullptr), VT(nullptr), nbRows(0), nbCols(0), level(0), rank(0), accuracy(FMath::pow(10.0,static_cast<FReal>(-ORDER))) {
}
// ctor
......@@ -104,8 +134,11 @@ public:
// Perform decomposition of rectangular block (jobu=O, jobvt=S => only first min(M,N) cols/rows of U/VT are stored)
computeSVD(nbRows, nbCols, block, S, U ,VT);
// TODO Determine numerical rank using prescribed accuracy
// ...
// Determine numerical rank using prescribed accuracy
computeNumericalRank(rank, S, accuracy);
//// display rank
//std::cout << "rank=" << rank << std::endl;
};
......
......@@ -23,6 +23,8 @@
#include "Utils/FParameters.hpp"
#include "Utils/FParameterNames.hpp"
#include "Utils/FTic.hpp"
#include <memory>
int main(int argc, char** argv){
......@@ -33,6 +35,11 @@ int main(int argc, char** argv){
FHelpDescribeAndExit(argc, argv, "Test the bisection.", DimParam, FParameterDefinitions::OctreeHeight);
////////////////////////////////////////////////////////////////////
/// Timers
FTic time;
time.tic();
const int dim = FParameters::getValue(argc, argv, DimParam.options, 100);
const int height = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 4);
......@@ -112,6 +119,9 @@ int main(int argc, char** argv){
std::cout << "Test Dense with partitions, Error = " << testDense << "\n";
}
double tOverall = time.tacAndElapsed();
std::cout << "... took @tOverall = "<< tOverall <<"\n";
return 0;
}
......@@ -19,15 +19,23 @@
#include "../Src/Containers/FStaticDiagonalBisection.hpp"
#include "../Src/Viewers/FMatDense.hpp"
#include "../Src/Blocks/FDenseBlock.hpp"
#include "../Src/Blocks/FSVDBlock.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FParameterNames.hpp"
#include "Utils/FTic.hpp"
#include <memory>
int main(int argc, char** argv){
FHelpDescribeAndExit(argc, argv, "Test the bisection.", FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight);
////////////////////////////////////////////////////////////////////
/// Timers
FTic time;
time.tic();
const char* filename = FParameters::getStr(argc, argv, FParameterDefinitions::InputFile.options, "../Addons/HMat/Data/unitCube1000_ONE_OVER_R.bin");
const int height = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 4);
......@@ -57,8 +65,12 @@ int main(int argc, char** argv){
}
{
typedef FDenseBlock<FReal> LeafClass;
typedef FDenseBlock<FReal> CellClass;
std::cout << "Test Dense:\n";
//typedef FDenseBlock<FReal> LeafClass;
//typedef FDenseBlock<FReal> CellClass;
typedef FSVDBlock<FReal> LeafClass;
typedef FSVDBlock<FReal> CellClass;
typedef FStaticDiagonalBisection<FReal, LeafClass, CellClass> GridClass;
GridClass grid(dim, height);
......@@ -67,16 +79,29 @@ int main(int argc, char** argv){
std::unique_ptr<FReal[]> resDense(new FReal[dim]);
FSetToZeros(resDense.get(), dim);
std::cout << " Perform GEMV ";
FTic timeGEMV;
timeGEMV.tic();
grid.gemv(resDense.get(), vec.get());
double tGEMV = timeGEMV.tacAndElapsed();
std::cout << "... took @tGEMV = "<< tGEMV <<"\n";
FMath::FAccurater<FReal> testDense(resTest.get(), resDense.get(), dim);
std::cout << "Test Dense, Error = " << testDense << "\n";
std::cout << " Error = " << testDense << "\n";
}
{
typedef FDenseBlock<FReal> LeafClass;
typedef FDenseBlock<FReal> CellClass;
std::cout << "Test Dense with partitions:\n";
//typedef FDenseBlock<FReal> LeafClass;
//typedef FDenseBlock<FReal> CellClass;
typedef FSVDBlock<FReal> LeafClass;
typedef FSVDBlock<FReal> CellClass;
typedef FStaticDiagonalBisection<FReal, LeafClass, CellClass> GridClass;
const int nbPartitions = FMath::pow2(height-1);
......@@ -96,13 +121,24 @@ int main(int argc, char** argv){
std::unique_ptr<FReal[]> resDense(new FReal[dim]);
FSetToZeros(resDense.get(), dim);
std::cout << " Perform GEMV ";
FTic timeGEMV;
timeGEMV.tic();
grid.gemv(resDense.get(), vec.get());
double tGEMV = timeGEMV.tacAndElapsed();
std::cout << "... took @tGEMV = "<< tGEMV <<"\n";
FMath::FAccurater<FReal> testDense(resTest.get(), resDense.get(), dim);
std::cout << "Test Dense with partitions, Error = " << testDense << "\n";
std::cout << " Error = " << testDense << "\n";
}
double tOverall = time.tacAndElapsed();
std::cout << "... took @tOverall = "<< tOverall <<"\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