cs_array.h

Go to the documentation of this file.

#ifndef CS_ARRAY_H
#define CS_ARRAY_H
 
/*============================================================================
 * Array handling utilities.
 *============================================================================*/
 
/*
  This file is part of code_saturne, a general-purpose CFD tool.
 
  Copyright (C) 1998-2026 EDF S.A.
 
  This program is free software; you can redistribute it and/or modify it under
  the terms of the GNU General Public License as published by the Free Software
  Foundation; either version 2 of the License, or (at your option) any later
  version.
 
  This program is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
  details.
 
  You should have received a copy of the GNU General Public License along with
  this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
  Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
 
/*----------------------------------------------------------------------------*/
 
/*----------------------------------------------------------------------------
 * Standard C and C++ library headers
 *----------------------------------------------------------------------------*/
 
#include <string.h>
 
/*----------------------------------------------------------------------------
 *  Local headers
 *----------------------------------------------------------------------------*/
 
#include "base/cs_defs.h"
#include "base/cs_base_accel.h"
#if defined (__CUDACC__)
#include "base/cs_array_cuda.h"
#include "base/cs_base_cuda.h"
#endif
#if defined (__HIPCC__)
#include "base/cs_array_hip.h"
#include "base/cs_base_hip.h"
#endif
 
#include "base/cs_dispatch.h"
#include "base/cs_mdspan.h"
#include "base/cs_parall.h"
 
#if defined(__cplusplus)
#include <limits>
#endif
/*----------------------------------------------------------------------------*/
 
/*=============================================================================
 * Macro definitions
 *============================================================================*/
 
/* Define the way to apply a subset. In case of a copy, the subset is related
   to the reference (input array) and/or the destination array (output
   array) */
 
#define CS_ARRAY_SUBSET_NULL  -1
#define CS_ARRAY_SUBSET_IN     0
#define CS_ARRAY_SUBSET_OUT    1
#define CS_ARRAY_SUBSET_INOUT  2
 
/*============================================================================
 * Type definitions
 *============================================================================*/
 
/*============================================================================
 * Global variables
 *============================================================================*/
 
/*=============================================================================
 * Public inline function prototypes
 *============================================================================*/
 
/*=============================================================================
 * Public function prototypes
 *============================================================================*/
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
template <typename T, size_t stride, typename... Arrays>
void
cs_arrays_set_value(const cs_lnum_t  n_elts,
                    const T         *ref_val,
                    Arrays&&...      arrays)
{
  /* Expand the parameter pack */
  T* array_ptrs[] = {arrays ... };
 
#if defined (__CUDACC__) || defined(__HIPCC__)
  bool is_available_on_device = cs_check_device_ptr(ref_val);
  for (T* array : array_ptrs)
    is_available_on_device =  is_available_on_device
                           && (cs_check_device_ptr(array)
                               == CS_ALLOC_HOST_DEVICE_SHARED);
 
  if (is_available_on_device) {
#if defined (__CUDACC__)
    cudaStream_t stream_ = cs_cuda_get_stream(0);
#else
    hipStream_t stream_ = cs_hip_get_stream(0);
#endif
    cs_arrays_set_value<T, stride>(stream_,
                                   false,
                                   n_elts,
                                   ref_val,
                                   arrays...);
    return;
  }
#endif
 
  auto set_value = [=](cs_lnum_t i_elt) {
    for (T* array : array_ptrs)
      memcpy(array + i_elt*stride, ref_val, stride*sizeof(T));
  };
 
  /* Loop through each index and assign values */
# ifdef _OPENMP
  int n_threads = cs_parall_n_threads(n_elts, CS_THR_MIN);
# pragma omp parallel for  num_threads(n_threads)
# endif
  for (cs_lnum_t i = 0; i < n_elts; i++)
    set_value(i);
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
template <typename T, size_t stride, typename... Arrays>
void
cs_arrays_set_value(const cs_lnum_t  n_elts,
                    const T          ref_val,
                    Arrays&&...      arrays)
{
  /* Expand the parameter pack */
  T* array_ptrs[] = {arrays ... };
 
#if defined (__CUDACC__) || defined (__HIPCC__)
  bool is_available_on_device = true;
  for (T* array : array_ptrs)
    is_available_on_device =  is_available_on_device
                           && (cs_check_device_ptr(array)
                               == CS_ALLOC_HOST_DEVICE_SHARED);
 
  if (is_available_on_device) {
#if defined (__CUDACC__)
    cudaStream_t stream_ = cs_cuda_get_stream(0);
#else
    hipStream_t stream_ = cs_hip_get_stream(0);
#endif
    cs_arrays_set_value<T, stride>(stream_,
                                   false,
                                   n_elts,
                                   ref_val,
                                   arrays...);
    return;
  }
#endif
 
  auto set_value = [=](cs_lnum_t i_elt) {
    for (T* array : array_ptrs)
      for (size_t k = 0; k < stride; k++) {
        array[i_elt*stride + k] = ref_val;
      }
  };
 
  /* Loop through each index and assign values */
# ifdef _OPENMP
  int n_threads = cs_parall_n_threads(n_elts, CS_THR_MIN);
# pragma omp parallel for  num_threads(n_threads)
# endif
  for (cs_lnum_t i = 0; i < n_elts; i++)
    set_value(i);
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
template <typename T, size_t stride, typename... Arrays>
void
cs_arrays_set_value(cs_dispatch_context  &ctx,
                    const cs_lnum_t       n_elts,
                    const T              *ref_val,
                    Arrays&&...           arrays)
{
  /* Expand the parameter pack */
  T* array_ptrs[] = {arrays ... };
 
#if defined (__CUDACC__) || defined (__HIPCC__)
  if (ctx.use_gpu()) {
    cs_stream_t stream_ = ctx.stream();
    cs_arrays_set_value<T, stride>(stream_,
                                   true,
                                   n_elts,
                                   ref_val,
                                   arrays...);
    return;
  }
#endif
 
  auto set_value = [=](cs_lnum_t i_elt) {
    for (T* array : array_ptrs)
      memcpy(array + i_elt*stride, ref_val, stride*sizeof(T));
  };
 
  /* Loop through each index and assign values */
# ifdef _OPENMP
  int n_threads = ctx.n_cpu_threads();
  if (n_threads < 0)
    n_threads = cs_parall_n_threads(n_elts, CS_THR_MIN);
# pragma omp parallel for  num_threads(n_threads)
# endif
  for (cs_lnum_t i = 0; i < n_elts; i++)
    set_value(i);
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
template <typename T, size_t stride, typename... Arrays>
void
cs_arrays_set_value(cs_dispatch_context  &ctx,
                    const cs_lnum_t       n_elts,
                    const T               ref_val,
                    Arrays&&...           arrays)
{
  /* Expand the parameter pack */
  T* array_ptrs[] = {arrays ... };
 
#if defined (__CUDACC__) || defined (__HIPCC__)
  if (ctx.use_gpu()) {
    cs_stream_t stream_ = ctx.stream();
    cs_arrays_set_value<T, stride>(stream_,
                                   true,
                                   n_elts,
                                   ref_val,
                                   arrays...);
    return;
  }
#endif
 
  auto set_value = [=](cs_lnum_t i_elt) {
    for (T* array : array_ptrs)
      for (size_t k = 0; k < stride; k++) {
        array[i_elt*stride + k] = ref_val;
      }
  };
 
  /* Loop through each index and assign values */
# ifdef _OPENMP
  int n_threads = ctx.n_cpu_threads();
  if (n_threads < 0)
    n_threads = cs_parall_n_threads(n_elts, CS_THR_MIN);
# pragma omp parallel for  num_threads(n_threads)
# endif
  for (cs_lnum_t i = 0; i < n_elts; i++)
    set_value(i);
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
template <typename T, size_t stride, typename... Arrays>
void
cs_arrays_set_zero(cs_dispatch_context  &ctx,
                   const cs_lnum_t       n_elts,
                   Arrays&&...           arrays)
{
  /* Expand the parameter pack */
  T* array_ptrs[] = {arrays ... };
 
#if defined (__CUDACC__) || defined (__HIPCC__)
  if (ctx.use_gpu()) {
    cs_stream_t stream_ = ctx.stream();
    cs_arrays_set_zero<T, stride>(stream_,
                                  true,
                                  n_elts,
                                  arrays...);
    return;
  }
#endif
 
  constexpr T c = static_cast<T>(0);
 
  cs_lnum_t n_vals = n_elts*stride;
 
# ifdef _OPENMP
  int n_threads = ctx.n_cpu_threads();
  if (n_threads < 0)
    n_threads = cs_parall_n_threads(n_vals, CS_THR_MIN);
# pragma omp parallel num_threads(n_threads)
# endif
  {
    for (T* array : array_ptrs) {
#     ifdef _OPENMP
#     pragma omp for nowait
#     endif
      for (cs_lnum_t i = 0; i < n_vals; i++)
        array[i] = c;
    }
  };
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
template <typename T, size_t stride, typename... Arrays>
void
cs_arrays_set_value_on_subset(const cs_lnum_t  n_elts,
                              const cs_lnum_t *elt_ids,
                              const T         *ref_val,
                              Arrays&&...      arrays)
{
  if (n_elts < 1)
    return;
 
  if (elt_ids == NULL)
    cs_arrays_set_value<T, stride>(n_elts, ref_val, arrays...);
  else {
    /* Expand the parameter pack */
    T* array_ptrs[] = {arrays ... };
 
    auto set_value = [=](cs_lnum_t i_elt) {
      for (T* array : array_ptrs)
        memcpy(array + i_elt*stride, ref_val, stride*sizeof(T));
    };
 
    /* Loop through each index on the subset and assign values */
#   ifdef _OPENMP
    int n_threads = cs_parall_n_threads(n_elts, CS_THR_MIN);
#   pragma omp parallel for  num_threads(n_threads)
#   endif
    for (cs_lnum_t i = 0; i < n_elts; i++)
      set_value(elt_ids[i]);
  }
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
template <typename T, size_t stride, typename... Arrays>
void
cs_arrays_set_value_on_subset(const cs_lnum_t  n_elts,
                              const cs_lnum_t *elt_ids,
                              const T          ref_val,
                              Arrays&&...      arrays)
{
  if (n_elts < 1)
    return;
 
  if (elt_ids == NULL)
    cs_arrays_set_value<T, stride>(n_elts, ref_val, arrays...);
  else {
 
    /* Expand the parameter pack */
    T* array_ptrs[] = {arrays ... };
 
    auto set_value = [=](cs_lnum_t i_elt) {
      for (T* array : array_ptrs)
        for (size_t k = 0; k < stride; k++) {
          array[i_elt*stride + k] = ref_val;
        }
    };
 
   /* Loop through each index on the subset and assign values */
#   ifdef _OPENMP
    int n_threads = cs_parall_n_threads(n_elts, CS_THR_MIN);
#   pragma omp parallel for  num_threads(n_threads)
#   endif
    for (cs_lnum_t i = 0; i < n_elts; i++)
      set_value(elt_ids[i]);
  }
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
template <typename T>
void
cs_array_copy(const cs_lnum_t  size,
              const T*         src,
              T*               dest)
{
#if defined (__CUDACC__) || defined (__HIPCC__)
  bool is_available_on_device =  cs_check_device_ptr(src)
                              && cs_check_device_ptr(dest);
 
  if (is_available_on_device) {
#if defined (__CUDACC__)
    cudaStream_t stream_ = cs_cuda_get_stream(0);
#else
    hipStream_t stream_ = cs_hip_get_stream(0);
#endif
    cs_array_copy<T>(stream_, false, size, src, dest);
    return;
  }
#endif
 
# ifdef _OPENMP
  int n_threads = cs_parall_n_threads(size, CS_THR_MIN);
# pragma omp parallel for  num_threads(n_threads)
# endif
  for (cs_lnum_t ii = 0; ii < size; ii++)
    dest[ii] = src[ii];
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
template <typename T>
void
cs_array_difference(cs_lnum_t   size,
                    const T    *x,
                    const T    *y,
                    T          *diff)
{
  cs_array_copy(size, x, diff);
 
# ifdef _OPENMP
  int n_threads = cs_parall_n_threads(size, CS_THR_MIN);
# pragma omp parallel for  num_threads(n_threads)
# endif
  for (cs_lnum_t ii = 0; ii < size; ii++)
    diff[ii] -= y[ii];
}
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign true to all elements of an array. Case of an array of booleans
 *
 * \param[in]      size    total number of elements to set
 * \param[in, out] a       array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_bool_fill_true(cs_lnum_t  size,
                        bool       a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign false to all elements of an array. Case of an array of booleans
 *
 * \param[in]      size    total number of elements to set
 * \param[in, out] a       array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_bool_fill_false(cs_lnum_t  size,
                         bool       a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign zero to all elements of an array. Case of a cs_flag_t array.
 *
 * \param[in]      size    total number of elements to set to zero
 * \param[in, out] a       array of flags to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_flag_fill_zero(cs_lnum_t  size,
                        cs_flag_t  a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign zero to all elements of an array. Case of a cs_lnum_t array.
 *
 * \param[in]      size    total number of elements to set to zero
 * \param[in, out] a       array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_lnum_fill_zero(cs_lnum_t  size,
                        cs_lnum_t  a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Copy values from an array of cs_lnum_t type to another of the same
 *        dimensions.
 *
 * \param[in]  size  number of elements * dimension
 * \param[in]  src   source array values (size: n_elts*dim)
 * \param[out] dest  destination array values (size: n_elts*dim)
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_lnum_copy(cs_lnum_t       size,
                   const cs_lnum_t src[],
                   cs_lnum_t       dest[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign the value "num" to all elements of an array. Case of a
 *        cs_lnum_t array.
 *
 * \param[in]      size    total number of elements to set
 * \param[in]      num     value to set
 * \param[in, out] a       array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_lnum_set_value(cs_lnum_t  size,
                        cs_lnum_t  num,
                        cs_lnum_t  a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign the value "num" to an array on a selected subset of elements.
 *        if elt_ids is null, then one recovers the function
 *        \ref cs_array_lnum_set_value
 *
 * \param[in]      n_elts   number of elements
 * \param[in]      elt_ids  list of ids in the subset or null (size: n_elts)
 * \param[in]      num      value to set
 * \param[in, out] a        array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_lnum_set_value_on_subset(cs_lnum_t        n_elts,
                                  const cs_lnum_t  elt_ids[],
                                  cs_lnum_t        num,
                                  cs_lnum_t        a[]);
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign zero to all elements of an array. Case of a int array.
 *
 * \param[in]      size    total number of elements to set to zero
 * \param[in, out] a       array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_int_fill_zero(cs_lnum_t  size,
                       int        a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign the value "num" to all elements of an array. Case of a
 *        int array.
 *
 * \param[in]      size    total number of elements to set
 * \param[in]      num     value to set
 * \param[in, out] a       array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_int_set_value(cs_lnum_t  size,
                       int        num,
                       int        a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign the value "num" to an array on a selected subset of elements.
 *        if elt_ids is null, then one recovers the function
 *        \ref cs_array_int_set_value
 *
 * \param[in]      n_elts   number of elements
 * \param[in]      elt_ids  list of ids in the subset or null (size: n_elts)
 * \param[in]      num      value to set
 * \param[in, out] a        array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_int_set_value_on_subset(cs_lnum_t        n_elts,
                                 const cs_lnum_t  elt_ids[],
                                 int              num,
                                 int              a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Copy an array ("ref") into another array ("dest") on possibly only a
 *        part of the array(s). Array with stride > 1 are assumed to be
 *        interlaced. The subset of elements on which working is defined by
 *        "elt_ids". The way to apply the subset is defined with the parameter
 *        "mode" as follows:
 *        - Only the "ref" array if mode = 0 (CS_ARRAY_SUBSET_IN)
 *        - Only the "dest" array if mode = 1 (CS_ARRAY_SUBSET_OUT)
 *        - Both "ref" and "dest" arrays if mode = 2 (CS_ARRAY_SUBSET_INOUT)
 *
 *        It elt_ids is null or mode < 0 (CS_ARRAY_SUBSET_NULL), then the
 *        behavior is as \ref cs_array_real_copy
 *
 *        One assumes that all arrays are allocated with a correct size.
 *
 * \param[in]      n_elts   number of elements in the array
 * \param[in]      stride   number of values for each element
 * \param[in]      mode     apply the subset ids to which array(s)
 * \param[in]      elt_ids  list of ids in the subset or null (size: n_elts)
 * \param[in]      ref      reference values to copy
 * \param[in, out] dest     array storing values after applying the indirection
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_copy_subset(cs_lnum_t         n_elts,
                          int               stride,
                          const cs_lnum_t   elt_ids[],
                          int               mode,
                          const cs_real_t   ref[],
                          cs_real_t         dest[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Copy real values from an array to another of the same dimensions.
 *
 * \param[in]   size    number of elements * dimension
 * \param[in]   src     source array values (size: n_elts*dim)
 * \param[out]  dest    destination array values (size: n_elts*dim)
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_copy(cs_lnum_t         size,
                   const cs_real_t   src[],
                   cs_real_t         dest[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Multiply each value by a scaling factor s.t. dest *= scaling_factor
 *        If elt_ids is non-null, one applies an indirection.
 *        A stride can also be applied. One assumes an interlaced array.
 *
 * \param[in]  n_elts          number of elements
 * \param[in]  stride          number of values for each element
 * \param[in]  elt_ids         list of ids in the subset or null (size: n_elts)
 * \param[in]  scaling_factor  value of the scaling factor
 * \param[out] dest            destination array values
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_scale(cs_lnum_t         n_elts,
                    int               stride,
                    const cs_lnum_t  *elt_ids,
                    cs_real_t         scaling_factor,
                    cs_real_t         dest[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Add in place an array s.t. r += l_add
 *
 * \param[in]  n_elts   number of elements
 * \param[in]  l_add    array to add
 * \param[out] r        destination array values
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_padd(cs_lnum_t       n_elts,
                   const cs_real_t l_add[],
                   cs_real_t       r[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a constant value of dim "stride" to an interlaced array
 *        sharing the same stride
 *
 * \param[in]      n_elts    number of elements
 * \param[in]      stride    number of values for each element
 * \param[in]      ref_val   list of values to assign (size: stride)
 * \param[in, out] a         array to set (size: n_elts*stride)
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_value(cs_lnum_t        n_elts,
                        int              stride,
                        const cs_real_t  ref_val[],
                        cs_real_t        a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a weighted constant value of dim "stride" to an interlaced
 *        array sharing the same stride. Apply a weight for each element. This
 *        weight is constant for each component of an element.
 *
 * \param[in]      n_elts    number of elements
 * \param[in]      stride    number of values for each element
 * \param[in]      ref_val   list of values to assign (size: stride)
 * \param[in]      weight    values of the weight to apply (size: n_elts)
 * \param[in, out] a         array to set (size: n_elts*stride)
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_wvalue(cs_lnum_t        n_elts,
                         int              stride,
                         const cs_real_t  ref_val[],
                         const cs_real_t  weight[],
                         cs_real_t        a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a constant value of dim "stride" to an interlaced array
 *        sharing the same stride. Only a subset of elements are considered.
 *        If elt_ids is null, then one recovers the function
 *        \ref cs_array_real_set_value
 *
 * \param[in]      n_elts    number of elements
 * \param[in]      stride    number of values for each element
 * \param[in]      elt_ids   list of ids in the subset or null (size: n_elts)
 * \param[in]      ref_val   list of values to assign (size: stride)
 * \param[in, out] a         array to set (size >= n_elts * stride)
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_value_on_subset(cs_lnum_t        n_elts,
                                  int              stride,
                                  const cs_lnum_t  elt_ids[],
                                  const cs_real_t  ref_val[],
                                  cs_real_t        a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a weighted constant value of dim "stride" to an interlaced
 *        array sharing the same stride. Only a subset of elements are
 *        considered.  If elt_ids is null, then one recovers the function \ref
 *        cs_array_real_set_wvalue Apply a weight for each element. This
 *        weight is constant for each component of an element.
 *
 * \param[in]      n_elts   number of elements
 * \param[in]      stride   number of values for each element
 * \param[in]      elt_ids  list of ids in the subset or null (size: n_elts)
 * \param[in]      ref_val  list of values to assign (size: stride)
 * \param[in]      weight   values of the weight to apply (size >= n_elts)
 * \param[in, out] a        array to set (size >= n_elts*stride)
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_wvalue_on_subset(cs_lnum_t        n_elts,
                                   int              stride,
                                   const cs_lnum_t  elt_ids[],
                                   const cs_real_t  ref_val[],
                                   const cs_real_t  weight[],
                                   cs_real_t        a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a constant scalar value to an array.
 *
 * \param[in]      n_elts   number of elements
 * \param[in]      ref_val  value to assign
 * \param[in, out] a        array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_scalar(cs_lnum_t  n_elts,
                         cs_real_t  ref_val,
                         cs_real_t  a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a weighted constant scalar value to an array.
 *        The weight array has the same size as the array "a".
 *
 * \param[in]      n_elts   number of elements
 * \param[in]      ref_val  value to assign
 * \param[in]      weight   values of the weight to apply
 * \param[in, out] a        array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_wscalar(cs_lnum_t        n_elts,
                          cs_real_t        ref_val,
                          const cs_real_t  weight[],
                          cs_real_t        a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a constant scalar value to an array on a selected subset of
 *        elements. If elt_ids is null, then one recovers the function
 *        cs_array_real_set_scalar
 *
 * \param[in]      n_elts   number of elements
 * \param[in]      elt_ids  list of ids in the subset or null (size: n_elts)
 * \param[in]      ref_val  value to assign
 * \param[in, out] a        array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_scalar_on_subset(cs_lnum_t        n_elts,
                                   const cs_lnum_t  elt_ids[],
                                   cs_real_t        ref_val,
                                   cs_real_t        a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a weighted constant scalar value to an array on a selected
 *        subset of elements. If elt_ids is null, then one recovers the function
 *        cs_array_real_set_wscalar
 *
 * \param[in]      n_elts   number of elements
 * \param[in]      elt_ids  list of ids in the subset or null (size: n_elts)
 * \param[in]      ref_val  value to assign
 * \param[in]      weight   values of weights to apply
 * \param[in, out] a        array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_wscalar_on_subset(cs_lnum_t        n_elts,
                                    const cs_lnum_t  elt_ids[],
                                    cs_real_t        ref_val,
                                    const cs_real_t  weight[],
                                    cs_real_t        a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a constant vector to an array of stride 3 which is interlaced
 *
 * \param[in]      n_elts   number of elements
 * \param[in]      ref_val  vector to assign
 * \param[in, out] a        array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_vector(cs_lnum_t         n_elts,
                         const cs_real_t   ref_val[3],
                         cs_real_t         a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a weighted constant vector value to an interlaced array (of
 *        stride 3). The array of weights has the same size as the array "a".
 *
 * \param[in]      n_elts   number of elements
 * \param[in]      ref_val  vector to assign
 * \param[in]      weight   values of the weight to apply
 * \param[in, out] a        array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_wvector(cs_lnum_t        n_elts,
                          const cs_real_t  ref_val[3],
                          const cs_real_t  weight[],
                          cs_real_t        a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a constant vector to an interlaced array (of stride 3) on a
 *        selected subset of elements. If elt_ids is null, then one recovers the
 *        function cs_array_real_set_vector
 *
 * \param[in]      n_elts   number of elements
 * \param[in]      elt_ids  list of ids in the subset or null (size: n_elts)
 * \param[in]      ref_val  vector to assign
 * \param[in, out] a        array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_vector_on_subset(cs_lnum_t        n_elts,
                                   const cs_lnum_t  elt_ids[],
                                   const cs_real_t  ref_val[3],
                                   cs_real_t        a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a weighted constant vector value to an interlaced array (of
 *        stride 3). The subset selection is given by elt_ids. If null, then
 *        one recovers the function \ref cs_array_real_set_wvector
 *        The array of weights has the same size as the array "a".
 *
 * \param[in]      n_elts   number of elements
 * \param[in]      elt_ids  list of ids in the subset or null (size: n_elts)
 * \param[in]      ref_val  vector to assign
 * \param[in]      weight   values of the weight to apply
 * \param[in, out] a        array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_wvector_on_subset(cs_lnum_t        n_elts,
                                    const cs_lnum_t  elt_ids[],
                                    const cs_real_t  ref_val[3],
                                    const cs_real_t  weight[],
                                    cs_real_t        a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a constant 3x3 tensor to an array (of stride 9) which is
 *        interlaced
 *
 * \param[in]      n_elts    number of elements
 * \param[in]      ref_tens  tensor to assign
 * \param[in, out] a         array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_tensor(cs_lnum_t        n_elts,
                         const cs_real_t  ref_tens[3][3],
                         cs_real_t        a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a constant 3x3 tensor to an interlaced array (of stride 9) on
 *        a subset of elements. If elt_ids is null, then one recovers the
 *        function cs_array_real_set_tensor
 *
 * \param[in]      n_elts    number of elements
 * \param[in]      elt_ids   list of ids defining the subset or nullptr
 * \param[in]      ref_tens  tensor to assign
 * \param[in, out] a         array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_set_tensor_on_subset(cs_lnum_t         n_elts,
                                   const cs_lnum_t   elt_ids[],
                                   const cs_real_t   ref_tens[3][3],
                                   cs_real_t         a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign zero to all elements of an array.
 *
 * \param[in]      size    total number of elements to set to zero
 * \param[in, out] a       array to set
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_real_fill_zero(cs_lnum_t  size,
                        cs_real_t  a[]);
 
/*----------------------------------------------------------------------------*/
/*
 * \brief Assign a constant value to an array (deprecated function).
 *
 * \param[in]   n_elts  number of associated elements
 * \param[in]   dim     associated dimension
 * \param[in]   v       value to assign
 * \param[out]  a       array values (size: n_elts*dim)
 */
/*----------------------------------------------------------------------------*/
 
void
cs_array_set_value_real(cs_lnum_t  n_elts,
                        cs_lnum_t  dim,
                        cs_real_t  v,
                        cs_real_t  a[]);
 
/*----------------------------------------------------------------------------*/
 
namespace cs {
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
template<class T, int N = 1, layout L = layout::right>
class array : public mdspan<T, N, L> {
 
private:
  using _span = mdspan<T,N,L>;
 
public:
  using mdspan<T,N,L>::copy_data;
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void *
  operator new
  (
    std::size_t count,   
#if (defined(__GNUC__) || defined(__clang__)) \
    && !defined(CS_AOCC_BUILTIN_WORKAROUND) \
    && __has_builtin(__builtin_LINE) && __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  {
    return cs_mem_malloc_hd(CS_ALLOC_HOST,
                            1,
                            count,
                            "cs::array",
                            file_name,
                            line_number);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void *
  operator new[]
  (
    std::size_t count,   
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
 __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  {
    return cs_mem_malloc_hd(CS_ALLOC_HOST,
                            1,
                            count,
                            "cs::array",
                            file_name,
                            line_number);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  operator delete(void *ptr)
  {
    CS_FREE(ptr);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  operator delete[](void *ptr)
  {
    CS_FREE(ptr);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  array()
  {
    _owner = true;
#if !defined(__CUDA_ARCH__) && \
    !defined(SYCL_LANGUAGE_VERSION) && \
    !defined(__HIP_DEVICE_COMPILE__)
    _mode = cs_alloc_mode;
#else
    // use default and avoid compiler warnings
    _mode = CS_ALLOC_HOST_DEVICE_SHARED;
#endif
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  array
  (
    cs_lnum_t     size, 
    cs_alloc_mode_t alloc_mode = cs_alloc_mode, 
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
 __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  :
    _owner(true),
    _mode(alloc_mode)
  {
    /* Only usable for array */
    static_assert(N == 1,
                  "Instantiation using only total size only possible for "
                  " cs::array<T,1> or cs::array<T>");
    _span::set_size_(size);
    allocate_(file_name, line_number);
#if defined(CS_ARRAY_DBG_INIT)
    set_to_nan_();
#endif
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  array
  (
    const cs_lnum_t(&dims)[N], 
    cs_alloc_mode_t alloc_mode = cs_alloc_mode, 
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
 __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  :
    _owner(true),
    _mode(alloc_mode)
  {
    _span::set_size_(dims);
    allocate_(file_name, line_number);
#if defined(CS_ARRAY_DBG_INIT)
    set_to_nan_();
#endif
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  array
  (
    cs_lnum_t     size1, 
    cs_lnum_t     size2, 
    cs_alloc_mode_t alloc_mode = cs_alloc_mode, 
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
 __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  :
    _owner(true),
    _mode(alloc_mode)
  {
    /* Only usable for array */
    static_assert(N == 2,
                  "Instantiation using (size1, size2) only possible for "
                  " cs::array<T,2>");
    cs_lnum_t tmp_size[N] = {size1, size2};
    _span::set_size_(tmp_size);
    allocate_(file_name, line_number);
#if defined(CS_ARRAY_DBG_INIT)
    set_to_nan_();
#endif
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  array
  (
    cs_lnum_t     size1, 
    cs_lnum_t     size2, 
    cs_lnum_t     size3, 
    cs_alloc_mode_t alloc_mode = cs_alloc_mode, 
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
 __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  :
    _owner(true),
    _mode(alloc_mode)
  {
    /* Only usable for array */
    static_assert(N == 3,
                  "Instantiation using (size1, size2, size3) only possible "
                  " for cs::array<T,3>");
    cs_lnum_t tmp_size[N] = {size1, size2, size3};
    _span::set_size_(tmp_size);
    allocate_(file_name, line_number);
#if defined(CS_ARRAY_DBG_INIT)
    set_to_nan_();
#endif
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  array
  (
    cs_lnum_t     size1, 
    cs_lnum_t     size2, 
    cs_lnum_t     size3, 
    cs_lnum_t     size4, 
    cs_alloc_mode_t alloc_mode = cs_alloc_mode, 
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
 __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  :
    _owner(true),
    _mode(alloc_mode)
  {
    /* Only usable for array */
    static_assert(N == 4,
                  "Instantiation using (size1, size2, size3, size4) only possible "
                  " for cs::array<T,4>");
    cs_lnum_t tmp_size[N] = {size1, size2, size3, size4};
    _span::set_size_(tmp_size);
    allocate_(file_name, line_number);
#if defined(CS_ARRAY_DBG_INIT)
    set_to_nan_();
#endif
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  template<typename... Args>
  CS_F_HOST_DEVICE
  array
  (
    T       *data,   
    Args...  indices 
  )
  :
    _span(data, indices...),
    _owner(false)
  {
#if !defined(__CUDA_ARCH__) && \
    !defined(SYCL_LANGUAGE_VERSION) && \
    !defined(__HIP_DEVICE_COMPILE__)
    _mode = cs_alloc_mode;
#else
    // use default and avoid compiler warnings
    _mode = CS_ALLOC_HOST_DEVICE_SHARED;
#endif
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  array
  (
    T              *data,      
    const cs_lnum_t(&dims)[N] 
  )
  :
    _span(data, dims),
    _owner(false)
  {
#if !defined(__CUDA_ARCH__) && \
    !defined(SYCL_LANGUAGE_VERSION) && \
    !defined(__HIP_DEVICE_COMPILE__)
    _mode = cs_alloc_mode;
#else
    // use default and avoid compiler warnings
    _mode = CS_ALLOC_HOST_DEVICE_SHARED;
#endif
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  array
  (
    const array& other 
  )
  :_span()
  {
    _span::set_size_(other._extent);
    _mode = other._mode;
    _owner = false;
    _span::set_data_ptr(other.data());
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  array
  (
    array&& other 
  )
  : array()
  {
    swap(*this, other);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  ~array()
  {
    clear();
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  friend void
  swap
  (
    array& first, 
    array& second 
  )
  {
    /* Swap the different members between the two references. */
    cs::swap_objects(first._span::_extent, second._span::_extent);
    cs::swap_objects(first._span::_offset, second._span::_offset);
    cs::swap_objects(first._span::_size, second._span::_size);
    cs::swap_objects(first._span::_data, second._span::_data);
    cs::swap_objects(first._owner, second._owner);
    cs::swap_objects(first._mode, second._mode);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  array& operator=(array other)
  {
    swap(*this, other);
 
    return *this;
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  void
  clear()
  {
    if (_owner) {
      // Avoid compiler warnings. If device, object cannot be owner...
#if   !defined(__CUDA_ARCH__) && \
      !defined(SYCL_LANGUAGE_VERSION) &&  \
      !defined(__HIP_DEVICE_COMPILE__)
      CS_FREE(_span::_data);
#endif
    }
    else {
      _span::set_data_ptr(nullptr);
    }
    _span::_size = 0;
    for (int i = 0; i < N; i++) {
      _span::_offset[i] = 0;
      _span::_extent[i] = 0;
    }
  }
 
  CS_F_HOST
  array<T, N, L>
  get_deep_copy
  (
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
  __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  ) const
  {
    array<T, N, L> new_array(this->_extent,
                             this->_mode,
                             file_name,
                             line_number);
    new_array.copy_data(*this);
 
    return new_array;
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  mdspan<T, N, L>
  view()
  {
    return mdspan<T,N,L>(_span::_data, _span::_extent);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  template<int _N_, layout _L_ = L>
  CS_F_HOST_DEVICE
  mdspan<T, _N_, _L_>
  get_mdspan
  (
    const cs_lnum_t(&dims)[_N_] 
  )
  {
    /* sanity check for total size */
    cs_lnum_t s = 1;
    for (int i = 0; i < _N_; i++)
      s *= dims[i];
 
    if (s != _span::_size) {
#if !defined(__CUDA_ARCH__) && \
    !defined(SYCL_LANGUAGE_VERSION) && \
    !defined(__HIP_DEVICE_COMPILE__)
      bft_error(__FILE__, __LINE__, 0,
                _("%s: requested span has total size of %d instead of %d "
                  "for this array.\n"),
                __func__, s, _span::_size);
#else
      return mdspan<T,_N_,_L_>();
#endif
    }
 
    return mdspan<T,_N_,_L_>(_span::_data, dims);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  template<layout _L_ = L, typename... Args>
  CS_F_HOST_DEVICE
  mdspan<T, sizeof...(Args), _L_>
  get_mdspan
  (
    Args... indices 
  )
  {
    check_non_empty_args_(indices...);
 
    constexpr int n_vals = sizeof...(Args);
    cs_lnum_t tmp_args[n_vals] = {indices...};
 
    return get_mdspan<n_vals, _L_>(tmp_args);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  void
  set_empty()
  {
    cs_lnum_t tmp[N] = {0};
    _span::set_size_(tmp);
    _span::set_data_ptr(nullptr);
    _owner = false;
#if !defined(__CUDA_ARCH__) && \
    !defined(SYCL_LANGUAGE_VERSION) && \
    !defined(__HIP_DEVICE_COMPILE__)
    _mode = cs_alloc_mode;
#else
    _mode = CS_ALLOC_HOST_DEVICE_SHARED;
#endif
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void reshape
  (
    cs_lnum_t    new_size,      
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
  __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  {
    static_assert(N == 1,
                  "Method reshape(size, ...) only possible for "
                  " cs::array<T> or cs::array<T,1>");
    assert(new_size >= 0);
 
    /* If same dimensions, nothing to do ... */
    if (new_size == _span::_size)
      return;
 
    if (_owner) {
      _span::set_size_(new_size);
      reallocate_(file_name, line_number);
    }
    else
      bft_error(__FILE__, __LINE__, 0,
                _("%s: array cannot be reshaped if non owner, "
                  "use mdspan instead.\n"), __func__);
 
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  reshape_and_copy
  (
    cs_lnum_t       new_size,     
    cs_lnum_t       size_to_keep, 
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
  __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  {
    static_assert(N==1,
                  "Method reshape_and_copy(size, ...) only possible for "
                  " cs::array<T> or cs::array<T,1>");
    assert(new_size >= 0);
    assert(size_to_keep <= new_size && size_to_keep <= _span::_size);
 
    /* If same dimensions, nothing to do ... */
    if (new_size == _span::_size)
      return;
 
    if (!(size_to_keep <= new_size && size_to_keep <= _span::_size))
      bft_error(__FILE__, __LINE__, 0,
                "%s: Data cannot be saved when new size is smaller "
                "than size to keep.\n",
                __func__);
 
    if (_owner) {
      /* If new size is larger, realloc is sufficient. */
      _span::set_size_(new_size);
      reallocate_(file_name, line_number);
    }
    else
      bft_error(__FILE__, __LINE__, 0,
                _("%s: array cannot be reshaped if non owner, "
                  "use mdspan instead.\n"), __func__);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void reshape
  (
    const cs_lnum_t(&dims)[N], 
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
  __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  {
    /* If same dimensions, nothing to do ... */
    bool same_dim = true;
    for (int i = 0; i < N; i++)
      if (dims[i] != _span::_extent[i])
        same_dim = false;
 
    if (same_dim)
      return;
 
    if (_owner) {
      _span::set_size_(dims);
      reallocate_(file_name, line_number);
    }
    else
      bft_error(__FILE__, __LINE__, 0,
                _("%s: array cannot be reshaped if non owner, "
                  "use mdspan instead.\n"), __func__);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void reshape_and_copy
  (
    const cs_lnum_t(&dims)[N],   
    cs_lnum_t size_to_keep = -1, 
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
  __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  {
    assert(size_to_keep <= dims[N-1] && size_to_keep <= _span::_extent[N-1]);
 
    /* If same dimensions, nothing to do ... */
    bool same_dim = true;
    for (int i = 0; i < N; i++)
      if (dims[i] != _span::_extent[i])
        same_dim = false;
 
    if (same_dim)
      return;
 
    if (_owner) {
      /* check that size_to_keep is either >0 or == -1 (keep all) */
      if (size_to_keep != 0) {
        /* Temporary copy */
        array<T,N,L> tmp = this->get_deep_copy();
 
        /* Update instance size */
        _span::set_size_(dims);
        reallocate_(file_name, line_number);
 
        /* Local copies for the context loop */
        cs_lnum_t new_o[N], new_e[N], old_o[N], old_e[N], max_e[N];
 
        for (int i = 0; i < N; i++) {
          new_o[i] = _span::offset(i);
          new_e[i] = _span::extent(i);
          old_o[i] = tmp.offset(i);
          old_e[i] = tmp.extent(i);
 
          max_e[i] = (new_e[i] > old_e[i]) ? old_e[i] : new_e[i];
        }
        if (max_e[N-1] > size_to_keep && size_to_keep > 0)
          max_e[N-1] = size_to_keep;
 
        T* old_data = tmp.data();
        T* new_data = _span::data();
 
        cs_lnum_t loop_size = tmp.size();
 
        /* Loop using dispatch */
        cs_dispatch_context ctx;
 
        ctx.parallel_for(loop_size, CS_LAMBDA (cs_lnum_t e_id) {
          cs_lnum_t idx[N];
          cs_lnum_t dummy = e_id;
 
          for (int i = N-1; i >= 1; i--) {
            idx[i] = dummy % old_e[i];
            dummy = (dummy - idx[i]) / old_e[i];
          }
          idx[0] = dummy;
 
          bool all_ok = true;
          for (int i = 0; i < N; i++) {
            if (idx[i] >= max_e[i]) {
              all_ok = false;
              break;
            }
          }
 
          if (all_ok) {
            cs_lnum_t o_id = 0;
            cs_lnum_t n_id = 0;
 
            for (int i = 0; i < N; i++) {
              o_id += idx[i] * old_o[i];
              n_id += idx[i] * new_o[i];
            }
 
            new_data[n_id] = old_data[o_id];
          }
        });
 
        ctx.wait();
      }
      else {
        reshape(dims, file_name, line_number);
      }
    }
    else
      bft_error(__FILE__, __LINE__, 0,
                _("%s: array cannot be reshaped if non owner, "
                  "use mdspan instead.\n"), __func__);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  reshape
  (
    const cs_lnum_t size1,       
    const cs_lnum_t size2,       
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
  __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  {
    static_assert(N == 2,
                  "Method reshape(size1, size2, ...) only possible for "
                  " cs::array<T,2>");
    cs_lnum_t tmp_size[N] = {size1, size2};
 
    reshape(tmp_size, file_name, line_number);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  reshape_and_copy
  (
    const cs_lnum_t size1,  
    const cs_lnum_t size2,  
    cs_lnum_t size_to_keep, 
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
  __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  {
    static_assert(N == 2,
                  "Method reshape_and_copy(size1, size2, ...) only possible for "
                  " cs::array<T,2>");
    cs_lnum_t tmp_size[N] = {size1, size2};
 
    reshape_and_copy(tmp_size, size_to_keep, file_name, line_number);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  reshape
  (
    const cs_lnum_t size1,       
    const cs_lnum_t size2,       
    const cs_lnum_t size3,       
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
  __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  {
    static_assert(N == 3,
                  "Method reshape(size1, size2, size3, ...) only possible for "
                  " cs::array<T,3>");
    cs_lnum_t tmp_size[N] = {size1, size2, size3};
 
    reshape(tmp_size, file_name, line_number);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  reshape_and_copy
  (
    const cs_lnum_t size1,  
    const cs_lnum_t size2,  
    const cs_lnum_t size3,  
    cs_lnum_t size_to_keep, 
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
  __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  {
    static_assert(N == 3,
                  "Method reshape_and_copy(size1, size2, size3, ...) only "
                  "possible for cs::array<T,3>");
    cs_lnum_t tmp_size[N] = {size1, size2, size3};
 
    reshape_and_copy(tmp_size, size_to_keep, file_name, line_number);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  reshape_and_copy
  (
    const cs_lnum_t size1,  
    const cs_lnum_t size2,  
    const cs_lnum_t size3,  
    const cs_lnum_t size4,  
    cs_lnum_t size_to_keep, 
#if (defined(__GNUC__) || defined(__clang__)) && \
  __has_builtin(__builtin_LINE) && \
  __has_builtin(__builtin_FILE)
    const char *file_name   = __builtin_FILE(), 
    const int   line_number = __builtin_LINE()  
#else
    const char *file_name   = __FILE__, 
    const int   line_number = __LINE__  
#endif
  )
  {
    static_assert(N == 4,
                  "Method reshape_and_copy(size1, size2, size3, size4, ...)"
                  " only possible for cs::array<T,4>");
    cs_lnum_t tmp_size[N] = {size1, size2, size3, size4};
 
    reshape_and_copy(tmp_size, size_to_keep, file_name, line_number);
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  void set_alloc_mode
  (
    cs_alloc_mode_t mode 
  )
  {
    if (_mode != mode) {
      /* Since allocation mode is changed, deallocate data */
      if (_span::_size != 0)
        clear();
 
      _mode = mode;
    }
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  inline
  bool
  owner() const
  {
    return _owner;
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST_DEVICE
  inline
  cs_alloc_mode_t
  mode() const
  {
    return _mode;
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  inline
  void
  copy_data
  (
    const array&    other,      
    const cs_lnum_t n_vals = -1 
  )
  {
    auto& ctx = cs::execution::default_context();
    copy_data(ctx, other, n_vals);
    ctx.wait();
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  copy_data
  (
    cs_dispatch_context &ctx,        
    const array         &other,      
    const cs_lnum_t      n_vals = -1 
  )
  {
    const cs_lnum_t loop_size = (n_vals == -1) ? _span::_size : n_vals;
 
    assert(loop_size <= _span::size());
    assert(loop_size <= other.size());
 
    // Explicit pointer, avoid passing internal member of class to the functor
    T* data_ptr = _span::data();
    T* o_data_ptr = other.data();
 
    ctx.parallel_for(loop_size, CS_LAMBDA (cs_lnum_t e_id) {
      data_ptr[e_id] = o_data_ptr[e_id];
    });
  }
 
private:
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  template<typename... Args>
  CS_F_HOST_DEVICE
  static inline
  void
  check_non_empty_args_
  (
    Args... 
  )
  {
    static_assert(sizeof...(Args) > 0, "At least one input argument is needed.");
    static_assert(cs::are_integral<Args...>::value, "Non integral input arguments.");
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  allocate_
  (
    const char *file_name,  
    const int   line_number 
  )
  {
    const char *_ptr_name = "cs::array._data";
    _span::_data = static_cast<T *>(cs_mem_malloc_hd(_mode,
                                                     _span::_size,
                                                     sizeof(T),
                                                     _ptr_name,
                                                     file_name,
                                                     line_number));
  }
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
 
  CS_F_HOST
  void
  reallocate_
  (
    const char *file_name,  
    const int   line_number 
  )
  {
    /* If not owner no-op */
    if (_owner) {
      const char *_ptr_name = "cs::array._data";
      _span::_data = static_cast<T *>(cs_mem_realloc_hd(_span::_data,
                                                        _mode,
                                                        _span::_size,
                                                        sizeof(T),
                                                        _ptr_name,
                                                        file_name,
                                                        line_number));
    }
  };
 
  /*--------------------------------------------------------------------------*/
  /*--------------------------------------------------------------------------*/
  CS_F_HOST
  inline
  void
  set_to_nan_()
  {
    auto& ctx = cs::execution::default_context();
    constexpr T tmp_nan = std::numeric_limits<T>::signaling_NaN();
    _span::set_to_val(ctx, tmp_nan);
    ctx.wait();
  }
  /*===========================================================================
   * Private members
   *==========================================================================*/
 
  bool            _owner;
  cs_alloc_mode_t _mode;
 
};
 
} /* namespace cs */
 
/*--------------------------------------------------------------------------*/
/* Useful aliases without namespace */
/*--------------------------------------------------------------------------*/
 
template<class T>
using cs_array = cs::array<T,1,cs::layout::right>;
 
template<class T, cs::layout L = cs::layout::right>
using cs_array_2d = cs::array<T,2,L>;
 
template<class T, cs::layout L = cs::layout::right>
using cs_array_3d = cs::array<T,3,L>;
 
template<class T, cs::layout L = cs::layout::right>
using cs_array_4d = cs::array<T,4,L>;
 
template<class T, int N, cs::layout L = cs::layout::right>
using cs_mdarray = cs::array<T,N,L>;
 
template<class T, int N>
using cs_mdarray_r = cs::array<T,N,cs::layout::right>;
 
template<class T, int N>
using cs_mdarray_l = cs::array<T,N,cs::layout::left>;
 
/*----------------------------------------------------------------------------*/
 
#endif /* CS_ARRAY_H */