Weather and climate modeling codes from Fortran to C++

41

I always saw the main point of Eigen as providing matrix and vector classes, which are otherwise missing in C++.

Looking at some of the showcases of the Eigen API, Fortran has equivalents for most of them:

// performing row/column operations on a matrix
m.row(i) += alpha * m.row(j);
m.col(i) = alpha * m.col(j) + beta * m.col(k);
m.row(i).swap(m.row(j));
m.col(i) *= factor;

I’d say the overloaded in-place arithmetic operators are the nice thing here. In Fortran you’d use array expressions (or alternatively reach for BLAS), and hopefully the compiler won’t generate temporary array expressions:

M(i,:) = M(i,:) + alpha*M(j,:)
M(:,i) = alpha * M(:,j) + beta * M(:,k)

tmp = M(j,:)
M(j,:) = M(i,:)
M(i,:) = tmp

M(:,i) = factor*M(:,i)

Swapping rows can be accomplished with a helper subroutine:

call swap(a(1,:),a(2,:))

contains
    elemental subroutine swap(a,b)
        real, intent(inout) :: a, b
        real :: t
        t = a
        a = b
        b = t
    end subroutine

Computing various reductions in Fortran is easier than Eigen:

// computing sums with Eigen
result = m.sum(); // returns the sum of all coefficients in m
result = m.row(i).sum();
result = m.rowwise().sum(); // returns a vector of the sums in each row

// sum of the cubes of the i-th column of a matrix:
result = m.col(i).array().cube().sum(); // array() returns an array-expression

! computing sums in Fortran
result = sum(m)
result = sum(m(i,:))
result = sum(m,dim=2)
result = sum(m(:,i)**3)

The factorization syntax and solvers in Eigen are very convenient:

result = m.lu().solve(right_hand_side);            // using partial-pivoting LU
result = m.fullPivLu().solve(right_hand_side);     // using full-pivoting LU
result = m.householderQr().solve(right_hand_side); // using Householder QR
result = m.ldlt().solve(right_hand_side);          // using LDLt Cholesky

Fortran does not allow chaining operators, so this would become a two-step process:

type(matrix) :: M
M = matrix([...])
call M%lu(inplace=.true.)
result = M%solve(right_hand_side)

Nevertheless, thanks to Eigen’s powerful template-based design, it’s relatively easily to use from Fortran too. Here’s how to do it using F2018:

// eigen.cpp
#include <ISO_Fortran_binding.h> // F2018
#include <Eigen/Dense>

using namespace Eigen;

// Convert descriptor to matrix
template<typename T>
auto fmat(CFI_cdesc_t *M) -> Map<Matrix<T,Dynamic,Dynamic>> {
    return {static_cast<T *>(M->base_addr), M->dim[0].extent, M->dim[1].extent};
}

// Convert descriptor to column vector
template<typename T>
auto fvec(CFI_cdesc_t *x) -> Map<Matrix<T,Dynamic,1>> {
    return {static_cast<T *>(x->base_addr), x->dim[0].extent};
}

template<typename T>
void solve_lu(CFI_cdesc_t *A, CFI_cdesc_t *b)
{
    auto A_ = fmat<T>(A);
    auto b_ = fvec<T>(b);
    b_ = A_.lu().solve(b_);
}

extern "C" {

void sslvlu(CFI_cdesc_t *A, CFI_cdesc_t *b) { solve_lu<float>(A, b); }
void dslvlu(CFI_cdesc_t *A, CFI_cdesc_t *b) { solve_lu<double>(A, b); }

} // extern "C"

! test_eigen.f90

module eigen
use, intrinsic :: iso_c_binding, only: c_double, c_float
implicit none
interface solve
   subroutine sslvlu(A,b) bind(c)
      import c_float
      real(c_float), intent(in), contiguous :: A(:,:)
      real(c_float), intent(inout), contiguous :: b(:)
   end subroutine
   subroutine dslvlu(A,b) bind(c)
      import c_double
      real(c_double), intent(in), contiguous :: A(:,:)
      real(c_double), intent(inout), contiguous :: b(:)
   end subroutine
end interface
end module

program main
use eigen
implicit none
real(c_float) :: A(3,3)
real(c_float) :: b(3), x(3)
A = reshape([1,-1,2,-2,3,-5,3,-1,5],[3,3])
b = [9,-6,17]
x = b
call solve(A,x)
print *, "Solution: ", x
print *, "Norm: ", norm2(matmul(A,x) - b)
end program

Commands to run the example:

$ g++-13 -c -O2 eigen.cpp -I/usr/include/eigen3/
$ gfortran -Wall -O2 test_eigen.f90 eigen.o -lstdc++
$ ./a.out
 Solution:    1.00000000      -1.00000000       2.00000000    
 Norm:    0.00000000

In this manner you can wrap all of Eigen’s dense matrix decompositions in an evening.

4 Likes