Commit 9d30ca9c authored by Cyrille Piacibello's avatar Cyrille Piacibello

Commit because I need to change branch

parent 7877e421
......@@ -149,16 +149,16 @@ public:
return norm2;
}
std::pair<Primary,Primary> getMinMaxNorm(){
std::pair<Primary,Primary> getMinMaxNorm(std::vector<Primary>& normB){
std::pair<Primary,Primary> minmax{100,0};
for(int i=0 ; i<SizeBlock ; ++i){
Primary norm = getNorm(i);
if(norm<minmax.first){
Primary norm = getNorm(i)/normB[i];
std::cout<<"NormRes["<<i<<"]\t"<<norm<<"\n";
if(norm > minmax.second){
minmax.second = norm;
}
if(norm < minmax.first){
minmax.first = norm;
}else{
if(norm>minmax.second){
minmax.second = norm;
}
}
}
return minmax;
......
......@@ -36,7 +36,7 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
//W and P will be stored in the same block
BlockWP<Scalar,Primary> WP(SizeBlock,B.getLeadingDim());
//Store Solution at each iteration
Block<Scalar,Primary> Y(SizeBlock,dim);
Block<Scalar,Primary> Y(SizeBlock,dim+SizeBlock);
//If no Restart, Restart corresponds to the Dimension/Size of Block
HessExtended<Scalar> L(SizeBlock,B.getLeadingDim()/SizeBlock);
......@@ -57,7 +57,7 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
R0.ComputeQRFacto(V1,R1);
//Fisrt block of base added
base.addBlockDatas(V1);
R1.displayBlock("First residual (Gamma)");
//We store the size of each block -- Could be delegated to a class ?
std::vector<int> sizeOfEachBlock;
sizeOfEachBlock.reserve(max_iter);
......@@ -71,7 +71,7 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
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");
//L.displayHessExtended("Hess Before ortho");
Scalar * toWrite = L.getNewStartingCol(WP.getSizeW()); //Where to write
base.ComputeMGSOrtho(WP,toWrite,L.getLeadingDim());
......@@ -90,7 +90,10 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
L.ComputeResidual(R1,Y,Res);
std::cout<<"Residual Computed !"<<std::endl;
//Do we really need that ?
std::pair<Primary,Primary> MinMax = Res.getMinMaxNorm();
std::pair<Primary,Primary> MinMax = Res.getMinMaxNorm(normB);
std::cout<<"Ite ["<<j<<"] Done\t"<<
"Min\t"<<MinMax.first<<"\t"<<
"Max\t"<<MinMax.second<<"\n";
if(MinMax.second < epsilon){
std::cout<<"Convergence ?"<<std::endl;
......@@ -106,7 +109,7 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
Block<Scalar,Primary> Directions(SizeBlock,SizeBlock);
//TO BE DISCARDED
Block<Scalar,Primary> tempRpart(U.getSizeBlock(),U.getSizeBlock());
U.displayBlock("U: Kept directions");
//U.displayBlock("U: Kept directions");
if(U.getSizeBlock() == SizeBlock){//Thus no inexact breakdown occured
......@@ -124,16 +127,13 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
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.displayHessExtended("hess before update bottom line");
//L_{j+1,} = | Directions_1^{H} | * H_j
//Gj = | Directions_2^{H} |
L.UpdateBottomLine(Directions,U.getSizeBlock());
L.displayHessExtended("hess after update bottom line");
//L.displayHessExtended("hess after update bottom line");
sizeOfEachBlock.push_back(U.getSizeBlock());
sumOfSizeOfEachBlock.push_back(U.getSizeBlock()+sumOfSizeOfEachBlock.back());
std::cout<<"Ite ["<<j<<"] Done\t"<<
"Min\t"<<MinMax.first<<"\t"<<
"Max\t"<<MinMax.second<<"\n";
}
//sol = Base*Y + X0
base.ComputeProduct(Y,Sol);
......
......@@ -38,7 +38,7 @@ int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
}
//Residual vectors needed for computing the solution
Block<Scalar,Primary> Yj(SizeBlock,dim);
Block<Scalar,Primary> Yj(SizeBlock,dim+SizeBlock);
//Add a loop for restart
int k=0,nbRestart=0, globalStep=0;
log.init_log();
......@@ -167,7 +167,7 @@ int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
base.addBlockDatas(Q);
std::cout<<"Base has been augmented\n";
globalStep++;
Hess.displayHess("Hess at the end of the loop");
//Hess.displayHess("Hess at the end of the loop");
}
k += nb_iter_restart;
nbRestart += 1;
......
......@@ -117,9 +117,6 @@ public:
r[ldr*i+j] = R1.getPtr(i)[j];
}
}
for(int i=0 ; i<ldr*SizeBlock ; ++i){
std::cout<<r[i]<<"\n";
}
Scalar*workingCopy = new Scalar[SizeBlock*SizeBlock*
(currentSize+1)*currentSize];
......@@ -131,7 +128,6 @@ public:
}
}
displayHess("Hess before gels");
callLapack_gels<Scalar>(ldr,currentSize*SizeBlock,SizeBlock,
workingCopy,ldr,
r,ldr);
......@@ -158,11 +154,11 @@ public:
Block<Scalar,Primary> blockRes(SizeBlock,(currentSize+1)*SizeBlock);
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(),
blockRes.getPtr(),blockRes.getLeadingDim(),
Scalar(1),Scalar(0));
call_gemm<Scalar>((currentSize+1)*SizeBlock,SizeBlock,SizeBlock*currentSize,
datas, ldh,
Y.getPtr(), Y.getLeadingDim(),
blockRes.getPtr(),blockRes.getLeadingDim(),
Scalar(1),Scalar(0));
//blockRes - R
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
......@@ -176,6 +172,7 @@ public:
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
auto norm = blockRes.getNorm(idxRHS)/normB[idxRHS];
//std::cout<<"Norm of vector["<<idxRHS<<"] is "<<norm<<"\n";
std::cout<<"NormRes["<<idxRHS<<"]\t"<<norm<<"\n";
if(norm > MinMax.second){
MinMax.second = norm;
}
......
......@@ -133,7 +133,7 @@ public:
std::cout<<"Gamma : is "<<gamma.getSizeBlock()<<" x "
<<gamma.getLeadingDim()<<"\n";
std::cout<<"Nb Line used : "<<nbLineUsed<<"\n";
displayHessExtended("Hess before gels");
int ldR = nbLineUsed + p;
Scalar * r = new Scalar[ldR*p];
memset(r,0,((nbLineUsed + p)*p)*sizeof(Scalar));
......@@ -142,9 +142,6 @@ public:
memcpy(&r[i*ldR],gamma.getPtr(i),
sizeof(Scalar)*(i+1));
}
for(int i=0 ; i<ldR*p ; ++i){
std::cout<<r[i]<<"\n";
}
//get a working copy of F
Scalar*workingCopy = new Scalar[ldR*sumNbVectInBlck.back()];
memset(workingCopy,0,ldR*sumNbVectInBlck.back()*sizeof(Scalar));
......@@ -166,8 +163,10 @@ public:
}
outputY.displayBlock("Y resulting from least square");
//Call gemm : outputRes = F*Y-gama
std::cout<<"About to compute Block Res:\t"<<nbLineUsed + p
<<"-"<<p<<"-"<<sumNbVectInBlck.back()<<"\n";
call_gemm<>(nbLineUsed + p,p,nbVectInBlck.back(),
call_gemm<>(nbLineUsed + p,p,sumNbVectInBlck.back(),
data,ldHess,
outputY.getPtr(),outputY.getLeadingDim(),
outputResidual.getPtr(),outputResidual.getLeadingDim(),
......
......@@ -149,14 +149,14 @@ int main(int ac, char ** av){
//Compute A starting Block
Block<Scalar,Primary> X0(SizeBlock,dim);
//Fill it!
for(int i=0 ; i<SizeBlock-1 ; ++i){ //last col not filled
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ;j<dim; ++j){
X0.getPtr(i)[j] = {dis(gen),dis(gen)};
}
}
//fill last col with one vect from solution
int last_indice = SizeBlock-1;
memcpy(X0.getPtr(last_indice),X_Exact.getPtr(0),
memcpy(X0.getPtr(last_indice),X_Exact.getPtr(last_indice),
sizeof(Scalar)*X0.getLeadingDim());
Block<Scalar,Primary> sol(SizeBlock,dim);
......@@ -218,7 +218,7 @@ int main(int ac, char ** av){
int nb2 = BlockGMRes<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,X0,5,
sol_ref,5,logger);
std::cout<<"Return value : "<<nb<<"\n";
std::cout<<"Return value : "<<nb2<<"\n";
//Check solution
{
//Compute what i need for my metrics
......
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