What mapping strategy should I use when solving many large linear systems of equations?

Question

I am working on a problem that involves solving many (thousands) of distinct linear systems of equations, each with thousands of variables. Let's assume that the size of each matrix is exactly the same, and has the same structure (i.e. same symmetry, sparsity, diagonal dominance, etc).

Assuming each linear system is solved by the same method, and each linear system is independent, I foresee three ways to approach solving all linear systems:

1. I can solve each linear system with all the processors at my disposal, one at a time.

2. I could allocate a certain number of linear systems to each processor, then solve each one sequentially.

3. A hybrid of the first two options... That is, if we have P processors, we can designate G groups of processors (assume P is a multiple of G for simplicity). Then, for N linear systems of equations, we allocate N/G systems to each group (assume N is a multiple of G for simplicity). Then each group can solve its N/G systems one by one in parallel using its P/G processors simultaneously.

I believe the mapping strategy that yields the solutions of all systems of equations in the fastest time depends both upon the size of each system, and the number of systems. Which strategy works best:

1. If the number of systems remains constant, but the size each matrix system increases?

2. If the size of each matrix system remains constant, but the number of systems increases?

Answer 1

Communication between nodes is (usually) much slower than accessing RAM local to that node. Strategy 2 is probably your best bet, unless you're systems are so big they can't fit into memory on one node any more.

Because communication between nodes is usually faster than communication with a hard disk, you should split your problem up between multiple nodes (strategy 3) to keep the whole thing in RAM if it won't fit on a single processor. If each problem is so large you need all your available computational power to fit it into memory, then use strategy 1.

To answer your questions:

1. If you're dealing with matrices big enough to fill up all of the memory on one node (you'll need more than "thousands of variables" to do this), go with either strategy 1 or 3, depending on how big your systems are. If they're not big enough to fill up memory, go with strategy 2.

2. If the matrices can all fit into main memory, go with strategy 2 regardless of how many systems you have.

Answer 2

I'd go for option two with a hybrid shared/distributed memory twist... I'm assuming that the 13 machines in your cluster all access the same shared file system and that each computer has enough memory to load 8 problems.

Your code could then look something like this:

int fd, next_task = 0, my_task, nr_tasks;

/* Open the file you want to use to synchronize the nodes. */
if ( ( fd = open( "sharedlockfile" , O_CREAT ) ) < 0 ) {
    printf( "failed to create lock file.\n" );
    abort();
    }

/* Initialize the curr_task in the file. */
if ( <I am the first node in this cluster!> )
    if ( write( fd , &next_task , sizeof(int) ) < 0 ) {
        printf( "failed to write to lock file." );
        abort();
        }

/* Main parallel loop. */
#pragma omp parallel shared(fd,next_task), private(mytask)
while ( next_task < nr_tasks ) {

    /* Try to get a hold of the lock file, grab a task, and write the new counter. */
    #pragma omp critical
    { 

        /* Get a hold of the file lock (this will block until the file is free). */
        if ( flock( fd , LOCK_EX ) < 0 ) {
            printf( "file locking failed.\n" );
            abort(); // not sure I'm allowed to abort in a parallel block...
            }

        /* Read the index of the next task. */
        if ( lseek( fd , 0 , SEEK_SET ) < 0 || read( fd , &next_task , sizeof(int) ) < 0 ) {
            printf( "error reading next task id.\n" );
            abort();
            }

        /* Remember my task and update the counter. */
        my_task = next_task;
        next_task += 1;

        /* Write the index of the next task back to the file. */
        if ( lseek( fd , 0 , SEEK_SET ) < 0 || write( fd , &next_task , sizeof(int) ) < 0 ) {
            printf( "error writing next task id.\n" );
            abort();
            }

        /* Let go of the file. */
        if ( flock( fd , LOCK_UN ) < 0 ) {
            printf( "file un-locking failed.\n" );
            abort();
            }

        } /* end of critical section */

    /* Did we get a valid task ID? */
    if ( my_task < nr_tasks ) {

        /*
         * Do whatever this task implies.
         */

        } /* check if valid task ID. */

    } /* Main loop. */

Please note that I haven't tested the above code! I've done similar things before and most of this is off the top of my head... Check out the manpage for flock to be sure.

If all the nodes in your cluster execute this, they should each spawn as many threads as they have cores (or as many as you specified in OMP_NUM_THREADS) and each core will try to acquire a task id. Access to next_task is controlled on two levels:

• By the #pragma omp critical so that no two threads on the same node try to mess with the next_task variable at the same time.
• By the exclusive file lock so that no two nodes read/update the next_task variable stored in the file at the same time.

This should give you a good, adaptive, dynamic scheduling of your computations.

Edit

Uhm, may not be as easy as I thought with the locking as it seems to be fd-dependent and thus may not work accross nodes. There seems to be a good solution here though.

Addendum

If solving each system on a separate core is too much of a strain for your memory bus, e.g. you notice that it's not scaling well at all, you might want to try MAGMA, which is like LAPACK, but scales very nicely on multicore architectures.