Commit 8ef263f3 authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

Innovation Work Test code updated, display now for the three methods the time spent.

parent 06e93feb
......@@ -34,7 +34,7 @@ struct Arnoldi_QRInc{
//Storage needed for computating
Base<Scalar,Primary> base(SizeBlock,max_iter+1,dim);
HessQR<Scalar> Hess(max_iter+1,SizeBlock);
HessQR<Scalar,Primary> Hess(max_iter+1,SizeBlock);
Block<Scalar,Primary> Yj(SizeBlock,dim+SizeBlock);
Block<Scalar,Primary> R0(SizeBlock,dim);
......
......@@ -19,7 +19,7 @@
*
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<typename Arn,class Matrix,class Scalar,class Primary=Scalar>
template<typename Arn,class Matrix,class Scalar,class Primary>
int BGMRes(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
int max_iter, int restart,
Block<Scalar,Primary>& Sol,Primary epsilon = 0.0000001){
......
......@@ -2,14 +2,14 @@
#include <iostream>
#include <fstream>
#include <stdint.h>
#include "../src/Block.hpp"
#include "../src/BlockGMRes.hpp"
#include "../src/LapackInterface.hpp"
#include "../src/Logger.hpp"
#include "../src/BlockArnoldiIb.hpp"
#include "../src/IBGMRes.hpp"
#include "../src/Arnoldi.hpp"
#include "../src/Arnoldi_IB.hpp"
#include "../src/Arnoldi_QRInc.hpp"
#include "../src/BGMRes.hpp"
/**
* Driver to read Innovation Works arrays
*/
......@@ -133,7 +133,8 @@ int main(int ac, char ** av){
std::cout << ag << "\n";
using Primary=double;
using Scalar=std::complex<double>;
using Scalar =std::complex<double>;
using Matrix = InputMatrix<Scalar,Primary>;
//Random engine used to generate RHS
std::random_device rd;
......@@ -141,7 +142,7 @@ int main(int ac, char ** av){
std::uniform_real_distribution<> dis(0, 2);
InputMatrix<Scalar,Primary> mat{};
Matrix mat{};
{
std::string location("../data/matrices_BEM/MatconeSpherePC_MAIN_MAIN_0");
std::ifstream myfile(location,std::ios::binary);
......@@ -174,77 +175,90 @@ int main(int ac, char ** av){
int dim = mat.size();
int SizeBlock = RHS.getSizeBlock();
int max_ite = dim/SizeBlock;
int maxiteIB = dim / SizeBlock + (SizeBlock);
Block<Scalar,Primary> X0(SizeBlock,dim);
// 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> logger;
// int nb = BlockArnoldiIb<InputMatrix<Scalar,Primary>,
// Scalar,Primary,5>(mat,RHS,
// X0,
// maxiteIB,
// sol,logger,
// 0.00001);
int nb = IBGMRes<InputMatrix<Scalar,Primary>,Scalar,Primary,5>(mat,RHS,X0,
3,maxiteIB,sol,
logger,0.00001);
std::string locationRes("../data/res/MatconeSpherePC_MAIN_MAIN_0.res");
std::ofstream fileRes(locationRes,std::ios::out);
fileRes<<"Ite\tSizeBlck\tMinRes\tMaxRes\tTime\tMatVect\tOrtho\tResidual\tSVD\n";
logger.forEachIte([&](int ite, int sizeB, Primary min, Primary max, double* time)
{
fileRes<<ite<<"\t"
<<sizeB<<"\t"
<<min<<"\t"
<<max<<"\t"
<<*time<<"\t";
for(int i=1 ; i<5 ; ++i){
fileRes<<time[i]<<"\t";
}
fileRes<<"\n";
});
std::cout<<"Return value / max_ite : "<<nb<<"/"<<maxiteIB<<"\n";
//Check solution
using Arnoldi1 = Arnoldi<Matrix,Scalar,Primary>;
using Arnoldi2 = Arnoldi_IB<Matrix,Scalar,Primary>;
using Arnoldi3 = Arnoldi_QRInc<Matrix,Scalar,Primary>;
double elapsed1,elapsed2,elapsed3;
{
//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<<"\n\tVersion Standard\n\n";
Block<Scalar,Primary> X0(SizeBlock,dim);
Block<Scalar,Primary> sol(SizeBlock,dim);
double tic = Logger<Primary>::GetTime();
int nb = BGMRes<Arnoldi1,Matrix,Scalar,Primary>(mat,RHS,X0,max_ite,max_ite,sol,
0.000001);
elapsed1 = Logger<Primary>::GetTime() - tic;
std::cout<<"Return value / max_ite : "<<nb<<"/"<<max_ite<<"\n";
}
{
std::cout<<"\tVersion with Inexact Breakdown\n\n";
Block<Scalar,Primary> X0(SizeBlock,dim);
Block<Scalar,Primary> sol(SizeBlock,dim);
double tic = Logger<Primary>::GetTime();
int nb = BGMRes<Arnoldi2,Matrix,Scalar,Primary>(mat,RHS,X0,maxiteIB,max_ite,sol,
0.000001);
elapsed2 = Logger<Primary>::GetTime() - tic;
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);
double tic = Logger<Primary>::GetTime();
int nb = BGMRes<Arnoldi3,Matrix,Scalar,Primary>(mat,RHS,X0,max_ite,max_ite,sol,
0.000001);
elapsed3 = Logger<Primary>::GetTime() - tic;
std::cout<<"Return value / max_ite : "<<nb<<"/"<<max_ite<<"\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";
// }
// }
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