Commit 5d73507d authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

Forget Utils.hpp

parent 49e17ea6
#ifndef UTILS_HPP
#define UTILS_HPP
#include "Base.hpp"
#include "Block.hpp"
/**
* @brief This struct will hold some output concerning an arnoldi
* procedure.
*/
struct ArnReturn{
int nbIteDone;
bool hasConverged;
};
/**
* @brief This function test if the current solution is accurate
* enough. It will works only if the RHS are scaled, (i.e each column
* vector of B has a norm of 1)
*
* @param A Matrix provided buy the user
* @param B RHS (each column vector of B must be of norm 1)
* @param X0 : initial guess, will be overwritten by
* solution if criteria is satisfied
* @param Base Current base
* @param Y Current solution of least square problem
* @param epsilon target accuracy
*/
template<class Matrix, class Scalar, class Primary=Scalar>
bool CheckCurrentSolution(Matrix&A, Block<Scalar,Primary>& B,
Block<Scalar,Primary>& X0,
Base<Scalar,Primary>& Base,
Block<Scalar,Primary>& Y, Primary epsilon){
int dim = A.size();
int SizeBlock = B.getSizeBlock();
//Create temporary block to store solution
Block<Scalar,Primary> Sol(SizeBlock,dim);
if(A.useRightPreCond()){
Block<Scalar,Primary> temp(SizeBlock,dim);
Base.ComputeProduct(Y,temp);
A.preCondBlockVect(temp,Sol);
}else{//No pre Cond used
Base.ComputeProduct(Y,Sol);
}
//Adding X0
for(int i=0 ; i<SizeBlock ; ++i){
for(int k=0 ; k<dim ; ++k){
Sol.getPtr(i)[k] += X0.getPtr(i)[k];
}
}
//Then ||A*Sol -B||
Block<Scalar,Primary> AxSol(SizeBlock,dim);
A.MatBlockVect(Sol,AxSol);
for(int i=0 ; i<SizeBlock ; ++i){
for(int j=0 ; j<dim ; ++j){
AxSol.getPtr(i)[j] -= B.getPtr(i)[j];
}
}
std::pair<Primary,Primary> RealMinMax = AxSol.getMinMaxNorm();
if(RealMinMax.second < epsilon){
std::cout<<"Convergence !!"<<std::endl;
X0.CopyBlock(Sol);
return true;
}
return false;
}
/**
* @brief This function compare two blocks of vectors of the same
* size, and display the min and max of frobenius norm overs the vectors
*
* @param block1
* @param block2
*/
template<class Scalar, class Primary>
std::pair<Primary,Primary>
CompareBlocks(Block<Scalar,Primary>& block1, Block<Scalar,Primary>& block2){
std::pair<Primary,Primary> MinMax{1000,0};
if(block1.getSizeBlock() != block2.getSizeBlock()){
std::cout<<"Size of blocks not compatible\nExiting\n";
exit(0);
}
if(block1.getLeadingDim() != block2.getLeadingDim()){
std::cout<<"Leading dim of blocks not compatible\nExiting\n";
exit(0);
}
for(int i=0 ; i<block1.getSizeBlock() ; ++i){
Primary norm{0};
for(int j=0 ; j<block1.getLeadingDim() ; ++j){
norm += module<Scalar,Primary>(block1.getPtr(i)[j] - block2.getPtr(i)[j]);
}
auto norm2 = std::sqrt(norm);
if(norm2 > MinMax.second){
MinMax.second = norm2;
}
if(norm2 < MinMax.first){
MinMax.first = norm2;
}
}
return MinMax;
}
#endif //UTILS_HPP
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