Commit dd4f1493 authored by Cyrille Piacibello's avatar Cyrille Piacibello

Qr incremental added, IB fully functionnal, API updated, changes is naming,...

Qr incremental added, IB fully functionnal, API updated, changes is naming, everything soon to be released
parent 3e81442a
......@@ -10,9 +10,10 @@ typedef void* IbGMResDr_handle;
/**
* @file Header for using the Ib-GMRes-Dr library from C.
*
* @brief The User must provide a MatrixBlockVector product (could be
* called gemm...). Then he can provide a Right preconditionner as
* another function with similar prototype.
* @brief This Library is Matrix free, meaning that the user must
* provide a method to compute a Matrix x Block of Vector
* product. Then, he can provide a Right PreConditionner x block of
* Vector product too.
*
*/
......@@ -67,7 +68,7 @@ typedef enum ibgmresdr_arithmetic{
* @brief Init the library.
* @param ARI : Arithmetic : integer between 0 and 3 corresponding to
* the arithmetic used.
* @param Matrix : generic ptr to the Matrix structur of the user,
* @param Matrix : generic ptr to the Matrix structure of the user,
* will be provided when Callback_MatBlockVect will be called.
* @param dim : Dimension of problem to be solved
* @param userEnv : generic ptr set by the user. Will be provided to
......@@ -76,7 +77,8 @@ typedef enum ibgmresdr_arithmetic{
* @return Generic ptr, must be provided to each call to Library's
* API function
*/
IbGMResDr_handle Init_ibgmresdr(ibgmresdr_Arithmetic ARI, void * Matrix, int dim, void * userEnv);
IbGMResDr_handle Init_ibgmresdr(ibgmresdr_Arithmetic ARI, void * Matrix, int dim,
void * userEnv);
/**
* @brief Free the Library
......@@ -123,7 +125,8 @@ void ibgmresdr_set_parameters(int max_ite, int restart, void * tolerance,
IbGMResDr_handle handle);
/**
* @brief Solve the system AX=B starting with X0
* @brief Solve the system AX=B starting with X0 using Block GMRes
* algorithm.
*
* @param nbRHS : number of right hand side (will be used as
* SizeBlock)
......@@ -133,6 +136,19 @@ void ibgmresdr_set_parameters(int max_ite, int restart, void * tolerance,
int ibgmresdr_solve(int nbRHS,void * RHS,void *X0,
IbGMResDr_handle handle);
/**
* @brief Solve the system AX=B starting with X0 using Block GMRes
* algorithm with detection of Inexact Breakdowns
*
* @param nbRHS : number of right hand side (will be used as
* SizeBlock)
* @param RHS : Right Hand Side (Colomn major mode)
* @param X0 : Starting block of vector
*/
int ibgmresdr_solve_IB(int nbRHS,void * RHS,void *X0,
IbGMResDr_handle handle);
/**
* @brief Retreive the block of solution computed. The solution is
* stored the same way RHS have been given in the first place.
......
......@@ -2,7 +2,8 @@
#define IBGMRESENGINE_HPP
#include "UserMatrix.hpp"
#include "../src/BlockGMRes.hpp"
#include "../src/BlockGMResQR.hpp"
#include "../src/IBGMRes.hpp"
/**
* @brief This class is used to dereference the handle given by the
* user without knowing in which arithmetic we are.
......@@ -25,6 +26,10 @@ public:
return 0;
};
virtual int solve_IB(int, void * , void *){
return 0;
};
virtual void * get_results(){
return nullptr;
};
......@@ -50,7 +55,7 @@ class IbGMResDrEngine : public EngineDispatcher{
//Solution storage
Block<T,S>* sol;
Logger<S>* log;
Logger<S,5>* log;
public:
IbGMResDrEngine(void * userMatrix, int dim, void * inUserEnv) :
matrix(nullptr),
......@@ -100,10 +105,24 @@ public:
B.InitBlock(reinterpret_cast<T*>(RHS));
//Init an empty block to store solution
sol = new Block<T,S>(nbRHS,matrix->size());
log = new Logger<S>();
log = new Logger<S,5>();
//Call to BlockGMRes with parameters
int res = BlockGMRes<UserMatrix<T,S>,T,S>(*matrix,B,X_init,max_iter,*sol,
restart,*log,tolerance);
int res = BlockGMResQR<UserMatrix<T,S>,T,S>(*matrix,B,X_init,max_iter,*sol,
restart,*log,tolerance);
return res;
}
int solve_IB(int nbRHS,void * RHS,void *X0)override{
//Convert raw data to Block
Block<T,S> X_init(nbRHS,matrix->size());
Block<T,S> B(nbRHS,matrix->size());
X_init.InitBlock(reinterpret_cast<T*>(X0));
B.InitBlock(reinterpret_cast<T*>(RHS));
//Init an empty block to store solution
sol = new Block<T,S>(nbRHS,matrix->size());
log = new Logger<S,5>();
int res = IBGMRes<UserMatrix<T,S>,T,S,5>(*matrix,B,X_init,
max_iter,restart,*sol,*log);
return res;
}
......@@ -149,6 +168,11 @@ extern "C" int ibgmresdr_solve(int nbRHS,void * RHS,void *X0,
return (reinterpret_cast<EngineDispatcher*>(handle))->solve(nbRHS,RHS,X0);
}
extern "C" int ibgmresdr_solve_IB(int nbRHS,void * RHS,void *X0,
IbGMResDr_handle handle){
return (reinterpret_cast<EngineDispatcher*>(handle))->solve_IB(nbRHS,RHS,X0);
}
extern "C" void ibgmresdr_dealloc(IbGMResDr_handle handle){
delete (reinterpret_cast<EngineDispatcher*>(handle));
}
......
......@@ -49,8 +49,8 @@ public:
~UserMatrix(){
}
void MatBlockVect(Block<T,S>& input, Block<T,S>& output){
void* toWrite = output.getPtr();
void MatBlockVect(Block<T,S>& input, Block<T,S>& output, int idxToWrite=0){
void* toWrite = output.getPtr(idxToWrite);
userProduct(userPtrMat,input.getSizeBlock(),input.getPtr(),&toWrite,
userEnv);
}
......
......@@ -38,12 +38,11 @@ void MatxBlockVect(void * mat, int nbVect, void * input, void **output, void * e
Env * my_env = (Env*)env;
int nbLine = my_mat->dim1;
int nbRHS = my_env->nbRHS;
int nbCol = my_mat->dim1;
double complex alpha = 1;
double complex beta = 0;
cblas_zgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,
nbLine,nbRHS,nbCol,&alpha,my_mat->datas,nbCol,
nbLine,nbVect,nbCol,&alpha,my_mat->datas,nbCol,
input,nbCol,&beta,*output,nbCol);
}
......@@ -52,11 +51,10 @@ void RightPreCondxBlockVect(void * rpc, int nbVect, void * input, void **output,
RightPreCond * my_rpc = (RightPreCond *) rpc;
Env * my_env = (Env*)env;
int nbRHS = my_env->nbRHS;
int size = my_rpc->dim;
double complex * inputBlock = input;
double complex * outputBlock = *output;
for(int i=0 ; i<nbRHS ; ++i){
for(int i=0 ; i<nbVect ; ++i){
for(int j=0 ; j<size ; ++j){
outputBlock[i*size+j] = inputBlock[i*size+j] * my_rpc->datas[j];
}
......@@ -132,13 +130,14 @@ int main(int ac , char ** av){
//Give to the library the PreConditionner and the method to call
ibgmresdr_set_rightprecond(RightPreCondxBlockVect,my_rpc,handle);
//Tolerance
double tol = 0.0001;
double tol = 0.000001;
//set param : nb_iter, restart, ptr to tolerance
ibgmresdr_set_parameters((mat->dim1)/userEnv->nbRHS, //nb max ite
((mat->dim1)/userEnv->nbRHS)/2, //size of restart
((mat->dim1)/userEnv->nbRHS), //size of restart
&tol, handle);
//Solve env
//Solve
int nb_iter_needed = ibgmresdr_solve(userEnv->nbRHS,RHS->datas,X0,handle);
//Won't be used (remark : we won't free the memory since it will
......@@ -146,16 +145,16 @@ int main(int ac , char ** av){
double complex * res = ibgmresdr_get_results(handle);
//Get the convergence history
double * historic = ibgmresdr_get_logs(handle);
//double * historic = ibgmresdr_get_logs(handle);
//Write down the convergence
FILE * fd3 = fopen("/home/cyrille/Data/matrices_BEM/ConvergenceHisto.txt",
"w");
fprintf(fd3,"Ite\tMin\tMax\tTime\n");
for(int i=0 ; i<nb_iter_needed ; ++i){
fprintf(fd3,"%d\t%e\t%e\t%e\n",i,historic[3*i],historic[3*i+1],historic[3*i+2]);
}
fclose(fd3);
/* FILE * fd3 = fopen("/home/cyrille/Data/matrices_BEM/ConvergenceHisto.txt", */
/* "w"); */
/* fprintf(fd3,"Ite\tMin\tMax\tTime\n"); */
/* for(int i=0 ; i<nb_iter_needed ; ++i){ */
/* fprintf(fd3,"%d\t%e\t%e\t%e\n",i,historic[3*i],historic[3*i+1],historic[3*i+2]); */
/* } */
/* fclose(fd3); */
//Free the handle
ibgmresdr_dealloc(handle);
......
......@@ -13,7 +13,7 @@ private:
int ldBase; //Leading dimension
Scalar * datas; //raw ptr
int totalNbBlock; //Number of block of vectors
int * sumNbVectInBlock; //Index of starting block (block size may change)
std::vector<int> sumNbVectInBlock; //Index of starting block (block size may change)
int currentBlockUsed; //nbblock used
......@@ -21,7 +21,7 @@ private:
//entry in sumNbVectInBlock
void addBlock(int inputSizeBlock){
currentBlockUsed++;
sumNbVectInBlock[currentBlockUsed] = sumNbVectInBlock[currentBlockUsed-1]+ inputSizeBlock;
sumNbVectInBlock.push_back(sumNbVectInBlock[currentBlockUsed-1]+ inputSizeBlock);
}
......@@ -31,14 +31,13 @@ public:
ldBase(sizeVector),
datas(nullptr),
totalNbBlock(restart),
sumNbVectInBlock(nullptr),
sumNbVectInBlock(0),
currentBlockUsed(0)
{
totalSize = ldBase*restart*SizeBlock;
datas = new Scalar[totalSize];
memset(datas,0,totalSize*sizeof(Scalar));
sumNbVectInBlock = new int[totalNbBlock+1];
memset(sumNbVectInBlock,0,sizeof(int)*(totalNbBlock+1));
sumNbVectInBlock.push_back(0);
// Initializing Base
std::cout<<"Init Base with :\n"
<<"\ttotalSize "<<totalSize<<"\n"
......@@ -49,10 +48,21 @@ public:
}
/**
* This method reset all the counters
*
* Semantically equivalent to ~Base(), Base(int,int,int),
*/
void reset(){
currentBlockUsed = 0;
memset(datas,0,totalSize*sizeof(Scalar));
sumNbVectInBlock.clear();
sumNbVectInBlock.push_back(0);
}
~Base(){
delete [] datas;
delete [] sumNbVectInBlock;
}
......@@ -109,8 +119,8 @@ public:
//Loop over different blocks
for(int i=0 ; i< currentBlockUsed ; ++i){
int sizeCurrBlock = sumNbVectInBlock[i+1]-sumNbVectInBlock[i];
std::cout<<"Ortho against block number "<<i<<" of size "
<<sizeCurrBlock<<"\n";
// std::cout<<"Ortho against block number "<<i<<" of size "
// <<sizeCurrBlock<<"\n";
//Compute hessenberg coeff
call_Tgemm<Scalar>(sizeCurrBlock,nbColinWj, ldBase,
getPtr(i),ldBase,
......@@ -132,10 +142,12 @@ public:
*
*/
void ComputeProduct(Block<Scalar,Primary> & input, Block<Scalar,Primary>& output){
//input.displayBlock("Y, before base.ComputeProduct");
int minLine = std::min(sumNbVectInBlock[currentBlockUsed],input.getLeadingDim());
std::cout<<"Computing Xm solution\t: gemm inputs : "<<ldBase
<<" "<<input.getSizeBlock()
<<" "<<minLine<<"\n";
//displayBase("Base before base.ComputeProduct");
call_gemm<Scalar>(ldBase,input.getSizeBlock(),minLine,
getPtr(),ldBase,
input.getPtr(),input.getLeadingDim(),
......
......@@ -55,7 +55,8 @@ public:
}
void displayBlock(std::string name = ""){
std::cout<<name<<"\n";
std::cout<<name<<"\t LeadingDim x SizeBlock : "<<getLeadingDim()
<<" x "<<getSizeBlock()<<"\n";
for(int a=0 ; a<SizeBlock ; ++a){
for(int b=0 ; b<ldb ; ++b){
std::cout<<datas[a*ldb+b]<<"\t";
......@@ -100,8 +101,11 @@ public:
std::cout<<"Problem on QR facto\n\t output="<<geqrf<<"\n";
}
//Copy triangular R part inside input R
for(int i=0; i<sizeQrBlock ; ++i){
memcpy(R.getPtr(i),&datas[(vStart+i)*ldb+hStart],(i+1)*sizeof(Scalar));
if(R.getSizeBlock() != 0){ //R isn't allocated, then we don't
//copy inside
for(int i=0; i<sizeQrBlock ; ++i){
memcpy(R.getPtr(i),&datas[(vStart+i)*ldb+hStart],(i+1)*sizeof(Scalar));
}
}
//Call to orgqr //orgqr does not exist in Lapack for complex,
//but Lapack has a complex counter part called ungqr, which is
......@@ -149,16 +153,16 @@ public:
Primary getNorm(int id){
//Starting point
Scalar * vect = getPtr(id);
auto norm = module<Primary>(vect[0]);
auto norm = module<Scalar,Primary>(vect[0]);
for(int i=1 ; i<getLeadingDim() ; ++i){
norm+=module<Primary>(vect[i]);
norm+=module<Scalar,Primary>(vect[i]);
}
auto norm2 = std::sqrt(norm);
return norm2;
}
std::pair<Primary,Primary> getMinMaxNorm(std::vector<Primary>& normB){
std::pair<Primary,Primary> minmax{100,0};
std::pair<Primary,Primary> minmax{100000,0};
for(int i=0 ; i<SizeBlock ; ++i){
Primary norm = getNorm(i)/normB[i];
if(norm > minmax.second){
......@@ -205,6 +209,7 @@ public:
//What I need:
char jobu = 'S'; //jobu : S --> Only first size_block vectors computed
char jobvt = 'N'; //jobvt : N --> no part of V^{H} is computed
Primary * sigma = new Primary[getSizeBlock()];
Scalar* Vt = nullptr; //Won't be used
int ldVt = 1; //Won't be used
......
......@@ -3,16 +3,20 @@
#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>
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,
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;
......@@ -34,7 +38,6 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
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());
......@@ -43,13 +46,15 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
//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 j=0 ; j<dim ; ++j){
R0.getPtr(i)[j] = B.getPtr(i)[j] - R0.getPtr(i)[j];
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
......@@ -68,26 +73,32 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
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=1 ; j<=max_iter ; ++j){
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;
......@@ -95,51 +106,33 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
//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";
// //Here we compute Base*Y +X0 and display res
// {
// //if no precond
// // X0 = X0 + Base*Yj
// Block<Scalar,Primary> Resultat{SizeBlock,dim};
// base.ComputeProduct(Y,Resultat);
// for(int i=0 ; i<SizeBlock ; ++i){
// for(int j=0 ; j<dim ; ++j){
// Resultat.getPtr(i)[j] += X0.getPtr(i)[j];
// }
// }
// //Resultat.displayBlock("Current Xm");
// //Compute A*Xm-B, and check norm
// Block<Scalar,Primary> AxSol(SizeBlock,dim);
// A.MatBlockVect(Resultat,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(int i=0 ; i<SizeBlock ; ++i){
// std::cout<<"AAA "<<i<<"\t"<<j<<"\t"<<AxSol.getNorm(i)<<"\n";
// }
// }
//Do we really need that ?
std::pair<Primary,Primary> MinMax = Res.getMinMaxNorm(normB);
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]
......@@ -151,15 +144,12 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
std::cout<<"### Inexact Breakdowns Occured ###!\t\t Ite : \t"<<j<<"\n";
//Create Block to store \W
Block<Scalar,Primary> Directions(SizeBlock,SizeBlock);
//TO BE DISCARDED
Block<Scalar,Primary> tempRpart(U.getSizeBlock(),U.getSizeBlock());
Block<Scalar,Primary> tempRpart;
//Qr facto of bottom block of U
U.ComputeQRFacto(Directions,tempRpart,sumOfSizeOfEachBlock.back());
//Directions.displayBlock("Ortho part of U after QR");
//Update
//V_{j+1} = [P_{j-1},~Wj] * Directions_1
WP.ComputeProduct(Directions,U.getSizeBlock(),
base.getLastColPtr(U.getSizeBlock()),base.getLeadingDim());
......@@ -173,14 +163,16 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
}
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 j=0 ; j<dim ; ++j){
Sol.getPtr(i)[j] += X0.getPtr(i)[j];
for(int k=0 ; k<dim ; ++k){
Sol.getPtr(i)[k] += X0.getPtr(i)[k];
}
}
......
......@@ -3,7 +3,7 @@
#include "LapackInterface.hpp"
#include "Hess.hpp"
#include "HessStandard.hpp"
#include "Base.hpp"
#include "Block.hpp"
#include "Logger.hpp"
......@@ -14,11 +14,11 @@
* Ref : Saad P.241
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<class Matrix,typename Scalar,typename Primary=Scalar>
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>& log,
Logger<Primary,5>& log,
Primary epsilon = 0.00001){
if(A.useRightPreCond()){
......@@ -41,10 +41,12 @@ int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
Block<Scalar,Primary> Yj(SizeBlock,dim+SizeBlock);
//Add a loop for restart
int k=0,nbRestart=0, globalStep=0;
log.init_log();
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
......@@ -93,32 +95,39 @@ int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
Wj.ComputeQRFacto(Q,R);
//Copy Block R inside Hess
Hess.addRPartToHess(R);
if(j>7){//last ite
std::cout<<"Iteration \t\t\t\t"<<j<<"§§§§§§§§§§§§§\n";
R.displayBlock("Last R ?");
}
std::cout<<"QR done"<<std::endl;
//Tests over criterion
Hess.ComputeBlockResidual(R1,Yj);
Yj.displayBlock("Y Computed");
//Here we compute Base*Y +X0 and display res
{
//if no precond
// X0 = X0 + Base*Yj
Block<Scalar,Primary> Resultat{SizeBlock,dim};
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];
}
}
Resultat.displayBlock("Current Xm");
for(int i=0 ; i<SizeBlock ; ++i){
std::cout<<"PREFIX "<<i<<"\t"<<j<<"\t"<<Resultat.getNorm(i)<<"\n";
}
}
// {
// //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;
......@@ -191,7 +200,6 @@ int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
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;
......@@ -209,6 +217,8 @@ int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
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];
......@@ -219,4 +229,4 @@ int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
return globalStep;
}
#endif
#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