-
-
Save ddemidov/1a67493145387cb7212c79e4015fe482 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#include <iostream> | |
#include <vector> | |
#include <atomic> | |
#include <chrono> | |
#include <cmath> | |
#include <omp.h> | |
double now() { | |
return omp_get_wtime(); | |
} | |
void TreeSolve(const int i, | |
const std::vector<int>& ptr, | |
std::vector<std::atomic<int>>& deg, | |
const std::vector<double>& rhs, | |
std::vector<double>& sum, | |
const std::vector<double>& val, | |
std::vector<double>& x, | |
const std::vector<int>& row) | |
{ | |
//note that there is not busy wait here. | |
//mark as negative deg[i] so it will not be executed later | |
#pragma omp task untied shared(ptr,deg,rhs,sum,val,x,row) | |
{ | |
deg[i]--; | |
int col_beg = ptr[i]; | |
int col_end = ptr[i+1]; | |
double xi = (rhs[i] - sum[i]) / val[col_beg]; | |
x[i] = xi; | |
for(int j = col_beg + 1; j < col_end; ++j) { | |
int r = row[j]; | |
double d = val[j] * xi; | |
if(r!=i) //other rows | |
{ | |
#pragma omp atomic | |
sum[r] += d; | |
std::cout << sum[r] << " "; | |
--deg[r]; | |
if(deg[r]==0) | |
{ | |
//std::cout << "spawning task for r = " << r << std::endl; | |
TreeSolve(r,ptr,deg,rhs,sum, val,x,row); | |
} | |
} | |
std::cout << std::endl; | |
} | |
} //task must terminate here so not to be kept alive | |
//std::cout << "finished task for row = " << i << std::endl; | |
} | |
int main() { | |
int n = 10000000; | |
std::vector<double> rhs(n, 1), x(n); | |
// Parallel | |
{ | |
std::vector<int> ptr, row; | |
std::vector<std::atomic<int>> deg(n); | |
std::vector<double> val, sum(n, 0.0); | |
ptr.reserve(n + 1); ptr.push_back(0); | |
row.reserve(n * 3); | |
val.reserve(n * 3); | |
for(int i = 0; i < n; ++i) { | |
deg[i] = 0; | |
if (i > 2) ++deg[i]; | |
if (i > 0) ++deg[i]; | |
row.push_back(i); | |
val.push_back(2); | |
if (i + 1 < n) { | |
row.push_back(i+1); | |
val.push_back(-1); | |
} | |
if (i + 3 < n) { | |
row.push_back(i+3); | |
val.push_back(-1); | |
} | |
ptr.push_back(val.size()); | |
} | |
double tic = now(); | |
// Parallel solve, column-wise matrix | |
{ | |
#pragma omp parallel | |
{ | |
#pragma omp single | |
for(int i = 0; i < n; ++i) | |
{ | |
if(deg[i] == 0) | |
{ | |
std::cout << "inside main loop for i = " << i << std::endl; | |
TreeSolve(i,ptr,deg,rhs,sum,val,x,row); | |
} | |
} | |
} | |
double toc = now(); | |
std::cout << "time: " << toc - tic << std::endl; | |
// Check residual | |
std::vector<double> res = rhs; | |
for(int i = 0; i < n; ++i) { | |
for(int j = ptr[i]; j < ptr[i+1]; ++j) { | |
res[row[j]] -= val[j] * x[i]; | |
} | |
} | |
double error = 0; | |
for(int i = 0; i < n; ++i) | |
error += res[i] * res[i]; | |
std::cout << "res: " << sqrt(error) << std::endl; | |
} | |
} | |
// Serial | |
{ | |
std::vector<int> ptr, col; | |
std::vector<double> val; | |
ptr.reserve(n + 1); ptr.push_back(0); | |
col.reserve(n * 3); | |
val.reserve(n * 3); | |
for(int i = 0; i < n; ++i) { | |
if (i > 2) { | |
col.push_back(i - 3); | |
val.push_back(-1); | |
} | |
if (i > 0) { | |
col.push_back(i-1); | |
val.push_back(-1); | |
} | |
col.push_back(i); | |
val.push_back(2); | |
ptr.push_back(val.size()); | |
} | |
double tic = now(); | |
for(int i = 0; i < n; ++i) { | |
double s = rhs[i]; | |
for(int j = ptr[i], e = ptr[i+1]; j < e; ++j) { | |
int c = col[j]; | |
double v = val[j]; | |
if (c < i) { | |
s -= v * x[c]; | |
} else { | |
x[i] = s / v; | |
} | |
} | |
} | |
double toc = now(); | |
std::cout << "time: " << toc - tic << std::endl; | |
// Check residual | |
double error = 0; | |
for(int i = 0; i < n; ++i) { | |
double res = rhs[i]; | |
for(int j = ptr[i]; j < ptr[i+1]; ++j) { | |
res -= val[j] * x[col[j]]; | |
} | |
error += res * res; | |
} | |
std::cout << "res: " << sqrt(error) << std::endl; | |
} | |
} | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment