Commit 98f18541 authored by Cyrille Piacibello's avatar Cyrille Piacibello

Start of dev of Inexact breakdown

parent caa3f1cb
#ifndef BLOCKARNOLDIIB_HPP
#define BLOCKARNOLDIIB_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 BlockArnoldiIb(Matrix& A,Block<Scalar>& B, Block<Scalar>& X0,
int max_iter, Block<Scalar>& Sol,
Primary epsilon = 0.00001){
//store norms
std::vector<Primary> normB;
for(int i=0 ; i< SizeBlock ; ++i){
normB.push_back(B.template getNorm<Primary>(i));
}
//Compute epsilon_R (step 1)
Primary epsilon_R = normB.front();
for(auto &a : nomrB){
if(a<epsilon_R){
epsilon_R = a;
}
}
epsilon_R *= epsilon;
//Compute first block residual (step 1)
Block<Scalar> 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];
}
}
//V1 will be first block of base, and R1, will be used for least
//square pb, in criteria verification
Block<Scalar> V1(SizeBlock,dim);
Block<Scalar> R1(SizeBlock,SizeBlock);
//Q-R Decomposition of (R0)
R0.ComputeQRFacto(V1,R1);
}
#endif //BLOCKARNOLDIIB_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