#include <flens/flens.h>
template <typename INT, typename T, typename VX, typename VY, typename VZ> void myComputationalRoutine(INT n, const T &alpha, const VX *x, INT incX, const VY *y, INT incY, VZ *z, INT incZ) { for (INT i=0, iX=0, iY=0, iZ=0; i<n; ++i, iX+=incX, iY+=incY, iZ+=incZ) { z[iZ] = alpha*(x[iX] + y[iY]); } } namespace flens { namespace blas { // // z = alpha*(x+y) // template <typename T, typename VX, typename VY> using MyCrazyVectorClosure = VectorClosure<OpMult, ScalarValue<T>, VectorClosure<OpAdd, DenseVector<VX>, DenseVector<VY> > >; template <typename T, typename VX, typename VY, typename VZ> void copy(const MyCrazyVectorClosure<T, VX, VY> &expr, DenseVector<VZ> &z) { const auto alpha = expr.left().value(); const auto x = expr.right().left(); const auto y = expr.right().right(); ASSERT(x.length()==z.length()); ASSERT(y.length()==z.length()); std::cerr << "calling myComputationalRoutine" << std::endl; myComputationalRoutine(z.length(), alpha, x.data(), x.stride(), y.data(), y.stride(), z.data(), z.stride()); } } } // namespace blas, flens #include <flens/flens.cxx> using namespace flens; using namespace std; int main() { typedef DenseVector<Array<double> > RealDenseVector; RealDenseVector x(5), y(5), z(5); x = 1; y = 2; z = 0; z = 3*(x+y); cout << "z = " << z << endl; } |