Suyeon Son collects some great tips about Googling for answers to programming questions:

I asked Jeremy Bowers, a news applications developer at NPR, who said he follows a specific pattern when Googling to achieve specific outcomes:

e.g. “JavaScript remove key from object”

I would also recommend Googling the exact string of an error message (minus any specific variable names that might appear in it). This is particularly useful for C compile errors.

Also, note the great reference to the original 2000 email to Yak shaving!

I have created a LinkedIn group for IDL Users. If you are on LinkedIn, look it up and invite other IDL users.

QR factorization of a matrix A is the process of determining an orthogonal matrix Q and upper triangular matrix R such that $$A = QR$$ Using a QR factorization is faster and much more numerically stable than finding an inverse matrix of A to solve the system of equations Ax = b. QR factorization is used in fundamental linear algebra IDL routines such as LA_EIGENPROBLEM, LA_LEAST_SQUARES, LA_SVD, and others. I was recently using MAGMA’s [SDCZ]GEQRF within GPULib to retrieve a QR factorization and found excellent performance improvements by using the GPU.

Also, the full decomposition into Q and R matrices is seldom needed, so the results are not returned in a naive manner. I will show how to retrieve Q and R for the GPULib routines (it’s slightly different than the CPU version).

Performance of MAGMA for QR factorization was excellent and it scaled quite well compared to the CPU as shown below.

For small matrices of size 500 x 500, I got about a 25x speedup, but on larger 5000 x 5000 element matrices about a 435x speedup.^[1]

Next, we’ll go through how to do perform QR factorization using GPULib and verify that we are really factoring A correctly. First, as always, we initialize GPULib:

IDL> gpuinit
GPULib Full 1.7.0 (Revision: 2876M)
Graphics card: Tesla C2070, compute capability: 2.0, memory: 1205 MB available, 1279 MB total
CUDA version: 5.5
MAGMA version: 1.4.0
Checking GPU memory allocation...cudaSuccess

Next, we make a simple test case and transfer it to the GPU:

IDL> a = [[9.0, 2.0, 6.0], [4.0, 8.0, 7.0]]
IDL> da = gpuputarr(a)

There are a few variables we need to setup:

IDL> dims = size(da, /dimensions)
IDL> m = dims[0]
IDL> n = dims[1]

The GPU routine we will use, GPUSGEQRF, is a call directly into the MAGMA C layer, so we need to setup a GPU workspace variable of the correct size as well as create a CPU variable to return the multipliers in tau:

IDL> info = 0L
IDL> tau = fltarr(m < n)
IDL> nb = gpu_get_sgeqrf_nb(m)
IDL> lwork = (2L * (m < n) + (n + 31L) /32L * 32L ) * nb
IDL> dworkspace = gpufltarr(lwork)

We are ready to call the routine to perform the QR factorization:

IDL> status = gpusgeqrf(m, n, da->_getHandle(), m, tau, dworkspace->_getHandle(), info)

Part of the result is done in place in da (the other portion of the result is in the CPU variable tau):

IDL> result = gpugetarr(da)

Q is decomposed into H reflector matrices such that: $$Q = H_1 H_2 . . . H_k$$ where k = min(m, n), $$H_i = I – \tau_i v v^T$$ and v is stored in da:

v[i + 1:m - 1L] = result[i + 1:m - 1L, i]

The code to reconstruct Q is therefore:

IDL> q = identity(m)
IDL> .run
- for i = 0L, n - 1L do begin
- v = fltarr(m)
- v[i] = 1.0
- v[i + 1:m - 1L] = result[i + 1:m - 1L, i]
- q_sub = (identity(m) - tau[i] * (v # v))
- q = q ## q_sub
- endfor
-
- end

R is stored in the other triangle of da:

IDL> r = fltarr(m, n)
IDL> for i = 0L, n - 1L do r[0L:i, i] = result[0L:i, i]

To verify the result, r is upper triangular (IDL and standard mathematical notation have different conventions for displaying matrices):

IDL> print, r
-11.0000 0.00000 0.00000
-8.54546 -7.48166 0.00000

Q is orthogonal:

IDL> print, transpose(q) # q
1.00000 0.00000 5.96046e-08
0.00000 1.00000 -2.98023e-08
5.96046e-08 -2.98023e-08 1.00000

And, $$A = QR$$

IDL> print, transpose(q) # r<br />9.00000 2.00000 6.00000<br />4.00000 8.00000 7.00000

Full disclosure: I work for Tech-X and I am the product manager for the GPULib.

Jake Vanderplas, writing on AstroBetter, discusses the licensing of scientific code, providing the following suggestions:

1. Always license your code. Unlicensed code is closed code, so any open license is better than none (but see #2).
2. Always use a GPL-compatible license. GPL-compatible licenses ensure broad compatibility for your code, and include GPL, new BSD, MIT, and others (but see #3).
3. Always use a permissive, BSD-style license. A permissive license such as new BSD or MIT is preferable to a copyleft license such as GPL or LGPL.

I agree. IDLdoc, mgunit, rIDL, my library, and all the other open source code I have on GitHub uses a BSD-style license.

As part of the process of getting IDL recognized by GitHub, I had to add IDL support to Pygments, a popular syntax highlighter supporting many languages. This allows generating syntax highlighting of IDL code.

Use pip to install Pygments:

$ pip install Pygments

Then use pygmentize to generate output:

$ pygmentize -O full,style=colorful -f html -o mg_repmat.html mg_repmat.pro

This produces the output below:

Here is the full output. This is obviously not perfect (I just wanted GitHub to recognize the language as IDL).

mpiDL is a set of IDL bindings for the Message Passing Interface (MPI). It is used for tasks where communication between processes is required, as opposed to the independent behavior of TaskDL workers. It can make use of the multiple cores of a single computer and/or multiple nodes of a cluster. mpiDL is supported on OS X and Linux for both OpenMPI and MPICH.

As an example of using mpiDL, I will present a simple probability-based computation of pi using many cores of a computer. If you are interesting in evaluating mpiDL or require more information about, please contact me.

Setting up for using mpiDL requires placing the mpiDL lib/ directory in your IDL path and DLM path, as well as setting the MPIDL_DIR environment variable to the root of your mpiDL distribution (make sure to do it in .bashrc to make it available to non-interactive processes). Then you are ready to run an mpiDL program, such as the parallel_pi example provided in the distribution:

examples$ runmpidl -np 8 ${PWD}/parallel_pi.sav
Running mpiexec -np 8 /home/research/mgalloy/software/mpiDL-2.4.0-Linux64/bin/mpidlstart
/home/research/mgalloy/software/mpidl-r431-par/examples/parallel_pi.sav...
******************************************************************
MPIDL Version 2.4.0 - a parallel implementation of IDL
(C) Copyright 2000 - 2014, Tech-X Corp.
All rights reserved.
******************************************************************
Process 1 gave: 3.143600
Process 2 gave: 3.131600
Process 3 gave: 3.137200
Process 4 gave: 3.116400
Process 5 gave: 3.140000
Process 6 gave: 3.172800
Process 7 gave: 3.157200
The final value of pi is: 3.142686

The -np 8 argument indicates that the example should run with 8 processes; the parallel_pi example code has decided that it will have one master process collecting results and 7 workers computing pi. The example can also run with a single process, where that process both computes pi and collects the result:

examples$ runmpidl -np 1 ${PWD}/parallel_pi.sav
Running mpiexec -np 1 /home/research/mgalloy/software/mpiDL-2.4.0-Linux64/bin/mpidlstart
/home/research/mgalloy/software/mpidl-r433-434M-par/examples/parallel_pi.sav...
Process 0 gave: 3.142200
The final value of pi is: 3.142200

Also, note that a .sav file is specified. This allows an each process to use a runtime IDL license instead of a full development license.

There are two routines in the source code: parallel_pi_calcpi, which does not know about MPI and just computes an estimate of pi, and parallel_pi itself which coordinates the work using the MPI interface. By the nature of MPI, the same routine is executed in each process, but the processes are given an identifier, called the “rank”, which lets the routine decide what it should be doing. In parallel_pi, we determine that with the following code:

rank = mpidl_comm_rank()
nprocs = mpidl_comm_size()

Here, nprocs is the total number of processes and rank is an identifier from 0 to nprocs - 1. From this information, parallel_pi can determine if it will be computing pi and sending the results back or if it is the receiver of the information:

am_a_sender = (nprocs eq 1L) or (rank gt 0L)
am_a_receiver = (nprocs eq 1L) or (rank eq 0L)

This is more complicated than might seem to be required in order to handle the case of only one process which would need to be both a sender and receiver.

If it is a sender, the process must compute pi and send it back to the receiver (the rank 0 process referenced by DEST in the mpidl_send routine below):

if (am_a_sender) then begin
  seedr = rank
  a = dblarr(1)
  a[0] = parallel_pi_calcpi(neval, seedr)
  mpidl_send, a, DEST=0
endif

If it is the receiver, it allocates an array to hold the result from each sender and then receives the result from each sender in turn with mpidl_recv. Finally, it computes the average of the values to give the result:

if (am_a_receiver) then begin
  mypi = dblarr(n_senders)
  for j = 0L, n_senders - 1L do begin
    mypi[j] = mpidl_recv(COUNT=1, SOURCE=senders[j], /DOUBLE)
    print, senders[j], mypi[j], format='(%"Process %d gave: %f")'
  endfor
  print, total(mypi, /preserve_type) / n_senders, $
         format='(%"The final value of pi is: %f")'
endif

Check out the source code for all the details. Also, download the users guide for more information.

Full disclosure: I work for Tech-X and I am the product manager for the FastDL suite which includes mpiDL.

There are IDL-specific editing modes for many of the more popular text editors:

IDLwave — This is the oldest and most complete mode for editing IDL code. JD Smith has been maintaining this for years. JD was looking for a new maintainer a while ago and moved IDLwave from its own website to GitHub, but everything seems to be still working.

BBEdit/TextWrangler — This mode includes a routine to scan IDL code and generate tags used for autocompletion.

Vim mode — Marshall Perrin has a Vim mode supporting IDL.

TextMate 2 — I used TextMate 2 on the Mac and decided to write this simple mode that supports IDL code and IDLdoc comments. Ethan Gutmann had a TextMate bundle which I liked, but hadn’t been updated for TextMate 2 or newer versions of IDL.

Please let me know if you know of IDL modes for other text editors.

UPDATE 3/5: Some I missed that Jeff N. pointed out:

Notepad++ — David Higgins maintains the IDL support for Notepad++.

Cream — Support for IDL is builtin!

TaskDL is a task farming library for IDL that allows you to farm out tasks using multiple cores of a single computer or even multiple computers. It is available on Linux, OS X, and Windows. Task farming is suitable for tasks which do not need to communicate with each other, i.e. “naturally” or “embarrassingly” parallel tasks, such as processing many files independently. For more complicated programs which required interprocess communication, mpiDL provides an interface to MPI (Message Passing Interface).

As an example of using TaskDL, I will present a program to compute some areas of the Mandelbrot set and create output files representing them. If you are interesting in evaluating TaskDL or require more information about, please contact me.

The first program needed when using TaskDL is the compute task program that will be called, in our case this will be called mandelbrot_compute.pro. It is a normal IDL program that typically does not need to know about TaskDL and normally just places output in files. Our mandelbrot_compute.pro example has the following interface:

pro mandelbrot_compute, x_range, y_range, nx, ny, $
                        max_iterations=max_iterations, $
                        bound=bound, $
                        color_table=color_table, $
                        image_file=image_file, $
                        data_file=data_file, $
                        uniform_color=uniform_color

The driver of this program, mandelbrot.pro, is in charge of setting up the task farm, creating the tasks, and sending them off to the workers. To begin, creating a TaskDL object and open a new session on a particular host and port:

oFarm = obj_new('TaskDL', _extra=e)
oFarm->open_session, host=server_host, port=server_port

To run locally, server_host would simply be localhost. Port is typically just any unused port.

TaskDL has optimizations for running locally, use the ::spawn_local_worker and ::spawn_worker methods as needed to create as many workers as required, typically matching the number of processing units (cores or nodes) as available:

for w = 0L, n_workers - 1L do begin
  if (keyword_set(local)) then begin
    ofarm->spawn_local_worker
  endif else begin
    ofarm->spawn_worker, host=_worker_host[w mod n_hosts]
  endelse
endfor

The commands, as strings, to be sent to the workers much be constructed. This construction and the ::add_task call would typically be done in a loop, in our example, over the number of zoom levels desired:

cmd_format = '(%"mandelbrot_compute, %s, %s, %s, %s, ' $
               + 'max_iterations=%s, ' $
               + 'image_file=''%s'', data_file=''%s''")'
cmd = string(x_range_str, $
             y_range_str, $
             nx_str, $
             ny_str, $
             max_iterations_str, $
             image_file, $
             data_file, $
             format=cmd_format)
ofarm->add_task, cmd, queueid=0, stage=1

Multiple queues can be created associated with specific workers, but in our simple example we use the default queue. Stages provide the ability to require work to progress in stages, i.e., all stage 1 tasks must complete before stage 2 tasks start, etc. Again, that is not needed for our example.

When done, close the TaskDL session:

ofarm->close_session

Output is placed in mandelbrot-[zoom_level].png and mandelbrot-[zoom_level].nc files.

Full disclosure: I work for Tech-X and I am the product manager for the FastDL suite which includes TaskDL.

UPDATE 3/24/2014: Here is the TaskDL Users Guide.

I’ve been annoyed for a long time that GitHub marks my IDL code as Prolog in the language statistics. After a quick edit to their Linguist package that determines the language, a long wait for the pull request to get accepted, and more waiting for GitHub to rerun the language statistics on my repos, I now see the true statistics for my projects!

For example, if you go to mglib and click on the bright orange bar, you will see that mglib is 86.7% IDL, 12.1% C, and 1.2% other languages (mostly CMake).

I was told that GitHub reruns language statistics on a push, so if you have a repo that contains IDL code in .pro files that is not getting recognized correctly, try that. I was still not seeing that, though, so I had to contact GitHub support.

I’ve updated Modern IDL with some of the new features^[1] of IDL 8.3.

If you purchased a PDF of Modern IDL in the past, you should have received an email (at the address for your PayPal account) from me already with a link to download your new version. If you haven’t received it yet, please contact me.