Efficient computation of kronecker products in C

Question

I'm fairly new to C, not having much need to anything faster than python for most of my research. However, it turns out that recent work I've been doing required the computation of fairly large vectors/matrices, and there therefore a C+MPI solution might be in order.

Mathematically speaking, the task is very simple. I have a lot of vectors of dimensionality ~40k and wish to compute the Kronecker Product of selected pairs of these vectors, and then sum these kronecker products.

The question is, how to do this efficiently? Is there anything wrong with the following structure of code, using for loops, or obtain the effect?

The function kron described below passes vectors A and B of lengths vector_size, and computes their kronecker product, which it stores in C, a vector_size*vector_size matrix.

void kron(int *A, int *B, int *C, int vector_size) {

    int i,j;

    for(i = 0; i < vector_size; i++) {
        for (j = 0; j < vector_size; j++) {
            C[i*vector_size+j] = A[i] * B[j];
        }
    }
    return;
}

This seems fine to me, and certainly (if I've not made some silly syntax error) produce the right result, but I have a sneaking suspicion that embedded for loops is not optimal. If there's another way I should be going about this, please let me know. Suggestions welcome.

I thank you for you patience and any advice you may have. Once again, I'm very inexperienced with C, but Googling around has brought me little joy for this query.

Answer 1

Since your loop bodies are all completely independent, there is certainly a way to accelerate this. Easiest would be already to take advantage of several cores before thinking of MPI. OpenMP should do quite fine on this.

#pragma omp parallel for
for(int i = 0; i < vector_size; i++) {
    for (int j = 0; j < vector_size; j++) {
        C[i][j] = A[i] * B[j];
    }
}

This is supported by many compilers nowadays.

You could also try to drag some common expressions out of the inner loop but decent compilers e.g gcc, icc or clang should do this quite well all by themselves:

#pragma omp parallel for
for(int i = 0; i < vector_size; ++i) {
    int const x = A[i];
    int * vec = &C[i][0];
    for (int j = 0; j < vector_size; ++j) {
        vec[j] = x * B[j];
    }
}

BTW, indexing with int is usually not the right thing to do. size_t is the correct typedef for everything that has to do with indexing and sizes of objects.

Answer 2

For double-precision vectors (single-precision and complex are similar), you can use the BLAS routine DGER (rank-one update) or similar to do the products one-at-a-time, since they are all on vectors. How many vectors are you multiplying? Remember that adding a bunch of vector outer products (which you can treat the Kronecker products as) ends up as a matrix-matrix multiplication, which BLAS's DGEMM can handle efficiently. You might need to write your own routines if you truly need integer operations, though.

Answer 3

If your compiler supports C99 (and you never pass the same vector as A and B), consider compiling in a C99-supporting mode and changing your function signature to:

void kron(int * restrict A, int * restrict B, int * restrict C, int vector_size);

The restrict keyword promises the compiler that the arrays pointed to by A, B and C do not alias (overlap). With your code as written, the compiler must re-load A[i] on every execution of the inner loop, because it must be conservative and assume that your stores to C[] can modify values in A[]. Under restrict, the compiler can assume that this will not happen.

Answer 4

Solution found (thanks to @Jeremiah Willcock): GSL's BLAS bindings seem to do the trick beautifully. If we're progressively selecting pairs of vectors A and B and adding them to some 'running total' vector/matrix C, the following modified version of the above kron function

void kronadd(int *A, int *B, int *C, int vector_size, int alpha) {

    int i,j;

    for(i = 0; i < vector_size; i++) {
        for (j = 0; j < vector_size; j++) {
            C[i*vector_size+j] = alpha * A[i] * B[j];
        }
    }
    return;
}

precisely corresponds to the BLAS DGER function (accessible as gsl_blas_dger), functionally speaking. The initial kron function is DGER with alpha = 0 and C being an uninitialised (zeroed) matrix/vector of the correct dimensionality.

It turns out, it might well be easier to simply use python bindings for these libraries, in the end. However, I think I've learned a lot while trying to figure this stuff out. There are some more helpful suggestions in the other responses, do check them out if you have the same sort of problem to deal with. Thanks everyone!