scal.cxx

// Stupidly simple scal & axpy example

#include <cstdlib>
#include <ctime>
#include <iostream>

#ifndef REAL
# define REAL double
#endif

#ifndef N
# define N 10000
#endif

typedef REAL real;

template<typename RealT>
void scal(RealT *x, RealT s, int n) {
    for (int i = 0; i < n; i++) {
        x[i] *= s;
    }
}

template<typename RealT>
void axpy(RealT a, RealT *x, RealT *y, int n) {
    for (int i = 0; i < n; i++) {
        x[i] = a * x[i] + y[i];
    }
}

template<typename T>
void init(T *x, int n) {
    for (int i = 0; i < n; i++) {
        x[i] = static_cast<T>(rand()) / RAND_MAX;
    }
}

int main() {
    srand(time(NULL));
    real *x = new real[N];
    real *y = new real[N];
    real a = static_cast<real>(rand()) / RAND_MAX;

    init(x, N);
    init(y, N);

    scal(x, a, N);
    axpy(a, x, y, N);

    std::cout << x[rand() % N] << std::endl; // avoid optimization

    delete [] x;
    delete [] y;
    return 0;
}