added advanced assignment for More Features session. Code tested working.

D3-OpenMP/03-OpenMP-MoreFeatures/exercise/README.md

-Fill with exercise instructions.
+Assignment 1
+============
+
+In this assignment we are optimizing code for computing a sparse matrix vector product.
+
+We will use the file spmv.cc, (c++ only).
+
+Remember to load the gcc module as follows:
+
+module load gcc
+
+and then compile the code with
+
+g++ spmv.cc -o spmv -fopenmp -O3
+
+At first, the code does not take advantage of OpenMP work-sharing constructs.  This assignment
+is to use them to optimize the bottleneck operation of the code, which is the function spMv,
+which implements a sparse matrix vector product.
+
+You do not need to really understand the function being optimized, just see it as a 'black box'
+in terms of loops, parallel sections etc.
+
+The program takes a number of arguments.  For example, run it with
+
+/spmv 100000 100 1 10 100 
+
+to create a matrix with 100000 rows, with 100 nonzero entries per row, using 1 thread, checking every 10 iterations for iteration time, up to a maximum of 100 iterations.
+
+You can experiment with different parameters but you should aim to add work sharing constructs so that e.g. a 40 thread run gets a much smaller iteration time than a 1 thread run of the same code.
+
+ 

D3-OpenMP/03-OpenMP-MoreFeatures/exercise/spmv.cc

+#include <stdio.h>
+#include <iostream>	// cout
+#include <vector>	// vector
+#include <omp.h>
+#include <iomanip>	// setprecision
+
+using std::cout;
+using std::endl;
+using std::string;
+using std::fixed;
+using std::setprecision;
+
+class spM {
+	private:
+	public:
+		int nrows;
+		double** rows;	// rows, jagged array
+		int** colj;	// column position of jth entry in row i
+		int* rw;		// row weight
+		
+		spM(int nrows1, int w, int maxcol) {
+			nrows = nrows1;
+			rw = new int[nrows];
+			rows = new double*[nrows];
+			colj = new int*[nrows];
+			for (int i = 0; i < nrows; i++) {
+				// read ith row weight
+				rw[i] = w;
+				rows[i] = new double[rw[i]];
+				colj[i] = new int[rw[i]];
+				for (int j = 0; j < rw[i]; j++) {
+					colj[i][j] = rand() % maxcol;
+					rows[i][j] = ((double)rand())/((double)RAND_MAX);  // float between 0 and 1
+				}
+			}
+		}
+		~spM() {
+			for (int i = 0; i < nrows; i++) {
+				delete[] rows[i];
+				delete[] colj[i];
+			}
+			delete[] colj;
+			delete[] rows;
+			delete[] rw;
+		}
+};
+
+void spMv(int nrows, spM* M, double* w0, double* v0, int t_total, double* u0);
+void vset(int nrows, double* v, double* w);
+
+int main(int argc, char** argv)
+{
+	if (argc != 6) {
+		cout << endl << "Usage: ./spmv nrows row_weight num_threads checkiters maxiters" << endl << endl;
+		return 0;
+	}
+
+	int nrows = atoi(argv[1]);
+	int row_weight = atoi(argv[2]);
+	int num_threads = atoi(argv[3]);
+	int checkiters = atoi(argv[4]);
+	int maxiters = atoi(argv[5]);
+
+	omp_set_num_threads(num_threads);
+
+	// initialize vectors
+	double* b = new double[nrows];
+	double* wi = new double[nrows];
+	double* vi = new double[nrows];
+
+	// initialize thread separation array
+	double* u0 = new double[nrows * num_threads];
+
+	// initialize matrix
+	spM* M = new spM(nrows, row_weight, nrows);
+
+	// compute random b
+	cout << "Generating random starting vector b...";
+	for (int i = 0; i < nrows; i++) {
+		b[i] = ((double)rand())/((double)RAND_MAX);
+	}
+	cout << "Done." << endl;
+
+	vset(nrows, wi, b);
+
+	cout << "Starting spMv product iteration" << endl;
+
+	// *********************************************************
+	// main iteration
+	// *********************************************************
+	int iter = 0;
+	double start = omp_get_wtime();
+	while(true) {
+		double spMvstart = omp_get_wtime();
+		spMv(nrows, M, wi, vi, num_threads, u0);	// vi = M~*M*wi (mod ell)
+		double spMvstop = omp_get_wtime();
+		double spMvtime = spMvstop - spMvstart;
+		if (iter >= maxiters) break;
+		//cout << iter << ":" << mpz_get_str(NULL, 10, t) << endl; 
+		
+		vset(nrows, wi, vi);
+		iter++;
+		if (iter % checkiters == 0) {
+			double stop = omp_get_wtime();
+			double itertime = (stop - start)/checkiters;
+			cout << iter << ":iteration time: " << itertime << " seconds, spMv time: "
+				<< spMvtime << " seconds" << endl;
+			start = stop;
+		}
+	}
+
+	// release memory
+	delete M;
+	delete[] u0;
+	delete[] vi;
+	delete[] wi;
+	delete[] b;
+}
+
+
+void spMv(int nrows, spM* M, double* w0, double* v0, int t_total, double* u0)
+{
+	// v0 = M~*M*w0 (mod ell)
+	{
+		// first v0 = M*w0
+		
+		for (int i = 0; i < nrows; i++) {
+			v0[i] = 0.0;
+			for (int j = 0; j < M->rw[i]; j++) {
+				v0[i] = M->rows[i][j] * w0[M->colj[i][j]];
+			}
+		}
+		// now u0 = M~*v0.  This is the gnarly part
+		// first set u0 = 0
+		
+		for (int i = 0; i < nrows * t_total; i++) u0[i] = 0.0;
+		// now matrix transpose vector product
+		
+		for (int i = 0; i < nrows; i++) {
+			int t = omp_get_thread_num();
+			for (int j = 0; j < M->rw[i]; j++) {
+				u0[nrows*t + M->colj[i][j]] += M->rows[i][j] * v0[i];
+			}
+		}
+		// collapse u0 to a single vector
+		
+		for (int i = 0; i < nrows; i++) {
+			v0[i] = 0.0;
+			for (int t = 0; t < t_total; t++) {
+				v0[i] = v0[i] + u0[nrows*t + i];
+			}
+		}
+	}
+}
+
+void vset(int nrows, double* v, double* w)
+{
+	for (int i = 0; i < nrows; i++) {
+		v[i] = w[i];
+	}
+}
+