Commit 3263cb10 authored by MIJIEUX Thomas's avatar MIJIEUX Thomas

Test factorization (begin)

parent 972e51d4
......@@ -135,22 +135,21 @@ public:
int solve(int nbRHS, void *RHS , void *X0) override
{
//Convert raw data to Block
Block<S,P> X_init(nbRHS,matrix->size());
Block<S,P> B(nbRHS,matrix->size());
ConvertStructures(RHS,X0,B,X_init);
Block<S,P> X_init{nbRHS,matrix->size()};
Block<S,P> B{nbRHS,matrix->size()};
ConvertStructures(RHS, X0, B, X_init);
//Init an empty block to store solution
sol = new Block<S,P>(nbRHS, matrix->size());
log = new Logger<P>();
//Compute
using Arnoldi = Arnoldi<UserMatrix<S,P>,S,P>;
int res = call_solve<Arnoldi>(B,X_init,*sol,*log);
int res = call_solve<Arnoldi>(B, X_init, *sol, *log);
return res;
}
//Solve method : RHS and X0 have the same size
int solve_QR(int nbRHS, void *RHS , void *X0) override
int solve_QR(int nbRHS, void *RHS, void *X0) override
{
//Convert raw data to Block
Block<S,P> X_init(nbRHS,matrix->size());
......@@ -162,8 +161,7 @@ public:
log = new Logger<P>();
//Compute
using Arnoldi = Arnoldi_QRInc<UserMatrix<S,P>,S,P>;
int res = call_solve<Arnoldi>(B,X_init,*sol,*log);
int res = call_solve<Arnoldi_QRInc>(B,X_init,*sol,*log);
return res;
}
......@@ -178,9 +176,7 @@ public:
sol = new Block<S,P>(nbRHS,matrix->size());
log = new Logger<P>();
//Compute
using Arnoldi = Arnoldi_IB<UserMatrix<S,P>,S,P>;
int res = call_solve<Arnoldi>(B,X_init,*sol,*log);
int res = call_solve<Arnoldi_IB>(B,X_init,*sol,*log);
return res;
}
......@@ -193,19 +189,19 @@ public:
B.InitBlock(reinterpret_cast<S*>(RHS));
}
template<class Arn>
int call_solve(Block<S,P>& B,
Block<S,P>& X0,
Block<S,P>& sol,
Logger<P>& log)
template<class ARNOLDI,
class Block,
class Primary = typename Block::primary_type
>
int call_solve(Block &B,
Block &X0,
Block &sol,
Logger<Primary>& log)
{
int res = BGMRes<Arn,UserMatrix<S,P>,S,P>(
*matrix,B,X0,
max_iter,restart,
sol,log,tolerance,ortho,
ArnOrtho::RUHE
);
return res;
return BGMRes<ARNOLDI>(*matrix, B, X0,
max_iter, restart,
sol, log, tolerance,
ortho, ArnOrtho::RUHE );
}
void* get_results() override
......@@ -235,54 +231,54 @@ FABULOUS_BEGIN_C_DECL
//Here, write the implementation of API functions
void fabulous_set_usergemm(Callback_MatBlockVect userProduct,
Fabulous_handle handle)
Fabulous_handle handle)
{
(reinterpret_cast<EngineDispatcher*>(handle))->set_usergemm(userProduct);
}
void fabulous_set_rightprecond(Callback_RightCond rightPreCond,
void *ptrRightPreCond,
Fabulous_handle handle)
void *ptrRightPreCond,
Fabulous_handle handle)
{
(reinterpret_cast<EngineDispatcher*>(handle))->set_rightprecond(rightPreCond,
ptrRightPreCond);
}
void fabulous_set_dotProduct(Callback_DotProduct dotProduct,
Fabulous_handle handle)
Fabulous_handle handle)
{
(reinterpret_cast<EngineDispatcher*>(handle))->set_dotProduct(dotProduct);
}
void fabulous_set_parameters(int max_ite, int restart, void *tolerance,
int nbTol,
Fabulous_handle handle)
int nbTol,
Fabulous_handle handle)
{
(reinterpret_cast<EngineDispatcher*>(handle))->set_parameters(max_ite,restart,
tolerance, nbTol);
}
void fabulous_set_ortho_process(enum fabulous_ortho_process value,
Fabulous_handle handle)
Fabulous_handle handle)
{
(reinterpret_cast<EngineDispatcher*>(handle))->set_ortho_process(value);
}
int fabulous_solve(int nbRHS,void *RHS,void *X0,
Fabulous_handle handle)
Fabulous_handle handle)
{
return (reinterpret_cast<EngineDispatcher*>(handle))->solve(nbRHS,RHS,X0);
}
int fabulous_solve_QR(int nbRHS,void *RHS,void *X0,
Fabulous_handle handle)
Fabulous_handle handle)
{
return (reinterpret_cast<EngineDispatcher*>(handle))->solve_QR(nbRHS,RHS,X0);
}
int fabulous_solve_IB(int nbRHS,void *RHS,void *X0,
Fabulous_handle handle)
Fabulous_handle handle)
{
return (reinterpret_cast<EngineDispatcher*>(handle))->solve_IB(nbRHS,RHS,X0);
}
......
......@@ -40,23 +40,20 @@ public:
}
//Return number of col
int size(){
return dim;
}
int size() const { return dim; }
UserMatrix(void * inUserEnv, void * userPtrMatrix, int inDim) :
UserMatrix(void *inUserEnv, void *userPtrMatrix, int inDim):
userProduct(nullptr),
dim(inDim),
userPtrMat(userPtrMatrix),
userPtrPreCond(nullptr),
userEnv(inUserEnv),
usingRightPreCond(false){
}
~UserMatrix(){
usingRightPreCond(false)
{
}
void MatBlockVect(Block<T,S>& input, Block<T,S>& output, int idxToWrite=0){
void MatBlockVect(Block<T,S>& input, Block<T,S>& output, int idxToWrite=0)
{
void* toWrite = output.getPtr(idxToWrite);
std::cout<<"MatBlockVect : Input is "<<input.getSizeBlock()<<" x "
<<input.getLeadingDim()<<"\n";
......@@ -66,45 +63,50 @@ public:
userEnv);
}
void MatBaseProduct(T * ptrToRead,int nbRHS,void * ptrToWrite){
void MatBaseProduct(T * ptrToRead,int nbRHS,void * ptrToWrite)
{
userProduct(userPtrMat,nbRHS,ptrToRead,&ptrToWrite,
userEnv);
}
void setMatBlockVect(Callback_MatBlockVect inUserProduct){
void setMatBlockVect(Callback_MatBlockVect inUserProduct)
{
userProduct = inUserProduct;
}
//Right Pre Conditionner
//Since it is set by the user, we use the same structure to hold it
void setRightPreCond(Callback_RightCond inRightPreCond, void * usrPtrRightpreCond){
void setRightPreCond(Callback_RightCond inRightPreCond, void *usrPtrRightpreCond)
{
rightPreCond = inRightPreCond;
userPtrPreCond = usrPtrRightpreCond;
usingRightPreCond = true;
}
//set Dot product
void setDotProduct(Callback_DotProduct dotProduct){
void setDotProduct(Callback_DotProduct dotProduct)
{
this->dotProduct = dotProduct;
}
bool useRightPreCond(){
return usingRightPreCond;
}
bool useRightPreCond() const { return usingRightPreCond; }
void preCondBlockVect(Block<T,S>& input, Block<T,S>& output){
void preCondBlockVect(Block<T,S>& input, Block<T,S>& output)
{
void* toWrite = output.getPtr();
rightPreCond(userPtrPreCond,input.getSizeBlock(),input.getPtr(),&toWrite,
userEnv);
}
void preCondBaseProduct(T * ptrToRead,Block<T,S>& output){
void preCondBaseProduct(T * ptrToRead,Block<T,S>& output)
{
void* toWrite = output.getPtr();
rightPreCond(userPtrPreCond,output.getSizeBlock(),ptrToRead,&toWrite,
userEnv);
}
void DotProduct(int size, int nbVect, T* vectA, T*vectB, T*res){
void DotProduct(int size, int nbVect, T* vectA, T*vectB, T*res)
{
dotProduct(size,nbVect, vectA,vectB,res,userEnv);
}
};
......
#ifndef FABULOUS_ARITHMETIC_H
#define FABULOUS_ARITHMETIC_H
namespace fabulous {
#define FABULOUS_ARITHMETIC_LIST(ARITHMETIC) \
ARITHMETIC(FLOAT, float, 0) \
ARITHMETIC(DOUBLE, double, 1) \
ARITHMETIC(COMPLEX_FLOAT, std::complex<float>, 2) \
ARITHMETIC(COMPLEX_DOUBLE, std::complex<double>, 3) \
/* **************** Arithmetic CLASS ********************** */
template<class U> class Arithmetik {};
#define FABULOUS_DEFINE_ARITHMETIC(name_, type_, value_) \
template<> \
struct Arithmetik<type_> { \
typedef type_ type; \
static const constexpr char *type_name = #type_; \
static const constexpr char *name = #name_; \
};
FABULOUS_ARITHMETIC_LIST(FABULOUS_DEFINE_ARITHMETIC);
/* **************** Arithmetic Enum ********************** */
#define FABULOUS_ARITHMETIC_TO_ENUM(name_, type_, value_) \
FABULOUS_##name_ = value_,
enum class Arithmetic {
FABULOUS_ARITHMETIC_LIST(FABULOUS_ARITHMETIC_TO_ENUM)
};
}; // end namespace fabulous
#endif // FABULOUS_ARITHMETIC_H
#ifndef IBBGMRESDR_ARNOLDI_HPP
#define IBBGMRESDR_ARNOLDI_HPP
#ifndef FABULOUS_ARNOLDI_HPP
#define FABULOUS_ARNOLDI_HPP
#include "Base.hpp"
#include "Block.hpp"
......@@ -16,7 +16,7 @@
* @param A : User Matrix
* @param X0 : Needed for computing real solution if convergence
* @param B : Needed for computing real residual if convergence
* @param MaxKSize : Maximum size of Krylov Search space.
* @param MaxKrylovSpaceSize : Maximum size of Krylov Search space.
* @param epsilon : tolerance for Residual
* @param log : logger
* @param DR : Struct holding what we need to deflate at restart (see DR_param)
......@@ -28,53 +28,54 @@
* @return Struct containing information about the procedure
*
* 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,
class KI=LapackKernI /* Kernel Interface Concept */
>
struct Arnoldi {
template<class Matrix,
class Block,
class DRParam,
class Scalar = typename Block::value_type,
class Primary = typename Block::primary_type,
class Base = typename DRParam::base_type
>
static ArnReturn Exec_Procedure(
Matrix& A,
Block<Scalar,Primary,KI>& X0,
Block<Scalar,Primary,KI>& B,
int MaxKSize,
std::vector<Primary>& epsilon,
Logger<Primary>& log,
DRParam<Scalar,Primary>& DR,
ArnOrtho ArnChoice,
OrthoChoice ortho = OrthoChoice::MGS)
Matrix &A,
Block &X0,
Block &B,
const int MaxKrylovSpaceSize,
std::vector<Primary> &epsilon,
Logger<Primary> &log,
DRParam &DR,
OrthoChoice ortho = OrthoChoice::MGS,
ArnOrtho ArnChoice = ArnOrtho::BLOCK)
{
std::cout<<"#######################################################\n";
std::cout<<"################# Arnoldi standard ##################\n";
std::cout<<"################# Mat Vect product scheduled "<<MaxKSize<<"\n";
//Initial Size of block
int SizeBlock = B.getSizeBlock();
//Size of matrix
int dim = A.size();
//Needed storage classes
Block<Scalar,Primary,KI> Yj(SizeBlock,MaxKSize/*+SizeBlock*/);
Base<Scalar,Primary,KI>& base = DR.base;
Block<Scalar,Primary,KI> R0;
double timestamp;
//Init hess
HessStandard<Scalar> hess{MaxKSize + 1,SizeBlock};
//V1 will be first block of base, and R1, will be used for least
//square pb, in criteria verification
Block<Scalar,Primary,KI> V1(SizeBlock,dim);
Block<Scalar,Primary,KI> R1;
std::cout<<"########### Mat Vect product scheduled "
<< MaxKrylovSpaceSize << " #############\n";
const int SizeBlock = B.getSizeBlock(); //Initial Size of block
const int dim = A.size(); //Size of matrix
// Needed storage classes
Block Yj{SizeBlock,MaxKrylovSpaceSize/*+SizeBlock*/};
Base &base = DR.base;
Block R0;
HessStandard<Scalar> hess{MaxKrylovSpaceSize + 1,SizeBlock}; //Init hess
// V1 will be first block of base, and R1, will be used for least
// square pb, in criteria verification
Block V1{SizeBlock, dim};
Block R1;
int newK = 0;
double timestamp = 0.0;
if (DR.dr && base.getCurrentBlock() != 0) {
//DR activated AND we already restart once
newK = DeflatedRestart(base,DR.hessSave,R1,DR.k,DR.LS);
newK = DeflatedRestart(base, DR.hessSave, R1, DR.k, DR.LS);
//Init hess with Hess Storage
for (int i=0; i<newK; ++i) {
memcpy(&hess.getPtr()[i*hess.getLeadingDim()],
......@@ -87,67 +88,62 @@ struct Arnoldi {
} else { //No DR, or first run of arnoldi
R0.initData(SizeBlock,dim);
//Fill matrix R0 with each vector r0_i = B[i] - A * X0[i]
A.MatBlockVect(X0,R0); // R0 = A*X0
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];
std::pair<Primary,Primary> minmaxR0 = R0.getMinMaxNorm(A);
timestamp = Logger<Primary>::GetTime();
double timestamp = Logger<Primary>::GetTime();
//Init Log with norm of R0 (only if log is empty)
if (log.getNbIteLogged() == 0) {
log.add_ite(-1,B.getSizeBlock(),
minmaxR0.first,minmaxR0.second,timestamp);
}
//Q-R Decomposition of (R0)
R1.initData(SizeBlock,SizeBlock);
R0.ComputeQRFacto_own(V1,R1,A);
R0.ComputeQRFacto_own(V1,R1,A); // QR Decomposition of (R0)
//First block of base added
// First block of base added
base.addBlockDatas(V1);
}
bool convergence = false;
bool residualCoherence = true;
bool convergence = false, residualCoherence = true;
int ite_meter = 0, j = newK;
/* If DR is used, Hess and Base already contains newK vectors
* (which is equal to k or k+1 in some real cases) */
int ite_meter = 0;
int j = newK; //If DR, Hess and Base already contains newK vectors (which is
//equal to k or k+1 in some real cases)
while (j < MaxKSize) {
while (j < MaxKrylovSpaceSize) {
ite_meter++;
timestamp = Logger<Primary>::GetTime();
std::cout<<"\t\tStart of Iteration\n";
//Compute block Wj
Block<Scalar,Primary,KI> Wj(SizeBlock,dim);
//Where to write in hess
Block Wj{SizeBlock, dim}; // Compute block Wj
// 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,KI> temp(SizeBlock,dim);
Block 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());
Wj.getSizeBlock(), Wj.getPtr());
}
std::cout<<"["<<j<<"]\tAbout to Arnoldi_Ortho(...)\n";
Arnoldi_Ortho(hess,base,Wj,ArnChoice,ortho,A);
//hess.template displayHessBitMap<Primary>("Hess Before Least Square");
//Tests over criterion
hess.ComputeLeastSquareSol(R1,Yj);
Arnoldi_Ortho(hess, base, Wj, ArnChoice, ortho, A);
hess.template ComputeBlockResidual<Primary>(Yj,R1,DR.LS);
// hess.displayHessBitMap("Hess Before Least Square");
// Tests over criterion
hess.ComputeLeastSquareSol(R1, Yj);
hess.ComputeBlockResidual(Yj, R1, DR.LS);
auto evalConv = DR.LS.CheckPrecision(epsilon);
if (evalConv.hasConverged) {//Residual converged !
if (evalConv.hasConverged) { // Residual converged !
if (CheckCurrentSolution(A,B,X0,base,Yj,hess.getNbColUsed(),epsilon)) {
log.add_ite(ite_meter,SizeBlock,evalConv.Min,evalConv.Max,
timestamp);
......@@ -167,13 +163,13 @@ struct Arnoldi {
} // End of main loop
if (!convergence) { //No convergence, we need to write current sol into X0
Block<Scalar,Primary,KI> Sol{SizeBlock,dim};
Block Sol{SizeBlock,dim};
if (A.useRightPreCond()) {
Block<Scalar,Primary,KI> temp(SizeBlock,dim);
base.ComputeProduct(Yj,temp);
A.preCondBlockVect(temp,Sol);
} else {//No pre Cond used
base.ComputeProduct(Yj,Sol);
Block 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)
......@@ -190,4 +186,4 @@ struct Arnoldi {
}
};
#endif // IBBGMRESDR_ARNOLDI_HPP
#endif // FABULOUS_ARNOLDI_HPP
This diff is collapsed.
This diff is collapsed.
#ifndef FABULOUS_ORTHO_IMPL_INCLUDE
#error "Never include this header directly; include 'Arnoldi_Ortho.hpp' instead"
#endif
/**
* @brief Arnoldi Block version
*/
template<class Hess,
class Base,
class Block,
class = typename std::enable_if<!Block::IB_supported::value>::type,
class = void
>
void ArnoldiBlock(Hess &H, Base &base, Block &Wj, OrthoChoice ortho)
{
const int SizeBlock = Wj.getSizeBlock();
const int dim = base.getLeadingDim();
auto *ptrToHess = H.getNewCol(0);
// Ortho against vectors inside base
base.ComputeOrtho(ortho, Wj, ptrToHess, H.getLeadingDim());
std::cout<<"Ortho done"<<std::endl;
Block Q{SizeBlock, dim};
Block R{SizeBlock, SizeBlock};
Wj.ComputeQRFacto(Q, R); // QR facto of orthogonalized Wj
H.addRPartToHess(R); // Copy Block R inside hess
base.addBlockDatas(Q); // Base augmentation
std::cout<<"Block Arnoldi done"<<std::endl;
}
#ifndef FABULOUS_ORTHO_IMPL_INCLUDE
#error "Never include this header directly; include 'Arnoldi_Ortho.hpp' instead"
#endif
/**
* @brief Arnoldi Block version for IB
*/
template<class Hess,
class Base,
class BlockWP,
class = typename std::enable_if<BlockWP::IB_supported::value>::type
>
void ArnoldiBlock(Hess &L, Base &base, BlockWP &WP, OrthoChoice ortho)
{
// We already test if there was room to augment Hessenberg
auto *ptrToHess = L.getNewStartingCol(WP.getSizeW());
//Ortho against vectors inside base
base.ComputeOrtho(ortho, WP, ptrToHess, L.getLeadingDim(), WP.getSizeP());
std::cout << "Ortho done" << std::endl;
// Orthogonalisation of Wj against Pj
WP.OrthoWP(L.getPtrToC(), L.getLeadingDim());
WP.QRofW(L.getPtrToD(),L.getLeadingDim()); // QR of w part
/* Warning : No base augmentation because we need to know if IB
* already happened at this stage */
}
#ifndef FABULOUS_ARNOLDI_ORTHO_ENUM_H
#define FABULOUS_ARNOLDI_ORTHO_ENUM_H
/**
* @brief Enum over the two different approach to ortho phase during Arnoldi.
*
*/
enum class ArnOrtho {
RUHE = 0, /*!< Arnoldi-Ruhe : otrtho vector by vector*/
BLOCK = 1, /*!< Block Arnoldi Ortho block wise followed by a a QR factorization*/
};
/**
* @brief Enum over different available Orthogonalization scheme.
*/
enum class OrthoChoice {
MGS = 0, /*!< Modified Gram Schmidt */
CGS = 1, /*!< Classical Gram Schmidt */
IMGS = 2, /*!< Iterated Modified Gram Schmidt */
ICGS = 3 /*!< Iterated Classical Gram Schmidt */
};
/* **************** STRING CONVERSION ******************** */
static std::unordered_map<ArnOrtho, std::string>
arnortho_names= {
{ ArnOrtho::BLOCK, "BLOCK" },
{ ArnOrtho::RUHE, "RUHE" },
};
static std::unordered_map<OrthoChoice, std::string>
ortho_choice_names = {
{ OrthoChoice::MGS, "MGS" },
{ OrthoChoice::CGS, "CGS" },
{ OrthoChoice::IMGS, "IMGS" },
{ OrthoChoice::ICGS, "ICGS" },
};
/* **************** OPERATOR<< OVERLOAD ********************* */
inline std::ostream& operator<<(std::ostream& os, const ArnOrtho ortho)
{
auto it = arnortho_names.find(ortho);
if (it == arnortho_names.end())
os << "No choice made between BLOCK and RUHE";
else
os << it->second;
return os;
}
inline std::ostream& operator<<(std::ostream& os, const OrthoChoice ortho)
{
auto it = ortho_choice_names.find(ortho);
if (it == ortho_choice_names.end())
os << "Orthogonalization process unknown";
else
os << it->second;
return os;
}
#endif // FABULOUS_ARNOLDI_ORTHO_ENUM_H
#ifndef FABULOUS_ORTHO_IMPL_INCLUDE
#error "Never include this header directly; include 'Arnoldi_Ortho.hpp' instead"
#endif
/* **************************** RUHE ******************************** */
/**
* @brief Arnoldi Ruhe version with CGS ortho
*/
template<class Hess,
class Scalar,
class Primary,
class Matrix,
class KI=LapackKernI /* KernelInterface Concept */
>
void ArnoldiRuhe_CGS(Hess& H,
Base<Scalar,Primary,KI> & base,
Block<Scalar,Primary,KI>& Wj,
Matrix& A)
{
//Cste
int SizeBlock = Wj.getSizeBlock();
int dim = base.getLeadingDim();
//Where to start writing into Hessenberg
Scalar *ptrToHess = H.getNewCol(0);
//Init a new block in base to filled
base.setANewBlock();
//Loop over vector in Wj block
for (int k=0; k<SizeBlock; ++k){
Scalar * ptrToHessK = ptrToHess + (k*H.getLeadingDim());
int nbVectInBase = base.getNbVectUsed();
A.DotProduct(dim,nbVectInBase,base.getPtr(),Wj.getPtr(k),ptrToHessK);
//Wj = Wj + (-1) * Vm * ptrToHessK
KI::gemm(dim,1,nbVectInBase,
base.getPtr(),base.getLeadingDim(),
ptrToHessK,H.getLeadingDim(),
Wj.getPtr(k),Wj.getLeadingDim(),
Scalar(-1),Scalar(1));
// Hess[...] = norm(wj_k)
auto norm = Wj.getNorm(k,A);
ptrToHessK[nbVectInBase] = norm;
//Wj[k] /= norm(Wj[k])
Primary invNorm = Primary{1.0} / norm;
for (int i=0; i<Wj.getLeadingDim(); ++i)
Wj.getPtr(k)[i] *= invNorm;
base.pushOneVect(Wj.getPtr(k));
}
H.incCursor(Wj.getSizeBlock());
std::cout<<"Block Arnoldi done"<<std::endl;
}
/**
* @brief Arnoldi Ruhe version with CGS ortho
*/
template<class Hess,
class Scalar,
class Primary,
class Matrix,
class KI=LapackKernI /* KernelInterface Concept */
>
void ArnoldiRuhe_ICGS(Hess& H,
Base<Scalar,Primary,KI>& base,
Block<Scalar,Primary,KI>& Wj,
Matrix& A)
{
//Cste
int SizeBlock = Wj.getSizeBlock();
int dim = base.getLeadingDim();
//Where to start writing into Hessenberg
Scalar *ptrToHess = H.getNewCol(0);
//Init a new block in base to be filled
base.setANewBlock();
//Loop over vector in Wj block
for (int k=0; k<SizeBlock; ++k){
Scalar * ptrToHessK = ptrToHess + (k*H.getLeadingDim());