Commit 858d578a authored by Oisin Robinson's avatar Oisin Robinson
Browse files

add solution for smpv.cc advanced assignment

parent a7521198
#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)
#pragma omp parallel
{
// first v0 = M*w0
#pragma omp for
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
#pragma omp for
for (int i = 0; i < nrows * t_total; i++) u0[i] = 0.0;
// now matrix transpose vector product
#pragma omp for
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
#pragma omp for
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];
}
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment