Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
// 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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment