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

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?

• What kind of system are you running this on? – Dan Feb 6 '12 at 23:42
• Are these shared-memory or distributed-memory processors? – Pedro Feb 6 '12 at 23:51
• I would have to say that it is a hybrid cluster: 13 distributed nodes, each with 8 cores sharing memory. – Paul Feb 7 '12 at 0:04
• @Paul: I have recently been faced with a similar question and the answer heavily depends on the scalability characteristics of your solver, and if/when you will run out of memory. Would you mind posting which solver you will be using and what your ideal problem size will be? Also, since you have asked an efficiency question, I would recommend looking into solving several right-hand sides in batches rather than independently, as it lowers the required number of memops. – Jack Poulson Feb 7 '12 at 16:43
• @JackPoulson: Since the linear systems are sparse, symmetric, diagonally dominant, and positive definite, I figured a conjugate gradient method would work very quickly. But I'm not sure what you mean by "solving several RHS's in batches"... could you elaborate on that? – Paul Feb 7 '12 at 16:51

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.

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.

• are there situations in which strategy 3 might be faster? – Paul Feb 6 '12 at 23:33
• Sure: when the problem gets too big to fit on one node. You probably shouldn't use strategy 1 unless need all (or at least most) of the nodes in order to fit the problem in RAM. – Dan Feb 6 '12 at 23:37
• @Paul: I suppose 3 is also a good idea when you have more processors than systems to solve. – Dan Feb 6 '12 at 23:40
• @Dan: Assuming there are more right-hand sides than processes, I completely agree; if he can perform approach two without running out of memory, it is always ideal. Otherwise, he should relax to approach three, and approach one would be the extreme case for very large systems. – Jack Poulson Feb 7 '12 at 17:26
• @Dan: I usually think of superlinear scaling as the result of being able to move the local computation higher up in the memory hierarchy. As an extreme example, suppose that on one process the problem cannot fit in memory and has to be solved out-of core: it is then very reasonable to get a superlinear speedup as soon as the problem can be solved in core. – Jack Poulson Feb 7 '12 at 17:45

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. */
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. */

/* 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.