The tensor ops functions typically require a single/combination of integer arguments to set the bounds for the function. Take for example:

template< std::ptrdiff_t ISIZE, typename VECTOR >
LVARRAY_HOST_DEVICE CONSTEXPR_WITHOUT_BOUNDS_CHECK inline
auto l2NormSquared( VECTOR const & vector )
{
  static_assert( ISIZE > 0, "ISIZE must be greater than zero." );
  internal::checkSizes< ISIZE >( vector );

  auto norm = vector[ 0 ] * vector[ 0 ];
  for( std::ptrdiff_t i = 1; i < ISIZE; ++i )
  {
    norm = norm + vector[ i ] * vector[ i ];
  }
  return norm;
}

The call looks something like:

l2NormSquared<3>(array)

where the N(3) is always required, and the type VECTOR is deduced. This is a little bit clunky imo. We could do something where we deduced both a size and type. Like so:

#include <tuple>
#include <iostream>
#include <type_traits>

template< typename _T > 
  struct HasMemberFunction_size_impl 
  { 
private: 
    template< typename CLASS > static constexpr auto test( int )->decltype( std::is_convertible< decltype( std::declval< CLASS >().size(  ) ), int >::value, int() ) 
    { return CLASS::size(); } 
    template< typename CLASS > static constexpr auto test( ... )->int 
    { return sizeof(CLASS)/sizeof(std::declval<CLASS>()[0]);; } 
public: 
    static constexpr int value = test< _T >( 0 );
  }; 
template< typename CLASS > 
static constexpr int HasMemberFunction_size = HasMemberFunction_size_impl< CLASS >::value;




template< int N >
struct Tensor
{
    static constexpr int size() 
    { return N; }

    double & operator[]( int const i )
    {
        return data[i];
    }

    double const & operator[]( int const i ) const
    {
        return data[i];
    }

    double data[N];
};



template< typename T, int N>
double func_impl( T const & array )
{
    double rvalue = array[0]*array[0];
    for( int i=1 ; i<N ; ++i )
    {
        rvalue += array[i]*array[i];
    }
    return rvalue;
}

template< typename T, int N = HasMemberFunction_size<T> >
typename std::enable_if<!std::is_pointer<T>::value, double>::type func( T const & array )
{
   std::cout<<"Calling reference version"<<std::endl;
   return func_impl<T,N>(array);
}

template< int N, typename T >
typename std::enable_if<std::is_pointer<T>::value, double>::type func( T const array )
{
   std::cout<<"Calling pointer version"<<std::endl;
   return func_impl<T,N>(array);
}



////////////////////////////////////////////

int main()
{
    Tensor<3> t;

    t[0] = 1;
    t[1] = 2;
    t[2] = 3;


    double a[3] = {1,2,3};
    double * const pa = a;


    std::cout<<func(t)<<std::endl;
    std::cout<<func(a)<<std::endl;
    std::cout<<func<3>(pa)<<std::endl;

}

Here it is in a compiler explorer:
https://godbolt.org/z/K1P4T9qdE

If we were to standardize the way we return compile time sizes, this works fairly well. With a raw pointer you would still require the specification of N. With multi-dimensions this would be more complicated...but I don't think it is too prohibitive, and the usability is certainly nicer.

@corbett5 @klevzoff Thoughts?

I get this on Ascent with gcc-8.1.1
Is this a legit error or gcc acting weird?

In file included from /autofs/nccsopen-svm1_home/klevtsov/GEOSX/src/coreComponents/cxx-utilities/src/src/ArrayView.hpp:29,
                 from /autofs/nccsopen-svm1_home/klevtsov/GEOSX/src/coreComponents/cxx-utilities/src/src/Array.hpp:30,
                 from /autofs/nccsopen-svm1_home/klevtsov/GEOSX/src/coreComponents/common/DataTypes.hpp:40,
                 from /autofs/nccsopen-svm1_home/klevtsov/GEOSX/src/coreComponents/dataRepository/xmlWrapper.hpp:27,
                 from /autofs/nccsopen-svm1_home/klevtsov/GEOSX/src/coreComponents/dataRepository/Group.hpp:27,
                 from /autofs/nccsopen-svm1_home/klevtsov/GEOSX/src/coreComponents/dataRepository/Group.cpp:15:
/autofs/nccsopen-svm1_home/klevtsov/GEOSX/src/coreComponents/cxx-utilities/src/src/ChaiVector.hpp: In instantiation of 'void LvArray::ChaiVector<T>::move(chai::ExecutionSpace, bool) [with T = std::__cxx11::basic_string<char>]':
/autofs/nccsopen-svm1_home/klevtsov/GEOSX/src/coreComponents/cxx-utilities/src/src/Array.hpp:426:5:   required from 'typename std::enable_if<(! LvArray::detail::is_array<U>::value), void>::type LvArray::Array<T, NDIM, INDEX_TYPE, DATA_VECTOR_TYPE>::move(chai::ExecutionSpace, bool) [with U = std::__cxx11::basic_string<char>; T = std::__cxx11::basic_string<char>; int NDIM = 1; INDEX_TYPE = long int; DATA_VECTOR_TYPE = LvArray::ChaiVector<std::__cxx11::basic_string<char> >; typename std::enable_if<(! LvArray::detail::is_array<U>::value), void>::type = void]'
/autofs/nccsopen-svm1_home/klevtsov/GEOSX/src/coreComponents/dataRepository/Wrapper.hpp:991:7:   required from 'static typename std::enable_if<geosx::dataRepository::Wrapper<U>::move_wrapper::has_memberfunction_move<TT>::value, void>::type geosx::dataRepository::Wrapper<U>::move_wrapper::move(U&, chai::ExecutionSpace, bool) [with U = LvArray::Array<std::__cxx11::basic_string<char>, 1, long int, LvArray::ChaiVector<std::__cxx11::basic_string<char> > >; T = LvArray::Array<std::__cxx11::basic_string<char>, 1, long int, LvArray::ChaiVector<std::__cxx11::basic_string<char> > >; typename std::enable_if<geosx::dataRepository::Wrapper<U>::move_wrapper::has_memberfunction_move<TT>::value, void>::type = void]'
/autofs/nccsopen-svm1_home/klevtsov/GEOSX/src/coreComponents/dataRepository/Wrapper.hpp:1001:30:   required from 'void geosx::dataRepository::Wrapper<U>::move(chai::ExecutionSpace, bool) [with T = LvArray::Array<std::__cxx11::basic_string<char>, 1, long int, LvArray::ChaiVector<std::__cxx11::basic_string<char> > >]'
/autofs/nccsopen-svm1_home/klevtsov/GEOSX/src/coreComponents/dataRepository/Wrapper.hpp:1000:16:   required from here
/autofs/nccsopen-svm1_home/klevtsov/GEOSX/src/coreComponents/cxx-utilities/src/src/ChaiVector.hpp:473:63: error: dereferencing type-punned pointer will break strict-aliasing rules [-Werror=strict-aliasing]
     else reinterpret_cast< chai::ManagedArray< T const > & >( m_array ).move( space );
                                                               ^~~~~~~

Most of the LvArray tests use TYPED_TEST but not all. Fix this. Also remove dependence on chai::forall and instead use RAJA directly, then refactor CUDA test cases to be templated on the execution policy and do, serial, host parallel, and device parallel. See testCRSMatrix for an example.

I could be wrong but I don't think anything in either SciPy or numpy has something similar to the ArrayOfArrays. Therefore I think it is best to create a python class implemented in C which will do the trick. This class needs needs to support

• Returning the number of inner arrays.
• Returning a 1D numpy array representing an inner array.
• Returning a specific value.
• A python iterator interface?

To create one object that can represent any combination of T and INDEX_TYPE might be tricky. One approach is to have the object contain three ndarrays representing the values, sizes and offsets. Then the type agnostic ndarray c-api can be used to access the entries and create slices. I think that would work but if we can only easily get away with having one object per INDEX_TYPE that'd be fine too.

Add tests to cover the contents of #109

LvArray has been released with BSD3. We need to change the license file, and the header for each file.

In the conversion to const user defined conversions we have code like this:

  operator Array<T const, NDIM, INDEX_TYPE> const & () const
  {
    return reinterpret_cast<Array<T const, NDIM, INDEX_TYPE> const &>(*this);
  }

While this is great from the brevity point of view, it may result in undefined behavior. We should keep this in mind, and perhaps at some point consider just shallow constructing new objects for as part of this type of user defined conversion.

We can improve the pylvarray test coverage with pure C++ tests as follows

template< typename T >
void testSortedArray()
{
  SortedArray< T > set;
  populate( set);
  SortedArray< T > copy( arr );
  PyObjectRef<> pySet = create( arr );
  mutate( pySet, copy );
  check( copy, arr, pySet);
}

instead of

void append( T * const ptr, INDEX_TYPE const size, T const * const values, INDEX_TYPE const n )

make it

void append( T * const ptr, INDEX_TYPE const size, InputIt fist, INDEX_TYPE const n )

and likewise for similar functions.

Also create a portable wrapper for std::move_iterator and use it here
https://github.com/GEOSX/cxx-utilities/blob/ac9b2469dccefc19e4adf5810944a8a050f12967/src/src/sortedArrayManipulation.hpp#L715

to move insert the values in buffer.

Without pragmas to suppress -Wstrict-aliasing around the reinterpret cast in Array<T>::move() (when T is an array) the following error occurs when building with gcc in release:

In file included from /usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/common/DataTypes.hpp:30,
                 from /usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/dataRepository/BufferOps.hpp:20,
                 from /usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/dataRepository/WrapperHelpers.hpp:28,
                 from /usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/dataRepository/Wrapper.hpp:23,
                 from /usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/dataRepository/Group.hpp:27,
                 from /usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/managers/ObjectManagerBase.hpp:22,
                 from /usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/mesh/ElementSubRegionBase.hpp:22,
                 from /usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/mesh/FaceElementSubRegion.hpp:22,
                 from /usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/mesh/FaceElementSubRegion.cpp:19:
/usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/cxx-utilities/src/Array.hpp: In instantiation of 'std::enable_if_t<isArray<U> > LvArray::Array<T, NDIM, PERMUTATION, INDEX_TYPE, DATA_VECTOR_TYPE>::move(chai::ExecutionSpace, bool) [with U = LvArray::Array<long int, 1, camp::int_seq<long int, 0>, long int, LvArray::ChaiBuffer>; T = LvArray::Array<long int, 1, camp::int_seq<long int, 0>, long int, LvArray::ChaiBuffer>; int NDIM = 1; PERMUTATION = camp::int_seq<long int, 0>; INDEX_TYPE = long int; BUFFER_TYPE = LvArray::ChaiBuffer; std::enable_if_t<isArray<U> > = void]':
/usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/dataRepository/Wrapper.hpp:531:7:   required from 'static typename std::enable_if<has_chai_move_method<U>, void>::type geosx::dataRepository::Wrapper<T>::move_wrapper::move(U&, chai::ExecutionSpace, bool) [with U = geosx::InterObjectRelation<LvArray::Array<LvArray::Array<long int, 1, camp::int_seq<long int, 0>, long int, LvArray::ChaiBuffer>, 1, camp::int_seq<long int, 0>, long int, LvArray::ChaiBuffer> >; T = geosx::InterObjectRelation<LvArray::Array<LvArray::Array<long int, 1, camp::int_seq<long int, 0>, long int, LvArray::ChaiBuffer>, 1, camp::int_seq<long int, 0>, long int, LvArray::ChaiBuffer> >; typename std::enable_if<has_chai_move_method<U>, void>::type = void]'
/usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/dataRepository/Wrapper.hpp:541:30:   required from 'void geosx::dataRepository::Wrapper<T>::move(chai::ExecutionSpace, bool) [with T = geosx::InterObjectRelation<LvArray::Array<LvArray::Array<long int, 1, camp::int_seq<long int, 0>, long int, LvArray::ChaiBuffer>, 1, camp::int_seq<long int, 0>, long int, LvArray::ChaiBuffer> >]'
/usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/dataRepository/Wrapper.hpp:540:16:   required from here
/usr/WS1/tobin6/geosx/GEOSX/src/coreComponents/cxx-utilities/src/Array.hpp:545:5: error: dereferencing type-punned pointer will break strict-aliasing rules [-Werror=strict-aliasing]
     reinterpret_cast< BUFFER_TYPE< typename T::ViewType > & >( m_dataBuffer ).move( space, touch );

SortedArrayView and by extension SortedArray should be able to create an ArraySlice< T const, 1, 0 >.

Make something similar to the RAJA MultiView. Here's what I imagine it would look like.

template< int NARRAYS, typename T, int NDIM, int USD, typename INDEX_TYPE, template< typename > class BUFFER_TYPE >
class MultiView
{
public:
  MultiView( CArray< ArrayView< T, NDIM, USD, INDEX_TYPE, BUFFER_TYPE > > && views ):
    m_dims{ views[ 0 ].m_dims },
    m_strides{ views[ 0 ].m_strides }
  {
    // Error checking

    for( int i = 0; i < NARRAYS; ++i )
    {
      m_dataBuffers[ i ] = std::move( views[ i ].m_dataBuffer );
    }
  }

  inline LVARRAY_HOST_DEVICE constexpr
  ArraySlice< T, NDIM, USD, INDEX_TYPE >
  toSlice( int const i ) const & noexcept
  { return ArraySlice< T, NDIM, USD, INDEX_TYPE >( m_dataBuffers[ i ].data(), m_dims.data, m_strides.data ); }

private:
  CArray< INDEX_TYPE, NDIM > m_dims;
  CArray< INDEX_TYPE, NDIM > m_strides;
  CArray< BUFFER_TYPE< T >, NARRAYS > m_dataBuffers;
};

Would need to be a bit more complicated to ensure the same copy-semantics as ArrayView.

With a little work it could be extended to support views with different types, which would be necessary if you want to mix T with T const.

Usage would look something like

void foo( ArrayView2d< int > const & xIn, ArrayView2d< int > const & yIn )
{
  MultiView2d< 2, int > const multiView( { xIn, yIn } );

  forAll< POLICY >( xIn.size( 0 ), [multiView] ( int const i )
  {
    ArraySlice2d< int > const x = multiView.toSlice( 0 );
    ArraySlice2d< int > const y = multiView.toSlice( 1 );

    bar( x, y );
  } );
}

Definitely since it can avoid a move. I think the best way to do this is have registerTouch allocate in the given space if it didn't already exist and then touch it there. That way all memcpy has to do is replace move( space, true ) with registerTouch( space ). I don't really want to add that to this PR though so I'll create an issue.

Originally posted by @corbett5 in #231 (comment)

in addition to these stream operators, it would be nice to have an output option that would dump using printf so that we could use it on device. For instance a printRow( INDEX_TYPE const row ).

Originally posted by @rrsettgast in #164

Since they are such an important backbone of the geosx infrastructure, it would be nice to have clear developer guidelines on proper use of array/view/slice classes. Example of topics to cover:

• when to use each, e.g. when accessing field data should I use array<T,N> &, array_view<T,N> & or array_view<T,N> to capture the value returned by getReference<array<T,N>>()
• which of the classes are safe to capture by value/reference in lambdas passed to RAJA
• which operations involve shallow/deep copies of data

A lot of it can be inferred by reading the source, but some of it requires going deep, e.g. knowing that ChaiVector's CC does shallow copies by default (and so does ArrayView). Or for example it is not clear why

array_view<T,N> v = obj->getReference<array<T,N>>(...)

is explicitly prohibited, but

array_view<T,N> v;
v = obj->getReference<array<T,N>>(...)

is allowed. Or why there is no assignment operator of any kind in ArraySlice (or size() for that matter).

In the CRSMatrix::addToRow methods there are multiple uses of LVARRAY_ASSERT. These should probably be LVARRAY_ERROR and only enabled if ARRAY_BOUNDS_CHECK is on. Same thing in other places.

@andrembcosta needs the spectral decompositions that used to live in R2SymTensor before it was removed from the code for phase field fracture.

Currently sorting::introsortLoop uses recursion, nvcc allocates space for the stack at compile time and it issues a warning that it can't statically determine the stack size. This isn't an pressing issue but it does bring up the possibility that we blow through the stack (unlikely though since it has a controlled recursion level). However if it can be rewritten without recursion and without sacrificing performance this would be nice.

The current tensorOps tests are to their benefit

• Exhaustive.

and to their detriment

• Complicated. They involve heavy usage of macros and are spread across multiple source files.
• Slow to compile, maybe a minute or two with NVCC and super long with Intel or XL (not that I really care about those).

This is proving to be the part of LvArray that is touched most frequently by other developers and as such I would like it to be easy to add tests. If there is an existing test that is similar to the new function copying works fairly well, but due the unique file structure it can be hard to find a specific tests, and if there isn't a similar test writing one from scratch is non-trivial.

Currently much of the complexity arises from my desire to test all of the template functions with many of the possible argument permutations. For example to test the trivial copy I test all nine permutations of ArraySlice< T, 1, 0 >, ArraySlice< T, 1 -1 > and T[ N ]. This test is then instantiated for three different combinations of T and N, and another 3 when using CUDA. Just looking at the code for copy it is pretty obvious that it is correct so this seems like overkill, the purpose of the tests is as much to verify that things compile as to verify runtime correctness. The template explosion leads to long compile times and I've spread the tests across multiple translation units to make it faster.

I'm starting to think that it might be worth decreasing the test coverage if it could significantly decrease the barriers to entry. What has been your experience adding tests, was there anything in particular that took some work to figure out?

Tagging people who have added tests:
@klevzoff @francoishamon @CusiniM @rrsettgast

pretty self explanatory.

Currently we can convert an ArrayView<T> object to ArrayView<const T> const &, which comes in handy when you need to use it in a non-modifying kernel, since const guarantees a touch is not registered in the target memory space, which helps avoid copying arrays back and forth.

It would also be useful to have similar conversion for nested arrays, i.e. ArrayView<ArrayView<T>> to
ArrayView<ArrayView<const T> const> const & (not sure if the middle const is needed).

This may be needed in FV kernels in GEOSX which don't loop over the subregions (e.g. flux computations), as they need to capture ElementViewAccessor objects, that are essentially Array<Array<ArrayView<T>>>, and we would like to avoid touching the innermost array.

Create something like this

namespace LvArray
{
namespace python
{
namespace internal
{

/**
 * @brief Create a numpy ndarray representing the given array. If T is const then the array should be immutable.
 */
template< typename T, typename INDEX_TYPE >
PyObject * create( T * const data, int ndim, INDEX_TYPE const * const dims, INDEX_TYPE const * const strides )
{
  ...
}

} // namespace internal

template< typename T, int NDIM, ... >
PyObject * create( ArraySlice< T, NDIM, ... > const & slice )
{ return internal::create( slice.data(), NDIM, slice.dims(), slice.strides() ); }

template< ... >
PyObject * create( ArrayView< ... > const & view )
{
  LVARRAY_ERROR_IF_NE( view.currentSpace(), MemorySpace::CPU );
  return create( view.toSlice() );
}

template< typename T, typename INDEX_TYPE, ... >
PyObject * create( SortedArrayView< T const, INDEX_TYPE, .... > const & view )
{
  INDEX_TYPE const dims = view.size();
  INDEX_TYPE const strides = 1;
  return internal::create( view.data(), 1, dims, strides );
}

} // namespace python
} // namespace LvArray

So that a user can use these methods in their own python module.

You should not be able to create an ArrayView or an ArraySlice from a Array && and you should not be able to create an ArraySlice from an ArrayView &&. Similarly for the other classes. This includes the creation of ArraySlice from the operator[] but not from CRSMatrixView::getEntries.

At the same time it would be nice to see if there's an easy way to disable the conversion from Array to ArrayView &. If this involves making ArrayView a protected base class then a separate issue should be made.

The moved-from object will segfault on destruction with the following stack trace:

LvArray::ArrayOfArraysView<long long, long, false>::destroyValues<>(long, long)::{lambda(auto:1&)#1}::operator()<LvArray::ChaiBuffer<long long> >(LvArray::ChaiBuffer<long long>&) const ArrayOfArraysView.hpp:828
for_each_arg<LvArray::ArrayOfArraysView<long long, long, false>::destroyValues<>(long, long)::{lambda(auto:1&)#1}, LvArray::ChaiBuffer<long long>&>(LvArray::ArrayOfArraysView<long long, long, false>::destroyValues<>(long, long)::{lambda(auto:1&)#1}&&, LvArray::ChaiBuffer<long long>&) templateHelpers.hpp:46
LvArray::ArrayOfArraysView<long long, long, false>::destroyValues<>(long, long) ArrayOfArraysView.hpp:823
LvArray::ArrayOfArraysView<long long, long, false>::free<>() ArrayOfArraysView.hpp:536
LvArray::SparsityPattern<long long, long>::~SparsityPattern SparsityPattern.hpp:93

Basically, the class invariant is broken after move, because m_numArrays is not reset to 0 and is no longer consistent with m_offsets and m_sizes (which are empty and destroyValues is trying to access them). Having to reset m_numArrays during move to ensure safe destruction, unfortunately, will prevent us from defaulting move constructor/assignment. Maybe there's an easier way?

The call to backtrace in stackTrace fails if used after an OpenMP parallel region but only with clang10 debug on Quartz. Strange.

This will help prevent combinatorial explosion to cover all the permutation types.

template <class T, class INDEX_TYPE, bool isSorted>
class FooView
{
  ...
  ChaiVector<INDEX_TYPE const> m_offsets;
  ChaiVector<INDEX_TYPE> m_sizes;
  ChaiVector<T> m_values;
};
  
template <class T, class INDEX_TYPE, bool isSorted>
class Foo : public FooView<T, INDEX_TYPE, isSorted>
{
  ...
};

Then we can do

template <class T, class INDEX_TYPE>
using CSArray2DView = FooView<T, INDEX_TYPE, false>;

template <class T, class INDEX_TYPE>
using CSArray2D = Foo<T, INDEX_TYPE, false>;

template <class T, class INDEX_TYPE>
using CSArrayOfSetsView = FooView<T, INDEX_TYPE, true>;

template <class T, class INDEX_TYPE>
using CSArrayOfSets = Foo<T, INDEX_TYPE, true>;

template <class T, class INDEX_TYPE>
SparsityPatternView : protected CSArrayOfSets<T, INDEX_TYPE>
{
 ...
};

We have an ArrayOfSetsView but no ArrayOfSets. Should be very similar to SparsityPattern.

https://github.com/LLNL/camp/blob/master/include/camp/resource/platform.hpp
https://github.com/LLNL/RAJA/blob/develop/include/RAJA/internal/get_platform.hpp
https://github.com/LLNL/RAJA/blob/develop/include/RAJA/util/camp_aliases.hpp

I am having compilation errors on Mac OSX with Xcode, Apple clang version 11.0.0 (clang-1100.0.33.8 in the following files:

src/src/ArrayView.hpp, line 284: Variable of non-literal type 'std::ostringstream' (aka 'basic_ostringstream<char>') cannot be defined in a constexpr function
in GEOS_ASSERT_MSG

Same error in:

• src/src/ArrayOfArraysView.hpp, line 198 in GEOS_ASSERT
• src/src/ArraySlice.hpp, line 202 in GEOS_ASSERT_MSG

It builds OK if I comment out the asserts.
Looks like it may be a question of standard.

@francoishamon saw some striking differences when he switched to the RAJA view, I'd like to add it as a benchmark here.
See LLNL/RAJA#1002

@corbett5 I started to include the call to move for the cases we discussed (resize, reserve), but this is actually a larger issue which we should have a unified approach for. It is really a situation where upon every type of array size/capacity modification that must be done on host, we should move back to host prior to performing the operation. This should be done across the board for all types (Array, ArrayOfArrays, CRSMatrix, SortedArray...).

Pros/Cons:

1. Pro. The developer can be assured that they will not perform even if operations are not called in a host policy raja loop. There is a lot of benefit to this.
2. Con. There are a lot of places that we will call move. It seems like we may call move multiple times with the way that the manipulation functions are setup. I don't know if this is a real issue, but I guess move is essentially a no-op if the buffer is already current on host? I think it is fine, but perhaps it isn't.
3. Is this always what the developer would want? Are there cases where this would be undesirable?

Split Macros.hpp up into two files, one with just the logging macros. That way you don't need to #include <iostream> just to get LVARRAY_HOST_DEVICE.

figure out what the best INDEX_TYPE for arrays is on each platform. I notice all the index types I see in cuda examples are "int". Does that mean that we should be using "int" indices when running cuda kernels? Perhaps this should be a platform dependent type?

When building standalone we should build with the same flags as GEOSX. For example we're not building with -Wall or -Werror.

One option is to have geosxOptions.cmake include cxx-utilitiesOptions.cmake.

Solution: Change the int in the the Array template parameter list to an std::integral_constant like object. Or just leave it. I don't believe there's a way to fix is_instantiation_of.

Also examine using std::fma.

The ArrayOfArrays family have a method resizeFromCapacities( INDEX_TYPE, INDEX_TYPE const * ) this should be rewritten as resizeFromCapacities( Array< T, 1, RAJA::PERM_I, INDEX_TYPE, BUFFER_TYPE > && ) where the ArrayOfArrays steals the buffer from the Array.

We need to get the TPL dependencies built as part of this repository for stand-alone builds. I think this is just:

• CHAI/Umpire
• RAJA/Camp

anything else?

These comments/issues were provided by steve glenn after some effort to integrate LvArray into Tribol.

1. git/config contains references to ssh://github.com/LLNL instead of https://github.com/LLNL, which causes failures if ssh is not configured.

• Fixed

2. Docs fail to build (hang) on Lassen.

• Fixed

3. Doc build cannot be disabled using ENABLE_DOCS=Off

• Fixed

4. Library and header exports are missing from install directory, so clients are forced to pull from the build directory. The install step emits a “Cannot create regular file” error.

• Fixed

5. The build should complete for all four combinations of ENABLE_CUDA=On/Off and ENABLE_OPENMP=On/Off.

• Fixed

6. LvArray cannot be built with xlc V16.1.1 on Lassen. According to LC documentation, https://hpc.llnl.gov/training/tutorials/using-lcs-sierra-system#Compilers, “… recommended and supported compilers are those delivered from IBM (XL and Clang ) and NVIDIA (NVCC).” A fix from IBM may also be needed for xlc to work with LvArray.

• Fixed

7. LvArray should be buildable against independent builds of recent versions of Axom, CHAI/Umpire and RAJA, not only those in /usr/gapps/GEOSX/thirdPartyLibs. This is important since applications may need to customize various options.
   For example, Tribol needs Slic, Spin, Quest, and Primal from Axom, but /usr/gapps/GEOSX/thirdPartyLibs/2020-05-23/install-lassen-clang@upstream-release only contains slic and lumberjack.

• Fixed. Hard dependency on RAJA, soft dependency on Chai and Caliper. No dependency on Axom.

8. The documentation needs better examples demonstrating use of two-dimensional Arrays, ArraySlices, and ArrayOfArrays. Currently, the unit tests are the only available examples.
   #159

9. Consider adding support for emplace_back(). In addition to making LvArray conform more closely to the std::vector interface, some applications may want this instead of push_back() for better performance.

• Fixed.

For exporting a CRSMatrix I think we want to use the scipy.sparse.csr_matrix object. This would involve verifying that the matrix is compressed, then creating three 1D ndarrays and calling the python method to construct the csr_matrix object.

There doesn't seem to be something equivalent to the SparsityPattern so we could either export it as an ArrayOfSets or create a ndarray of true representing the values and then create a crs_matrix object.

https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html#scipy.sparse.csr_matrix
https://stackoverflow.com/questions/3286448/calling-a-python-method-from-c-c-and-extracting-its-return-value

Describe the bug
Funny things happen when cxx_utilities::stringToArray is called with a string that has missing commas, e.g. "{1 2 3}". The first pass over the input string detects dimension 0 (because no commas are encountered) and resizes the array to 0. The second pass then reads as many values as the string contains into unallocated memory and corrupts the heap. The function completes fine, but some time later another heap allocation throws a glibc error and completely hangs the program (has to be kill -9ed, maybe it's doing a core dump, but I didn't wait long enough). Perhaps different systems will behave differently. Adding commas fixes the problem.

To Reproduce
Steps to reproduce the behavior:

1. Run GEOSX with an input file that has an array input without commas.

Expected behavior
Some sort of user-friendly error. Perhaps the second pass should check that the total number of values read does not exceed size() of the array.

Platform (please complete the following information):

• Ubuntu 19.04
• gcc-8.3
• debug or release mode

Additional context
Output from geosx before it hangs (last line is a glibc error):

real64 is alias of double
localIndex is alias of long
globalIndex is alias of long long
== CALIPER: (0): Initialized
== CALIPER: (0): Config check: No services enabled. Caliper will not record data.
GEOS must be configured to use Python to use parameters, symbolic math, etc. in input files
Adding Solver of type LaplaceFEM, named laplace
Adding Mesh: InternalMesh, mesh1
Adding Event: PeriodicEvent, solverApplications
Adding Event: PeriodicEvent, outputs
Adding Event: PeriodicEvent, restarts
   TableFunction: timeFunction
Adding Output: Silo, siloOutput
corrupted size vs. prev_size

GDB stack trace after error:

#0  __GI_raise (sig=sig@entry=6) at ../sysdeps/unix/sysv/linux/raise.c:50
#1  0x00007fffee7f4535 in __GI_abort () at abort.c:79
#2  0x00007fffee85b726 in __libc_message (action=action@entry=do_abort, fmt=fmt@entry=0x7fffee981952 "%s\n") at ../sysdeps/posix/libc_fatal.c:181
#3  0x00007fffee86259a in malloc_printerr (str=str@entry=0x7fffee97f9b0 "corrupted size vs. prev_size") at malloc.c:5352
#4  0x00007fffee862841 in unlink_chunk (p=<optimized out>, av=0x7fffee9b3c40 <main_arena>) at malloc.c:1464
#5  0x00007fffee863a18 in malloc_consolidate (av=av@entry=0x7fffee9b3c40 <main_arena>) at malloc.c:4510
#6  0x00007fffee865ca0 in _int_malloc (av=av@entry=0x7fffee9b3c40 <main_arena>, bytes=bytes@entry=1280) at malloc.c:3711
#7  0x00007fffee867caa in __GI___libc_malloc (bytes=1280) at malloc.c:3073
#8  0x00007ffff7a405f8 in google::dense_hashtable<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const, long>, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::hash<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, google::dense_hash_map<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, long, std::hash<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::equal_to<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, google::libc_allocator_with_realloc<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const, long> > >::SelectKey, google::dense_hash_map<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, long, std::hash<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::equal_to<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, google::libc_allocator_with_realloc<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const, long> > >::SetKey, std::equal_to<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, google::libc_allocator_with_realloc<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const, long> > >::set_empty_key(std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const, long> const&) ()
   from /home/klevtsov/work/GEOSX/build-clang8-debug/lib/libgeosx_core.so
#9  0x00007ffff7a404f4 in google::dense_hash_map<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, long, std::hash<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, std::equal_to<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > >, google::libc_allocator_with_realloc<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const, long> > >::set_empty_key(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&) ()
   from /home/klevtsov/work/GEOSX/build-clang8-debug/lib/libgeosx_core.so
#10 0x00007ffff7a3fa0b in axom::sidre::MapCollection<axom::sidre::Group>::insertItem(axom::sidre::Group*, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&) ()
   from /home/klevtsov/work/GEOSX/build-clang8-debug/lib/libgeosx_core.so
#11 0x00007ffff7a37a91 in axom::sidre::Group::createGroup(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&) ()
   from /home/klevtsov/work/GEOSX/build-clang8-debug/lib/libgeosx_core.so
#12 0x00007ffff6f4363d in geosx::dataRepository::ManagedGroup::setSidreGroup (name="siloOutput", parent=0x6037d0) at /home/klevtsov/work/GEOSX/src/coreComponents/dataRepository/ManagedGroup.cpp:58
#13 0x00007ffff6f43732 in geosx::dataRepository::ManagedGroup::ManagedGroup (this=0x649910, name="siloOutput", parent=0x6037d0)
    at /home/klevtsov/work/GEOSX/src/coreComponents/dataRepository/ManagedGroup.cpp:73
#14 0x00007ffff718ae2a in geosx::ExecutableGroup::ManagedGroup (this=0x649910) at /home/klevtsov/work/GEOSX/src/coreComponents/dataRepository/ExecutableGroup.hpp:46
#15 0x00007ffff71958c3 in geosx::OutputBase::OutputBase (this=0x649910, name="siloOutput", parent=0x6037d0) at /home/klevtsov/work/GEOSX/src/coreComponents/managers/Outputs/OutputBase.cpp:34
#16 0x00007ffff7196900 in geosx::SiloOutput::SiloOutput (this=0x649910, name="siloOutput", parent=0x6037d0) at /home/klevtsov/work/GEOSX/src/coreComponents/managers/Outputs/SiloOutput.cpp:36
#17 0x00007ffff7197a90 in std::make_unique<geosx::SiloOutput, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, geosx::dataRepository::ManagedGroup* const&> (
    __args=@0x7fffffffbec8: 0x6037d0, __args=@0x7fffffffbec8: 0x6037d0) at /usr/bin/../lib/gcc/x86_64-linux-gnu/8/../../../../include/c++/8/bits/unique_ptr.h:831
#18 0x00007ffff71979e3 in cxx_utilities::CatalogEntry<geosx::OutputBase, geosx::SiloOutput, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, geosx::dataRepository::ManagedGroup* const>::Allocate (this=0x45d270, args=0x6037d0, args=0x6037d0) at /home/klevtsov/work/GEOSX/src/coreComponents/cxx-utilities/src/src/ObjectCatalog.hpp:222
#19 0x00007ffff71949bd in cxx_utilities::CatalogInterface<geosx::OutputBase, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, geosx::dataRepository::ManagedGroup* const>::Factory (objectTypeName="Silo", args=0x6037d0, args=0x6037d0) at /home/klevtsov/work/GEOSX/src/coreComponents/cxx-utilities/src/src/ObjectCatalog.hpp:139
#20 0x00007ffff7194740 in geosx::OutputManager::CreateChild (this=0x6037d0, childKey="Silo", childName="siloOutput") at /home/klevtsov/work/GEOSX/src/coreComponents/managers/Outputs/OutputManager.cpp:47
#21 0x00007ffff6f44bc6 in geosx::dataRepository::ManagedGroup::ProcessInputFileRecursive (this=0x6037d0, targetNode=...)
    at /home/klevtsov/work/GEOSX/src/coreComponents/dataRepository/ManagedGroup.cpp:171
#22 0x00007ffff6f44c72 in geosx::dataRepository::ManagedGroup::ProcessInputFileRecursive (this=0x7fffffffd0b0, targetNode=...)
    at /home/klevtsov/work/GEOSX/src/coreComponents/dataRepository/ManagedGroup.cpp:178
#23 0x00007ffff71c58cd in geosx::ProblemManager::ParseInputFile (this=0x7fffffffd0b0) at /home/klevtsov/work/GEOSX/src/coreComponents/managers/ProblemManager.cpp:641
#24 0x000000000040e870 in main (argc=3, argv=0x7fffffffdfc8) at /home/klevtsov/work/GEOSX/src/main/main.cpp:59

The unit tests are currently kind of difficult to read for newcomers. We need to generate a set of examples that illustrate common usage for each of the classes. Also we should document the examples in sphinx. The sphinx is probably sufficient, but it would be nice to have the usage examples coded up and run with the unit tests s.t. any changes will be caught. Then we can just include the examples into the sphinx directly.

Henri Calandra had issues compiling GEOSX with GCC 9.2.1. The problem were asserts of the form static_assert( !typeManipulation ::always_true<T>, … ). The issue appears to be in the implementation of always_true whose purpose is to obfuscate the conditional so that the static_assert is only triggered if the function is used, not when it's parsed. It appears that the current implementation is too simple.

Add extensive documentation for the interplay with RAJA and CHAI.

Add documentation on when you can take a slice of an object, and what the lifetime of that slice must then be.

Currently, if we have const object that has an Array<double> m_array, our access to that array will be through the interface of Array const. This implies that you can't modify the values of m_array. However if you have an Array<double> const & array = m_array then you can get an ArrayView<double> const & arrayView = array. Now you can edit the values of array through arrayView. There is no way to stop this because Array derives from ArrayView. However, if Array did not derive from ArrayView, then we could control the conversion from Array -> ArrayView and disallow conversion of an Array<double> const to an ArrayView<double> const by only providing the user defined conversion Array<T>::operator ArrayView<T const>() and not Array<T>::operator ArrayView<T>(). Of course this UDC would not trigger a touch in the underlying chai::ManagedArray.

@corbett5 @klevzoff any thoughts?