Commit 6d089831 authored by Cyrille Piacibello's avatar Cyrille Piacibello

Changes from Block<T> to Block<T,S=T> in order to deal with complex

parent 46985994
......@@ -15,17 +15,23 @@ public:
virtual ~EngineDispatcher(){};
virtual void set_rightprecond(Callback_RightCond rightPreCond,void * ptr){};
virtual void set_rightprecond(Callback_RightCond /*rightPreCond*/,void * /*ptr*/){};
virtual void set_usergemm(Callback_MatBlockVect userProduct){};
virtual void set_usergemm(Callback_MatBlockVect /*userProduct*/){};
virtual void set_parameters(int inMaxIte,int restart, void * inTolerance){};
virtual void set_parameters(int /*inMaxIte*/,int /*restart*/, void * /*inTolerance*/){};
virtual int solve(int, void * , void *){};
virtual int solve(int, void * , void *){
return 0;
};
virtual void * get_results(){};
virtual void * get_results(){
return nullptr;
};
virtual double * get_logs(){};
virtual double * get_logs(){
return nullptr;
};
};
......@@ -34,7 +40,7 @@ template<typename T,typename S=T>
class IbGMResDrEngine : public EngineDispatcher{
//User parameters
UserMatrix<T>* matrix;
UserMatrix<T,S>* matrix;
void * userEnvPtr;
//Solve parameters
......@@ -43,7 +49,7 @@ class IbGMResDrEngine : public EngineDispatcher{
int restart;
//Solution storage
Block<T>* sol;
Block<T,S>* sol;
Logger<S>* log;
public:
IbGMResDrEngine(void * userMatrix, int dim, void * inUserEnv) :
......@@ -55,7 +61,7 @@ public:
sol(nullptr),
log(nullptr)
{
matrix = new UserMatrix<T>(userEnvPtr,userMatrix,dim);
matrix = new UserMatrix<T,S>(userEnvPtr,userMatrix,dim);
}
~IbGMResDrEngine(){
......@@ -88,15 +94,15 @@ public:
//Solve method : RHS and X0 have the same size
int solve(int nbRHS, void * RHS , void * X0)override{
//Convert raw data to Block
Block<T> X_init(nbRHS,matrix->size());
Block<T> B(nbRHS,matrix->size());
Block<T,S> X_init(nbRHS,matrix->size());
Block<T,S> B(nbRHS,matrix->size());
X_init.InitBlock(reinterpret_cast<T*>(X0));
B.InitBlock(reinterpret_cast<T*>(RHS));
//Init an empty block to store solution
sol = new Block<T>(nbRHS,matrix->size());
sol = new Block<T,S>(nbRHS,matrix->size());
log = new Logger<S>();
//Call to BlockGMRes with parameters
int res = BlockGMRes<UserMatrix<T>,T,S>(*matrix,B,X_init,max_iter,*sol,
int res = BlockGMRes<UserMatrix<T,S>,T,S>(*matrix,B,X_init,max_iter,*sol,
restart,*log,tolerance);
return res;
}
......@@ -140,7 +146,7 @@ extern "C" void ibgmresdr_set_parameters(int max_ite, int restart, void * tolera
}
extern "C" int ibgmresdr_solve(int nbRHS,void * RHS,void *X0,
IbGMResDr_handle handle){
(reinterpret_cast<EngineDispatcher*>(handle))->solve(nbRHS,RHS,X0);
return (reinterpret_cast<EngineDispatcher*>(handle))->solve(nbRHS,RHS,X0);
}
extern "C" void ibgmresdr_dealloc(IbGMResDr_handle handle){
......
......@@ -9,7 +9,7 @@
*
* Note : the matrix is supposed to be square
*/
template<typename T>
template<typename T,typename S=T>
class UserMatrix{
private:
typedef void (*Callback_MatBlockVect)(void * ,int , void * ,void **,void*);
......@@ -49,13 +49,13 @@ public:
~UserMatrix(){
}
void MatBlockVect(Block<T>& input, Block<T>& output){
void MatBlockVect(Block<T,S>& input, Block<T,S>& output){
void* toWrite = output.getPtr();
userProduct(userPtrMat,input.getSizeBlock(),input.getPtr(),&toWrite,
userEnv);
}
void MatBaseProduct(T * ptrToRead,Block<T>& output){
void MatBaseProduct(T * ptrToRead,Block<T,S>& output){
void* toWrite = output.getPtr();
userProduct(userPtrMat,output.getSizeBlock(),ptrToRead,&toWrite,
userEnv);
......@@ -76,13 +76,13 @@ public:
return (userPtrPreCond!=nullptr)&&(rightPreCond!=nullptr);
}
void preCondBlockVect(Block<T>& input, Block<T>& output){
void preCondBlockVect(Block<T,S>& input, Block<T,S>& output){
void* toWrite = output.getPtr();
rightPreCond(userPtrPreCond,input.getSizeBlock(),input.getPtr(),&toWrite,
userEnv);
}
void preCondBaseProduct(T * ptrToRead,Block<T>& output){
void preCondBaseProduct(T * ptrToRead,Block<T,S>& output){
void* toWrite = output.getPtr();
rightPreCond(userPtrPreCond,output.getSizeBlock(),ptrToRead,&toWrite,
userEnv);
......
......@@ -52,7 +52,7 @@ option(DEBUG_MODE "Set to On to compile with debug info" OFF)
option(BUILD_SHARED_LIBS "Build shared libraries" OFF)
if(DEBUG_MODE)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g3")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g3 -W -Wall -fdiagnostics-color=always")
endif()
if(BUILD_SHARED_LIBS)
......
......@@ -5,7 +5,7 @@
#include "LapackInterface.hpp"
template<class Scalar>
template<class Scalar,typename Primary=Scalar>
class Base{
private:
int SizeBlock;
......@@ -72,7 +72,7 @@ public:
return currentBlockUsed+1;
}
void addBlockDatas(Block<Scalar>& inBlock){
void addBlockDatas(Block<Scalar,Primary>& inBlock){
//Check
if(inBlock.getLeadingDim() != ldBase){std::cout<<"Problem in addBlockDatas\n";}
......@@ -90,7 +90,7 @@ public:
* orthogonalization coeffs
* @param ldHess : Leading dimension of Hessenberg
*/
void ComputeMGSOrtho(Block<Scalar> & Wj, Scalar * HessToWrite, int ldHess){
void ComputeMGSOrtho(Block<Scalar,Primary> & Wj, Scalar * HessToWrite, int ldHess){
Scalar * ptrOfHess = HessToWrite;
//Loop over different blocks
......@@ -117,7 +117,7 @@ public:
* @brief Compute the product Base*input -> output
*
*/
void ComputeProduct(Block<Scalar> & input, Block<Scalar>& output){
void ComputeProduct(Block<Scalar,Primary> & input, Block<Scalar,Primary>& output){
int minLine = std::min(nbVectInBlock[currentBlockUsed],input.getLeadingDim());
call_gemm<Scalar>(ldBase,input.getSizeBlock(),minLine,
getPtr(),ldBase,
......
......@@ -2,13 +2,14 @@
#define BLOCK_HPP
#include <string.h>
#include "LapackInterface.hpp"
#include "numerical.hpp"
#include "LapackInterface.hpp"
/**
* @brief Will handle Wj
*
*/
template<class Scalar>
template<class Scalar,typename Primary=Scalar>
class Block{
private:
int SizeBlock;
......@@ -92,7 +93,6 @@ public:
* Return the norm of id vector
*
*/
template<typename Primary>
Primary getNorm(int id){
//Starting point
Scalar * vect = getPtr(id);
......@@ -107,7 +107,7 @@ public:
/**
* @brief Copy input Block datas inside
*/
void CopyBlock(Block<Scalar>& in){
void CopyBlock(Block<Scalar,Primary>& in){
//Check
if(in.getSizeBlock() != SizeBlock){
std::cout<<"SizeBlock different between input block and current block\n";
......@@ -129,6 +129,18 @@ public:
memcpy(getPtr(),rawDatas, sizeof(Scalar)*ldb*SizeBlock);
}
/**
* @brief Compute SVD of current block, and write down U inside
* output, and Singular Values inside sigma.
*/
void DecompositionSVD(Block<Scalar,Primary>& output, Block<Scalar,Primary> sigma){
//What I need:
//jobu : A --> all the U part is computed
//jobvt : N --> no part of V^{H} is computed
}
//Methods for Pj
/**
......@@ -137,6 +149,23 @@ public:
int getUsedVects()const{
return nbVectUsed;
}
void OrthoBlock(Block<Scalar,Primary>& input,Scalar * toWrite, int ldToWrite){
Fill the coeff
call_Tgemm<>(getUsedVects(),input.getSizeBlock(),getLeadingDim(),
getPtr(),getLeadingDim(),
input.getPtr(),input.getLeadingDim(),
toWrite,ldToWrite,Scalar(1),Scalar(0));
//Read again the coeffs in order to modify input
//Input = (-1) * Pj * toWrite + Wj;
call_gemm<Scalar>(getLeadingDim(),input.getSizeBlock(),getUsedVects(),
getPtr(i),getLeadingDim(),
toWrite, ldToWrite,
Wj.getPtr(),Wj.getLeadingDim(),
Scalar(-1), Scalar(1));
}
};
......
......@@ -7,8 +7,8 @@
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<class Matrix,typename Scalar,typename Primary=Scalar>
int BlockArnoldiIb(Matrix& A,Block<Scalar>& B, Block<Scalar>& X0,
int max_iter, Block<Scalar>& Sol,
int BlockArnoldiIb(Matrix& A,Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
int max_iter, Block<Scalar,Primary>& Sol,
Primary epsilon = 0.00001){
......@@ -31,12 +31,12 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar>& B, Block<Scalar>& X0,
Base<Scalar> base(SizeBlock,max_iter+1,dim);
//Pj is a matrix ortho built with the abandonned directions, thus
//it's cannot be larger than the range(B).
Block<Scalar> Pj(SizeBlock,dim);
Block<Scalar,Primary> Pj(SizeBlock,dim);
Hessenberg<Scalar> L(nb_iter_restart,SizeBlock);
//Compute first block residual (step 1)
Block<Scalar> R0(SizeBlock,dim);
Block<Scalar,Primary> R0(SizeBlock,dim);
//Fill matrix R0 with each vector r0_i = B[i] - A * X0[i]
A.MatBlockVect(X0,R0); //R0 = A*X0
for(int i=0 ; i<SizeBlock ; ++i){ //R0 = B - R0
......@@ -46,8 +46,8 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar>& B, Block<Scalar>& X0,
}
//V1 will be first block of base, and R1, will be used for least
//square pb, in criteria verification
Block<Scalar> V1(SizeBlock,dim);
Block<Scalar> R1(SizeBlock,SizeBlock);
Block<Scalar,Primary> V1(SizeBlock,dim);
Block<Scalar,Primary> R1(SizeBlock,SizeBlock);
//Q-R Decomposition of (R0)
R0.ComputeQRFacto(V1,R1);
//Fisrt block of base added
......@@ -61,7 +61,7 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar>& B, Block<Scalar>& X0,
//Main Loop
for(int j=1 ; j<max_iter ; ++j){
//Compute block Wj
Block<Scalar> Wj(SizeBlock,dim);
Block<Scalar,Primary> Wj(SizeBlock,dim);
A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),Wj);
//Ortho against V and store coeff inside new block column of Hess L
Scalar * toWrite = L.getPtr(j);
......@@ -69,7 +69,7 @@ int BlockArnoldiIb(Matrix& A,Block<Scalar>& B, Block<Scalar>& X0,
std::cout<<"Ortho against V done"<<std::endl;
//Same against Pj
Pj.OrthoBlock(Wj,);
}
}
......
......@@ -15,9 +15,9 @@
* Ref : IB-BGMRES-DR (Giraud Jing Agullo): p.21
*/
template<class Matrix,typename Scalar,typename Primary=Scalar>
int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
int BlockGMRes(Matrix& A, Block<Scalar,Primary>& B, Block<Scalar,Primary>& X0,
//Matrix preCond Deal with that later
int max_iter, Block<Scalar>& Sol, int nb_iter_restart,
int max_iter, Block<Scalar,Primary>& Sol, int nb_iter_restart,
Logger<Primary>& log,
Primary epsilon = 0.00001){
......@@ -34,18 +34,18 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
//store norms
std::vector<Primary> normB;
for(int i=0 ; i< SizeBlock ; ++i){
normB.push_back(B.template getNorm<Primary>(i));
normB.push_back(B.getNorm(i));
}
//Residual vectors needed for computing the solution
Block<Scalar> Yj(SizeBlock,dim);
Block<Scalar,Primary> Yj(SizeBlock,dim);
//Add a loop for restart
int k=0,nbRestart=0, globalStep=0;
log.init_log();
while(k<max_iter){
std::cout<<"\t\t\tRestarting number "<<nbRestart<<" ("<<nb_iter_restart<<")\n";
Block<Scalar> R0(SizeBlock,dim);
Block<Scalar,Primary> R0(SizeBlock,dim);
//Fill matrix R0 with each vector r0_i = B[i] - A * X0[i]
A.MatBlockVect(X0,R0); //R0 = A*X0
for(int i=0 ; i<SizeBlock ; ++i){ //R0 = B - R0
......@@ -55,12 +55,12 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
}
//Storage needed for computating
Base<Scalar> base(SizeBlock,nb_iter_restart+1,dim);
Base<Scalar,Primary> base(SizeBlock,nb_iter_restart+1,dim);
Hessenberg<Scalar> Hess(nb_iter_restart,SizeBlock);
//V1 will be first block of base, and R1, will be used for least
//square pb, in criteria verification
Block<Scalar> V1(SizeBlock,dim);
Block<Scalar> R1(SizeBlock,SizeBlock);
Block<Scalar,Primary> V1(SizeBlock,dim);
Block<Scalar,Primary> R1(SizeBlock,SizeBlock);
//Q-R Decomposition of (R0)
R0.ComputeQRFacto(V1,R1);
//Fisrt block of base added
......@@ -69,10 +69,10 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
std::cout<<"Starting main loop"<<std::endl;
for(int j=0 ; (j<nb_iter_restart) && (j+k)<= max_iter ; ++j){
//Compute block Wj
Block<Scalar> Wj(SizeBlock,dim);
Block<Scalar,Primary> Wj(SizeBlock,dim);
if(A.useRightPreCond()){//Wj = A*preCond*V[j] if preCond
Block<Scalar> temp(SizeBlock,dim);
Block<Scalar,Primary> temp(SizeBlock,dim);
A.preCondBaseProduct(base.getPtr(base.getCurrentBlock()-1),temp);
A.MatBlockVect(temp,Wj);
}else{//Wj = A*V[j] without preCond
......@@ -86,8 +86,8 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
std::cout<<"Ortho done"<<std::endl;
//QR facto of orthogonalized Wj
Block<Scalar> Q(SizeBlock,dim);
Block<Scalar> R(SizeBlock,SizeBlock);
Block<Scalar,Primary> Q(SizeBlock,dim);
Block<Scalar,Primary> R(SizeBlock,SizeBlock);
Wj.ComputeQRFacto(Q,R);
//Copy Block R inside Hess
Hess.addRPartToHess(R);
......@@ -110,9 +110,9 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
//Generate block Sol = Base*Yj
//if Precond : block Sol = preCond*Base*yj
Block<Scalar> blockSol(SizeBlock,dim);
Block<Scalar,Primary> blockSol(SizeBlock,dim);
if(A.useRightPreCond()){
Block<Scalar> temp(SizeBlock,dim);
Block<Scalar,Primary> temp(SizeBlock,dim);
base.ComputeProduct(Yj,temp);
A.preCondBlockVect(temp,blockSol);
}else{
......@@ -126,7 +126,7 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
}
}
//Compute A*sol computed
Block<Scalar> AxSol(SizeBlock,dim);
Block<Scalar,Primary> AxSol(SizeBlock,dim);
A.MatBlockVect(blockSol,AxSol);
//Compute real residual
......@@ -137,7 +137,7 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
}
//For each vector, compute norm
for(int i = 0 ; i<SizeBlock ; ++i){
auto norm = AxSol.template getNorm<Primary>(i)/normB[i];
auto norm = AxSol.getNorm(i)/normB[i];
if(norm > RealMinMax.second){
RealMinMax.second = norm;
}
......@@ -173,9 +173,9 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
// X0 = X0 + Base*Yj
//else
// X0 = X0 + preCond*Base*yj
Block<Scalar> newX0(SizeBlock,dim);
Block<Scalar,Primary> newX0(SizeBlock,dim);
if(A.useRightPreCond()){
Block<Scalar> temp(SizeBlock,dim);
Block<Scalar,Primary> temp(SizeBlock,dim);
base.ComputeProduct(Yj,temp);
A.preCondBlockVect(temp,newX0);
}else{
......
......@@ -78,7 +78,8 @@ public:
}
//Here add the last triang block
void addRPartToHess(Block<Scalar> & R){
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);
......@@ -102,7 +103,8 @@ public:
}
//Call to lapack here
void ComputeBlockResidual(Block<Scalar>& R1, Block<Scalar> & Yj){
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];
......@@ -143,10 +145,10 @@ public:
//Compute the Hess*input set of vectors, and return min and max norm
template<typename Primary>
std::pair<Primary,Primary> displayBlockResidual(Block<Scalar>& Y,Block<Scalar>& R,
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> blockRes(SizeBlock,(currentSize+1)*SizeBlock);
Block<Scalar,Primary> blockRes(SizeBlock,(currentSize+1)*SizeBlock);
call_gemm<Scalar>((currentSize+1)*SizeBlock,SizeBlock,SizeBlock*currentSize,
datas, ldh,
Y.getPtr(), Y.getLeadingDim(),
......@@ -163,7 +165,7 @@ public:
//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.template getNorm<Primary>(idxRHS)/normB[idxRHS];
auto norm = blockRes.getNorm(idxRHS)/normB[idxRHS];
//std::cout<<"Norm of vector["<<idxRHS<<"] is "<<norm<<"\n";
if(norm > MinMax.second){
MinMax.second = norm;
......
......@@ -311,10 +311,18 @@ int callLapack_orgqr(int m, int p,
* char jobvt, lapack_int m, lapack_int n, lapack_complex_float* a,
* lapack_int lda, float* s, lapack_complex_float* u, lapack_int ldu,
* lapack_complex_float* vt, lapack_int ldvt, float* superb );
*
* Ref : https://software.intel.com/en-us/node/521150
*/
template<typename Scalar>
int callLapack_gesvd(){
int callLapack_gesvd(int ,int ,
Scalar *, int,
Scalar *, int){
std::cout<<"Only complx float and compl double supported\n";
exit(0);
return 0;
}
#endif //LapackInterface.hpp
......@@ -55,7 +55,7 @@ public:
double * writeDownArray(){
array = new double[3 * iterations.size()];
for(int i =0 ; i<iterations.size() ; ++i){
for(unsigned int i =0 ; i<iterations.size() ; ++i){
array[3*i] = iterations[i].min;
array[3*i+1] = iterations[i].max;
array[3*i+2] = iterations[i].timestamp;
......
#ifndef NUMERICAL_HPP
#define NUMERICAL_HPP
#include <vector>
#include <complex>
#include <cmath>
#include <vector>
#include <array>
#include <algorithm>
//template function to get the module of a cmplx (or the square of
......
#include "../src/Block.hpp"
#include <vector>
#include <random>
int main(int ac, char ** av){
std::vector<std::string > args(av,av+ac);
for(auto ag : args)
std::cout << ag << "\n";
using Scalar=double;
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<> dis(0, 2);
......
#include <random>
#include "../src/Block.hpp"
#include "../src/BlockGMRes.hpp"
#include "../src/LapackInterface.hpp"
#include "../src/Logger.hpp"
/**
* Matrix Class : need to implement MatBlockVect(), MatBaseproduct(),
* useRightPreCond(), preCondBlockVect(),preCondBaseProduct() and
......@@ -142,7 +144,6 @@ int main(int ac, char ** av){
//Check solution
{
//Compute what i need for my metrics
//MatxSol : Mat*Sol_computed
Block<Scalar> MatxSol(SizeBlock,dim);
mat.MatBlockVect(sol,MatxSol);
......@@ -161,16 +162,16 @@ int main(int ac, char ** av){
//Then display for each vect
for(int idxRHS=0 ; idxRHS<SizeBlock ; ++idxRHS){
std::cout<<"Metrics on "<<idxRHS<<"\n";
Primary ForwardError = X_Diff.getNorm<Primary>(idxRHS);
Primary ForwardError = X_Diff.getNorm(idxRHS);
std::cout <<"["<<idxRHS<<"]\t"
<<"Forward Error (on X) \t" << ForwardError <<"\n";
Primary SecondMemberError = MatxSolMinusB.getNorm<Primary>(idxRHS);
Primary SecondMemberError = MatxSolMinusB.getNorm(idxRHS);
std::cout <<"["<<idxRHS<<"]\t"
<<"Second Member Error \t" << SecondMemberError <<"\n";
Primary ForwardErrorNormalized = ForwardError/(X_Exact.getNorm<Primary>(idxRHS));
Primary ForwardErrorNormalized = ForwardError/(X_Exact.getNorm(idxRHS));
std::cout <<"["<<idxRHS<<"]\t"
<<"Forward Error Norm. \t" << ForwardErrorNormalized <<"\n";
Primary normB = RHS.getNorm<Primary>(idxRHS);
Primary normB = RHS.getNorm(idxRHS);
Primary NormwiseBackwardError = SecondMemberError / normB;
std::cout <<"["<<idxRHS<<"]\t"
<<"Normwise Backward Error\t" << NormwiseBackwardError <<"\n\n";
......
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