Commit a768410d authored by BRAMAS Berenger's avatar BRAMAS Berenger
Browse files
parents bb04a165 649e4168
......@@ -58,7 +58,41 @@ 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];
FReal SqrtSumSigma2 = std::sqrt(sumSigma2);
// set rank to 1
rank = 1;
// increase
FReal sumSigma2r = S[0]*S[0];
while(std::sqrt(sumSigma2r)<(FReal(1.)-epsilon)*SqrtSumSigma2 && rank<maxRank){
sumSigma2r+=S[rank]*S[rank];
rank++;
}
//std::cout << "std::sqrt(sumSigma2r)=" << std::sqrt(sumSigma2r) << std::endl;
//std::cout << "std::sqrt(sumSigma2)=" << std::sqrt(sumSigma2) << std::endl;
//std::cout << "R/S=" << (std::sqrt(sumSigma2)-std::sqrt(sumSigma2r))/std::sqrt(sumSigma2) << std::endl;
}
template <class FReal, int ORDER = 14>
class FSVDBlock{
protected:
// members
......@@ -70,12 +104,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
......@@ -96,16 +131,32 @@ public:
// SVD specific (col major)
rank = std::min(nbRows,nbCols);
U = new FReal[nbRows*nbCols]; // Call to computeSVD() copies block into U
S = new FReal[rank];
VT = new FReal[rank*nbCols];
FReal* _U = new FReal[nbRows*nbCols]; // Call to computeSVD() copies block into U
FReal* _VT = new FReal[rank*nbCols];
FBlas::setzero(int(rank), S);
FBlas::setzero(int(rank*nbCols),VT);
FBlas::setzero(int(rank*nbCols),_VT);
// 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);
computeSVD(nbRows, nbCols, block, S, _U ,_VT);
// Determine numerical rank using prescribed accuracy
computeNumericalRank(rank, S, accuracy);
// TODO Determine numerical rank using prescribed accuracy
// ...
//// display rank
//std::cout << "rank=" << rank << std::endl;
// Resize U and VT
U = new FReal[nbRows*rank]; // Call to computeSVD() copies block into U
VT = new FReal[rank*nbCols];
for(int idxRow = 0 ; idxRow < nbRows ; ++idxRow)
for(int idxCol = 0 ; idxCol < rank ; ++idxCol)
U[idxCol*nbRows+idxRow] = _U[idxCol*nbRows+idxRow];
for(int idxRow = 0 ; idxRow < rank ; ++idxRow)
for(int idxCol = 0 ; idxCol < nbCols ; ++idxCol)
VT[idxCol*rank+idxRow] = _VT[idxCol*std::min(nbRows,nbCols)+idxRow];
delete [] _U;
delete [] _VT;
};
......@@ -131,17 +182,23 @@ public:
}
void gemv(FReal res[], const FReal vec[], const FReal scale = FReal(1.)) const {
//// Apply (dense) block
//FBlas::gemva(nbRows, nbCols, scale, const_cast<FReal*>(block), const_cast<FReal*>(vec), res);
//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[nbCols];
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 S
for(int idxS = 0 ; idxS < rank ; ++idxS)
VTvec[idxS]*=S[idxS];
VTvec[idxS]*=S[idxS];
// Apply U
FBlas::gemva(nbRows, rank, scale, const_cast<FReal*>(U), const_cast<FReal*>(VTvec), res);
......
......@@ -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,7> LeafClass;
typedef FSVDBlock<FReal,7> 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,7> LeafClass;
typedef FSVDBlock<FReal,7> 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;
}
......
......@@ -43,9 +43,12 @@
##########
### List of vendors (BLA_VENDOR) valid in this module
########## List of vendors (BLA_VENDOR) valid in this module
## Goto,ATLAS PhiPACK,CXML,DXML,SunPerf,SCSL,SGIMATH,IBMESSL,Intel10_32 (intel mkl v10 32 bit),Intel10_64lp (intel mkl v10 64 bit,lp thread model, lp64 model),
## Open (for OpenBlas), Eigen (for EigenBlas), Goto, ATLAS PhiPACK,
##  CXML, DXML, SunPerf, SCSL, SGIMATH, IBMESSL,
## Intel10_32 (intel mkl v10 32 bit), Intel10_64lp (intel mkl v10 64 bit,lp thread model, lp64 model),
## Intel10_64lp_seq (intel mkl v10 64 bit,sequential code, lp64 model),
## Intel( older versions of mkl 32 and 64 bit), ACML,ACML_MP,ACML_GPU,Apple, NAS, Generic
## Intel( older versions of mkl 32 and 64 bit),
##  ACML, ACML_MP, ACML_GPU, Apple, NAS, Generic
# C/CXX should be enabled to use Intel mkl
###
# We handle different modes to find the dependency
......@@ -719,6 +722,7 @@ if (BLA_VENDOR STREQUAL "Goto" OR BLA_VENDOR STREQUAL "All")
endif (BLA_VENDOR STREQUAL "Goto" OR BLA_VENDOR STREQUAL "All")
# OpenBlas
if (BLA_VENDOR STREQUAL "Open" OR BLA_VENDOR STREQUAL "All")
if(NOT BLAS_LIBRARIES)
......@@ -736,6 +740,36 @@ if (BLA_VENDOR STREQUAL "Open" OR BLA_VENDOR STREQUAL "All")
endif (BLA_VENDOR STREQUAL "Open" OR BLA_VENDOR STREQUAL "All")
# EigenBlas
if (BLA_VENDOR STREQUAL "Eigen" OR BLA_VENDOR STREQUAL "All")
if(NOT BLAS_LIBRARIES)
# eigenblas (http://eigen.tuxfamily.org/index.php?title=Main_Page)
check_fortran_libraries(
BLAS_LIBRARIES
BLAS
sgemm
""
"eigen_blas"
""
)
endif()
if(NOT BLAS_LIBRARIES)
# eigenblas (http://eigen.tuxfamily.org/index.php?title=Main_Page)
check_fortran_libraries(
BLAS_LIBRARIES
BLAS
sgemm
""
"eigen_blas_static"
""
)
endif()
endif (BLA_VENDOR STREQUAL "Eigen" OR BLA_VENDOR STREQUAL "All")
if (BLA_VENDOR STREQUAL "ATLAS" OR BLA_VENDOR STREQUAL "All")
if(NOT BLAS_LIBRARIES)
......
......@@ -43,7 +43,7 @@
# add a cache variable to let the user specify the BLAS vendor
set(BLA_VENDOR "" CACHE STRING "list of possible BLAS vendor:
Goto, ATLAS PhiPACK, CXML, DXML, SunPerf, SCSL, SGIMATH, IBMESSL,
Open, Eigen, Goto, ATLAS PhiPACK, CXML, DXML, SunPerf, SCSL, SGIMATH, IBMESSL,
Intel10_32 (intel mkl v10 32 bit),
Intel10_64lp (intel mkl v10 64 bit, lp thread model, lp64 model),
Intel10_64lp_seq (intel mkl v10 64 bit, sequential code, lp64 model),
......
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