triangle.cpp

#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;
    }

}