parallelizing linear algebra using gsl library

user_2811074 0
Viewed 184

In my c++ scripts, I have many for loops to compute linear algebra operations. I am wondering what is the best way to make the loops parallel? One example is the following function which computes the kronecker product of two matrices.

void Kronecker(const gsl_matrix *K, const gsl_matrix *V, gsl_matrix *H) 
{
    for (size_t i=0; i<K->size1; i++) {
        for (size_t j=0; j<K->size2; j++) {
            gsl_matrix_view H_sub=gsl_matrix_submatrix (H, i*V->size1, j*V->size2, V->size1, V->size2);
            gsl_matrix_memcpy (&H_sub.matrix, V);
            gsl_matrix_scale (&H_sub.matrix, gsl_matrix_get (K, i, j));
        }
    }
    return;
}

How can I improve the computation time of my code, when I have for loops which can be parallel?

gsl linear-algebra parallel-processing c++
2 Answers

Without knowing the memory layout, allocations, syscalls, and potential side-effects in your underlying gsl calls, a really easy way to get parallelization is via OpenMP. That of course introduces a dependency and requires compiler support, but it's particularly effective on simple loops like yours. Untested and probably needs a bit more to ensure H is written properly, but something like:

#pragma omp parallel for private(i, j)
for (size_t i=0; i<K->size1; i++) {
    for (size_t j=0; j<K->size2; j++) {
        gsl_matrix_view H_sub=gsl_matrix_submatrix (H, i*V->size1, j*V->size2, V->size1, V->size2);
        gsl_matrix_memcpy (&H_sub.matrix, V);
        gsl_matrix_scale (&H_sub.matrix, gsl_matrix_get (K, i, j));
    }
}

See https://curc.readthedocs.io/en/latest/programming/OpenMP-C.html for more details.

If you don't want to introduce the dependency or have some other constraint (e.g., OpenMP can be problematic in library code), you can always do it yourself by having the inner for loop in a thread, kicking off N threads at the start and joining at the end. That's of course assuming you have enough work, which seems like you might if the matrices are big enough.

user_379290

Not sure if this will be much help but I do have an old example of using the pthread.h library to compute Gauss elimination with partial pivoting matrices.

In summary, the highlights are:

  • Create an array of threads pthread_t threads[N];
  • initlize the halting barrier for threads to run to pthread_barrier_init(&barrier, NULL, numThreads);
  • Set your barriers in the function you're trying to multithread so that it will wait until each function has the dependencies needed to continue. Add pthread_barrier_wait(&barrier); at your points
  • Start your threads
    for (i = 0; i < nthreads; i++)
    {
        pthread_create(&threads[i], NULL, functionWithThreading, (void *)i);
    }
  • Finally, wait for all the threads to finish and join them up
    for (i = 0; i < nthreads; i++)
    {
        pthread_join(threads[i], NULL);
    }

I know this might not be the exact solution you're looking for but I hope the example might help

user_123798840
Related