WLS and RBF-FD based stencil algorithms

[ivan-pi]

For stencil algorithms we can compute the WLS and RBF-FD weights in advance (using Medusa) and directly insert them into the streaming kernels in Fortran.

Alternatively, we can implement the streaming in C++, but we need to figure a way to pass Fortran allocated memory to C++ and wrap it in an Eigen vector. This can be done with the Eigen Map class.

[ivan-pi]

Turns out to be easier than I expected:

/* multiply_by_two.cc */

#include<Eigen/Dense>

#ifdef __cplusplus
extern "C" {
#endif

// Multiplies a C or Fortran array by 2
//
//  Inputs: 
//     nx, ny - dimensions of the array
//     a - pointer to an array of floats
//
void multiply_by_two(int nx, int ny, float* a) {

    Eigen::Map<Eigen::Matrix<float,Eigen::Dynamic,Eigen::Dynamic>> mf(a,nx,ny);
    mf *= 2;    

}

#ifdef __cplusplus
}
#endif

The Fortran driver

program pass_to_eigen

  use iso_c_binding, only: c_float, c_int
  implicit none

  interface
    subroutine multiply_by_two(nx,ny,a) bind(c,name="multiply_by_two")
      import c_float, c_int
      integer(c_int), intent(in), value :: nx, ny
      real(c_float), intent(inout) :: a(nx, ny)
    end subroutine
  end interface

  real(c_float) :: a(3,2)

  a = reshape([real(c_float) :: 1, 2, 3, 4, 5, 6], [3, 2])

  call multiply_by_two(3,2,a)

  print *, a

end program

After I installed Eigen using sudo apt install libeigen3-dev, I was able to run the Fortran main program using the commands:

$ g++ -c multiply_by_two.cc -I /usr/include/eigen3/
$ gfortran pass_to_eigen.f90 multiply_by_two.o -lstdc++
$ ./a.out
   2.00000000       4.00000000       6.00000000       8.00000000       10.0000000       12.0000000    

[ivan-pi]

It's also possible to pass a CSR or CSC matrix allocated in Fortran to Eigen. This way we could wrap the many solvers and preconditioners that are available in Eigen. The matrices we would allocate in Fortran are "compressed" according to Eigen terminology. We could also use Eigen as a backend for a Sparse BLAS like library with solver extensions.

The entry point is the class Map<SparseMatrix<>> which accepts the raw buffers of a CSR or CSC matrix format. An older (and deprecated) solution is described in this Stack Overflow thread.

Another interesting use of Eigen are the Matrix-free solvers. It should be possible to define a linear operator in Fortran, and use it in a matrix-free solver like BiCGSTAB.

[ivan-pi]

BiCGSTAB example:

#include <type_traits>
#include <iostream>

#include <Eigen/Dense>  
#include <Eigen/SparseCore>
#include <Eigen/IterativeLinearSolvers>

template<typename T, typename idx = int>
using CsrMat = Eigen::SparseMatrix<typename std::remove_const<T>::type, Eigen::RowMajor, idx>;

template<typename T>
using Vec = Eigen::Map<
    typename std::conditional<std::is_const<T>::value, 
    const Eigen::Matrix<typename std::remove_const<T>::type, Eigen::Dynamic, 1>, 
    Eigen::Matrix<T, Eigen::Dynamic, 1> >::type >;

template<typename T, typename idx = int>
using CsrMatrix = Eigen::Map< typename std::conditional<std::is_const<T>::value,
    const CsrMat<T,idx>, CsrMat<T,idx>>::type>;

extern "C" {

// Solve sparse system Ax = b, where A is a CSR matrix
void solve_csr_bicgstab(const int nr, const int nc, const int nnz, 
                        const int ia[], const int ja[], const double a[],
                        const double b[], double x[],
                        const int* maxIters,
                        const double* tolerance) {

    Vec<const double> b_{b,nr};
    Vec<double> x_{x,nr};

    CsrMatrix<const double> A_{nr,nc,nnz,ia,ja,a,nullptr};

    Eigen::BiCGSTAB<Eigen::SparseMatrix<double>> solver{A_};

    if (maxIters)
        solver.setMaxIterations((*maxIters));

    if (tolerance)
        solver.setTolerance((*tolerance));

    x_ = solver.solveWithGuess(b_,x_);

    std::cout << "#iterations:     " << solver.iterations() << std::endl;
    std::cout << "estimated error: " << solver.error()      << std::endl;

// To check for success:
// if (stat)
//      *stat = static_cast<int>(solver.info());

}

} // extern "C"

[ivan-pi]

One problem with using Eigen solvers from Fortran is they expect 0-based indexing. This means it would be necessary to subtract and add one before and after each call:

ia = ia - 1
ja = ja - 1
call solve_csr_bicgstab(nr,nc,nnz,ia,ja,a,b,x[,max_iters,tol,stat])
ia = ia + 1
ja = ja + 1

If we are calling the solver infrequently, this might be affordable. Otherwise we should probably keep a boolean flag and convert only when needed.

Potentially, it's easier to just re-implement the solver in Fortran: https://gitlab.com/libeigen/eigen/-/blob/master/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h