Commit 9c4e5ef6 authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

Old stuff removed

parent bee2c8ea
#ifndef BLOCKARNOLDIIB_HPP
#define BLOCKARNOLDIIB_HPP
#include "BlockWP.hpp"
#include "HessExtended.hpp"
#include "Logger.hpp"
#include "Base.hpp"
/**
* @brief : Arnoldi iteration with inexact breakdowns
*
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<class Matrix,typename Scalar,typename Primary=Scalar,int N=5>
int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
int max_iter, Block<Scalar,Primary>& Sol,Logger<Primary,N>& log,
Primary epsilon = 0.0000001){
double tic; //main tic
double local_tic;
//For monitoring purpose
int nbDeflDir = 0;
int SizeBlock = B.getSizeBlock();
int dim = A.size();
//store norms
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;
//Here will the restart loop
Base<Scalar,Primary> base(SizeBlock,max_iter+1,dim);
//W and P will be stored in the same block
BlockWP<Scalar,Primary> WP(SizeBlock,B.getLeadingDim());
//Store Solution at each iteration
Block<Scalar,Primary> Y(SizeBlock,dim+SizeBlock);
//If no Restart, Restart corresponds to the Dimension/Size of Block
HessExtended<Scalar> L(SizeBlock,B.getLeadingDim()/SizeBlock);
log.init_log(SizeBlock);
//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];
}
}
//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);
//We store the size of each block -- Could be delegated to a class ?
std::vector<int> sizeOfEachBlock;
sizeOfEachBlock.reserve(max_iter);
sizeOfEachBlock.push_back(SizeBlock);
std::vector<int> sumOfSizeOfEachBlock;
sumOfSizeOfEachBlock.push_back(SizeBlock);
int j;
//Add first ite
std::pair<Primary,Primary> MinMax = R0.getMinMaxNorm(normB);
log.add_ite(-1,MinMax.first,MinMax.second,B.getSizeBlock());
//Main Loop
for(j=0 ; j<max_iter ; ++j){
tic=Logger<Primary,N>::GetTime();
std::cout<<"################## Start of Iteration #################\n";
//Compute block Wj inside WP
local_tic=Logger<Primary,N>::GetTime();
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),WP.getSizeW(),
WP.getW());
log.add_step(1,local_tic);
//Ortho against V and store coeff inside new block column of Hess L
//L.displayHessExtended("Hess Before ortho");
Scalar * toWrite = L.getNewStartingCol(WP.getSizeW()); //Where to write
local_tic=Logger<Primary,N>::GetTime();
base.ComputeMGSOrtho(WP,toWrite,L.getLeadingDim(),WP.getSizeP());
log.add_step(2,local_tic);
std::cout<<"Ortho against V and P"<<std::endl;
//Orthogonalisation of Wj against Pj
WP.OrthoWP(L.getPtrToC(),L.getLeadingDim());
//QR of Wj --> {Wj~,Dj}
WP.QRofW(L.getPtrToD(),L.getLeadingDim());
//L.template displayHessExtendedBitMap<Primary>("L after Orthogonalization");
std::cout<<"QR of W done"<<std::endl;
//Create Block and Residual
Block<Scalar,Primary> Res(SizeBlock,L.getNbLineUsed() + SizeBlock);
//Compute residual
local_tic=Logger<Primary,N>::GetTime();
L.ComputeResidual(R1,Y,Res);
log.add_step(3,local_tic);
std::cout<<"\nResidual Computed ...\n";
//Do we really need that ?
MinMax = Res.getMinMaxNorm(normB);
std::cout<<"Ite ["<<j<<"] Done\t"<<
"Min\t"<<MinMax.first<<"\t"<<
"Max\t"<<MinMax.second<<"\n";
log.add_ite(j,MinMax.first,MinMax.second,WP.getSizeW());
if(MinMax.second < epsilon){
std::cout<<"Convergence !!"<<std::endl;
break;
}
//Block to Store U
Block<Scalar,Primary> U;
local_tic=Logger<Primary,N>::GetTime();
Res.DecompositionSVD(U,epsilon_R);
log.add_step(4,local_tic);
std::cout<<"Decomposition SVD done !\n\t\tSingular Values kept "
<<U.getSizeBlock()<<" !"<<std::endl;
//Here : If no inexact breakdowns occured, updates procedures are easiers
if(U.getSizeBlock() == SizeBlock){
//Augment base : [P_{j-1},~Wj] = [~Wj]
base.addBlockDatas(WP);
}else{
if(U.getSizeBlock() + nbDeflDir != SizeBlock){
nbDeflDir++;
}
std::cout<<"### Inexact Breakdowns Occured ###!\t\t Ite : \t"<<j<<"\n";
//Create Block to store \W
Block<Scalar,Primary> Directions(SizeBlock,SizeBlock);
Block<Scalar,Primary> tempRpart;
//Qr facto of bottom block of U
U.ComputeQRFacto(Directions,tempRpart,sumOfSizeOfEachBlock.back());
//Update
//V_{j+1} = [P_{j-1},~Wj] * Directions_1
WP.ComputeProduct(Directions,U.getSizeBlock(),
base.getLastColPtr(U.getSizeBlock()),base.getLeadingDim());
//Pj = [P_{j-1},~Wj] * Directions_2
WP.ComputeProduct(Directions,U.getSizeBlock());
//L.displayHessExtended("hess before update bottom line");
//L_{j+1,} = | Directions_1^{H} | * H_j
//Gj = | Directions_2^{H} |
L.UpdateBottomLine(Directions);
//L.template displayHessExtendedBitMap<Primary>("L after update");
}
sizeOfEachBlock.push_back(U.getSizeBlock());
sumOfSizeOfEachBlock.push_back(U.getSizeBlock()+sumOfSizeOfEachBlock.back());
std::cout<<"################## End of Iteration ##################\n"
<<"Base contains exactly :: "<<base.getNbVectUsed()<<"\t!!\n";
log.add_step(0,tic);
}
//sol = Base*Y + X0
base.ComputeProduct(Y,Sol);
for(int i=0 ; i<SizeBlock ; ++i){
for(int k=0 ; k<dim ; ++k){
Sol.getPtr(i)[k] += X0.getPtr(i)[k];
}
}
return j;
}
#endif //BLOCKARNOLDIIB_HPP
#ifndef IBARNOLDI_HPP
#define IBARNOLDI_HPP
#include "BlockWP.hpp"
#include "HessExtended.hpp"
#include "Logger.hpp"
/**
* @brief : Arnoldi iteration with inexact breakdowns
*
* @param A : User Matrix
* @param R1 : First residual
* @param X0 : Needed for computing real solution if convergence
* @param B : Needed for computing real residual if convergence
* @param Base : Base filled only with V1
* @param L : Hessenberg, empty
* @param normB : vector containing the norm of each vector inside block B
* @param epsilon : tolerance for Residual
* @param epsilon_R : tolerance for SVD of Residual
* @param max_ite : number of iteration to do (at max)
* @param Sol : Block to store current solution if convergence or at exit
*
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<class Matrix,typename Scalar,typename Primary=Scalar,int N=5>
int IBArnoldi(Matrix& A,Block<Scalar,Primary>& R1, //Needed for each step
Block<Scalar,Primary>& X0, Block<Scalar,Primary>& B,
Base<Scalar,Primary>& base, //Base has only one block
HessExtended<Scalar>& L, //Hessenberg is empty
std::vector<Primary>& normB,//Provided in order to avoid
//multiple calculations
Primary epsilon,Primary epsilon_R,int max_iter, //Control parameters
Block<Scalar,Primary>& Sol,
Logger<Primary,N>& log){ //Output
//Initial Size of block
int SizeBlock = R1.getSizeBlock();
//Size of matrix
int dim = A.size();
//We store the size of each block -- Could be delegated to a class ?
std::vector<int> sizeOfEachBlock;
sizeOfEachBlock.reserve(max_iter);
sizeOfEachBlock.push_back(SizeBlock);
std::vector<int> sumOfSizeOfEachBlock;
sumOfSizeOfEachBlock.push_back(SizeBlock);
int j;
//Add first ite
std::pair<Primary,Primary> MinMax = R1.getMinMaxNorm(normB);
log.add_ite(-1,MinMax.first,MinMax.second,B.getSizeBlock());
//Create tic :
double tic,local_tic;
//Create Block to Store W and P (if IB happens)
BlockWP<Scalar,Primary> WP(SizeBlock,dim);
//Create Block to Store Y
Block<Scalar,Primary> Y(SizeBlock,dim+SizeBlock);
//Nb directions discarded
int nbDeflDir = 0;
//Main Loop
for(j=0 ; j<max_iter ; ++j){
tic=Logger<Primary,N>::GetTime();
std::cout<<"################## Start of Iteration #################\n";
//Compute block Wj inside WP
local_tic=Logger<Primary,N>::GetTime();
if(A.useRightPreCond()){//Wj = A*preCond*V[j] if preCond
Block<Scalar,Primary> temp(SizeBlock-nbDeflDir,dim);
A.preCondBaseProduct(base.getPtr(base.getCurrentBlock()-1),temp);
A.MatBlockVect(temp,WP,WP.getSizeP());
}else{//Wj = A*V[j] without preCond
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),WP.getSizeW(),
WP.getW());
}
log.add_step(1,local_tic);
//Ortho against V and store coeff inside new block column of Hess L
//L.displayHessExtended("Hess Before ortho");
Scalar * toWrite = L.getNewStartingCol(WP.getSizeW()); //Where to write
local_tic=Logger<Primary,N>::GetTime();
base.ComputeMGSOrtho(WP,toWrite,L.getLeadingDim(),WP.getSizeP());
log.add_step(2,local_tic);
std::cout<<"Ortho against V and P"<<std::endl;
//Orthogonalisation of Wj against Pj
WP.OrthoWP(L.getPtrToC(),L.getLeadingDim());
//QR of Wj --> {Wj~,Dj}
WP.QRofW(L.getPtrToD(),L.getLeadingDim());
std::cout<<"QR of W done"<<std::endl;
//Create Block and Residual
Block<Scalar,Primary> Res(SizeBlock,L.getNbLineUsed() + SizeBlock);
//Compute residual
local_tic=Logger<Primary,N>::GetTime();
L.ComputeResidual(R1,Y,Res);
log.add_step(3,local_tic);
std::cout<<"\nResidual Computed ...\n";
//Do we really need that ?
MinMax = Res.getMinMaxNorm(normB);
std::cout<<"Ite ["<<j<<"] Done\t"<<
"Min\t"<<MinMax.first<<"\t"<<
"Max\t"<<MinMax.second<<"\n";
log.add_ite(j,MinMax.first,MinMax.second,WP.getSizeW());
if(MinMax.second < epsilon){
//here, add a test on the real residual (By computing
//(A*sol -B) with
//sol = Base*Y + X0
if(A.useRightPreCond()){
Block<Scalar,Primary> temp(SizeBlock,dim);
base.ComputeProduct(Y,temp);
A.preCondBlockVect(temp,Sol);
}else{//No pre Cond used
base.ComputeProduct(Y,Sol);
}
for(int i=0 ; i<SizeBlock ; ++i){
for(int k=0 ; k<dim ; ++k){
Sol.getPtr(i)[k] += X0.getPtr(i)[k];
}
}
//Then ||A*Sol -B||
Block<Scalar,Primary> AxSol(SizeBlock,dim);
A.MatBlockVect(Sol,AxSol);
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ; j<dim ; ++j){
AxSol.getPtr(i)[j] -= B.getPtr(i)[j];
}
}
std::pair<Primary,Primary> RealMinMax
= AxSol.getMinMaxNorm(normB);
if(RealMinMax.second < epsilon){
std::cout<<"Convergence !!"<<std::endl;
}
break;
}
//Block to Store U
Block<Scalar,Primary> U;
local_tic=Logger<Primary,N>::GetTime();
Res.DecompositionSVD(U,epsilon_R);
log.add_step(4,local_tic);
std::cout<<"Decomposition SVD done !\n\t\tSingular Values kept "
<<U.getSizeBlock()<<" !"<<std::endl;
//Here : If no inexact breakdowns occured, updates procedures are easiers
if(U.getSizeBlock() == SizeBlock){
//Augment base : [P_{j-1},~Wj] = [~Wj]
base.addBlockDatas(WP);
}else{
if(U.getSizeBlock() + nbDeflDir != SizeBlock){
nbDeflDir = SizeBlock-U.getSizeBlock();
}
std::cout<<"### Inexact Breakdowns Occured ###!\t\t Ite : \t"<<j<<"\n";
//Create Block to store \W
Block<Scalar,Primary> Directions(SizeBlock,SizeBlock);
Block<Scalar,Primary> tempRpart;
//Qr facto of bottom block of U
U.ComputeQRFacto(Directions,tempRpart,sumOfSizeOfEachBlock.back());
//Update
//V_{j+1} = [P_{j-1},~Wj] * Directions_1
WP.ComputeProduct(Directions,U.getSizeBlock(),
base.getLastColPtr(U.getSizeBlock()),base.getLeadingDim());
//Pj = [P_{j-1},~Wj] * Directions_2
WP.ComputeProduct(Directions,U.getSizeBlock());
//L.displayHessExtended("hess before update bottom line");
//L_{j+1,} = | Directions_1^{H} | * H_j
//Gj = | Directions_2^{H} |
L.UpdateBottomLine(Directions);
}
//L.template displayHessExtendedBitMap<Primary>("L after end iteration");
sizeOfEachBlock.push_back(U.getSizeBlock());
sumOfSizeOfEachBlock.push_back(U.getSizeBlock()+sumOfSizeOfEachBlock.back());
std::cout<<"################## End of Iteration ##################\n\n";
log.add_step(0,tic);
}
//Test if convergence occured, if so, solution doesn't need to be
//computed
if(j==max_iter){
//sol = Base*Y + X0
if(A.useRightPreCond()){
Block<Scalar,Primary> temp(SizeBlock,dim);
base.ComputeProduct(Y,temp);
A.preCondBlockVect(temp,Sol);
}else{//No pre Cond used
base.ComputeProduct(Y,Sol);
}
for(int i=0 ; i<SizeBlock ; ++i){
for(int k=0 ; k<dim ; ++k){
Sol.getPtr(i)[k] += X0.getPtr(i)[k];
}
}
}
return j;
}
#endif //IBARNOLDI_HPP
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