How to do Y = aX + Y in C++

Question

I am trying to do matrix operations using C++ STL containers. There are two vectors of sizes say Y and X of sizes m,n(m>n). I want to multiply X with a scalar and add it from a given index in Y. In this process I don't want to transform X (don't want to use std::transform). In fact the X's are columns in a matrix DPtr. One version I tried is given below.

    std::vector<double> D(&DPtr[index1],&DPtr[index1]+size_C);

    std::transform(D.begin(), D.end(), D.begin(),std::bind2nd(std::multiplies<double>(), val1*val2));

    std::transform (D.begin(), D.end(), YPtr.begin()+index2, YPtr.begin()+index2, std::plus<double>());

I am trying to copy the column in to a temporary vector and do operations on it.

Can some one help me in rewriting the code in a much simpler manner where I need not copy columns into another vector?

I am guessing I have to use std::for_each and lamda expression or a function call? But I am new to C++?

Just to give a lead, I want to write it as

std::for_each(YPtr.begin()+index2,YPtr.begin()+index2+(size_c-1),scalarAdd);

using a function scalarAdd or any lamda expression where I can access DPtr directly. Also can I write

YPtr.begin()+index2+(size_c-1)
     . 

as the second argument.Is it valid? Also Imagine I made the matrix as a C++-vector where all columns of DPtr matrix are stored in one single C++ vector D.

Visual Representation of my question

Answer 1

May I suggest you use a dedicated linear-algebra library like Eigen? With that you could simply write Y += X * a, and the library+compiler will figure out the best implementation for you.

---

Since you are using C++11, you could use std::transform together with a lambda (std::for_each is not recommended since you need to transform your Y).

std::transform(DPtr.begin() + index1, DPtr.begin() + index1 + size_c, // X
               YPtr.begin() + index2,                                 // Y_in
               YPtr.begin() + index2,                                 // Y_out
               [](double x, double y_in) { return a*x + y_in; });
               // y_out = a*x + y_in, for every entry.