diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 7e2aa4eea3a76a54115311cc3becea0cf3146e2e..3a084d2c00d37925d0f709ff097ba1a022e7c916 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -3,14 +3,10 @@ cmake_minimum_required(VERSION 2.8)
#Inside variable, is all my sources
set(GMRES_TEST_SRC
testBlockClass.cpp
- testBlockGMRes.cpp
testBlockSvd.cpp
testBlockWP.cpp
- testBlockArnoldiIb.cpp
testMatrixMarket.cpp
testMatrixIW.cpp
- testIbBGMRes.cpp
- testQRIncrementalGMRes.cpp
testBGMRes.cpp
testMatrixIW_noRHS.cpp
)
diff --git a/test/testBlockGMRes.cpp b/test/testBlockGMRes.cpp
deleted file mode 100644
index a4e33d5da8b98422545fe2063a0cea8fbe218fd4..0000000000000000000000000000000000000000
--- a/test/testBlockGMRes.cpp
+++ /dev/null
@@ -1,189 +0,0 @@
-#include <random>
-#include "../src/Block.hpp"
-#include "../src/BlockGMRes.hpp"
-#include "../src/LapackInterface.hpp"
-#include "../src/Logger.hpp"
-
-/**
- * Matrix Class : need to implement MatBlockVect(), MatBaseproduct(),
- * useRightPreCond(), preCondBlockVect(),preCondBaseProduct() and
- * size() member function in order to works with BlockGMRes algorithm
- *
- * An abstract class could be used to enforce that.
- */
-template<class Scalar>
-class InputMatrix{
- int dim;
- Scalar * data;
-public:
- InputMatrix(int inSize) : dim(inSize),data(nullptr){
- data = new Scalar[dim*dim];
- }
-
- ~InputMatrix(){
- delete [] data;
- }
-
- int size(){
- return dim;
- }
-
- void initMatrix(std::random_device& rd){
- std::mt19937 gen(/*rd()*/0);
- std::uniform_real_distribution<> dis(0, 2);
- for(int i=0 ; i<dim ; ++i){
- for(int j=0 ; j<dim ; ++j){
- data[dim*i+j] = dis(gen);
- }
- }
- // Scalar diag = dim;
- // Scalar no_diag = 1;
- // for(int i=0 ; i<dim ; ++i){
- // for(int j=0 ; j<dim ; ++j){
- // if(i==j)
- // data[i*dim + j] = diag+j;
- // else
- // data[i*dim + j] = no_diag;
- // }
- // }
- }
-
- void MatBlockVect(Block<Scalar>& in, Block<Scalar>& out){
- call_gemm<Scalar>(dim,in.getSizeBlock(),dim,
- data,dim,
- in.getPtr(), in.getLeadingDim(),
- out.getPtr(),out.getLeadingDim(),
- Scalar(1),Scalar(0));
- }
-
- void MatBaseProduct(Scalar * toRead,int nbRHS, Scalar * ptrToWrite){
- call_gemm<Scalar>(dim,nbRHS, dim,
- data,dim,
- toRead,dim,
- ptrToWrite,dim,
- Scalar(1),Scalar(0));
- }
-
- bool useRightPreCond(){
- return false;
- }
- void preCondBlockVect(Block<Scalar>& input, Block<Scalar>& output){
- for(int j=0 ; j<output.getSizeBlock() ; ++j){
- for(int i=0 ; i<dim ; ++i){
- output.getPtr(j)[i] = input.getPtr(j)[i]*(1/data[i*dim+i]);
- }
- }
- }
- void preCondBaseProduct(Scalar * ptr,Block<Scalar>& output){
- for(int j=0 ; j<output.getSizeBlock() ; ++j){
- for(int i=0 ; i<dim ; ++i){
- output.getPtr(j)[i] = ptr[j*dim+i]*1/(data[i*dim+i]);
- }
- }
- }
-};
-
-
-int main(int ac, char ** av){
- std::vector<std::string > args(av,av+ac);
- for(auto ag : args)
- std::cout << ag << "\n";
-
- int SizeBlock=2,nbBlock=20,restart=15;
-
- if (args.size()>1){
- SizeBlock = std::stoi(args[1]);
- }
- if (args.size()>2){
- nbBlock = std::stoi(args[2]);
- }
- restart = nbBlock;
- if (args.size()>3){
- restart = std::stoi(args[3]);
- }
- int dim = nbBlock*SizeBlock;
-
- using Scalar=double;
- using Primary=Scalar;
-
- //Random engine used to generate matrix
- std::random_device rd;
- std::mt19937 gen(/*rd()*/0);
- std::uniform_real_distribution<> dis(0, 2);
-
- //Create matrix
- InputMatrix<Scalar> mat(dim);
- mat.initMatrix(rd);
-
- //Create X_Exact
- Block<Scalar> X_Exact(SizeBlock,dim);
- //Fill it!
- for(int i=0 ; i<SizeBlock ; ++i){
- for(int j=0 ;j<dim; ++j){
- X_Exact.getPtr(i)[j] = dis(gen);
- }
- }
-
- //Compute RHS as Mat*X_Exact
- Block<Scalar> RHS(SizeBlock,dim);
- mat.MatBlockVect(X_Exact,RHS);
-
- //Compute a starting Block
- Block<Scalar> X0(SizeBlock,dim);
- //Fill it!
- for(int i=0 ; i<SizeBlock ; ++i){
- for(int j=0 ;j<dim; ++j){
- X0.getPtr(i)[j] = dis(gen);
- }
- }
-
- Block<Scalar> sol(SizeBlock,dim);
-
- Logger<Primary,5> log;
-
- //Three steps
- //With Restart
- int nb_iter = BlockGMRes<InputMatrix<Scalar>,Scalar>(mat,RHS,
- X0,nbBlock,sol,
- restart,log,
- 0.000001);
-
- //Check solution
- {
- //Compute what i need for my metrics
- //MatxSol : Mat*Sol_computed
- Block<Scalar> MatxSol(SizeBlock,dim);
- mat.MatBlockVect(sol,MatxSol);
- //MatxSolMinusB : Mat*Sol_computed-B
- Block<Scalar> MatxSolMinusB(SizeBlock,dim);
- //X_Diff : X_Exact - Sol_computed
- Block<Scalar> X_Diff(SizeBlock,dim);
- MatxSolMinusB.CopyBlock(MatxSol);
- for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
- for(int j=0 ; j<dim; ++j){
- MatxSolMinusB.getPtr(idxRHS)[j] -= RHS.getPtr(idxRHS)[j];
- X_Diff.getPtr(idxRHS)[j] = X_Exact.getPtr(idxRHS)[j] -
- sol.getPtr(idxRHS)[j];
- }
- }
- //Then display for each vect
- for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
- std::cout<<"Metrics on "<<idxRHS<<"\n";
- Primary ForwardError = X_Diff.getNorm(idxRHS);
- std::cout <<"["<<idxRHS<<"]\t"
- <<"Forward Error (on X) \t" << ForwardError <<"\n";
- Primary SecondMemberError = MatxSolMinusB.getNorm(idxRHS);
- std::cout <<"["<<idxRHS<<"]\t"
- <<"Second Member Error \t" << SecondMemberError <<"\n";
- Primary ForwardErrorNormalized = ForwardError/(X_Exact.getNorm(idxRHS));
- std::cout <<"["<<idxRHS<<"]\t"
- <<"Forward Error Norm. \t" << ForwardErrorNormalized <<"\n";
- Primary normB = RHS.getNorm(idxRHS);
- Primary NormwiseBackwardError = SecondMemberError / normB;
- std::cout <<"["<<idxRHS<<"]\t"
- <<"Normwise Backward Error\t" << NormwiseBackwardError <<"\n\n";
-
- }
- }
- return 0;
-}
diff --git a/test/testMatrixIW.cpp b/test/testMatrixIW.cpp
index 716309e0f3a875b3f7589081d8515a6cf0f8b2cf..ad493845df0a5ad8fe5f2e8c9e6b62a10df30039 100644
--- a/test/testMatrixIW.cpp
+++ b/test/testMatrixIW.cpp
@@ -99,7 +99,7 @@ public:
}
bool useRightPreCond(){
- return false;
+ return true;
}
void preCondBlockVect(Block<Scalar,Primary>& input, Block<Scalar,Primary>& output){
for(int j=0 ; j<output.getSizeBlock() ; ++j){
@@ -183,82 +183,127 @@ int main(int ac, char ** av){
using Arnoldi3 = Arnoldi_QRInc<Matrix,Scalar,Primary>;
double elapsed1,elapsed2,elapsed3;
+ int nb1,nb2,nb3;
{
std::cout<<"\n\tVersion Standard\n\n";
Block<Scalar,Primary> X0(SizeBlock,dim);
Block<Scalar,Primary> sol(SizeBlock,dim);
-
+ Logger<Primary> log;
double tic = Logger<Primary>::GetTime();
- int nb = BGMRes<Arnoldi1,Matrix,Scalar,Primary>(mat,RHS,X0,max_ite,max_ite,sol,
- 0.000001);
+ nb1 = BGMRes<Arnoldi1,Matrix,Scalar,Primary>(mat,RHS,X0,max_ite,max_ite,sol,log,
+ 0.00001);
elapsed1 = Logger<Primary>::GetTime() - tic;
- std::cout<<"Return value / max_ite : "<<nb<<"/"<<max_ite<<"\n";
+ std::cout<<"Return value / max_ite : "<<nb1<<"/"<<max_ite<<"\n";
+ std::string locationRes("../data/res/MatCone.res");
+ std::ofstream fileRes(locationRes,std::ios::out);
+ log.forEachIte([&](int ite, int nbMatVect, int currP,
+ Primary min, Primary max, double time)
+ {
+ fileRes<<ite<<"\t"
+ <<nbMatVect<<"\t"
+ <<currP<<"\t"
+ <<min<<"\t"
+ <<max<<"\t"
+ <<time<<"\n";
+ });
+
}
{
std::cout<<"\tVersion with Inexact Breakdown\n\n";
Block<Scalar,Primary> X0(SizeBlock,dim);
Block<Scalar,Primary> sol(SizeBlock,dim);
-
+ Logger<Primary> log;
double tic = Logger<Primary>::GetTime();
- int nb = BGMRes<Arnoldi2,Matrix,Scalar,Primary>(mat,RHS,X0,maxiteIB,max_ite,sol,
- 0.000001);
+ nb2 = BGMRes<Arnoldi2,Matrix,Scalar,Primary>(mat,RHS,X0,maxiteIB,maxiteIB,
+ sol,log,
+ 0.00001);
elapsed2 = Logger<Primary>::GetTime() - tic;
+ //Check solution
+
+ std::cout<<"Return value / max_ite : "<<nb2<<"/"<<maxiteIB<<"\n";
+ std::string locationRes("../data/res/MatCone_IB.res");
+ std::ofstream fileRes(locationRes,std::ios::out);
+ log.forEachIte([&](int ite, int nbMatVect, int currP,
+ Primary min, Primary max, double time)
+ {
+ fileRes<<ite<<"\t"
+ <<nbMatVect<<"\t"
+ <<currP<<"\t"
+ <<min<<"\t"
+ <<max<<"\t"
+ <<time<<"\n";
+ });
- std::cout<<"Return value / max_ite : "<<nb<<"/"<<maxiteIB<<"\n";
}
{
std::cout<<"\tVersion with Optimized Least Square\n\n";
Block<Scalar,Primary> X0(SizeBlock,dim);
Block<Scalar,Primary> sol(SizeBlock,dim);
-
+ Logger<Primary> log;
double tic = Logger<Primary>::GetTime();
- int nb = BGMRes<Arnoldi3,Matrix,Scalar,Primary>(mat,RHS,X0,max_ite,max_ite,sol,
- 0.000001);
+ nb3 = BGMRes<Arnoldi3,Matrix,Scalar,Primary>(mat,RHS,X0,max_ite,max_ite,sol,log,
+ 0.00001);
elapsed3 = Logger<Primary>::GetTime() - tic;
- std::cout<<"Return value / max_ite : "<<nb<<"/"<<max_ite<<"\n";
+ std::cout<<"Return value / max_ite : "<<nb3<<"/"<<max_ite<<"\n";
+ std::string locationRes("../data/res/MatCone_QR.res");
+ std::ofstream fileRes(locationRes,std::ios::out);
+ log.forEachIte([&](int ite, int nbMatVect, int currP,
+ Primary min, Primary max, double time)
+ {
+ fileRes<<ite<<"\t"
+ <<nbMatVect<<"\t"
+ <<currP<<"\t"
+ <<min<<"\t"
+ <<max<<"\t"
+ <<time<<"\n";
+ });
}
+ // {
+ // //Compute what i need for my metrics
+ // //MatxSol : Mat*Sol_computed
+ // Block<Scalar,Primary> MatxSol(SizeBlock,dim);
+ // mat.MatBlockVect(sol,MatxSol);
+ // //MatxSolMinusB : Mat*Sol_computed-B
+ // Block<Scalar,Primary> MatxSolMinusB(SizeBlock,dim);
+ // MatxSolMinusB.CopyBlock(MatxSol);
+ // for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
+ // for(int j=0 ; j<dim; ++j){
+ // MatxSolMinusB.getPtr(idxRHS)[j] -= RHS.getPtr(idxRHS)[j];
+ // }
+ // }
+ // //Then display for each vect
+ // for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
+ // std::cout<<"Metrics on "<<idxRHS<<"\n";
+ // Primary SecondMemberError = MatxSolMinusB.getNorm(idxRHS);
+ // std::cout <<"["<<idxRHS<<"]\t"
+ // <<"Second Member Error \t" << SecondMemberError <<"\n";
+ // Primary normB = RHS.getNorm(idxRHS);
+ // Primary NormwiseBackwardError = SecondMemberError / normB;
+ // std::cout <<"["<<idxRHS<<"]\t"
+ // <<"Normwise Backward Error\t" << NormwiseBackwardError <<"\n\n";
+
+ // }
+ // }
+
+
//Display timers
std::cout<<"Elapsed time for each method::\n";
- std::cout<<"Elapsed time for Standard Version\t"<<elapsed1<<"\n";
- std::cout<<"Elapsed time for Inexact Breakdown\t"<<elapsed2<<"\n";
- std::cout<<"Elapsed time for Optimized Least Square\t"<<elapsed3<<"\n";
- //Check solution
- // {
- // //Compute what i need for my metrics
- // //MatxSol : Mat*Sol_computed
- // Block<Scalar,Primary> MatxSol(SizeBlock,dim);
- // mat.MatBlockVect(sol,MatxSol);
- // //MatxSolMinusB : Mat*Sol_computed-B
- // Block<Scalar,Primary> MatxSolMinusB(SizeBlock,dim);
- // MatxSolMinusB.CopyBlock(MatxSol);
- // for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
- // for(int j=0 ; j<dim; ++j){
- // MatxSolMinusB.getPtr(idxRHS)[j] -= RHS.getPtr(idxRHS)[j];
- // }
- // }
- // //Then display for each vect
- // for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
- // std::cout<<"Metrics on "<<idxRHS<<"\n";
- // Primary SecondMemberError = MatxSolMinusB.getNorm(idxRHS);
- // std::cout <<"["<<idxRHS<<"]\t"
- // <<"Second Member Error \t" << SecondMemberError <<"\n";
- // Primary normB = RHS.getNorm(idxRHS);
- // Primary NormwiseBackwardError = SecondMemberError / normB;
- // std::cout <<"["<<idxRHS<<"]\t"
- // <<"Normwise Backward Error\t" << NormwiseBackwardError <<"\n\n";
-
- // }
- // }
+ std::cout<<"\tElapsed time for Standard Version\t\t"<<elapsed1
+ <<" seconds for "<<nb1<<"/"<<max_ite<<" iterations\n";
+ std::cout<<"\tElapsed time for Inexact Breakdown\t\t"<<elapsed2
+ <<" seconds for "<<nb2<<"/"<<max_ite<<" iterations\n";
+ std::cout<<"\tElapsed time for Optimized Least Square\t\t"<<elapsed3
+ <<" seconds for "<<nb3<<"/"<<max_ite<<" iterations\n";
return 0;
}
diff --git a/test/testMatrixMarket.cpp b/test/testMatrixMarket.cpp
index d6bee364a9ab80681ac60674b5d03b5273f2e95d..5113f3114202e85d5e0fda2cd86c44aca4be9d6e 100644
--- a/test/testMatrixMarket.cpp
+++ b/test/testMatrixMarket.cpp
@@ -2,12 +2,13 @@
#include <iostream>
#include <fstream>
#include "../src/Block.hpp"
-#include "../src/BlockGMRes.hpp"
-#include "../src/BlockGMResQR.hpp"
#include "../src/LapackInterface.hpp"
#include "../src/Logger.hpp"
-#include "../src/BlockArnoldiIb.hpp"
-#include "../src/IBGMRes.hpp"
+#include "../src/BGMRes.hpp"
+#include "../src/Arnoldi.hpp"
+#include "../src/Arnoldi_IB.hpp"
+#include "../src/Arnoldi_QRInc.hpp"
+
/**
* Matrix Class : need to implement MatBlockVect(), MatBaseproduct(),
@@ -68,7 +69,7 @@ public:
}
bool useRightPreCond(){
- return true;
+ return false;
}
void preCondBlockVect(Block<Scalar,Primary>& input, Block<Scalar,Primary>& output){
for(int j=0 ; j<output.getSizeBlock() ; ++j){
@@ -85,7 +86,7 @@ public:
}
}
- void importMatrix(std::ifstream& file){
+ void importMatrixCmplx(std::ifstream& file){
if(file.is_open()){
bool commentLine = true;
std::string line;
@@ -125,6 +126,45 @@ public:
exit(0);
}
}
+ void importMatrixReal(std::ifstream& file){
+ if(file.is_open()){
+ bool commentLine = true;
+ std::string line;
+ while(commentLine){
+ getline(file,line);
+ std::cout<<line<<std::endl;
+ if(line[0] != '%'){
+ commentLine = false;
+ }
+ }
+ std::stringstream streamLine(line);
+ std::vector<std::string> headerLine;
+ std::string item;
+ while(std::getline(streamLine,item,' ')){
+ headerLine.push_back(item);
+ }
+ //Check square ?
+ if(headerLine[0] != headerLine[1]){
+ std::cout<<"line[0] != line[1] :: "<<line[0]<<" != "<<line[1]<<"\n";
+ std::cout<<"Matrix is not square !\nExiting\n";
+ exit(0);
+ }
+ this->dim = std::stoi(headerLine[0]);
+ this->data = new Scalar[dim*dim];
+ int nbValues = std::stoi(headerLine[2]);
+ headerLine.clear();
+ //Get the values
+ for(int i = 0; i<nbValues ; ++i){
+ std::string currLine[3];
+ file >> currLine[0] >> currLine[1] >> currLine[2];
+ int indx = (std::stoi(currLine[0])-1)*dim + std::stoi(currLine[1])-1;
+ data[indx] = std::stod(currLine[2]);
+ }
+ }else{
+ std::cout<<"File not open !\nExiting\n";
+ exit(0);
+ }
+ }
};
@@ -134,110 +174,198 @@ int main(int ac, char ** av){
for(auto ag : args)
std::cout << ag << "\n";
- using Primary=double;
- using Scalar=std::complex<double>;
+ int SizeBlock = 16;
+ if(ac>1){
+ SizeBlock = std::stoi(av[1]);
+ }
+ using Primary= double;
+ using Scalar = std::complex<double>;
+ using Matrix = InputMatrix<Scalar,Primary>;
//Random engine used to generate RHS
std::random_device rd;
std::mt19937 gen(/*rd()*/0);
std::uniform_real_distribution<> dis(0, 2);
- InputMatrix<Scalar,Primary> mat{};
- std::string location("../data/young1c.mtx");
+ Matrix mat{};
+ std::string location("../data/conf5_0_00l4x4_1000.mtx");
std::ifstream myfile(location,std::ios::in);
- mat.importMatrix(myfile);
+ mat.importMatrixCmplx(myfile);
int dim = mat.size();
- int SizeBlock = 15;
-
- int maxiteIB = dim/SizeBlock + (SizeBlock-1);
- Block<Scalar,Primary> XSol(SizeBlock,dim);
+ Block<Scalar,Primary> XExact(SizeBlock,dim);
Block<Scalar,Primary> RHS(SizeBlock,dim);
- Block<Scalar,Primary> X0(SizeBlock,dim);
- //Fill XSol
+
+ //Fill XExact
//Fill it!
for(int i=0 ; i<SizeBlock ; ++i){
- for(int j=0 ;j<dim; ++j){
- XSol.getPtr(i)[j] = {dis(gen),dis(gen)};
- X0.getPtr(i)[j] = {dis(gen),dis(gen)};
+ for(int j=0 ; j<XExact.getSizeBlock() ; ++j){
+ XExact.getPtr(i)[j] = dis(gen);
}
}
- mat.MatBlockVect(XSol,RHS);
+ mat.MatBlockVect(XExact,RHS);
- Block<Scalar,Primary> sol(SizeBlock,dim);
- Logger<Primary,5> logger;
+ int nbIte = (dim/SizeBlock)*2;
+ int nbIteIB = nbIte + SizeBlock-1;
+ int restart = mat.size();
- int nbIte = dim/SizeBlock;
- int nbIteIB = dim/SizeBlock+SizeBlock-1;
- int restart = nbIte;
+ using Arnoldi1 = Arnoldi<Matrix,Scalar,Primary>;
+ using Arnoldi2 = Arnoldi_IB<Matrix,Scalar,Primary>;
+ using Arnoldi3 = Arnoldi_QRInc<Matrix,Scalar,Primary>;
- int nb = BlockGMResQR<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,
- X0,nbIte,sol,
- restart,logger,
- 0.000001);
+ double elapsed1,elapsed2,elapsed3;
+ int nb1,nb2,nb3;
+ std::pair<Primary,Primary> minmax1,minmax2,minmax3;
+ // {
+ // std::cout<<"\n\tVersion Standard\n\n";
+ // Block<Scalar,Primary> X0(SizeBlock,dim);
+ // Block<Scalar,Primary> sol(SizeBlock,dim);
+ // Logger<Primary> log;
- // int nb = IBGMRes<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,X0,nbIteIB,
- // restart,sol,logger,
- // 0.0001);
- // int nb = BlockArnoldiIb<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,
- // X0,
- // maxiteIB,
- // sol,logger,
- // 0.0001);
+ // double tic = Logger<Primary>::GetTime();
- std::string locationRes("../data/res/young1c.res");
- std::ofstream fileRes(locationRes,std::ios::out);
- logger.forEachIte([&](int ite, int sizeB, Primary min, Primary max, double* time)
- {
- fileRes<<ite<<"\t"
- <<sizeB<<"\t"
- <<min<<"\t"
- <<max<<"\t"
- <<*time<<"\n";
- });
+ // nb1 = BGMRes<Arnoldi1,Matrix,Scalar,Primary>(mat,RHS,X0,nbIte,nbIte,sol,log,
+ // 0.00001);
+ // elapsed1 = Logger<Primary>::GetTime() - tic;
- std::cout<<"Return value / max_ite : "<<nb<<"/"<<maxiteIB<<"\n";
- //Check solution
+ // std::cout<<"Return value / max_ite : "<<nb1<<"/"<<nbIte<<"\n";
+ // std::string locationRes("../data/res/conf5_0_00l4x4_1000.res");
+ // std::ofstream fileRes(locationRes,std::ios::out);
+ // log.forEachIte([&](int ite, int nbMatVect, int currSize,
+ // Primary min, Primary max, double time)
+ // {
+ // fileRes<<ite<<"\t"
+ // <<nbMatVect<<"\t"
+ // <<currSize<<"\t"
+ // <<min<<"\t"
+ // <<max<<"\t"
+ // <<time<<"\n";
+ // });
+ // minmax1 = CompareBlocks<Scalar,Primary>(sol,XExact);
+ // }
{
- //Compute what i need for my metrics
- //MatxSol : Mat*Sol_computed
- Block<Scalar,Primary> MatxSol(SizeBlock,dim);
- mat.MatBlockVect(sol,MatxSol);
- //MatxSolMinusB : Mat*Sol_computed-B
- Block<Scalar,Primary> MatxSolMinusB(SizeBlock,dim);
- //X_Diff : XSol - Sol_computed
- Block<Scalar,Primary> X_Diff(SizeBlock,dim);
- MatxSolMinusB.CopyBlock(MatxSol);
- for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
- for(int j=0 ; j<dim; ++j){
- MatxSolMinusB.getPtr(idxRHS)[j] -= RHS.getPtr(idxRHS)[j];
- X_Diff.getPtr(idxRHS)[j] = XSol.getPtr(idxRHS)[j] -
- sol.getPtr(idxRHS)[j];
- }
- }
- //Then display for each vect
- for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
- std::cout<<"Metrics on "<<idxRHS<<"\n";
- Primary ForwardError = X_Diff.getNorm(idxRHS);
- std::cout <<"["<<idxRHS<<"]\t"
- <<"Forward Error (on X) \t" << ForwardError <<"\n";
- Primary SecondMemberError = MatxSolMinusB.getNorm(idxRHS);
- std::cout <<"["<<idxRHS<<"]\t"
- <<"Second Member Error \t" << SecondMemberError <<"\n";
- Primary ForwardErrorNormalized = ForwardError/(XSol.getNorm(idxRHS));
- std::cout <<"["<<idxRHS<<"]\t"
- <<"Forward Error Norm. \t" << ForwardErrorNormalized <<"\n";
- Primary normB = RHS.getNorm(idxRHS);
- Primary NormwiseBackwardError = SecondMemberError / normB;
- std::cout <<"["<<idxRHS<<"]\t"
- <<"Normwise Backward Error\t" << NormwiseBackwardError <<"\n\n";
-
- }
+ std::cout<<"\tVersion with Inexact Breakdown\n\n";
+ Block<Scalar,Primary> X0(SizeBlock,dim);
+ Block<Scalar,Primary> sol(SizeBlock,dim);
+ Logger<Primary> log;
+ double tic = Logger<Primary>::GetTime();
+
+ nb2 = BGMRes<Arnoldi2,Matrix,Scalar,Primary>(mat,RHS,X0,nbIteIB,restart,sol,log,
+ 0.00001);
+
+ elapsed2 = Logger<Primary>::GetTime() - tic;
+
+ std::cout<<"Return value / max_ite : "<<nb2<<"/"<<nbIteIB<<"\n";
+ std::string locationRes("../data/res/conf5_0_00l4x4_1000_IB.res");
+ std::ofstream fileRes(locationRes,std::ios::out);
+ log.forEachIte([&](int ite, int nbMatVect, int currSize,
+ Primary min, Primary max, double time)
+ {
+ fileRes<<ite<<"\t"
+ <<nbMatVect<<"\t"
+ <<currSize<<"\t"
+ <<min<<"\t"
+ <<max<<"\t"
+ <<time<<"\n";
+ });
+ minmax2 = CompareBlocks<Scalar,Primary>(sol,XExact);
}
+ {
+ std::cout<<"\tVersion with Optimized Least Square\n\n";
+ Block<Scalar,Primary> X0(SizeBlock,dim);
+ Block<Scalar,Primary> sol(SizeBlock,dim);
+ Logger<Primary> log;
+ double tic = Logger<Primary>::GetTime();
+
+ nb3 = BGMRes<Arnoldi3,Matrix,Scalar,Primary>(mat,RHS,X0,nbIte,restart,sol,log,
+ 0.00001);
+
+ elapsed3 = Logger<Primary>::GetTime() - tic;
+
+ std::cout<<"Return value / max_ite : "<<nb3<<"/"<<nbIte<<"\n";
+ std::string locationRes("../data/res/sherman4_QR.res");
+ std::ofstream fileRes(locationRes,std::ios::out);
+ log.forEachIte([&](int ite, int nbMatVect, int currSize,
+ Primary min, Primary max, double time)
+ {
+ fileRes<<ite<<"\t"
+ <<nbMatVect<<"\t"
+ <<currSize<<"\t"
+ <<min<<"\t"
+ <<max<<"\t"
+ <<time<<"\n";
+ });
+ minmax3 = CompareBlocks<Scalar,Primary>(sol,XExact);
+ }
+
+
+
+ // std::string locationRes("../data/res/young1c.res");
+ // std::ofstream fileRes(locationRes,std::ios::out);
+ // logger.forEachIte([&](int ite, int sizeB, Primary min, Primary max, double* time)
+ // {
+ // fileRes<<ite<<"\t"
+ // <<sizeB<<"\t"
+ // <<min<<"\t"
+ // <<max<<"\t"
+ // <<*time<<"\n";
+ // });
+
+
+ //Check solution
+ // {
+ // //Compute what i need for my metrics
+ // //MatxSol : Mat*Sol_computed
+ // Block<Scalar,Primary> MatxSol(SizeBlock,dim);
+ // mat.MatBlockVect(sol,MatxSol);
+ // //MatxSolMinusB : Mat*Sol_computed-B
+ // Block<Scalar,Primary> MatxSolMinusB(SizeBlock,dim);
+ // //X_Diff : XExact - Sol_computed
+ // Block<Scalar,Primary> X_Diff(SizeBlock,dim);
+ // MatxSolMinusB.CopyBlock(MatxSol);
+ // for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
+ // for(int j=0 ; j<dim; ++j){
+ // MatxSolMinusB.getPtr(idxRHS)[j] -= RHS.getPtr(idxRHS)[j];
+ // X_Diff.getPtr(idxRHS)[j] = XExact.getPtr(idxRHS)[j] -
+ // sol.getPtr(idxRHS)[j];
+ // }
+ // }
+ // //Then display for each vect
+ // for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
+ // std::cout<<"Metrics on "<<idxRHS<<"\n";
+ // Primary ForwardError = X_Diff.getNorm(idxRHS);
+ // std::cout <<"["<<idxRHS<<"]\t"
+ // <<"Forward Error (on X) \t" << ForwardError <<"\n";
+ // Primary SecondMemberError = MatxSolMinusB.getNorm(idxRHS);
+ // std::cout <<"["<<idxRHS<<"]\t"
+ // <<"Second Member Error \t" << SecondMemberError <<"\n";
+ // Primary ForwardErrorNormalized = ForwardError/(XExact.getNorm(idxRHS));
+ // std::cout <<"["<<idxRHS<<"]\t"
+ // <<"Forward Error Norm. \t" << ForwardErrorNormalized <<"\n";
+ // Primary normB = RHS.getNorm(idxRHS);
+ // Primary NormwiseBackwardError = SecondMemberError / normB;
+ // std::cout <<"["<<idxRHS<<"]\t"
+ // <<"Normwise Backward Error\t" << NormwiseBackwardError <<"\n\n";
+
+ // }
+ // }
+
+ //Display timers
+ std::cout<<"Elapsed time for each method::\n";
+ std::cout<<"\tElapsed time for Standard Version\t\t"<<elapsed1
+ <<" seconds for "<<nb1<<"/"<<nbIte<<" iterations\n"
+ <<"Forward error : Min\t"<<minmax1.first<<" Max\t"<<minmax1.second<<"\n";
+ std::cout<<"\tElapsed time for Inexact Breakdown\t\t"<<elapsed2
+ <<" seconds for "<<nb2<<"/"<<nbIteIB<<" iterations\n"
+ <<"Forward error : Min\t"<<minmax2.first<<" Max\t"<<minmax2.second<<"\n";
+ std::cout<<"\tElapsed time for Optimized Least Square\t\t"<<elapsed3
+ <<" seconds for "<<nb3<<"/"<<nbIte<<" iterations\n"
+ <<"Forward error : Min\t"<<minmax3.first<<" Max\t"<<minmax3.second<<"\n";
+
return 0;
diff --git a/test/testQRIncrementalGMRes.cpp b/test/testQRIncrementalGMRes.cpp
deleted file mode 100644
index 28014938e24d4995852a17a78a79799d2f62df8d..0000000000000000000000000000000000000000
--- a/test/testQRIncrementalGMRes.cpp
+++ /dev/null
@@ -1,186 +0,0 @@
-#include <vector>
-#include <complex>
-#include <random>
-#include "../src/LapackInterface.hpp"
-#include "../src/numerical.hpp"
-#include "../src/Block.hpp"
-#include "../src/BlockGMResQR.hpp"
-#include "../src/BlockGMRes.hpp"
-
-/**
- * Matrix Class : need to implement MatBlockVect(), MatBaseproduct(),
- * useRightPreCond(), preCondBlockVect(),preCondBaseProduct() and
- * size() member function in order to works with BlockGMRes algorithm
- *
- * An abstract class could be used to enforce that.
- */
-template<class Scalar,class Primary>
-class InputMatrix{
- int dim;
- Scalar * data;
-public:
- InputMatrix(int inSize) : dim(inSize),data(nullptr){
- data = new Scalar[dim*dim];
- }
-
- ~InputMatrix(){
- delete [] data;
- }
-
- int size(){
- return dim;
- }
-
- void initMatrix(std::random_device& rd){
- std::mt19937 gen(/*rd()*/0);
- std::uniform_real_distribution<> dis(0, 2);
- for(int i=0 ; i<dim ; ++i){
- for(int j=0 ; j<dim ; ++j){
- data[dim*i+j] = {dis(gen),dis(gen)};
- }
- }
- // Scalar diag = dim;
- // Scalar no_diag = 1;
- // for(int i=0 ; i<dim ; ++i){
- // for(int j=0 ; j<dim ; ++j){
- // if(i==j)
- // data[i*dim + j] = diag+j;
- // else
- // data[i*dim + j] = no_diag;
- // }
- // }
- }
-
- void MatBlockVect(Block<Scalar,Primary>& in, Block<Scalar,Primary>& out){
- call_gemm<Scalar>(dim,in.getSizeBlock(),dim,
- data,dim,
- in.getPtr(), in.getLeadingDim(),
- out.getPtr(),out.getLeadingDim(),
- Scalar(1),Scalar(0));
- }
-
- void MatBaseProduct(Scalar * toRead,int nbRHS, Scalar * ptrToWrite){
- call_gemm<Scalar>(dim,nbRHS, dim,
- data,dim,
- toRead,dim,
- ptrToWrite,dim,
- Scalar(1),Scalar(0));
- }
-
- bool useRightPreCond(){
- return false;
- }
- void preCondBlockVect(Block<Scalar,Primary>& input, Block<Scalar,Primary>& output){
- // for(int j=0 ; j<output.getSizeBlock() ; ++j){
- // for(int i=0 ; i<dim ; ++i){
- // output.getPtr(j)[i] = input.getPtr(j)[i]*(1/data[i*dim+i]);
- // }
- // }
- }
- void preCondBaseProduct(Scalar * ptr,Block<Scalar,Primary>& output){
- // for(int j=0 ; j<output.getSizeBlock() ; ++j){
- // for(int i=0 ; i<dim ; ++i){
- // output.getPtr(j)[i] = ptr[j*dim+i]*1/(data[i*dim+i]);
- // }
- // }
- }
-};
-
-
-int main(int ac, char ** av){
- std::vector<std::string > args(av,av+ac);
- for(auto ag : args)
- std::cout << ag << "\n";
-
- int SizeBlock=3,nbBlock=4,restart=30;
-
- if (args.size()>1){
- SizeBlock = std::stoi(args[1]);
- }
- if (args.size()>2){
- nbBlock = std::stoi(args[2]);
- }
- restart = nbBlock;
- int dim = SizeBlock*nbBlock;
-
- //Random engine used to generate matrix
- std::random_device rd;
- std::mt19937 gen(/*rd()*/0);
- std::uniform_real_distribution<> dis(0, 2);
-
- using Primary = double;
- using Scalar = std::complex<Primary>;
-
- //Create matrix
- InputMatrix<Scalar,Primary> mat(dim);
- mat.initMatrix(rd);
-
- //Create X_Exact
- Block<Scalar,Primary> X_Exact(SizeBlock,dim);
- //Fill it!
- for(int i=0 ; i<SizeBlock ; ++i){
- for(int j=0 ;j<dim; ++j){
- X_Exact.getPtr(i)[j] = dis(gen);
- }
- }
-
- //Compute RHS as Mat*X_Exact
- Block<Scalar,Primary> RHS(SizeBlock,dim);
- mat.MatBlockVect(X_Exact,RHS);
-
- //Compute a starting Block
- Block<Scalar,Primary> X0(SizeBlock,dim);
-
- //Fill it!
- for(int i=0 ; i<SizeBlock ; ++i){
- for(int j=0 ;j<dim; ++j){
- X0.getPtr(i)[j] = {dis(gen),dis(gen)};
- }
- }
- Block<Scalar,Primary> sol(SizeBlock,dim);
- Logger<Primary,5> log;
- BlockGMResQR<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,
- X0,nbBlock,sol,
- restart,log,
- 0.000001);
-
- {
- //Compute what i need for my metrics
- //MatxSol : Mat*Sol_computed
- Block<Scalar,Primary> MatxSol(SizeBlock,dim);
- mat.MatBlockVect(sol,MatxSol);
- //MatxSolMinusB : Mat*Sol_computed-B
- Block<Scalar,Primary> MatxSolMinusB(SizeBlock,dim);
- //X_Diff : X_Exact - Sol_computed
- Block<Scalar,Primary> X_Diff(SizeBlock,dim);
- MatxSolMinusB.CopyBlock(MatxSol);
- for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
- for(int j=0 ; j<dim; ++j){
- MatxSolMinusB.getPtr(idxRHS)[j] -= RHS.getPtr(idxRHS)[j];
- X_Diff.getPtr(idxRHS)[j] = X_Exact.getPtr(idxRHS)[j] -
- sol.getPtr(idxRHS)[j];
- }
- }
- //Then display for each vect
- for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
- std::cout<<"Metrics on "<<idxRHS<<"\n";
- Primary ForwardError = X_Diff.getNorm(idxRHS);
- std::cout <<"["<<idxRHS<<"]\t"
- <<"Forward Error (on X) \t" << ForwardError <<"\n";
- Primary SecondMemberError = MatxSolMinusB.getNorm(idxRHS);
- std::cout <<"["<<idxRHS<<"]\t"
- <<"Second Member Error \t" << SecondMemberError <<"\n";
- Primary ForwardErrorNormalized = ForwardError/(X_Exact.getNorm(idxRHS));
- std::cout <<"["<<idxRHS<<"]\t"
- <<"Forward Error Norm. \t" << ForwardErrorNormalized <<"\n";
- Primary normB = RHS.getNorm(idxRHS);
- Primary NormwiseBackwardError = SecondMemberError / normB;
- std::cout <<"["<<idxRHS<<"]\t"
- <<"Normwise Backward Error\t" << NormwiseBackwardError <<"\n\n";
-
- }
- }
-
-
-
-}