Commit 46985994 authored by Cyrille Piacibello's avatar Cyrille Piacibello

So far so good

parent a505c7e5
......@@ -15,11 +15,13 @@ private:
int TotalSize;
int ldb;
Scalar * datas;
int nbVectUsed;
public:
Block(int sizeBlock, int sizeVector) : SizeBlock(sizeBlock),
TotalSize(0),
ldb(sizeVector),
datas(nullptr)
datas(nullptr),
nbVectUsed(0)
{
TotalSize = ldb * SizeBlock;
datas = new Scalar[TotalSize];
......@@ -126,6 +128,15 @@ public:
void InitBlock(Scalar * rawDatas){
memcpy(getPtr(),rawDatas, sizeof(Scalar)*ldb*SizeBlock);
}
//Methods for Pj
/**
* @brief get the number of vectors used
*/
int getUsedVects()const{
return nbVectUsed;
}
};
......
......@@ -25,6 +25,16 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar>& B, Block<Scalar>& X0,
}
}
epsilon_R *= epsilon;
//Here will the restart loop
Base<Scalar> base(SizeBlock,max_iter+1,dim);
//Pj is a matrix ortho built with the abandonned directions, thus
//it's cannot be larger than the range(B).
Block<Scalar> Pj(SizeBlock,dim);
Hessenberg<Scalar> L(nb_iter_restart,SizeBlock);
//Compute first block residual (step 1)
Block<Scalar> R0(SizeBlock,dim);
//Fill matrix R0 with each vector r0_i = B[i] - A * X0[i]
......@@ -40,8 +50,27 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar>& B, Block<Scalar>& X0,
Block<Scalar> R1(SizeBlock,SizeBlock);
//Q-R Decomposition of (R0)
R0.ComputeQRFacto(V1,R1);
//Fisrt block of base added
base.addBlockDatas(V1);
//We store the size of each block
std::vector<int> sizeOfEachBlock;
sizeOfEachBlock.reserve(max_iter);
sizeOfEachBlock.push_back(B.getSizeBlock());
//Main Loop
for(int j=1 ; j<max_iter ; ++j){
//Compute block Wj
Block<Scalar> Wj(SizeBlock,dim);
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),Wj);
//Ortho against V and store coeff inside new block column of Hess L
Scalar * toWrite = L.getPtr(j);
base.ComputeMGSOrtho(Wj,toWrite,L.getLeadingDim());
std::cout<<"Ortho against V done"<<std::endl;
//Same against Pj
}
}
#endif //BLOCKARNOLDIIB_HPP
......@@ -304,5 +304,17 @@ int callLapack_orgqr(int m, int p,
return LAPACKE_zungqr(LAPACK_COL_MAJOR, m, p, p, mat, ldMat,tau);
}
/**
* Computation of SVD
*
* Proto : lapack_int LAPACKE_cgesvd( int matrix_layout, char jobu,
* char jobvt, lapack_int m, lapack_int n, lapack_complex_float* a,
* lapack_int lda, float* s, lapack_complex_float* u, lapack_int ldu,
* lapack_complex_float* vt, lapack_int ldvt, float* superb );
*/
template<typename Scalar>
int callLapack_gesvd(){
}
#endif //LapackInterface.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