First commit

Author: RaulPPelaez

.gitignore

@@ -0,0 +1,6 @@
+*~
+a.out
+example
+*so
+*o
+compile_flags.txt

Makefile

@@ -0,0 +1,51 @@
+
+
+NVCC=nvcc
+CXX=g++
+
+#Uncomment for a GPU enabled library
+#CUDA_ENABLED=-DCUDA_ENABLED
+
+#Uncomment to compile in double precision mode
+#DOUBLE_PRECISION=-DDOUBLE_PRECISION
+
+
+INCLUDEFLAGS=-Iinclude/
+NVCCLDFLAGS= -lcublas
+LDFLAGS= -llapacke -lcblas
+
+LIBNAME=liblanczos.so
+
+
+CXXFLAGS=-fPIC -w -O3 -g -std=c++14 $(INCLUDEFLAGS) $(DOUBLE_PRECISION) $(CUDA_ENABLED)
+NVCCFLAGS=-ccbin=$(CXX) -Xcompiler "$(CXXFLAGS)" -std=c++14 -O3 $(INCLUDEFLAGS) $(DOUBLE_PRECISION) $(CUDA_ENABLED)
+
+ifndef CUDA_ENABLED
+COMPILER=$(CXX)
+CXXFLAGS:=$(CXXFLAGS)  -xc++
+else
+COMPILER=$(NVCC) 
+LDFLAGS:=$(LDFLAGS) $(NVCCLDFLAGS)
+CXXFLAGS:=$(NVCCFLAGS) -I$(CUDA_ROOT)/include
+endif
+
+all: shared $(patsubst %.cu, %, $(wildcard *.cu)) $(patsubst %.cpp, %, $(wildcard *.cpp))
+
+$(LIBNAME): $(wildcard include/*.cu)
+	$(COMPILER) -DSHARED_LIBRARY_COMPILATION -shared $(CXXFLAGS) $^ -o $@ $(LDFLAGS)
+
+
+shared: $(LIBNAME)
+
+
+%: %.cu Makefile
+	$(COMPILER) $(CXXFLAGS) $<  -o $@ $(LDFLAGS)
+
+
+
+#clean: $(patsubst include/%.cu, %.clean, $(wildcard include/*.cu))
+
+#%.clean:
+#rm -f $(@:.clean=.so)
+clean:
+	rm -rf include/*.o $(LIBNAME) example

README.md

@@ -0,0 +1,90 @@
+# Lanczos solver,
+  Computes the matrix-vector product sqrt(M)·v using a recursive algorithm.  
+  For that, it requires a functor in which the () operator takes an output real* array and an input real* (both in device memory if compiled in CUDA mode or host memory otherwise) as:  
+  ```c++ 
+  inline operator()(real* out_Mv, real * a_v);
+  ```  
+  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 
+
+## 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.  
+
+Let us go through the remaining one, a GPU-only example.  
+
+Create the module:
+```c++
+	real tolerance = 1e-6;
+    lanczos::Solver lanczos(tolerance);
+```
+Write a functor that computes the product between the original matrix and a given vector, "v":
+```c++
+//A functor that will return the result of multiplying a certain matrix times a given vector
+struct MatrixDot{
+  int size;
+  MatrixDot(int size): size(size){}
+  
+  void operator()(real* v, real* Mv){
+    //An example diagonal matrix
+    for(int i=0; i<size; i++){
+      Mv[i] = 2*v[i];
+    }
+  }
+
+};
+
+```
+
+Provide the solver with an instance of the functor and the target vector:  
+
+```c++
+    int size = 10;
+    //A vector filled with 1.
+    //Lanczos defines this type for convenience. It will be a thrust::device_vector if CUDA_ENABLED is defined and an std::vector otherwise
+    thrust::device_vector<real> v(size);
+    thrust::fill(v.begin(), v.end(), 1);
+    //A vector to store the result of sqrt(M)*v
+    thrust::device_vector<real> result(size);
+    //A functor that multiplies by the identity matrix times two
+    MatrixDot dot(size);
+    //Call the solver
+    real* d_result = thrust::raw_pointer_cast(result.data());
+    real* d_v = thrust::raw_pointer_cast(v.data());
+    int numberIterations = lanczos.solve(dot, d_result, d_v, size);
+```
+The solve function returns the number of iterations that were needed to achieve the requested accuracy.
+
+## 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:  
+```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:  
+
+**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.  
+
+## References:  
+
+  [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:  
+
+  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.  

example.cpp

@@ -0,0 +1,64 @@
+/*Raul P. Pelaez 2022. Lanczos solver CPU example
+This code will compute and print sqrt(M)*v using the iterative Krylov solver in [1].
+
+In this usage example, the matrix M is a diagonal matrix and v is filled with ones.
+
+This code is equivalent to example.cu, but can only be compiled in CPU mode.
+
+You can compile it with:
+g++ -std=c++14 -Iinclude example.cpp -llapacke -lcblas
+
+References:
+  [1] Krylov subspace methods for computing hydrodynamic interactions in Brownian dynamics simulations
+  J. Chem. Phys. 137, 064106 (2012); doi: 10.1063/1.4742347
+ */
+#include<iostream>
+#include <LanczosAlgorithm.h>
+
+//Using this floating point type (either float or double) will make the code be compiled with the same
+// precision as the library.
+using real = lanczos::real;
+//A functor that will return the result of multiplying a certain matrix times a given vector
+struct MatrixDot{
+  int size;
+  MatrixDot(int size): size(size){}
+  
+  void operator()(real* v, real* Mv){
+    //an example diagonal matrix
+    for(int i=0; i<size; i++){
+      Mv[i] = (2+i/10.0)*v[i];
+    }
+  }
+
+};
+
+
+int main(){
+  {
+    //Initialize the solver
+    real tolerance = 1e-6;
+    lanczos::Solver lanczos(tolerance);
+    int size = 10;
+    //A vector filled with 1.
+    //Lanczos defines this type for convenience. It will be a thrust::device_vector if CUDA_ENABLED is defined and an std::vector otherwise
+    std::vector<real> v(size);
+    std::fill(v.begin(), v.end(), 1);
+    //A vector to store the result of sqrt(M)*v
+    std::vector<real> result(size);
+    //A functor that multiplies by the identity matrix times two
+    MatrixDot dot(size);
+    //Call the solver
+    int numberIterations = lanczos.solve(dot, result.data(), v.data(), size);
+    std::cout<<"Solved after "<<numberIterations<< " iterations"<<std::endl;
+    //Now result is filled with sqrt(M)*v = sqrt(2)*[1,1,1...1]
+    std::cout<<"Result: ";for(int i = 0; i<10; i++) std::cout<<result[i]<<" "; std::cout<<std::endl;
+    //Compute error
+    std::cout<<"Error: ";
+    for(int i = 0; i<10; i++){
+      real truth = sqrt(2+i/10.0);
+      std::cout<<abs(result[i]-truth)/truth<<" ";
+    }
+    std::cout<<std::endl;
+  }
+  return 0;
+}

example.cu

@@ -0,0 +1,67 @@
+/*Raul P. Pelaez 2022. Lanczos solver CPU/GPU example
+This code will compute and print sqrt(M)*v using the iterative Krylov solver in [1].
+
+In this usage example, the matrix M is a diagonal matrix and v is filled with ones.
+
+The code below uses a series of utilities available in the library to make it work in either the CPU or GPU.
+
+If this code is compiled with CUDA_ENABLED defined, the code will run in a GPU. Otherwise it will be CPU only.
+
+References:
+  [1] Krylov subspace methods for computing hydrodynamic interactions in Brownian dynamics simulations
+  J. Chem. Phys. 137, 064106 (2012); doi: 10.1063/1.4742347
+ */
+#include<iostream>
+#include <LanczosAlgorithm.h>
+
+//Using this floating point type (either float or double) will make the code be compiled with the same
+// precision as the library.
+using real = lanczos::real;
+
+//A functor that will return the result of multiplying a certain matrix times a given vector
+struct MatrixDot{
+  int size;
+  MatrixDot(int size): size(size){}
+  
+  void operator()(real* v, real* Mv){
+    //An example diagonal matrix
+    for(int i=0; i<size; i++){
+      Mv[i] = (2+i/10.0)*v[i];
+    }
+  }
+
+};
+
+
+int main(){
+  {
+    //Initialize the solver
+    real tolerance = 1e-6;
+    lanczos::Solver lanczos(tolerance);
+    int size = 10;
+    //A vector filled with 1.
+    //Lanczos defines this type for convenience. It will be a thrust::device_vector if CUDA_ENABLED is defined and an std::vector otherwise
+    lanczos::device_container<real> v(size);
+    lanczos::detail::device_fill(v.begin(), v.end(), 1);
+    //A vector to store the result of sqrt(M)*v
+    lanczos::device_container<real> result(size);
+    //A functor that multiplies by the identity matrix times two
+    MatrixDot dot(size);
+    //Call the solver
+    real* d_result = lanczos::detail::getRawPointer(result);
+    real* d_v = lanczos::detail::getRawPointer(v);
+    int numberIterations = lanczos.solve(dot, d_result, d_v, size);
+    std::cout<<"Solved after "<<numberIterations<< " iterations"<<std::endl;
+    //Now result is filled with sqrt(M)*v = sqrt(2)*[1,1,1...1]
+    std::cout<<"Result: ";for(int i = 0; i<10; i++) std::cout<<result[i]<<" "; std::cout<<std::endl;
+    //Compute error
+    std::cout<<"Error: ";
+    for(int i = 0; i<10; i++){
+      real truth = sqrt(2+i/10.0);
+      std::cout<<abs(result[i]-truth)/truth<<" ";
+    }
+    std::cout<<std::endl;
+    
+  }
+  return 0;
+}