// x <- A * d * A' * b.  x is allowed to be the same as b.
template <class T, int N>
static inline void eigmult(const T (&A)[N][N],
                           const T d[N],
                           const T b[N],
                           T x[N])
{
    T e[N];
    for (int i = 0; i < N; i++) {
        T sum = 0;
        for (int j = 0; j < N; j++)
            sum += A[j][i] * b[j];
        e[i] = sum * d[i];
    }
    for (int i = 0; i < N; i++) {
        T sum = 0;
        for (int j = 0; j < N; j++)
            sum += A[i][j] * e[j];
        x[i] = sum;
    }
}


he calls the code in his program like this:

// Solve system
eigmult<double,6>(A, einv, b);