Commit 5869ac58 authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

Refactoring almost complete, lack of logger class, BGMres deal with restart...

Refactoring almost complete, lack of logger class, BGMres deal with restart and scaling mechanisme, then call the arnoldi procedure given in parameters. Arnoldi can be classical arnoldi, arnoldi with Incexact breakdown and Arnoldi with least square optimized.
parent 4b47a240
#ifndef ARNOLDI_HPP
#define ARNOLDI_HPP
#include "Base.hpp"
#include "Block.hpp"
#include "HessStandard.hpp"
#include "Utils.hpp"
/**
* @brief : Arnoldi iterations
*
* @param A : User Matrix
* @param X0 : Needed for computing real solution if convergence
* @param B : Needed for computing real residual if convergence
* @param max_ite : number of iteration to do (at max)
* @param epsilon : tolerance for Residual
* @return Number of iterations done
*
* Solution will be stored inside X0 at the end of computation.
*
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<class Matrix,class Scalar,class Primary=Scalar>
struct Arnoldi{
static int Exec_Procedure(Matrix& A,Block<Scalar,Primary>& X0,
Block<Scalar,Primary>& B,
int max_iter,Primary epsilon){
std::cout<<"################# Arnoldi standard ##################\n";
//Initial Size of block
int SizeBlock = B.getSizeBlock();
//Size of matrix
int dim = A.size();
//Storage needed for computating
Base<Scalar,Primary> base(SizeBlock,max_iter+1,dim);
HessStandard<Scalar> Hess(max_iter,SizeBlock);
Block<Scalar,Primary> Yj(SizeBlock,dim+SizeBlock);
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];
}
}
//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);
int j;
for(j=0 ; j<max_iter ; ++j){
std::cout<<"\t\tStart of Iteration\n";
//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);
std::cout<<"["<<j<<"]\t"<<"Residual computed"<<std::endl;
{//Display residual to see convergence
std::pair<Primary,Primary> MinMax
= Hess.template displayBlockResidual<Primary>(Yj,R1);
std::cout<<"["<<j<<"]\t"<<"Min "<<MinMax.first
<<"\tMax "<<MinMax.second<<"\n";
if(MinMax.second <= epsilon){
if(CheckCurrentSolution(A,B,X0,base,Yj,epsilon)){
break;
}else{
std::cout<<"Residual is small enough, but not the real's one \n";
}
}
}
//Base augmentation
base.addBlockDatas(Q);
std::cout<<"Base has been augmented\n";
}
if(j==max_iter){ //No convergence, we need to write current sol into X0
Block<Scalar,Primary> Sol{SizeBlock,dim};
if(A.useRightPreCond()){
Block<Scalar,Primary> temp(SizeBlock,dim);
base.ComputeProduct(Yj,temp);
A.preCondBlockVect(temp,Sol);
}else{//No pre Cond used
base.ComputeProduct(Yj,Sol);
}
for(int i=0 ; i<SizeBlock ; ++i){
for(int k=0 ; k<dim ; ++k){
X0.getPtr(i)[k] += Sol.getPtr(i)[k];
}
}
}
return j;
}
};
#endif //ARNOLDI_HP
#ifndef ARNOLDI_IB_HPP
#define ARNOLDI_IB_HPP
#include "Base.hpp"
#include "BlockWP.hpp"
#include "HessExtended.hpp"
#include "Utils.hpp"
/**
* @brief : Arnoldi iterations with inexact breakdowns
*
* @param A : User Matrix
* @param X0 : Needed for computing real solution if convergence
* @param B : Needed for computing real residual if convergence
* @param max_ite : number of iteration to do (at max)
* @param epsilon : tolerance for Residual
* @return Number of iterations done
*
* Solution will be stored inside X0 at the end of computation.
*
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<class Matrix,class Scalar,class Primary=Scalar>
struct Arnoldi_IB{
static int Exec_Procedure(Matrix& A,Block<Scalar,Primary>& X0,
Block<Scalar,Primary>& B,
int max_iter,Primary epsilon){
std::cout<<"########### Arnoldi with Inexact Breakdown ############\n";
//Initial Size of block
int SizeBlock = B.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);
//Counter
int j;
//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,max_iter+1,dim);
base.addBlockDatas(V1);
//Set the hessenberg
HessExtended<Scalar> L(SizeBlock,max_iter);
//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){
std::cout<<"\t\tStart of Iteration\n";
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());
}
//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
base.ComputeMGSOrtho(WP,toWrite,L.getLeadingDim(),WP.getSizeP());
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
L.ComputeResidual(R1,Y,Res);
std::cout<<"\nResidual Computed ...\n";
//Do we really need that ?
std::pair<Primary,Primary> MinMax = Res.getMinMaxNorm();
std::cout<<"Ite ["<<j<<"] Done\t"<<
"Min\t"<<MinMax.first<<"\t"<<
"Max\t"<<MinMax.second<<"\n";
if(MinMax.second < epsilon){
//here, add a test on the real residual
if(CheckCurrentSolution(A,B,X0,base,Y,epsilon)){
break;
}
}
//Block to Store U
Block<Scalar,Primary> U;
Res.DecompositionSVD(U,epsilon);
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);
}
sizeOfEachBlock.push_back(U.getSizeBlock());
sumOfSizeOfEachBlock.push_back(U.getSizeBlock()+sumOfSizeOfEachBlock.back());
std::cout<<"################## End of Iteration ##################\n\n";
}
if(j==max_iter){ //No convergence, we need to write current sol into X0
Block<Scalar,Primary> Sol{SizeBlock,dim};
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){
X0.getPtr(i)[k] += Sol.getPtr(i)[k];
}
}
}
return j;
}
};
#endif //ARNOLDI_IB_HPP
#ifndef ARNOLDI_QRINC_HPP
#define ARNOLDI_QRINC_HPP
#include "Base.hpp"
#include "Block.hpp"
#include "HessQR.hpp"
#include "Utils.hpp"
/**
* @brief : Arnoldi iterations
*
* @param A : User Matrix
* @param X0 : Needed for computing real solution if convergence
* @param B : Needed for computing real residual if convergence
* @param max_ite : number of iteration to do (at max)
* @param epsilon : tolerance for Residual
* @return Number of iterations done
*
* Solution will be stored inside X0 at the end of computation.
*
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<class Matrix,class Scalar,class Primary=Scalar>
struct Arnoldi_QRInc{
static int Exec_Procedure(Matrix& A,Block<Scalar,Primary>& X0,
Block<Scalar,Primary>& B,
int max_iter,Primary epsilon){
std::cout<<"################# Arnoldi with QR Inc ################\n";
//Initial Size of block
int SizeBlock = B.getSizeBlock();
//Size of matrix
int dim = A.size();
//Storage needed for computating
Base<Scalar,Primary> base(SizeBlock,max_iter+1,dim);
HessQR<Scalar> Hess(max_iter+1,SizeBlock);
Block<Scalar,Primary> Yj(SizeBlock,dim+SizeBlock);
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];
}
}
//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 the Gamma, that will be used as RHS in
//leastsquare calls.
Hess.initRHS(R1);
int j;
for(j=0 ; j<max_iter ; ++j){
std::cout<<"\t\tStart of Iteration\n";
//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,epsilon);
std::cout<<"["<<j<<"]\t"<<"Residual computed\n"
<<"\t\tMin\t"<<MinMax.first<<"\tMax"<<MinMax.second
<<std::endl;
if(MinMax.second <= epsilon){
if(CheckCurrentSolution(A,B,X0,base,Yj,epsilon)){
std::cout<<"There\n";
break;
}else{
std::cout<<"Residual is small enough, but not the real's one \n";
}
}
//Base augmentation
base.addBlockDatas(Q);
std::cout<<"Base has been augmented\n";
}
if(j==max_iter){ //No convergence, we need to write current sol into X0
Block<Scalar,Primary> Sol{SizeBlock,dim};
if(A.useRightPreCond()){
Block<Scalar,Primary> temp(SizeBlock,dim);
base.ComputeProduct(Yj,temp);
A.preCondBlockVect(temp,Sol);
}else{//No pre Cond used
base.ComputeProduct(Yj,Sol);
}
for(int i=0 ; i<SizeBlock ; ++i){
for(int k=0 ; k<dim ; ++k){
X0.getPtr(i)[k] += Sol.getPtr(i)[k];
}
}
}
return j;
}
};
#endif //ARNOLDI_QRINC_HPP
#ifndef BGMRES_HPP
#define BGMRES_HPP
#include <algorithm>
#include "Block.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 epsilon : Target accuracy
*
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<typename Arn,class Matrix,class Scalar,class Primary=Scalar>
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){
//Cste Size
int SizeBlock = B.getSizeBlock();
int dim = A.size();
std::cout<<"#######################################################\n";
std::cout<<"################## 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";
if(A.useRightPreCond()){
std::cout<<"\nRight Pre Conditionner used\n\n";
}else{
std::cout<<"\nNO - Right Pre Conditionner used\n\n";
}
//Norm of B
std::vector<Primary> normB;
B.Normalize(normB);
//transform from normB to 1/normB
std::vector<Primary> invNormB;
std::transform(normB.begin(),normB.end(),std::back_inserter(invNormB),
[&](const Primary &a)-> Primary{
return 1/a;});
X0.ScaleBlock(invNormB);
// Global step 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 Arnoldi's procedure specified in template
int nb = Arn::Exec_Procedure(A,X0,B,arnoldi_ite,epsilon);
//Check if Convergence
if (nb == arnoldi_ite){//No convergence
k+= restart;
nb_restart++;
}else{ //Convergence !
k+=nb;
break;
}
}
Sol.CopyBlock(X0);
Sol.ScaleBlock(normB);
B.ScaleBlock(normB);
return k;
}
#endif //BGMRES_HPP
......@@ -175,6 +175,49 @@ public:
return minmax;
}
std::pair<Primary,Primary> getMinMaxNorm(){
std::pair<Primary,Primary> minmax{100000,0};
for(int i=0 ; i<SizeBlock ; ++i){
Primary norm = getNorm(i);
if(norm > minmax.second){
minmax.second = norm;
}
if(norm < minmax.first){
minmax.first = norm;
}
}
return minmax;
}
/**
* @brief Normalize all the column vectors of block and store the
* norm inside input norms.
*
*/
void Normalize(std::vector<Primary>& norms){
for(int i=0 ; i< SizeBlock ; ++i){
Primary currNorm = getNorm(i);
norms.push_back(currNorm);
for(int j=0 ; j< getLeadingDim() ; ++j){
getPtr(i)[j] /= currNorm;
}
}
}
/**
* @brief Apply a primary coeff different for each vector inside
* Block.
*/
void ScaleBlock(std::vector<Primary>& coeffs){
for(int i=0 ; i< SizeBlock ; ++i){
for(int j=0 ; j< getLeadingDim() ; ++j){
getPtr(i)[j] *= coeffs[i];
}
}
}
/**
* @brief Copy input Block datas inside
*/
......
......@@ -20,7 +20,7 @@
* by the restart parameter.
*/
template<typename Scalar,typename Primary = Scalar>
class Hessenberg_2{
class HessQR{
//Constants
int SizeBlock;
int restart;
......@@ -38,7 +38,7 @@ class Hessenberg_2{
int idxCurrBlock;
public:
Hessenberg_2(int inRestart, int inSizeBlock) : SizeBlock(inSizeBlock),
HessQR(int inRestart, int inSizeBlock) : SizeBlock(inSizeBlock),
restart(inRestart),
totalSize(0),
data(nullptr),ldHess(0),
......@@ -63,7 +63,7 @@ public:
<<"\tldHess "<<ldHess<<"\n";
}
~Hessenberg_2(){
~HessQR(){
if(data){
delete [] data;
data = nullptr;
......@@ -211,7 +211,7 @@ public:
* inside Y
*/
std::pair<Primary,Primary> LeastSquare(Block<Scalar,Primary>& Y,
std::vector<Primary> & normB){
Primary epsilon){
//Five Steps
//-1 Loop over each block in last colomn to apply Householders
......@@ -267,7 +267,7 @@ public:
for(int j=0 ; j<SizeBlock ; ++j){
acc += module<Scalar,Primary>(start[j]);
}
Primary norm = std::sqrt(acc)/normB[i];
Primary norm = std::sqrt(acc);
if(norm>res.second){
res.second = norm;
}
......@@ -276,16 +276,27 @@ public:
}
}
//-5 Solve a triangular system
int nbLines = idxCurrBlock*SizeBlock;
//First, copy inside Y what's in RHS
for(int i=0 ; i<SizeBlock ; ++i){
memcpy(Y.getPtr(i),&RHS[ldRHS*i],nbLines*sizeof(Scalar));
//Test if Residual is small enough
if(res.second < epsilon){
std::cout<<"Residual is small enough\nComputing Y\n";
//-5 Solve a triangular system
int nbLines = idxCurrBlock*SizeBlock;
//First, copy inside Y what's in RHS
for(int i=0 ; i<SizeBlock ; ++i){
memcpy(Y.getPtr(i),&RHS[ldRHS*i],nbLines*sizeof(Scalar));
}
int trtrs = callLapack_trtrs(nbLines,SizeBlock,
getPtr(),ldHess,
Y.getPtr(),Y.getLeadingDim());
//Check return value
if(trtrs !=0){
std::cout<<"WARNING :: Return values of back solve is "
<<trtrs<<"\tExiting\n";
exit(0);
}
}else{
std::cout<<"Residual is not small enough,\nNo need to compute Y\n";
}
int trtrs = callLapack_trtrs(nbLines,SizeBlock,
getPtr(),ldHess,
Y.getPtr(),Y.getLeadingDim());
return res;
}
......
......@@ -4,7 +4,7 @@
#include "LapackInterface.hpp"
#include "Block.hpp"
/**
* This class is supposed to handle Hessenberg in order to provide
* @brief This class is supposed to handle Hessenberg in order to provide
* directly a ptr to Lapack...
*
* Format:
......@@ -17,7 +17,7 @@
* by the restart parameter.
*/
template<typename Scalar>
class Hessenberg{
class HessStandard{
//Depends on restart
int totalSize;
int restart;
......@@ -29,7 +29,7 @@ class Hessenberg{
Scalar * datas;
public:
Hessenberg(int restart, int SizeBlock) : totalSize(0),restart(restart),
HessStandard(int restart, int SizeBlock) : totalSize(0),restart(restart),
SizeBlock(SizeBlock), currentSize(0),
nbVectInBlock(nullptr),
ldh(SizeBlock*(restart+2)),
......@@ -48,7 +48,7 @@ public:
<<"\tldh "<<ldh<<"\n";
}
~Hessenberg(){
~HessStandard(){
delete [] datas;
delete [] nbVectInBlock;
datas = nullptr;
......@@ -126,8 +126,6 @@ public:
sizeof(Scalar)*ldr);
}
}
std::cout<<"ldr:\t"<<ldr<<" currentSize*SizeBlock:\t"<<currentSize*SizeBlock
<<"\n";
callLapack_gels<Scalar>(ldr,currentSize*SizeBlock,SizeBlock,
workingCopy,ldr,
r,ldr);
......@@ -148,12 +146,12 @@ public:
//Compute the Hess*input set of vectors, and return min and max norm
template<typename Primary>
std::pair<Primary,Primary> displayBlockResidual(Block<Scalar,Primary>& Y,Block<Scalar,Primary>& R,
std::vector<Primary> normB){
std::pair<Primary,Primary> displayBlockResidual(Block<Scalar,Primary>& Y,
Block<Scalar,Primary>& R){
//The Idea is to compute a matrix vector product...
Block<Scalar,Primary> blockRes(SizeBlock,(currentSize+1)*SizeBlock);
std::cout<<"About to compute Block Res:\t"<<(currentSize+1)*SizeBlock
<<"-"<<SizeBlock<<"-"<<SizeBlock*currentSize<<"\n";
// std::cout<<"About to compute Block Res:\t"<<(currentSize+1)*SizeBlock
// <<"-"<<SizeBlock<<"-"<<SizeBlock*currentSize<<"\n";
call_gemm<Scalar>((currentSize+1)*SizeBlock,SizeBlock,SizeBlock*currentSize,
datas, ldh,
Y.getPtr(), Y.getLeadingDim(),
......@@ -170,7 +168,7 @@ public:
//Loop to get the min/max of residuals
std::pair<Primary,Primary> MinMax = std::make_pair<>(1000,0);
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
auto norm = blockRes.getNorm(idxRHS)/normB[idxRHS];
auto norm = blockRes.getNorm(idxRHS);
if(norm > MinMax.second){
MinMax.second = norm;
}
......
#include <random>
#include "../src/Block.hpp"
#include "../src/BGMRes.hpp"
#include "../src/Arnoldi_IB.hpp"
#include "../src/Arnoldi.hpp"
#include "../src/Arnoldi_QRInc.hpp"
#include "../src/LapackInterface.hpp"
#include "../src/Logger.hpp"
/**
* Matrix Class : need to implement MatBlockVect(), MatBaseproduct(),
* useRightPreCond(), preCondBlockVect(),preCondBaseProduct() and
* size() member function in order to works with BlockGMRes algorithm
*
* An abstract class could be used to enforce that.
*/
template<class Scalar>
class InputMatrix{
int dim;
Scalar * data;
public: