Commit bee2c8ea authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

Old versions of IB, Standard and QR algo suppressed

parent ed50c6fa
#ifndef BLOCKGMRES_HPP
#define BLOCKGMRES_HPP
#include "LapackInterface.hpp"
#include "HessStandard.hpp"
#include "Base.hpp"
#include "Block.hpp"
#include "Logger.hpp"
/**
* @brief Computation of the GMRES Block.
*
* Ref : Saad P.241
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<class Matrix,typename Scalar,typename Primary=Scalar,int N=5>
int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
//Matrix preCond Deal with that later
int max_iter, Block<Scalar,Primary>& Sol, int nb_iter_restart,
Logger<Primary,5>& log,
Primary epsilon = 0.00001){
if(A.useRightPreCond()){
std::cout<<"Right Pre Conditionner used\n";
}else{
std::cout<<"NO - Right Pre Conditionner used\n";
}
//Global variables
int dim = A.size();
int SizeBlock = B.getSizeBlock();
//store norms
std::vector<Primary> normB;
for(int i=0 ; i< SizeBlock ; ++i){
normB.push_back(B.getNorm(i));
}
//Residual vectors needed for computing the solution
Block<Scalar,Primary> Yj(SizeBlock,dim+SizeBlock);
//Add a loop for restart
int k=0,nbRestart=0, globalStep=0;
log.init_log(SizeBlock);
while(k<max_iter){
std::cout<<"\t\t\tRestarting number "<<nbRestart<<" ("<<nb_iter_restart<<")\n";
// X0.displayBlock("New X0");
Block<Scalar,Primary> R0(SizeBlock,dim);
//Fill matrix R0 with each vector r0_i = B[i] - A * X0[i]
A.MatBlockVect(X0,R0); //R0 = A*X0
for(int i=0 ; i<SizeBlock ; ++i){ //R0 = B - R0
for(int j=0 ; j<dim ; ++j){
R0.getPtr(i)[j] = B.getPtr(i)[j] - R0.getPtr(i)[j];
}
}
//Storage needed for computating
Base<Scalar,Primary> base(SizeBlock,nb_iter_restart+1,dim);
Hessenberg<Scalar> Hess(nb_iter_restart,SizeBlock);
//V1 will be first block of base, and R1, will be used for least
//square pb, in criteria verification
Block<Scalar,Primary> V1(SizeBlock,dim);
Block<Scalar,Primary> R1(SizeBlock,SizeBlock);
//Q-R Decomposition of (R0)
R0.ComputeQRFacto(V1,R1);
//Fisrt block of base added
base.addBlockDatas(V1);
std::cout<<"Starting main loop"<<std::endl;
for(int j=0 ; (j<nb_iter_restart) && (j+k)<= max_iter ; ++j){
//Compute block Wj
Block<Scalar,Primary> Wj(SizeBlock,dim);
if(A.useRightPreCond()){//Wj = A*preCond*V[j] if preCond
Block<Scalar,Primary> temp(SizeBlock,dim);
A.preCondBaseProduct(base.getPtr(base.getCurrentBlock()-1),temp);
A.MatBlockVect(temp,Wj);
}else{//Wj = A*V[j] without preCond
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),
Wj.getSizeBlock(),Wj.getPtr());
}
//Where to write in hess
Scalar * toWrite = Hess.getPtr(j);
base.ComputeMGSOrtho(Wj,toWrite,Hess.getLeadingDim());
std::cout<<"Ortho done"<<std::endl;
//QR facto of orthogonalized Wj
Block<Scalar,Primary> Q(SizeBlock,dim);
Block<Scalar,Primary> R(SizeBlock,SizeBlock);
Wj.ComputeQRFacto(Q,R);
//Copy Block R inside Hess
Hess.addRPartToHess(R);
std::cout<<"QR done"<<std::endl;
//Tests over criterion
Hess.ComputeBlockResidual(R1,Yj);
//Here we compute Base*Y +X0 and display res
// {
// //if no precond
// // X0 = X0 + Base*Yj
// Block<Scalar,Primary> Resultat{SizeBlock,dim};
// //Yj.displayBlock("Yj at current iter");
// base.ComputeProduct(Yj,Resultat);
// for(int i=0 ; i<SizeBlock ; ++i){
// for(int j=0 ; j<dim ; ++j){
// Resultat.getPtr(i)[j] += X0.getPtr(i)[j];
// }
// }
// std::cout<<"#########\t"<<j<<"\t#########\n";
// //Resultat.displayBlock("Current Xm");
// //Then we compute real residual
// Block<Scalar,Primary> AxSol(SizeBlock,dim);
// A.MatBlockVect(Resultat,AxSol);
// for(int i=0 ; i<SizeBlock ; ++i){
// for(int j=0 ; j<dim ; ++j){
// AxSol.getPtr(i)[j] -= B.getPtr(i)[j];
// }
// }
// for(int i=0 ; i<SizeBlock ; ++i){
// std::cout<<"PREFIX "<<i<<"\t"<<j<<"\t"<<AxSol.getNorm(i)<<"\n";
// }
// }
std::cout<<"["<<j<<"]\t"<<"Residual computed"<<std::endl;
{//Display residual to see convergence
std::pair<Primary,Primary> MinMax
= Hess.template displayBlockResidual<Primary>(Yj,R1,normB);
log.add_ite(j,MinMax.first,MinMax.second);
std::cout<<"["<<j<<"]\t"<<"Min "<<MinMax.first
<<"\tMax "<<MinMax.second<<"\n";
if(MinMax.second <= epsilon){
//Compare to the real residual
std::pair<Primary,Primary> RealMinMax
= std::make_pair<Primary,Primary>(100,0);
//Generate block Sol = Base*Yj
//if Precond : block Sol = preCond*Base*yj
Block<Scalar,Primary> blockSol(SizeBlock,dim);
if(A.useRightPreCond()){
Block<Scalar,Primary> temp(SizeBlock,dim);
base.ComputeProduct(Yj,temp);
A.preCondBlockVect(temp,blockSol);
}else{
base.ComputeProduct(Yj,blockSol);
}
//Compute addition with X0
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ; j<dim ; ++j){
blockSol.getPtr(i)[j] += X0.getPtr(i)[j];
}
}
//Compute A*sol computed
Block<Scalar,Primary> AxSol(SizeBlock,dim);
A.MatBlockVect(blockSol,AxSol);
//Compute real residual
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ; j<dim ; ++j){
AxSol.getPtr(i)[j] -= B.getPtr(i)[j];
}
}
//For each vector, compute norm
for(int i = 0 ; i<SizeBlock ; ++i){
auto norm = AxSol.getNorm(i)/normB[i];
if(norm > RealMinMax.second){
RealMinMax.second = norm;
}
if(norm < RealMinMax.first){
RealMinMax.first = norm;
}
}//If convergence
if(RealMinMax.second <= epsilon){
std::cout<<"Convergence !\t Step : "<<j+k<<std::endl;
std::cout<<"Min "<<RealMinMax.first
<<"\tMax "<<RealMinMax.second<<"\n";
//Modify k in order to stop all iterations
k = max_iter;
//Base augmentation
base.addBlockDatas(Q);
break;
}else{
std::cout
<<"Residual is small enough, but not the real's one \nMin2"
<<RealMinMax.first<<"\tMax2"<<RealMinMax.second<<"\n";
}
}
}
//Base augmentation
base.addBlockDatas(Q);
std::cout<<"Base has been augmented\n";
globalStep++;
}
k += nb_iter_restart;
nbRestart += 1;
//Compute new starting block X0
//if no precond
// X0 = X0 + Base*Yj
//else
// X0 = X0 + preCond*Base*yj
Block<Scalar,Primary> newX0(SizeBlock,dim);
if(A.useRightPreCond()){
Block<Scalar,Primary> temp(SizeBlock,dim);
base.ComputeProduct(Yj,temp);
A.preCondBlockVect(temp,newX0);
}else{
base.ComputeProduct(Yj,newX0);
}
//Copy back
//End of use of Yj, reset it
Yj.reset();
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ; j<dim ; ++j){
X0.getPtr(i)[j] = X0.getPtr(i)[j] + newX0.getPtr(i)[j];
}
}
}
Sol.CopyBlock(X0);
return globalStep;
}
#endif //BLOCKGMRES_HPP
#ifndef BLOCKGMRESQR_HPP
#define BLOCKGMRESQR_HPP
#include "LapackInterface.hpp"
#include "HessQR.hpp"
#include "Base.hpp"
#include "Block.hpp"
#include "Logger.hpp"
/**
* @brief Computation of the GMRES Block.
*
* Ref : Saad P.241
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<class Matrix,typename Scalar,typename Primary=Scalar,int N=5>
int BlockGMResQR(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
//Matrix preCond Deal with that later
int max_iter, Block<Scalar,Primary>& Sol, int nb_iter_restart,
Logger<Primary,5>& log,
Primary epsilon = 0.00001){
if(A.useRightPreCond()){
std::cout<<"Right Pre Conditionner used\n";
}else{
std::cout<<"NO - Right Pre Conditionner used\n";
}
//Global variables
int dim = A.size();
int SizeBlock = B.getSizeBlock();
//store norms
std::vector<Primary> normB;
for(int i=0 ; i< SizeBlock ; ++i){
normB.push_back(B.getNorm(i));
}
//Residual vectors needed for computing the solution
Block<Scalar,Primary> Yj(SizeBlock,dim+SizeBlock);
//Add a loop for restart
int k=0,nbRestart=0, globalStep=0;
log.init_log(SizeBlock);
while(k<max_iter){
std::cout<<"\t\t\tRestarting number "<<nbRestart<<" ("<<nb_iter_restart<<")\n";
// X0.displayBlock("New X0");
Block<Scalar,Primary> R0(SizeBlock,dim);
//Fill matrix R0 with each vector r0_i = B[i] - A * X0[i]
A.MatBlockVect(X0,R0); //R0 = A*X0
for(int i=0 ; i<SizeBlock ; ++i){ //R0 = B - R0
for(int j=0 ; j<dim ; ++j){
R0.getPtr(i)[j] = B.getPtr(i)[j] - R0.getPtr(i)[j];
}
}
//Storage needed for computating
Base<Scalar,Primary> base(SizeBlock,nb_iter_restart+1,dim);
Hessenberg_2<Scalar,Primary> Hess(nb_iter_restart+1,SizeBlock);
//V1 will be first block of base, and R1, will be used for least
//square pb, in criteria verification
Block<Scalar,Primary> V1(SizeBlock,dim);
Block<Scalar,Primary> R1(SizeBlock,SizeBlock);
//Q-R Decomposition of (R0)
R0.ComputeQRFacto(V1,R1);
//Fisrt block of base added
base.addBlockDatas(V1);
//Give to Hessenberg Gamma
Hess.initRHS(R1);
std::cout<<"Starting main loop"<<std::endl;
for(int j=0 ; (j<nb_iter_restart) && (j+k)<= max_iter ; ++j){
//Compute block Wj
Block<Scalar,Primary> Wj(SizeBlock,dim);
if(A.useRightPreCond()){//Wj = A*preCond*V[j] if preCond
Block<Scalar,Primary> temp(SizeBlock,dim);
A.preCondBaseProduct(base.getPtr(base.getCurrentBlock()-1),temp);
A.MatBlockVect(temp,Wj);
}else{//Wj = A*V[j] without preCond
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),
Wj.getSizeBlock(),Wj.getPtr());
}
//Where to write in hess
Scalar * toWrite = Hess.getNewCol();
base.ComputeMGSOrtho(Wj,toWrite,Hess.getLeadingDim());
std::cout<<"Ortho done"<<std::endl;
//QR facto of orthogonalized Wj
Block<Scalar,Primary> Q(SizeBlock,dim);
Block<Scalar,Primary> R(SizeBlock,SizeBlock);
Wj.ComputeQRFacto(Q,R);
//Copy Block R inside Hess
Hess.addRPartToHess(R);
std::cout<<"QR done"<<std::endl;
//Tests over criterion
std::pair<Primary,Primary> MinMax = Hess.LeastSquare(Yj,normB);
std::cout<<"["<<j<<"]\t"<<"Residual computed"<<std::endl;
//Here we compute Base*Y +X0 and display res
// {
// //if no precond
// // X0 = X0 + Base*Yj
// Block<Scalar,Primary> Resultat{SizeBlock,dim};
// //Yj.displayBlock("Yj at current iter");
// base.ComputeProduct(Yj,Resultat);
// for(int i=0 ; i<SizeBlock ; ++i){
// for(int j=0 ; j<dim ; ++j){
// Resultat.getPtr(i)[j] += X0.getPtr(i)[j];
// }
// }
// std::cout<<"#########\t"<<j<<"\t#########\n";
// //Resultat.displayBlock("Current Xm");
// //Then we compute real residual
// Block<Scalar,Primary> AxSol(SizeBlock,dim);
// A.MatBlockVect(Resultat,AxSol);
// for(int i=0 ; i<SizeBlock ; ++i){
// for(int j=0 ; j<dim ; ++j){
// AxSol.getPtr(i)[j] -= B.getPtr(i)[j];
// }
// }
// for(int i=0 ; i<SizeBlock ; ++i){
// std::cout<<"PREFIX "<<i<<"\t"<<j<<"\t"<<AxSol.getNorm(i)<<"\n";
// }
// }
{//Display residual to see convergence
log.add_ite(j,MinMax.first,MinMax.second);
std::cout<<"["<<j<<"]\t"<<"Min "<<MinMax.first
<<"\tMax "<<MinMax.second<<"\n";
if(MinMax.second <= epsilon){
//Compare to the real residual
std::pair<Primary,Primary> RealMinMax
= std::make_pair<Primary,Primary>(100,0);
//Generate block Sol = Base*Yj
//if Precond : block Sol = preCond*Base*yj
Block<Scalar,Primary> blockSol(SizeBlock,dim);
if(A.useRightPreCond()){
Block<Scalar,Primary> temp(SizeBlock,dim);
base.ComputeProduct(Yj,temp);
A.preCondBlockVect(temp,blockSol);
}else{
base.ComputeProduct(Yj,blockSol);
}
//Compute addition with X0
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ; j<dim ; ++j){
blockSol.getPtr(i)[j] += X0.getPtr(i)[j];
}
}
//Compute A*sol computed
Block<Scalar,Primary> AxSol(SizeBlock,dim);
A.MatBlockVect(blockSol,AxSol);
//Compute real residual
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ; j<dim ; ++j){
AxSol.getPtr(i)[j] -= B.getPtr(i)[j];
}
}
//For each vector, compute norm
for(int i = 0 ; i<SizeBlock ; ++i){
auto norm = AxSol.getNorm(i)/normB[i];
if(norm > RealMinMax.second){
RealMinMax.second = norm;
}
if(norm < RealMinMax.first){
RealMinMax.first = norm;
}
}//If convergence
if(RealMinMax.second <= epsilon){
std::cout<<"Convergence !\t Step : "<<j+k<<std::endl;
std::cout<<"Min "<<RealMinMax.first
<<"\tMax "<<RealMinMax.second<<"\n";
//Modify k in order to stop all iterations
k = max_iter;
//Base augmentation
base.addBlockDatas(Q);
break;
}else{
std::cout
<<"Residual is small enough, but not the real's one \nMin2"
<<RealMinMax.first<<"\tMax2"<<RealMinMax.second<<"\n";
}
}
}
//Base augmentation
base.addBlockDatas(Q);
std::cout<<"Base has been augmented\n";
globalStep++;
//Hess.displayHess("Hess at the end of the loop");
}
k += nb_iter_restart;
nbRestart += 1;
//Compute new starting block X0
//if no precond
// X0 = X0 + Base*Yj
//else
// X0 = X0 + preCond*Base*yj
Block<Scalar,Primary> newX0(SizeBlock,dim);
if(A.useRightPreCond()){
Block<Scalar,Primary> temp(SizeBlock,dim);
base.ComputeProduct(Yj,temp);
A.preCondBlockVect(temp,newX0);
}else{
base.ComputeProduct(Yj,newX0);
}
//Copy back
//Yj.displayBlock("Yj at last iter before restart");
//End of use of Yj, reset it
Yj.reset();
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ; j<dim ; ++j){
X0.getPtr(i)[j] = X0.getPtr(i)[j] + newX0.getPtr(i)[j];
}
}
}
Sol.CopyBlock(X0);
return globalStep;
}
#endif //BLOCKGMRESQR_HPP
#ifndef IBGMRES_HPP
#define IBGMRES_HPP
#include "BlockWP.hpp"
#include "HessExtended.hpp"
#include "Logger.hpp"
#include "IBArnoldi.hpp"
/**
* @brief : GMRes with Inexact Breakdowns, restart and preconditionner
*
* @param A : Matrix, need to provide MatBlockVect(), size(),
* MatBaseProduct(), useRightPreCond(), preCondBlockVect(), preCondBaseProduct()
*
* @param B : Block containing Right Hand Side
* @param X0 : Initial guess
* @param max_iter : Total iteration
* @param restart : max number of iteration before restarting
* @param Sol : Block to store the solution
* @param log : instance of Logger class, to store information
* gathered at each step.
* @param epsilon : Target accuracy
*
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<class Matrix,typename Scalar,typename Primary=Scalar,int N=5>
int IBGMRes(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
int max_iter, int restart,
Block<Scalar,Primary>& Sol,Logger<Primary,N>& log,
Primary epsilon = 0.0000001){
if(A.useRightPreCond()){
std::cout<<"\nRight Pre Conditionner used\n\n";
}else{
std::cout<<"\nNO - Right Pre Conditionner used\n\n";
}
//Cste Size
int SizeBlock = B.getSizeBlock();
int dim = A.size();
std::cout<<"#######################################################\n";
std::cout<<"################## I.B. Block GMRes ###################\n";
std::cout<<"#######################################################\n";
std::cout<<"Dimension of Problem :\t"<<dim<<"\n";
std::cout<<"Number of Right Hand Side :\t"<<SizeBlock<<"\n";
std::cout<<"Maximum nb of Iterations:\t"<<max_iter<<"\n";
std::cout<<"Restart IBGMRes each:\t"<<restart<<" iterations.\n";
std::cout<<"#######################################################\n";
//Norm of B
std::vector<Primary> normB;
for(int i=0 ; i< SizeBlock ; ++i){
normB.push_back(B.getNorm(i));
}
//Compute epsilon_R (step 1)
Primary epsilon_R = normB.front();
for(auto &a : normB){
if(a<epsilon_R){
epsilon_R = a;
}
}
epsilon_R *= epsilon;
//Compute first block residual (step 1)
Block<Scalar,Primary> R0(SizeBlock,dim);
//Fill matrix R0 with each vector r0_i = B[i] - A * X0[i]
A.MatBlockVect(X0,R0); //R0 = A*X0
for(int i=0 ; i<SizeBlock ; ++i){ //R0 = B - R0
for(int k=0 ; k<dim ; ++k){
R0.getPtr(i)[k] = B.getPtr(i)[k] - R0.getPtr(i)[k];
}
}
Block<Scalar,Primary> V1(SizeBlock,dim);
Block<Scalar,Primary> R1(SizeBlock,SizeBlock);
//Q-R Decomposition of (R0)
R0.ComputeQRFacto(V1,R1);
//Set the base
Base<Scalar,Primary> base(SizeBlock,restart+1,dim);
base.addBlockDatas(V1);
//Set the hessenberg
HessExtended<Scalar> L(SizeBlock,restart);
// Global sstep is total number of iteration (thus is <max_iter)
int k=0,nb_restart=0;
while(k<max_iter){
//Compute nb_iter to give to Arnoldi procedure
int arnoldi_ite = std::min(restart,max_iter-k);
std::cout<<"\n################## Restart number "<<nb_restart
<<" #################\n";
//Call Block Arnoldi's procedure with IB
int nb = IBArnoldi<Matrix,Scalar,Primary>(A,R1,X0,B,
base,L,normB,
epsilon,epsilon_R,
arnoldi_ite,Sol,log);
//Check if Convergence
if (nb == arnoldi_ite){//No convergence
k+= restart;
nb_restart++;
//Compute new starting block from solution
X0.CopyBlock(Sol);
//Compute new residual
//Fill matrix R0 with each vector r0_i = B[i] - A * X0[i]
A.MatBlockVect(X0,R0); //R0 = A*X0
for(int i=0 ; i<SizeBlock ; ++i){ //R0 = B - R0
for(int k=0 ; k<dim ; ++k){
R0.getPtr(i)[k] = B.getPtr(i)[k] - R0.getPtr(i)[k];
}
}
//Decomposition QR
R0.ComputeQRFacto(V1,R1);
//Reset datas structure
base.reset();
L.reset();
base.addBlockDatas(V1);
}else{ //Convergence !
k+=nb;
break;
}
}
return k;
}
#endif //IBGMRES_HPP
This diff is collapsed.
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