Parallel version of lu_rec with task groups

lu_rec_grp_tasks.cpp

/*parallhlopoihsh sthn lower solve*/

void lower_solve(double ** a, int x1, int y1, double ** l, int x2, int y2, int N)
{

    int hn;
    tbb::task_group g;

    /* Check base case. */
    if (N <= BLOCK) {
        block_lower_solve(a,x1,y1,l,x2,y2,N);
        return;
    }
         
    hn = N / 2;

       /* Solve with recursive calls. */
	g.run([=]{aux_lower_solve(a,x1,y1,a,x1+hn,y1,l,x2,y2,hn);});
	g.run([=]{aux_lower_solve(a,x1,y1+hn,a,x1+hn,y1+hn,l,x2,y2,hn);});

    g.wait();
}

/*parallhlopoihsh sthn upper solve*/

void upper_solve(double ** a, int x1, int y1, double ** u, int x2, int y2, int N)
{
    int hn;
    tbb::task_group g;

    /* Check base case. */
     if (N <= BLOCK) {
        block_upper_solve(a,x1,y1,u,x2,y2,N);
	return;
 	}
    /* Break matrices into 4 pieces. */
 	hn = N / 2;

    /* Solve with recursive calls. */
	g.run([=]{aux_upper_solve(a,x1,y1,a,x1,y1+hn,u,x2,y2,hn);});
    g.run([=]{aux_upper_solve(a,x1+hn,y1,a,x1+hn,y1+hn,u,x2,y2,hn);});
    
    g.wait();
}

/*parallhlopoihsh sthn schur*/

void schur(double ** a, int x1, int y1, double ** v, int x2, int y2, double ** w, int x3, int y3, int N)
{

    int hn;
    tbb::task_group g;

    /* Check base case. */
    if (N <= BLOCK) {
        block_schur(a, x1, y1, v, x2, y2, w, x3, y3, N);
        return;
    }
    hn = N / 2;
	

    /* Form Schur complement with recursive calls. */
	g.run([=]{schur(a,x1,y1,v,x2,y2,w,x3,y3,hn);});
	g.run([=]{schur(a,x1,y1+hn,v,x2,y2,w,x3,y3+hn,hn);});
	g.run([=]{schur(a,x1+hn,y1,v,x2+hn,y2,w,x3,y3,hn);});
	g.run([=]{schur(a,x1+hn,y1+hn,v,x2+hn,y2,w,x3,y3+hn,hn);});


	g.run([=]{schur(a,x1,y1,v,x2,y2+hn,w,x3+hn,y3,hn);});
	g.run([=]{schur(a,x1,y1+hn,v,x2,y2+hn,w,x3+hn,y3+hn,hn);});
	g.run([=]{schur(a,x1+hn,y1,v,x2+hn,y2+hn,w,x3+hn,y3,hn);});
	g.run([=]{schur(a,x1+hn,y1+hn,v,x2+hn,y2+hn,w,x3+hn,y3+hn,hn);});

    g.wait();
}

/*parallhlopoihsh sthn lu*/
void lu(double ** a, int xs, int ys, int N)
{
    int hn;
    tbb::task_group g;

    /* Check base case. */
    if (N <= BLOCK) {
    block_lu(a, xs, ys,N);
    return;
    }

    hn = N / 2;

    /***** Compute LU decomposition on upper left tile A00 recursively*****/
    lu(a,xs,ys,hn);

    /***** Compute LU decomposition on upper right A01 and lower left A10 tiles via substitution recursively*****/
    g.run([=]{lower_solve(a,xs,ys+hn,a,xs,ys,hn);});
    g.run([=]{upper_solve(a,xs+hn,ys,a,xs,ys,hn);});
    
    g.wait();

    /***** Update lower right tile A11 recursively*****/
    schur(a,xs+hn,ys+hn,a,xs+hn,ys,a,xs,ys+hn,hn);

    /***** Compute LU decomposition on lower right tile A11 recursively*****/
    lu(a,xs+hn,ys+hn,hn);

}