Resistance to modernization

13

At least for BLAS, it may be possible to build an assumed-shape variant based upon the future C++ stdBLAS and mdspan container in combination the the C-Fortran array descriptor available in ISO_Fortran_binding.h.

Using GEMV interface from the Fortran 95 Thin BLAS as an example,

interface gemv
   subroutine cpp_dgemv(alpha,a,op_a,b,beta,c) bind(c)
      import c_double, blas_key
      real(c_double), intent(in), optional :: alpha, beta
      real(c_double), intent(in) :: a(:,:), b(:)
      integer(blas_key), optional :: op_a          ! blas_trans, or blas_notrans
      real(c_double), intent(inout) :: c(:)
   end subroutine
   ! ... other real kinds
end interface

I think it will be possible to define a family of routines for different kinds in C++ as shown in the draft below. (Disclaimer: the code below contains future language features and is meant only for demonstration purposes; I have not attempted any verification or validation whatsoever).

#include <mdspan> // P0009r16 - MDSPAN
#include <linalg> // P1673R9 - A free function linear algebra interfaces based on the BLAS

#include <ISO_Fortran_binding.h>

enum blas_key
{
    blas_notrans,
    blas_trans
};

template<typename T>
auto as_mdspan_2d(CFI_cdesc_t *obj) {

    using ext_t = std::extents<std::size_t,std::dynamic_extent,std::dynamic_extent>;
    ext_t extent{obj->dim[0].extent,obj->dim[1].extent};

    using mat_t = std::mdspan<std::size_t,ext_t,std::layout_left>;
    mat_t a{static_cast<T *>(obj->base_addr),extent};

    return a;
};

template<typename T>
auto as_mdspan(CFI_cdesc_t *obj) {

    using ext_t = std::extents<size_t,std::dynamic_extent,1>;
    ext_t extent{obj->dim[0].extent};

    using vec_t = std::mdspan<size_t,ext_t>;
    vec_t a{static_cast<T *>(obj->base_addr),extent};

    return a;
};

template<typename T>
void cpp_Tgemv(const T *alpha, const CFI_cdesc_t *A, const blas_key *op_a,
    const CFI_cdesc_t *b, const T *beta, CFI_cdesc_t *c) {

    using std::linalg::scaled;
    using std::linalg::transposed;
    using std::linalg::matrix_vector_product;

    T _alpha = (T) 1.0
    if (alpha) _alpha = *alpha;

    blas_key _op_a = blas_notrans;
    if (op_a) _op_a = *op_a;

    T _beta = (T) 0.0;
    if (beta) _beta = *beta;

    // assume a, b, c contiguous for now; for non-contiguous or strided
    // array submdspan is needed

    auto _A = as_mdspan_2d<T>(A);
    auto _b =    as_mdspan<T>(b);
    auto _c =    as_mdspan<T>(c);

    matrix_vector_product(
        scaled(_alpha, _op_a == blas_notrans ? _A : transposed(_A)),
        _b,scaled(_beta,_c),_c);
}

/* ... instantiate function templates for different 
       real kinds: float, c_double std::complex, ...
*/

I’m not sure if there is any communication between the C++ and Fortran committees on this topic.

(For brevity, handling of optional arguments can be deferred to a small function template.)

2 Likes