Commit ae53c28a authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

IB functionnal, with detyection of IB at restart

parent c4fdc197
......@@ -39,6 +39,8 @@ struct Arnoldi_IB{
ArnReturn returned_struct;
returned_struct.sizeKSpace = 0;
// B.displayBlock("RHS");
//Set the hessenberg
HessExtended<Scalar,Primary> L(MaxKSize,SizeBlock);
Base<Scalar,Primary> base(SizeBlock,MaxKSize,dim);
......@@ -58,9 +60,13 @@ struct Arnoldi_IB{
}
Block<Scalar,Primary> q0(SizeBlock,dim);
Block<Scalar,Primary> r0(SizeBlock,SizeBlock);
// R0.displayBlock("R0");
//Q-R Decomposition of (R0)
R0.ComputeQRFacto(q0,r0);
bool IB_happened = false;
//here, insert inexact breakdown on R0
{
//Copy r0 because needed to compute Lambda, but is
......@@ -72,8 +78,13 @@ struct Arnoldi_IB{
Block<Scalar,Primary> U{r0.getSizeBlock(),r0.getLeadingDim()};
//SVD
int nbKDir = r0_cpy.DecompositionSVD(U,epsilon);
// U.displayBlock("U...");
std::cout<<"U.getSizeBlock() : "<<U.getSizeBlock()
<<"\nnbKDir : "<<nbKDir<<"\n";
//Test if needed to launch algorithm again
if(nbKDir == 0){
IB_happened = true;
std::cout<<"\t\tINEXACT BREAKDOWN on R0\n";
std::cout<<"\t\tNb Kept Dir :"<<nbKDir<<"/"<<SizeBlock<<std::endl;
returned_struct.sizeKSpace = 0;
......@@ -88,6 +99,7 @@ struct Arnoldi_IB{
base.addBlockDatas(q0);
L.initLambda(r0);
}else{
IB_happened = true;
std::cout<<"\t\tINEXACT BREAKDOWN on R0\n";
std::cout<<"\t\tNb Kept Dir :"<<nbKDir<<"/"<<SizeBlock<<std::endl;
//Block empty that will contains the first base vectors
......@@ -108,9 +120,27 @@ struct Arnoldi_IB{
Scalar(1),Scalar(0));
WP.incSizeP(SizeBlock-nbKDir);
base.addBlockDatas(V1ToBeAdded);
Block<Scalar,Primary> Temp{nbKDir,WP.getSizeP()};
{ //Check Ortho of P against V1
call_Tgemm<>(SizeBlock-nbKDir,nbKDir,dim,
WP.getP(),WP.getLeadingDim(),
V1ToBeAdded.getPtr(),V1ToBeAdded.getLeadingDim(),
Temp.getPtr(),Temp.getLeadingDim(),
Scalar(1),Scalar(0));
// Temp.displayBlock("P transpose * V1");
}
// WP.displayBlock("WP after P has been initialized");
std::cout<<"Norm of P : "<<WP.getNorm(0)<<"\n";
base.addBlockDatas(V1ToBeAdded);
//std::cout<<"Init Base with "<<V1ToBeAdded.getSizeBlock()<<" vectors\n";
// V1ToBeAdded.displayBlock("V1ToBeAdded");
//Create a Lambda1 = [U1,U2]^{H} * R0;
//r0.displayBlock("r0");
//q0.displayBlock("q0");
//U.displayBlock("U1U2");
Block<Scalar,Primary> Lambda1{SizeBlock,SizeBlock};
call_Tgemm<>(SizeBlock,SizeBlock,SizeBlock,
......@@ -145,7 +175,7 @@ struct Arnoldi_IB{
if(nullptr == toWrite){ //No more room for this iteration
break;
}
//WP.displayBlock("WP before first ite (but after setup R0)");
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);
......@@ -154,15 +184,20 @@ struct Arnoldi_IB{
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),WP.getSizeW(),
WP.getW());
}
//WP.displayBlock("WP after first multiplication");
returned_struct.sizeKSpace += WP.getSizeW();
//Ortho against V and store coeff inside new block column of Hess L
base.ComputeMGSOrtho(WP,toWrite,L.getLeadingDim(),WP.getSizeP());
std::cout<<"Ortho against V and P"<<std::endl;
std::cout<<"Ortho against V done"<<std::endl;
//Orthogonalisation of Wj against Pj
// L.displayHessExtendedBitMap("Before adding C part");
WP.OrthoWP(L.getPtrToC(),L.getLeadingDim());
//WP.CheckOrthoPW();
//QR of Wj --> {Wj~,Dj}
// L.displayHessExtendedBitMap("After adding C part");
WP.QRofW(L.getPtrToD(),L.getLeadingDim());
// L.displayHessExtendedBitMap("After adding D part");
std::cout<<"QR of W done"<<std::endl;
//Create Block and Residual
......@@ -194,9 +229,11 @@ struct Arnoldi_IB{
int nbDirKept = Res.DecompositionSVD(U,epsilon);
std::cout<<"Decomposition SVD done !\n\t\tSingular Values kept "
<<nbDirKept<<" !"<<std::endl;
if(nbDirKept != SizeBlock){
IB_happened = true;
}
//Here : If no inexact breakdowns occured, updates procedures are easier
if(U.getSizeBlock() == SizeBlock){
if(!IB_happened){
//Augment base : [P_{j-1},~Wj] = [~Wj]
base.addBlockDatas(WP);
}else{
......@@ -232,6 +269,7 @@ struct Arnoldi_IB{
//Gj = | Directions_2^{H} |
L.UpdateBottomLine(Directions);
}
//base.checkOrtho();
//sizeOfEachBlock.push_back(U.getSizeBlock());
sumOfSizeOfEachBlock.push_back(U.getSizeBlock()+sumOfSizeOfEachBlock.back());
......
......@@ -119,6 +119,7 @@ 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;
......@@ -141,7 +142,7 @@ struct Arnoldi_QRInc{
}
}
std::cout<<"################# Iterations done ... "<<j<<"\n";
returned_struct.sizeKSpace = (j+1)*SizeBlock;
returned_struct.sizeKSpace = j*SizeBlock;
returned_struct.nbIteDone = j;
returned_struct.hasConverged = convergence;
return returned_struct;
......
......@@ -139,6 +139,7 @@ public:
*/
void ComputeMGSOrtho(Block<Scalar,Primary> & Wj, Scalar * HessToWrite, int ldHess,
int idxToStart = 0){
std::cout<<"Idx To Start :: "<<idxToStart<<"\n";
Scalar * ptrOfHess = HessToWrite;
int nbColinWj = Wj.getSizeBlock() - idxToStart;
//Loop over different blocks
......@@ -146,13 +147,20 @@ public:
int sizeCurrBlock = sumNbVectInBlock[i+1]-sumNbVectInBlock[i];
// std::cout<<"Ortho against block number "<<i<<" of size "
// <<sizeCurrBlock<<"\n";
//Compute hessenberg coeff
call_Tgemm<Scalar>(sizeCurrBlock,nbColinWj, ldBase,
getPtr(i),ldBase,
Wj.getPtr(idxToStart),Wj.getLeadingDim(),
ptrOfHess, ldHess, 1, 0);
ptrOfHess, ldHess,
Scalar(1), Scalar(0));
//Modify Wj
//Wj = (-1) * V_i * Hij + Wj;
// std::cout<<"Params : \n"
// <<"ldBase\t:: "<<ldBase<<"\n"
// <<"nbColinWj\t:: "<<nbColinWj<<"\n"
// <<"sizeCurrBlock\t:: "<<sizeCurrBlock<<"\n"
// <<"Wj.getLeadingDim()\t:: "<<Wj.getLeadingDim()<<"\n";
call_gemm<Scalar>(ldBase,nbColinWj,sizeCurrBlock,
getPtr(i),ldBase,
ptrOfHess, ldHess,
......@@ -168,10 +176,12 @@ public:
*/
void ComputeProduct(Block<Scalar,Primary> & input, Block<Scalar,Primary>& output){
//input.displayBlock("Y, before base.ComputeProduct");
int minLine = std::min(sumNbVectInBlock[currentBlockUsed],input.getLeadingDim());
std::cout<<"Computing Xm solution\t: gemm inputs : "<<ldBase
<<" "<<input.getSizeBlock()
<<" "<<minLine<<"\n";
<<"min between "<<sumNbVectInBlock[currentBlockUsed]
<<" and "<< input.getLeadingDim() <<"\n";
int minLine = std::min(sumNbVectInBlock[currentBlockUsed],input.getLeadingDim());
//displayBase("Base before base.ComputeProduct");
call_gemm<Scalar>(ldBase,input.getSizeBlock(),minLine,
getPtr(),ldBase,
......
......@@ -57,10 +57,10 @@ public:
void displayBlock(std::string name = ""){
std::cout<<name<<"\t LeadingDim x SizeBlock : "<<getLeadingDim()
<<" x "<<getSizeBlock()<<"\n";
for(int a=0 ; a<SizeBlock ; ++a){
for(int a=0 ; a<ldb ; ++a){
std::cout<<"["<<a<<"] ";
for(int b=0 ; b<ldb ; ++b){
std::cout<<datas[a*ldb+b]<<"\t";
for(int b=0 ; b<getSizeBlock() ; ++b){
std::cout<<datas[b*ldb+a]<<"\t";
}
std::cout<<std::endl;
}
......@@ -278,6 +278,13 @@ public:
std::cout<<"\n";
}
}
//DisplaySingular values
// std::cout<<"DisplaySingular values\n";
// for(int i=0 ; i<getSizeBlock() ; ++i){
// std::cout<<sigma[i]<<"\n";
// }
//Count the number of Singular Value superior to the threshold
//epsilon
int id = 0;
......
......@@ -59,6 +59,7 @@ public:
toWrite,ldToWrite,
getW(),Parent::getLeadingDim(),
Scalar(-1),Scalar(1));
std::cout<<"WP :: Ortho done\n";
}
}
......@@ -116,10 +117,15 @@ public:
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(sizeInputBlock,Parent::getLeadingDim());
std::cout<<"Updating P part :\n"
<<"nbKeptDir :: "<<nbKeptDir<<"\n"
<<"sizeInputBlock :: "<<sizeInputBlock<<"\n";
call_gemm<>(Parent::getLeadingDim(),sizeInputBlock,Parent::getSizeBlock(),
Parent::getPtr(),Parent::getLeadingDim(),
input.getPtr(nbKeptDir),input.getLeadingDim(),
......@@ -135,7 +141,20 @@ public:
//Not necessary, should be written over at the next iteration
//but still...
memset(getW(),0,getSizeW()*Parent::getLeadingDim()*sizeof(Scalar));
//Parent::displayBlock("WP at the end of iteration");
}
//Check othogonality between p and W
void CheckOrthoPW(){
Block<Scalar,Primary> temp{getSizeW(),getSizeP()};
call_Tgemm(getSizeP(),getSizeW(),Parent::getLeadingDim(),
getP(),Parent::getLeadingDim(),
getW(),Parent::getLeadingDim(),
temp.getPtr(),temp.getLeadingDim(),
Scalar(1),Scalar(0));
temp.displayBlock("Ortho between P and W");
}
};
......
......@@ -78,6 +78,7 @@ public:
}
Lambda1 = new Block<Scalar,Primary>{p,p};
memcpy(Lambda1->getPtr(),lambda.getPtr(),p*p*sizeof(Scalar));
//Lambda1->displayBlock("lambda1 at init");
}
......@@ -92,6 +93,7 @@ public:
for(int i=0 ; i<Lambda1->getSizeBlock() ; ++i){
Phi->getPtr(i)[i] = Scalar(1);
}
// Phi->displayBlock("Phi at init");
}else{
std::cout<<"Should not be there !!\nExiting\n";
exit(0);
......@@ -246,6 +248,13 @@ public:
Scalar(1),Scalar(0));
std::cout<<"Computing RHS done at ite "
<< sumNbVectInBlck.size() <<"\n";
// std::cout<<"Display RHS : \n";
// for(int i=0 ; i<ldRHS ; i++){
// for(int j=0 ; j<p ; ++j){
// std::cout<<RHS[j*ldRHS + i]<<" ";
// }
// std::cout<<"\n";
// }
}
/**
......@@ -258,8 +267,8 @@ public:
//Complete gamma
int ldR = nbLineUsed + p;
Scalar * r = new Scalar[ldR*p];
memset(r,0,ldR*p*sizeof(Scalar));
if(!Phi){
memset(r,0,((nbLineUsed + p)*p)*sizeof(Scalar));
for(int i=0 ; i<p ; ++i){
memcpy(&r[i*ldR],Lambda1->getPtr(i),
sizeof(Scalar)*(i+1));
......@@ -267,14 +276,29 @@ public:
}else{
ComputeRHS(r,ldR);
}
Block<Scalar,Primary> RHS{p,outputResidual.getLeadingDim()};
for(int i=0 ; i<outputResidual.getSizeBlock();++i){
for(int j=0 ; j<outputResidual.getLeadingDim() ; ++j){
RHS.getPtr(i)[j] = r[i*ldR+j];
}
}
//get a working copy of F
// displayHessExtendedBitMap("Before calling gels");
//displayHessExtended("Same");
Scalar*workingCopy = new Scalar[ldR*sumNbVectInBlck.back()];
memset(workingCopy,0,ldR*sumNbVectInBlck.back()*sizeof(Scalar));
for(int i=0 ; i<sumNbVectInBlck.back() ; ++i){
memcpy(&workingCopy[i*ldR],&data[i*ldHess],
sizeof(Scalar)*ldR);
}
// displayHessExtendedBitMap<Primary>("Residual");
// displayHessExtendedBitMap("Residual");
std::cout<<"Computing gels ::\n"
<<"\t ldr :: "<<ldR<<"\n"
<<"\t sumNbVectInBlck.back() :: "<<sumNbVectInBlck.back()<<"\n"
<<"\t p :: "<<p<<"\n";
//Call least square
int gels = callLapack_gels<>(ldR,sumNbVectInBlck.back(),p,
......@@ -288,6 +312,8 @@ public:
memcpy(outputY.getPtr(i),&r[i*ldR],sizeof(Scalar)*sumNbVectInBlck.back());
}
// outputY.displayBlock("Y after least square");
//Call gemm : outputRes = F*Y-gama
call_gemm<>(nbLineUsed + p,p,sumNbVectInBlck.back(),
data,ldHess,
......@@ -295,15 +321,20 @@ public:
outputResidual.getPtr(),outputResidual.getLeadingDim(),
Scalar(1),Scalar(0));
// outputResidual.displayBlock("F * Y");
std::cout<<"Compute F*Y done\n";
//Add gamma to residual
for(int i=0 ; i<outputResidual.getSizeBlock();++i){
for(int j=0 ; j<i+1 ; ++j){
outputResidual.getPtr(i)[j] -= Lambda1->getPtr(i)[j];
for(int j=0 ; j<outputResidual.getLeadingDim() ; ++j){
outputResidual.getPtr(i)[j] -= RHS.getPtr(i)[j];
}
}
// outputResidual.displayBlock("Residu after least square");
delete [] workingCopy;
delete [] r;
}
......
......@@ -166,7 +166,9 @@ void call_gemm(int m, int n, int k,
double alpha_raw[2] = {alpha.real(),alpha.imag()};
double beta_raw[2] = {beta.real(),beta.imag()};
cblas_zgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,
m,n,k,alpha_raw,base,ldBase,block,ldBlock,beta_raw,
m,n,k,
alpha_raw,base,ldBase,
block,ldBlock,beta_raw,
res,ldRes);
}
......
......@@ -32,12 +32,16 @@ public:
return dim;
}
Scalar * getPtr(){
return data;
}
void initMatrix(std::random_device& rd){
std::mt19937 gen(/*rd()*/0);
std::uniform_real_distribution<> dis(0, 2);
std::uniform_real_distribution<> dis(0, 6);
for(int i=0 ; i<dim ; ++i){
for(int j=0 ; j<dim ; ++j){
data[dim*i+j] = dis(gen);
data[dim*i+j] = std::floor(dis(gen));
}
}
// Scalar diag = dim;
......@@ -52,6 +56,17 @@ public:
// }
}
void displayMat(){
std::cout<<"Matrice A ::\n";
std::cout<<"{";
for(int i=0 ; i<dim ; ++i){
for(int j=0 ; j< dim ; ++j){
std::cout<<data[j*dim+i]<<" ";
}
std::cout<<";\n";
}
}
void MatBlockVect(Block<Scalar>& in, Block<Scalar>& out,
int idxToWrite=0){
call_gemm<Scalar>(dim,in.getSizeBlock(),dim,
......@@ -93,7 +108,7 @@ int main(int ac, char ** av){
for(auto ag : args)
std::cout << ag << "\n";
int SizeBlock=2,nbBlock=20,restart=15;
int SizeBlock=3,nbBlock=10,restart=15;
if (args.size()>1){
SizeBlock = std::stoi(args[1]);
......@@ -119,39 +134,51 @@ int main(int ac, char ** av){
//Create matrix
Matrix mat(dim);
mat.initMatrix(rd);
mat.displayMat();
//Create X_Exact
Block<Scalar> X_Exact(SizeBlock,dim);
//Fill it!
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ;j<dim; ++j){
X_Exact.getPtr(i)[j] = dis(gen);
//Fill the SizeBLock - 1 first vectors with random
for(int i=0 ; i<SizeBlock-1 ; ++i){
X_Exact.getPtr(i)[i] = Scalar(1);
}
Scalar epsilon = 0.0001;
//Fill the last col with the sum of the others
for(int i=0 ; i< dim ; ++i){
for(int j=0 ; j<SizeBlock-1 ; ++j){
X_Exact.getPtr(SizeBlock-1)[i] += X_Exact.getPtr(j)[i];
}
}
//Compute RHS as Mat*X_Exact
Block<Scalar> RHS(SizeBlock,dim);
mat.MatBlockVect(X_Exact,RHS);
X_Exact.displayBlock("X_Exact");
RHS.displayBlock("RHS");
//Compute a starting Block
Block<Scalar> X0(SizeBlock,dim);
//Fill it!
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ;j<dim; ++j){
X0.getPtr(i)[j] = dis(gen);
}
}
//Copy inside first col of X0 the first col of B
memcpy(X0.getPtr(),X_Exact.getPtr(0),sizeof(Scalar)*dim);
// for(int i=0 ; i<SizeBlock ; ++i){
// for(int j=0 ;j<dim; ++j){
// X0.getPtr(i)[j] = dis(gen);
// }
// }
Block<Scalar> sol(SizeBlock,dim);
Logger<Primary> log;
//Using for the arnoldi procedure choosen
using Arnoldi1 = Arnoldi_QRInc<InputMatrix<Scalar>,Scalar,Primary>;
using Arnoldi2 = Arnoldi<InputMatrix<Scalar>,Scalar,Primary>;
using Arnoldi3 = Arnoldi_IB<InputMatrix<Scalar>,Scalar,Primary>;
int nb = BGMRes<Arnoldi3,Matrix,Scalar,Primary>(mat,RHS,X0,nbBlock+4,restart,sol,log,
0.000001);
epsilon);
std::cout<<"Total Nb iter done::\t"<<nb<<"\n";
//Check solution
......
......@@ -99,7 +99,7 @@ public:
}
bool useRightPreCond(){
return true;
return false;
}
void preCondBlockVect(Block<Scalar,Primary>& input, Block<Scalar,Primary>& output){
for(int j=0 ; j<output.getSizeBlock() ; ++j){
......@@ -171,7 +171,7 @@ int main(int ac, char ** av){
Block<Scalar,Primary> RHS{SizeBlock,dim};
Block<Scalar,Primary> XExact{SizeBlock,dim};
for(int i=0 ; i<SizeBlock ; ++i){
XExact.getPtr(i)[i] = Primary(1);
XExact.getPtr(i)[i] = Scalar(1);
}
mat.MatBlockVect(XExact,RHS);
......@@ -257,35 +257,35 @@ int main(int ac, char ** av){
});
minmax2 = CompareBlocks<Scalar,Primary>(sol,XExact);
}
{
std::cout<<"\tVersion with Optimized Least Square\n\n";
Block<Scalar,Primary> X0(SizeBlock,dim);
Block<Scalar,Primary> sol(SizeBlock,dim);
Logger<Primary> log;
double tic = Logger<Primary>::GetTime();
// {
// std::cout<<"\tVersion with Optimized Least Square\n\n";
// Block<Scalar,Primary> X0(SizeBlock,dim);
// Block<Scalar,Primary> sol(SizeBlock,dim);
// Logger<Primary> log;
// double tic = Logger<Primary>::GetTime();
nb3 = BGMRes<Arnoldi3,Matrix,Scalar,Primary>(mat,RHS,X0,max_ite,restart,sol,log,
0.0001);
// nb3 = BGMRes<Arnoldi3,Matrix,Scalar,Primary>(mat,RHS,X0,max_ite,restart,sol,log,
// 0.0001);
elapsed3 = Logger<Primary>::GetTime() - tic;
// elapsed3 = Logger<Primary>::GetTime() - tic;
std::cout<<"Return value / max_ite : "<<nb3<<"/"<<max_ite<<"\n";
std::string locationRes("../data/res/MatCone_QR.res");
std::ofstream fileRes(locationRes,std::ios::out);
log.forEachIte([&](int ite, int nbMatVect, int currSize,
Primary min, Primary max, double time)
{
fileRes<<ite<<"\t"
<<nbMatVect<<"\t"
<<currSize<<"\t"
<<min<<"\t"
<<max<<"\t"
<<time<<"\n";
});
minmax3 = CompareBlocks<Scalar,Primary>(sol,XExact);
}
// std::cout<<"Return value / max_ite : "<<nb3<<"/"<<max_ite<<"\n";
// std::string locationRes("../data/res/MatCone_QR.res");
// std::ofstream fileRes(locationRes,std::ios::out);
// log.forEachIte([&](int ite, int nbMatVect, int currSize,
// Primary min, Primary max, double time)
// {
// fileRes<<ite<<"\t"
// <<nbMatVect<<"\t"
// <<currSize<<"\t"
// <<min<<"\t"
// <<max<<"\t"
// <<time<<"\n";
// });
// minmax3 = CompareBlocks<Scalar,Primary>(sol,XExact);
// }
// {
// {
// //Compute what i need for my metrics
// //MatxSol : Mat*Sol_computed
// Block<Scalar,Primary> MatxSol(SizeBlock,dim);
......
......@@ -225,34 +225,34 @@ int main(int ac, char ** av){
double elapsed1,elapsed2,elapsed3;
int nb1,nb2,nb3;
std::pair<Primary,Primary> minmax1,minmax2,minmax3;
// {
// std::cout<<"\n\tVersion Standard\n\n";
// Block<Scalar,Primary> X0(SizeBlock,dim);
// Block<Scalar,Primary> sol(SizeBlock,dim);
// Logger<Primary> log;
// double tic = Logger<Primary>::GetTime();
// nb1 = BGMRes<Arnoldi1,Matrix,Scalar,Primary>(mat,RHS,X0,nbIte,nbIte,sol,log,
// 0.00001);
// elapsed1 = Logger<Primary>::GetTime() - tic;
// std::cout<<"Return value / max_ite : "<<nb1<<"/"<<nbIte<<"\n";
// std::string locationRes("../data/res/conf5_0_00l4x4_1000.res");
// std::ofstream fileRes(locationRes,std::ios::out);
// log.forEachIte([&](int ite, int nbMatVect, int currSize,
// Primary min, Primary max, double time)
// {
// fileRes<<ite<<"\t"
// <<nbMatVect<<"\t"
// <<currSize<<"\t"
// <<min<<"\t"
// <<max<<"\t"
// <<time<<"\n";
// });
// minmax1 = CompareBlocks<Scalar,Primary>(sol,XExact);
// }
{
std::cout<<"\n\tVersion Standard\n\n";
Block<Scalar,Primary> X0(SizeBlock,dim);
Block<Scalar,Primary> sol(SizeBlock,dim);
Logger<Primary> log;
double tic = Logger<Primary>::GetTime();
nb1 = BGMRes<Arnoldi1,Matrix,Scalar,Primary>(mat,RHS,X0,nbIte,nbIte,sol,log,
0.00001);
elapsed1 = Logger<Primary>::GetTime() - tic;
std::cout<<"Return value / max_ite : "<<nb1<<"/"<<nbIte<<"\n";
std::string locationRes("../data/res/conf5_0_00l4x4_1000.res");
std::ofstream fileRes(locationRes,std::ios::out);
log.forEachIte([&](int ite, int nbMatVect, int currSize,
Primary min, Primary max, double time)
{
fileRes<<ite<<"\t"
<<nbMatVect<<"\t"
<<currSize<<"\t"
<<min<<"\t"
<<max<<"\t"
<<time<<"\n";
});
minmax1 = CompareBlocks<Scalar,Primary>(sol,XExact);
}
{
std::cout<<"\tVersion with Inexact Breakdown\n\n";
Block<Scalar,Primary> X0(SizeBlock,dim);
......
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