frame1.png
frame2.png
frame3.png
frame4.png
frame5.png
frame6.png
frame7.png
frame8.png
frame9.png
frame10.png

and so on, right up to frame750.png. The plan was to turn these .png files into an uncompressed movie using mencoder via the following command (original source)

mencoder mf://*.png -mf w=720:h=720:fps=25:type=png -ovc raw -oac copy -o output.avi

but I ended up with a movie that jumped all over the place since the frames were in an odd order. In the following order in fact

frame0.csv
frame100.csv
frame101.csv
frame102.csv
frame103.csv
frame104.csv
frame105.csv
frame106.csv
frame107.csv
frame108.csv
frame109.csv
frame10.csv
frame110.csv

This is because globbing expansion (the *.png bit) is alphabetical in bash rather than numerical.

One way to get the frames in the order that I want is to zero-pad them. In other words I replace file1.png with file001.png and file20.png with file020.png and so on. Here’s how to do that in bash

#!/bin/bash
num=`expr match "$1" '[^0-9]*\([0-9]\+\).*'`
paddednum=`printf "%03d" $num`
echo ${1/$num/$paddednum}

Save the above to a file called zeropad.sh and then do the following command to make it executable

chmod +x ./zeropad.sh

You can then use the zeropad.sh script as follows

./zeropad.sh frame1.png

which will return the result

frame001.png

All that remains is to use this script to rename all of the .png files in the current directory such that they are zeropadded.

for i in *.png;do mv $i `./zeropad.sh $i`; done

You may want to change the number of digits used in each filename from 3 to 5 (say). To do this just change %03d in zeropad.sh to %05d

Let me know if you find this useful or have an alternative solution you’d like to share (in another language maybe?)

These days I am not a full-time researcher (although I do assist many of them on a regular basis) but I still like to read journal articles.  The type of journal I read, however, has changed rather a lot since I’m doing it for fun and personal interest rather than for my salary.  I sometimes take what I read in these articles and turn them into blog posts and/or Wolfram Demonstrations.

Here’s a list of some of my favourites

• European Journal of Physics: According to their website “The primary mission of European Journal of Physics is to assist in maintaining and improving the standard of taught physics in universities and other institutes of higher education.”  There’s lots of cool stuff to be found with past articles including The Kaye effect, Mean Free Path in Soccer and Gases and Cooling and warming laws: an exact analytical solution
• The College Mathematics Journal: According to their website “The College Mathematics Journal is designed to enhance classroom learning and stimulate thinking regarding undergraduate mathematics.”
• The Mathematical Gazette: I could spend all day browsing through the archives of this one – it’s a math nerds dream.  For example, when you read an article written in 1894 with the title ‘Some old Text-Books‘ you know that they are going to be really old! It’s also fun to compare the ‘problems and solutions’ from 1900 to those from 2000.  From the website “The Mathematical Gazette is the original journal of the Mathematical Association and it is now over a century old”
• Mathematics Magazine: I’ve had one or two ideas for Wolfram Demonstrations from this magazine.  From the website “Mathematics Magazine presents articles and notes on undergraduate mathematical topics”

What journals do you recommend for fun and/or teaching purposes in areas such as physics, mathematics and statistics?

Update The following are recommendations from readers in the comments section. Thanks to everyone who responded.  Feel free to let me know if your favourite isn’t on this list.

• The Mathematical Intelligencer
• Notices of the AMS
• Journal of Recreational Mathematics
• American Journal of Physics

??? License checkout failed.
License Manager Error -4
Maximum number of users for Statistics_Toolbox reached.
Try again later.
To see a list of current users use the lmstat utility or contact your License Administrator.

An alternative to the Statistics Toolbox is the NAG Toolbox for MATLAB for which we have an unlimited number of licenses. Here’s how to replace ranksum with the NAG routine g08ah.

Original MATLAB / Statistics Toolbox code

x = [0.8147;0.9058;0.1270;0.9134;0.6324;0.0975;0.2785;0.5469;0.9575;0.9649];
y=  [0.4076;1.220;1.207;0.735;1.0502;0.3918;0.671;1.165;1.0422;1.2094;0.9057;0.285;1.099;1.18;0.928];
p = ranksum(x,y)

The result is p = 0.0375

Code using the NAG Toolbox for MATLAB

x = [0.8147;0.9058;0.1270;0.9134;0.6324;0.0975;0.2785;0.5469;0.9575;0.9649];
y =  [0.4076;1.220;1.207;0.735;1.0502;0.3918;0.671;1.165;1.0422;1.2094;0.9057;0.285;1.099;1.18;0.928];
tail = 'T';
[u, unor, p, ties, ranks, ifail] = g08ah(x, y, tail);

The value for p is the same as that calculated by ranksum: p = 0.0375

NAG’s g08ah routine returns a lot more than just the value p but, for this particular example, we can just ignore it all. In fact, if you have MATLAB 2009b or above then you could call g08ah like this

tail = 'T';
[~, ~, p, ~, ~, ~] = g08ah(x, y, tail);

Which explicitly indicates that you are not going to use any of the outputs other than p.

People at Manchester are using the NAG toolbox for MATLAB more and more; not only because we have a full site license for it but because it can sometimes be very fast.  Here’s some more articles on the NAG toolbox you may find useful.

• A faster version of MATLAB’s interp1
• An alternative to binopdf using the NAG toolbox for MATLAB
• Mark 22 of the NAG toolbox for MATLAB released

I work in the IT department of the University of Manchester and we are currently developing a Condor Pool which is basically a method of linking together hundreds of desktop machines to produce a high-throughput computing resource.  A MATLAB user recently submitted some jobs to our pool and complained that all of them gave identical results which is stupid because his code used MATLAB’s rand command to mix things up a bit.

I was asked if I knew why this should happen to which I replied ‘yes.’  I was then asked to advise the user how to fix the problem and I did so.  The next request was for me to write some recommendations and tutorials on how users should use random numbers in MATLAB (and Mathematica and possibly Python while I was at it) along with our Condor Pool and I went uncharacteristically quiet for a while.

It turned out that I had a lot to learn about random numbers.  This is the first of a series of (probably 2) posts that will start off by telling you what I knew and move on to what I have learned.  It’s as much a vehicle for getting the concepts straight in my head as it is a tutorial.

Ask MATLAB for 10 Random Numbers

Before we go on, I’d like you to try something for me. You have to start on a system that doesn’t have MATLAB running at all so if MATLAB is running then close it before proceeding. Now, open up MATLAB and before you do anything else issue the following command

rand(10,1)

As many of you will know, the rand command produces random numbers from the uniform distribution between 0 and 1 and the command above is asking for 10 such numbers. You may reasonably expect that the 10 random numbers that you get will be different from the 10 random numbers that I get; after all, they are random right? Well, I got the following numbers when running the above command on MATLAB 2009b running on Linux.

ans =
0.8147
0.9058
0.1270
0.9134
0.6324
0.0975
0.2785
0.5469
0.9575
0.9649

Look familiar?

Now I’ve done this experiment with a number of people over the last few weeks and the responses can be roughly split into two different camps as follows:

1. Oh yeah, I know all about that – nothing to worry about. It’s pretty obvious why it happens isn’t it?
2. It’s a bug. How can the numbers be random if MATLAB always returns the same set?

What does random mean anyway?

If you are new to the computer generation of random numbers then there is something that you need to understand and that is that, strictly speaking, these numbers (like all software generated ‘random’ numbers) are not ‘truly’ random.  Instead they are pseudorandom – my personal working definition of which is “A sequence of numbers generated by some deterministic algorithm in such a way that they have the same statistical properties of ‘true’ random numbers”.  In other words, they are not random they just appear to be but the appearance is good enough most of the time.

Pseudorandom numbers are generated from deterministic algorithms with names like Mersenne Twister, L’Ecuyer’s mrg32k3a [1]  and Blum Blum Schub whereas ‘true’ random numbers come from physical processes such as radioactive decay or atmospheric noise (the website www.random.org uses atmospheric noise for example).

For many applications, the distinction between ‘truly random’ and ‘pseudorandom’ doesn’t really matter since pseudorandom numbers are ‘random enough’ for most purposes.  What does ‘random enough’ mean you might ask?  Well as far as I am concerned it means that the random number generator in question has passed a set of well defined tests for randomness – something like Marsaglia’s Diehard tests or, better still, L’Ecuyer and Simard’s TestU01 suite will do nicely for example.

The generation of random numbers is a complicated topic and I don’t know enough about it to do it real justice but I know a man who does.  So, if you want to know more about the theory behind random numbers then I suggest that you read Pierre L’Ecuyer’s paper simply called ‘Random Numbers’ (pdf file).

Back to MATLAB…

Always the same seed

So, which of my two groups are correct?  Is there a bug in MATLAB’s random number generator or not?

There is nothing wrong with MATLAB’s random number generator at all. The reason why the command rand(10,1) will always return the same 10 numbers if executed on startup is because MATLAB always uses the same seed for its pseudorandom number generator (which at the time of writing is a Mersenne Twister) unless you tell it to do otherwise.

Without going into details, a seed is (usually) an integer that determines the internal state of a random number generator.  So, if you initialize a random number generator with the same seed then you’ll always get the same sequence of numbers and that’s what we are seeing in the example above.

Sometimes, this behaviour isn’t what we want.  For example, say I am doing a Monte Carlo simulation and I want to run it several times to verify my results.  I’m going to want a different sequence of random numbers each time or the whole exercise is going to be pointless.

One way to do this is to initialize the random number generator with a different seed at startup and a common way of achieving this is via the system clock.  The following comes straight out of the current MATLAB documentation for example

RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(100*clock)));

Do this once per MATLAB session and you should be good to go (there is usually no point in doing it more than once per session by the way….your numbers won’t be any ‘more random’ if you so.  In fact, there is a chance that they will become less so!).

Condor and ‘random’ random seeds

Sometimes the system clock approach isn’t good enough either.  For example, at my workplace, Manchester University, we have a Condor Pool of hundreds of desktop machines which is perfect for people doing Monte Carlo simulations.  Say a single simulation takes 5 hours and it needs to be run 100 times in order to get good results.  On one machine that’s going to take about 3 weeks but on our Condor Pool it can take just 5 hours since all 100 simulations run at the same time but on different machines.

If you don’t think about random seeding at all then you end up with 100 identical sets of results using MATLAB on Condor for the reasons I’ve explained above.  Of course you know all about this so you switch to using the clock seeding method, try again and….get 100 identical sets of results[2].

The reason for this is that the time on all 100 machines is synchronized using internet time servers.  So, when you start up 100 simultaneous jobs they’ll all have the same timestamp and, therefore, have the same random seed.

It seems that what we need to do is to guarantee (as far as possible) that every single one of our condor jobs gets a unique seed in order to provide a unique random number stream and one way to do this would be to incorporate the condor process ID into the seed generation in some way and there are many ways one could do this.  Here, however, I’m going to take a different route.

On Linux machines it is possible to obtain small numbers of random numbers using the special files /dev/random and /dev/urandom which are interfaces to the kernel’s random number generator.  According to the documentation ‘The random number generator gathers environmental noise from device drivers and other sources into an entropy pool. The generator also keeps an estimate of the number of bit of the noise in the entropy pool.  From this entropy pool random numbers are created.’

This kernel generator isn’t suitable for simulation purposes but it will do just fine for generating an initial seed for MATLAB’s pseudorandom number generator.  Here’s the MATLAB commands

[status seed] = system('od /dev/urandom --read-bytes=4 -tu | awk ''{print $2}''')
seed=str2double(seed)
RandStream.setDefaultStream(RandStream('mt19937ar','seed',seed));

Put this at the beginning of the MATLAB script that defines your condor job and you should be good to go.  Don’t do it more than once per MATLAB session though – you won’t gain anything!

The end or the end of the beginning?

If you asked me the question ‘How do I generate a random seed for a pseudorandom number generator?’ then I think that the above solution answers it quite well.  If, however, you asked me ‘What is the best way to generate multiple independent random number streams that would be good for thousands of monte-carlo simulations?‘ then we need to rethink somewhat for the following reasons.

Seed collisions: The Mersenne twister algorithm currently used as the default random number generator for MATLAB uses a 32bit integer seed which means that it can take on 2^32 or 4,294,967,296 different values – which seems a lot!  However, by considering a generalisation of the birthday problem it can be seen that if you select such a seed at random then you have a 50% chance choosing two identical seeds after only 65,536 runs.  In other words, if you perform 65,536 simulations then there is a 50% chance that two such simulations will produce identical results.

Bad seeds: I have read about (but never experienced) the possibility of ‘bad seeds’; that is some seeds that produce some very strange, non-random results – for the first few thousand iterates at least.  This has led to some people advising that you should ‘warm-up’ your random number generator by asking for, and throwing away, a few thousand random numbers before you start using them. Does anyone know of any such bad seeds?

Overlapping or correlated sequences: All pseudorandom number generators are periodic (at least, all the ones that I know about are) – which means that after N iterations the sequence repeats itself.  If your generator is good then N is usually large enough that you don’t need to worry about this.  The Mersenne Twister used in MATLAB, for instance, has a huge period of (2^19937 – 1)/2 (half of the standard 32bit implementation).

The point I want to make is that you don’t get several different streams of random numbers, you get just one, albeit a very big one.  Now, when you choose a seed you are essentially choosing a random point in this stream and there is no guarantee how far apart these two points are.  They could be separated by a distance of trillions of points or they could be right next to each other – we simply do not know – and this leads to the possibility of overlapping sequences.

Now, one could argue that the possibility of overlap is very small in a generator such as the Mersenne Twister and I do not know of any situation where it has occurred in practice but that doesn’t mean that we shouldn’t worry about it.  If your work is based on the assumption that all of your simulations have used independent, uncorrelated random number streams then there is a possibility that your assumptions could be wrong which means that your conclusions could be wrong.  Unlikely maybe, but still no way to do science.

Next Time

Next time I’ll be looking at methods for generating guaranteed independent random number streams using MATLAB’s in-built functions as well as methods taken from the NAG Toolbox for MATLAB.  I’ll also be including explicit examples that use this stuff in Condor.

Ask the audience

I assume that some of you will be in the business of performing Monte-Carlo simulations and so you’ll probably know much more about all of this than me.  I have some questions

• Has anyone come across any ‘bad seeds’ when dealing with MATLAB’s Mersenne Twister implementation?
• Has anyone come across overlapping sequences when using MATLAB’s Mersenne Twister implementation?
• How do YOU set up your random number generator(s).

I’m going to change my comment policy for this particular post in that I am not going to allow (most) anonymous comments.  This means that you will have to give me your email address (which, of course, I will not publish) which I will use once to verify that it really was you that sent the comment.

Notes and References

[1] L’Ecuyer P (1999) Good parameter sets for combined multiple recursive random number generators Operations Research 47:1 159–164

[2] More usually you’ll get several different groups of results.  For example you might get 3 sets of results, A B C, and get 30 sets of A, 50 sets of B and 20 sets of C.  This is due to the fact that all 100 jobs won’t hit the pool at precisely the same instant.

[3] Much of this stuff has already been discussed by The Mathworks and there is an excellent set of articles over at Loren Shure’s blog – Loren onThe Art of MATLAB.

• New ways with random numbers, Part 1
• New ways with random numbers, Part 2

Erroneous write during file extend. write -1 instead of 4096
Probably out of disk space.
Erroneous write during file extend. write -1 instead of 4096
Probably out of disk space.
Write error in NtrExt1
Write error in NtrExt1: Bad address
Segmentation fault

A bit of googling suggested that the following might work

sudo echo 0 > /proc/sys/kernel/randomize_va_space

but this will result in permission denied (explanation here). The command you really want to use is

sudo bash -c "echo 0 > /proc/sys/kernel/randomize_va_space"

Once this was done, Gaussian worked as advertised. Maybe this post will help a googler sometime in the future.

X=rand(10,2);
dists=pdist(X,'euclidean');

It’s a nice function but the problem with it is that it is part of the Statistics Toolbox and that costs extra. I was recently approached by a user who needed access to the pdist function but all of the statistics toolbox license tokens on our network were in use and this led to the error message

??? License checkout failed.
License Manager Error -4
Maximum number of users for Statistics_Toolbox reached.
Try again later.
To see a list of current users use the lmstat utility or contact your License Administrator

One option, of course, is to buy more licenses for the statistics toolbox but there is another way. You may have heard of GNU Octave, a free,open-source MATLAB-like program that has been in development for many years.  Well, Octave has a sister project called Octave-Forge which aims to provide a set of free toolboxes for Octave.  It turns out that not only does Octave-forge contain a statistics toolbox but that toolbox contains an pdist function.  I wondered how hard it would be to take Octave-forge’s pdist function and modify it so that it ran on MATLAB.

Not very!  There is a script called oct2mat that is designed to automate some of the translation but I chose not to use it – I prefer to get my hands dirty you see.  I named the resulting function octave_pdist to help clearly identify the fact that you are using an Octave function rather than a  MATLAB function.  This may matter if one or the other turns out to have bugs.  It appears to work rather nicely:

dists_oct = octave_pdist(X,'euclidean');
% let's see if it agrees with the stats toolbox version
all( abs(dists_oct-dists)<1e-10)

ans =
     1

Let’s look at timings on a slightly bigger problem.

>> X=rand(1000,2);
>> tic;matdists=pdist(X,'euclidean');toc
Elapsed time is 0.018972 seconds.
>> tic;octdists=octave_pdist(X,'euclidean');toc
Elapsed time is 6.644317 seconds.

Uh-oh! The Octave version is 350 times slower (for this problem) than the MATLAB version. Ouch! As far as I can tell, this isn’t down to my dodgy porting efforts, the original Octave pdist really does take that long on my machine (Ubuntu 9.10, Octave 3.0.5).

This was far too slow to be of practical use and we didn’t want to be modifying algorithms so we ditched this function and went with the NAG Toolbox for MATLAB instead (routine g03ea in case you are interested) since Manchester effectively has an infinite number of licenses for that product.

If,however, you’d like to play with my MATLAB port of Octave’s pdist then download it below.

• octave_pdist.m makes use of some functions in the excellent NaN Toolbox so you will need to download and install that package first.

Some time later I followed this up with another blog post and a Wolfram Demonstration.

Well, over at Stack Overflow, some people have been rendering this cool equation using MATLAB. Here’s the first version

x = linspace(-3,3,50);
y = linspace(-5,5,50);
[X Y]=meshgrid(x,y);
Z = exp(-X.^2-Y.^2/2).*cos(4*X) + exp(-3*((X+0.5).^2+Y.^2/2));
Z(Z>0.001)=0.001;
Z(Z<-0.001)=-0.001;
surf(X,Y,Z);
colormap(flipud(cool))
view([1 -1.5 2])

and here’s the second.

[x y] = meshgrid( linspace(-3,3,50), linspace(-5,5,50) );
z = exp(-x.^2-0.5*y.^2).*cos(4*x) + exp(-3*((x+0.5).^2+0.5*y.^2));
idx = ( abs(z)>0.001 );
z(idx) = 0.001 * sign(z(idx)); 

figure('renderer','opengl')
patch(surf2patch(surf(x,y,z)), 'FaceColor','interp');
set(gca, 'Box','on', ...
    'XColor',[.3 .3 .3], 'YColor',[.3 .3 .3], 'ZColor',[.3 .3 .3], 'FontSize',8)
title('$e^{-x^2 - \frac{y^2}{2}}\cos(4x) + e^{-3((x+0.5)^2+\frac{y^2}{2})}$', ...
    'Interpreter','latex', 'FontSize',12) 

view(35,65)
colormap( [flipud(cool);cool] )
camlight headlight, lighting phong

Do you have any cool graphs to share?

SpaceTime Mathematics

The Rolls Royce of mobile mathematical applications and one that I have been using since my days as a Windows Mobile user.  The iPad version was one of the first apps I bought when I received my device and it is just beautiful!  Symbolic algebra and calculus, 2 and 3D interactive plotting, scripting, fractals linear algebra…the list of functions just goes on and on.  I would have loved to have access to this app when I was in high school or early university.

If you want to get an idea of the quality of SpaceTime’s graphical capabilities then check out the free demo, Graphbook, but be aware that there is a lot more to SpaceTime than just graphics.

Regular readers of Walking Randomly will know that I am a big fan of Wolfram’s Demonstration project which is made possible by Mathematica’s Manipulate function.  Well, SpaceTime has a similar, albeit simplified, version of Manipulate – a function called Scroll.  Interactive Fourier Series on the iPad anyone?

Something else that I like about SpaceTime is the fact that it is cross-platform with versions for Linux and Windows available in addition to iPhone, iPad and Windows Mobile.  So, students could use it in a classroom setting on PCs and use what they have learned on their own iPad/iPhone version.

If you only buy one mathematical application for iPad then this should be it.  It’s relatively expensive for an iPad app at £11.99 (at the time of writing) but is worth every penny and I bought it without hesitation – so should you!

• SpaceTime iTunes Link
• Main website for SpaceTime

PocketCAS Pro

PocketCAS Pro is a computer algebra system that started out life as a Windows Mobile app and is now available for iPhone and iPad.  I haven’t had chance to try it out yet so I can’t comment on its quality but it has a lot of features including symbolic algebra and calculus, 2D plotting, numerical solution of equations and more.

At the time of writing, it is the same price as SpaceTime mathematics – £11.99 – and yet my first impression is that it has less functionality.   No 3d plotting for example.  I’ll know more when I buy a copy next month.

There is a free lite version available which includes some of the functionality of the main product to allow you to try it out.

• PocketCAS Pro on iTunes.
• PocketCAS lite – free, cut down version of PocketCAS Pro
• PocketCAS main website

fxIntegrator

My favourite operating system is Linux where there is a philosophy of “Write programs that do one thing and do it well”.  fxIntegrator does one thing -the numerical solution of 1d integrals – but does it do it well?

Well, it’s not bad.  You enter the function you want to integrate using the nice, specially designed keyboard, then you enter the limits and press the = button to get the result.  Couldn’t get any easier and I like it.  The equation editor is very nice resulting in well formatted integrands but I did manage to confuse it once or twice.  FxIntegrator is also very cheap at only 59p – a real bargain!

I tried a few straightforward integrals on it and it gave the correct answer in all cases.  Then I got nasty and tried the following which has an algebraic-logarithmic singularity at the origin (original source for this integral).

I wasn’t expecting fxIntegrator to cope and it didn’t. Rather than giving the answer I just got an unhappy face indicating that it couldn’t compute the solution.  This isn’t a criticism though! I like the fact that rather than giving numerical garbage, fxIntegator simply said ‘I can’t do that’.

There are some niggles, however.  First of all, the list of elementary functions available is rather limited as it only includes square roots, powers,the trigonometric functions sin,cos and tan, the natural logarithm function ln and basic arithmetic.  Even when I was in high school I would have wanted more such as inverse trig functions.

Another problem with it is that although you can use the customised keyboard to enter the integrand, when you try to enter limits the standard iPad keyboard pops up.

These niggles aside, however, this is a nice little app for 59p and I hope the author continues to develop it.  If he does then here are some suggestions for functionality I’d like to see.

• Add a few more functions.  Inverse trig for a start.  If possible then maybe things such as Bessel functions.
• Help turn this into a better teaching and learning tool by implementing a range of numerical methods for computing the integrand and allow the user to choose between them.  Methods such as the rectangle rule, trapezoidal rule and simpson’s rule along with the ability to change the sub-divison size.  The more methods the better
• Perhaps add some tutorial notes on each numerical method.
• Give the calculation time for the result in seconds along with the number of evaluations of the integrand.  This will help students compare the trade off in speed/accuracy of each method.
• Add the ability to plot the integrand along with the limits.  Allow the user to change limits by moving them on the graph as well as by direct input.  Once the calculation is performed, show the points on the curve where the algorithm sampled the function.

This good little app could be turned into a great little app.

• iTunes page for fxIntegrator
• Main wesbite for fxIntegrator

More articles from Walking Randomly on mobile Mathematics

• PocketCAS – a new computer algebra system for iPhone
• Detexify answers ‘What’s the LaTeX code for this symbol?’
• Review of Gnuplot 4.2 for Windows Mobile

All this got me thinking….how fast is my mobile phone compared to these computational behemoths?

The key to answering this question lies with the Linpack benchmarks developed by Jack Dongarra back in 1979.  Wikipedia explains:

‘they [The Linpack benchmarks] measure how fast a computer solves a dense N by N system of linear equations Ax = b, which is a common task in engineering. The solution is obtained by Gaussian elimination with partial pivoting, with 2/3·N3 + 2·N2 floating point operations. The result is reported in millions of floating point operations per second (MFLOP/s).’

People have been using the Linpack benchmarks to measure the speed of computers for decades and so we can use the historical results to see just how far computers have come over the last thirty years or so.  Back in 1979, for example, the fastest computer on the block according to the N=100 Linpack benchmark was the Cray 1 supercomputer which had a measured speed of 3.4 Mflop/s per processor.

More recently, a Java version of the Linpack benchmark was developed and this was used by GreeneComputing to produce an Android version of the benchmark.

I installed the benchmark onto my trusty T-Mobile G2 (a rebadged HTC Hero, currently running Android 1.5) and on firing it up discovered that it tops out at around 2.3 Mflop/s which makes it around 66% as fast as a single processor on a 1979 Cray 1 supercomputer.  OK, so maybe that’s not particularly impressive but the very latest crop of Android phones are a different matter entirely.

According to the current Top 10 Android Linpack results, a tweaked Motorola Droid is capable of scoring 52 Mflop/s which is over 15 times faster than the 1979 Cray 1 CPU.  Put another way, if you transported that mobile phone back to 1987 then it would be on par with the processors in one of the fastest computers in the world of the time, the ETA 10-E, and they had to be cooled by liquid nitrogen.

Like all benchmarks, however, you need to take this one with a pinch of salt.  As explained on the Java Linpack page ‘This [the Java version of the] test is more a reflection of the state of the Java systems than of the floating point performance of the underlying processors.’ In other words, the underlying processors of our mobile phones are probably faster than these Java based tests imply.

How about you though?  Will this new application be useful for the way you work?

More articles from Walking Randomly on mobile mathematics software

• Mathematics on iPad #1
• Mathematics on the Pocket PC
• PocketCAS – a new computeger algebra system for the iPhone

1031,-948,-76
507,635,-1148
-1031,948,750
-507,-635,114

The following, very simple command, is all that you need

>> veryclean = csvread('very_clean.txt')

veryclean =

        1031        -948         -76
         507         635       -1148
       -1031         948         750
        -507        -635         114

In the real world, however, your data is rarely this nice and clean. One of the most common problems faced by MATLABing data importers is that of header lines. Take the file quite_clean.txt for instance. This is identical to the previous example apart from the fact that it contains two header lines

These are some data that I made using my hand-crafted code
Date:12th July 1996
1031,-948,-76
507,635,-1148
-1031,948,750
-507,-635,114

This is all too much for the csvread command

>> data=csvread('quite_clean.txt')
??? Error using ==> dlmread at 145
Mismatch between file and format string.
Trouble reading number from file (row 1, field 1) ==> This

Error in ==> csvread at 52
    m=dlmread(filename, ',', r, c);

Not to worry, we can just use the more capable importdata command instead

>> quiteclean = importdata('quite_clean.txt')

quiteclean = 

        data: [4x3 double]
    textdata: {2x1 cell}

The result above is a two element structure array and our numerical values are contained in a field called data. Here’s how you get at it.

>> quiteclean.data

ans =

        1031        -948         -76
         507         635       -1148
       -1031         948         750
        -507        -635         114

So far so good. How do we handle a file like messy_data.txt though?

header 1;
header 2;
1031,-948,-76, ,"12"
507,635,-1148, ,"34"
-1031,948,750, ,"45"
-507,-635,114, ,"67"

This is the kind of file encountered by Walking Randomly reader ‘reen’ and it contains exactly the same numerical values as the previous two examples. Unfortunately, it also contains some cruft that makes life more difficult for us. Let’s bring out the big-guns!

Using textscan to import csv files in MATLAB

When the going gets tough, the tough use textscan.  Here’s the incantation for importing messy_data.txt

fid=fopen('messy_data.txt');
data = textscan(fid,'%f %f %f %*s %*s','HeaderLines',2,'Delimiter',',','CollectOutput',1);
fclose(fid)

The result is a one-element cell array that contains an array of doubles.  Let’s get the array of doubles out of the cell

>> data=data{1}
data =
        1031        -948         -76
         507         635       -1148
       -1031         948         750
        -507        -635         114

If the importdata command is a chauffeur then textscan is a Ferrari and I don’t know about you but I’d much rather be driving my own Ferrari than being chauffeured around (John Cook over at The Endeavour has more to say on Ferraris and Chauffeurs).

Let’s de-construct the above set of commands.  The first thing to notice is that, unlike csvread and importdata, you have to explicitly open and close your file when using the textscan command.  So, you open your file using fopen and give it a file ID (which in this example is fid).

fid=fopen('messy_data.txt');

The first argument to textscan is just this file ID, fid. Next you need to supply a conversion specifier which in this case is

'%f %f %f %*s %*s'

The conversion specifier tells textscan what you want each row in your csv file to be converted to. %f means “64 bit double” and %s means “string” so ‘%f %f %f %s %s’ means “3 doubles followed by 2 strings” (we’ll get onto the asterisks in the original specifier later). You can use all sorts of data types in a conversion specifier such as integers, quoted strings and pattern matched strings among others. Check out the MATLAB documentation for textscan for the full list but an abbreviated list is shown below:

%d signed 32bit integer
%u unsigned 32bit integer
%f 64bit double (you'll want this most of the time when using MATLAB)
%s string

Now, in the command I used to import messy_data.txt the conversion specifier contained some asterisks such as %*s so what do these mean?  Quite simply, the asterisk just means ‘ignore’ so %*s means ‘ignore the string in this field’.  So, the full meaning of my conversion specifier ‘%f %f %f %*s %*s’ is “read 3 doubles and ignore 2 strings” and textscan will do this for every row.

The rest of the command is pretty self explanatory but I’ll explain it anyway for the sake of completeness

'HeaderLines',2

The file has 2 headerlines which should be ignored

'Delimiter',','

The fields are delimited (a posh word for separated) by a comma

'CollectOutput',1

If you supply a 1 (which stands for True) to the CollectOutput option then textscan will join consecutive output cells with the same data type into a single array. Since I want all of my doubles to be in a single array then this is the behaviour I went for.

Finally, once you have finished textscanning, don’t forget to close your file

fclose(fid)

That’s pretty much it for this mini-tutorial – I hope you find it useful.

• Accelerated linear algebra using graphics cards via NVIDIA’s CUDA.    Maple’s advertising blurb says that they have implemented matrix multiplication and that this will help speed up many linear algebra routines since this is such a fundamental operation.  I think that Maple are the first of the big general purpose mathematical packages to offer direct CUDA integration out of the box and this development is well worth watching.  The speed-ups that are possible from CUDA technology are nothing short of astonishing – hundreds of times in some cases.  However, Maplesoft are going to need to add a lot more than matrix multiplication in order for this to be truly useful.  A set of fast random number generators would be nice for example (I’m thinking superfast Monte-carlo simulations – the finance people would love it).
• Maple uses Intel’s Math Kernel Library (MKL) for many of its low-level numerical linear algebra routines and this has been updated to version 10.0 in Maple 14.  For 32bit windows users this has sped certain operations up quite a lot but it is 64bit Windows users who will really see the benefit since 64bit Maple 13 only used a set of generic BLAS routines.  The practical upshot is that certain basic linear algebra routines, such as Matrix Multiplication, can be around 10 times faster in 64bit windows Maple 14 compared to the previous version.  I couldn’t find mention of the Linux version.
• A shed load of updates to their differential equation solvers including a new numerical routine called the Cash-Karp pair.
• The Maple toolbox for MATLAB is no longer a separate product and is now included with Maple itself.  This is great news if you, like me, tend to work with several mathematical packages simultaneously.  Of course you need to have a copy of MATLAB installed to make use of this functionality – you don’t get a copy of MATLAB for free
• You can now import MATLAB binary files (compressed and uncompressed) directly into Maple using the ImportMatrix command.
• Another product, The Maple-NAG connector, has also been integrated with Maple itself.  This allows you to easily call the NAG C library directly from Maple but, similar to the MATLAB toolbox, you’ll have to purchase the NAG C library separately to make use of this.

As you can see, I tend to favour new features that lead to improved performance or better interoperability with other software packages in the first instance.  New mathematics and usability features take a little longer to sink in (for me at least).

I’ve not got a copy of Maple 14 yet but will try to write more if I upgrade (finances permitting).

For a full list of changes check out Maple’s online help section.

More on Maple from Walking Randomly

• Japanese firm buys Maplesoft
• Simulating Santa using Maple
• Parallel Programming in Maple

a=int64(10);
b=int64(20);
a+b
??? Undefined function or method 'plus' for input arguments of type 'int64'.

It doesn’t like any of the other operators either

>> a-b
??? Undefined function or method 'minus' for input arguments of type 'int64'.

>> a*b
??? Undefined function or method 'mtimes' for input arguments of type 'int64'.

>> a/b
??? Undefined function or method 'mrdivide' for input arguments of type 'int64'.

At first I thought that there was something wrong with my MATLAB installation but it turns out that this behaviour is expected and documented. At the time of writing, the MATLAB documentation contains the lines

Note   Only the lower order integer data types support math operations.
Math operations are not supported for int64 and uint64.

So, there you have it. When can’t MATLAB add up? When you ask it to add 64 bit integers!

Update: Just had an email from someone who points out that Octave (Free MATLAB Clone) can handle 64bit integers just fine

octave:1> a=int64(10);
octave:2> b=int64(20);
octave:3> a+b
ans = 30
octave:4>

Update 2: Since I got slashdotted, people have been asking why I (or, more importantly, the user I was helping) needed 64bit integers. Well, we were using the NAG Toolbox for MATLAB which is a MATLAB interface to the NAG Fortran library. Some of its routines require integer arguments and on a 64bit machine these must be 64bit integers. We just needed to do some basic arithmetic on them before passing them to the toolbox and discovered that we couldn’t. The work-around was simple – use int32s for the arithmetic and then convert to int64 at the end so no big deal really. I was simply surprised that arithmetic wasn’t supported for int64s directly – hence this post.

Update 3: In the comments, someone pointed out that there is a package on the File Exchange that adds basic arithmetic support for 64bit integers.  I’ve not tried it myself though so can’t comment on its quality.

Update 4: Someone has made me aware of an interesting discussion concerning this topic on MATLAB Central a few years back: http://www.mathworks.com/matlabcentral/newsreader/view_thread/154955

I also spent a few minutes talking about Python including why I thought it was a good language for numerical computing and how to interface it with the NAG C library using my tutorials and PyNAG code.  I finished off with a simple graphical demonstration that minimised the Rosenbrock function using Python, Matplotlib and the NAG C-library running on Ubuntu Linux.

“So!”, I asked, “Any questions?”

The very first reply was “Does your Python-NAG interface work on Windows machines?” to which the answer at the time was “No!”  I took the opportunity to ask the audience how many of them did their numerical computing in Windows (Most of the room of around 50+ people), how many people did it using Mac OS X (A small group at the back), and how many people did it in Linux (about 3).

So, if I wanted the majority of that audience to use PyNAG then I had to get it working on Windows.  Fortunately, thanks to the portability of Python and the consistency of the NAG library across platforms, this proved to be rather easy to do and the result is now available for download on the main PyNAG webpage.

Let’s look at the main differences between the Linux and Windows versions

Loading the NAG Library

In the examples given in my Linux NAG-Python tutorials, the cllux08dgl NAG C-library was loaded using the line

libnag = cdll.LoadLibrary("/opt/NAG/cllux08dgl/lib/libnagc_nag.so.8")

In Windows the call is slightly different. Here it is for CLDLL084ZL (Mark 8 of the C library for Windows)

C_libnag = ctypes.windll.CLDLL084Z_mkl

If you are modifying any of the codes shown in my NAG-Python tutorials then you’ll need to change this line yourself. If you are using PyNAG, however, then it will already be done for you.

Callback functions

For my full introduction to NAG callback functions on Linux click here. The main difference between using callbacks on Windows compared to Linux is that on Linux we have

C_FUNC = CFUNCTYPE(c_double,c_double )

but on Windows we have

C_FUNC = WINFUNCTYPE(c_double,c_double )

There are several callback examples in the PyNAG distribution.

That’s pretty much it I think.  Working through all of the PyNAG examples and making sure that they ran on Windows uncovered one or two bugs in my codes that didn’t affect Linux for one reason or another so it was a useful exercise all in all.

So, now you head over to the main PyNAG page and download the Windows version of my Python/NAG interface which includes a set of example codes.  I also took the opportunity to throw in a couple of extra features and so upgraded PyNAG to version 0.16, check out the readme for more details.  Thanks to several employees at NAG for all of their help with this including Matt, John, David, Marcin and Sorin.

Yes, we are big users of NAG at Manchester but then that stands to reason because NAG and Manchester have a collaborative history stretching back 40 years to NAG’s very inception.  Manchester takes a lot of NAG’s products but for reasons that are lost in the mists of time, we have never (to my knowledge at least) had a site license for their SMP library (more recently called The NAG Library for SMP & multicore).  Very recently, that changed!

SMP stands for Symmetric Multi-Processor which essentially means ‘two or more CPUs sharing the same bit of memory.’  Five years ago, it would have been rare for a desktop user to have an SMP machine but these days they are everywhere.  Got a dual-core laptop or a quad-core desktop?  If the answer’s ‘yes’ then you have an SMP machine and you are probably wondering how to get the most out of it.

‘How can I use all my cores (without any effort)’

One of the most common newbie questions I get asked these days goes along the lines of ‘Program X is only using 50%/25%/12.5% of my CPU – how can I fix this?’ and, of course, the reason for this is that the program in question is only using a single core of their multicore machine.  So, the problem is easy enough to explain but not so easy to fix because it invariably involves telling the user that they are going to have to learn how to program in parallel.

Explicit parallel programming is a funny thing in that sometimes it is trivial and other  times it is pretty much impossible – it all depends on the problem you see.  Sometimes all you need to do is drop a few OpenMP pragmas here and there and you’ll get a 2-4 times speed up.  Other times you need to completely rewrite your algorithm from the ground up to get even a modest speed up.  Yet more times you are simply stuck because your problem is inherently non-parallelisable.  It is even possible to slow your code down by trying to parallelize it!

If you are lucky, however, then you can make use of all those extra cores with no extra effort at all!  Slowly but surely, mathematical software vendors have been re-writing some of their functions to ensure that they work efficiently on parallel processors in such a way that it is transparent to the user.  This is sometimes referred to as implicit parallelism.

Take MATLAB, for example, ever since release 2007a more and more built in MATLAB functions have been modified to allow them to make full use of multi-processor systems.  If your program makes extensive use of these functions then you don’t need to spend extra money or time on the parallel computing toolbox, just run your code on the latest version of MATLAB, sit back and enjoy the speed boost.  It doesn’t work for all functions of course but it does for an ever increasing subset.

The NAG SMP Library – zero effort parallel programming

For users of NAG routines, the zero-effort approach to making better use of your multicore system is to use their SMP library.  According to NAG’s advertising blurb you don’t need to rewrite your code to get a speed boost – you just need to link to the SMP library instead of the Fortran library at compile time.

Just like MATLAB, you won’t get this speed boost for every function, but you will get it for a significant subset of the library (around 300+ functions as of Mark 22 – the full list is available on NAG’s website).  Also, just like MATLAB, sometimes the speed-up will be massive and other times it will be much more modest.

I wanted to test NAG’s claims for myself on the kind of calculations that researchers at Manchester tend to perform so I asked NAG for a trial of the upcoming Mark 22 of the SMP library and, since I am lazy, I also asked them to provide me with simple demonstration programs and I’d like to share the results with you all.  These are not meant to be an exhaustive set of benchmarks (I don’t have the time to do such a thing) but should give you an idea of what you can expect from NAG’s new library.

System specs

All of these timings were made on a 3Ghz Intel Core2 Quad Q9650 CPU desktop machine with 8Gb of RAM running Ubuntu 9.10.  The serial version of the NAG library used was fll6a22dfl  and the SMP version of the library was fsl6a22dfl.  The version of gfortran used was 4.4.1 and the cpu-frequency governor was switched off as per a previous blog post of mine.

Example 1 – Nearest correlation matrix

The first routine I looked at was one that calculated the nearest correlation matrix. In other words ‘Given a symmetric matrix, what is the nearest correlation matrix—that is, the nearest symmetric positive semidefinite matrix with unit diagonal?’^[1] This is a problem that pops up a lot in Finance – option pricing, risk management – that sort of thing.

The NAG routine that calculates the nearest correlation matrix is G02AAF which is based on the algorithm developed by Qi and Sun^[2].  The example program I was sent first constructs a N x N tridiagonal matrix A of the form A(j,j)=2, A(j-1,j)=-1.1 and A(j,j-1)=1.0.  It then times how long it takes G02AAF to find the nearest correlation matrix to A.  You can download this example, along with all supporting files, from the links below.

To compile this benchmark against the serial library I did

gfortran ./bench_g02aaf.f90 shtim.c wltime.f /opt/NAG/fll6a22dfl/lib/libnag_nag.a  \
/opt/NAG/fll6a22dfl/acml/libacml_mv.a -o serial_g02aaf

To compile the parallel version I did

gfortran -fopenmp ./bench_g02aaf.f90 shtim.c wltime.f /opt/NAG/fsl6a22dfl/lib/libnagsmp.a \
  /opt/NAG/fsl6a22dfl/acml4.3.0/lib/libacml_mp.a -o smp_g02aaf

I didn’t need to modify the source code in any way when going from a serial version to a parallel version of this benchmark. The only work required was to link to the SMP library instead of the serial library – so far so good. Let’s run the two versions and see the speed differences.

First things first, let’s ensure that there is no limit to the stack size of my shell by doing the following in bash

ulimit -s unlimited

Also, for the parallel version, I need to set the number of threads to equal the number of processors I have on my machine by setting the OMP_NUM_THREADS environment variable.

export OMP_NUM_THREADS=4

This won’t affect the serial version at all. So, onto the program itself. When you run it, it will ask you for two inputs – an integer, N, and a boolean, IO. N gives the size of the starting matrix and IO determines whether or not to output the result.

Here’s an example run of the serial version with N=1000 and IO set to False. (In Fortran True is represented as .t. and False is represented as .f. )

./serial_g02aaf
 G02AAF: Enter N, IO
1000 .f.
 G02AAF: Time, IFAIL =    37.138997077941895                0
 G02AAF: ITER, FEVAL, NRMGRD =            4           5  4.31049822255544465E-012

The only output I am really interested in here is the time and 37.13 seconds doesn’t seem too bad for a 1000 x 1000 matrix at first glance. Move to the parallel version though and you get a very nice surprise

./smp_g02aaf
 G02AAF: Enter N, IO
1000 .f.
 G02AAF: Time, IFAIL =    5.1906139850616455                0
 G02AAF: ITER, FEVAL, NRMGRD =            4           5  4.30898281428799420E-012

The above times were typical and varied little from run to run (although the SMP version varied by a bigger percentage than the serial version). However I averaged over 10 runs to make sure and got 37.14 s for the serial version and 4.88 s for the SMP version which gives a speedup of 7.6 times!

Now, this is rather impressive. Usually, when one parallelises over N cores then the very best you can expect in an ideal word is a speed up of just less than N times, so called ‘linear scaling’. Here we have N=4 and a speedup of 7.6 implying that NAG have achieved ’super-linear scaling’ which is usually pretty rare.

I dropped them a line to ask what was going on. It turns out that when they looked at parallelising this routine they worked out a way of speeding up the serial version as well. This serial speedup will be included in the serial library in its next release, Mark 23 but the SMP library got it as of Mark 22.

So, some of that 7.6 times speedup is as a result of serial speedup and the rest is parallelisation. By setting OMP_NUM_THREADS to 1 we can force the SMP version of the benchmark to only run on only one core and thus see how much of the speedup we can attribute to parallelisation:

export OMP_NUM_THREADS=1
./smp_g02aaf
 G02AAF: Enter N, IO
1000 .f.
 G02AAF: Time, IFAIL =    12.714214086532593                0
 G02AAF: ITER, FEVAL, NRMGRD =            4           5  4.31152424170390294E-012

Recall that the 4 core version took an average of 4.88 seconds so the speedup from parallelisation alone is 2.6 times – much closer to what I expected to see.  Now, it is probably worth mentioning that there is an extra level of complication (with parallel programming there is always an extra level of complication) in that some of this parallelisation speedup comes from extra work that NAG has done in their algorithm and some of it comes from the fact that they are linking to parallel versions of the BLAS and LAPACK libraries.  We could go one step further and determine how much of the speed up comes from NAG’s work and how much comes from using parallel BLAS/LAPACK but, well, life is short.

The practical upshot is that if you come straight from the Mark 22 serial version then you can expect a speed-up of around 7.6 times.  In the future, when you compare the Mark 22 SMP library to the Mark 23 serial library then you can expect a speedup of around 2.6 times on a 4 core machine like mine.

Example 2 – Quasi random number generation

Quasi random numbers (also referred to as ‘Low discrepancy sequences’) are extensively used in Monte-Carlo simulations which have applications in areas such as finance, chemistry and computational physics.  When people need a set of quasi random numbers they usually need a LOT of them and so the faster they can be produced the better.  The NAG library contains several quasi random number generators but the one considered here is the routine g05ymf.  The benchmark program I used is called bench_g05ymf.f90 and the full set of files needed to compile it are available at the end of this section.

The benchmark program requires the user to input 4 numbers and a boolean as follows.

• The first number is the quasi random number generator to use with the options being:
  1. NAG’s newer Sobol generator (based on the 2008 Joe and Kuo algorithm^[3] )
  2. NAG’s older Sobol generator
  3. NAG’s Niederreiter generator
  4. NAG’s Faure generator
• The second number is the order in which the generated values are returned (The parameter RCORD as referred to in the documentation for g05ymf). Say that the matrix of returned values is called QUAS then if RCORD=1, QUAS(i,j) holds the jth value for the ith dimension, otherwise QUAS(i,j) holds the ith value for the jth dimension.
• The third number is the number of dimensions required.
• The fourth number is the number of number of quasi-random numbers required.
• The boolean (either .t. or .f.) determines whether or not the output should be verbose or not.  A value of .t. will output the first 100 numbers in the sequence.

To compile the benchmark program against the serial library I did:

gfortran ./bench_g05ymf.f90 shtim.c wltime.f /opt/NAG/fll6a22dfl/lib/libnag_nag.a  \
/opt/NAG/fll6a22dfl/acml/libacml_mv.a -o serial_g02aaf

As before, the only thing needed to turn this into a parallel program was to compile against the SMP library and add the -fopenmp switch

gfortran -fopenmp ./bench_g05ymf.f90 shtim.c wltime.f /opt/NAG/fsl6a22dfl/lib/libnagsmp.a  \
/opt/NAG/fsl6a22dfl/acml4.3.0/lib/libacml_mp.a -o smp_g02aaf

The first set of parameters I used was

1 2 900 1000000 .f.

Which uses NAG’s new Sobol generator with RCORD=2 to generate and store 1,000,000 numbers over 900 dimensions with no verbose output.  Averaged over 10 runs the times were 5.17 seconds for the serial version and 2.14 seconds for the parallel version giving a speedup of 2.4x on a quad-core machine.  I couldn’t push number of dimensions much higher because the benchmark stores all of the numbers in one large array and I was starting to run out of memory.

If you only want a relatively small sequence then switching to the SMP library is actually slower thanks to the overhead involved in spawning extra threads.  For example if you only want 100,000 numbers over 8 dimensions:

1 2 8 100000 .f.

then the serial version of the code takes an average of  0.0048 seconds compared to 0.0479 seconds for the parallel version so the parallel version is almost 10 times slower when using 4 threads for small problems. Setting OMP_NUM_THREADS to 1 gives exactly the same speed as the serial version as you might expect.

NAG have clearly optimized this function to ensure that you get a good speedup for large problem sets which is where it really matters so I am very pleased with it.  However, the degradation in performance for smaller problems concerns me a little.  I think I’d prefer it if NAG were to modify this function so that it works serially for small sequences and to automatically switch to parallel execution if the user requests a large sequence.

Conclusions

In the old days we could speedup our programs simply by buying a new computer.  The new computer would have a faster clock speed than the old one and so our code would run faster with close to zero effort on our part.  Thanks to the fact that clock speeds have stayed more or less constant for the last few years those days are over.  Now, when we buy a new computer we get more processors rather than faster ones and this requires a change in thinking.  Products such as the NAG Library for SMP & multicore help us to to get the maximum benefit out of our machines with the minimum amount of effort on our part.  If switching to a product like this doesn’t give you the performance boost you need then the next thing for you to do is to learn how to program in parallel.  The free ride is over.

In summary:

• You don’t need to modify your code if it already uses NAG’s serial library.  Just recompile against the new library and away you go.  You don’t need to know anything about parallel programming to make use of this product.
• The SMP Library works best with big problems.  Small problems don’t work so well because of the inherent overheads of parallelisation.
• On a quad-core machine you can expect speed-ups around 2-4 times compared to the serial library.  In exceptional circumstances you can expect speed-up as large as 7 times or more.
• You should notice a speed increase in over 300 functions compared to the serial library.
• Some of this speed increase comes from fundamental libraries such as BLAS and LAPACK, some of it comes from NAG directly parallelising their algorithms and yet more comes from improvements to the underlying serial code.
• I’m hoping that this Fortran version is just the beginning.  I’ve always felt that NAG program in Fortran so I don’t have to and I’m hoping that they will go on to incorporate their SMP improvements in their other products, especially their C library and MATLAB toolbox.

Acknowledgements

Thanks to several staff at NAG who suffered my seemingly endless stream of emails during the writing of this article.  Ed, in particular, has the patience of a saint.  Any mistakes left over are all mine!

Links

• The Nearest Correlation Matrix benchmark program used here – bench_g02aaf.f90 – along with the source code for the timing routine.
• NAG’s documentation for the g02aaf routine.
• The Quasi Random Number benchmark program used here – bench_g05ymf.f90 – along with the source code for the timing routine.
• An Introduction to quasi random numbers (written by a member of NAG’s staff)
• Official NAG documentation for g05ymf – the quasi random number generator considered here.

References

[1] – Higham N, Computing the nearest correlation matrix—a problem from finance, IMA Journal of Numerical Analysis 22 (3): 329–343

[2] – Qi H and Sun D (2006), A Quadratically Convergent Newton Method for Computing the Nearest
Correlation Matrix, SIAM J. Matrix AnalAppl 29(2) 360–385

[3] – Joe S and Kuo F Y (2008) Constructing Sobol sequences with better two-dimensional projects SIAM J. Sci. Comput. 30 2635–2654

Dislaimer: This is a personal blog and the opinions expressed on it are mine and do not necessarily reflect the policies and opinions of my employer, The University of Manchester.