Update: Progress on OpenMP and MEX on OS X (EDITED)

I've made some progress on getting things working.  By making the following changes to the mexpots.sh in my home home directory, I got things to compile and link directory.  First, change

CC='xcrun  -sdk macosx10.7  clang'

to

CC='llvm-gcc'

and make a likewise changes for CXX, using llvm-g++.  The other neccessary change is to ensure the lines

CFLAGS="$CFLAGS  -fopenmp"
CXXFLAGS="$CXXFLAGS -fopenmp"
LDFLAGS="$LDFLAGS -fopenmp"

are present.

Once that's taken care of, you should be able to compile a mex file that include OpenMP directives without any additional work; for straightforward C code that doesn't use any external libraries

mex your_omp_code.c -I${MATLABPATH}/extern/include

should be sufficient.  As an example, consider the following code, that also makes use of GSL Library.

This code can then be compiled as:

mex omp_reduction_mex.c -I${GSLPATH}/include -I${MATLABPATH}/extern/include -L${GSLPATH}/lib -lgsl

The only subtlety I've found in this is that to set the number of threads, you need to open a shell, run

export OMP_NUM_THREADS=2 (or 4 or 8, etc.)

and then from within the same terminal, launch MATLAB. If you just launch MATLAB, it seems to grab the maximum number of possible threads. On my machine, which is dual core with hyperthreading, it gets 4.

NOTE: I was having some weird lock ups with my code. I found that switching from the llvm compliers with gcc frontend to just the MacPorts installed gcc (I have 4.7.3 on my machine) seemed to clear things up.

Current Gripes on MATLAB, MEX, OpenMP and OS X

I write a lot of C code around which I then bind to MEX so that I can use MATLAB to manage the pre/post processing phases of my simulations. Some of my recent work has a lot of loops for computing averages that are ripe for OpenMP, so it'd be great if I could get OpenMP and MEX to play together nicely. Currently, I'm facing a couple of obstacles, many of which seem to revolve around OS X 10.8. First, the only automated MEX configuration available for OS X 10.8 uses clang, which does not support OpenMP (yet). My next thought was to write pure C, but use MATLAB's .mat interface to handle the pre/post processing script. Unfortunately, to use the libraries needed for that, they need to be set in the DLYD_LIBRARY_PATH, and OS X 10.8 seems to have dyld issues. I know I could always resort to just reading/writing ASCII files, but this all makes me grumpy about OS X as a development environment for scientific computing. I moved to OS X about 10 years ago specifically thinking that it'd give me well maintained OS with the ability to leverage all sorts of *nix software.

Getting Good Figures out of MATLAB

If you're working with MATLAB and need to produce a figure that can be included in a publication, getting the font sizes and line weights right is essential. Moreover, you should first output the figure in a vector format, either eps or pdf, and then convert to something else, such as png, if needed.

Here is an example script that produces the subsequent figure. Note the use of latex in both the figure axes and in the legend and title.

% Construct two curves
x = linspace(0,4);
u = exp(-x.^2);
v = sin(x);

% Plot the curves
plot(x, u, x, v,'linewidth',2);
set(gca, 'fontsize',16);

% Annotate the figure
xlabel('$x$','Interpreter', 'Latex');
lgnd = legend('$u$', '$v$');
set(lgnd,'interpreter','latex');
title('Example Figure','interpreter','latex');

% Save figure, note the argument 'epsc2', which saves color information
saveas(gcf,'fig1.eps','epsc2');

And here's an PNG conversion of the EPS output:

Cleaning Up Mac Ports

This is a good tip, as stuff tends to get installed for one reason or another, and then never removed.  I also don't entirely understand why macports doesn't remove the inactive ports when upgrades are installed.

Checking Progress on Asynchronous Jobs

I recently discovered the wc command in *nix environments, which is a handy solution for when I run asynchronously processed jobs.  I often need to run ensembles of simulations.  If I process them serially, and put an output command somewhere in the script, I can immediately see how far along I am by checking a log file.  But if the data is processed asynchronously, it can be a bit tougher to see what's going on.

Running the following on my log file tells me the number of lines so far:

wc -l log.out

Hence, if the log file has n lines, and I have a total of N cases, then, assuming each case produces a single output statement, n/N % of my job is done.

Graphics Grievances

When refereeing papers with plots of numerical simulations, I very often have to write back to the authors telling them that figures, while readable on the screen, where they can be blown up, are completely illegible when printed out in black and white.

When preparing images (MATLAB and matplotlib images anyways), unless you intend the figure to take up all of a letter sized paper, fonts and line widths must be checked for readability.  For the typical size figure that you'll see in a publication, a font size of at least 16, and often 18, is needed for the axes labels, tick marks, and other elements.  Line widths should generally be at least 2.  These things can easily be set in MATLAB using the exportfig script.  For matplotlib, consider editing your ~/.matplotlib/matplotlibrc file. You can set the default values for all of these graphical items as desired.

Pseudo-Spectral Methods

I was previously affiliated with the University of Toronto, where I wrote some pages for the math department's wiki on how to get pseudo-spectral algorithms up and running for nonlinear wave equations such as Nonlinear Schrödinger (NLS) and Korteweg - de Vries (KdV).  Much of this was based on Trefethen's Spectral Methods in MATLAB.

While that's still an excellent book that I recommend to people who are starting out, it has a very serious flaw in its discussion of nonlinear problems.  Trefethen's p27.m example, which simulates KdV solitons does not include a de-aliasing step.  Early in my career, this led me to think that de-aliasing was unnecessary.  That was, until, I started looking at problems with cubic nonlinearities; the p27.m example is only quadratic.

Consider focusing, cubic NLS on the torus, \(\mathbb{T}\),
\[
i u_t + u_{xx} + |u|^2 u = 0, \quad u: \mathbb{T} \times \mathbb{R} \to \mathbb{C}.
\]
When posed on the whole real line, instead of the torus, it admits soliton solutions,
\[
\Phi(x,t) = e^{i \lambda^2 t} \sqrt{2} \lambda {\rm sech}(\lambda x), \quad \lambda \neq 0.
\]
Since they are exponentially localized, if we take \( \lambda \) sufficiently large, we can approximately simulate it on the finite domain size.

Modifying p27.m for NLS, we have:

% Set discretization parameters
N = 512;
dt = .1/N;
x = (2*pi/N)*(-N/2:N/2-1)';
% Set initial condition
lam = 10;
u = sqrt(2) * lam * sech(lam * x);
v = fft(u);

% Construct integration operators
k = [0:N/2-1 0 -N/2+1:-1]';
E = exp(-.5i * dt * k.^2);
E2 = E.^2;
g = 1i * dt;

% Set simulation parameters
tmax = 5; nplt = floor((tmax/50)/dt); nmax = round(tmax/dt);
udata = u; tdata = 0;

% Integrate using the Integrating Factor method
for n = 1:nmax
  t = n*dt;
  a = g.*fft(abs( ifft(     v   ) ).^2 .* ifft(v));
  b = g.*fft(abs( ifft(E.*(v+a/2)) ).^2 .* ifft(E.*(v+a/2)) );
  c = g.*fft(abs( ifft(E.*v + b/2) ).^2 .* ifft(E.*v + b/2) );
  d = g.*fft(abs( ifft(E2.*v+E.*c) ).^2 .* ifft(E2.*v+E.*c));
  v = E2.*v + (E2.*a + 2*E.*(b+c) + d)/6;
  if mod(n,nplt) == 0
    u = ifft(v);
    udata = [udata u]; tdata = [tdata t];
  end
end

% Plot and save
figure;
for n = 1:length(tdata)
  plot(x, abs(udata(:,n)),...
    x, real(udata(:,n)), x, imag(udata(:,n)),...
    x,abs(udata(:,1)), 'linewidth',2);
  xlim([-pi, pi]);
  ylim([-16,16]);

  set(gca,'fontsize',16);
  xlabel('x');
  legend('|u|', 'Re. u', 'Im. u', '|u_0|');
  title(sprintf(' t = %g', tdata(n)));
  drawnow;
  saveas(gcf,sprintf('frame%i.png', n));
end

This generates the following output:

Notice the appearance of very high frequency oscillations.  These are aliasing errors. De-aliasing them, by zeroing out the upper half of the frequencies in Fourier, produces the following, desired, output:

This is the result of the following code.  Notice that I doubled the number of grid points to compensate for the de-aliasing.

% Set discretization parameters
N = 1024;
dt = .1/N;
x = (2*pi/N)*(-N/2:N/2-1)';
% Set initial condition
lam = 10;
u = sqrt(2) * lam * sech(lam * x);
v = fft(u);
v(N/4+2:end - N/4+1) = 0;

% Construct integration operators
k = [0:N/2-1 0 -N/2+1:-1]';
E = exp(-.5i * dt * k.^2);
E2 = E.^2;
g = 1i * dt;

% Set simulation parameters
tmax = 5; nplt = floor((tmax/50)/dt); nmax = round(tmax/dt);
udata = u; tdata = 0;

% Integrate using the Integrating Factor method
for n = 1:nmax
  t = n*dt;
  a = g.*fft(abs( ifft(     v   ) ).^2 .* ifft(v));
  a(N/4+2:end - N/4+1) = 0;
  b = g.*fft(abs( ifft(E.*(v+a/2)) ).^2 .* ifft(E.*(v+a/2)) );
  b(N/4+2:end - N/4+1) = 0;
  c = g.*fft(abs( ifft(E.*v + b/2) ).^2 .* ifft(E.*v + b/2) );
  c(N/4+2:end - N/4+1) = 0;
  d = g.*fft(abs( ifft(E2.*v+E.*c) ).^2 .* ifft(E2.*v+E.*c));
  d(N/4+2:end - N/4+1) = 0;
  v = E2.*v + (E2.*a + 2*E.*(b+c) + d)/6;
  if mod(n,nplt) == 0
    u = ifft(v);
    udata = [udata u]; tdata = [tdata t];
  end
end

% Plot and save
figure;
for n = 1:length(tdata)
  plot(x, abs(udata(:,n)),...
    x, real(udata(:,n)), x, imag(udata(:,n))...
    , 'linewidth',2);
  xlim([-pi, pi]);
  ylim([-16,16]);

  set(gca,'fontsize',16);
  xlabel('x');
  legend('|u|', 'Re. u', 'Im. u');
  title(sprintf(' t = %g', tdata(n)));
  drawnow;
  saveas(gcf,sprintf('frame%i.png', n));
end

While I still recommend Trefethen's book to people, I would supplement it with another text that addresses the aliasing problem. Two books I like are Boyd's, Chebyshev and Fourier Spectral Methods, and Canuto et al's, Spectral Methods: Fundamentals in Single Domains.