|
23 | 23 | #ifndef _H_GRB_DENSEREF_BLAS1 |
24 | 24 | #define _H_GRB_DENSEREF_BLAS1 |
25 | 25 |
|
| 26 | +#include <graphblas/backends.hpp> |
| 27 | +#include <graphblas/config.hpp> |
| 28 | +#include <graphblas/rc.hpp> |
| 29 | + |
26 | 30 | namespace grb { |
27 | 31 |
|
| 32 | + /** |
| 33 | + * Calculates the dot product, \f$ \alpha = (x,y) \f$, under a given additive |
| 34 | + * monoid and multiplicative operator. |
| 35 | + * |
| 36 | + * @tparam descr The descriptor to be used (descriptors::no_operation |
| 37 | + * if left unspecified). |
| 38 | + * @tparam Ring The semiring type to use. |
| 39 | + * @tparam OutputType The output type. |
| 40 | + * @tparam InputType1 The input element type of the left-hand input vector. |
| 41 | + * @tparam InputType2 The input element type of the right-hand input vector. |
| 42 | + * |
| 43 | + * @param[in,out] z The output element \f$ z + \alpha \f$. |
| 44 | + * @param[in] x The left-hand input vector. |
| 45 | + * @param[in] y The right-hand input vector. |
| 46 | + * @param[in] addMonoid The additive monoid under which the reduction of the |
| 47 | + * results of element-wise multiplications of \a x and |
| 48 | + * \a y are performed. |
| 49 | + * @param[in] anyop The multiplicative operator under which element-wise |
| 50 | + * multiplications of \a x and \a y are performed. This can |
| 51 | + * be any binary operator. |
| 52 | + * |
| 53 | + * By the definition that a dot-product operates under any additive monoid and |
| 54 | + * any binary operator, it follows that a dot-product under any semiring can be |
| 55 | + * trivially reduced to a call to this version instead. |
| 56 | + * |
| 57 | + * @return grb::MISMATCH When the dimensions of \a x and \a y do not match. All |
| 58 | + * input data containers are left untouched if this exit |
| 59 | + * code is returned; it will be as though this call was |
| 60 | + * never made. |
| 61 | + * @return grb::SUCCESS On successful completion of this call. |
| 62 | + * |
| 63 | + * \parblock |
| 64 | + * \par Performance semantics |
| 65 | + * -# This call takes \f$ \Theta(n/p) \f$ work at each user process, where |
| 66 | + * \f$ n \f$ equals the size of the vectors \a x and \a y, and |
| 67 | + * \f$ p \f$ is the number of user processes. The constant factor |
| 68 | + * depends on the cost of evaluating the addition and multiplication |
| 69 | + * operators. A good implementation uses vectorised instructions |
| 70 | + * whenever the input domains, output domain, and the operators used |
| 71 | + * allow for this. |
| 72 | + * |
| 73 | + * -# This call takes \f$ \mathcal{O}(1) \f$ memory beyond the memory used |
| 74 | + * by the application at the point of a call to this function. |
| 75 | + * |
| 76 | + * -# This call incurs at most |
| 77 | + * \f$ n( \mathit{sizeof}(\mathit{D1}) + \mathit{sizeof}(\mathit{D2}) ) + \mathcal{O}(p) \f$ |
| 78 | + * bytes of data movement. |
| 79 | + * |
| 80 | + * -# This call incurs at most \f$ \Theta(\log p) \f$ synchronisations |
| 81 | + * between two or more user processes. |
| 82 | + * |
| 83 | + * -# A call to this function does result in any system calls. |
| 84 | + * \endparblock |
| 85 | + * |
| 86 | + * \note This requires an implementation to pre-allocate \f$ \Theta(p) \f$ |
| 87 | + * memory for inter-process reduction, if the underlying communication |
| 88 | + * layer indeed requires such a buffer. This buffer may not be allocated |
| 89 | + * (nor freed) during a call to this function. |
| 90 | + * |
| 91 | + * \parblock |
| 92 | + * \par Valid descriptors |
| 93 | + * -# grb::descriptors::no_operation |
| 94 | + * -# grb::descriptors::no_casting |
| 95 | + * -# grb::descriptors::dense |
| 96 | + * \endparblock |
| 97 | + * |
| 98 | + * If the dense descriptor is set, this implementation returns grb::ILLEGAL if |
| 99 | + * it was detected that either \a x or \a y was sparse. In this case, it shall |
| 100 | + * otherwise be as though the call to this function had not occurred (no side |
| 101 | + * effects). |
| 102 | + * |
| 103 | + * \note The standard, in contrast, only specifies undefined behaviour would |
| 104 | + * occur. This implementation goes beyond the standard by actually |
| 105 | + * specifying what will happen. |
| 106 | + */ |
| 107 | + template< |
| 108 | + Descriptor descr = descriptors::no_operation, |
| 109 | + class AddMonoid, class AnyOp, |
| 110 | + typename OutputType, typename InputType1, typename InputType2, |
| 111 | + typename InputView1, typename InputView2, |
| 112 | + typename InputStorage1, typename InputStorage2, |
| 113 | + typename InputCoords1, typename InputCoords2 |
| 114 | + > |
| 115 | + RC dot( OutputType &z, |
| 116 | + const VectorView< InputType1, InputView1, InputStorage1, reference_dense, InputCoords1 > &x, |
| 117 | + const VectorView< InputType2, InputView2, InputStorage2, reference_dense, InputCoords2 > &y, |
| 118 | + const AddMonoid &addMonoid = AddMonoid(), |
| 119 | + const AnyOp &anyOp = AnyOp(), |
| 120 | + const typename std::enable_if< !grb::is_object< OutputType >::value && |
| 121 | + !grb::is_object< InputType1 >::value && |
| 122 | + !grb::is_object< InputType2 >::value && |
| 123 | + grb::is_monoid< AddMonoid >::value && |
| 124 | + grb::is_operator< AnyOp >::value, |
| 125 | + void >::type * const = NULL |
| 126 | + ) { |
| 127 | + (void)z; |
| 128 | + (void)x; |
| 129 | + (void)y; |
| 130 | + (void)addMonoid; |
| 131 | + (void)anyOp; |
| 132 | + // static sanity checks |
| 133 | + // NO_CAST_ASSERT( ( !( descr & descriptors::no_casting ) || std::is_same< InputType1, typename AnyOp::D1 >::value ), "grb::dot", |
| 134 | + // "called with a left-hand vector value type that does not match the first " |
| 135 | + // "domain of the given multiplicative operator" ); |
| 136 | + // NO_CAST_ASSERT( ( !( descr & descriptors::no_casting ) || std::is_same< InputType2, typename AnyOp::D2 >::value ), "grb::dot", |
| 137 | + // "called with a right-hand vector value type that does not match the second " |
| 138 | + // "domain of the given multiplicative operator" ); |
| 139 | + // NO_CAST_ASSERT( ( !( descr & descriptors::no_casting ) || std::is_same< typename AddMonoid::D3, typename AnyOp::D1 >::value ), "grb::dot", |
| 140 | + // "called with a multiplicative operator output domain that does not match " |
| 141 | + // "the first domain of the given additive operator" ); |
| 142 | + // NO_CAST_ASSERT( ( !( descr & descriptors::no_casting ) || std::is_same< OutputType, typename AddMonoid::D2 >::value ), "grb::dot", |
| 143 | + // "called with an output vector value type that does not match the second " |
| 144 | + // "domain of the given additive operator" ); |
| 145 | + // NO_CAST_ASSERT( ( !( descr & descriptors::no_casting ) || std::is_same< typename AddMonoid::D3, typename AddMonoid::D2 >::value ), "grb::dot", |
| 146 | + // "called with an additive operator whose output domain does not match its " |
| 147 | + // "second input domain" ); |
| 148 | + // NO_CAST_ASSERT( ( !( descr & descriptors::no_casting ) || std::is_same< OutputType, typename AddMonoid::D3 >::value ), "grb::dot", |
| 149 | + // "called with an output vector value type that does not match the third " |
| 150 | + // "domain of the given additive operator" ); |
| 151 | + |
| 152 | + // dynamic sanity check |
| 153 | + // const size_t n = internal::getCoordinates( y ).size(); |
| 154 | + // if( internal::getCoordinates( x ).size() != n ) { |
| 155 | + // return MISMATCH; |
| 156 | + // } |
| 157 | + |
| 158 | + // // cache nnzs |
| 159 | + // const size_t nnzx = internal::getCoordinates( x ).nonzeroes(); |
| 160 | + // const size_t nnzy = internal::getCoordinates( y ).nonzeroes(); |
| 161 | + |
| 162 | + // // catch trivial case |
| 163 | + // if( nnzx == 0 && nnzy == 0 ) { |
| 164 | + // return SUCCESS; |
| 165 | + // } |
| 166 | + |
| 167 | + // // dot will be computed out-of-place here. A separate field is needed because |
| 168 | + // // of possible multi-threaded computation of the dot. |
| 169 | + // OutputType oop = addMonoid.template getIdentity< OutputType >(); |
| 170 | + |
| 171 | + // if descriptor says nothing about being dense... |
| 172 | + RC ret = SUCCESS; |
| 173 | + // if( !( descr & descriptors::dense ) ) { |
| 174 | + // // check if inputs are actually dense... |
| 175 | + // if( nnzx == n && nnzy == n ) { |
| 176 | + // // call dense implementation |
| 177 | + // ret = internal::dot_generic< descr | descriptors::dense >( oop, x, y, addMonoid, anyOp ); |
| 178 | + // } else { |
| 179 | + // // pass to sparse implementation |
| 180 | + // ret = internal::dot_generic< descr >( oop, x, y, addMonoid, anyOp ); |
| 181 | + // } |
| 182 | + // } else { |
| 183 | + // // descriptor says dense, but if any of the vectors are actually sparse... |
| 184 | + // if( nnzx < n || nnzy < n ) { |
| 185 | + // return ILLEGAL; |
| 186 | + // } else { |
| 187 | + // // all OK, pass to dense implementation |
| 188 | + // ret = internal::dot_generic< descr >( oop, x, y, addMonoid, anyOp ); |
| 189 | + // } |
| 190 | + // } |
| 191 | + |
| 192 | + // fold out-of-place dot product into existing input, and exit |
| 193 | + // ret = ret ? ret : foldl( z, oop, addMonoid.getOperator() ); |
| 194 | + return ret; |
| 195 | + } |
| 196 | + |
| 197 | + /** |
| 198 | + * Provides a generic implementation of the dot computation on semirings by |
| 199 | + * translating it into a dot computation on an additive commutative monoid |
| 200 | + * with any multiplicative operator. |
| 201 | + * |
| 202 | + * For return codes, exception behaviour, performance semantics, template |
| 203 | + * and non-template arguments, @see grb::dot. |
| 204 | + */ |
| 205 | + template< |
| 206 | + Descriptor descr = descriptors::no_operation, class Ring, |
| 207 | + typename IOType, typename InputType1, typename InputType2, |
| 208 | + typename InputView1, typename InputView2, |
| 209 | + typename InputStorage1, typename InputStorage2, |
| 210 | + Backend backend, typename Coords1, typename Coords2 |
| 211 | + > |
| 212 | + RC dot( IOType &x, |
| 213 | + const VectorView< InputType1, InputView1, InputStorage1, backend, Coords1 > &left, |
| 214 | + const VectorView< InputType2, InputView2, InputStorage2, backend, Coords2 > &right, |
| 215 | + const Ring &ring = Ring(), |
| 216 | + const typename std::enable_if< |
| 217 | + !grb::is_object< InputType1 >::value && |
| 218 | + !grb::is_object< InputType2 >::value && |
| 219 | + !grb::is_object< IOType >::value && |
| 220 | + grb::is_semiring< Ring >::value, |
| 221 | + void >::type * const = NULL |
| 222 | + ) { |
| 223 | + // return grb::dot< descr >( x, |
| 224 | + return grb::dot( x, |
| 225 | + left, right, |
| 226 | + ring.getAdditiveMonoid(), |
| 227 | + ring.getMultiplicativeOperator() |
| 228 | + ); |
| 229 | + } |
| 230 | + |
| 231 | + /** |
| 232 | + * Provides a generic implementation of the 2-norm computation. |
| 233 | + * |
| 234 | + * Proceeds by computing a dot-product on itself and then taking the square |
| 235 | + * root of the result. |
| 236 | + * |
| 237 | + * This function is only available when the output type is floating point. |
| 238 | + * |
| 239 | + * For return codes, exception behaviour, performance semantics, template |
| 240 | + * and non-template arguments, @see grb::dot. |
| 241 | + * |
| 242 | + * @param[out] x The 2-norm of \a y. The input value of \a x will be ignored. |
| 243 | + * @param[in] y The vector to compute the norm of. |
| 244 | + * @param[in] ring The Semiring under which the 2-norm is to be computed. |
| 245 | + * |
| 246 | + * \warning This function computes \a x out-of-place. This is contrary to |
| 247 | + * standard ALP/GraphBLAS functions that are always in-place. |
| 248 | + * |
| 249 | + * \warning A \a ring is not sufficient for computing a two-norm. This |
| 250 | + * implementation assumes the standard <tt>sqrt</tt> function |
| 251 | + * must be applied on the result of a dot-product of \a y with |
| 252 | + * itself under the supplied semiring. |
| 253 | + */ |
| 254 | + template< |
| 255 | + Descriptor descr = descriptors::no_operation, class Ring, |
| 256 | + typename InputType, typename OutputType, |
| 257 | + typename InputView, |
| 258 | + typename InputStorage, |
| 259 | + Backend backend, typename Coords |
| 260 | + > |
| 261 | + RC norm2( OutputType &x, |
| 262 | + const VectorView< InputType, InputView, InputStorage, backend, Coords > &y, |
| 263 | + const Ring &ring = Ring(), |
| 264 | + const typename std::enable_if< |
| 265 | + std::is_floating_point< OutputType >::value, |
| 266 | + void >::type * const = NULL |
| 267 | + ) { |
| 268 | + RC ret = grb::dot< descr >( x, y, y, ring ); |
| 269 | + if( ret == SUCCESS ) { |
| 270 | + x = sqrt( x ); |
| 271 | + } |
| 272 | + return ret; |
| 273 | + } |
| 274 | + |
28 | 275 | } // end namespace ``grb'' |
29 | 276 |
|
30 | 277 | #endif // end ``_H_GRB_DENSEREF_BLAS1'' |
|
0 commit comments