Skip to content

GitLab

fabulous
fabulous
Commit 49e17ea6 authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille
Browse files
Options

Cleaning on its way : tests updated (removed old one, and update others, while... 

Cleaning on its way : tests updated (removed old one, and update others, while maintaining test/CMakelists right)
parent 9c4e5ef6
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 300 additions and 506 deletions
test/CMakeLists.txt
View file @ 49e17ea6
......@@ -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
)
......
test/testBlockGMRes.cpp deleted 100644 → 0
View file @ 9c4e5ef6
#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;
}
test/testMatrixIW.cpp
View file @ 49e17ea6
......@@ -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;
}
test/testMatrixMarket.cpp
View file @ 49e17ea6
This diff is collapsed.
test/testQRIncrementalGMRes.cpp deleted 100644 → 0
View file @ 9c4e5ef6
#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";
}
}
}
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