MPI and lapack

I've been upgrading my code to use mpi and since I use lapack in C for the diagonalization and some other procedures I needed to make a doublecomplex type for mpi and I didn't find anyone else who has posted this. It's not hard but when you are starting out with mpi one less thing to worry about is nice so I'll post the relevant few lines.

MPI_Datatype MPI_DOUBLE_COMPLEX;
MPI_Type_contiguous(2, MPI_DOUBLE, &MPI_DOUBLE_COMPLEX);
MPI_Type_commit(&MPI_DOUBLE_COMPLEX);
MPI_Type_free(&MPI_DOUBLE_COMPLEX);

This works by making MPI_DOUBLE_COMPLEX from two contiguous doubles since doublecomplex is a struct with two doubles r and i.

More details on openmp

There are a few little things I've learned that I want to share in hopes others using the openmp will not spend the 6 hours "debugging" the same types of mistakes I did.
(This is meant to clear up some stuff that I found unclear on other pages but is in no way complete).

The first thing to try is compiling without openmp (usually by not including the -openmp flag) and see if there are any normal mistakes. If the code is mistakeless (no hanging parenthesis or undeclared variables) then it is still compilable without openmp it will just run as a serial program.

The next thing is that directs for improving loops is incredibly sensitive and they really mean the loop must be on the very line. For instance
#pragma omp parallel for private(i)
norm =1; for(i=0 ;i <N ; i++) {...}

gave me errors of: syntax error before '{' token but on a line that was 10 or so lines preceding this code. That took me a long time to figure out. It doesn't complain about the code itself but says there are errors in parenthesis in preceding lines (at least that's what happened for me). But just putting the norm = 1; before the #pragma directive it all works. Putting parenthesis encapsulating the for loop also seem to be a source of error. So it wants
#pragma omp parallel for private(i)
for(i=0 ;i <N ; i++) {...}
NOT
#pragma omp parallel for private(i)
{for(i=0 ;i <N ; i++) {...} }

Also there seems to be no reason to do everything under a giant #pragma omp parallel code segment with subsections for #pragma omp for or #pragma omp sections. In fact it seems much easier to do it in small forks. As in having a segment:
...
#pragma omp parallel sections
{
#pragma omp section
{ code; code(A,B); etc... }
#pragma omp section
{ code; code(B,A); etc... }
}
...
printf(...); other regular code...
...
#pragma omp parallel for schedule(static) private (i,j)
for(i=0 ;i <N ; i++) {...}

Well I hope this save someone time and effort.

Threading a scientific program

I wanted one of my codes to have a section where it would diagonalize 4 different matrices at the same time (since this is my bottleneck) and the computers that I'm running on (sharcnet) are mostly quad processors. I looked in to MPI and it's apparent processing and it is hard!! It is clearly not meant for what I wanted to do. I have a code that executes in serial and at a bottleneck point I wanted it to do many things on different CPU's. Then I found openMP. This was what I wanted, a way to parallelize a section of my code without a complete rewrite! I found the SECTIONS and SECTION commands especially easy. For instance say you have 2 matrices A and B and you want to do something to them with the function diag(). In a normal C code you'd have
...
diag(A); //long wait
diag(B); //also long wait
...

With openmp you can fork this section of the code with
#pragma omp parallel sections
{ }
This created an area that is forked and then joined at the end. So once this region finishes executing it goes back to serial-type execution. In the SECTIONS part you now (optimally) make a number of SECTIONs that equals the number of processors (not rquired but makes sense) and each section will execute on a different CPU at the same time but once the SECTIONS section ends it's regular execution again. So for matrices A and B you'd have
...
#pragma omp parallel sections
{
#pragma omp section
{ diag(A)}

#pragma omp section
{ diag(B)}
}
... //this will be normal C or C++ code. There is plenty more but I don't need it right now so that's all I know. I can say this is defiantly easier than using MPI send/receive! Some good sites are AMD's short one and a more complete one here. Hope this helps someone but leave me some CPU's please. ;-)