Submit parallel mkl openmp jobs (blupf90, yams) in genotoul cluster

The blupf90 software and my own one includes parallelised code using Intel fortran MKL and OPENMP facilities. These use several threads, usually witihin a single node (as they share memory). I was getting nuts to make this work in bioinfo.genotoul.fr until I went to see the Support people. One of the problems for me was that there are actually several grid engines so google does not always find the good solution.

Short message: use this:

alegarra@genotoul2 ~/work $ cat test_omp.sh

#!/bin/bash

#$ -pe parallel_smp 10

export OMP_NUM_THREADS=10

your_program

alegarra@genotoul2 ~/work $ qsub test_omp.sh

where 10 is the number of threads that you want.

Long explanation:

I have a simple example that just creates a big matrix and inverts it. Let's try a 1000 x 1000 matrix, and the program runs smoothly without any particular command or request for parallel:

$cat test_omp.sh

#!/bin/bash

echo "1000 1000" | /save/alegarra/progs/parallel_inverse/a.out

alegarra@genotoul2 ~/work $ qsub test_omp.sh

$cat test_omp.sh.o8160505
nanim,nsnpp

        1000        1000

 Z=rnd()

 GG=ZZ

    Dgemm MKL #threads=    20   40 Elapsed omp_get_time:     0.1530

 GG=GG/(tr/n)

 GG=0.95GG+0.05I

 GG^-1

    Inverse LAPACK MKL dpotrf/i #threads=   40   40 Elapsed omp_get_time:     0.1071

Epilog : job finished at Thu Mar  9 15:45:53 CET 2017

However, this fails if I try a larger 10000 x 10000 matrix:

OMP: Error #34: System unable to allocate necessary resources for OMP thread:

OMP: System error #11: Resource temporarily unavailable

OMP: Hint: Try decreasing the value of OMP_NUM_THREADS.

OMP: Error #178: Function pthread_getattr_np failed:

OMP: System error #12: Cannot allocate memory

forrtl: error (76): Abort trap signal

The reason is that every thread allocates memory, so that this is too much. 

Then I try by requesting explicitly 4 threads and 1000 x 1000 matrix. This is doable by requesting the adequate environment:

-pe parallel_smp demands X cores on the same node (multi-thread, OpenMP) as explained in the docs.

Do not forget putting the SAME number in -pe parallel_smp and in export OMP_NUM_THREADS

alegarra@genotoul2 ~/work $ cat test_omp.sh

#!/bin/bash

export OMP_NUM_THREADS=4 

echo "1000 1000" | /save/alegarra/progs/parallel_inverse/a.out

alegarra@genotoul2 ~/work $ qsub -pe parallel_smp 4 test_omp.sh

and it does work:

alegarra@genotoul2 ~/work $ cat test_omp.sh.o8160541 

 nanim,nsnpp

        1000        1000

 Z=rnd()

 GG=ZZ

    Dgemm MKL #threads=     4    4 Elapsed omp_get_time:     0.0320

 GG=GG/(tr/n)

 GG=0.95GG+0.05I

 GG^-1

    Inverse LAPACK MKL dpotrf/i #threads=    4    4 Elapsed omp_get_time:     0.0230

Epilog : job finished at Thu Mar  9 16:01:30 CET 2017

now, try with 10000 x 10000:

alegarra@genotoul2 ~/work $ cat test_omp.sh.o8160558 

 nanim,nsnpp

       10000       10000

 Z=rnd()

 GG=ZZ

    Dgemm MKL #threads=     4    4 Elapsed omp_get_time:    45.8333

 GG=GG/(tr/n)

 GG=0.95GG+0.05I

 GG^-1

    Inverse LAPACK MKL dpotrf/i #threads=    4    4 Elapsed omp_get_time:    29.2803

Epilog : job finished at Thu Mar  9 16:05:11 CET 2017

and it also works. I can also put the -pe within the script:

alegarra@genotoul2 ~/work $ cat test_omp.sh

#!/bin/bash

#$ -pe parallel_smp 4

export OMP_NUM_THREADS=4

echo "10000 10000" | /save/alegarra/progs/parallel_inverse/a.out

alegarra@genotoul2 ~/work $ qsub test_omp.sh

Now, I try with 10 threads and it also works. Then I will try using yams:


A couple of tricks more. To know the parallel environments in the cluster and some characteristics:

alegarra@genotoul2 ~/work/mtr_2015/mf/methodR $ qconf -spl

parallel_10

parallel_20

parallel_40

parallel_48

parallel_fill

parallel_fill_amd

parallel_rr

parallel_rr_amd

parallel_smp

parallel_testq

alegarra@genotoul2 ~/work/mtr_2015/mf/methodR $ qconf -sp parallel_smp

pe_name            parallel_smp

slots              9999

user_lists         NONE

xuser_lists        NONE

start_proc_args    NONE

stop_proc_args     NONE

allocation_rule    $pe_slots

control_slaves     TRUE

job_is_first_task  FALSE

urgency_slots      min

accounting_summary TRUE

Julia syntax in vim

I try to put julia syntax in Macvim. I found this github repository: https://github.com/JuliaEditorSupport/julia-vim and I blindly follow the instructions in https://github.com/JuliaEditorSupport/julia-vim#manually , starting from ~/.vim . It does work:

There are some cool tricks, for instance I can write Unicode :
\alpha + Tab = α

tests with Julia

So, after having installed Julia months ago I gave it some tries. The speedup seems larger than an order of magnitude in both cases

looping through the same matrix

R

system.time({val=0; for (i in 1:2000){for (j in 1:2000){val=val+a[i,j]}}; val})
   user  system elapsed 
  1.862   0.015   1.868 

Julia

julia> tic(); val=0; for i in 1:2000; for j in 1:2000; val=val+a[i,j]; end; end; println(val); toc()
5.142821063269366
elapsed time: 0.327145074 seconds
0.327145074




Singular value decomposition of a 2000 x 2000 matrix

> a=matrix(rnorm(4000000),2000)
> system.time(svd(a))
   user  system elapsed 
 27.419   0.218  27.452 

Julia

julia> tic(); a=randn(2000,2000); toc()
elapsed time: 0.673701438 seconds
0.673701438

julia> tic(); svd(a); toc()
elapsed time: 5.236212822 seconds
5.236212822

transposing marker data

Quite frequently one may find marker data ordered in this way:

alegarra@genotoul2 ~/save/progs $ cat ex_rachel

 1 a b

 1 b b

 1 c c

 2 b b

 2 c c

 2 d d

 3 a b

 3 b b

 3 c c

e.g. there are three animals (1 to 3) and three markers. Rachel and I would like them to be formatted one animal per line and markers one after each other, in this way:

alegarra@genotoul2 ~/save/progs $ cat out

         1 ab bb cc 

         2 bb cc dd 

         3 ab bb cc 

This is conceptually simple if animals are sorted:

1. Read a line
2. If the animal is the same as the old one, print the markers after the previous one
3. If the animal is different, start a new line, print the animal and the markers.
4. Add special cases for the first and last line

Here is an awk implementation

alegarra@genotoul2 ~/save/progs $ cat SNPcol2line.awk 

#! /bin/awk -f

# this program reads genotypes in one line per locus

# then puts them as

# individual allele1 allele2 allele1 allele 2

#

BEGIN{

idold=0

i=0

}

{

id=$1

# if new animal

if(id!=idold){

if(idold!=0){

# close previous line

printf("\n")

}

# write new id

printf("%10s%1s",id," ")

idold=id

}

printf("%1s%1s%1s",$2,$3," ")

}

END{

# last individual

printf("\n")

}

which works:

alegarra@genotoul2 ~/save/progs $ ./SNPcol2line.awk ex_rachel 

         1 ab bb cc 

         2 bb cc dd 

         3 ab bb cc 

Sires in validation and Sires in training

I want to test metafounders for genomic selection. So I have these MTR sires that have daughters in the validation , in the training, or in both of them, and I want a file with these numbers. I have two files:

T is:
juan 10
pepe 5

V is:

luis 6
pepe 4

so luis has no daughters on Training, and I want this to appear. Then I use this pipeline from Toni Reverter from CSIRO:

awk '{print $1}' $1 $2 | sort -u | join -a1 - $1 | awk '{print $1, (NF==2?$2:0)}' | join -a1 - $2 | awk '{print $1, $2, (NF==3?$3:0)}'

which gives

putterri:yarp andres$ ./joinToni.sh V T

juan 0 10

luis 6 0

pepe 4 5

Install GNU coreutils in the Mac

The MacOSX join is poorer than GNU one. So I try toinstall it using macports. First I have to update the data base of port:

putterri:mf andres$ port -v selfupdate

But this does not work. Then I try

putterri:mf andres$ sudo port -d selfupdate

This does work. Join is actually part of coreutils:

putterri:mf andres$ port search --name coreutils

coreutils @8.25 (sysutils)

    GNU File, Shell, and Text utilities

xml-coreutils @0.8.1_1 (textproc, xml)

    Command line tools for XML processing

Found 2 ports.

So I do install coreutils:

putterri:mf andres$ sudo port install coreutils

The most important information is at the end:

The tools provided by GNU coreutils are prefixed with the character 'g' by default to distinguish them from the BSD commands.

For example, cp becomes gcp and ls becomes gls.

If you want to use the GNU tools by default, add this directory to the front of your PATH environment variable:

    /opt/local/libexec/gnubin/

Plotting correlations

I found this nice R package to plot correlations : corrplot . Imagine that you want to present five genetic correlations to your buddies. Here's how to visualize them nicely.

#create a covariance matrix
set.seed(1234)
V=rWishart(1,10,diag(5))[,,1]
cov2cor(V)
            [,1]       [,2]       [,3]       [,4]        [,5]
[1,]  1.00000000  0.1169636  0.2236465 -0.2667448  0.04351447
[2,]  0.11696358  1.0000000 -0.2260988 -0.3425368  0.64837326
[3,]  0.22364648 -0.2260988  1.0000000 -0.1883763 -0.22484960
[4,] -0.26674477 -0.3425368 -0.1883763  1.0000000 -0.50315007
[5,]  0.04351447  0.6483733 -0.2248496 -0.5031501  1.00000000
#give it names
colnames(V)=c("MY","gain","longevity","SCS","pietin")
rownames(V)=c("MY","gain","longevity","SCS","pietin")

#plot
require(corrplot)
corrplot.mixed(cor(V),col=gray.colors(10))