Commit 3e81442a authored by Cyrille Piacibello's avatar Cyrille Piacibello

Seems to be working, need to be tested with Real matrix

parent 094906ea
......@@ -13,15 +13,15 @@ private:
int ldBase; //Leading dimension
Scalar * datas; //raw ptr
int totalNbBlock; //Number of block of vectors
int * nbVectInBlock; //Index of starting block (block size may change)
int * sumNbVectInBlock; //Index of starting block (block size may change)
int currentBlockUsed; //nbblock used
//Use to modify the cursor currentBlockUsed and fill corresponding
//entry in nbVectInBlock
//entry in sumNbVectInBlock
void addBlock(int inputSizeBlock){
currentBlockUsed++;
nbVectInBlock[currentBlockUsed] = nbVectInBlock[currentBlockUsed-1]+ inputSizeBlock;
sumNbVectInBlock[currentBlockUsed] = sumNbVectInBlock[currentBlockUsed-1]+ inputSizeBlock;
}
......@@ -31,14 +31,14 @@ public:
ldBase(sizeVector),
datas(nullptr),
totalNbBlock(restart),
nbVectInBlock(nullptr),
sumNbVectInBlock(nullptr),
currentBlockUsed(0)
{
totalSize = ldBase*restart*SizeBlock;
datas = new Scalar[totalSize];
memset(datas,0,totalSize*sizeof(Scalar));
nbVectInBlock = new int[totalNbBlock+1];
memset(nbVectInBlock,0,sizeof(int)*(totalNbBlock+1));
sumNbVectInBlock = new int[totalNbBlock+1];
memset(sumNbVectInBlock,0,sizeof(int)*(totalNbBlock+1));
// Initializing Base
std::cout<<"Init Base with :\n"
<<"\ttotalSize "<<totalSize<<"\n"
......@@ -52,13 +52,13 @@ public:
~Base(){
delete [] datas;
delete [] nbVectInBlock;
delete [] sumNbVectInBlock;
}
Scalar * getPtr(int idBlock = 0){
//starting point is in nbVectInBlock[idBlock]
return &datas[nbVectInBlock[idBlock]*ldBase];
//starting point is in sumNbVectInBlock[idBlock]
return &datas[sumNbVectInBlock[idBlock]*ldBase];
}
int getLeadingDim(){
......@@ -71,7 +71,9 @@ public:
int getNbBlockUsed(){
return currentBlockUsed+1;
}
int getNbVectUsed(){
return sumNbVectInBlock[currentBlockUsed];
}
/**
* @brief return the last block col
*/
......@@ -98,28 +100,28 @@ public:
* @param HessToWrite Array of Scalar to write the
* orthogonalization coeffs
* @param ldHess : Leading dimension of Hessenberg
* @param idxToStart : number of vectors to discard in Ortho process
*/
void ComputeMGSOrtho(Block<Scalar,Primary> & Wj, Scalar * HessToWrite, int ldHess){
void ComputeMGSOrtho(Block<Scalar,Primary> & Wj, Scalar * HessToWrite, int ldHess,
int idxToStart = 0){
Scalar * ptrOfHess = HessToWrite;
int nbColinWj = Wj.getSizeBlock() - idxToStart;
//Loop over different blocks
for(int i=0 ; i< currentBlockUsed ; ++i){
int sizeCurrBlock = nbVectInBlock[i+1]-nbVectInBlock[i];
int sizeCurrBlock = sumNbVectInBlock[i+1]-sumNbVectInBlock[i];
std::cout<<"Ortho against block number "<<i<<" of size "
<<sizeCurrBlock<<"\n"
<<"param to gemm : \t"<<sizeCurrBlock<<"\t"<<Wj.getSizeBlock()
<<"\t"<<ldBase<<"\n";
<<sizeCurrBlock<<"\n";
//Compute hessenberg coeff
call_Tgemm<Scalar>(sizeCurrBlock,Wj.getSizeBlock(), ldBase,
call_Tgemm<Scalar>(sizeCurrBlock,nbColinWj, ldBase,
getPtr(i),ldBase,
Wj.getPtr(),Wj.getLeadingDim(),
Wj.getPtr(idxToStart),Wj.getLeadingDim(),
ptrOfHess, ldHess, 1, 0);
//Modify Wj
//Wj = (-1) * V_i * Hij + Wj;
call_gemm<Scalar>(ldBase,Wj.getSizeBlock(),sizeCurrBlock,
call_gemm<Scalar>(ldBase,nbColinWj,sizeCurrBlock,
getPtr(i),ldBase,
ptrOfHess, ldHess,
Wj.getPtr(),Wj.getLeadingDim(),
Wj.getPtr(idxToStart),Wj.getLeadingDim(),
Scalar(-1), Scalar(1));
ptrOfHess += sizeCurrBlock;
}
......@@ -130,7 +132,10 @@ public:
*
*/
void ComputeProduct(Block<Scalar,Primary> & input, Block<Scalar,Primary>& output){
int minLine = std::min(nbVectInBlock[currentBlockUsed],input.getLeadingDim());
int minLine = std::min(sumNbVectInBlock[currentBlockUsed],input.getLeadingDim());
std::cout<<"Computing Xm solution\t: gemm inputs : "<<ldBase
<<" "<<input.getSizeBlock()
<<" "<<minLine<<"\n";
call_gemm<Scalar>(ldBase,input.getSizeBlock(),minLine,
getPtr(),ldBase,
input.getPtr(),input.getLeadingDim(),
......@@ -140,7 +145,7 @@ public:
void displayBase(std::string name = ""){
std::cout<<name<<"\t"<<"currentBlockUsed "<<currentBlockUsed<<"\n";
for(int i = 0 ; i< nbVectInBlock[currentBlockUsed] ; ++i){
for(int i = 0 ; i< sumNbVectInBlock[currentBlockUsed] ; ++i){
for(int j=0 ; j< getLeadingDim() ; ++j){
std::cout<<datas[i*ldBase+j]<<"\t";
}
......@@ -148,6 +153,18 @@ public:
}
}
void checkOrtho(){
int nbVect = sumNbVectInBlock[currentBlockUsed];
Block<Scalar,Primary> A(nbVect,nbVect);
call_Tgemm(nbVect,nbVect, ldBase,
getPtr(),ldBase,
getPtr(),ldBase,
A.getPtr(),A.getLeadingDim(),Scalar(1),Scalar(0));
std::cout<<"Base*Base^{H} : \n";
for(int i=0 ; i<nbVect ; ++i){
std::cout<<"\tNorm of ["<<i<<"]\t"<<A.getNorm(i)<<"\n";
}
}
};
......
......@@ -85,6 +85,7 @@ public:
//QR facto
void ComputeQRFacto(Block& Q, Block& R, //input are already allocated
int hStart=0, int vStart=0){ //if QR of part of block
//Sizes needed
int sizeQrBlock = SizeBlock - vStart;
......@@ -93,11 +94,6 @@ public:
memset(tau,0,sizeof(Scalar)*sizeQrBlock);
//Qr facto can be done in place
std::cout<<"Parameters to call_geqrf : "
<<"\nldb-hStart:\t"<<ldb-hStart
<<"\nsizeQrBlock:\t"<<sizeQrBlock
<<"\nStarting point:\t"<<hStart+vStart*ldb
<<std::endl;
int geqrf = callLapack_geqrf(ldb-hStart,sizeQrBlock,
&datas[hStart+vStart*ldb],ldb,tau);
if(geqrf != 0){
......@@ -115,22 +111,34 @@ public:
//use the full QR
int orgqr;
if(Q.getSizeBlock() == Q.getLeadingDim()){
std::cout<<"Generating full Q after QR\n";
std::cout<<"Generating full Q after QR\n"<<std::endl;
//First copy inside Q becaus current block does not have
//necessary space for a full QR facto
for(int i=0 ; i<sizeQrBlock ; ++i){
memcpy(Q.getPtr(i),&datas[hStart+(vStart+i)*ldb],
(ldb-hStart)*sizeof(Scalar));
}
orgqr = callLapack_orgqr(ldb-hStart,sizeQrBlock,
&datas[hStart+vStart*ldb],ldb,tau,true);
Q.getPtr(),Q.getLeadingDim(),
tau,true);
if(orgqr != 0){
std::cout<<"Problem on Q generation after QR facto\n";
}
}else{ //Else, only leading p directions are computed
std::cout<<"Generating partial Q after QR\n";
std::cout<<"Generating partial Q after QR\n"<<std::endl;
orgqr = callLapack_orgqr(ldb-hStart,sizeQrBlock,
&datas[hStart+vStart*ldb],ldb,tau,false);
if(orgqr != 0){
std::cout<<"Problem on Q generation after QR facto\n";
}
//Copy output to Q
for(int i=0 ; i<Q.getSizeBlock() ; ++i){
memcpy(Q.getPtr(i),&datas[hStart+(vStart+i)*ldb],
(ldb-hStart)*sizeof(Scalar));
}
}
if(orgqr != 0){
std::cout<<"Problem on Q generation after QR facto\n";
}
//Copy output to Q
for(int i=0 ; i<Q.getSizeBlock() ; ++i){
memcpy(Q.getPtr(i),&datas[hStart+(vStart+i)*ldb],
(ldb-hStart)*sizeof(Scalar));
}
delete [] tau;
}
......@@ -153,7 +161,6 @@ public:
std::pair<Primary,Primary> minmax{100,0};
for(int i=0 ; i<SizeBlock ; ++i){
Primary norm = getNorm(i)/normB[i];
std::cout<<"NormRes["<<i<<"]\t"<<norm<<"\n";
if(norm > minmax.second){
minmax.second = norm;
}
......@@ -224,10 +231,10 @@ public:
id++;
}
//Display all the singular values
for(int i=0 ; i<getSizeBlock() ; ++i){
std::cout<<"sigma["<<i<<"]\t : "<<sigma[i]<<"\n";
}
std::cout<<"Number of sv sup. to "<<epsilon<<"\t: "<<id<<"\n";
// for(int i=0 ; i<getSizeBlock() ; ++i){
// std::cout<<"sigma["<<i<<"]\t : "<<sigma[i]<<"\n";
// }
//std::cout<<"Number of sv sup. to "<<epsilon<<"\t: "<<id<<"\n";
output.initData(id,getLeadingDim());
memcpy(output.getPtr(),U.getPtr(),sizeof(Scalar)*id*U.getLeadingDim());
delete [] sigma;
......
......@@ -12,10 +12,13 @@
template<class Matrix,typename Scalar,typename Primary=Scalar>
int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
int max_iter, Block<Scalar,Primary>& Sol,
Primary epsilon = 0.00001){
Primary epsilon = 0.0000001){
//For monitoring purpose
int nbDeflDir = 0;
int SizeBlock = B.getSizeBlock();
int dim = A.size();
//store norms
std::vector<Primary> normB;
for(int i=0 ; i< SizeBlock ; ++i){
......@@ -57,7 +60,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);
......@@ -67,28 +70,59 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
int j;
//Main Loop
for(j=1 ; j<=max_iter ; ++j){
std::cout<<"################## Start of Iteration #################\n";
//Compute block Wj inside WP
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");
Scalar * toWrite = L.getNewStartingCol(WP.getSizeW()); //Where to write
base.ComputeMGSOrtho(WP,toWrite,L.getLeadingDim());
base.ComputeMGSOrtho(WP,toWrite,L.getLeadingDim(),WP.getSizeP());
std::cout<<"Ortho against V done"<<std::endl;
std::cout<<"Ortho against V and P"<<std::endl;
//Orthogonalisation of Wj against Pj
WP.OrthoWP(L.getPtrToC(),L.getLeadingDim());
std::cout<<"Ortho against of W against Pj"<<std::endl;
//QR of Wj --> {Wj~,Dj}
WP.QRofW(L.getPtrToD(),L.getLeadingDim());
//L.template displayHessExtendedBitMap<Primary>("L after Orthogonalization");
std::cout<<"QR of W done"<<std::endl;
//Create Block and Residual
Block<Scalar,Primary> Res(SizeBlock,L.getNbLineUsed() + SizeBlock);
//Compute residual
L.ComputeResidual(R1,Y,Res);
std::cout<<"Residual Computed !"<<std::endl;
std::cout<<"\nResidual Computed ...\n";
// //Here we compute Base*Y +X0 and display res
// {
// //if no precond
// // X0 = X0 + Base*Yj
// Block<Scalar,Primary> Resultat{SizeBlock,dim};
// base.ComputeProduct(Y,Resultat);
// for(int i=0 ; i<SizeBlock ; ++i){
// for(int j=0 ; j<dim ; ++j){
// Resultat.getPtr(i)[j] += X0.getPtr(i)[j];
// }
// }
// //Resultat.displayBlock("Current Xm");
// //Compute A*Xm-B, and check norm
// Block<Scalar,Primary> AxSol(SizeBlock,dim);
// A.MatBlockVect(Resultat,AxSol);
// //Compute real residual
// for(int i=0 ; i<SizeBlock ; ++i){
// for(int j=0 ; j<dim ; ++j){
// AxSol.getPtr(i)[j] -= B.getPtr(i)[j];
// }
// }
// for(int i=0 ; i<SizeBlock ; ++i){
// std::cout<<"AAA "<<i<<"\t"<<j<<"\t"<<AxSol.getNorm(i)<<"\n";
// }
// }
//Do we really need that ?
std::pair<Primary,Primary> MinMax = Res.getMinMaxNorm(normB);
std::cout<<"Ite ["<<j<<"] Done\t"<<
......@@ -96,7 +130,7 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
"Max\t"<<MinMax.second<<"\n";
if(MinMax.second < epsilon){
std::cout<<"Convergence ?"<<std::endl;
std::cout<<"Convergence !!"<<std::endl;
break;
}
......@@ -105,13 +139,16 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
Res.DecompositionSVD(U,epsilon_R);
std::cout<<"Decomposition SVD done !\n\t\tSingular Values kept "
<<U.getSizeBlock()<<" !"<<std::endl;
//U.displayBlock("U: Kept directions");
//Here : If no inexact breakdowns occured, updates procedures are easiers
if(U.getSizeBlock() == SizeBlock){
//Augment base and move on
//Augment base : [P_{j-1},~Wj] = [~Wj]
base.addBlockDatas(WP);
}else{
if(U.getSizeBlock() + nbDeflDir != SizeBlock){
nbDeflDir++;
}
std::cout<<"### Inexact Breakdowns Occured ###!\t\t Ite : \t"<<j<<"\n";
//Create Block to store \W
Block<Scalar,Primary> Directions(SizeBlock,SizeBlock);
//TO BE DISCARDED
......@@ -119,23 +156,25 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0
//Qr facto of bottom block of U
U.ComputeQRFacto(Directions,tempRpart,sumOfSizeOfEachBlock.back());
Directions.displayBlock("Ortho part of U after QR");
//Directions.displayBlock("Ortho part of U after QR");
//Update
//V_{j+1} = [P_{j-1},~Wj] * Directions_1
WP.ComputeProduct(Directions,U.getSizeBlock(),
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_{j+1,} = | Directions_1^{H} | * H_j
//Gj = | Directions_2^{H} |
L.UpdateBottomLine(Directions,U.getSizeBlock());
L.UpdateBottomLine(Directions);
//L.template displayHessExtendedBitMap<Primary>("L after update");
}
sizeOfEachBlock.push_back(U.getSizeBlock());
sumOfSizeOfEachBlock.push_back(U.getSizeBlock()+sumOfSizeOfEachBlock.back());
std::cout<<"################## End of Iteration ##################\n"
<<"Base contains exactly :: "<<base.getNbVectUsed()<<"\t!!\n";
}
//sol = Base*Y + X0
base.ComputeProduct(Y,Sol);
......
......@@ -93,10 +93,34 @@ int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
Wj.ComputeQRFacto(Q,R);
//Copy Block R inside Hess
Hess.addRPartToHess(R);
if(j>7){//last ite
std::cout<<"Iteration \t\t\t\t"<<j<<"§§§§§§§§§§§§§\n";
R.displayBlock("Last R ?");
}
std::cout<<"QR done"<<std::endl;
//Tests over criterion
Hess.ComputeBlockResidual(R1,Yj);
Yj.displayBlock("Y Computed");
//Here we compute Base*Y +X0 and display res
{
//if no precond
// X0 = X0 + Base*Yj
Block<Scalar,Primary> Resultat{SizeBlock,dim};
base.ComputeProduct(Yj,Resultat);
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ; j<dim ; ++j){
Resultat.getPtr(i)[j] += X0.getPtr(i)[j];
}
}
Resultat.displayBlock("Current Xm");
for(int i=0 ; i<SizeBlock ; ++i){
std::cout<<"PREFIX "<<i<<"\t"<<j<<"\t"<<Resultat.getNorm(i)<<"\n";
}
}
std::cout<<"["<<j<<"]\t"<<"Residual computed"<<std::endl;
{//Display residual to see convergence
......
......@@ -41,6 +41,8 @@ public:
*/
void OrthoWP(Scalar * toWrite, int ldToWrite){
if(cursor>0){
std::cout<<"Ortho of "<<getSizeW()<<" vectors in W against "<<getSizeP()
<<" vectors in P\n";
//Generate coef inside C
call_Tgemm<>(getSizeP(),getSizeW(),Parent::getLeadingDim(),
getP(),Parent::getLeadingDim(),
......
......@@ -111,7 +111,6 @@ public:
memset(r,0,R1.getSizeBlock()*SizeBlock*(currentSize+1)*sizeof(Scalar));
int ldr = (currentSize+1)*SizeBlock;
std::cout<<"R : is "<<R1.getSizeBlock()<<" x "<<ldr<<"...\n";
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0; j<=i ; ++j){
r[ldr*i+j] = R1.getPtr(i)[j];
......@@ -140,7 +139,7 @@ public:
}
}
Yj.displayBlock("Y resulting of Gels");
//Yj.displayBlock("Y resulting of Gels");
delete [] workingCopy;
delete [] r;
}
......@@ -171,8 +170,6 @@ public:
std::pair<Primary,Primary> MinMax = std::make_pair<>(1000,0);
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;
}
......
......@@ -76,6 +76,26 @@ public:
}
std::cout<<std::endl;
}
template<class Primary>
void displayHessExtendedBitMap(std::string name = ""){
std::cout<<name<<"\n";
for(int j=0 ; j<nbLineUsed+p ; ++j){
int col = 0;
for(int i=0 ; i<nbVectInBlck.size() ; ++i){
for(int t=0; t<nbVectInBlck[i] ; ++t){
if(module<Primary,Scalar>(data[col*ldHess+j]) != Primary(0)){
std::cout<<"."<<" ";
}else{
std::cout<<" "<<" ";
}
col++;
}
std::cout<<"| ";
}
std::cout<<std::endl;
}
std::cout<<std::endl;
}
/**
* @brief return the ptr to the start of new blockcol, and augment
......@@ -87,7 +107,6 @@ public:
sumNbVectInBlck.push_back(ref+inSizeBlock);
nbBlckColUsed++;
nbLineUsed += nbVectInBlck.back();
std::cout<<"New Starting col is in data["<<ref*ldHess<<"]\n";
return &data[ref*ldHess];
}
......@@ -103,7 +122,7 @@ public:
* @param sizeToRead : limite between W1 and W2
*/
template<typename Primary>
void UpdateBottomLine(Block<Scalar,Primary>& W,int sizeToRead){
void UpdateBottomLine(Block<Scalar,Primary>& W){
int nj = sumNbVectInBlck.back();
//Save H in order to write directly inside L
Block<Scalar,Primary> Hj(nj,p);
......@@ -129,15 +148,9 @@ public:
Block<Scalar,Primary>& outputY,
Block<Scalar,Primary>& outputResidual){
//Complete gamma
std::cout<<"About to compute Residual...\n";
std::cout<<"Gamma : is "<<gamma.getSizeBlock()<<" x "
<<gamma.getLeadingDim()<<"\n";
std::cout<<"Nb Line used : "<<nbLineUsed<<"\n";
int ldR = nbLineUsed + p;
Scalar * r = new Scalar[ldR*p];
memset(r,0,((nbLineUsed + p)*p)*sizeof(Scalar));
std::cout<<"R : is "<<p<<" x "<<ldR<<"...\n";
for(int i=0 ; i<p ; ++i){
memcpy(&r[i*ldR],gamma.getPtr(i),
sizeof(Scalar)*(i+1));
......@@ -161,11 +174,8 @@ public:
for(int i=0 ; i<p ; ++i){
memcpy(outputY.getPtr(i),&r[i*ldR],sizeof(Scalar)*ldR);
}
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 : outputRes = F*Y-gama
call_gemm<>(nbLineUsed + p,p,sumNbVectInBlck.back(),
data,ldHess,
outputY.getPtr(),outputY.getLeadingDim(),
......@@ -177,7 +187,7 @@ public:
outputResidual.getPtr(i)[j] -= gamma.getPtr(i)[j];
}
}
outputResidual.displayBlock("Residual obtained from LeastSquare");
delete [] workingCopy;
delete [] r;
}
......@@ -194,7 +204,7 @@ public:
*
*/
Scalar * getPtrToG(){
// std::cout<<"\tCall to getPtrToG():\n\t return data["<<nbLineUsed<<"]\n";
//std::cout<<"\tCall to getPtrToG():\n\t return data["<<nbLineUsed<<"]\n";
return &data[nbLineUsed];
}
......
......@@ -326,8 +326,6 @@ int callLapack_orgqr(int m, int p,
std::complex<double> * mat, int ldMat,
std::complex<double> * tau, bool full){
if(full){
std::cout<<"There !\n"<<
"Dans l'ordre : m "<<m<<" p "<<p<<"\n";;
return LAPACKE_zungqr(LAPACK_COL_MAJOR, m, m, p, mat, ldMat,tau);
}else{
return LAPACKE_zungqr(LAPACK_COL_MAJOR, m, p, p, mat, ldMat,tau);
......
......@@ -7,6 +7,7 @@ set(GMRES_TEST_SRC
testBlockSvd.cpp
testBlockWP.cpp
testBlockArnoldiIb.cpp
testMatrixMarket.cpp
)
#set(LIBS_FOR_TESTS ibgmresdr)
......
......@@ -120,6 +120,7 @@ int main(int ac, char ** av){
}
int dim = nbBlock*SizeBlock;
int maxiteIB = nbBlock + (SizeBlock-1);
using Primary=double;
using Scalar=std::complex<double>;
......@@ -162,17 +163,19 @@ int main(int ac, char ** av){
Block<Scalar,Primary> sol(SizeBlock,dim);
Block<Scalar,Primary> sol_ref(SizeBlock,dim);
display_matrix_matlab(mat.getPtr(),dim,dim);
std::cout<<"\n\n";
display_matrix_matlab(RHS.getPtr(),dim,SizeBlock);
std::cout<<"\n\n";
display_matrix_matlab(X0.getPtr(),dim,SizeBlock);
// display_matrix_matlab(mat.getPtr(),dim,dim);
// std::cout<<"\n\n";
// display_matrix_matlab(RHS.getPtr(),dim,SizeBlock);
// std::cout<<"\n\n";
// display_matrix_matlab(X0.getPtr(),dim,SizeBlock);
std::cout<<"\n\n";
Logger<Primary> logger;
int nb = BlockArnoldiIb<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,
X0,5,sol);
X0,
maxiteIB,
sol);
std::cout<<"Return value : "<<nb<<"\n";
......@@ -213,49 +216,61 @@ 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 : "<<nb2<<"\n";
//Check solution
{
//Compute what i need for my metrics
//MatxSol : Mat*Sol_computed
Block<Scalar,Primary> MatxSol(SizeBlock,dim);
mat.MatBlockVect(sol_ref,MatxSol);
//MatxSolMinusB : Mat*Sol_computed-B
Block<Scalar,Primary> MatxSolMinusB(SizeBlock,dim);
//X_Diff : X_Exact - Sol_computed
Block<Scalar,Primary> X_Diff(SizeBlock,dim);
MatxSolMinusB.CopyBlock(MatxSol);
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
for(int j=0 ; j<dim; ++j){
MatxSolMinusB.getPtr(idxRHS)[j] -= RHS.getPtr(idxRHS)[j];
X_Diff.getPtr(idxRHS)[j] = X_Exact.getPtr(idxRHS)[j] -
sol_ref.getPtr(idxRHS)[j];
}
}
//Then display for each vect
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
std::cout<<"Metrics on "<<idxRHS<<"\n";
Primary ForwardError = X_Diff.getNorm(idxRHS);
std::cout <<"["<<idxRHS<<"]\t"
<<"Forward Error (on X) \t" << ForwardError <<"\n";
Primary SecondMemberError = MatxSolMinusB.getNorm(idxRHS);
std::cout <<"["<<idxRHS<<"]\t"
<<"Second Member Error \t" << SecondMemberError <<"\n";
Primary ForwardErrorNormalized = ForwardError/(X_Exact.getNorm(idxRHS));
std::cout <<"["<<idxRHS<<"]\t"
<<"Forward Error Norm. \t" << ForwardErrorNormalized <<"\n";
Primary normB = RHS.getNorm(idxRHS);
Primary NormwiseBackwardError = SecondMemberError / normB;
std::cout <<"["<<idxRHS<<"]\t"
<<"Normwise Backward Error\t" << NormwiseBackwardError <<"\n\n";
}
}
// int nb2 = BlockGMRes<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,X0,9,
// sol_ref,9,
// logger);
// std::cout<<"Return value : "<<nb2<<"\n";
// //Check solution
// {
// //Compute what i need for my metrics
// //MatxSol : Mat*Sol_computed
// Block<Scalar,Primary> MatxSol(SizeBlock,dim);
// mat.MatBlockVect(sol_ref,MatxSol);
// //MatxSolMinusB : Mat*Sol_computed-B
// Block<Scalar,Primary> MatxSolMinusB(SizeBlock,dim);
// //X_Diff : X_Exact - Sol_computed
// Block<Scalar,Primary> X_Diff(SizeBlock,dim);
// MatxSolMinusB.CopyBlock(MatxSol);
// for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
// for(int j=0 ; j<dim; ++j){
// MatxSolMinusB.getPtr(idxRHS)[j] -= RHS.getPtr(idxRHS)[j];
// X_Diff.getPtr(idxRHS)[j] = X_Exact.getPtr(idxRHS)[j] -
// sol_ref.getPtr(idxRHS)[j];
// }
// }
// //Then display for each vect
// for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
// std::cout<<"Metrics on "<<idxRHS<<"\n";
// Primary ForwardError = X_Diff.getNorm(idxRHS);
// std::cout <<"["<<idxRHS<<"]\t"
// <<"Forward Error (on X) \t" << ForwardError <<"\n";
// Primary SecondMemberError = MatxSolMinusB.getNorm(idxRHS);
// std::cout <<"["<<idxRHS<<"]\t"
// <<"Second Member Error \t" << SecondMemberError <<"\n";
// Primary ForwardErrorNormalized = ForwardError/(X_Exact.getNorm(idxRHS));
// std::cout <<"["<<idxRHS<<"]\t"
// <<"Forward Error Norm. \t" << ForwardErrorNormalized <<"\n";
// Primary normB = RHS.getNorm(idxRHS);
// Primary NormwiseBackwardError = SecondMemberError / normB;
// std::cout <<"["<<idxRHS<<"]\t"
// <<"Normwise Backward Error\t" << NormwiseBackwardError <<"\n\n";
// }
// }
// std::cout<<"Difference betweens solutions found\n";
// Block<Scalar,Primary> Sol_Diff(SizeBlock,dim);
// for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
// for(int j=0 ; j<dim; ++j){
// Sol_Diff.getPtr(idxRHS)[j] = sol_ref.getPtr(idxRHS)[j] -
// sol.getPtr(idxRHS)[j];
// }
// }
// for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
// std::cout<<"Norm sol-sol_res["<<idxRHS<<"]\t"
// <<Sol_Diff.getNorm(idxRHS)<<"\n";
// }
return 0;
}
......@@ -89,7 +89,7 @@ int main(int ac, char ** av){
for(auto ag : args)
std::cout << ag << "\n";
int SizeBlock=1,nbBlock=4;
int SizeBlock=4,nbBlock=5;
if (args.size()>1){
SizeBlock = std::stoi(args[1]);
......
#include <random>
#include <iostream>
#include <fstream>
#include "../src/Block.hpp"
#include "../src/BlockGMRes.hpp"
#include "../src/LapackInterface.hpp"
#include "../src/Logger.hpp"
#include "../src/BlockArnoldiIb.hpp"
/**
* Matrix Class : need to implement MatBlockVect(), MatBaseproduct(),
* useRightPreCond(), preCondBlockVect(),preCondBaseProduct() and
* size() member function in order to works with BlockGMRes algorithm
*
* An abstract class could be used to enforce that.
*/
template<class Scalar,class Primary=Scalar>
class InputMatrix{
int dim;
Scalar * data;
public:
InputMatrix() : dim(0),data(nullptr){
}
InputMatrix(int inSize) : dim(inSize),data(nullptr){