Commit 53b5d030 authored by Cyrille Piacibello's avatar Cyrille Piacibello

Work in progress

parent 96b8debc
......@@ -18,6 +18,13 @@ private:
Scalar * datas;
int nbVectUsed;
public:
Block() : SizeBlock(0),
TotalSize(0),
ldb(0),
datas(nullptr),
nbVectUsed(0){
}
Block(int sizeBlock, int sizeVector) : SizeBlock(sizeBlock),
TotalSize(0),
ldb(sizeVector),
......@@ -29,8 +36,22 @@ public:
memset(datas,0,sizeof(Scalar)*TotalSize);
}
/**
* Will be used if default ctor is called, thus if we need to
* create the instance before knowing the size
*/
void initData(int sizeBlock, int sizeVector){
SizeBlock = sizeBlock;
ldb = sizeVector;
TotalSize = ldb * SizeBlock;
datas = new Scalar[TotalSize];
memset(datas,0,sizeof(Scalar)*TotalSize);
}
~Block(){
delete [] datas;
if(datas){
delete [] datas;
}
}
void displayBlock(std::string name = ""){
......@@ -62,30 +83,38 @@ public:
//Lapack part
//QR facto
void ComputeQRFacto(Block& Q, Block& R){ //input are already allocated
void ComputeQRFacto(Block& Q, Block& R, //input are already allocated
int vStart=0, int hStart=0){ //if QR of part of block
//Sizes needed
int sizeQrBlock = SizeBlock - hStart;
//Create Tau
Scalar * tau = new Scalar[SizeBlock];
memset(tau,0,sizeof(Scalar)*SizeBlock);
Scalar * tau = new Scalar[sizeQrBlock];
memset(tau,0,sizeof(Scalar)*sizeQrBlock);
//Qr facto can be done in place
int geqrf = callLapack_geqrf(ldb,SizeBlock,datas,ldb,tau);
int geqrf = callLapack_geqrf(ldb-hStart,sizeQrBlock,
&datas[vStart+hStart*ldb],ldb,tau);
if(geqrf != 0){
std::cout<<"Problem on QR facto\n";
}
//Copy triangular R part inside input R
for(int i=0; i<SizeBlock ; ++i){
memcpy(R.getPtr(i),&datas[i*ldb],(i+1)*sizeof(Scalar));
for(int i=0; i<sizeQrBlock ; ++i){
memcpy(R.getPtr(i),&datas[(hStart+i)*ldb+vStart],(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
//called if Scalar is a complex type
int orgqr = callLapack_orgqr(ldb,SizeBlock,datas,ldb,tau);
int orgqr = callLapack_orgqr(ldb-hStart,sizeQrBlock,
&datas[vStart+hStart*ldb],ldb,tau);
if(orgqr != 0){
std::cout<<"Problem on Q generation after QR facto\n";
}
//Copy output to Q
memcpy(Q.getPtr(),datas,TotalSize*sizeof(Scalar));
for(int i=0 ; i<sizeQrBlock ; ++i){
memcpy(Q.getPtr(i),&datas[vStart+(hStart+i)*ldb],
(ldb-vStart)*sizeof(Scalar));
}
delete [] tau;
}
......@@ -104,6 +133,21 @@ public:
return norm2;
}
std::pair<Primary,Primary> getMinMaxNorm(){
std::pair<Primary,Primary> minmax{100,0};
for(int i=0 ; i<SizeBlock ; ++i){
Primary norm = getNorm(i);
if(norm<minmax.first){
minmax.first = norm;
}else{
if(norm>minmax.second){
minmax.second = norm;
}
}
}
return minmax;
}
/**
* @brief Copy input Block datas inside
*/
......@@ -134,18 +178,21 @@ public:
* @brief Compute SVD of current block, and write down U inside
* output, and Singular Values inside sigma.
*/
void DecompositionSVD(Block<Scalar,Primary>& output, Block<Primary>& sigma){
void DecompositionSVD(Block<Scalar,Primary>& output, Primary epsilon){
//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
//Needed in case Lapack did not converge
Primary * superb = new Primary[output.getSizeBlock()-1];
Primary * superb = new Primary[getSizeBlock()-1];
//Temporary block, will contains the left singular vectors
Block<Scalar,Primary> U(getSizeBlock(),getLeadingDim());
int res = callLapack_gesvd<>(jobu,jobvt,getLeadingDim(),getSizeBlock(),
getPtr(),getLeadingDim(),
sigma.getPtr(),
output.getPtr(),output.getLeadingDim(),
sigma,
U.getPtr(),U.getLeadingDim(),
Vt,ldVt, superb);
if(res!=0){
std::cout<<"Problem ... in SVD, did not converge\n";
......@@ -154,9 +201,17 @@ public:
std::cout<<"\n";
}
}
//Count the number of Singular Value superior to the threshold
//epsilon
int id = 0;
while(sigma[id]>epsilon){
id++;
}
output.initData(id,getLeadingDim());
memcpy(output.getPtr(),U.getPtr(),sizeof(Scalar)*id*U.getLeadingDim());
delete [] sigma;
delete [] superb;
}
};
......
......@@ -11,10 +11,11 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
int max_iter, Block<Scalar,Primary>& Sol,
Primary epsilon = 0.00001){
int SizeBlock = B.getSizeBlock();
//store norms
std::vector<Primary> normB;
for(int i=0 ; i< SizeBlock ; ++i){
for(int i=0 ; i< p ; ++i){
normB.push_back(B.getNorm(i));
}
//Compute epsilon_R (step 1)
......@@ -28,12 +29,12 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
//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,Primary> Pj(SizeBlock,dim);
Base<Scalar> base(p,max_iter+1,dim);
//W and P will be stored in the same block
BlockWP<Scalar,Primary> WP(SizeBlock,B.getLeadingDim());
Hessenberg<Scalar> L(nb_iter_restart,SizeBlock);
//If no Restart, Restart corresponds to the Dimension/Size of Block
HessExtended<Scalar> L(SizeBlock,B.getLeadingDim()/SizeBlock);
//Compute first block residual (step 1)
Block<Scalar,Primary> R0(SizeBlock,dim);
......@@ -53,27 +54,57 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
//Fisrt block of base added
base.addBlockDatas(V1);
//We store the size of each block
//We store the size of each block -- Could be delegated to a class ?
std::vector<int> sizeOfEachBlock;
sizeOfEachBlock.reserve(max_iter);
sizeOfEachBlock.push_back(B.getSizeBlock());
sizeOfEachBlock.push_back(SizeBlock);
std::vector<int> sumOfSizeOfEachBlock;
sumOfSizeOfEachBlock.push_back(SizeBlock);
//Main Loop
for(int j=1 ; j<max_iter ; ++j){
//Compute block Wj
Block<Scalar,Primary> Wj(SizeBlock,dim);
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),Wj);
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),WP);
//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;
//Orthogonalisation de Wj par rapport à Pj
//Stockage des coeffs de l'ortho dans la partie en haut à droite de H|_j
//QR de Wj --> Wj~ et Dj
//Reconstitution de la matrice F ? ou bien tout sera fait in place
//
//Orthogonalisation of Wj against Pj
WP.OrthoWP(L.getPtrToC(),L.getLeadingDim());
//QR of Wj --> {Wj~,Dj}
WP.QRofW(L.getPtrToD(),L.getLeadingDim());
//Create Block for Y and Residual
Block<Scalar,Primary> Y(SizeBlock,L.getNbLineUsed() + SizeBlock);
Block<Scalar,Primary> Res(SizeBlock,L.getNbLineUsed() + SizeBlock);
//Compute residual
L.ComputeResidual(R1,Y,Res);
//Do we really need that ?
std::pair<Primary,Primary> MinMax = getMinMaxNorm();
//Block to Store U
Block<Scalar,Primary> U;
Res.DecompositionSVD(U,epsilon_R);
//Create Block to store W
Block<Scalar,Primary> Directions(U.getSizeBlock(),SizeBlock);
Block<Scalar,Primary> tempRpart; //TO BE DISCARDED
U.ComputeQRFacto(Directions,tempRpart,sumOfSizeOfEachBlock.back());
//Update
//V_{j+1} = [P_{j-1},~Wj] * Directions_1
//Pj = [P_{j-1},~Wj] * Directions_2
//L_{j+1,} = Directions_1^{H} * H_j
//Gj = Directions_2^{H} * H_j
sizeOfEachBlock.push_back(U.getSizeBlock());
sumOfSizeOfEachBlock.push_back(U.getSizeBlock()+sumOfSizeOfEachBlock.back());
}
}
......
......@@ -89,8 +89,10 @@ public:
*
*/
void ComputeProduct(Block<Scalar,Primary>& input,
int nbKeptDir, //Correspond to the number of
//kept directions
Scalar toWrite,int ldToWrite){
call_gemm<>(Parent::getLeadingDim(),input.getSizeBlock(),Parent::getSizeBlock(),
call_gemm<>(Parent::getLeadingDim(),nbKeptDir,Parent::getSizeBlock(),
Parent::getPtr(),Parent::getLeadingDim(),
input.getPtr(),input.getLeadingDim(),
toWrite,ldToWrite,
......@@ -101,13 +103,16 @@ public:
* @brief Compute the product between [P,W] and input block and
* write the result in place of P. (obviously used to update P)
*/
void ComputeProduct(Block<Scalar,Primary>& input){
void ComputeProduct(Block<Scalar,Primary>& input,
int nbKeptDir){ //Correspond to the start
//of discarded directions
int sizeInputBlock = input.getSizeBlock() - nbKeptDir;
//Create temporary blck to store result,
Block<Scalar,Primary> temp(input.getSizeBlock(),Parent::getLeadingDim());
Block<Scalar,Primary> temp(sizeInputBlock,Parent::getLeadingDim());
call_gemm<>(Parent::getLeadingDim(),input.getSizeBlock(),Parent::getSizeBlock(),
call_gemm<>(Parent::getLeadingDim(),sizeInputBlock,Parent::getSizeBlock(),
Parent::getPtr(),Parent::getLeadingDim(),
input.getPtr(),input.getLeadingDim(),
input.getPtr(nbKeptDir),input.getLeadingDim(),
temp.getPtr(),temp.getLeadingDim(),
Scalar(1),Scalar(0));
//temp will be copied back inside P part, and cursor will be
......
......@@ -52,6 +52,10 @@ public:
return ldHess;
}
int getNbLineUsed(){
return nbLineUsed;
}
/**
* @brief return the ptr to the start of new blockcol, and augment
* the cursors
......@@ -71,11 +75,11 @@ public:
* @param W1 Block
*/
template<typename Primary>
void UpdateBottomLine(Block<Scalar,Primary>& W1){
void UpdateBottomLine(Block<Scalar,Primary>& W1,int sizeToRead){
int nj = sumNbVectInBlck.back();
//Create temp block
Block<Scalar,Primary> tempRes(nj,W1.getSizeBlock());
call_Tgemm<>(W1.getSizeBlock(),nj,W1.getLeadingDim(),
call_Tgemm<>(sizeToRead,nj,W1.getLeadingDim(),
W1.getPtr(),W1.getLeadingDim(),
getPtrToG(),ldHess,
tempRes.getPtr(),tempRes.getLeadingDim(),
......@@ -90,7 +94,7 @@ public:
sizeof(Scalar)*diff);
}
//Update nbLine
nbLineUsed += W1.getSizeBlock();
nbLineUsed += sizeToRead;
}
/**
......
......@@ -33,7 +33,7 @@ int main(int ac, char ** av){
Block<Scalar> _Wj(SizeBlock,sizeVect);
memcpy(_Wj.getPtr(), Wj.getPtr(), SizeBlock*sizeVect*sizeof(Scalar));
Wj.ComputeQRFacto(V,R);
Wj.ComputeQRFacto(V,R,1,1);
V.displayBlock("V");
R.displayBlock("R");
......
......@@ -26,15 +26,16 @@ int main(int ac, char ** av){
}
}
//Store Sigma
Block<Primary> Sigma(1,SizeBlock);
//Threshold
Primary epsilon = 2;
//Store U
Block<Scalar,Primary> U(SizeBlock,sizeVect);
Block<Scalar,Primary> U;
B.DecompositionSVD(U,Sigma);
B.DecompositionSVD(U,epsilon);
Sigma.displayBlock("Singular values");
std::cout<<"Number of singular values kept over max number : "
<<U.getSizeBlock()<<"/"<<SizeBlock<<std::endl;
return 0;
}
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