flintlib / arb

Arb has been merged into FLINT -- use https://github.com/flintlib/flint/ instead
http://arblib.org/
GNU Lesser General Public License v2.1
455 stars 137 forks source link

Is it possible to index acb_t matrices? #235

Open rudolph-git-acc opened 6 years ago

rudolph-git-acc commented 6 years ago

An acb_mat type matrix (say 100x100 cells) is the output of a summation done in a long (n=1..10^12) loop. To speed up the processing, the loop is split up in equal 'chunks' based on the number of threads chosen and these threads are then run in parallel. This all works fine and delivers the expected speed gains.

There is a simple final task to add up all the individualacb_mat matrices that each thread produces. This is easy to do in a 'hard-coded' manner by assigning a unique name to each thread output matrix (using the thread-id), but this is not very elegant nor flexible for scaling up. Is there a way to 'index' an acb_mat matrix, so the final addition process of the sub-matrices could be properly coded (or any other idea on how this final step could be done)?

fredrik-johansson commented 6 years ago

I don't understand the problem.

rudolph-git-acc commented 6 years ago

Let me try to explain differently.

Question in a nutshell: It is allowed to put arb_poly_t or acb_poly_t variables (and others) in an array using a ptr. I need to make a similar array of arb_mat or acb_mat matrices. Is this possible and/or is there a workaround?

Question in the context of the application: A long loop 1...10^12 is split into various smaller loops each starting at startand ending at stop (these are parameters passed through to the thread). Each thread produces a unique acb_matrix as output. This is the simplified thread code:

void* storedsums_calc(void *voidData)
{
    struct ThreadData* data=voidData;

    acb_mat_t resacb;

    //recover the data passed to this specific thread
    slong expterms=data->expterms;
    slong taylorterms=data->taylorterms;
    slong prec=data->prec;
    slong start=data->start;
    slong stop=data->stop;
    slong id=data->threadid;

    acb_mat_init(resacb, expterms, taylorterms);

    //produce the stored sums matrix(expterms, taylorterms)
    for (n = start; n < stop; n++)
    {
        perform all calculations to fill the 'resacb' matrix for this n
        establish 'nexpo'

       //this is the current way I make a ' finmat_x' output matrix unique for each thread (3 in this case)
        if (threadid == 0)
            acb_mat_scalar_addmul_acb(finmat0, resacb, nexpo, prec);
        if (threadid == 1)
            acb_mat_scalar_addmul_acb(finmat1, resacb, nexpo, prec);
        if (threadid == 2)
            acb_mat_scalar_addmul_acb(finmat2, resacb, nexpo, prec);
        etc.
    }
    return(NULL);
}

and then I add all output matrices from the threads together like this:

(...)
    //wait for all threads to complete
    for (i = 0; i < numthreads; i++)
    {
    pthread_join(thread[i], NULL);
    }

    //sum all thread outputs (stored in globally defined matrices) together
    acb_mat_add(finalmat, finmat0, finmat1, prec);
    acb_mat_add(finalmat, finalmat, finmat2, prec);
    etc.
(...)

This code works for three threads, but I'd like to make this dependent on any number of threads chosen (numthreads). Therefore a finmat matrix that is stored in an array, indexed by threadid, is what I'm after. Is this possible in ARB? If not, is there maybe an alternative way to do it?

fredrik-johansson commented 6 years ago

You can create arrays of acb_mat_t's the same way you create arrays of any other type in C.

acb_mat_t * vec;
vec = malloc(sizeof(acb_mat_t) * len);
/* or: acb_mat_t vec[len]; if len is a fixed small number */

for (i = 0; i < len; i++)
    acb_mat_init(vec[i], n, n);

for (i = 0; i < len; i++)
    acb_mat_add(vec[i], ...);

for (i = 0; i < len; i++)
    acb_mat_clear(vec[i]);
free(vec);

This does not require special library support.