RK4 w/ OMP

rk4_omp.cpp

#include <cstdlib>
#include <iostream>
#include <vector>
#include <ctime>
#include <omp.h>

using namespace std;

#define dim 2
#define USE_DEFAULT 1			// skip user input
#define DO_HARLEM_SHAKE 1		// show the result
#ifdef _OPENMP					// patch
	#ifndef _OMP_H
		#define _OMP_H _OPENMP
	#endif // _OMP_H
#endif // _OPENMP

/***********************************************
* comment the following line out to use OpenMP *
***********************************************/
#undef _OMP_H

vector<double> f(double t)
{
	return vector<double>(dim, 1.);		// here you should return whatever f(t) should be equal to
}

vector<double> operator*(const vector<double>& v, double s)
{
	vector<double> _(v);
	unsigned i, n;
	n = _.size();
	#ifdef _OMP_H
		#pragma omp parallel for private(i)
	#endif
	for(i=0; i<n; ++i)
		_[i]*=s;
	return _;
}

vector<double> operator/(const vector<double>& v, double s)
{
	vector<double> _(v);
	unsigned i, n;
	n = _.size();
	#ifdef _OMP_H
		#pragma omp parallel for private(i)
	#endif
	for(i=0; i<n; ++i)
		_[i]/=s;
	return _;
}

vector<double> operator+(const vector<double>& v1, const vector<double>& v2)
{
	vector<double> _(v1);
	unsigned i, n;
	n = _.size();
	if(v2.size()!=n)
		throw("wrong size");
	#ifdef _OMP_H
		#pragma omp parallel for private(i)
	#endif
	for(i=0; i<n; ++i)
		_[i]=v1[i]+v2[i];
	return _;
}

vector<double> &operator+=(vector<double>& v1, const vector<double>& v2)
{
	unsigned i, n;
	n = v1.size();
	if(v2.size()!=n)
		throw("wrong size");
	#ifdef _OMP_H
		#pragma omp parallel for private(i)
	#endif
	for(i=0; i<n; ++i)
		v1[i]+=v2[i];
	return v1;
}

vector<double> operator*(const vector<vector<double> >& m, const vector<double>& v)
{
	vector<double> _(v);
	unsigned i, j, n;
	n = _.size();
	if(m.size()!=n)
		throw("wrong size");
	#ifdef _OMP_H
		#pragma omp parallel for private(i,j)
	#endif
	for(i=0; i<n; ++i)
	{
		_[i] = 0.;
		if(m[i].size()!=n)
			throw("wrong size");
		for(j=0; j<n; ++j)
			_[i]+=m[i][j]*v[j];
	}
	return _;
}

int main()
{
	vector<double> y(dim), k1(dim), k2(dim), k3(dim), k4(dim);
	int n, i, j;
	double t, h, h_2, ti, tf;
	vector<vector<double> > A(dim);
    cout << "Enter the initial value of t: ";
	#if USE_DEFAULT
		ti = 0.;
		cout << ti << endl;
	#else
		cin >> ti;
	#endif // USE_DEFAULT
	cout << "Enter the final value of t: ";
	#if USE_DEFAULT
		tf = 1.;
		cout << tf << endl;
	#else
		cin >> tf;
	#endif // USE_DEFAULT
	cout << "Enter the value of h (step): ";
	#if USE_DEFAULT
		h = .0001;
		cout << h << endl;
	#else
		cin >> h;
	#endif // USE_DEFAULT
	for(i=1; i<=dim; ++i)
	{
		cout << "Enter the initial value of the " << i << ((((i%10)==1)&&((i%100)!=11))?"st":((((i%10)==2)&&((i%100)!=12))?"nd":"th")) << " component of y: ";
		A[i-1].resize(dim);			// a user reads slowly enough to perform the task
		#if USE_DEFAULT
			y[i-1] = i;
			cout << y[i-1] << endl;
		#else
			cin >> y[i-1];
		#endif // USE_DEFAULT
	}
	cout << endl;
	// initialize A
	#ifdef _OMP_H
		#pragma omp parallel for private(i,j)
	#endif
	for(i=0; i<dim; ++i)
		for(j=0; j<dim; ++j)
			A[i][j]=(i==j)*(1-2*(i&1));

	n = int((tf-ti)/h);
	h_2 = h/2.;
	#ifdef _OMP_H
		cout << "multiple thread execution (max " << omp_get_max_threads() << " threads on " << omp_get_num_procs() << " procs.)" << endl;
	#else
		cout << "single thread execution" << endl;
	#endif
	#if DO_HARLEM_SHAKE
		cout << 't';
		for(j=0; j<dim; ++j)
			cout << "\ty(" << j+1 << ')';
		cout << endl;
	#endif // DO_HARLEM_SHAKE
	// tic
	time_t tic = time(NULL);
	for(i=0; i<=n; ++i)
	{
		t = ti + h*i;
		#if DO_HARLEM_SHAKE
			cout << t;
			for(j=0; j<dim; ++j)
				cout << '\t' << y[j];
			cout << endl;
		#endif // DO_HARLEM_SHAKE
		k1 = (A* y        + f(t    ))*h;
		k2 = (A*(y+k1/2.) + f(t+h_2))*h;
		k3 = (A*(y+k2/2.) + f(t+h_2))*h;
		k4 = (A*(y+k3   ) + f(t+h  ))*h;
		y += (k1+k4)/6.+(k2+k3)/3.;
	}
	// toc
	time_t toc = time(NULL);
	cout << "the calculation took " << difftime(toc, tic) << " seconds" << endl;
	cout << "done" << endl;
	return 0;
}