Multithreading pthread program

Question

My multithreaded version of code works slowly than singlethreaded. Can you help me to understand why?

#include <iostream>
#include <vector>
#include <getopt.h>
#include <stdlib.h>
#include <cmath>
#include <pthread.h>
#include <unistd.h>
#include <time.h>

std::vector<std::vector<double> > phi, I;
std::vector<int> condition;
static int THREADS, X_SIZE, Y_SIZE;
pthread_barrier_t bar, syn;
pthread_mutex_t lock;
static double R = 10, C = 1e-6, dt, T;
long long clock_time() {
    struct timespec tp;
    clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &tp);
    return (long long)(tp.tv_nsec + (long long)tp.tv_sec * 1000000000ll);
}
void print_matrix(std::vector<std::vector<double> > v) {

    for(int i = 0; i < v.size(); i++) {
        for(int j = 0; j < v[i].size(); j++) {
            std::cout << v[i][j] << " ";
        }
        std::cout << std::endl;
    }
}

void print_matrix(std::vector<double> v) {
    for(int i = 0; i < v.size(); i++){ 
        std::cout << v[i] << std::endl;
    }
}

std::vector<double> solve_progon(std::vector<std::vector<double> > A, std::vector <double> B) {
    std::vector<double> P, Q, ans;
    int N = B.size();
    P.resize(N);
    Q.resize(N);
    ans.resize(N);
    P[0] = - A[0][1] / A[0][0];
    Q[0] = B[0]/A[0][0];
    //Pr9moy xod
    for(int i = 1; i < N; i++) {
        P[i] = A[i][i + 1]/(-A[i][i] - A[i][i - 1]*P[i-1]);
        Q[i] = (A[i][i - 1] * Q[i - 1] - B[i]) / (-A[i][i] - A[i][i-1] * P[i - 1]);
    }

    //Obratniy xod
    ans[N - 1] = (A[N-1][N-2] * Q[N-2] - B[N-1]) / (-A[N-1][N-1] - A[N-1][N-2] * P[N-2]);
    for(int i = N - 2; i >= 0; --i) {
        ans[i] = P[i] * ans[i+1] + Q[i];
    }
    return ans;
}

std::vector<double> check_solution(std::vector<std::vector<double> >A , std::vector<double> X, std::vector<double> B) {
    std::vector<double> diff;
    diff.resize(B.size());
    int size = B.size();
    for(int i = 0; i < size; i++) {
        double temp = 0;
        for(int j = 0; j < size; j++) {
            temp += A[i][j] * X[j];
        }
        diff[i] = temp - B[i];
    }
    return diff;
}

std::vector<double> solve(std::vector<double> left, std::vector<double> mid, std::vector<double> right, int j) {
    int size = left.size();
    std::vector<double> B;
    std::vector<std::vector<double> > A;
    A.resize(size);
    B.resize(size);
    //std::cout << left.size() << " " << mid.size() << " " << right.size() << std::endl;
    for(int i = 0; i < A.size(); i++) {
        A[i].resize(size);
    }
    for(int i = 0; i < A.size(); i++) {
        if(i - 1 >=0)
            A[i][i-1] = -1;
        A[i][i] = 2;

        if(i + 1 < A.size())
            A[i][i+1] = -1;
        B[i] = (left[i] - 2 * mid[i] + right[i]) * (dt / (R * C) - 1) + I[j][i] * dt / C ;
    }
    //print_matrix(A);
    std::vector<double> s = solve_progon(A, B);
    std::vector<double> check = check_solution(A, s, B);
    for (int i = 0; i < check.size(); i++) {
        if (fabs(check[i]) > 1e-10)
            std::cout << "ERROR TOO BIG" << check[i] << std::endl;
    }
    return s;
}

void * calc(void *thread) {
    long t = (long) thread;
    int start_index = t * (X_SIZE / THREADS); 
    int end_index = (t != THREADS - 1)?(t + 1) * (X_SIZE / THREADS) - 1: X_SIZE - 1;
    std::vector<std::vector<double> > local, next;
    std::vector<double> zeros;
    for (int i = 0; i < Y_SIZE; i++) {
        zeros.push_back(0);
    }
    double cur_time = 0;
    std::cout << "I am thread " << t << ". My start index: " << start_index << ". My end index: " << end_index << std::endl;
    while (cur_time < T) {
        /*if (start_index == 0) {
            local.push_back(zeros);
        } else {
            local.push_back(phi[start_index - 1]);
        }
        //pthread_mutex_lock(&lock);
        for (int i = start_index; i <= end_index; i++) {
            local.push_back(phi[i]);
        }
        //pthread_mutex_unlock(&lock);
        if (end_index == X_SIZE - 1) {
            local.push_back(zeros);
        }
        else {
            local.push_back(phi[end_index + 1]);
        }*/
        //print_matrix(local);
        //std::cout << local.size() << std::endl;
        for (int i = start_index; i <= end_index; i ++) {
            if (i == start_index && start_index == 0) 
                next.push_back(solve(zeros, phi[i], phi[i + 1], i - start_index));
            else if (i == end_index && end_index == X_SIZE - 1) 
                next.push_back(solve(phi[i - 1], phi[i], zeros, i - start_index));
            else 
                next.push_back(solve(phi[i - 1], phi[i], phi[i + 1], i - start_index));
        }
        cur_time += dt;
//        pthread_barrier_wait(&bar);
        for (int i = start_index; i <=end_index; i++) {
            phi[i] = next[i - start_index];
        }
        next.clear();
        pthread_barrier_wait(&syn);
    }
    pthread_exit(NULL);
}

int main(int argc, char **argv) {
    int time_to_modulate;
    int num_threads;
    int x_grid, y_grid;
    const char* short_options = "tpxy";
    int c;
    T = 2;
    dt = 1e-5;
    pthread_attr_t attr;
    pthread_attr_init(&attr);
    pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
    x_grid = atoi(argv[1]);
    y_grid = atoi(argv[2]);
    num_threads = atoi(argv[3]);
    THREADS = num_threads;
    time_to_modulate = 2;
    X_SIZE = x_grid;
    Y_SIZE = y_grid;
    phi.resize(x_grid);
    for(int i = 0; i < phi.size(); i++) {
        phi[i].resize(y_grid);
        for (int j = 0; j < phi[i].size(); j++) 
            phi[i][j] = (double) rand() / RAND_MAX;
    }
    I.resize(x_grid);
    for(int i = 0; i < I.size(); i++) {
        I[i].resize(y_grid);
    }
    I[1][1] = 10;
    pthread_barrier_init(&bar, NULL, THREADS);
    pthread_barrier_init(&syn, NULL, THREADS);
    pthread_mutex_init(&lock, NULL);
    /*std::vector<std::vector<double> > A;
    std::vector<double> B, X;
    int size = atoi(argv[1]);;
    A.resize(size);
    for(int i = 0; i < A.size(); i++) {
        A[i].resize(size);
    }
    B.resize(size);
    */
    pthread_barrier_init(&bar, NULL, THREADS);
    pthread_t *threads = new pthread_t[THREADS];
    unsigned long long start = clock_time();
    for (long i = 0; i < THREADS; i++) {
        if (pthread_create(&threads[i], NULL, calc, (void *)i) != 0) {
            std::cout << "Can't create thread " << i << std::endl;
        }
    }
    for (int i = 0; i < THREADS; i++) {
        pthread_join(threads[i], NULL);
    }
    std::cout << "It takes " << (double)(clock_time() - start) / 1e9 << std::endl;
    /*for(int i = 0; i < A.size(); i++) {
        if(i - 1 >=0)
            A[i][i-1] = (double)rand() / RAND_MAX;
        A[i][i] = (double)rand() / RAND_MAX;

        if(i + 1 < A.size())
            A[i][i+1] = (double)rand() / RAND_MAX;
        B[i] = (double)rand() / RAND_MAX;

    }
    //print_matrix(A);
    //print_matrix(B);
    X = solve_progon(A, B);
    std::vector<double> check = check_solution(A,X,B);
    for (int i = 0; i < check.size(); i++) {
        if (fabs(check[i]) > 1e-10)
            std::cout << "ERROR TOO BIG" << check[i] << std::endl;
    }*/
    //print_matrix(check_solution(A, X, B));
    /*while (1) {
        const struct option long_options[] = {
            {"delta_t", required_argument, 0, 0},
            {"threads", required_argument, 0, 0},
            {"x_size", required_argument, 0, 0},
            {"y_size", required_argument, 0, 0},
            {NULL, 0, NULL, 0}
        };
        int index = 0;
        c = getopt_long(argc, argv, "t:p:x:y:", long_options, &index);
        if (c == -1)
            break;
        std::cout << c << " " << optarg << std::endl;
        switch (c) {
            case 't':
                delta_t = atoi(optarg);
                break;
            case 'p':
                num_threads = atoi(optarg);
                break;
            case 'x':
                x_grid = atoi(optarg);
                break;
            case 'y':
                y_grid = atoi(optarg);
                break;
        }
    }
    std::cout << delta_t << " " << num_threads << " " << x_grid << " " << y_grid<< std::endl;*/

    return 0;
}

Answer 1

C++, unlike a lot of other languages, has what are known as value semantics. What this means is that when you pass a parameter such as a std::vector<T>, a full copy of that vector is made. Here, for every single parameter that takes a vector, you'll be doing a lot of unnecessary work.

Changing every function to take its parameters by const reference will likely make your program run much faster, single or multi-threaded.

void print_matrix(const std::vector<std::vector<double> >& v)

If you're using a vector of vectors as a matrix, I'd suggest you use a typedef (or in C++11, a using declaration) to cut down on the boilerplate:

// C++11
using matrix = std::vector<std::vector<double>>;

// C++98/03
typedef std::vector<std::vector<double> > matrix;

You might want to consider something similar for std::vector<double> as well.

A lot of the resize operations you perform on vectors can simply be moved into the constructor:

std::vector<double> solve(const std::vector<double>& left, const std::vector<double>& mid, 
                          const std::vector<double>& right, int j)
{
    // Note that std::vector<T>.size() doesn't actually return an int
    // It defines a size_type that should be used
    using size_type = std::vector<double>::size_type;
    const size_type size = left.size();

    // In C++11 the above can simply be:
    const auto size = left.size();

    std::vector<double> B(size);
    // Initialize each row of A with a std::vector<double> of the correct size.
    matrix A(size, std::vector<double>(size));

    ...
}

You might want to consider overloading operator<< to perform your printing:

std::ostream& operator<<(std::ostream& os, const matrix& m)
{
    // Code from print matrix goes here
    return os;
}

This way you print a matrix the same way you'd print any basic type in C++:

matrix m;
std::cout << m << "\n";

You have two functions called print_matrix, one that operates on matrices, and one that looks like it operates on row/column vectors. The second should be given a correct name, like print_row (or again, overloading operator<< for std::vector<double>).

If you aren't tied to pthreads, I'd suggest using a higher-level threading library (such as boost, or with C++11, the threading facilities that come as part of the language). Not only is pthreads less portable, it is also a much uglier API to work with.