Commit 5b9746b3 authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

do a few indentation matters and memcpy element size bugfix in LapackKernelInterface

parent eb5f4f34
cmake .. -DCMAKE_INSTALL_PREFIX=/home/tonton/usr -DCMAKE_BUILD_TYPE=Debug -DIBBGMRESDR_DEBUG_MODE=ON -DSANITIZE_ADDRESS=ON -DSANITIZE_UNDEFINED=ON
......@@ -14,10 +14,12 @@ set(CMAKE_C_FLAGS "-std=c99")
#Ajout à la liste des repertoires dans lesquels CMAKE cherche ses modules.
list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake_modules/morse/find")
list(APPEND CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake_modules/")
find_package(CBLAS COMPONENTS BLASEXT REQUIRED)
find_package(LAPACKE COMPONENTS LAPACKEXT REQUIRED)
find_package(MPI REQUIRED)
find_package(Sanitizers)
find_package(PkgConfig)
pkg_check_modules(CHAMELEON chameleon REQUIRED)
pkg_check_modules(STARPU starpu-1.1 REQUIRED)
......
......@@ -57,7 +57,7 @@ struct Arnoldi {
//Needed storage classes
Block<Scalar,Primary,KI> Yj(SizeBlock,MaxKSize/*+SizeBlock*/);
Base<Scalar,Primary,KI> & base = DR.base;
Base<Scalar,Primary,KI>& base = DR.base;
Block<Scalar,Primary,KI> R0;
double timestamp;
......
......@@ -511,7 +511,6 @@ void ArnoldiRuhe_CGS(Hess& H,
base.pushOneVect(Wj.getPtr(k));
}
// base.checkOrtho();
H.incCursor(Wj.getSizeBlock());
std::cout<<"Block Arnoldi done"<<std::endl;
}
......@@ -565,7 +564,6 @@ void ArnoldiRuhe_ICGS(Hess& H,
for (int i=0 ; i<Wj.getLeadingDim(); ++i)
Wj.getPtr(k)[i] /= norm;
base.pushOneVect(Wj.getPtr(k));
// base.checkOrtho();
}
H.incCursor(Wj.getSizeBlock());
std::cout<<"Block Arnoldi done"<<std::endl;
......@@ -588,17 +586,15 @@ void ArnoldiBlock(Hess& H,
int SizeBlock = Wj.getSizeBlock();
int dim = base.getLeadingDim();
Scalar * ptrToHess = H.getNewCol(0);
Scalar *ptrToHess = H.getNewCol(0);
//Ortho against vectors inside base
base.ComputeOrtho(ortho,Wj,ptrToHess,H.getLeadingDim());
std::cout<<"Ortho done"<<std::endl;
//QR facto of orthogonalized Wj
Block<Scalar,Primary,KI> Q(SizeBlock,dim);
Block<Scalar,Primary,KI> R(SizeBlock,SizeBlock);
Wj.ComputeQRFacto(Q,R);
//Copy Block R inside hess
......@@ -642,8 +638,6 @@ void ArnoldiBlock(Hess& L,
//already happened at this stage
}
template<class Matrix,
class Hess,
class Scalar,
......@@ -715,7 +709,7 @@ template<class Matrix,
class KI=LapackKernI /* Kernel Interface Concept */
>
void Arnoldi_Ortho(Hess& H,
Base<Scalar,Primary,KI> base,
Base<Scalar,Primary,KI>& base,
Block<Scalar,Primary,KI>& Wj,
ArnOrtho choice,
OrthoChoice ortho,
......
......@@ -76,6 +76,11 @@ private:
sumNbVectInBlock.push_back(sumNbVectInBlock[currentBlockUsed-1]+ inputSizeBlock);
}
Base& operator=(const Base&) = delete;
Base& operator=(Base&&) = delete;
Base(const Base&) = delete;
Base(Base&&) = delete;
public:
/**
* @brief ctor
......@@ -126,6 +131,7 @@ public:
~Base()
{
delete [] datas;
datas = nullptr;
}
/**
......
......@@ -21,6 +21,14 @@ private:
Scalar * datas;/*!< Raw datas*/
int nbVectUsed;/*!< number of vectors set*/
/**
* @brief Move operator is deleted in order to avoid implicit moves
*/
Block& operator=(const Block&) = delete;
Block& operator=(Block&&) = delete;
Block(const Block&) = delete;
Block(Block&&) = delete;
public:
typedef Scalar value_type;
typedef Primary primary_type;
......@@ -61,10 +69,6 @@ public:
memset(datas,0,sizeof(Scalar)*TotalSize);
}
/**
* @brief Move operator is deleted in order to avoid implicit moves
*/
Block<Scalar,Primary,KI> & operator=(const Block<Scalar,Primary,KI> &) = delete;
/**
* @brief Will be used if default ctor is called, thus if we need
......@@ -94,10 +98,12 @@ public:
/**
* @brief DTor
*/
~Block()
virtual ~Block()
{
if (datas)
if (datas) {
delete [] datas;
datas = nullptr;
}
}
/**
......@@ -223,7 +229,7 @@ public:
template<class Matrix>
void ComputeQRFacto_own(Block& Q, // inputs are already allocated
Block& R, // inputs are already allocated
Matrix & A)
Matrix& A)
{
int n = getLeadingDim();
int m = getSizeBlock();
......@@ -514,9 +520,8 @@ public:
getPtr(),getLeadingDim(),
getPtr(),getLeadingDim(),
A.getPtr(),A.getLeadingDim(),Scalar(1),Scalar(0));
for (int i = 0; i < getSizeBlock(); ++i) {
for (int i = 0; i < getSizeBlock(); ++i)
A.getPtr(i)[i] -= Scalar(1);
}
std::pair<Primary,Primary> minmax = A.getMinMaxNorm();
std::cout<<"Ortho of Block ::\n";
std::cout<<"A - In = \t\tMin :: "<<minmax.first<<"\tMax :: "<<minmax.second<<"\n";
......
......@@ -14,9 +14,16 @@ template<typename Scalar,
>
class BlockWP : public Block<Scalar,Primary, KI>
{
private:
int cursor;
typedef Block<Scalar,Primary,KI> Parent;
typedef Block<Scalar,Primary,KI> BlockSPK;
using Parent = Block<Scalar,Primary,KI>;
using BlockSPK = Parent;
BlockWP& operator=(const BlockWP&) = delete;
BlockWP& operator=(BlockWP&&) = delete;
BlockWP(const BlockWP&) = delete;
BlockWP(BlockWP&&) = delete;
public:
BlockWP(int sizeBlock, int sizeVector):
......
......@@ -14,7 +14,8 @@ template<typename Scalar,
typename Primary=Scalar,
typename KI=LapackKernI /* KernelInterface Concept */
>
class HessExtended : public HessGeneric<Scalar>{
class HessExtended : public HessGeneric<Scalar> {
private:
//Alias
using Parent = HessGeneric<Scalar>;
......@@ -33,8 +34,13 @@ class HessExtended : public HessGeneric<Scalar>{
//Storage of what's needed to compte the rhs for Least square
//Ref : IB BGMRes Dr : Page 10.
Block<Scalar,Primary,KI> * Phi;
Block<Scalar,Primary,KI> * Lambda1;
Block<Scalar,Primary,KI> *Phi;
Block<Scalar,Primary,KI> *Lambda1;
HessExtended& operator=(const HessExtended&) = delete;
HessExtended& operator=(HessExtended&&) = delete;
HessExtended(const HessExtended&) = delete;
HessExtended(HessExtended&&) = delete;
public:
/**
......@@ -48,11 +54,15 @@ public:
sumNbVectInBlck.push_back(0);
}
~HessExtended() {
if (Lambda1)
virtual ~HessExtended() {
if (Lambda1) {
delete Lambda1;
if (Phi)
Lambda1 = nullptr;
}
if (Phi) {
delete Phi;
Phi = nullptr;
}
}
void initLambda(Block<Scalar,Primary,KI>& lambda) {
......
......@@ -7,9 +7,16 @@
* class. The aim of this interface is to have a generic pointer to be
* used to read the hessenberg during Deflation at Restart.
*
* THIS IS A BASE CLASS TO DO POLYMORPHISM STUFF !!! MUST HAVE VIRTUAL DESTUCTOR
*/
template<typename Scalar>
class HessGeneric{
private:
HessGeneric& operator=(const HessGeneric&) = delete;
HessGeneric& operator=(HessGeneric&&) = delete;
HessGeneric(const HessGeneric&) = delete;
HessGeneric(HessGeneric&&) = delete;
protected:
int totalSize;
int SizeBlock;
......@@ -44,8 +51,10 @@ public:
virtual ~HessGeneric()
{
if (datas)
if (datas) {
delete [] datas;
datas = nullptr;
}
}
/**
......
......@@ -24,6 +24,7 @@ template<typename Scalar,
typename KI=LapackKernI /* KernelInterface Concept */
>
class HessQR : public HessGeneric<Scalar>{
private:
//Alias
using Parent = HessGeneric<Scalar>;
......@@ -37,6 +38,11 @@ class HessQR : public HessGeneric<Scalar>{
//Current state
int idxCurrBlock;
HessQR& operator=(const HessQR&) = delete;
HessQR& operator=(HessQR&&) = delete;
HessQR(const HessQR&) = delete;
HessQR(HessQR&&) = delete;
public:
HessQR(int inMaxKSize, int inSizeBlock) : Parent(inMaxKSize,inSizeBlock),
RHS(nullptr),
......
......@@ -18,35 +18,39 @@
* by the restart parameter.
*/
template<typename Scalar,
typename KI=LapackKernI /* KernelInterface Concept */
typename KI=LapackKernI /* Kernel Interface Concept */
>
class HessStandard : public HessGeneric<Scalar>{
//Alias
private:
using Parent = HessGeneric<Scalar>;
int currentSize; //Nb of block Column used (if regular gmres,
//number of column)
/** Nb of block Column used (if regular gmres: number of column) */
int currentSize;
std::vector<int> sumNbVectInBlck;
int nb_ev;
public:
HessStandard(int MaxKSize, int SizeBlock,
int nb_ev=0 ) : Parent(MaxKSize,SizeBlock),
currentSize(0),
sumNbVectInBlck(0),
nb_ev(nb_ev) {
HessStandard& operator=(const HessStandard&) = delete;
HessStandard& operator=(HessStandard&&) = delete;
HessStandard(const HessStandard&) = delete;
HessStandard(HessStandard&&) = delete;
public:
HessStandard(int MaxKSize, int SizeBlock, int nb_ev=0 ):
Parent(MaxKSize,SizeBlock),
currentSize(0),
sumNbVectInBlck(0),
nb_ev(nb_ev)
{
sumNbVectInBlck.push_back(0);
if (0 != nb_ev) {
if (nb_ev != 0) {
sumNbVectInBlck.push_back(nb_ev);
currentSize += 1;
}
}
~HessStandard() {
sumNbVectInBlck.clear();
}
void reserveDR(int inNb_ev) {
void reserveDR(int inNb_ev)
{
nb_ev = inNb_ev;
sumNbVectInBlck.push_back(nb_ev);
currentSize += 1;
......@@ -55,26 +59,20 @@ public:
/**
* Return the raw ptr
*/
Scalar * getPtr(int idBlock = 0) {
Scalar* getPtr(int idBlock = 0)
{
return &(Parent::getPtr()[sumNbVectInBlck[idBlock]*Parent::getLeadingDim()]);
}
int getLeadingDim() {
return Parent::getLeadingDim();
}
int getCurrentSize() {
return currentSize;
}
int getNbColUsed() {
return sumNbVectInBlck[currentSize];
}
int getLeadingDim() const { return Parent::getLeadingDim(); }
int getCurrentSize() const { return currentSize; }
int getNbColUsed() const { return sumNbVectInBlck[currentSize]; }
/**
* @brief Increment Cursor, to be used only by Ruhe variant
*/
void incCursor(int inSizeBlock) {
void incCursor(int inSizeBlock)
{
currentSize++;
//Sum prefix
sumNbVectInBlck.push_back(inSizeBlock + sumNbVectInBlck[currentSize-1]);
......@@ -82,14 +80,15 @@ public:
template<typename Primary>
void displayHessBitMap(std::string name = "") {
void displayHessBitMap(std::string name = "")
{
std::cout<<name<<"\tnbLines :: "
<<sumNbVectInBlck[currentSize]+Parent::SizeBlock<<"\t"
<<"nbCols :: "<<sumNbVectInBlck[currentSize]<<"\n";
for(int j=0 ; j<sumNbVectInBlck[currentSize]+Parent::SizeBlock ; ++j) {
for(int i=0 ; i<sumNbVectInBlck[currentSize] ; ++i) {
if (module<Scalar,Primary>(Parent::getPtr()[i*Parent::getLeadingDim()+j])
!= Primary(0)) {
!= Primary(0)) {
std::cout<<"."<<" ";
} else {
std::cout<<" "<<" ";
......@@ -105,7 +104,8 @@ public:
* @brief Function will return ptr of the new col if there is room
* available, and nullptr if not
*/
Scalar * getNewCol(int inSizeBlock) {
Scalar* getNewCol(int inSizeBlock)
{
if (sumNbVectInBlck[currentSize] + inSizeBlock > Parent::MaxKSize) {
std::cout<<"Max reached\n";
std::cout<<"sumNbVectInBlck[currentSize]:: "<<sumNbVectInBlck[currentSize]
......@@ -121,13 +121,14 @@ public:
//Here add the last triang block
template<typename Primary>
void addRPartToHess(Block<Scalar,Primary,KI> & R) {
void addRPartToHess(Block<Scalar,Primary,KI> & R)
{
int inSizeBlock = R.getSizeBlock();
Scalar* toWrite = (getPtr(currentSize)
+ sumNbVectInBlck[currentSize]
+ inSizeBlock );
Scalar * toWrite = getPtr(currentSize)
+ sumNbVectInBlck[currentSize] + inSizeBlock;
for(int i=0 ; i< inSizeBlock ; ++i) {
for (int i=0 ; i< inSizeBlock ; ++i) {
memcpy(&toWrite[i*Parent::getLeadingDim()],R.getPtr(i),
inSizeBlock*sizeof(Scalar));
}
......@@ -136,8 +137,8 @@ public:
sumNbVectInBlck.push_back(inSizeBlock + sumNbVectInBlck[currentSize-1]);
}
void displayHess(std::string name = "") {
void displayHess(std::string name = "")
{
std::cout<<name<<" ^{H} \t"<<"Current Size "<<currentSize<<"\n";
for(int i = 0 ; i < sumNbVectInBlck[currentSize] ; ++i) { //Loop over columns
std::cout<<"Col"<<i<<" ";
......@@ -150,8 +151,9 @@ public:
//Call to lapack here
template<typename Primary>
void ComputeLeastSquareSol(Block<Scalar,Primary,KI>& R1, Block<Scalar,Primary,KI> & Yj) {
void ComputeLeastSquareSol(Block<Scalar,Primary,KI>& R1,
Block<Scalar,Primary,KI>& Yj)
{
int nbCols = sumNbVectInBlck[currentSize];
int nbLineUsed = nbCols + Parent::SizeBlock;
......@@ -209,7 +211,8 @@ public:
template<typename Primary>
void ComputeBlockResidual(Block<Scalar,Primary,KI>& Y,
Block<Scalar,Primary,KI>& R,
Block<Scalar,Primary,KI>& outputRes) {
Block<Scalar,Primary,KI>& outputRes)
{
//Cstes
int nbCols = sumNbVectInBlck[currentSize];
int nbLineUsed = nbCols + Parent::SizeBlock;
......@@ -226,11 +229,9 @@ public:
Scalar(-1),Scalar(0));
//outputRes - R
for(int idxRHS=0 ; idxRHS<Parent::SizeBlock ; ++idxRHS) {
for(int i = 0 ; i<R.getLeadingDim() ; ++i) {
for(int idxRHS=0 ; idxRHS<Parent::SizeBlock ; ++idxRHS)
for(int i = 0 ; i<R.getLeadingDim() ; ++i)
outputRes.getPtr(idxRHS)[i] += R.getPtr(idxRHS)[i];
}
}
}
};
......
......@@ -16,6 +16,11 @@ private:
int nbTotCols;
Scalar *datas;
HessStorage& operator=(const HessStorage&) = delete;
HessStorage& operator=(HessStorage&&) = delete;
HessStorage(const HessStorage&) = delete;
HessStorage(HessStorage&&) = delete;
public:
HessStorage():
ldH(0),
......@@ -24,7 +29,13 @@ public:
{
}
~HessStorage() { delete [] datas; }
virtual ~HessStorage()
{
if (datas) {
delete [] datas;
datas = nullptr;
}
}
template<class inH>
void StoreHess(inH& H)
......
......@@ -302,6 +302,7 @@ void LapackKernI::gemm(int m, int n, int k,
cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,
m,n,k,alpha,base,ldBase,block,ldBlock,beta,
res,ldRes);
}
template<> //complex float
......@@ -442,7 +443,7 @@ int LapackKernI::orgqr(int m, int p, float *A, int lda,
bool full)
{
for (int i = 0; i < p; ++i)
memcpy(Q+i*ldq, A+i*lda, m*sizeof(float));
memcpy(Q+i*ldq, A+i*lda, m*sizeof*Q);
if (full)
return LAPACKE_sorgqr(LAPACK_COL_MAJOR, m, m, p, Q, ldq, tau.ptr());
......@@ -458,7 +459,7 @@ int LapackKernI::orgqr(int m, int p,
bool full)
{
for (int i = 0; i < p; ++i)
memcpy(Q+i*ldq, A+i*lda, m*sizeof(float));
memcpy(Q+i*ldq, A+i*lda, m*sizeof*Q);
if (full)
return LAPACKE_dorgqr(LAPACK_COL_MAJOR, m, m, p, Q, ldq, tau.ptr());
......@@ -475,7 +476,7 @@ int LapackKernI::orgqr(int m, int p,
bool full)
{
for (int i = 0; i < p; ++i)
memcpy(Q+i*ldq, A+i*lda, m*sizeof(float));
memcpy(Q+i*ldq, A+i*lda, m*sizeof*Q);
if (full)
return LAPACKE_cungqr(LAPACK_COL_MAJOR, m, m, p, Q, ldq, tau.ptr());
......@@ -491,7 +492,7 @@ int LapackKernI::orgqr(int m, int p,
bool full)
{
for (int i = 0; i < p; ++i)
memcpy(Q+i*ldq, A+i*lda, m*sizeof(float));
memcpy(Q+i*ldq, A+i*lda, m*sizeof*Q);
if (full)
return LAPACKE_zungqr(LAPACK_COL_MAJOR, m, m, p, Q, ldq, tau.ptr());
......
......@@ -106,7 +106,7 @@ bool CheckCurrentSolution(Matrix&A,
for(int j=0; j < dim; ++j)
AxSol.getPtr(i)[j] = B.getPtr(i)[j] - AxSol.getPtr(i)[j];
auto evalRes = AxSol.CheckPrecision(epsilon, A);
std::cout<<"\t\\t\tReal Residual : Min\t"<<evalRes.Min<<" Max\t"<<evalRes.Max<<"\n";
std::cout<<"\t\t\tReal Residual : Min\t"<<evalRes.Min<<" Max\t"<<evalRes.Max<<"\n";
if (evalRes.hasConverged) {
std::cout<<"Convergence !\n";
X0.CopyBlock(Sol);
......
......@@ -49,6 +49,7 @@ foreach(_test ${GMRES_TEST_SRC})
# )
# set_property(TARGET ${_name_exe} PROPERTY LINKER_LANGUAGE Fortran)
target_link_libraries(${_name_exe} ${LIBS_FOR_TESTS})
add_sanitizers(${_name_exe})
install(TARGETS ${_name_exe}
DESTINATION ${CMAKE_INSTALL_PREFIX}/bin/.)
......
......@@ -17,8 +17,6 @@ split(std::string line, char sep = ' ')
return res;
}
template<class MatrixStorage,
class E = typename std::enable_if< // complex types only:
std::is_same<
......@@ -40,7 +38,6 @@ void load_line(MatrixStorage &M,const std::vector<std::string> &v)
M.setValue(line, column, CPLX{real, imag});
}
template<class MatrixStorage,
class E = typename std::enable_if< // real types only:
std::is_same<typename MatrixStorage::value_type, float>::value ||
......
......@@ -93,7 +93,7 @@ int main(int ac, char *av[])
using Arnoldi2 = Arnoldi<Matrix,Scalar,Primary>;
using Arnoldi3 = Arnoldi_IB<Matrix,Scalar,Primary>;
int nb = BGMRes<Arnoldi3,Matrix,Scalar,Primary>(mat,RHS,X0,dim+SizeBlock,restart,
int nb = BGMRes<Arnoldi2,Matrix,Scalar,Primary>(mat,RHS,X0,dim+SizeBlock,restart,
sol,log,
epsilon,OrthoChoice::ICGS,
ArnOrtho::RUHE);
......
......@@ -18,7 +18,7 @@ int main(int ac, char *av[])
for(auto ag : args)
std::cout << ag << "\n";
int SizeBlock=3,nbBlock=10,restart=15;
int SizeBlock=3,nbBlock=10,restart=40;
if (args.size()>1)
SizeBlock = std::stoi(args[1]);
......@@ -69,7 +69,7 @@ int main(int ac, char *av[])
Scalar target(0);
//Logs
Logger<Primary> log,log2;
Logger<Primary> log(true),log2(true);
//Using for the arnoldi procedure choosen
using Arn1 = Arnoldi<Matrix,Scalar,Primary>;
......@@ -77,12 +77,12 @@ int main(int ac, char *av[])
{
//Compute a starting Block
Block<Scalar,Primary,KI> X0(SizeBlock,dim);
Block<Scalar,Primary,KI> sol(SizeBlock,dim);
int nb = BGMRes<Arn1,Matrix,Scalar,Primary>(mat,RHS,X0,maxMVP,restart,
sol,log, epsilon,
OrthoChoice::MGS);
int nb = BGMRes<Arn1>(mat, RHS,X0,maxMVP,restart,
sol,log, epsilon,
OrthoChoice::MGS,
ArnOrtho::BLOCK);
std::cout<<"Total Nb iter done::\t"<<nb<<"\n";
//Check solution
......@@ -125,13 +125,13 @@ int main(int ac, char *av[])
{
//Compute a starting Block
Block<Scalar,Primary,KI> X0(SizeBlock,dim);
Block<Scalar,Primary,KI> sol(SizeBlock,dim);
int nb = BGMResDR<Arn1,Matrix,Scalar,Primary>(mat,RHS,X0,maxMVP,restart,
sol,log2, epsilon,
k,target,
OrthoChoice::MGS);
int nb = BGMResDR<Arn1>(mat,RHS,X0,maxMVP,restart,
sol,log2, epsilon,
k,target,
OrthoChoice::MGS,
ArnOrtho::BLOCK );
std::cout<<"Total Nb iter done::\t"<<nb<<"\n";
//Check solution
......
......@@ -22,39 +22,37 @@ using KI = LapackKernI;
using Matrix = UserInputMatrix<Scalar,Primary,KI>;
using BlockSP = Block<Scalar,Primary,KI>;
const int MaxNumberMatrixVectorProduct = 2000;
struct TestRet {
double elapsed;
int nb;
std::pair<Primary, Primary> minmax;
};
template<class T>
template<class ARNOLDI>
TestRet runTest(Matrix &mat,
BlockSP &RHS,
BlockSP &XExact,
int restart,
int SizeBlock, int dim,
int SizeBlock,
int dim,
int maxMVP,
std::vector<Primary> &epsilon,
ArnOrtho ortho)
{
std::cout<<"\n\tVersion Inexact Breakdown Ruhe\n\n";
BlockSP X0(SizeBlock,dim);
BlockSP sol(SizeBlock,dim);
Logger<Primary> log(true);
Logger<Primary> log;
double tic = Logger<Primary>::GetTime();
int nb = BGMRes<T,Matrix,Scalar,Primary>(
mat, RHS, X0, MaxNumberMatrixVectorProduct,
restart, sol, log, epsilon,
OrthoChoice::CGS, ortho
);
int nb = BGMRes<ARNOLDI>(mat, RHS, X0, maxMVP,
restart, sol, log, epsilon,
OrthoChoice::CGS, ortho );
double elapsed = Logger<Primary>::GetTime() - tic;
std::cout<<"Return value / max_ite : "
<< nb << "/" << MaxNumberMatrixVectorProduct << std::endl;
<< nb << "/" << maxMVP << std::endl;
std::stringstream ss;
ss << "../data/res/young1c_IB_"
......@@ -120,30 +118,33 @@ int main(int argc, char *argv[])
restart = std::stoi(args[2]);
else
restart = dim+1;
int maxMVP = 2000;
if (args.size() > 3)
nb_ev = std::stoi(args[3]);
(void) nb_ev;
std::vector<Primary> epsilon{1e-5};
using Arnoldi2 = Arnoldi_IB<Matrix,Scalar,Primary>;
TestRet r1 = runTest<Arnoldi2>(mat, RHS, XExact,
restart, SizeBlock, dim,
epsilon, ArnOrtho::BLOCK);
TestRet r2 = runTest<Arnoldi2>(mat, RHS, XExact,
restart, SizeBlock, dim,
epsilon, ArnOrtho::RUHE);
using ARN_IB = Arnoldi_IB<Matrix,Scalar,Primary>;
TestRet r1 = runTest<ARN_IB>(mat, RHS, XExact,
restart, SizeBlock, dim, maxMVP,
epsilon, ArnOrtho::BLOCK);
// TestRet r2 = runTest<ARN_IB>(mat, RHS, XExact,
// restart, SizeBlock, dim, maxMVP,
// epsilon, ArnOrtho::RUHE);
//Display timers
std::cout<<"Elapsed time for each method::\n";
std::cout<<"\tElapsed time for IB- Arnoldi BLOCK\t\t"<<r1.elapsed
<<" seconds for "<<r1.nb<<"/"<<MaxNumberMatrixVectorProduct<<" iterations\n"
<<" seconds for "<<r1.nb<<"/"<<maxMVP<<" iterations\n"
<<"Forward error : Min\t"<<r1.minmax.first<<" Max\t"<<r1.minmax.second<<"\n";
std::cout<<"\tElapsed time for IB Arnoldi RUHE\t\t"<<r2.elapsed
<<" seconds for "<<r2.nb<<"/"<<MaxNumberMatrixVectorProduct<<" iterations\n"
<<"Forward error : Min\t"<<r2.minmax.first<<" Max\t"<<r2.minmax.second<<"\n";
// std::cout<<"\tElapsed time for IB Arnoldi RUHE\t\t"<<r2.elapsed
// <<" seconds for "<<r2.nb<<"/"<<maxMVP<<" iterations\n"
// <<"Forward error : Min\t"<<r2.minmax.first<<" Max\t"<<r2.minmax.second<<"\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