reference and reference_omp grb::dot do not handle aliasing properly

[anyzelman]

If the two inputs are the same vector, then we should not take restricted pointers from both-- yet the reference and reference_omp implementations presently do do that.

Possibly relatedly in the context of #379 , the internal CI emits a warning that may be indirectly related:

In file included from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas/internalops.hpp:29,
                 from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas/ops.hpp:30,
                 from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas.hpp:519,
                 from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/tests/smoke/conjugate_gradient.cpp:28:
/builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas/base/internalops.hpp: In function 'grb::RC grb::internal::dot_generic(OutputType&, const grb::Vector<InputType, grb::reference, Coords>&, const grb::Vector<MaskType, grb::reference, Coords>&, const AddMonoid&, const AnyOp&, const grb::Phase&) [with unsigned int descr = 0; AddMonoid = grb::Monoid<grb::operators::add<double, double, double, grb::hyperdags>, grb::identities::zero>; AnyOp = grb::operators::conjugate_left_mul<double, double, double, grb::hyperdags>; OutputType = double; InputType1 = double; InputType2 = double; Coords = grb::internal::Coordinates<grb::reference>]':
/builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas/base/internalops.hpp:2772:60: warning: 'yy[0]' may be used uninitialized in this function [-Wmaybe-uninitialized]
 2772 |        *c = grb::utils::is_complex< IN1 >::conjugate( *a ) * *b;
      |             ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~
In file included from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas/blas1.hpp:29,
                 from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas.hpp:538,
                 from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/tests/smoke/conjugate_gradient.cpp:28:
/builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas/reference/blas1.hpp:9400:27: note: 'yy[0]' was declared here
 9400 |        typename AnyOp::D2 yy[ AnyOp::blocksize ];
      |                           ^~
In file included from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas/internalops.hpp:29,
                 from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas/ops.hpp:30,
                 from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas.hpp:519,
                 from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/tests/smoke/conjugate_gradient.cpp:28:
/builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas/base/internalops.hpp:2772:60: warning: 'xx[0]' may be used uninitialized in this function [-Wmaybe-uninitialized]
 2772 |        *c = grb::utils::is_complex< IN1 >::conjugate( *a ) * *b;
      |             ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~
In file included from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas/blas1.hpp:29,
                 from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas.hpp:538,
                 from /builds/zrc-von-neumann-lab/spatial-computing/graphblas/tests/smoke/conjugate_gradient.cpp:28:
/builds/zrc-von-neumann-lab/spatial-computing/graphblas/include/graphblas/reference/blas1.hpp:9399:27: note: 'xx[0]' was declared here
 9399 |        typename AnyOp::D1 xx[ AnyOp::blocksize ];
      |                           ^~

The resolution/suppression of this warning I'll leave to this issue also

[GiovaGa]

It is not clear to me what we want to do with this issue. After this issue is solved, what is the behavior we want?
Do we want:

• To assert that the inputs in the dot function are different (and so disallow this operation, by throwing an error or failed assertion), or
• To check if the inputs are the same, in which case we use a different algorithm to evaluate the dot?
  Or is this second option for a future release?

[djelovina]

I am not sure if the warning and the issue are directly related

1. This seems to be related to the issue of restricted pointers usage

   const InputType1 * __restrict__ a = internal::getRaw( x );
   const InputType2 * __restrict__ b = internal::getRaw( y );

   if x and y are the same vector we have violation of __restrict__

2. the warning

   warning: 'yy[0]' may be used uninitialized in this function
   

   and and others, come from

   								for( size_t k = 0; k < AnyOp::blocksize; ++k ) {
   									if( mask[ k ] ) {
   										apply( zz[ k ], xx[ k ], yy[ k ], anyOp ); 
   

   which is false positive since compiler cannot conclude that

   							for( size_t k = 0; k < AnyOp::blocksize; ++k, ++i ) {
   								if( mask[ k ] ) {
   									xx[ k ] = static_cast< typename AnyOp::D1 >( a[ coors_r.index( i ) ] );
   									yy[ k ] = static_cast< typename AnyOp::D2 >( b[ coors_r.index( i ) ] );
   								}
   							}
   

   properly initializes xx and yy because mask[ k ] = coors_q.assigned( coors_r.index( i ) );

Is this reasoning correct @anyzelman ?

[anyzelman]

What you write seems correct to me, and apologies to Giovanni -- I had missed his earlier comment. The MR that you have however doesn't close the issue, because an ALP program could call grb::dot( alpha, x, x, ring ); in which case we violate __restrict__, as you said at the start of the thought chain?

For Giovanni's comments:

To assert that the inputs in the dot function are different (and so disallow this operation, by throwing an error or failed assertion)

I think that would definitely be good to do if we don't already

To check if the inputs are the same, in which case we use a different algorithm to evaluate the dot?

Yes, I think such an alternative code path is needed.

[djelovina]

Yes, I think such an alternative code path is needed.

This PR does both

1. provides the alternative path by exposing 1-vector versions or norm2() which is not internally caling 2-vector function like dot(),
2. asserts that dot() forbids aliased calls.

It also shows that that we can express complex operations as a single-vector-to-scalar reduction by defining a wrapper type and monoid, which allows us to avoid the aliasing issue entirely, issue.