Commit 62a4cdbf authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

New version for each arnoldi procedure is functionnal, to use with testBGMRes

parent 8ef263f3
......@@ -5,6 +5,7 @@
#include "Block.hpp"
#include "HessStandard.hpp"
#include "Utils.hpp"
#include "Logger.hpp"
/**
* @brief : Arnoldi iterations
......@@ -13,7 +14,10 @@
* @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 MaxKSize : Maximum size of Krylov Search space.
* @param epsilon : tolerance for Residual
* @param log : logger
*
* @return Number of iterations done
*
* Solution will be stored inside X0 at the end of computation.
......@@ -22,18 +26,23 @@
*/
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){
static ArnReturn Exec_Procedure(Matrix& A,Block<Scalar,Primary>& X0,
Block<Scalar,Primary>& B,
int max_iter,int MaxKSize,Primary epsilon,
Logger<Primary>& log){
std::cout<<"################# Arnoldi standard ##################\n";
std::cout<<"################# Iterations scheduled "<<max_iter<<"\n";
//Initial Size of block
int SizeBlock = B.getSizeBlock();
//Size of matrix
int dim = A.size();
//Return struct
ArnReturn returned_struct;
//Storage needed for computating
Base<Scalar,Primary> base(SizeBlock,max_iter+1,dim);
HessStandard<Scalar> Hess(max_iter,SizeBlock);
Base<Scalar,Primary> base(SizeBlock,MaxKSize,dim);
HessStandard<Scalar> Hess(MaxKSize,SizeBlock);
Block<Scalar,Primary> Yj(SizeBlock,dim+SizeBlock);
Block<Scalar,Primary> R0(SizeBlock,dim);
......@@ -54,13 +63,20 @@ struct Arnoldi{
//Fisrt block of base added
base.addBlockDatas(V1);
int j;
for(j=0 ; j<max_iter ; ++j){
bool convergence = false;
int j = 0;
while(j<max_iter){
double timestamp = Logger<Primary>::GetTime();
std::cout<<"\t\tStart of Iteration\n";
//Compute block Wj
Block<Scalar,Primary> Wj(SizeBlock,dim);
//Where to write in hess
Scalar * toWrite = Hess.getNewCol(Wj.getSizeBlock());
if(nullptr == toWrite){
break;
}
if(A.useRightPreCond()){//Wj = A*preCond*V[j] if preCond
Block<Scalar,Primary> temp(SizeBlock,dim);
A.preCondBaseProduct(base.getPtr(base.getCurrentBlock()-1),temp);
......@@ -70,10 +86,7 @@ struct Arnoldi{
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
......@@ -87,13 +100,16 @@ struct Arnoldi{
//Tests over criterion
Hess.ComputeBlockResidual(R1,Yj);
std::cout<<"["<<j<<"]\t"<<"Residual computed"<<std::endl;
std::pair<Primary,Primary> MinMax;
{//Display residual to see convergence
std::pair<Primary,Primary> MinMax
= Hess.template displayBlockResidual<Primary>(Yj,R1);
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)){
//Adding last ite
log.add_ite(j,SizeBlock,MinMax.first,MinMax.second,timestamp);
convergence = true;
break;
}else{
std::cout<<"Residual is small enough, but not the real's one \n";
......@@ -103,8 +119,10 @@ struct Arnoldi{
//Base augmentation
base.addBlockDatas(Q);
std::cout<<"Base has been augmented\n";
log.add_ite(j,SizeBlock,MinMax.first,MinMax.second,timestamp);
j++;
}
if(j==max_iter){ //No convergence, we need to write current sol into X0
if(!convergence){ //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);
......@@ -119,7 +137,10 @@ struct Arnoldi{
}
}
}
return j;
std::cout<<"################# Iterations done ... "<<j<<"\n";
returned_struct.nbIteDone = j;
returned_struct.hasConverged = convergence;
return returned_struct;
}
};
#endif //ARNOLDI_HP
......@@ -5,6 +5,7 @@
#include "BlockWP.hpp"
#include "HessExtended.hpp"
#include "Utils.hpp"
#include "Logger.hpp"
/**
* @brief : Arnoldi iterations with inexact breakdowns
......@@ -14,6 +15,8 @@
* @param B : Needed for computing real residual if convergence
* @param max_ite : number of iteration to do (at max)
* @param epsilon : tolerance for Residual
* @param log : logger
*
* @return Number of iterations done
*
* Solution will be stored inside X0 at the end of computation.
......@@ -22,22 +25,34 @@
*/
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){
static ArnReturn Exec_Procedure(Matrix& A,Block<Scalar,Primary>& X0,
Block<Scalar,Primary>& B,
int max_iter,int MaxKSize, Primary epsilon,
Logger<Primary>& log){
std::cout<<"########### Arnoldi with Inexact Breakdown ############\n";
//Initial Size of block
int SizeBlock = B.getSizeBlock();
//Size of matrix
int dim = A.size();
//Return struct
ArnReturn returned_struct;
//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> sizeOfEachBlock;
// sizeOfEachBlock.reserve(max_iter);
// sizeOfEachBlock.push_back(SizeBlock);
std::vector<int> sumOfSizeOfEachBlock;
sumOfSizeOfEachBlock.push_back(SizeBlock);
//Counter
int j;
//Set the hessenberg
HessExtended<Scalar> L(MaxKSize,SizeBlock);
Base<Scalar,Primary> base(SizeBlock,MaxKSize,dim);
//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);
//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]
......@@ -51,22 +66,24 @@ struct Arnoldi_IB{
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);
//First block of base added
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){
bool convergence = false;
//Counter
int j=0;
while(j<max_iter){ //Main Loop
double timestamp = Logger<Primary>::GetTime();
std::cout<<"\t\tStart of Iteration\n";
//Get new starting point, and test if restart is triggered
Scalar * toWrite = L.getNewStartingCol(WP.getSizeW()); //Where to write
if(nullptr == toWrite){ //No more room for this iteration
break;
}
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);
......@@ -76,8 +93,6 @@ struct Arnoldi_IB{
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;
......@@ -101,6 +116,9 @@ struct Arnoldi_IB{
if(MinMax.second < epsilon){
//here, add a test on the real residual
if(CheckCurrentSolution(A,B,X0,base,Y,epsilon)){
//Adding last ite
log.add_ite(j,WP.getSizeW(),MinMax.first,MinMax.second,timestamp);
convergence = true;
break;
}
}
......@@ -129,9 +147,14 @@ struct Arnoldi_IB{
//Update
//V_{j+1} = [P_{j-1},~Wj] * Directions_1
WP.ComputeProduct(Directions,U.getSizeBlock(),
base.getLastColPtr(U.getSizeBlock()),
base.getLeadingDim());
Scalar * baseToWrite = base.getLastColPtr(U.getSizeBlock());
if(baseToWrite){
WP.ComputeProduct(Directions,U.getSizeBlock(),
baseToWrite,
base.getLeadingDim());
}else{
break;
}
//Pj = [P_{j-1},~Wj] * Directions_2
WP.ComputeProduct(Directions,U.getSizeBlock());
//L.displayHessExtended("hess before update bottom line");
......@@ -139,12 +162,14 @@ struct Arnoldi_IB{
//Gj = | Directions_2^{H} |
L.UpdateBottomLine(Directions);
}
sizeOfEachBlock.push_back(U.getSizeBlock());
//sizeOfEachBlock.push_back(U.getSizeBlock());
sumOfSizeOfEachBlock.push_back(U.getSizeBlock()+sumOfSizeOfEachBlock.back());
std::cout<<"################## End of Iteration ##################\n\n";
log.add_ite(j,WP.getSizeW(),MinMax.first,MinMax.second,timestamp);
++j;
}
if(j==max_iter){ //No convergence, we need to write current sol into X0
if(!convergence){ //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);
......@@ -159,7 +184,10 @@ struct Arnoldi_IB{
}
}
}
return j;
std::cout<<"################# Iterations done ... "<<j<<"\n";
returned_struct.nbIteDone = j;
returned_struct.hasConverged = convergence;
return returned_struct;
}
};
......
......@@ -13,7 +13,10 @@
* @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 MaxKSize : Maximum size of Krylov Search space.
* @param epsilon : tolerance for Residual
* @param log : logger
*
* @return Number of iterations done
*
* Solution will be stored inside X0 at the end of computation.
......@@ -22,9 +25,10 @@
*/
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){
static ArnReturn Exec_Procedure(Matrix& A,Block<Scalar,Primary>& X0,
Block<Scalar,Primary>& B,
int max_iter,int MaxKSize, Primary epsilon,
Logger<Primary>& log){
std::cout<<"################# Arnoldi with QR Inc ################\n";
//Initial Size of block
......@@ -32,9 +36,12 @@ struct Arnoldi_QRInc{
//Size of matrix
int dim = A.size();
//Return struct
ArnReturn returned_struct;
//Storage needed for computating
Base<Scalar,Primary> base(SizeBlock,max_iter+1,dim);
HessQR<Scalar,Primary> Hess(max_iter+1,SizeBlock);
Base<Scalar,Primary> base(SizeBlock,MaxKSize,dim);
HessQR<Scalar,Primary> Hess(MaxKSize,SizeBlock);
Block<Scalar,Primary> Yj(SizeBlock,dim+SizeBlock);
Block<Scalar,Primary> R0(SizeBlock,dim);
......@@ -52,19 +59,27 @@ struct Arnoldi_QRInc{
Block<Scalar,Primary> R1(SizeBlock,SizeBlock);
//Q-R Decomposition of (R0)
R0.ComputeQRFacto(V1,R1);
//Fisrt block of base added
//First 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){
bool convergence = false;
int j = 0;
while(j<max_iter){
double timestamp = Logger<Primary>::GetTime();
std::cout<<"\t\tStart of Iteration\n";
//Compute block Wj
Block<Scalar,Primary> Wj(SizeBlock,dim);
//Where to write in hess
Scalar * toWrite = Hess.getNewCol(Wj.getSizeBlock());
if(nullptr == toWrite){ //No more room for this iteration
break;
}
if(A.useRightPreCond()){//Wj = A*preCond*V[j] if preCond
Block<Scalar,Primary> temp(SizeBlock,dim);
A.preCondBaseProduct(base.getPtr(base.getCurrentBlock()-1),temp);
......@@ -74,8 +89,6 @@ struct Arnoldi_QRInc{
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;
......@@ -83,7 +96,6 @@ struct Arnoldi_QRInc{
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;
......@@ -91,12 +103,14 @@ struct Arnoldi_QRInc{
//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
<<"\t\tMin\t"<<MinMax.first<<"\tMax\t"<<MinMax.second
<<std::endl;
if(MinMax.second <= epsilon){
if(CheckCurrentSolution(A,B,X0,base,Yj,epsilon)){
std::cout<<"There\n";
//Adding last ite
log.add_ite(j,SizeBlock,MinMax.first,MinMax.second,timestamp);
convergence = true;
break;
}else{
std::cout<<"Residual is small enough, but not the real's one \n";
......@@ -105,8 +119,12 @@ struct Arnoldi_QRInc{
//Base augmentation
base.addBlockDatas(Q);
std::cout<<"Base has been augmented\n";
log.add_ite(j,SizeBlock,MinMax.first,MinMax.second,timestamp);
++j;
}
if(j==max_iter){ //No convergence, we need to write current sol into X0
if(!convergence){ //No convergence, we need to write current sol into X0
//First, call the triangular solver to get the solution
Hess.ComputeCurrentSol(Yj);
Block<Scalar,Primary> Sol{SizeBlock,dim};
if(A.useRightPreCond()){
Block<Scalar,Primary> temp(SizeBlock,dim);
......@@ -121,7 +139,10 @@ struct Arnoldi_QRInc{
}
}
}
return j;
std::cout<<"################# Iterations done ... "<<j<<"\n";
returned_struct.nbIteDone = j;
returned_struct.hasConverged = convergence;
return returned_struct;
}
};
......
......@@ -3,6 +3,7 @@
#include <algorithm>
#include "Block.hpp"
#include "Utils.hpp"
/**
* @brief : GMRes with Inexact Breakdowns, restart and preconditionner
......@@ -13,7 +14,7 @@
* @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 MaxKSize : maximum size of Krylov space
* @param Sol : Block to store the solution
* @param epsilon : Target accuracy
*
......@@ -21,8 +22,9 @@
*/
template<typename Arn,class Matrix,class Scalar,class Primary>
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){
int max_iter, int MaxKSize,
Block<Scalar,Primary>& Sol,Logger<Primary>& log,
Primary epsilon = 0.0000001){
//Cste Size
int SizeBlock = B.getSizeBlock();
......@@ -34,7 +36,7 @@ int BGMRes(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
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<<"Max Size Krylov space before restart :\t"<<MaxKSize<<"\n";
std::cout<<"#######################################################\n";
if(A.useRightPreCond()){
......@@ -56,28 +58,35 @@ int BGMRes(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
// Global step is total number of iteration (thus is <max_iter)
int k=0,nb_restart=0;
int nbIteDone = 0;
while(k<max_iter){
//Compute nb_iter to give to Arnoldi procedure
int arnoldi_ite = std::min(restart,max_iter-k);
int arnoldi_ite = max_iter - nbIteDone;
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;
ArnReturn ret = Arn::Exec_Procedure(A,X0,B,arnoldi_ite,MaxKSize,epsilon,log);
std::cout<<"\n################## Restart Done "<<nb_restart
<<" #################\n";
nbIteDone+=ret.nbIteDone;
//Three case :
//--1 Arnoldi converged
//--2 Arnoldi reached Max Krylov Size allowed
//--3 Arnoldi reached arnoldi_ite itartions
if(ret.hasConverged){//Convergence
break;
}else{
k += nbIteDone;//This will out us from while loop
nb_restart++;
}
}
Sol.CopyBlock(X0);
Sol.ScaleBlock(normB);
B.ScaleBlock(normB);
return k;
return nbIteDone;
}
......
......@@ -12,10 +12,9 @@ private:
int totalSize; //allocated space
int ldBase; //Leading dimension
Scalar * datas; //raw ptr
int totalNbBlock; //Number of block of vectors
std::vector<int> sumNbVectInBlock; //Index of starting block (block size may change)
int currentBlockUsed; //nbblock used
int MaxKSize;
//Use to modify the cursor currentBlockUsed and fill corresponding
//entry in sumNbVectInBlock
......@@ -26,15 +25,14 @@ private:
public:
Base(int sizeBlock, int restart, int sizeVector) : SizeBlock(sizeBlock),
totalSize(0),
ldBase(sizeVector),
datas(nullptr),
totalNbBlock(restart),
sumNbVectInBlock(0),
currentBlockUsed(0)
{
totalSize = ldBase*restart*SizeBlock;
Base(int sizeBlock, int MaxKSize, int sizeVector) : SizeBlock(sizeBlock),
totalSize(0),
ldBase(sizeVector),
datas(nullptr),
sumNbVectInBlock(0),
currentBlockUsed(0),
MaxKSize(MaxKSize){
totalSize = ldBase*MaxKSize;
datas = new Scalar[totalSize];
memset(datas,0,totalSize*sizeof(Scalar));
sumNbVectInBlock.push_back(0);
......@@ -58,6 +56,7 @@ public:
memset(datas,0,totalSize*sizeof(Scalar));
sumNbVectInBlock.clear();
sumNbVectInBlock.push_back(0);
MaxKSize(0);
}
......@@ -85,9 +84,17 @@ public:
return sumNbVectInBlock[currentBlockUsed];
}
/**
* @brief return the last block col
* @brief return the last block col or nullptr if no room
* available. One does not need to check the value, because no
* room left problem should have triggered when asking for new
* room inside Hess.
*/
Scalar * getLastColPtr(int sizeNewBlock){
//Check if there is room for the new block of vect
if(sumNbVectInBlock.back()+sizeNewBlock > MaxKSize){
std::cout<<"WARNING :: No room left for augmenting the base\n";
return nullptr;
}
Scalar * ptr = getPtr(currentBlockUsed);
addBlock(sizeNewBlock);
return ptr;
......@@ -96,10 +103,12 @@ public:
void addBlockDatas(Block<Scalar,Primary>& inBlock){
//Check
if(inBlock.getLeadingDim() != ldBase){std::cout<<"Problem in addBlockDatas\n";}
Scalar * toWrite = getPtr(currentBlockUsed);
memcpy(toWrite,inBlock.getPtr(),inBlock.getSizeBlock()*ldBase*sizeof(Scalar));
addBlock(inBlock.getSizeBlock());
if(sumNbVectInBlock.back()+ inBlock.getSizeBlock() <= MaxKSize){
Scalar * toWrite = getPtr(currentBlockUsed);
memcpy(toWrite,inBlock.getPtr(),
inBlock.getSizeBlock()*ldBase*sizeof(Scalar));
addBlock(inBlock.getSizeBlock());
}
}
......
......@@ -15,7 +15,7 @@ class HessExtended{
int p; //Size of first block
//Relative to total size
int nbTotalCol;
int MaxKSize;
int ldHess;
//Relative to current state for columns
......@@ -35,19 +35,19 @@ class HessExtended{
Scalar * data;
public:
HessExtended(int inP,int inRestart) : p(inP),
nbTotalCol(p*inRestart),
ldHess(p*(inRestart+1)),
nbBlckColUsed(0),
nbLineUsed(0),data(nullptr){
data = new Scalar[ldHess*nbTotalCol];
memset(data,0,sizeof(Scalar)*ldHess*nbTotalCol);
HessExtended(int inMaxKSize,int inP) : p(inP),
MaxKSize(inMaxKSize),
ldHess(inMaxKSize+inP),
nbBlckColUsed(0),
nbLineUsed(0),data(nullptr){
data = new Scalar[ldHess*MaxKSize];
memset(data,0,sizeof(Scalar)*ldHess*MaxKSize);
sumNbVectInBlck.push_back(0);
std::cout<<"Init Hess Extended with :\n"
<<"\ttotalSize "<<nbTotalCol*ldHess<<"\n"
<<"\ttotalSize "<<MaxKSize*ldHess<<"\n"
<<"\tldHess "<<ldHess<<"\n"
<<"\tSize First block "<<p<<"\n"
<<"\tnbCol "<<nbTotalCol<<"\n";
<<"\tnbCol "<<MaxKSize<<"\n";
}
~HessExtended(){
......@@ -55,7 +55,7 @@ public:
}
void reset(){
memset(data,0,sizeof(Scalar)*ldHess*nbTotalCol);
memset(data,0,sizeof(Scalar)*ldHess*MaxKSize);
sumNbVectInBlck.clear();
nbVectInBlck.clear();
nbBlckColUsed = 0;
......@@ -63,10 +63,6 @@ public:
sumNbVectInBlck.push_back(0);
}
Scalar * getPtr(int idx){
return data;
}
int getLeadingDim() const {
return ldHess;
}
......@@ -111,11 +107,12 @@ public:
* the cursors
*/
Scalar * getNewStartingCol(int inSizeBlock){
std::cout<<"Requesting "<<inSizeBlock<<" new vect col\n";
//Check if there is room
if(sumNbVectInBlck.back() + inSizeBlock > MaxKSize){
return nullptr;
}
int ref = sumNbVectInBlck.back();
std::cout<<"sumNbVectInBlck.back() "<<ref<<"\n";
nbVectInBlck.push_back(inSizeBlock);
std::cout<<"nbVectInBlck.back() "<<nbVectInBlck.back()<<"\n";
sumNbVectInBlck.push_back(ref+inSizeBlock);
nbBlckColUsed++;
nbLineUsed += nbVectInBlck.back();
......
......@@ -23,7 +23,7 @@ template<typename Scalar,typename Primary = Scalar>
class HessQR{
//Constants
int SizeBlock;
int restart;
int MaxKSize;
int totalSize;
......@@ -38,27 +38,29 @@ class HessQR{
int idxCurrBlock;
public:
HessQR(int inRestart, int inSizeBlock) : SizeBlock(inSizeBlock),
restart(inRestart),
totalSize(0),
data(nullptr),ldHess(0),
tau(nullptr),RHS(nullptr),
ldRHS(0),idxCurrBlock(0){
HessQR(int inMaxKSize, int inSizeBlock) : SizeBlock(inSizeBlock),
MaxKSize(inMaxKSize),
totalSize(0),
data(nullptr),ldHess(inMaxKSize+SizeBlock),
tau(nullptr),RHS(nullptr),
ldRHS(0),idxCurrBlock(0){
//First set memory
totalSize = restart * SizeBlock * (restart+1) * SizeBlock;
totalSize = MaxKSize * (MaxKSize+SizeBlock);
//max nb of iteration
int maxNbIte = MaxKSize / SizeBlock + 1;
data = new Scalar[totalSize];
ldHess = (restart+1) * SizeBlock;
tau = new Scalar[restart*SizeBlock];
RHS = new Scalar[restart*SizeBlock*SizeBlock];
ldRHS = restart*SizeBlock;
tau = new Scalar[maxNbIte*SizeBlock];
RHS = new Scalar[maxNbIte*SizeBlock*SizeBlock];
ldRHS = maxNbIte*SizeBlock;
//Set memory to zero
memset(data,0,totalSize*sizeof(Scalar));
memset(tau,0,restart*SizeBlock*sizeof(Scalar));
memset(RHS,0,restart*SizeBlock*SizeBlock*sizeof(Scalar));
memset(tau,0,maxNbIte*SizeBlock*sizeof(Scalar));
memset(RHS,0,maxNbIte*SizeBlock*SizeBlock*sizeof(Scalar));
// Initializing Hessenberg
std::cout<<"Init Hess with :\n"
<<"\ttotalSize "<<totalSize<<"\n"
<<"\trestart "<<restart<<"\n"
<<"\tMaxKSize "<<MaxKSize<<"\n"
<<"\tSizeBlock "<<SizeBlock<<"\n"
<<"\tldHess "<<ldHess<<"\n";
}
......@@ -164,7 +166,11 @@ public:
/**
* Return a ptr to the start of the last column block
*/
Scalar * getNewCol(){
Scalar * getNewCol(int inSizeBlock){
//check if room for it
if(idxCurrBlock*SizeBlock + inSizeBlock > MaxKSize){
return nullptr;
}
int ref = idxCurrBlock;
idxCurrBlock++;
return &data[ref*SizeBlock*ldHess];
......@@ -280,25 +286,30 @@ public:
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);
}
ComputeCurrentSol(Y);
}else{
std::cout<<"Residual is not small enough,\nNo need to compute Y\n";
}
return res;
}
void ComputeCurrentSol(Block<Scalar,Primary>& Y){
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));