# Weather and climate modeling codes from Fortran to C++

**URL:** <https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531>\
**Category:** Uncategorized\
**Created:** [July 9, 2021, 2:30am UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531 "2021-07-09T02:30:55Z")\
**Posts on this page:** 1\
**Showing post:** 41

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [March 24, 2024, 7:52pm UTC](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531/41 "2024-03-24T19:52:41Z")

</div>

> [@jkd2022](#):
>
> To quote from Eigen’s FAQ:
> 
> > Eigen has an incomparably better API than BLAS and LAPACK

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](https://eigen.tuxfamily.org/index.php?title=API_Showcase), Fortran has equivalents for most of them:

```cpp
// 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:

```fortran
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:

```fortran
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:

```cpp
// 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

```

```fortran
! 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:

```cpp
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:

```fortran
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:

```cpp
// 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"

```

```auto
! 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:

```txt
$ 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.

---

_[View the full topic](https://fortran-lang.discourse.group/t/weather-and-climate-modeling-codes-from-fortran-to-c/1531)._