Commit a16196b9 authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

Corresponding Tests added

parent dd4f1493
#ifndef HESS_HPP
#define HESS_HPP
#include "LapackInterface.hpp"
#include "Block.hpp"
/**
* This class is supposed to handle Hessenberg in order to provide
* directly a ptr to Lapack...
*
* Format:
* | X X X |
* | X X X |
* | X X |
* | X |
*
* @attention : the whole matrix is allocated, and its size is driven
* by the restart parameter.
*/
template<typename Scalar>
class Hessenberg{
//Depends on restart
int totalSize;
int restart;
int SizeBlock;
int currentSize; //Nb of block Column used (if regular gmres,
//number of column)
int * nbVectInBlock; //Index of starting block (block size may change)
int ldh;
Scalar * datas;
public:
Hessenberg(int restart, int SizeBlock) : totalSize(0),restart(restart),
SizeBlock(SizeBlock), currentSize(0),
nbVectInBlock(nullptr),
ldh(SizeBlock*(restart+2)),
datas(nullptr){
totalSize = (restart+1)*(restart+2)*SizeBlock*SizeBlock;
datas = new Scalar[totalSize];
memset(datas,0,totalSize*sizeof(Scalar));
nbVectInBlock = new int[restart+2];
memset(nbVectInBlock,0,sizeof(int)*(restart+2));
// Initializing Hessenberg
std::cout<<"Init Hess with :\n"
<<"\ttotalSize "<<totalSize<<"\n"
<<"\trestart "<<restart<<"\n"
<<"\tSizeBlock "<<SizeBlock<<"\n"
<<"\tldh "<<ldh<<"\n";
}
~Hessenberg(){
delete [] datas;
delete [] nbVectInBlock;
datas = nullptr;
}
/**
* Return the raw ptr
*/
Scalar * getPtr(int idBlock = 0){
return &datas[nbVectInBlock[idBlock]*ldh];
}
int getLeadingDim(){
return ldh;
}
Scalar & operator[](int ind){
return datas[ind];
}
int getCurrentSize(){
return currentSize;
}
int getNbColUsed(){
return nbVectInBlock[currentSize];
}
//Here add the last triang block
template<typename Primary>
void addRPartToHess(Block<Scalar,Primary> & R){
int inSizeBlock = R.getSizeBlock();
//Get the ptr where we currently are
Scalar * toWrite = getPtr(currentSize) + inSizeBlock*(currentSize+1);
for(int i=0 ; i< inSizeBlock ; ++i){
memcpy(&toWrite[i*ldh],R.getPtr(i),inSizeBlock*sizeof(Scalar));
}
currentSize++;
//Sum prefix
nbVectInBlock[currentSize] = inSizeBlock + nbVectInBlock[currentSize-1];
}
void displayHess(std::string name = ""){
std::cout<<name<<"\t"<<"Current Size "<<currentSize<<"\n";
for(int i=0 ; i<currentSize*SizeBlock ; ++i){
for(int j=0; j<(currentSize+1)*SizeBlock ; ++j){
std::cout<<getPtr()[i*ldh+j]<<"\t";
}
std::cout<<"\n";
}
}
//Call to lapack here
template<typename Primary>
void ComputeBlockResidual(Block<Scalar,Primary>& R1, Block<Scalar,Primary> & Yj){
//First format R...
Scalar * r = new Scalar[R1.getSizeBlock() * (currentSize+1)*SizeBlock];
memset(r,0,R1.getSizeBlock()*SizeBlock*(currentSize+1)*sizeof(Scalar));
int ldr = (currentSize+1)*SizeBlock;
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0; j<=i ; ++j){
r[ldr*i+j] = R1.getPtr(i)[j];
}
}
Scalar*workingCopy = new Scalar[SizeBlock*SizeBlock*
(currentSize+1)*currentSize];
for(int i=0 ; i<currentSize ; ++i){
for(int j=0 ; j<SizeBlock ; ++j){
memcpy(&workingCopy[i*SizeBlock*ldr + j*ldr], //Dest
&getPtr()[i*SizeBlock*ldh + j*ldh], //Src
sizeof(Scalar)*ldr);
}
}
callLapack_gels<Scalar>(ldr,currentSize*SizeBlock,SizeBlock,
workingCopy,ldr,
r,ldr);
//Build return Matrix
//copy inside Yj
for(int i=0 ; i<Yj.getSizeBlock() ; ++i){
for(int j=0 ; j < ldr ; ++j){
Yj.getPtr(i)[j] = r[ldr*i + j];
}
}
//Yj.displayBlock("Y resulting of Gels");
delete [] workingCopy;
delete [] r;
}
//Compute the Hess*input set of vectors, and return min and max norm
template<typename Primary>
std::pair<Primary,Primary> displayBlockResidual(Block<Scalar,Primary>& Y,Block<Scalar,Primary>& R,
std::vector<Primary> normB){
//The Idea is to compute a matrix vector product...
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));
//blockRes - R
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
for(int i = 0 ; i<R.getLeadingDim() ; ++i){
blockRes.getPtr(idxRHS)[i] -= R.getPtr(idxRHS)[i];
}
}
blockRes.displayBlock("Block residual");
//Loop to get the min/max of residuals
std::pair<Primary,Primary> MinMax = std::make_pair<>(1000,0);
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
auto norm = blockRes.getNorm(idxRHS)/normB[idxRHS];
if(norm > MinMax.second){
MinMax.second = norm;
}
if(norm < MinMax.first){
MinMax.first = norm;
}
}
return MinMax;
}
};
#endif //HESS_HPP
......@@ -65,6 +65,19 @@ public:
}
}
void initDiagMatrix(int nbDiagEl){
int j=0;
Scalar El{1,1};
int nbSameEl = dim/nbDiagEl;
while(j<dim){
for(int i=0 ; (i<nbSameEl) && (j<dim) ; ++i){
getPtr()[j*dim+j] = El;
j++;
}
El = Scalar{El.imag()+1,El.real()+1};
}
}
void MatBlockVect(Block<Scalar,Primary>& in, Block<Scalar,Primary>& out){
call_gemm<Scalar>(dim,in.getSizeBlock(),dim,
data,dim,
......@@ -118,7 +131,8 @@ int main(int ac, char ** av){
if (args.size()>2){
nbBlock = std::stoi(args[2]);
}
std::cout<<"Solving AX=B with A: "<<nbBlock*SizeBlock<<"x"<<nbBlock*SizeBlock
<<" and B containing: "<<SizeBlock<<"\n";
int dim = nbBlock*SizeBlock;
int maxiteIB = nbBlock + (SizeBlock-1);
......@@ -156,26 +170,26 @@ int main(int ac, char ** av){
}
}
//fill last col with one vect from solution
int last_indice = SizeBlock-1;
memcpy(X0.getPtr(last_indice),X_Exact.getPtr(0),
sizeof(Scalar)*X0.getLeadingDim());
// int last_indice = SizeBlock-1;
// memcpy(X0.getPtr(last_indice),X_Exact.getPtr(0),
// sizeof(Scalar)*X0.getLeadingDim());
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(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;
Logger<Primary,4> logger;
int nb = BlockArnoldiIb<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,
X0,
maxiteIB,
sol);
int nb = BlockArnoldiIb<InputMatrix<Scalar,Primary>,Scalar,Primary,4>(mat,RHS,
X0,
maxiteIB,
sol,logger);
std::cout<<"Return value : "<<nb<<"\n";
......@@ -216,6 +230,7 @@ int main(int ac, char ** av){
}
}
logger.resetLogs();
// int nb2 = BlockGMRes<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,X0,9,
// sol_ref,9,
// logger);
......
......@@ -29,23 +29,23 @@ public:
}
void initMatrix(std::random_device& rd){
// std::mt19937 gen(rd());
// std::uniform_real_distribution<> dis(0, 2);
// for(int i=0 ; i<dim ; ++i){
// for(int j=0 ; j<dim ; ++j){
// data[dim*i+j] = dis(gen);
// }
// }
Scalar diag = dim;
Scalar no_diag = 1;
std::mt19937 gen(/*rd()*/0);
std::uniform_real_distribution<> dis(0, 2);
for(int i=0 ; i<dim ; ++i){
for(int j=0 ; j<dim ; ++j){
if(i==j)
data[i*dim + j] = diag+j;
else
data[i*dim + j] = no_diag;
data[dim*i+j] = dis(gen);
}
}
// Scalar diag = dim;
// Scalar no_diag = 1;
// for(int i=0 ; i<dim ; ++i){
// for(int j=0 ; j<dim ; ++j){
// if(i==j)
// data[i*dim + j] = diag+j;
// else
// data[i*dim + j] = no_diag;
// }
// }
}
void MatBlockVect(Block<Scalar>& in, Block<Scalar>& out){
......@@ -65,7 +65,7 @@ public:
}
bool useRightPreCond(){
return true;
return false;
}
void preCondBlockVect(Block<Scalar>& input, Block<Scalar>& output){
for(int j=0 ; j<output.getSizeBlock() ; ++j){
......@@ -89,7 +89,7 @@ int main(int ac, char ** av){
for(auto ag : args)
std::cout << ag << "\n";
int SizeBlock=4,nbBlock=5;
int SizeBlock=2,nbBlock=20,restart=15;
if (args.size()>1){
SizeBlock = std::stoi(args[1]);
......@@ -97,7 +97,10 @@ int main(int ac, char ** av){
if (args.size()>2){
nbBlock = std::stoi(args[2]);
}
restart = nbBlock;
if (args.size()>3){
restart = std::stoi(args[3]);
}
int dim = nbBlock*SizeBlock;
using Scalar=double;
......@@ -105,7 +108,7 @@ int main(int ac, char ** av){
//Random engine used to generate matrix
std::random_device rd;
std::mt19937 gen(rd());
std::mt19937 gen(/*rd()*/0);
std::uniform_real_distribution<> dis(0, 2);
//Create matrix
......@@ -125,7 +128,7 @@ int main(int ac, char ** av){
Block<Scalar> RHS(SizeBlock,dim);
mat.MatBlockVect(X_Exact,RHS);
//Compute A starting Block
//Compute a starting Block
Block<Scalar> X0(SizeBlock,dim);
//Fill it!
for(int i=0 ; i<SizeBlock ; ++i){
......@@ -136,11 +139,15 @@ int main(int ac, char ** av){
Block<Scalar> sol(SizeBlock,dim);
Logger<Scalar> log;
Logger<Primary,5> log;
//Three steps
//With Restart
int nb_iter = BlockGMRes<InputMatrix<Scalar>,Scalar>(mat,RHS,
X0,nbBlock,sol,nbBlock,log,
0.0000000001);
X0,nbBlock,sol,
restart,log,
0.000001);
//Check solution
{
//Compute what i need for my metrics
......@@ -178,6 +185,5 @@ int main(int ac, char ** av){
}
}
return 0;
}
This diff is collapsed.
#include <random>
#include <iostream>
#include <fstream>
#include <stdint.h>
#include "../src/Block.hpp"
#include "../src/BlockGMRes.hpp"
#include "../src/LapackInterface.hpp"
#include "../src/Logger.hpp"
#include "../src/BlockArnoldiIb.hpp"
#include "../src/IBGMRes.hpp"
/**
* Driver to read Innovation Works arrays
*/
void readHeader(std::ifstream& file, int & dim1, int& dim2, int & sizeEl){
if(file.is_open()){
int32_t header[5];
file.read(reinterpret_cast<char*>(header),sizeof(int)*5);
std::cout<<"Header : "<<header[0]<<"\t"<<header[1]
<<"\t"<<header[2]<<"\t"<<header[3]<<"\t"<<header[4]<<"\n";
dim1 = header[1];
dim2 = header[2];
switch(header[0]){
case 0:
sizeEl = sizeof(float);
break;
case 1:
sizeEl = sizeof(double);
break;
case 2:
sizeEl = sizeof(std::complex<float>);
break;
case 3:
sizeEl = sizeof(std::complex<double>);
break;
}
}
}
/**
* 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){
data = new Scalar[dim*dim];
}
~InputMatrix(){
delete [] data;
}
int size(){
return dim;
}
Scalar * getPtr(){
return data;
}
void initMatrix(std::random_device& rd){
std::mt19937 gen(/*rd()*/0);
std::uniform_real_distribution<> dis(0, 2);
for(int i=0 ; i<dim ; ++i){
for(int j=0 ; j<dim ; ++j){
data[dim*i+j] = {dis(gen),dis(gen)};
}
}
}
void MatBlockVect(Block<Scalar,Primary>& in, Block<Scalar,Primary>& out,
int idxToWrite=0){
call_gemm<Scalar>(dim,in.getSizeBlock(),dim,
data,dim,
in.getPtr(), in.getLeadingDim(),
out.getPtr(idxToWrite),out.getLeadingDim(),
Scalar(1),Scalar(0));
}
void MatBaseProduct(Scalar * toRead, int nbRHS, Scalar* ptrToWrite){
call_gemm<Scalar>(dim,nbRHS, dim,
data,dim,
toRead,dim,
ptrToWrite,dim,
Scalar(1),Scalar(0));
}
bool useRightPreCond(){
return false;
}
void preCondBlockVect(Block<Scalar,Primary>& input, Block<Scalar,Primary>& output){
for(int j=0 ; j<output.getSizeBlock() ; ++j){
for(int i=0 ; i<dim ; ++i){
output.getPtr(j)[i] = input.getPtr(j)[i]*(Scalar(1)/data[i*dim+i]);
}
}
}
void preCondBaseProduct(Scalar * ptr,Block<Scalar,Primary>& output){
for(int j=0 ; j<output.getSizeBlock() ; ++j){
for(int i=0 ; i<dim ; ++i){
output.getPtr(j)[i] = ptr[j*dim+i]*(Scalar(1)/(data[i*dim+i]));
}
}
}
void initDim(int inDim){
this->dim = inDim;
if(data){
std::cout<<"Error : Data already allocated\nExiting\n";
exit(0);
}
this->data = new Scalar[dim*dim];
}
};
int main(int ac, char ** av){
std::vector<std::string > args(av,av+ac);
for(auto ag : args)
std::cout << ag << "\n";
using Primary=double;
using Scalar=std::complex<double>;
//Random engine used to generate RHS
std::random_device rd;
std::mt19937 gen(/*rd()*/0);
std::uniform_real_distribution<> dis(0, 2);
InputMatrix<Scalar,Primary> mat{};
{
std::string location("../data/matrices_BEM/MatconeSpherePC_MAIN_MAIN_0");
std::ifstream myfile(location,std::ios::binary);
int dim1,dim2,sizeEl;
readHeader(myfile,dim1,dim2,sizeEl);
//check size of element
if(sizeEl != sizeof(Scalar)){
std::cout<<"Different Size between Scalar and sizeEl\nCheck what you doing\nExiting\n";
exit(0);
}
mat.initDim(dim1);
myfile.read(reinterpret_cast<char*>(mat.getPtr()),sizeEl*dim1*dim1);
}
Block<Scalar,Primary> RHS{};
{
std::string location("../data/matrices_BEM/coneSpherePC_RHS_MAIN.0.res");
std::ifstream myfile(location,std::ios::binary);
int dim1,dim2,sizeEl;
readHeader(myfile,dim1,dim2,sizeEl);
//check size of element
if(sizeEl != sizeof(Scalar)){
std::cout<<"Different Size between Scalar and sizeEl\nCheck what you doing\nExiting\n";
exit(0);
}
RHS.initData(dim2,dim1);
myfile.read(reinterpret_cast<char*>(RHS.getPtr()),sizeEl*dim1*dim2);
}
int dim = mat.size();
int SizeBlock = RHS.getSizeBlock();
int maxiteIB = dim / SizeBlock + (SizeBlock);
Block<Scalar,Primary> X0(SizeBlock,dim);
// for(int i=0 ; i<SizeBlock ; ++i){
// for(int j=0 ;j<dim; ++j){
// X0.getPtr(i)[j] = {dis(gen),dis(gen)};
// }
// }
Block<Scalar,Primary> sol(SizeBlock,dim);
Logger<Primary,5> logger;
// int nb = BlockArnoldiIb<InputMatrix<Scalar,Primary>,
// Scalar,Primary,5>(mat,RHS,
// X0,
// maxiteIB,
// sol,logger,
// 0.00001);
int nb = IBGMRes<InputMatrix<Scalar,Primary>,Scalar,Primary,5>(mat,RHS,X0,
3,maxiteIB,sol,
logger,0.00001);
std::string locationRes("../data/res/MatconeSpherePC_MAIN_MAIN_0.res");
std::ofstream fileRes(locationRes,std::ios::out);
fileRes<<"Ite\tSizeBlck\tMinRes\tMaxRes\tTime\tMatVect\tOrtho\tResidual\tSVD\n";
logger.forEachIte([&](int ite, int sizeB, Primary min, Primary max, double* time)
{
fileRes<<ite<<"\t"
<<sizeB<<"\t"
<<min<<"\t"
<<max<<"\t"
<<*time<<"\t";
for(int i=1 ; i<5 ; ++i){
fileRes<<time[i]<<"\t";
}
fileRes<<"\n";
});
std::cout<<"Return value / max_ite : "<<nb<<"/"<<maxiteIB<<"\n";
//Check solution
{
//Compute what i need for my metrics
//MatxSol : Mat*Sol_computed
Block<Scalar,Primary> MatxSol(SizeBlock,dim);
mat.MatBlockVect(sol,MatxSol);
//MatxSolMinusB : Mat*Sol_computed-B
Block<Scalar,Primary> MatxSolMinusB(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];
}
}
//Then display for each vect
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
std::cout<<"Metrics on "<<idxRHS<<"\n";
Primary SecondMemberError = MatxSolMinusB.getNorm(idxRHS);
std::cout <<"["<<idxRHS<<"]\t"
<<"Second Member Error \t" << SecondMemberError <<"\n";
Primary normB = RHS.getNorm(idxRHS);
Primary NormwiseBackwardError = SecondMemberError / normB;
std::cout <<"["<<idxRHS<<"]\t"
<<"Normwise Backward Error\t" << NormwiseBackwardError <<"\n\n";
}
}
return 0;
}
......@@ -3,10 +3,11 @@
#include <fstream>
#include "../src/Block.hpp"
#include "../src/BlockGMRes.hpp"
#include "../src/BlockGMResQR.hpp"
#include "../src/LapackInterface.hpp"
#include "../src/Logger.hpp"
#include "../src/BlockArnoldiIb.hpp"
#include "../src/IBGMRes.hpp"
/**
* Matrix Class : need to implement MatBlockVect(), MatBaseproduct(),
......@@ -49,11 +50,12 @@ public:
}
}
void MatBlockVect(Block<Scalar,Primary>& in, Block<Scalar,Primary>& out){
void MatBlockVect(Block<Scalar,Primary>& in, Block<Scalar,Primary>& out,
int idxToWrite=0){
call_gemm<Scalar>(dim,in.getSizeBlock(),dim,
data,dim,
in.getPtr(), in.getLeadingDim(),
out.getPtr(),out.getLeadingDim(),
out.getPtr(idxToWrite),out.getLeadingDim(),
Scalar(1),Scalar(0));
}
......@@ -66,7 +68,7 @@ public:
}
bool useRightPreCond(){
return false;
return true;
}
void preCondBlockVect(Block<Scalar,Primary>& input, Block<Scalar,Primary>& output){
for(int j=0 ; j<output.getSizeBlock() ; ++j){
......@@ -126,6 +128,7 @@ public:
};
int main(int ac, char ** av){
std::vector<std::string > args(av,av+ac);
for(auto ag : args)
......@@ -146,7 +149,7 @@ int main(int ac, char ** av){
mat.importMatrix(myfile);
int dim = mat.size();
int SizeBlock = 29;
int SizeBlock = 15;
int maxiteIB = dim/SizeBlock + (SizeBlock-1);
......@@ -157,19 +160,45 @@ int main(int ac, char ** av){
//Fill it!
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ;j<dim; ++j){
XSol.getPtr(i)[j] = {dis(gen),/*dis(gen)*/0};
X0.getPtr(i)[j] = {dis(gen),/*dis(gen)*/0};
XSol.getPtr(i)[j] = {dis(gen),dis(gen)};
X0.getPtr(i)[j] = {dis(gen),dis(gen)};
}
}
mat.MatBlockVect(XSol,RHS);
Block<Scalar,Primary> sol(SizeBlock,dim);
Logger<Primary,5> logger;
int nbIte = dim/SizeBlock;
int nbIteIB = dim/SizeBlock+SizeBlock-1;
int restart = nbIte;
int nb = BlockGMResQR<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,
X0,nbIte,sol,
restart,logger,
0.000001);
// int nb = IBGMRes<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,X0,nbIteIB,
// restart,sol,logger,
// 0.0001);
// int nb = BlockArnoldiIb<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,
// X0,
// maxiteIB,
// sol,logger,
// 0.0001);
std::string locationRes("../data/res/young1c.res");
std::ofstream fileRes(locationRes,std::ios::out);
logger.forEachIte([&](int ite, int sizeB, Primary min, Primary max, double* time)
{
fileRes<<ite<<"\t"
<<sizeB<<"\t"
<<min<<"\t"
<<max<<"\t"
<<*time<<"\n";
});
int nb = BlockArnoldiIb<InputMatrix<Scalar,Primary>,Scalar,Primary>(mat,RHS,
X0,
maxiteIB,
sol);
std::cout<<"Return value / max_ite : "<<nb<<"/"<<maxiteIB<<"\n";
//Check solution
......
#include <vector>
#include <complex>
#include <random>
#include "../src/LapackInterface.hpp"