Right PreCond is available and working

Api/src/Ib-GMRes-Dr.h

@@ -101,11 +101,12 @@ void ibgmresdr_set_usergemm(Callback_MatBlockVect userProduct,
  *
  * @param userPreCond : callback, refer to documentation to get the
  * prototype
-
+ * @param ptrRightPreCond : ptr to what user need to compute right
+ * preconditionning. Won't be read by library.
  * @Attention : If this function is never called, then the library
  * will assume no PreConditionner is used.
  */
-void ibgmresdr_set_rightprecond(Callback_RightCond rightPreCond,
+void ibgmresdr_set_rightprecond(Callback_RightCond rightPreCond,void * ptrRightPreCond,
                                 IbGMResDr_handle handle);
 
 

Api/src/IbGMResDrEngine.hpp

@@ -15,7 +15,7 @@ public:
 
     virtual ~EngineDispatcher(){};
 
-    virtual void set_rightprecond(Callback_RightCond rightPreCond){};
+    virtual void set_rightprecond(Callback_RightCond rightPreCond,void * ptr){};
 
     virtual void set_usergemm(Callback_MatBlockVect userProduct){};
 
@@ -70,7 +70,9 @@ public:
      * @brief This function allow to set a right pre conditionner
      *
      */
-    void set_rightprecond(Callback_RightCond rightPreCond)override{
+    void set_rightprecond(Callback_RightCond rightPreCond,
+                          void * usrPtrRightPreCond)override{
+        matrix->setRightPreCond(rightPreCond,usrPtrRightPreCond);
     }
 
     void set_usergemm(Callback_MatBlockVect userProduct)override{
@@ -125,8 +127,10 @@ extern "C" void ibgmresdr_set_usergemm(Callback_MatBlockVect userProduct,
 }
 
 extern "C" void ibgmresdr_set_rightprecond(Callback_RightCond rightPreCond,
+                                           void * ptrRightPreCond,
                                            IbGMResDr_handle handle){
-    (reinterpret_cast<EngineDispatcher*>(handle))->set_rightprecond(rightPreCond);
+    (reinterpret_cast<EngineDispatcher*>(handle))->set_rightprecond(rightPreCond,
+                                                                    ptrRightPreCond);
 }
 
 extern "C" void ibgmresdr_set_parameters(int max_ite, int restart, void * tolerance,

Api/src/UserMatrix.hpp

@@ -13,15 +13,23 @@ template<typename T>
 class UserMatrix{
 private:
     typedef void (*Callback_MatBlockVect)(void * ,int , void * ,void **,void*);
+    typedef void (*Callback_RightCond)(void * ,int , void * ,void ** , void *);
+
     Callback_MatBlockVect userProduct;
+    Callback_RightCond rightPreCond;
+
     int dim;
     void* userPtrMat;
-    void * userEnv;
+    void* userPtrPreCond;
+    void* userEnv;
+
 public:
 
     UserMatrix() : userProduct(nullptr),
+                   rightPreCond(nullptr),
                    dim(0), //Do nothing
                    userPtrMat(nullptr),
+                   userPtrPreCond(nullptr),
                    userEnv(nullptr){
     }
 
@@ -34,6 +42,7 @@ public:
         userProduct(nullptr),
         dim(inDim),
         userPtrMat(userPtrMatrix),
+        userPtrPreCond(nullptr),
         userEnv(inUserEnv){
     }
 
@@ -55,6 +64,29 @@ public:
     void setMatBlockVect(Callback_MatBlockVect inUserProduct){
         userProduct = inUserProduct;
     }
+
+    //Right Pre Conditionner
+    //Since it is set by the user, we use the same structure to hold it
+    void setRightPreCond(Callback_RightCond inRightPreCond, void * usrPtrRightpreCond){
+        rightPreCond = inRightPreCond;
+        userPtrPreCond = usrPtrRightpreCond;
+    }
+
+    bool useRightPreCond(){
+        return (userPtrPreCond!=nullptr)&&(rightPreCond!=nullptr);
+    }
+
+    void preCondBlockVect(Block<T>& input, Block<T>& output){
+        void* toWrite = output.getPtr();
+        rightPreCond(userPtrPreCond,input.getSizeBlock(),input.getPtr(),&toWrite,
+                    userEnv);
+    }
+
+    void preCondBaseProduct(T * ptrToRead,Block<T>& output){
+        void* toWrite = output.getPtr();
+        rightPreCond(userPtrPreCond,output.getSizeBlock(),ptrToRead,&toWrite,
+                     userEnv);
+    }
 };
 
 

Api/tests/testApiGMres.c

 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
+#include <time.h>
 #include "Ib-GMRes-Dr.h"
 
 

Api/tests/testApiGMresUserCase.c

@@ -24,6 +24,13 @@ typedef struct matrix{
 }Matrix;
 
 
+//Struct for storing the Right PreCond
+typedef struct rightPreCond{
+    int dim;
+    double complex * datas;
+}RightPreCond;
+
+
 //Define callbacks
 //Matrix block of vector products
 void MatxBlockVect(void * mat, int nbVect, void * input, void **output, void * env){
@@ -40,6 +47,22 @@ void MatxBlockVect(void * mat, int nbVect, void * input, void **output, void * e
                 input,nbCol,&beta,*output,nbCol);
 }
 
+void RightPreCondxBlockVect(void * rpc, int nbVect, void * input, void **output,
+                            void * env){
+    RightPreCond * my_rpc = (RightPreCond *) rpc;
+    Env * my_env = (Env*)env;
+
+    int nbRHS = my_env->nbRHS;
+    int size = my_rpc->dim;
+    double complex * inputBlock = input;
+    double complex * outputBlock = *output;
+    for(int i=0 ; i<nbRHS ; ++i){
+        for(int j=0 ; j<size ; ++j){
+            outputBlock[i*size+j] = inputBlock[i*size+j] * my_rpc->datas[j];
+        }
+    }
+}
+
 void load_datas(Matrix * mat,FILE * fd){
     int ari;
     fread(&ari,sizeof(int),1,fd);
@@ -90,17 +113,29 @@ int main(int ac , char ** av){
     userEnv->nbRHS = RHS->dim2;
 
     //Create Starting Block
-    double * X0 = malloc(sizeof(double)* userEnv->nbRHS * mat->dim1);
-
+    double * X0 = malloc(sizeof(double complex)* userEnv->nbRHS * mat->dim1);
+    memset(X0,0,sizeof(double complex)*userEnv->nbRHS * mat->dim1);
     //Init library
     IbGMResDr_handle handle = Init_ibgmresdr(CMPLX_DBL,mat,mat->dim1,userEnv);
     //Set my product Matrix block of vectors
     ibgmresdr_set_usergemm(MatxBlockVect,handle);
+
+    //Create right precond
+    RightPreCond * my_rpc = malloc(sizeof(RightPreCond));
+    my_rpc->datas = malloc(sizeof(double complex)*mat->dim1);
+    my_rpc->dim = mat->dim1;
+
+    //Fill with the diag of Mat
+    for(int j=0 ; j<mat->dim1 ; ++j){
+        my_rpc->datas[j] = mat->datas[j*(mat->dim1)+j];
+    }
+    //Give to the library the PreConditionner and the method to call
+    ibgmresdr_set_rightprecond(RightPreCondxBlockVect,my_rpc,handle);
     //Tolerance
-    double tol = 0.01;
+    double tol = 0.0001;
     //set param : nb_iter, restart, ptr to tolerance
-    ibgmresdr_set_parameters((mat->dim1)/userEnv->nbRHS,
-                             ((mat->dim1)/userEnv->nbRHS)/2,
+    ibgmresdr_set_parameters((mat->dim1)/userEnv->nbRHS,     //nb max ite
+                             ((mat->dim1)/userEnv->nbRHS)/2, //size of restart
                              &tol, handle);
 
     //Solve env
@@ -125,7 +160,8 @@ int main(int ac , char ** av){
     //Free the handle
     ibgmresdr_dealloc(handle);
 
-
+    free(my_rpc->datas);
+    free(my_rpc);
     free(X0);
     free(userEnv);
     free(RHS->datas);

CMakeLists.txt

@@ -54,6 +54,8 @@ if(DEBUG_MODE)
   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g3")
 endif()
 
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
+
 #Ajout d'un repertoire src
 add_subdirectory(src)
 #Ajout du meme repertoire pour les includes des tests

src/BlockGMRes.hpp

@@ -21,11 +21,11 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
                Logger<Primary>& log,
                Primary epsilon = 0.00001){
 
-    // if(preCond.size() != 0){
-    //     std::cout<<"Right Pre Conditionner used\n";
-    // }else{
-    //     std::cout<<"NO - Right Pre Conditionner used\n";
-    // }
+    if(A.useRightPreCond()){
+        std::cout<<"Right Pre Conditionner used\n";
+    }else{
+        std::cout<<"NO - Right Pre Conditionner used\n";
+    }
 
     //Global variables
     int dim = A.size();
@@ -71,8 +71,13 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
             //Compute block Wj
             Block<Scalar> Wj(SizeBlock,dim);
 
-            //A*V[j]
-            A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),Wj);
+            if(A.useRightPreCond()){//Wj = A*preCond*V[j] if preCond
+                Block<Scalar> temp(SizeBlock,dim);
+                A.preCondBaseProduct(base.getPtr(base.getCurrentBlock()-1),temp);
+                A.MatBlockVect(temp,Wj);
+            }else{//Wj = A*V[j] without preCond
+                A.MatBaseProduct(base.getPtr(base.getCurrentBlock()-1),Wj);
+            }
 
             //Where to write in hess
             Scalar * toWrite = Hess.getPtr(j);
@@ -102,9 +107,17 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
                     //Compare to the real residual
                     std::pair<Primary,Primary> RealMinMax
                         = std::make_pair<Primary,Primary>(100,0);
-                    //Generate block Sol
+
+                    //Generate block Sol = Base*Yj
+                    //if Precond : block Sol = preCond*Base*yj
                     Block<Scalar> blockSol(SizeBlock,dim);
-                    base.ComputeProduct(Yj,blockSol);
+                    if(A.useRightPreCond()){
+                        Block<Scalar> temp(SizeBlock,dim);
+                        base.ComputeProduct(Yj,temp);
+                        A.preCondBlockVect(temp,blockSol);
+                    }else{
+                        base.ComputeProduct(Yj,blockSol);
+                    }
 
                     //Compute addition with X0
                     for(int i=0 ; i<SizeBlock ; ++i){
@@ -112,7 +125,6 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
                             blockSol.getPtr(i)[j] += X0.getPtr(i)[j];
                         }
                     }
-
                     //Compute A*sol computed
                     Block<Scalar> AxSol(SizeBlock,dim);
                     A.MatBlockVect(blockSol,AxSol);
@@ -152,15 +164,23 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
             //Base augmentation
             base.addBlockDatas(Q);
             std::cout<<"Base has been augmented\n";
-
             globalStep++;
         }
         k += nb_iter_restart;
         nbRestart += 1;
         //Compute new starting block X0
-        //X0 = Base*Yj
+        //if no precond
+        //    X0 = X0 + Base*Yj
+        //else
+        //    X0 = X0 + preCond*Base*yj
         Block<Scalar> newX0(SizeBlock,dim);
-        base.ComputeProduct(Yj,newX0);
+        if(A.useRightPreCond()){
+            Block<Scalar> temp(SizeBlock,dim);
+            base.ComputeProduct(Yj,temp);
+            A.preCondBlockVect(temp,newX0);
+        }else{
+            base.ComputeProduct(Yj,newX0);
+        }
         //Copy back
         for(int i=0 ; i<SizeBlock ; ++i){
             for(int j=0 ; j<dim ; ++j){
@@ -169,7 +189,6 @@ int BlockGMRes(Matrix& A, Block<Scalar>& B, Block<Scalar>& X0,
         }
     }
     Sol.CopyBlock(X0);
-
     return globalStep;
 }