task-for-matrix-vector-multiplication-and-adding-by-openmp-parallel

# include <Eigen/Dense>
# include <iostream>

// # include <mpc_omp_task_label.h>

# define EPSILON 10e-6

using Eigen::MatrixXd;
using Eigen::VectorXd;
using std::cout;
using std::cerr;
using std::endl;

int
main(int argc, char ** argv)
{
    if (argc != 3)
    {
        cerr << "usage: " << argv[0] << " [n] [m]" << endl;
        return 1;
    }
    int n = atoi(argv[1]);
    int m = atoi(argv[2]);
    # pragma omp parallel
    {
        # pragma omp single
        {
            // random matrices, vector
            MatrixXd A   = MatrixXd::Random(n, n);
            MatrixXd C   = MatrixXd::Random(n, n);
            MatrixXd R   = MatrixXd::Random(n, n);
            MatrixXd D1  = MatrixXd::Random(n, n);
            MatrixXd D2  = MatrixXd::Random(n, n);
            MatrixXd CT  = C.transpose();
            MatrixXd RT  = R.transpose();
            MatrixXd D1T = D1.transpose();
            MatrixXd D2T = D2.transpose();
            VectorXd x   = VectorXd::Random(n);

            // sequential result storage
            VectorXd y_seq(n);
            double t0 = omp_get_wtime();
            // MPC_OMP_TASK_SET_LABEL("sequential");
            # pragma omp task if(0) default(shared)
            {
                y_seq = (A + C + CT + R + RT + D1 + D1T + D2 + D2T) * x;
            }
            double t1 = omp_get_wtime();
            cout << "sequential version took " << (t1 - t0) << " s." << endl;

            // taskwait result storage
            VectorXd y_tw(n);
            VectorXd yA(n),yC(n), yCt(n), yR(n), yRt(n), yD1(n), yD1t(n), yD2(n), yD2t(n);

            double t2 = omp_get_wtime();

            #pragma omp task default(shared)
            yA = A * x;

            #pragma omp task default(shared)
            yC = C * x;

            #pragma omp task default(shared)
            yCt = CT * x;

            #pragma omp task default(shared)
            yR = R * x;

            #pragma omp task default(shared)
            yRt = RT * x;

            #pragma omp task default(shared)
            yD1 = D1 * x;

            #pragma omp task default(shared)
            yD1t = D1T * x;

            #pragma omp task default(shared)
            yD2 = D2 * x;

            #pragma omp task default(shared)
            yD2t = D2T * x;

            #pragma omp taskwait
            y_tw = yA + yC + yCt + yR + yRt + yD1 + yD1t + yD2 + yD2t;
            double t3 = omp_get_wtime();

            cout << "taskwait version took " << (t3 - t2) << " s.";
            cout << '(' << ((y_seq - y_tw).norm() < EPSILON ? "SUCCESS" : "FAILURE") << ')';
            cout << endl;

            // with task dependences and block operations
            MatrixXd M(n, n);
            VectorXd y_task(n);

            // 1 memory address per block
            int b = n / m;
            char Md[m][m];  (void)Md; // 'fake' data dependences
            char yd[m];     (void)yd; // 'fake' data dependences
            int m1, m2;

            double t4 = omp_get_wtime();
            for (m1 = 0 ; m1 < m ; ++m1)
            {
                for (m2 = 0 ; m2 < m ; ++m2)
                {
                    // MPC_OMP_TASK_SET_LABEL("add(%d, %d)", m1, m2);
                    # pragma omp task default(shared) firstprivate(m1, m2) depend(out: Md[m1][m2]) priority(1)
                    {
                        int b1 = m1 * b;
                        int b2 = m2 * b;
                        int i1 = b1;
                        int i2 = ((b1 + b < n) ? (b1 + b) : n);
                        int j1 = b2;
                        int j2 = ((b2 + b < n) ? (b2 + b) : n);
                        int si = i2 - i1;
                        int sj = j2 - j1;
                        M.block(i1, j1, si, sj) = ( A.block  (i1, j1, si, sj) +
                                                    C.block  (i1, j1, si, sj) +
                                                    CT.block (i1, j1, si, sj) +
                                                    R.block  (i1, j1, si, sj) +
                                                    RT.block (i1, j1, si, sj) +
                                                    D1.block (i1, j1, si, sj) +
                                                    D1T.block(i1, j1, si, sj) +
                                                    D2.block (i1, j1, si, sj) +
                                                    D2T.block(i1, j1, si, sj));
                    }
                }
            }
            for (m1 = 0 ; m1 < m ; ++m1)
            {
                // MPC_OMP_TASK_SET_LABEL("mult(%d)", m1);
                # pragma omp task default(shared) firstprivate(m1) depend(iterator(m2=0:m), in: Md[m1][m2]) depend(out: yd[m1])
                {
                    int b1 = m1 * b;
                    int i1 = b1;
                    int i2 = ((b1 + b < n) ? (b1 + b) : n);
                    int s  = i2 - i1;
                    y_task.segment(i1, s) = M.block(i1, 0, s, n) * x;
                }
            }

            //MPC_OMP_TASK_SET_LABEL("norm");
            # pragma omp task default(shared) depend(iterator(m2=0:m), in: yd[m2])
            {
                double t5 = omp_get_wtime();
                cout << "task version took " << (t5 - t4) << " s.";
                cout << '(' << ((y_seq - y_task).norm() < EPSILON ? "SUCCESS" : "FAILURE") << ')';
                cout << endl;
            }
        }
    }
    return 0;
}