Change dot function to operator() overload, so a lambda function can be passed

Author: RaulPPelaez

README.md

@@ -1,22 +1,22 @@
 # Lanczos solver,
-  Computes the matrix-vector product sqrt(M)·v using a recursive algorithm.  
-  For that, it requires a functor with a "dot" function that takes an output real* array and an input real* (both in device memory if compiled in CUDA mode or host memory otherwise) as:  
-  ```c++ 
-	virtual void dot(real* in_v, real * out_Mv) override;
-  ```  
-  The functor must inherit ```lanczos::MatrixDot``` (see example.cu).
-  This function must fill "out" with the result of performing the M·v dot product- > out = M·a_v.  
-  If M has size NxN and the cost of the dot product is O(M). The total cost of the algorithm is O(m·M). Where m << N.  
-  If M·v performs a dense M-V product, the cost of the algorithm would be O(m·N^2).  
+  Computes the matrix-vector product sqrt(M)·v using a recursive algorithm.
+  For that, it requires a functor that takes an output real* array and an input real* (both in device memory if compiled in CUDA mode or host memory otherwise) as:
+  ```c++
+	void operator()(real* in_v, real * out_Mv) override;
+  ```
+  The functor can inherit ```lanczos::MatrixDot``` (see example.cu).
+  This function must fill "out" with the result of performing the M·v dot product- > out = M·a_v.
+  If M has size NxN and the cost of the dot product is O(M). The total cost of the algorithm is O(m·M). Where m << N.
+  If M·v performs a dense M-V product, the cost of the algorithm would be O(m·N^2).
 
 This is a header-only library, although a shared library can be compiled instead.
 
-## Usage:  
+## Usage:
 
-See example.cu for an usage example that can be compiled to work in GPU or CPU mode instinctively.  
-See example.cpp for a CPU only example.  
+See example.cu for an usage example that can be compiled to work in GPU or CPU mode instinctively.
+See example.cpp for a CPU only example.
 
-Let us go through the remaining one, a GPU-only example.  
+Let us go through the remaining one, a GPU-only example.
 
 Create the module:
 ```c++
@@ -30,8 +30,8 @@ Write a functor that computes the product between the original matrix and a give
 struct DiagonalMatrix: public lanczos::MatrixDot{
   int size;
   DiagonalMatrix(int size): size(size){}
-  
-  void dot(real* v, real* Mv) override{
+
+  void operator()(real* v, real* Mv) override{
     //An example diagonal matrix
     for(int i=0; i<size; i++){
       Mv[i] = 2*v[i];
@@ -42,7 +42,7 @@ struct DiagonalMatrix: public lanczos::MatrixDot{
 
 ```
 
-Provide the solver with an instance of the functor and the target vector:  
+Provide the solver with an instance of the functor and the target vector:
 
 ```c++
     int size = 10;
@@ -55,7 +55,7 @@ Provide the solver with an instance of the functor and the target vector:
 	DiagonalMatrix dot(size);
     //Call the solver
     real* d_result = thrust::raw_pointer_cast(result.data());
-    real* d_v = thrust::raw_pointer_cast(v.data()); 
+    real* d_v = thrust::raw_pointer_cast(v.data());
 	real tolerance = 1e-6;
     int numberIterations = lanczos.run(dot, d_result, d_v, tolerance, size);
 	int iterations = 100;
@@ -64,24 +64,24 @@ Provide the solver with an instance of the functor and the target vector:
 The run function returns the number of iterations that were needed to achieve the requested accuracy.
 The runIterations returns the residual after the requested iterations.
 
-## Other functions:  
+## Other functions:
 
-After a certain number of iterations, if convergence was not achieved, the module will give up and throw an error.  
-To increase this threshold you can use this function:  
+After a certain number of iterations, if convergence was not achieved, the module will give up and throw an error.
+To increase this threshold you can use this function:
 ```c++
 lanczos::Solver::setIterationHardLimit(int newlimit);
 ```
-## Compilation:  
-This library requires lapacke and cblas (can be replaced by MKL). In GPU mode CUDA is also needed.  
-Note, however, that the heavy-weight of this solver comes from the Matrix-vector multiplication that must provided by the user. The main benefit of the CUDA mode is not an increased performance of the internal library code, but the fact that the input/output arrays will live in the GPU (saving potential memory copies).  
-## Optional macros:  
+## Compilation:
+This library requires lapacke and cblas (can be replaced by MKL). In GPU mode CUDA is also needed.
+Note, however, that the heavy-weight of this solver comes from the Matrix-vector multiplication that must provided by the user. The main benefit of the CUDA mode is not an increased performance of the internal library code, but the fact that the input/output arrays will live in the GPU (saving potential memory copies).
+## Optional macros:
 
-**CUDA_ENABLED**: Will compile a GPU enabled shared library, the solver expects input/output arrays to be in the GPU and most of the computations will be carried out in the GPU. Requires a working CUDA environment.  
-**DOUBLE_PRECISION**: The library is compiled in single precision by default. This macro switches to double precision, making ```lanczos::real``` be a typedef to double.  
-**USE_MKL**: Will include mkl.h instead of lapacke and cblas. You will have to modify the compilation flags accordingly.  
-**SHARED_LIBRARY_COMPILATION**: The Makefile uses this macro to compile a shared library. By default, this library is header only.  
+**CUDA_ENABLED**: Will compile a GPU enabled shared library, the solver expects input/output arrays to be in the GPU and most of the computations will be carried out in the GPU. Requires a working CUDA environment.
+**DOUBLE_PRECISION**: The library is compiled in single precision by default. This macro switches to double precision, making ```lanczos::real``` be a typedef to double.
+**USE_MKL**: Will include mkl.h instead of lapacke and cblas. You will have to modify the compilation flags accordingly.
+**SHARED_LIBRARY_COMPILATION**: The Makefile uses this macro to compile a shared library. By default, this library is header only.
 
-See the Makefile for further instructions.  
+See the Makefile for further instructions.
 
 ## Python interface
 
@@ -91,13 +91,13 @@ See python/example.py for more information.
 The root folder's Makefile will try to compile the python library as well. It expects pybind11 to be placed under the extern/ folder. Pybind11 is included as a submodule, so make sure to clone this repository with --recursive.
 Note that the python wrapper can only be compiled in CPU mode.
 
-## References:  
+## References:
+
+  [1] Krylov subspace methods for computing hydrodynamic interactions in Brownian dynamics simulations  J. Chem. Phys. 137, 064106 (2012); doi: 10.1063/1.4742347
 
-  [1] Krylov subspace methods for computing hydrodynamic interactions in Brownian dynamics simulations  J. Chem. Phys. 137, 064106 (2012); doi: 10.1063/1.4742347  
-  
-## Some notes:  
+## Some notes:
 
-  From what I have seen, this algorithm converges to an error of ~1e-3 in a few steps (<5) and from that point a lot of iterations are needed to lower the error.  
-  It usually achieves machine precision in under 50 iterations.  
+  From what I have seen, this algorithm converges to an error of ~1e-3 in a few steps (<5) and from that point a lot of iterations are needed to lower the error.
+  It usually achieves machine precision in under 50 iterations.
 
-  If the matrix does not have a sqrt (not positive definite, not symmetric...) it will usually be reflected as a nan in the current error estimation. In this case an exception will be thrown.  
+  If the matrix does not have a sqrt (not positive definite, not symmetric...) it will usually be reflected as a nan in the current error estimation. In this case an exception will be thrown.

example.cpp

@@ -22,8 +22,8 @@ using real = lanczos::real;
 struct DiagonalMatrix: public lanczos::MatrixDot{
   int size;
   DiagonalMatrix(int size): size(size){}
-  
-  void dot(real* v, real* Mv) override{
+
+  void operator()(real* v, real* Mv) override{
     //an example diagonal matrix
     for(int i=0; i<size; i++){
       Mv[i] = (2+i/10.0)*v[i]*2;

example.cu

@@ -23,7 +23,7 @@ struct DiagonalMatrix: public lanczos::MatrixDot{
   int size;
   DiagonalMatrix(int size): size(size){}
 
-  void dot(real* v, real* Mv) override{
+  void operator()(real* v, real* Mv) override{
     //An example diagonal matrix
     for(int i=0; i<size; i++){
       Mv[i] = (2+i/10.0)*v[i];

include/LanczosAlgorithm.cu

@@ -147,7 +147,7 @@ namespace lanczos{
 	auto d_w = detail::getRawPointer(w);	
 	/*w = D·vi*/
 	dot->setSize(N);
-	dot->dot(d_V+N*i, d_w);
+	dot->operator()(d_V+N*i, d_w);
 	if(i>0){
 	  /*w = w-h[i-1][i]·vi*/
 	  real alpha = -hsup[i-1];

include/LanczosAlgorithm.h

@@ -22,7 +22,8 @@ Some notes:
 #include"utils/defines.h"
 #include"utils/device_blas.h"
 #include<vector>
-#include<memory>
+#include <memory>
+#include<functional>
 #include"utils/device_container.h"
 #include"utils/MatrixDot.h"
 namespace lanczos{
@@ -39,13 +40,21 @@ namespace lanczos{
     int run(MatrixDot &dot, real *Bv, const real* v, real tolerance, int N){
       return run(&dot, Bv, v, tolerance, N);
     }
+    int run(std::function<void(real*, real*)> dot, real *Bv, const real* v, real tolerance, int N){
+      auto lanczos_dot = createMatrixDotAdaptor(dot);
+      return run(lanczos_dot, Bv, v, tolerance, N);
+    }
     //Given a Dotctor that computes a product M·v (where M is handled by Dotctor ), computes Bv = sqrt(M)·v
     //Returns the residual after numberIterations iterations
     //B = sqrt(M)
     real runIterations(MatrixDot *dot, real *Bz, const real*z, int numberIterations, int N);
     real runIterations(MatrixDot &dot, real *Bv, const real* v, int numberIterations, int N){
       return runIterations(&dot, Bv, v, numberIterations, N);
     }
+    real runIterations(std::function<void(real*, real*)> dot, real *Bv, const real* v, int numberIterations, int N){
+      auto lanczos_dot = createMatrixDotAdaptor(dot);
+      return runIterations(lanczos_dot, Bv, v, numberIterations, N);
+    }
 
     void setIterationHardLimit(int newLimit){this->iterationHardLimit = newLimit;}
 

include/utils/MatrixDot.h

@@ -1,13 +1,29 @@
 #ifndef LANCZOS_MATRIX_DOT_H
 #define LANCZOS_MATRIX_DOT_H
-#include"defines.h"
+#include "defines.h"
+#include<functional>
 namespace lanczos{
 
   struct MatrixDot{
     void setSize(int newsize){this->m_size = newsize;}
-    virtual void dot(real* v, real*Mv) = 0;
+    //virtual void dot(real* v, real*Mv) = 0;
+    virtual void operator()(real* v, real*Mv) = 0;
   protected:
     int m_size;
   };
+
+  //Transforms any callable into a MatrixDot valid to use with Lanczos
+  template<class Foo>
+  struct MatrixDotAdaptor: public lanczos::MatrixDot{
+    Foo& foo;
+    MatrixDotAdaptor(Foo &&foo):foo(foo){}
+    void operator()(real* v, real* Mv) override{foo(v,Mv);}    
+  };
+
+  template<class Foo>
+  auto createMatrixDotAdaptor(Foo &&foo){
+    return MatrixDotAdaptor<Foo>(foo);
+  }
+
 }
 #endif

python/python_wrapper.cpp

@@ -15,8 +15,12 @@ using namespace pybind11::literals;
 //This class allows to inherit lanczos::MatrixDot from python
 struct MatrixDotTrampoline: public lanczos::MatrixDot{  
   using MatrixDot::MatrixDot;
-  
-  void dot(real* v, real* Mv) override{
+
+  void dot(real* v, real* Mv){
+    this->operator()(v, Mv);
+  }
+
+  void operator()(real* v, real* Mv) override{
     pybind11::gil_scoped_acquire gil;  // Acquire the GIL while in this scope.
     // Try to look up the overridden method on the Python side.
     pybind11::function overridef = pybind11::get_override(this, "dot");
@@ -60,7 +64,7 @@ class PyLanczos{
 PYBIND11_MODULE(LANCZOS_PYTHON_NAME, m){
   py::class_<lanczos::MatrixDot, MatrixDotTrampoline>(m, "MatrixDot", "The virtual class required by the Lanczos solver").
     def(py::init<>()).
-    def("dot", &lanczos::MatrixDot::dot,
+    def("dot", &lanczos::MatrixDot::operator(),
 	"Given a result (Mv) and a vector (v), this method must write in Mv the result of multiplying the target matrix and v.",
 	"v"_a, "The input vector", "Mv"_a, "The output result vector");