The NAG Blog

Cloud computing or HPC? Finding trends.

2012-01-19T18:22:00.001Z

Enable innovation and efficiency in product design and manufacture by using more powerful simulations. Apply more complex models to better understand and predict the behaviour of the world around us. Process datasets faster and with more advance analyses to extract more reliable and previously hidden insights and opportunities.

All ambitions that will probably resonate with those seeking scientific advances, commercial innovation, industrial growth and more cost-effective research. Underpinning all of the above is the use of more powerful computing methods and technologies. Faster and more capable computers - but equally important - more advanced and better performing algorithms and software implementations.

It's a pretty convincing story for those who take the time to listen - whether business leaders, governments, or research funders. Even in these challenging economic times, it has led to investments from industry and governments for this reason - the potential return is well documented and significant. It is even enticing enough to interest the media and the public - especially when we use emotive descriptions like "world's fastest supercomputer", "international competitiveness in digital economy", "personal supercomputing", and so on.

And it is this last thought that cause me to diverge from the grand theme to explore names and attention. I will come back to the main theme later (a future blog), as it is both important and timely. But on to my side topic.

Anybody using, selling or funding the technologies and methods described above will know what I am talking about. But the names and labels applied to it can vary significantly across the diverse audience. High performance computing (HPC), supercomputing, computational science and engineering, technical computing, advanced computer modelling, advanced research computing, etc. The range of names/labels and the diversity of the audience involved mean that what is a common everyday term for many (e.g. HPC) is an unrecognised meaningless acronym to others - even though they are doing "HPC".

This can create a barrier to engaging politicians, companies that could benefit, the media, and people in search of solutions for their day-to-day modelling/simulation/data processing challenges.

Let's play.

Most of us who see this as part of our daily life use the terms HPC or supercomputing. How do these stack up with the wider world? Let's turn to Google Trends as an arbitrary tool of statistics.

The following graph shows the search popularity of these terms ("supercomputer" and "HPC") over the last few years. Clearly "HPC" is a more common keyword.

[Plot 1: blue = supercomputer, red = HPC]

But what if we add in that term that is so often considered just a buzzword by seasoned HPC professionals - "cloud computing"?

[Plot 2: blue = supercomputer, red = HPC, orange = cloud computing]

We see that in the last few years "cloud computing" has soared above the traditional names in usage. This means a wider audience - and thus more possibilities for that ambitious opening paragraph of mine.

Adding some more technical terms ("parallel computing", "multicore") hardly register in comparative popularity.

[Plot 3: blue = supercomputer, red = HPC, orange = cloud computing, green = parallel computing, purple = multicore]

Interestingly, adding a domain specific term ("CFD") tracks the popularity of HPC rather than cloud computing.

[Plot 4: blue = CFD, red = HPC, orange = cloud computing]

You can play the same game with key technologies of the supercomputing world - e.g. [MPI, OpenMP, CUDA, OpenCL, Fortran] - and discover more interesting trends, but as this blog is already getting long - that is for another day.

I'll just leave you with this one, which might be interpreted as speaking volumes to the challenges faced in delivering the promise of my opening paragraph - [computer, software, programmer, algorithm].

[Plot 5: red = computer, blue = software, orange = programmer, green = algorithm]

What interesting related trends can you find and analyze?

Self-Improvement

2012-01-05T18:09:00.000Z

Winston Churchill once said "The pessimist sees difficulty in every opportunity. The optimist sees the opportunity in every difficulty." (for you confirmed pedants, it may have been L.P. Jacks)

My custom is to use the time away from work at the end of the year to think about what I want to do differently in the year ahead. Among the topics that came up was e-mail, the bane of my life and perhaps yours as well. I get hundreds every day (and that doesn't include the SPAM).

Being a lifelong optimist I've decided to make e-mail less of a pain in my life, both work and personal. Being an occasional realist, I recognize that I have a limited number of options and they must focus on what I can do.

So, here's my plan for 2012 (with acknowledgement to Scott Belsky and Stever Robbins who supplied some of the ideas and got me thinking). See Disrupt Your Inbox and What you should never say in an e-mail

Experiment with three-sentence emails when I need an answer from someone. (improve the likelihood that someone will actually read and answer).
Start e-mails with action I want, Don't leave the reader guessing until the end.
Use subject lines that intrigue the reader and actually invite them to open and read the e-mail.
Take disagreements offline. There are volumes that could be written on this.
Don't "reply all" unless everyone needs to be involved. How do you feel when you are one of 92 people copied on an e-mail that doesn't interest or pertain to you? Resist the urge, it's probably illegal.
If you need to make several points and expect a response, use numbers for reference to reduce length, opportunity for confusion.

What's your plan?

More on managing e-mail in a future installment.

Question one: Where's my phone?

2012-01-03T17:20:00.000Z

The end of the year at NAG is always celebrated with a rather splendid Christmas lunch, which is generously paid for by the company as an acknowledgement of the hard work its employees have put in over the previous twelve months. Accordingly, it provides an occasion for some much-needed relaxation and refreshment before the Christmas break. It's also the time for the NAG Christmas Quiz which, whilst not necessarily contributing to the participants' relaxation, usually provides some entertainment or diversion for those who care about such arcane matters as the number of hearts an octopus has or Paul McCartney's middle name. One of those tortured souls is the present author, who cleverly realized some time ago that the only way to be sure of knowing all the answers to the questions was to set them.

Setting quiz questions in the connected age - particularly for a technical-savvy band such as the employees of NAG - can present a few challenges, however. For example, anyone with internet access (via, say, a smartphone's web browser) would be able to find the answers to the questions indicated above in a matter of moments. More direct questions, such as

Who wrote "A Child's Christmas In Wales"?
What kind of logs did Good King Wenceslas ask for?
When was "Merry Xmas Everybody" number 1 in the UK?

are even easier to answer, though some skill and judgement may still be required in the selection of the correct response to an ambiguous or ill-posed question such as the first one on this list (does it refer to Dylan Thomas's prose piece or John Cale's song?). Whilst this is clearly a valid and imaginative use of technology, I wondered whether it would provide an unfair advantage over those participants who wouldn't be using their phones in this fashion (or texting more knowledgeable friends for answers) and started wondering about ways to obviate their effectiveness.

I realized fairly quickly that confiscating all phones on entry to the restaurant didn't fall within my powers as quizmaster (chiefly because I didn't have any) and, for a similar reason, there were no funds in the quizmaster's budget for the purchase of a mobile phone jammer. My suggestion that we should move the location of the Christmas lunch to the interior of a Faraday cage in order to attenuate the phone signal wasn't looked on too kindly by senior management (or the restaurant) either. Accordingly, I began thinking about using questions that were more indirect, which might make searching for an answer more difficult. For example: what's the connection between these three things?

Queen Jezebel
The Regents of Prague
Chopin's piano

This type of query is - very roughly - analogous to a so-called inverse problem in science, in which we're asked to use observed data or results to deduce something about an underlying system or model. It can be harder to answer because there may be several models that are consistent with the observations - for example, one thing that connects the three things is that each contains the letter 'n', but that's not necessarily the right answer (which is ever-so-vaguely computer related, lest I be accused of straying too far off-topic).

Another type of problem that might present more challenges in the search for an answer is image recognition. Whilst mobile tools such as Google Goggles are already trying to make this easier for specific examples (e.g. in the identification of labels and landmarks), it's still either very difficult or impossible to recognize arbitrary objects. In the context of the quiz, this means questions such as: identify these people:

or name these films:

I'm not sure how effective my preparations were, but the quiz appeared to provide the usual amount of stimulation, along with a certain degree of exaltation and frustration (only for those who care about this sort of thing, naturally). Some of the participants have even started talking to me again, but that might be just the supervening effect of the Christmas break, and all the goodwill-to-all associated with that happy season.

PS If you'd like to answer any of the questions above, please feel free to add a comment below, although I regret to say that the unbelievably fabulous prizes have long been consumed by the winners of the quiz. Don't use your phone, though.

Coffee and Filters

2011-12-01T13:42:00.001Z

Given the origins of NAG and our mission, it’s natural for us to take the “long view” in giving back to the communities in which we operate. In the US, we’ve just finished our second year as sponsor of the DemandTec Retail Challenge (DTRC) scholarship competition for high school students in the Chicago area. DemandTec is one of our earliest and strongest software company partners who incorporate NAG components into their software products, providing demand management software for major retailers around the world.

The DTRC puts 2-person teams of high school seniors in the role of category managers for a retail store in a a 2-week computer simulation where each day represents a week in the “real world”. The students are given many weeks of data showing the price, inventory, unit sales, promotions used and the profit earned on each product they are managing. In the contest, they are responsible for two brands of coffee, one brand of tea and coffee filters. During the competition, they must analyze prior results to set the new price of each item, decide whether to run promotions and decide how much inventory to purchase. The simulation creates an interaction both with consumers and with the other teams in the competition. Ultimately, the three teams with the highest profit at the end of two weeks advance to the regional finals here in Chicago.

At the Chicago Regional finals held November 10th, teams from Glenbard West High School and Wheaton Warrenville South High School gave presentations of their data analysis, strategy formulation and how they worked as a team as they adapted their strategy in the course of the contest. The judges for the contest were experienced professionals in retail analytics from DemandTec and other local companies. Our winners were team “Price Lords” consisting of Peter Ericksen and Greg Grabarek. Each of them received $2500 toward their colleges expenses next year from NAG and our co-sponsor software company Informatica. Our third sponsor was the Network of Executive Women (NEW), which supports the education and advancement of women in leadership roles.

This contest is a labor of love for those of us involved. We have cultivated relationships with teachers of mathematics and statistics at local schools and have made presentations on the contest to a number of high school classes. It not only engages us with the students and teachers but it also helps us educate the next generation of applied scientists who might one day be NAG users. We wish Peter and Greg good luck as they compete in the national semi-finals in early December and are already making plans or next year’s contest.

Calling NAG routines from R

2011-11-23T11:21:00.001Z

R is a widely-used environment for statistical computing and data analysis. It is one of the modern implementations of the S programming language (i.e. much code written in S will run unaltered in R) although its underlying semantics are derived from Scheme. R is Free Software.

The capabilities of R can be extended through the use of add-on packages and a large number of these are available from the Comprehensive R Archive Network (CRAN). Some users have expressed an interest in calling NAG Library routines from within R; accordingly, we have recently created a package which provides access to some of NAG's numerical functionality. More specifically, the NAGFWrappers R package contains the local and global optimization chapters of the NAG Fortran Library, together with a few nearest correlation matrix solvers and some other simpler routines. The package incorporates documentation in the so-called Rdoc style that will automatically produce help pages within the R system (see figure below), and also in HTML and PDF - for example, here is the full list of NAG routines that are contained in the package.

Help for the NAG routine e04uc, as displayed in R

For completeness, and to help R users further, we have also published more general instructions about how to use the R extension mechanisms to access any NAG routine from within R.

The original version of NAGFWrappers has been available since mid-2011; we have just updated it to use Mark 23 of the Fortran Library, and are releasing R binary packages for Windows 32 bit and Windows 64 bit, along with the R source package which can be used on other platforms (for example, we have built and run it on 64 bit Linux).

It should perhaps be noted that this is a preview release of the package, which is aimed at obtaining user feedback. Although it has been built and run on the platforms mentioned above, it is not a NAG product. We are keen to receive user feedback, and will respond to technical queries and problem reports via support@nag.co.uk so that we can further refine this package and make it still more useful to the R community.

SC11 diary catch up

2011-11-01T18:41:00.000Z

I posted here a week or two ago about my diary leading up to the year's biggest supercomputing event - SC11 in Seattle. I though it would be handy to give a quick summary of the diary entries so far for those who haven't been reading along.

If you recall, I said: "On my twitter stream (@hpcnotes), I described it as: "the world of #HPC in one week & one place but hard to find time for all of technical program + exhibition + customer/partner/supplier meetings + social meetings + sleep!" To follow the news about SC11 on twitter, follow @supercomputing and/or the #SC11 hashtag."

"Any hope of "live" blogging or actively tweeting during the week of SC11 itself is almost nil - the week is just too busy with the day job. Even simply consuming the relevant news, comment and gossip is a challenge."

"So instead I am going to try to write a diary of the lead up to SC11."

If you've missed them, here are the 8 SC11 blogs so far:

Along the way, I have briefly alluded to a few things NAG will be doing at SC11. One of my colleagues will be along shortly to post here about our activities at SC11, but in the meantime, plan to visit us on booth 2622, or get in touch to arrange a conversation.

SC11 will be here soon

2011-10-20T16:52:00.002+01:00

The SC11 conference, or just "supercomputing", will be held in Seattle this November. For many in the high performance computing community, including NAG, SC is the big event of the year. Certainly it is the one that attracts the most press (and press releases), the most attendees, the biggest exhibition, and absorbs the most amount of time in planning before we even get there. It is the event where we get to meet with many of our customers, most of our potential suppliers, and many friends and collaborators.

SC11 will generate a lot of news, articles and other media in the coming weeks. You can keep up with SC11 on twitter by following @supercomputing and/or the #SC11 hashtag. NAG's twitter team (@NAGtalk) will also keep you up to date with NAG's activities and stories at SC11 and beyond.

The NAG team will also post here on The NAG Blog over the coming days and weeks with news about what NAG will be doing at SC11.

I will be adding to all this with a series of posts on my personal blog - a diary of the lead up and planning for SC11 - the first entry was posted today.

See you in Seattle.

Calling routines from the NAG Fortran and C libraries from within LabVIEW

2011-10-06T15:02:00.024+01:00

This is a follow-on from my colleague Sorin Serban's recent post, in which he looked at invoking methods from the NAG Library for .NET from within the LabVIEW programming environment. The motivation for this work was the realization that many LabVIEW users want to enhance their applications by making use of some of NAG's numerical routines. Recently, I've been using our Fortran and C libraries to do the same thing, and some of this work is discussed here.

The mechanism used within LabVIEW to call a routine from either of these libraries is different from that used in the case of the .NET Library. This is because that library is a so-called .NET assembly, whilst the C and Fortran libraries are (on the Windows platform) dynamic link libraries (DLLs). The .NET assembly uses the common language runtime (CLR) and the .NET framework to manage assembly functions, and to export information about classes, methods, properties and events. The practical implication of this is that when a method from the NAG Library for .NET is loaded into LabVIEW, information about its function arguments is automatically loaded as well. This is not the case with a DLL, as we shall see.

As previously, we begin building a LabVIEW application by creating a blank virtual instrument (VI) using, for example, the Blank VI option on LabVIEW’s Getting Started splash screen, and go to the block diagram, which is where our routine will be loaded. Once again, we build our example around the NAG routine g01aa, which performs simple statistical calculations on a set of ungrouped data; besides being implemented as a method in the .NET Library, it is also available in both the Fortran and the C libraries.

To load the routine, we first bring up the Functions Palette, using the option on the View menu on the block diagram window. Next, we open the Connectivity collection from the palette, and then the Libraries & Executables item. This palette contains a variety of functions for calling code from libraries and executing system commands, and we select the Call Library Function item – which calls a DLL directly – and drag it onto the block diagram in order to initialize it. To specify the library we want to call, we right-click on the node and select Configure.... This brings up the Call Library Function dialog box, and we enter the name of the DLL in the Library name or path box. In general, either the name or the path is acceptable, although the choice between between them has consequences if the VI is to be distributed or relocated. In our example, we're running LabVIEW under Windows, and are using the implementation of the NAG Fortran Library contained in FLDLL224M_nag.dll and the NAG C Library implementation in CLW3209DA_nag.dll.

Having specified the library, we next select the name of the routine we wish to call by picking it from the drop-down list in the Function name dialog box; in our example, this is G01AAF if we're using the NAG Fortran Library, or g01aac if it's the NAG C Library that we're calling. In either case, we specify the calling convention as stdcall (WINAPI).

The next step is to tell LabVIEW about the argument list for our routine via the Parameters tab on the Call Library Function dialog box. For each parameter, we specify its name, type and how it is to be passed (as a value or a pointer to a value) to LabVIEW when it calls the routine. As these details are entered, LabVIEW updates the Function prototype display at the bottom of the dialog box with the C prototype of the routine (LabVIEW is written in C, and it is from this language that the external routine is called). This table shows the parameters for G01AAF, including their Fortran type, their LabVIEW specification and the C type that LabVIEW generates for the function prototype:

Fortran variable name	NAG Fortran type	LabVIEW Type	LabVIEW Constant?	LabVIEW Data type	LabVIEW Pass	LabVIEW Array Format	LabVIEW C type
N	INTEGER	Numeric	Y	Signed 32-bit Integer	Pointer to Value	-	const int32_t*
X	double precision array	Array	Y	8-byte Double	-	Array Data Pointer	const double*
IWT	INTEGER	Numeric	N	Signed 32-bit Integer	Pointer to Value	-	int32_t*
WT	double precision array	Array	N	8-byte Double	-	Array Data Pointer	double*
XMEAN	double precision	Numeric	N	8-byte Double	Pointer to Value	-	double*
S2	double precision	Numeric	N	8-byte Double	Pointer to Value	-	double*
S3	double precision	Numeric	N	8-byte Double	Pointer to Value	-	double*
S4	double precision	Numeric	N	8-byte Double	Pointer to Value	-	double*
XMIN	double precision	Numeric	N	8-byte Double	Pointer to Value	-	double*
XMAX	double precision	Numeric	N	8-byte Double	Pointer to Value	-	double*
WTSUM	double precision	Numeric	N	8-byte Double	Pointer to Value	-	double*
IFAIL	INTEGER	Numeric	N	Signed 32-bit Integer	Pointer to Value	-	int32_t*

Although there are several parameters, many of them have the same characteristics. Firstly we recall that, since Fortran passes arguments by reference (not by value), each parameter must be passed as a pointer when calling the Fortran routine from LabVIEW (i.e. from C). Next, the documentation for this implementation of the NAG Fortran Library tells us that double precision means DOUBLE PRECISION - that is, an 8-byte floating point number, whilst the Fortran INTEGER type is a 32-bit integer. Hence, these two types can be respectively mapped to double and int32_t, via LabVIEW's 8-byte Double and Signed 32-bit Integer types. Finally, it can be seen how designating a variable as a constant (because it's only used on input to the routine) results in LabVIEW applying the const qualifier to the C type.

Similar remarks apply, for the most part, to the use of the g01aac routine (note that the parameter list is slightly different for this version of the routine):

C variable name	NAG C type	LabVIEW Type	LabVIEW Constant?	LabVIEW Data type	LabVIEW Pass	LabVIEW Array Format	LabVIEW C type
n	Integer	Numeric	Y	Signed 32-bit Integer	Value	-	const int32_t
x	const double array	Array	Y	8-byte Double	-	Array Data Pointer	const double*
wt	const double array	Array	Y	8-byte Double	-	Array Data Pointer	const double*
nvalid	Integer*	Numeric	N	Signed 32-bit Integer	Pointer to Value	-	int32_t*
xmean	double*	Numeric	N	8-byte Double	Pointer to Value	-	double*
xsd	double*	Numeric	N	8-byte Double	Pointer to Value	-	double*
xskew	double*	Numeric	N	8-byte Double	Pointer to Value	-	double*
xkurt	double*	Numeric	N	8-byte Double	Pointer to Value	-	double*
xmin	double*	Numeric	N	8-byte Double	Pointer to Value	-	double*
xmax	double*	Numeric	N	8-byte Double	Pointer to Value	-	double*
wsum	double*	Numeric	N	8-byte Double	Pointer to Value	-	double*
fail	NagError*	Numeric	N	Signed 32-bit Integer	Value	-	int32_t

Since we're now calling a C routine from C, there is a closer correspondence between the NAG type and the LabVIEW C type, although a couple of additional explanatory remarks might be necessary. Firstly, the NAG C Library uses a data type called Integer; for this implementation, this is defined as long, which is a 32-bit integer on Windows systems. Secondly, the fail variable is a pointer to a C structure named NagError; we access this using a 32-bit integer (not a pointer to an integer), which is the same size as the address of the structure.

More information about how to create an interface between LabVIEW and external libraries (including advice for more complicated cases than those used here) is available here. Having created the interface to the NAG routine, the final step is to add LabVIEW controls and indicators to handle the input and output of our application. This is exactly the same procedure as was described previously, and we conclude this post with Figure 1, which shows the block diagram of our application, centered around the Call Library Function, which is where the NAG routine (in this case, g01aac) is being called.

Figure 1. The g01aac routine in the LabVIEW block diagram

NAG communications

2011-09-08T14:35:00.003+01:00

We know that it can be hard for you to keep track of all the new stuff that is happening at NAG, such as the regular Library updates, new ways to use routines from various computing environments, events and conferences, training days, talks at customer sites, collaborations with numerical specialists, work on new technology platforms, and so on. The list of new activities and developments can seem endless, but we try to directly communicate the information that seems likely to be of most interest to you.

For example:

If you are a NAG user, we let you know when a new release of the Library is available for your system, so that you have the quickest access to new routines.
When we know your field of interest, we send you information about events and highlight relevant examples in order to help you get the most from NAG.
We periodically release brief case studies that contain novel ways of using NAG routines; some of these might help trigger new ideas for your problem solving.

However, one of the best ways to keep up with all the information is via NAGnews. This e-newsletter contains articles on new product development and technical advancements along with developer tips and hints and event information. More than 7,000 subscribers already receive NAGNews, which is distributed approximately every 6 weeks. Registration is free of charge, and NAG will not share your details with any third party.

NAGNews has been going for some time – in fact, the next issue will be the hundredth. Why not take a look here if you’d like to consider asking to receive a copy?

How do you use NAG in Excel?

2011-08-22T15:27:00.005+01:00

For years NAG has considered Excel users as a group with certain needs for good numerical algorithms. Nowadays Excel is the industry-standard for spreadsheets, but it is not a statistical package nor is it a numerical analysis tool. That’s why Excel users come to NAG and ask how our accurate, robust and efficient NAG Library can be utilised in Excel.

There are a few ways in which the Library can be called in Excel:

xls/xlsm spreadsheet and VBA

The most popular, widely used way: NAG and Excel page describes how to do it with the help of technical papers and downloadable examples. Basically you create a new spreadsheet and add VBA code to it. NAG has created VBA Declaration Statements which provide the bridge between NAG and Excel- the user just has to add the declaration file to a VBA module and the NAG DLL becomes visible to Excel & VBA. From this moment on you can create macros that will, for example, read data from an external data file, process it in VBA and then display the results on the spreadsheet. Or you can create a UDF (User Defined Function) that is called straight from the spreadsheet, but inside uses NAG to compute some statistics, find a zero of a function, or find a value of a European call option using Black-Scholes Model or Heston’s Stochastic Volatility Model. This is the great thing about using NAG in conjunction with Excel- you can do it in many ways!

xla/xlam add-in and VBA

Very similar to the previous point. The difference is that the VBA code for the macros and UDFs is in an Excel add-in. As a result the user’s data and the user’s code may be in two separate files. So the user needs to load the xla add-in while using his xls file to have the functions available. http://www.cpearson.com/excel/createaddin.aspx nicely describes how an xla add-in works and how to create it.

xll add-in

An xll is a different type of add-in, it actually allows you to use your C/C++, Fortran, or VB code from Excel. This is actually a Dynamic Link Library (DLL) file that is compiled and accessible via Excel API. The benefit of having an xll is that it’s compiled and will be faster than using VBA. However it requires very cautious programming due to memory issues and is generally harder to program.

VSTO

Visual Studio Tools for Office is a set of development tools that allows creating Visual Studio projects in C# or VB.NET in the form of Excel add-ins that extend Excel’s features or customise the user interface. Since this is a .NET environment, it’s advisable to use the NAG Library for .NET.

We look forward to hearing from our users and anyone that is interested in using Excel in conjunction with NAG. How do you use Excel? What example programs do you think we should have on our website? Currently we have many examples on using NAG via VBA, but not many on the other methods.

Tell us what you would find useful!

What is this HPC thing?

2011-08-19T22:38:00.003+01:00

I’m sure something like this is familiar to many readers of this blog. The focus here is HPC, but there is a similar story for mathematicians, numerical software engineeers, etc.

You've just met an old acquaintance. Or a family member is curious. Or at social events (when social means talking to real people not twitter/facebook). We see that question coming. We panic. Then the family/friend/stranger, asks it. We freeze. How to reply? Can I get a meaningful, ideally interesting, answer out before they get bored? What if I fail to get the message across correctly? Oops, this pause before answering has gone on too long. Now they are looking at me strangely. They are thinking the answer is embarrassing or weird. This is not a good start.

The question? “What do you do then?” Followed by: “Oh! So what exactly is supercomputing then?”

The problem is that it usually takes a several slides or a few minutes of explanation to give a decent overview to a scientist who is new to HPC. Your questioner is almost certainly not a scientist, and maybe even thinks using email, facebook or twitter is a big computing achievement.

So, where do you start? The big computer? The science it enables? The money involved? You want to quickly convey the scale, the enormous range of user disciplines, the specialist nature, the benefits to everyday life. You want to distinguish it from laptops, tablet computers, and corporate servers. You want to avoid mentioning nuclear weapons labs. You’d even like to come across as cool, but it’s way too late for that.

It’s a question that I have seen leaders in the field of HPC struggle to answer effectively - either talking rubbish, stalling, or starting a monologue on computer science theory. But really, especially in the current tough funding climate, it is something we should all learn to answer properly. Be able to explain quickly, simply – and relevantly – what HPC is all about. Make it sound the interesting, high paced, leading edge field that it is.

Why can’t HPC have the same positive media grab that particle physics and astronomy have? If the mass media feel they can talk about particle physics, radio telescopes, space telescopes, or sun storms, then we should be able to make supercomputing amenable to the same audience.

We know that many areas of research and industry can benefit hugely from High Performance Computing. And so it is important to the economy and society to have a strong HPC future. But for a strong future, HPC needs millions of HPC-aware students and public – maybe more than it needs thousands of extra HPC skilled programmers. The latter will follow naturally if the former is done well. It is more essential to ensure the benefits of HPC are recognised by end user organisations for future generations.

So, next time you are asked by family or strangers what you do, avoid the pause followed by a shy “computer stuff” – instead answer “supercomputing” and then quickly explain what that is and why it matters.

And, if you can get the message across positively to family and strangers, then as well as the long term benefits above, you will find you have the skills to get it across positively to your funders, politicians, senior managers, board, etc. when the situation arises.

Seeing the good for the trees

2011-08-05T09:43:00.013+01:00

We put a lot of effort into trying to ensure that our numerical libraries are available on the platforms that are most popular with our users. For each library, this leads to a proliferation of implementations, each of which is targeted at a specific combination of compute processor, operating system and compiler. The details of current implementations are on our website - for example, here is the list for the NAG Fortran Library, which also includes download links for each implementation. Although this list - which currently features forty-nine implementations - is an impressive array (and the fruits of a great deal of careful work on the part of our implementation team), its presentation could perhaps be viewed as being tricky to navigate by those who are searching for the appropriate implementation for their particular system.

Recently I was asked if I could produce a more intelligible representation of pages like this. I was reminded initially of a decision tree, mainly because we already make use of these tools in our documentation as a way to help users select the most appropriate algorithm (or NAG routine) for the solution of their problem. For example, figure 1 shows one of the decision trees which is part of the documentation for chapter D01 of the NAG Fortran Library, which deals with numerical integration, or quadrature.

Figure 1. Decision tree to help users select the appropriate NAG quadrature routine

This chapter currently contains twenty-nine routines, each of which is tailor-made for a different type of problem (e.g. integrands defined as a set of points or a function, functions that contain singularities, integration intervals that are finite, semi-infinite or infinite, etc), and the decision tree guides the user to a correct choice by asking questions about their problem.

Recasting the implementation list as a decision tree (with questions such as "what platform?", "what operating system?", "what compiler?") seemed, then, to be the way forward. However, I was conscious that I didn't want something which was as as set in stone as figure 1; for one thing, I was unsure about the initial layout of my tree. I also knew that the tree would require maintenance as new implementations became available - for example, the latest release of the NAG Fortran Library is currently being made available on more platforms, and the picture (and the list) will have to be periodically updated to reflect this. In addition, because the tree only contains branching nodes (points at which the path diverges) and leaf nodes (points at which the path terminates) but no converging nodes (points at which several paths come together), it could quickly become very large, and difficult to draw by hand.

Fortunately, there are a range of options for the display of decision trees. Some of these - for example, this one and this one - also incorporate the calculation of the tree itself (for use in machine learning applications), but this isn't necessary in the present example because we already have the structure of our tree. The display of trees and graphs is an active research area (see, for example, here) and there is a variety of freely-available applications (for example this one and this one) which can be used. In the end, I chose prefuse, which is a Java-based toolkit for building interactive applications in information visualization, chiefly because their treeview demo appeared to be close to what I wanted to do.

I downloaded the prefuse distribution and, after some only limited success with the supplied build script, followed the advice of the documentation and built it from within the Eclipse IDE (more specifically, their advice was to use Eclipse "if you are a Java novice", which certainly rang a bell with me). Adapting the treeview demo to my purposes was very straightforward because the tree structure and node names are read in as an XML file whose format is very simple (there are essentially just three elements: tree, branch and leaf; these are labelled via attributes). Figure 2 is a screen shot of the treeview demo running as a Java applet in the browser, displaying the implementation information for the NAG Fortran Library (i.e., based on the information presented here).

Figure 2. Screenshot of the treeview demo, running as an applet in the browser

The instructions at the bottom of the web page give the details of the applet's interface: how the user can pan, zoom and also change the orientation of the tree from its default left-to-right ordering. Clicking on any of the nodes expands the tree from that point (and collapses any expanded node at the same level) which facilitates navigation. Figure 3 shows the applet being used to determine the appropriate implementation of the NAG Fortran Library for applications built on machines having the Sun X64 architecture, and using the Sun Studio 11 compiler; its NAG product code is FLSA622DCL.

Figure 3. Using the applet to select an appropriate version of the NAG Fortran Library for Sun

In a similar fashion, Figure 4 is a partial screenshot of the applet being used to select the correct implementation of the NAG Fortran Library for applications built on machines having the AMD64 architecture, and using the GNU Fortran compiler with the ILP64 data model; its NAG product code is FLL6A22DHL. Note that, in this screenshot, the orientation of the tree has been switched to a top-to-bottom ordering, which might be easier to navigate for some users.

Figure 4. As figure 3, but now for Linux

The applet is currently being tested internally at NAG, but initial feedback seems to indicate that its presentation of the implementation data represents a helpful supplement to the list format which has been in use up until now.

Types in The NAG Toolbox for Matlab®

2011-07-14T20:26:00.011+01:00

I was recently asked why the NAG Toolbox for Matlab requires users to be explicit about the types of parameters passed to NAG routines, for example:

[vr, info] = f08nf(int64(1), int64(4), vr, tau);

and similarly why output parameters like info in this example have class int64 as opposed to double.

The original reason for this was efficiency. The NAG Toolbox consists of wrappers to NAG Library routines and there are always issues in transforming data into the correct format when calling between different languages. As well as the top-level call (like the one shown here), we also have to transform data when calling Matlab functions in "call backs", for example when evaluating the objective function during the solution of an optimisation problem. In such cases there may be many thousands of calls and so we want to make the data transformation as efficient as possible. It was for this reason that we decided to put the onus on the user to do any coercions and checks on the validity of those coercions. When we beta tested our original prototype the feedback we got was that this wasn't a big deal.

Back then, all the integers used in the NAG Toolbox were 32 bit, even on 64 bit machines. This was because the version of LAPACK that Matlab used was 32 bit and we needed to be compatible with that since we were calling those routines ourselves. In Matlab, 32 bit integers behave quite nicely: you can perform arithmetic on them and use them as the bounds in loops so users were unlikely to notice the types of objects that we returned. However when Matlab switched to 64 bit integer versions of LAPACK we had to follow, and that's when things started to get harder. Back then it wasn't possible to do arithmetic on 64 bit integers in Matlab or to loop over them, and so programs to call our routines and their associated call-backs started to look quite clumsy with lots of casts between int64 and double. At this point we started to consider whether to accept and return doubles and do the coercions transparently inside our wrappers.

Unfortunately there is a problem with this approach, in that not all 64 bit integers can be represented as doubles. For example consider the following Matlab session:

>> x=2^53
x =
 9.0072e+015
>> x+1-x
ans =
    0
>> n=int64(2)^53
n =
    9007199254740992
>> n+1-x
ans =
                   1

Since there are only 53 bits available to represent the mantissa in a double precision number, once we get beyond 2⁵³ we can't represent every integer exactly. For example in this case we can represent x and x+2, but not x+1. In many cases integers returned by NAG routines will never be big enough for this to be a problem but there are places where we cannot be sure of this. The good news is that recent versions of Matlab support arithmetic on 64 bit integers although you still cannot loop over them, so the problem has been alleviated but hasn't quite gone away.

There is another issue with having versions of the Toolbox that accept different kinds of integers, which is that it's harder to write portable code. In the next version of the Toolbox we will provide two functions nag_int and nag_int_name to allow users to perform operations like:

[vr, info] = f08nf(nag_int(1), nag_int(4), vr, tau);
a = zeros(100, nag_int_name);

which will work correctly on both 32- and 64-bit toolboxes.

As explained above, we can't always coerce a 64 bit integer to a double so we have no option but to pass integers to call-backs and return integers from calls to NAG routines. However we could give users the option of passing doubles in the place of integer parameters to NAG routines, and returning doubles from call-back functions. With this approach the user could avoid explicit coercions (except when using a 64 bit integer as the end point of a loop) if they wished, but they would incur a performance penalty, which would be noticeable if arrays were involved (both in terms of execution speed and memory usage). While we aren't planning this for the release later this year we would consider it for a future version if users wanted it. I'd love to hear your thoughts.

Why migrate from legacy systems?

2011-07-12T14:09:00.001+01:00

I was tempted by a colleague into posting an answer to a question on LinkedIn. I thought it might be useful to reproduce that here for a potentially different audience ...

Q: (from LinkedIn)
"Why is it important to migrate from legacy systems? Is there anyone who can simplify the migration process?"

A: (my answer, there are others too)

Why migrate from legacy systems? There are two aspects to consider - protection and opportunity.

On the protection side, it is a balance of risks - the future maintainability risk of staying with an older system against the inevitable pains of migration. The old saying "if it ain't broke, don't fix it" is relevant, but the catch is that the further your legacy system falls behind current practice, when it finally does break the cost of the eventual fix is likely to be much greater.

On the opportunity side, legacy systems may well be constrained in terms of performance, capabaility, features, etc. A more modern system may thus provide significant business benefits from better performance, more powerful capabilities, or more relevant features.

In the end, if you are using a different system (e.g. legacy) to most of your industry colleagues/competitors, then you are either operating at an advantage (e.g. cost, established solution, ...) or disadvantage (e.g. cost, performance, capability, ...). Make sure you know which.

As to who can help - it depends on what migration you need (from which system? migrate the application or recreate the functionaility? etc.)

Readers of this blog will know that NAG provides services, training and consulting to support migration of application software (especially science & engineering, or involving parallel computing - multicore, supercomputers, etc).

You can also follow NAG on LinkedIn by joining our group.

Using NAG .NET methods in LabVIEW

2011-07-08T14:01:00.025+01:00

Recently we were asked by users at a few European universities about how to call the NAG Library from within the popular LabVIEW programming environment. They're interested in doing this because they wanted to supplement the functionality of their LabVIEW applications using some of the numerical routines offered by NAG. We already have some preliminary material on our website which describes a few aspects of using the NAG Fortran Library in LabVIEW, but alternative approaches (which could be viewed by some as reflecting more modern programming paradigms such as object orientation) are also possible. Specifically, we’ve been taking a look at calling routines from the NAG Library for .NET from LabVIEW, and we describe some of this work here.

We begin building a LabVIEW application by creating a blank virtual instrument (VI) using, for example, the Blank VI option on LabVIEW’s Getting Started splash screen. This brings up two windows: a front panel and a block diagram. The former contains the user interface of the VI which the user creates using controls and indicators, each of which is selected from a palette of options. Code to control the objects in the front panel is added by the user to the block diagram; this is a manifestation of LabVIEW’s visual programming paradigm, which sees an application as being made up of a collection of nodes that are connected together in a way which reflects the flow of data between them.

Our simple example application is built around the NAG method g01aa, which calculates the mean, standard deviation and various other statistical quantities from a set of ungrouped data (details about the documentation for the NAG Library for .NET is available here, while the documentation for the g01aa routine in the NAG Fortran library is here). To load the method, we first bring up the Functions Palette, using the option on the View menu on the block diagram window. Next, we open the Connectivity collection from the palette, and then the .NET item. This palette contains a variety of functions for creating .NET objects, setting properties or calling methods on those objects. We select the Invoke Node item – which invokes methods or performs actions on an object – and drag it onto the block diagram in order to initialize it.

To specify the class on which to execute the property, we right-click on the node and select Select Class >> .NET >> Browse.... This brings up the Select Object From Assembly dialog box. If the NAG Library is not listed in the Assembly component, then we find it (typically in a location like C:/Program Files/NAG/NAG Library for .NET/) using the Browse... component, choosing either NAGLibrary32.dll or NAGLibrary64.dll, depending on our machine architecture. Double-clicking on the NagLibrary item in the Objects component on the dialog box then shows the objects which are available in this assembly; in our example we select G01, followed by OK to load it into the Property Node. To select the desired method from within the G01 class, we right-click once again on the node and select Create >> Method for NAGLibrary.G01 class >> g01aa (note that the items on the method menu explicitly show the function argument list, including the .NET types; they’re also prefixed with [S] if they’re a static .NET method). This creates a node corresponding to a new instance of the g01aa method in the block diagram (see figure 1).

Figure 1. The g01aa method in the LabVIEW block diagram

The next step in building our demo is to add LabVIEW controls and indicators to handle, respectively, the input and output of our application. These get placed on the front panel and corresponding objects appear on the block diagram, where they’re then connected to other objects to specify the flow of data and control through the program. In our application (see figure 2), we use a Numeric Indicator to display each output quantity from g01aa, and a Numeric Control within an Array shell control to allow the user to enter the data values for which the routine will calculate the statistics (see here for more information on creating and using array controls in LabVIEW).

Figure 2. Front panel of our demo application, with the (input) data array on the left, and the (output) statistics on the right

The block diagram of our application is shown in Figure 3. Along with the objects corresponding to the controls and indicators described above, we’ve also used the Array Size function to determine the number of elements in the Data values array. The output from this is used as input for g01aa's input parameter n (i.e., the length of the array x, which contains the data values), and also as the dimension size input parameter for the Initialize Array function. The other input parameter for this function is element, which specifies the value to which all elements of the array are to be set; we use a DBL Numeric Constant to set this to zero, and the resultant array is used as input for g01aa's input parameter wt, which is the array used to weight the data values in the in the calculation. In fact, in this example, all values are weighted equally; we specify this by using a Numeric Constant to set g01aa's input parameter iwt to zero (see the g01aa documentation for more details).

Figure 3. Block diagram for our demo application. Inputs to g01aa are on the left, and outputs on the right

Although simple in application, this example has hopefully given a flavour of the way in which the NAG library for .NET can be used within LabVIEW. Other methods which are offered by the library address the more complicated problems of local and global optimization of multivariate functions, pseudo and quasi random number generation, wavelet transforms, least-squares and eigenvalue determination, curve and surface fitting, summation of series and interpolation; more details about its functionality are available here.

If you have any questions about this work (or want to know more about using the NAG Fortran Library or the NAG C Library with LabVIEW), then please get in touch. I’d really love to hear from users (or prospective users) of NAG and/or LabVIEW, and want to build new examples that will be useful to researchers in this area.

Finally, I’m grateful to my colleague Jeremy Walton for his help in preparing this post, and his interest in this work.

ISC11 Review

2011-06-27T16:28:00.004+01:00

Several of NAG's staff attended the International Supercomputing Conference (ISC) in Hamburg, Germany last week. We took part in several aspects of the show ("Talk to us in Hamburg"). I also posted a ISC11 preview with comments on I would be looking out for. So, what were my impressions after the show?

First, as ever, ISC is a (the?) great mid-season HPC event - breathing renewed life into the HPC community half-way around the world and half-way around the year from the biggest show of the HPC calandar ("SC" in USA each November).

In my preview, I said I was watching out for three battles:

GPU vs MIC vs Fusion The fight for top voice in manycore/GPU world will be one interesting theme of ISC11. Will this be the year that the GPU/manycore theme really means more than just NVidia and CUDA? AMD has opened the lid on Fusion in recent weeks and has sparked some real interest. Intel's MIC (or Knights) is probably set for some profile at ISC11 now the Knights Ferry program has been running a while. How will NVidia react to no longer being the loudest (only?) noise in GPU/manycore land? Or will NVidia's early momentum carry through?

Review: None of this is definitive, but my gut reaction is that MIC won this battle. GPU lost. Fusion didn't play again. My feeling from talking to attendees was that MIC was second only to the K story, in terms of what people were talking about (and asking NAG - as collaborators in the MIC programme - what we thought). Partly because of the MIC hype, and the K success (performance and power efficient without GPUs), GPUs took a quieter role than recent years. Fusion, disappointingly, once again seemed to have a quiet time in terms of people talking about it (or not). Result? As I thought, manycore is now realistically meaning more than just NVidia/CUDA.

Exascale vs Desktop HPC Both the exascale vision/race/distraction (select according to your preference) and the promise of desktop HPC (personal supercomputing?) have space on the agenda and exhibit floor at ISC11. Which will be the defining scale of the show? Will most attendees be discussing exascale and the research/development challenges to get there? Or will the hopes and constraints of "HPC for the masses" have people talking in the aisles? Will the lone voices trying to link the two extremes be heard? (technology trickle down, market solutions to efficient parallel programming etc.) What about the "missing middle"?

Review: Exascale won this one hands down, I think. Some lone voices still tried to talk about desktop HPC, missing middles, mass usage of HPC and so-on. But exascale got the hype again (not necessarily wrong for one of the year's primary "supercomputing" shows!)

Software vs Hardware The biggie for me. Will this be the year that software really gets as much attention as hardware? Will the challenges and opportunities of major applications renovation get the profile it deserves? Will people just continue to say "and software too". Or will the debate - and actions - start to follow? The themes above might (should) help drive this (porting to GPU, new algorithms for manycore, new paradigms for exascale, etc). Will people trying to understand where to focus their budget get answers? Balance of hardware vs software development vs new skills? Balance of "protect legacy investment" against opportunity of fresh look at applications?

Review: Hardware still got more attention than software. Top500, MIC, etc. Although ease-of-programming for MIC was a common question too. I did miss lots of talks, so perhaps there was more there focusing on applications and software challenges than I caught. But the chat in the corridors was still hardware dominated I thought.

The rest? What have I not listed? National flag waving. I'm not sure I will be watching too closely whether USA, Japan, China, Russia or Europe get the most [systems|petaflops|press releases|whatever]. Nor the issue of cloud vs traditional HPC. I'm not saying those two don't matter. But I am guessing the three topics above will have more impact on the lives of HPC users and technology developers - both next week and for the next year once back at work.

Review: Well, I got those two wrong! Flags were out in force, with Japan (K, Fujitsu, Top500, etc) and France (Bull keynote) waving strongly among others. And clouds were seemingly the question to be asked at every panel! But in a way, I was still right - flags and clouds do matter and will get people talking - but I maintain that manycore, exascale vs desktop, and the desperation of software all matter more.

What did you learn? What stood out for you? Please add your comments and thoughts below ...

Storing Banded Matrices for Speed

2011-06-23T17:29:00.019+01:00

Introduction

A real n-by-n matrix A will, in general, be defined by its n² elements that would need to be stored for computations to be performed using A. If n is very large, this can amount to a lot of storage.

Fortunately, many real world applications use sparse matrices that are comprised mainly of zero-valued elements; the number of defining elements that need to be stored is typically some fixed multiple of n (30n, say). A special case is the banded matrix:

*Illustration 1: A banded matrix A, with k_l sub- and k_u super-diagonals*

This banded matrix has elements defined on its main diagonal (the dashed line), on k_l sub-diagonals (blue region) and on k_usuper-diagonals (green region). The total number of non-zero elements is (k_l+k_u+1)n minus a few corner points. What we would really like to do is solve problems using this matrix, like AX = Y, but storing only the non-zero elements. Various storage formats emerged from various packages, but the most familiar to people solving such problems is that adopted by the linear algebra package LAPACK (and also used extensively in the NAG Libraries). This works by storing the non-zero elements of A in a rectangular, (k_l+k_u+1) by n, matrix B, such that the diagonals of A become rows of B while the columns of A remain as columns of B:

*Illustration 2: The matrix A stored in LAPACK banded storage format as matrix B.*

While this is appealing in terms of being able to reconstruct readily in our mind's eye the original matrix, it is less well suited to solving problems (such as AX = B) as quickly as possible. The problem is that most matrix computations rely, for their speed, on very fast kernel operations (written by the likes of Intel or AMD) working on smallish matrix blocks; this requires matrix blocks to be stored as you would expect for a general matrix, while in B above the block elements have become higgledy-piggledy. What is required is a storage format that keeps some of the block structure of the original matrix A

The Simplest Case - to ease you in

To begin to describe an alternative storage scheme that preserves block structure, it is probably best to begin with the simplest scenario ... a symmetric matrix with bandwidth k and dimension n = (k+1)m (that is, n is an integer multiple of k+1)

*Illustration 3: Symmetric banded matrix with lower half partitioned into triangular blocks.*

Once you partition a lower band (including main diagonal) into triangles as above, the possibilities for improved storage forms become immediately clear. The B triangular blocks can either be slid up or slid across to the right to complete a square (k+1 by k+1) block with a neighbouring A triangular block. To make it easier to convert from LAPACK storage format, I will go for sliding the B blocks up, giving:

*Illustration 4: Storage of symmetric band matrix preserving block structure.*

The zero paranoiacs might want to remove the trailing zeros above by storing the last A₃ triangle in Rectangular Full Packed (RFP) format [see http://www.nag.co.uk/numeric/fl/nagdoc_fl23/xhtml/F07/f07intro.xml#recomm_32a], but this may be a little over-the-top.

What can you do with these blocks?

Those in touch with their linear algebra side will now wonder what can be done with the banded matrix stored in this form. In particular, they will whisper the words Cholesky decomposition. Don't worry, it is just the first of two steps in solving AX = Y, a set of simultaneous equations. Without going into the nitty-gritty, I will just say that the lower triangular banded L formed in this decomposition can be partitioned exactly as above, and it can work with each triangular block in turn, overwriting them as it goes. The second stage of solution, backsubstitution, then works with the triangular blocks of L to return the solution X

How do I put A into this strange storage format?

You couldn't be expected to do this yourself. On the other hand it would be easy to write a utility routine that converted from the LAPACK style format to the above triangular block style; and this can be done in-place since columns are preserved. The mapping between elements of the blocked storage matrix (C, say) and the elements of A is also quite straightforward:

C_qj = A_ij for j = 1,...,n and i = j,...,j+k, where q = mod(i-1,k+1) + 1.

This makes life easy when A_ij is a function of i and j.

The next simplest case -- to lessen your ease

In the simplest case everything fitted just nicely, but for more general n there are going to be some leftovers at the bottom-right of A. A picture tells the story ...

*Illustration 5: Symmetric banded matrix (lower half) with left-overs.*

And the mapping above tells us where to put the left-overs ...

*Illustration 6: Block storage pattern for symmteric banded matrix for general n.*

Next?

This stuff is all very well when we have symmetric matrices, but what about non-symmetric matrices like that in Illustration 1? I'll leave that till the next blog (very soon). Lawrence

Talk to us in Hamburg - ISC11

2011-06-15T10:56:00.003+01:00

NAG will be at the International Supercomputing Conference (ISC) in Hamburg, Germany next week (June 19-23). Come and see us at booth number 830 in the main exhibition hall. We will be demonstrating how our product, services and training can help your HPC applications.

Products on show will include the NAG Library for SMP & Multicore, the NAG Parallel Library, numerical routines for GPUs, the NAG Toolbox for MATLAB®.

We will also be able to tell you more about how NAG’s HPC training courses, impartial technology advice and HPC software engineering services are helping HPC users around the world.

We will have the latest news and code performance success stories from the Computational Science and Engineering Support Service provided by NAG as part of the UK’s HECToR National Supercomputing Service.

NAG is even contributing to the technical program. We are joint hosts (with Intel and RRZE Erlangen) of BoF session 7: “Data‐Parallelism/SIMD for HPC in the Future” on Wednesday. And NAG’s Andrew Jones is one of the four panel members on the conference closing “Analyst Crossfire”, hosted by Addison Snell of InterSect360 on Thursday.

Finally, NAG is expanding our HPC team in the USA and are looking for HPC programmers or CSE experts to be based at a customer site in Houston – so please come and talk to us if you are interested or can recommend anyone.

Now hear this

2011-06-14T11:46:00.004+01:00

Recently I’ve been doing some work producing a web-based seminar (or webinar, to use the fashionable neologism) which demonstrates some aspects of NAG’s optimization routines and their use in the finance industry (specifically for the optimization of financial portfolios). We originally scheduled a live broadcast of the seminar using the GoToWebinar package, which works very well; we’ve also used Webex’s offering, with similar results. Unfortunately, in spite of a lot of preparation and rehearsal, some technical problems (unconnected with the package) arose during the broadcast which rendered the recording of the session useless. Reasoning that better results might be obtained through simplification (our broadcast was a joint effort between two sites, which contributed to the problems), we decided to try just recording the seminar and put it up as a video.

And so it was last week that I found myself talking to my computer in an isolated room whilst deftly manipulating slides and demos. It was just as if I was giving a real live seminar - albeit one that nobody had turned up to (come to think of it, that wasn’t the first time such a thing had happened). I used CamStudio 2.0 to record the session (an earlier attempt using the more up-to-date CamStudio 2.6b failed on this Windows 7 machine for some unknown reason) and Windows Live Movie Maker to do some very minimal editing of the resultant video file. The results can be seen and heard here and here (there's also some explanatory text about the seminar on this page for the curious). I was quite pleased with the way they turned out, although I’d be interested to hear from anyone who knows of a tool that can automatically remove the large number of um’s and err’s from the soundtrack.

.NET examples: C++ and using Object fields

2011-06-02T15:58:00.004+01:00

Recently we had a question from a user of the NAG Library for .NET. He wondered whether it was possible to use methods from objects rather than static methods as callback parameters to the NAG methods. The library was designed such that this should be possible, but we were made aware that there were not many examples of this usage in the documentation. Also the user was coding in C++/CLR and currently we have no examples of this, all examples being C#, F# or VB.Net.

We took the existing C# example for the optimisation routine e04uc and modified it in two ways. Firstly converting to C++/CLR, and then modifying the example so that the objective function being minimised is a method of an object (of locally declared class NagDelegate) and to demonstrate two ways in which data may be passed between successive calls to this callback parameter.

Note that C++/CLR has extended syntax (notably ^ and gcnew in this example) which declare managed data that will be managed and potentially garbage collected by the .NET runtime system.
Once the class and its method are defined, the only thing that needs to be done to pass this method to the e04uc method is to declare the delegate object using the syntax

E04::E04UC_OBJFUN^ objfunE  =
    gcnew E04::E04UC_OBJFUN(nagdelegate, &NagDelegate::objfun);

Unlike the situation in C#, you need to specify both the class and the particular method as two separate parameters when declaring the delegate.

As mentioned above, the method NagDelegate::objfun demonstrates two ways in which data may be accessed, in both cases here the data is a counter reporting how many times the objective function was evaluated by the optimisation algorithm. static_counter is a static variable and counter is a local field of the NagDelegate object. As may be seen by the results when the program is run, both methods work as expected, returning the same value.

  static_counter++;
  counter++;

The C++ example program is shown in full at the end of this post. To run it using Visual Studio, just start a new project and replace the default program by the code below. You must then select the project references and add a reference to the NagLibrary32 or NagLibrary64 .NET assemblies. The program should then run and produce the output shown.

// csharpfromcpp.cpp : main project file.
  
 #include "stdafx.h"
  
 using namespace System;
 using namespace NagLibrary;
  
 static int static_counter;
                           
 public ref class NagDelegate
 {
 public:
     int counter;
     void objfun(int% mode, int n, array <double>^ x, double% objf,
                 array <double>^ objgrd, int nstate)
     {
         //      Routine to evaluate objective function
         //      and its 1st derivatives.
         static_counter++;
         counter++;
         double one = 1.00e0;
         double two = 2.00e0;
         if ((mode == 0) || (mode == 2))
             {
                 objf = x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2];
             }
         //
         if ((mode == 1) || (mode == 2))
             {
                 objgrd[0] = x[3] * (two * x[0] + x[1] + x[2]);
                 objgrd[1] = x[0] * x[3];
                 objgrd[2] = x[0] * x[3] + one;
                 objgrd[3] = x[0] * (x[0] + x[1] + x[2]);
             }
         //
     }
 };
  
 void confun( int% mode, int ncnln, int n, array< int >^ needc,
               array< double >^ x, array <double>^ c,
                     array <double,2>^ cjac, int nstate)
     {
       //      Routine to evaluate the nonlinear constraints
       //       and their 1st derivatives.
       double zero = 0.00e0;
       double two = 2.00e0;
       int i,  j;
       if (nstate == 1)
       {
         //  First call to confun.  Set all Jacobian elements to zero.
         //  Note that this will only work when 'Derivative Level = 3'
         //  (the default; see Section 11.2).
         for (j = 1; j <= n; j++)
         {
           for (i = 1; i <= ncnln; i++)
           {
             cjac[i - 1, j - 1] = zero;
           }
         }
       }
       //
       if (needc[0] > 0)
       {
         if ((mode == 0) || (mode == 2))
         {
           c[0] = (x[0]) * (x[0]) +
                          (x[1]) * (x[1]) + (x[2]) * (x[2]) + (x[3]) * (x[3]);
         }
         if ((mode == 1) || (mode == 2))
         {
           cjac[0, 0] = two * x[0];
           cjac[0, 1] = two * x[1];
           cjac[0, 2] = two * x[2];
           cjac[0, 3] = two * x[3];
         }
       }
       //
       if (needc[1] > 0)
       {
         if ((mode == 0) || (mode == 2))
         {
           c[1] = x[0] * x[1] * x[2] * x[3];
         }
         if ((mode == 1) || (mode == 2))
         {
           cjac[1, 0] = x[1] * x[2] * x[3];
           cjac[1, 1] = x[0] * x[2] * x[3];
           cjac[1, 2] = x[0] * x[1] * x[3];
           cjac[1, 3] = x[0] * x[1] * x[2];
         }
       }
       //
     }
  
  
 int main(array<System::String ^> ^args)
 {
         double objf = 0.0;
         int i,  ifail,  iter,  j,  n,  nclin,  ncnln;
         Console::WriteLine(L"e04uc Example Program Results");
         //      Skip heading in data file
         n = 4;
         nclin = 1;
         ncnln = 2;
         array< double, 2 >^ a = gcnew array< double, 2 >( nclin, n);
         a[0, 0]=1;
         a[0, 1]=1;
         a[0, 2]=1;
         a[0, 3]=1;
  
         array< double >^ bl = gcnew array< double >(n + nclin +ncnln);
         bl[0] = 1.0;
         bl[1] = 1.0;
         bl[2] = 1.0;
         bl[3] = 1.0;
         bl[4] = -1.0e+25;
         bl[5] = -1.0e+25;
         bl[6] = 25.0;
  
         array< double >^ bu = gcnew array< double >(n + nclin +ncnln);
         bu[0] = 5.0;
         bu[1] = 5.0;
         bu[2] = 5.0;
         bu[3] = 5.0;
         bu[4] = 20.0;
         bu[5] = 40.0;
         bu[6] = 1.0e25;
  
  
         array< double >^ c = gcnew array< double >(ncnln);
         array< double, 2 >^ cjac = gcnew array< double, 2 >(ncnln, n);
         array< double >^ clamda = gcnew array< double >(n + nclin +ncnln);
         array< double >^ objgrd = gcnew array< double >(n);
         array< double, 2 >^ r = gcnew array< double, 2 >(n, n);
         array< double >^ x = gcnew array< double >(n);
         x[0] = 1.0;        
         x[1] = 5.0;        
         x[2] = 5.0;        
         x[3] = 1.0;        
  
         array< int >^ istate = gcnew array< int >(n + nclin +ncnln);
  
              E04::e04ucOptions^ options = gcnew E04::e04ucOptions() ;
              E04::E04UC_CONFUN^ confunE = gcnew  E04::E04UC_CONFUN(confun);
  

              NagDelegate^ nagdelegate = 
                        gcnew NagDelegate();
              E04::E04UC_OBJFUN^ objfunE  =
                        gcnew  E04::E04UC_OBJFUN(nagdelegate,  &NagDelegate::objfun);
              E04::e04uc(n, nclin, ncnln, a, bl, bu, confunE, objfunE, iter, istate, c,
                         cjac, clamda, objf, objgrd, r, x, options, ifail);
  
              Console::WriteLine("The objective value on return is {0}. ",
                      objf);
              Console::WriteLine("Objfun has been  called {0} times, using class member",
                      nagdelegate->counter);
              Console::WriteLine("Objfun has been  called {0} times, using a static variable",
                      static_counter);
  
  
  
     return 0;
 }

When run the following output is produced:

e04uc Example Program Results
The objective value on return is 17.0140172891347.
Objfun has been  called 7 times, using class member
Objfun has been  called 7 times, using a static variable

Programming gnuplot

2011-05-23T15:03:00.020+01:00

One of the side effects of having the word 'visualisation' in your job title is that people expect you to know something about the subject. Accordingly, every now and again I get contacted by one of my colleagues who needs to draw a graph, or reformat an image, or incorporate a visualisation into a document, and is - usually - too busy with other (probably more important) activities to spend a lot of time figuring out how to make a graphics package behave itself and do the right thing.

I received the latest of these queries last week, from a colleague who wanted to reproduce Figure 1 (below) as part of the documentation for the NAG Parallel Library. It's a schematic illustration of the way in which, for certain routines in the Library, elements of a matrix are distributed between compute processors. More specifically, this example shows the distribution of a two-dimensional matrix consisting of twelve rows and columns (each numbered from 1 to 12) across a two by three grid of processors.

Figure 1. Data distribution across a grid of processors

The request was to reproduce the figure using gnuplot, a plotting program that can be used to create a wide variety of two- and three-dimensional plots of data, functions and approximations to data. It runs on a variety of platforms, and is freely distributed. Gnuplot's interface is command-driven - i.e. the user either enters commands one at a time at its interactive prompt, or puts the commands into a script, which the program then loads and runs. Several example scripts (and their associated output) are available at the gnuplot site - for example, these illustrate the range of the program's functionality for plotting analytical functions. We use the program a lot for the display of the results from the example programs which form an important part of NAG Library documentation; thus, Figure 2 is a plot of the output from the c06ab routine, which is used in the numerical summation of series.

Figure 2. Results from the c06ba example program

Whilst producing this sort of picture in gnuplot is comparatively straightforward, and can be done with a handful of commands, the diagram in Figure 1 presents a few more challenges. However, the elementary components aren't hard to find: thus, gnuplot has a set object rectangle command which creates a rectangle, given the coordinates of its diagonal corners, a set arrow command which creates an arrow (or a line) between two endpoints, and a set label command which adds text at a specified location. So I started writing the script in this fashion:

# draw box.
set object 1 rect from 0,0 to 6,8.2 

# add label.
set label "{1,0}" at 3,4.1 center

# draw y boundaries.
set arrow 2 from 1.5,0 to 1.5,8.2 nohead
set arrow 3 from 3.0,0 to 3.0,8.2 nohead
set arrow 4 from 4.5,0 to 4.5,8.2 nohead
...

but I hadn't got very far before realising that hardwiring all of the coordinate values wasn't a good idea if - say - the size of the boxes, or the spacing between them, needed changing. And I was getting tired of keeping track of the running number of objects, and working out in my head the locations of the boundaries. Fortunately, I recalled that a computer is good at this sort of thing (the clue is in the name, after all). I also discovered that gnuplot allows user-defined variables to be declared, and supports mathematical operations on them. This made writing the script a lot easier, but there's a large degree of repetition in the code which - although straightforward to produce by cutting and pasting - looks inelegant and verbose. In a more fully-featured programming environment, this drawback would be overcome by using control statements such as loops or branches, or by using subprograms. However, these features aren't supported by gnuplot (or if they are, I couldn't find them - and the call of more important activities was becoming too strident to ignore), so we're left with a lengthy script, which is reproduced below in all its resplendence for the curious.

Buried within it (almost at the end) is one technical point that might be of wider interest: because gnuplot is, when all's said and done, a plotting package, each script must contain at least one plot command, otherwise nothing will appear. More specifically in our example, we've used the set command to create items to be displayed by specifying values for its object, arrow and label options, but these won't be drawn on the graph until the plot command is issued. This has a compulsory argument, but having specified everything to be displayed, we don't have anything else left to be explicitly plotted. Other users who have used the command in similar circumstances (for example, here and here) have used the expression 1/0 for this argument; gnuplot evaluates this to NaN, which is interpreted as an undefined point, and is quietly ignored by the plot command when it's displaying the other objects in the figure.

# gnuplot file for figure (data distribution from the 
# processors' point of view) from 'essential introduction 
# to nag parallel library'

# NB: delete '#EPS' or '#PNG' from the following set commands to 
# generate individual image files
#EPSset term postscript eps enhanced monochrome font "Times" 
#EPSset output "fdessint2_fig.ps"
#PNGset term png enhanced medium font "TIMES"
#PNGset output "fdessint2_fig.png"

# set values for total x and y dimensions.
xLength = 25
yLength = 25

# set values for overall offset at left and top, to make room 
# for labels.
xLabelOffset = 1.0
yLabelOffset = 1.0

# set values for x and y box sizes.
xBox = 6.0
yBox = 8.2

# set values for number of cells within each box.
xNumCell = 4
yNumCell = 6

# calcuate offset between each box.  There are 3 boxes along x, 
# and 2 along y.
xBoxOffset = (xLength - xLabelOffset - 3*xBox)/2.0
yBoxOffset = yLength - yLabelOffset - 2*yBox

# calculate spacing between each cell.
xCellSpace = xBox/xNumCell
yCellSpace = yBox/yNumCell

# Specify the extent of the plotting space.
set xrange [0:xLength]
set yrange [0:yLength] 

# specify the aspect ratio of the plot.
set size ratio 0.8

# turn off the key, the axis tics and the border around the plot.
unset key
unset xtics
unset ytics
unset border

# set the fill style (for the boxes).
set style fill empty border rgb "black"

# initialize counter for objects.
nobj = 0

############
# initialize location of box (bottom row). 
rx = xLabelOffset
ry = 0

# draw box.
nobj = nobj+1
set object nobj rect from rx,ry to rx+xBox,ry+yBox 

# add label.
set label "{1,0}" at rx+0.5*xBox,ry+0.5*yBox center

# initialize coordinates for cell boundaries. 
ax = rx
ay = ry

# draw y boundaries.
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead

# draw x boundaries.
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead

# draw next box.
rx = rx+xBox+xBoxOffset 
nobj = nobj+1
set object nobj rect from rx,ry to rx+xBox,ry+yBox 

# add label.
set label "{1,1}" at rx+0.5*xBox,ry+0.5*yBox center

# initialize coordinates for cell boundaries. 
ax = rx
ay = ry

# draw y boundaries.
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead

# draw x boundaries.
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead

# draw next box.
rx = rx+xBox+xBoxOffset 
nobj = nobj+1
set object nobj rect from rx,ry to rx+xBox,ry+yBox 

# add label.
set label "{1,2}" at rx+0.5*xBox,ry+0.5*yBox center

# initialize coordinates for cell boundaries. 
ax = rx
ay = ry

# draw y boundaries.
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead

# draw x boundaries.
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead

############
# initialize location of box (top row). 
rx = xLabelOffset
ry = yBox + yBoxOffset

# draw next box.
nobj = nobj+1
set object nobj rect from rx,ry to rx+xBox,ry+yBox 

# add label.
set label "{0,0}" at rx+0.5*xBox,ry+0.5*yBox center

# initialize coordinates for cell boundaries. 
ax = rx
ay = ry

# draw y boundaries.
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead

# draw x boundaries.
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead

# draw next box.
rx = rx+xBox+xBoxOffset 
nobj = nobj+1
set object nobj rect from rx,ry to rx+xBox,ry+yBox 

# add label.
set label "{0,1}" at rx+0.5*xBox,ry+0.5*yBox center

# initialize coordinates for cell boundaries. 
ax = rx
ay = ry

# draw y boundaries.
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead

# draw x boundaries.
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead

# draw next box.
rx = rx+xBox+xBoxOffset 
nobj = nobj+1
set object nobj rect from rx,ry to rx+xBox,ry+yBox 

# add label.
set label "{0,2}" at rx+0.5*xBox,ry+0.5*yBox center

# initialize coordinates for cell boundaries. 
ax = rx
ay = ry

# draw y boundaries.
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead
ax = ax+xCellSpace
nobj = nobj+1; set arrow nobj from ax,ry to ax,ry+yBox nohead

# draw x boundaries.
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead
ay = ay+yCellSpace
nobj = nobj+1; set arrow nobj from rx,ay to rx+xBox,ay nohead

############
# Now draw the arrows between boxes.  
# First, in y direction (between rows).
ax = xLabelOffset + 0.5*xBox

# shave off a bit of the arrow so it's not touching the box 
# (looks nicer).
eps = 0.1*yBoxOffset

nobj = nobj+1 
set arrow nobj from ax,yBox+eps to ax,yBox+yBoxOffset-eps heads
ax = ax + xBox + xBoxOffset 
nobj = nobj+1
set arrow nobj from ax,yBox+eps to ax,yBox+yBoxOffset-eps heads
ax = ax + xBox + xBoxOffset 
nobj = nobj+1
set arrow nobj from ax,yBox+eps to ax,yBox+yBoxOffset-eps heads

# Second, in x direction (between columns - bottom row).
ax = xLabelOffset + xBox
ay = 0.5*yBox
eps = 0.1*xBoxOffset 

nobj = nobj+1
set arrow nobj from ax+eps,ay to ax+xBoxOffset-eps,ay heads
ax = ax + xBox + xBoxOffset 
nobj = nobj+1
set arrow nobj from ax+eps,ay to ax+xBoxOffset-eps,ay heads

# Third, in x direction (between columns - top row).
ax = xLabelOffset + xBox
ay = ay + yBox + yBoxOffset
eps = 0.1*xBoxOffset 

nobj = nobj+1
set arrow nobj from ax+eps,ay to ax+xBoxOffset-eps,ay heads
ax = ax + xBox + xBoxOffset 
nobj = nobj+1
set arrow nobj from ax+eps,ay to ax+xBoxOffset-eps,ay heads

############
# Add labels.  First, left of bottom row.
lx = 0.5*xLabelOffset
ly = 0.5*yCellSpace
set label "12" at lx,ly center  
ly = ly+yCellSpace
set label "10" at lx,ly center
ly = ly+yCellSpace
set label "8" at lx,ly center
ly = ly+yCellSpace
set label "6" at lx,ly center
ly = ly+yCellSpace
set label "4" at lx,ly center
ly = ly+yCellSpace
set label "2" at lx,ly center

# Second, left of top row.
lx = 0.5*xLabelOffset
ly = yBox + yBoxOffset + 0.5*yCellSpace
set label "11" at lx,ly center
ly = ly+yCellSpace
set label "9" at lx,ly center
ly = ly+yCellSpace
set label "7" at lx,ly center
ly = ly+yCellSpace
set label "5" at lx,ly center
ly = ly+yCellSpace
set label "3" at lx,ly center
ly = ly+yCellSpace
set label "1" at lx,ly center

# Third, top of the three columns.
rx = xLabelOffset
ly = yLength - yLabelOffset + 0.5*yCellSpace
lx = rx + 0.5*xCellSpace
set label "1" at lx,ly center
lx = lx+xCellSpace
set label "4" at lx,ly center
lx = lx+xCellSpace
set label "7" at lx,ly center
lx = lx+xCellSpace
set label "10" at lx,ly center

rx = rx+xBox+xBoxOffset 
lx = rx + 0.5*xCellSpace
set label "2" at lx,ly center
lx = lx+xCellSpace
set label "5" at lx,ly center
lx = lx+xCellSpace
set label "8" at lx,ly center
lx = lx+xCellSpace
set label "11" at lx,ly center

rx = rx+xBox+xBoxOffset 
lx = rx + 0.5*xCellSpace
set label "3" at lx,ly center
lx = lx+xCellSpace
set label "6" at lx,ly center
lx = lx+xCellSpace
set label "9" at lx,ly center
lx = lx+xCellSpace
set label "12" at lx,ly center

############
# call plot (with a dummy argument) to display the objects 
# we've defined.
plot 1/0

pause -1

Good Business Analytics Conference

2011-05-04T10:28:00.000+01:00

I have not yet written a note about the good experience I had at the 2011 INFORMS Conference on Business Analytics and Operations Research in Chicago, where NAG had a booth with demonstrations of our library routines. ....I guess 'now' is as good a time as any to do so. The conference was very well attended, up around 30% over last year by many accounts, and included a wide variety of analytical types; academics and practitioners with many different specialities. It was a good opportunity to understand how the best analytical techniques are applied in different environments. During the three days my colleagues and I managed to talk with lots of people. Some were users of NAG routines and some who did not know us at all. It was a pleasure to meet all of you. My note is just to say thanks to fellow attendees, as well as the organizers, for a stimulating 3 days – hope to speak with many of you again soon

Aging Fortran Codes - A Crisis in Slow Motion

2011-04-11T16:38:00.000+01:00

Fortran isn't dead but some of the critical scientific codes written in it will be on life-support soon if research sponsoring organizations don't react soon. It was just a little over four years years ago that we read about the death of John Backus, the leader of the IBM that invented the Fortran language in 1954. It was 54 years ago that the first Fortran compiler was released commercially.

In that span of time millions of lines of Fortran code have been written in Fortran IV, Fortran 66 and 77 as well as variants such as Watfor. In listening to managers of codes at DOE, NASA universities and private companies (including independent software vendors) over the past several years, I am constantly amazed at how much code is in the early versions of Fortran.

People readily acknowledge the amount of code and the concern about the aging "expertise base" (i.e., people) who are charged with maintaining and improving the code. They recognize that much of this code is "mission-critical" to research in physics, chemistry, aerospace and alternative energy and worry about the long certification process for updating modernized versions of the code. And there the "disconnect" starts.

When I ask how much resource they are putting into modernizing the codes to preserve the code and help it cope in a multi/many-core world the answer I often get is, effectively, "that's not one of our budget priorities right now".

This quite unsatisfying situation suggests two responses. For us at NAG who get called upon to port, modernize and scale codes as part of our work with HECToR users and other customers it suggests that we need to continue to find more cost-effective tools and processes to convert, modernize and scale codes. The bigger challenge is for organizations such as those above and other research sponsoring organizations such as the NSF in the US and various research councils in the UK to not only recognize the issue but to find the resources to support the effort before we lose the base of talent and institutional "memory" needed to do the job.

Turning off logic

2011-04-04T14:10:00.004+01:00

Everybody knows about the two halves of our brain. One half deals with logical thought, facts etc. while the other deals with feelings, imagination and other things. They are linked at the base of the brain by the corpus collosum which is the sole place that they communicate with each other. The two halves work together as they take in information from the world through all of your senses, painting a "complete" picture of what is going on around you. For example the left side of your brain deals with the names of objects, whereas the right side deals with the objects function, so that when presented with a computer you are able to name it and describe what it does.

Artists and programmers are probably at two extremes of the spectrum when it comes down to which parts of our brain we use for our work. Artists almost certainly use much of their right brain, whereas I imagine that programmers use a lot of their left brains.

Whatever parts of our brain we use get exercised and will be become more adept at what they do. So for programmers their logic should improve; for artists their artistic skills should improve also. But what happens to those parts of the brain that we don't use? Are they getting worse? Better? Staying the same? What would life be like if we made more of an effort to spread the work load over our whole brain?

I don't know what you do to relax, but usually I find that doing something completely different to my every day activities leads to me feeling very relaxed, even if that something involves physical activity. I guess purposefully not using the part of your brain that gets the most daily use, means that it gets to relax while another function of the brain is being used.

As a programmer I find that walking in the country side, enjoying the beauty of nature, gardening, creating something artistic are incredibly relaxing for me. I imagine conversely that if I was doing gardening or art for my job I would really enjoy programming as a relaxing thing to do.
All of this is just speculatory as I only have my personal experience to go on and haven't researched into any scientific evidence.

So my question is then, do you find the same thing? Do you need to relax by doing something completely different, whether by gardening, or doing something arty? If you do need to, but you find that you don't have the time, do you experience any negative consequences from not relaxing?
I know that I tend to start treating my family like programs, and expecting everything to be logical if I don't relax. This feeling of logicalness then becomes slightly detrimental when dealing with the emotions of a stressed wife or a frustrated toddler.

So do you relax from a days programming? And if so how?

Sticking to coding guidelines

2011-03-21T10:55:00.002Z

Here at NAG we have a lot of rules governing how we write our code. These include things like which constructs we're allowed to use (not all compilers support complete specifications of languages such as Fortran 2003), how we allocate and deallocate memory (we insist on explicit deallocation in the same program scope as the original allocation to avoid memory leaks), how we name and declare variables etc. For many of these we have tools which can verify that the rules are being adhered to, some of which are based on our own compiler. For the rest, however, we rely on an internal peer review process where a colleague checks our work independently to ensure that it conforms to our standards.

Recently a case arose where a relatively new colleague was writing a rather complicated piece of code and he came up against a problem. Being a clever chap he came up with an elegant solution which worked perfectly in the context that he was working in and required only a minimal change to the user interface of his routine. Strictly speaking that change was a subtle deviation from our standards since it involved some arrays being declared with explicit lengths rather than being assumed size. What he didn't realise was that, when this code came to be incorporated into some of our other products, his solution would cause problems with the way the parameters were being used there, requiring a lot of effort to work around.

Luckily he was paired with one of our most eagle-eyed peer reviewers who queried the fact that these arrays did not conform to our normal practice. This in turn led to the original problem and the issues with applying his solution in our other products being aired. In the short term a fix conforming with our standards is being implemented which will work correctly in all our products, and in the longer term we will adopt an improved mechanism, better to address his original problem.

One lesson from this is that, particularly in a multi-developer environment like ours, its important to stick to agreed standards since its not always clear what the consequences of breaking them will be. But we all make mistakes and it would be nice if we had tools which could validate stylistic programming rules alongside the more mechanistic ones. It might be possible to do this with formal specification languages but that seems like using a sledgehammer to crack a nut. We already have tools which check the consistency of our documentation and our software, e.g. to ensure that the types and dimensions of parameters match. Maybe we should develop this idea further to include more complicated checks.

Performance and Results

2011-03-18T16:38:00.002Z

What's in a catch phrase?

As you will hopefully know, NAG's strapline is "Results Matter. Trust NAG".

What matters to you, our customers, is results. Correct results that you can rely on. Our strapline invites you to trust NAG - our people and our software products - to deliver that for you.

When I joined NAG to help develop the High Performance Computing (HPC) services and consulting business, one of the early discussions raised the possibility of using a new version of this strapline for our HPC business, reflecting the performance emphasis of the increased HPC activity. Probably the best suggestion was "Performance Matters. Trust NAG." Close second was "Productivity Matters. Trust NAG."

Our (then newly proactive but not new) HPC consulting business offers to help customers extract better performance from parallel computing technologies - whether multicore workstations, clusters with hundreds of cores, or multi-thousand node supercomputers. In many cases, application software must be evolved, optimized, or occasionally substantially re-architected with new algorithms, to effectively exploit the processing power of multicore based workstations and supercomputers. We do this through collaborative development of your applications to introduce parallel processing, or to improve scalability, or to develop new algorithms suited for the next generation of computing technology. We also offer training courses in the relevant skills and technologies (e.g. OpenMP, MPI, Fortran). And some insights on HPC via the NAG blog.

However, we also recognised that as well as our comprehensive HPC skills base, NAG brings an added value to the market - experience of developing portable, robust, production quality software that users can rely on for accurate results.

HPC and supercomputing often focuses on getting the fastest performance. But that isn't much use unless you can trust those answers too. The converse is true also - answers must be delivered in a useful timescale - performance matters. I have written about this performance/robustness balance several times - on the NAG blog and in my ZDNet UK columns. Good enough performance must always be balanced with accuracy and reliability.

And so, we concluded that keeping the same strapline, "Results Matter. Trust NAG.", was just right for the HPC aspects of NAG's business too.

Because, in the end, Performance is part of Results - in fact, while the Results out of the code are obviously important, it is the Results out of your business or research that really Matters - and scalable application Performance could be the factor that drives your research, design, or even economic Results.

The Addictive Allure of Supercomputing

2011-03-17T10:26:00.002Z

The European Medical Device Technology (EMDT) magazine interviewed me recently. InsideHPC also has pointed to the interview here.

The interview discusses false hopes of users: "Computers will always get faster – I just have to wait for the next processor and my application will run faster."

We still see this so often - managers, researchers, programmers even - all waiting for the silver bullet that will make multicore processors run their application faster with no extra effort from them. There is nothing now or coming soon that will do that except for a few special cases. Getting performance from multicore processors means evolving your code for parallel processing. Tools and parallelized library plugins can help - but in many cases they won't be a substitute for re-writing key parts of the code using multithreading or similar techniques.

Of course, NAG can help ...

NAG out and about

2011-03-10T14:00:00.016Z

The NAG website has a section called "Meet our experts - NAG out and about", which gives a list of upcoming events worldwide that NAG experts will be attending or presenting at.

The page also notes: "We regularly organise and participate in conferences, seminars and training days with our customers and partners. If you would like to talk to us about hosting a NAG seminar at your organisation or any training requirements you might have email us at sales@nag.co.uk".

In my own focus of high performance computing (HPC), I have previously written (for ZDNet UK) about some key supercomputing events. For those of you interested in meeting up with HPC experts (especially from NAG!), I have set up a survey of HPC events - please let us know which events you plan to attend in 2011 - and see which events other readers of The NAG Blog are attending.

Calling NAG routines from IDL

2011-02-15T10:25:00.026Z

The collection of numerical algorithms contained in the NAG Library has been a valuable and widely-used resource for many years, and it's something that we're continually updating and maintaining lest the expression "has been" becomes an uncomfortably relevant adjective. To this end, we're - for example - developing and testing new routines for the next Mark of the Library, have just released a new version of the Library for the .NET environment, and are continuing to try and make our algorithms accessible from as many languages and packages as possible. Besides .NET, these include MATLAB^®, Excel, Maple, Java, Python and many others (some of which are listed on this page). In spite of our best efforts in this area however, there will always be environments that we've yet to cover fully.

Fortunately, our users - besides possessing a commendably high level of discernment in their selection of numerical libraries - are talented and resourceful individuals. For example, we've recently been in contact with a group at Culham Centre for Fusion Energy, which is the UK's national fusion research laboratory. They're users of the popular IDL programming language and visualization system, and were keen to enhance IDL's functionality by calling NAG Library routines from within that environment. They found that they were able to do this by writing a C wrapper around each Fortran routine call, and then passing the definition of the wrapper routine to IDL via its Dynamically Loadable Module (DLM) packaging method. We highlight their work - which they've been kind enough to share with us - using the E02BEF routine as an example; this calculates the cubic spline approximation to a set of data points (and has previously formed the basis of one of the demos of the NAG Toolbox for MATLAB^®)).

We first give the wrapper function, which checks its input parameters and then calls the NAG routine. The function prototypes for all the routines in the NAG Fortran Library are contained in the header file nagmk22.h; this is part of the Fortran Library distribution, and makes it easier to call any of the routines from within C or C++. We note in passing that the wrapper is only required because we're using the NAG Fortran Library. IDL is written in C, so DLMs must also be written in C; if the NAG C Library was being used instead, this code would be simpler.

//---------------------------------------------------------------
void IDL_CDECL idl_e02bef(int argc, IDL_VPTR argv[], char *argk){
// IDL wrapper function for NAG routine E02BEF.
//
// Arguments: 
//
//  0 start  char*   input  "C" or "W": cold or warm start
//  1 NP     int     input  number of data points
//  2 x      double* input  independent (abscissa) values      
//  3 y      double* input  dependent (ordinate) values
//  4 w      double* input  weights
//  5 s      double  input  smoothing factor 
//  6 nest   int     input  over-estimate of the no. of knots
//  7 n      int     input  prior no. of knots for a warm start
//                   output total number of spline knots      
//  8 lambda double* input  prior knot positions
//                   output knot positions
//  9 c      double* output spline coefficients
// 10 fp     double  output sum of square weighted residuals 
// 11 wrk    double* input  prior work array for a warm start
//                   output work array
// 12 lwrk   int     input  size of work array: 4*NP+16*nest+41
// 13 iwrk   int*    input  prior integer work array (warm start)
//                   output integer work array
// 14 ifail  int     input  must be set to 0, -1 or 1
//                   output 0 unless an error occurred
//
//      call E02BEF(start, NP, x, y, w, s, nest, n, 
//                  lamda, c, fp, wrk, lwrk, iwrk, ifail) 
//
//---------------------------------------------------------------
   char start[2];
   int i, NP, nest, lwrk, n, *iwrk, ifail, length_start;
   double s, fp, *x, *y, *w, *lambda, *c, *wrk;
   
// Import all arguments, and check their values.
   
// First, the input parameters only.      
   IDL_ENSURE_STRING(argv[0]);      
   IDL_ENSURE_SCALAR(argv[0]);
   strncpy(start, IDL_STRING_STR(&argv[0]->value.str), 1);
   start[1] = '\0';
   if(strlen(IDL_STRING_STR(&argv[0]->value.str)) != 1) 
      fprintf(stderr,"Start String Length Incorrect\n");

   IDL_ENSURE_SCALAR(argv[1]);
   NP = (int)IDL_LongScalar(argv[1]); 
   
   IDL_ENSURE_ARRAY(argv[2]);
   x = (double *)argv[2]->value.arr->data;   
   if(argv[2]->value.arr->n_elts < NP) 
      fprintf(stderr,"X Array Too Small\n");
      
   IDL_ENSURE_ARRAY(argv[3]);
   y = (double *)argv[3]->value.arr->data;   
   if(argv[3]->value.arr->n_elts < NP) 
      fprintf(stderr,"Y Array Too Small\n");  
    
   IDL_ENSURE_ARRAY(argv[4]);
   w = (double *)argv[4]->value.arr->data;   
   if(argv[4]->value.arr->n_elts < NP) 
      fprintf(stderr,"W Array Too Small\n");  
   
   IDL_ENSURE_SCALAR(argv[5]);
   s = (double)IDL_DoubleScalar(argv[5]); 
   
   IDL_ENSURE_SCALAR(argv[6]); 
   nest = (int)IDL_LongScalar(argv[6]); 
   
   IDL_ENSURE_SCALAR(argv[12]); 
   lwrk = (int)IDL_LongScalar(argv[12]); 
   
// Second, the input & output parameters.  
   IDL_ENSURE_SCALAR(argv[7]);
   IDL_EXCLUDE_EXPR(argv[7]);      
   n = (int)IDL_LongScalar(argv[7]);
        
   IDL_ENSURE_ARRAY(argv[8]);
   IDL_EXCLUDE_EXPR(argv[8]);
   lambda = (double *)argv[8]->value.arr->data;   
   if(argv[8]->value.arr->n_elts < nest) 
      fprintf(stderr,"Lambda Array Too Small\n");
           
   IDL_ENSURE_ARRAY(argv[9]);
   IDL_EXCLUDE_EXPR(argv[9]);
   c = (double *)argv[9]->value.arr->data;
   if(argv[9]->value.arr->n_elts < nest) 
      fprintf(stderr,"C Array Too Small\n");
   
   IDL_ENSURE_SCALAR(argv[10]);
   IDL_EXCLUDE_EXPR(argv[10]); 
   fp = (double)IDL_DoubleScalar(argv[10]);
   
   IDL_ENSURE_ARRAY(argv[11]);
   IDL_EXCLUDE_EXPR(argv[11]);
   wrk = (double *)argv[11]->value.arr->data; 
   if(argv[11]->value.arr->n_elts < lwrk) 
      fprintf(stderr,"Work Array Too Small\n");
   
   IDL_ENSURE_ARRAY(argv[13]);
   IDL_EXCLUDE_EXPR(argv[13]);
   iwrk = (int *)argv[13]->value.arr->data;
   if(argv[13]->value.arr->n_elts < lwrk) 
      fprintf(stderr,"IWork Array Too Small\n");
   
   IDL_ENSURE_SCALAR(argv[14]);
   IDL_EXCLUDE_EXPR(argv[14]);      
   ifail = (int)IDL_LongScalar(argv[14]);        

//---------------------------------------------------------------
// Call the NAG routine.  It was written in Fortran; here, we 
// assume we're running under Linux or Unix, where most 
// compilers convert the name to lower case and prepend an 
// underscore.  The function prototype is defined in nagmk22.h 
// (above).
//
// Because one of the parameters is a character string (and 
// we're calling Fortran from C), we also have to pass its 
// length.  Under this calling convention, this goes at the end 
// of the argument list.
   length_start = strlen(start); 
   e02bef_(start, &NP, x, y, w, &s, &nest, &n, lambda, c, &fp, 
      wrk, &lwrk, iwrk, &ifail, length_start); 

// Check for unsuccessful call of NAG routine (see documentation)
// and output warning if necessary.
   if(ifail != 0) 
      fprintf(stderr,
      "Warning: e02bef returned with ifail = %d\n", ifail);
                            
//---------------------------------------------------------------
// Return the scalar data to IDL (we don't have to explicitly 
// return the arrays, as they're automatically passed by 
// pointer/reference).
   IDL_StoreScalar(argv[7],  IDL_TYP_LONG,   (IDL_ALLTYPES*)&n);
   IDL_StoreScalar(argv[10], IDL_TYP_DOUBLE, (IDL_ALLTYPES*)&fp);
   IDL_StoreScalar(argv[14], IDL_TYP_LONG, 
      (IDL_ALLTYPES*)&ifail);
}

Each NAG routine that we want to use would require a separate wrapper, constructed along similar lines to those illustrated in this example.

Here's the rest of the code (which should precede the specification of the wrapper function in our file - which could be called, say, idlnag.c). Along with the NAG function prototypes, we also include idl_export.h, which defines the IDL API. Next, after defining a useful macro and the prototype for the wrapper function, we set up the DLM by first defining a structure idl_nagprocedures that will be used to pass information about our routine to IDL, and then using that in the specification of a function called IDL_Load. Every DLM must have this function, which is called when the DLM is first loaded; it calls IDL_SysRtnAdd to load the routine definition, and then defines the exit handler (which does nothing in our case).

#include <stdio.h>
#include <string.h>
#include "idl_export.h" // IDL API Header
#include "nagmk22.h"    // NAG C/Fortran Prototypes

// Useful macro.  
#define ARRLEN(arr) (sizeof(arr)/sizeof(arr[0]))

//---------------------------------------------------------------
// Procedure prototypes - one for each routine in this module.
extern void IDL_CDECL idl_e02bef(int argc, 
   IDL_VPTR argv[], char *argk);

//---------------------------------------------------------------
// Structure used to pass information about our routine to IDL.  
// Other routines would require similar entries in this table.
static IDL_SYSFUN_DEF2 idlnag_procedures[] = {
   {(IDL_FUN_RET) idl_e02bef, // name of routine in here.
   "IDL_E02BEF",              // name of the IDL procedure.
   15,            // Min no. of arguments for routine.
   15,            // Max "   "      "      "     " 
   0,             // Further options which we don't use.
   0              // Reserved for IDL use - always set to 0.
   },
};

//---------------------------------------------------------------
// This function is called when IDL is shut down.
void idlnag_exit_handler(void) {
   // Nothing to do.
}

//---------------------------------------------------------------
int IDL_Load(void) {
// This function is called when IDL loads this module.  

// Call IDL_SysRtnAdd to register the routine contained in this 
// module.  Note that this routine is what IDL calls a 
// procedure - i.e. it doesn't return a value (cf. what IDL 
// calls a function, which does).  This is specified by setting 
// the second parameter in the call to FALSE (it's TRUE for 
// functions). 
   if (!IDL_SysRtnAdd(idlnag_procedures, FALSE, 
      ARRLEN(idlnag_procedures))) return IDL_FALSE;

// Register the exit handler.
   IDL_ExitRegister(idlnag_exit_handler);

   return(IDL_TRUE);
}

idlnag.c must then be compiled and linked into a shared library (called idlnag.so on Unix) which should be placed in a location on the directory path specified by the IDL_DLM_PATH environment variable. The other file that must be in that location is the so-called DLM module description file. In our example, its contents would be:

MODULE idlnag                    # Name of the DLM
DESCRIPTION IDL Interface to NAG # Optional one-line description
VERSION 1.0                      # Optional version number
PROCEDURE IDL_E02BEF 15 15 # Routine name, min, max no. of args

Here, the name of the DLM must be the same as the name of the shared library we've created, whilst the routine name and the minimum and maximum number of arguments must be the same as that specified in the IDL_SYSFUN_DEF2 structure - i.e. idlnag_procedures in our case.

Having defined the new function in this fashion, it will be loaded the next time that IDL is started, and it can be invoked within IDL as

IDL_E02BEF, start, NP, x, y, w, s, nest, n, lambda, c, fp, $
            wrk, lwrk, iwrk, ifail

where start is an IDL variable of type string, NP is a long, x, y and w are double arrays of size NP, s is a double, nest and n are longs, lambda and c are double arrays of size nest, fp is a double, lwrk is a long, wrk is a double array of size lwrk, iwrk is a long array of size nest and ifail is a long. The next step would be to call E02BBF, which evaluates the spline curve using the parameters calculated by E02BEF - specifically, the knot positions lambda(nest), and the spline coefficients c(nest) - see this demo for an illustration (albeit using a different package).

A more complete description of how to call external routines from within IDL is available here. That account also illustrates the use of IDL keywords to modify the behaviour of the function; their use has been eschewed in this example in order to retain the exact correspondence between the argument list for the NAG Fortran Library routine and the way in which it's called in IDL. We've also opted for a simple error-handling mechanism in the wrapper; alternatives which are better-integrated with IDL are possible. Finally, we note that more information on IDL DLMs is available from the IDL documentation, which can be found in various locations (such as here, for example).

Wandering Precision

2011-02-09T16:31:00.008Z

If you call the same numerical library routine twice with exactly the same input data, you expect to get identical results - don't you?

Well, it all depends. If you're calling a routine that somehow uses random numbers, for example a solver based on a monte carlo method, then the results might be different - but then you'd probably accept that (assuming that you'd been warned in the documentation).

But what if you're calling a numerical routine that in theory contains no randomness, but you still get slightly different results? And what if those slightly different results cause you problems - for example if you depend on repeatability for testing purposes? Then you might waste a lot of time hunting for the cause.

Well, this problem has turned up more than once recently for us at NAG, and in the cases we've looked at it has turned out to be a problem with the way the machine does its arithmetic - or rather, in the way that the compiler tells the machine to do its arithmetic.

Let me give you an example. Suppose you want to compute the inner product of two vectors x and y, each of length n. Simple - you just multiply each element of x by the corresponding element of y, then add all the products together to get the answer. No apparent randomness there. What could go wrong?

Well, if the hardware you are running on is equipped with Streaming SIMD Extension (SSE) instructions - and almost all modern hardware is - here's what.

SSE instructions enable low-level parallelism of floating-point arithmetic operations. For example, you can hold four single precision numbers at the same time in a 128-bit register, and operate on them all at the same time. This leads to massive time savings when working on large amounts of data.

But this may come at a price. Efficient use of SSE instructions can sometimes depend on exactly how the memory used to store vectors x and y is aligned. If it's aligned nicely - by which I mean, in the inner product example, that the addresses of the first elements of the arrays x and y are multiples of 16 bytes - then that's good. The hardware can efficiently move numbers from memory into registers to work on them, using instructions that depend on that alignment. So for our inner product, with a good optimizing compiler, we'd load numbers four at a time, multiply them together four at a time, and accumulate the results as we go along into our final result.

But if the memory is not nicely aligned - and there's a good chance it may not be - the compiler needs to generate a different code path to deal with the situation. Here the result will take longer to get because the numbers have to be accumulated one at a time. At run time, the code checks whether it can take the fast path or not, and works appropriately.

The problem is that by messing with the order of the accumulations, you are quite possibly changing the final result, simply due to rounding differences when working with finite precision computer arithmetic.

Instead of getting

  s = x1*y1 + x2*y2 + x3*y3 + ...

you get

  s = (x1*y1 + x5*y5 + x9*y9 + x13*y13) +
      (x2*y2 + x6*y6 + x10*y10 + x14*y14) + ...

Chances are that the result will be just as accurate either way - but it's different by a tiny amount. And if that tiny difference leads to a different decision made by the code that called the inner product routine, the difference can be magnified.

In some cases, for example when solving a mathematical optimization problem using a local minimizer, the final result can be completely different - and yet still valid - because a local minimizer might converge to any one of several different minimum points, if the function does not have a unique minimum.

And all because of how a piece of memory happens to be aligned. If you allocated the memory yourself you might be able to do something about it, but on the other hand the alignment might be completely outside your control. NAG users have reported cases where running a program from the command line would give consistently different results from running the same program by clicking on its name in Windows Explorer.

Is there anything we can do about problems like this? The fact is that NAG users want their code to run fast, so using SSE instructions makes sense. We can lobby the compiler writers to ask them to be careful how they optimize the code they produce, but they won't necessarily take notice. Compiler writers could also help by making their memory allocation routines only return 16-byte aligned memory, specifically to help with the use of SSE instructions. In the past, though, I had no success trying to convince the gfortran compiler folks to do that. In any case, even if they did, it wouldn't always help - if the first element of a single or double precision array is aligned on a 16 byte boundary, the second element will definitely not be, so if you want to operate on a vector starting at that location you've no chance.

We could choose not to use SSE instructions at all. But, apart from efficiency reasons, the "legacy" extended precision 80-bit register arithmetic which started in the 1980s with the 8087 numeric co-processor had its own "wobbling precision" problems for numerical code.

As new hardware with AVX instructions and 256 bit registers comes along, even more numerical work can be done in parallel. So - it seems that for the foreseeable future we're just going to have to live with this, and try to minimize the effects on NAG users by means of documentation.

The joy of creating new solvers...

2011-01-27T11:34:00.002Z

It’s an exciting time at the moment with a whole range of new routines being completed ready for the next Library release. The development group here are a highly motivated set of hard working folks, many with post-doc experience in a plethora of different areas of mathematics and computer science. While bringing their own individual knowledge to the code, they also get together as a well orchestrated team as all the methods are cross checked and fine tuned and documented for inclusion as new additions to the NAG library.

The next Library release will include more optimization routines, additional regression methods, extended wavelet functions, new surface fitting routines, further ODE solvers, enhanced random number generators... to mention only some of the content. So, particularly at the moment, it’s a real pleasure to have conversations with the individual developers and hear how some of the subtle design decisions that are made during the development process provide such large benefits, in usability, accuracy and performance, to users of our routines. For example one of the new optimization routines provides a new approach to solving bound constrained problems without requiring derivatives – making this solver very efficient for a range of large dimensional problems.

All these recent conversations have reminded me that we certainly have the expertise and the right attitude to help you. So do get in touch, at any time, if you want to discuss which of the hundreds of NAG routines you can use to solve your specific problem.

True Customer Engagement (when social is not so social)

2011-01-21T11:45:00.001Z

About a year ago NAG embarked on a new journey into the realms of social media as a new way of engaging with our customers and potential customers. For many years we’d relied on communicating via the website, email, telephone and letter (all of which we still do of course), but we felt social media outlets could provide a great new way of improving communication.

First we created the NAG Blog. The vision for the blog is to provide an outlet for NAG staff to share their thoughts with the outside world in an informal way. We also created a space for NAG on Twitter and a group on Linkedin. So far, we’ve been pleased with how the groups have grown with followers and readers, but feel what might be lacking with all three places is true interaction with our audience.

The statistics show that people are reading and following, yet there's little two way dialogue. I’ve been trying to understand why this is, I've concluded that it's probably down to a host of reasons. Here are my thoughts:

1. Maybe our audience don’t have time or inclination to communicate with us through social media?

2. Or, maybe we’re not asking the right questions of our audience, or in fact any questions?

3. Maybe our audience still prefer more traditional methods of communication?

During 2011 we’ll be thinking of ways to further enhance dialogue with our customers. We really do want to hear how you’re using our software and learn about your successes and failures. This helps us grow and improve as an organisation, not only by making product offerings more customer centric, but improving our customer service.

Please do get in touch if you have ideas or thoughts of your own on this. I’d love to hear from you.

NAG Life Service Recognition Award

2011-01-14T13:54:00.003Z

You don’t often feel like singing about being 40 years old except when you are a software company with such a rich history and as many good friends, colleagues and controbutors to remember. So it would seem very fitting in NAG’s 40th anniversary year to create the NAG Life Service Recognition Award. The award aims to recognise valuable contribution to the company over an extended period of time. An invitation went out last week to all members of NAG to submit nominations for this award. The award will be made at AGM35 this September. It’s a shame that we have lost touch with quite a few good friends, colleagues and contributors, I would love to hear from anyone who is interested in what we are up to now.

ParaView, VTK files and endianness

2010-12-22T10:32:00.012Z

Recently, I've been doing some work with ParaView, a rather nice open-source visualization application that's aimed at the display and analysis of large data sets. It has a wide range of visualization functionality, an interface that allows for interaction via its graphical user interface, or through scripting, and a distributed architecture. Executable versions are available for a variety of platforms, or you can download the source and compile it yourself. I'm currently running the Windows 64bit distribution on my laptop, and have also been building the source on HECToR, with a view to trying the application out on that machine, and using it for the analysis of some of the large data sets generated there.

ParaView uses the open-source Visualization Toolkit (VTK) for data processing and rendering; this consists of a core C++ class library and interpreted interface layers for Tcl/Tk, Java and Python. Much of VTK's functionality - specifically, its visualization techniques and modelling methods - is available from within ParaView, and it's also possible to extend the application to support other VTK classes by (for example) providing an XML description of the interface.

The first thing the user of any visualization system wants to do is to read their own data into it (of course, this will also be the last thing they ever do with the system, if it proves to be too recalcitrant). Doing this requires some understanding of the type of data that the system can process, and the translation of the different components of the user's data into those types. Generally speaking, it's usually most straightforward for the user to write their data to a file in a format which is supported by the system (the alternative approach - more useful when, for example, the the user has a lot of data files in a particular format - is to extend the system to support the reading of that format).

ParaView can read files in several formats, including so-called VTK Legacy files (so named because the format was introduced in an earlier version of VTK; since then, it's been supplemented by a more flexible XML-based format). Its structure is pretty straightforward, and well-documented. Here's the first part of an example Legacy file, which contains scalar data values located on a regular 3D mesh:

# vtk DataFile Version 2.0
3D example
ASCII
DATASET STRUCTURED_POINTS
DIMENSIONS 61 61 101
ORIGIN 0.0 0.0 0.0
SPACING 0.05 0.05 0.05
POINT_DATA 375821
SCALARS volume_scalars float 1
LOOKUP_TABLE default

The coordinates of the nodes on this kind of mesh (which VTK calls a Structured Point dataset) are fully specified by providing values for the parameters given by the keywords DIMENSION, ORIGIN and SPACING on each axis; the other file keywords are described elsewhere, along with an account of the other types of mesh and the other data structures which VTK supports. More specifically here, the keyword ASCII specifies that the data values (which come next in the file) are written in this format. Selecting this option ensures that the entire file is human-readable and portable from one machine to another, but it might not be the best choice if the data file is very large (for the usual reasons of space and speed). Switching to BINARY makes the file more compact and allows it to be written and read more quickly.

Here's where I ran into a (slight) problem: although ParaView read the the ASCII version of my file and displayed it correctly, this broke when I switched to the binary version (more accurately, ParaView read the file in happily but the values which it found there were wrong). After some poking around in the documentation, I realised the problem was that VTK - and hence, ParaView - always expects the binary data in the Legacy file to be stored with the most significant byte first - also known as the big-endian representation. (To be strictly honest, I couldn't find this point specified unambiguously in the otherwise-excellent VTK file formats guide, but it was emphasised more explicitly elsewhere). However, my file had been output by a program in which the data was stored with the least significant byte first - i.e., in the so-called little-endian fashion. The representation used by the program to store the data depends on both the operating system and hardware architecture for which it's been compiled (for example, it's little-endian for Windows on x64 - which is what I'm currently using - but big-endian for Solaris on SPARC), and this representation is preserved when the data is output as a raw stream of bytes - as is done, for example, by the C fwrite() function.

Having realised the problem, the solution was clear: swap the order of the bytes around in each data value before writing it out. There are many ways to do this, but I found this document helpful. It proposes a framework for handling several data types which is intended to be portable across both little- and big-endian architectures (note also that the framework receives much attention here as alternatives which are possibly more robust, general and/or portable are suggested). Since my data is stored as floats, I only needed a small part of the framework; here's the C function I used:

float FloatSwap( float f )
{    
   union
   {
      float f;
      byte b[4];
   } dat1, dat2;

   dat1.f = f;
   dat2.b[0] = dat1.b[3];
   dat2.b[1] = dat1.b[2];
   dat2.b[2] = dat1.b[1];
   dat2.b[3] = dat1.b[0];
   return dat2.f;
}

and here's how I used it when writing out the data values:

float *data, val;
FILE *fp;

/* Generate the data values, open the file and write the header.  */
...

/* Output the data values.  */
index = 0;
for( k=0; k<nz; k++ )
   for( j=0; j<ny; j++ )
      for( i=0; i<nx; i++ ) {
         /* Swap the bytes in this data value, write it out.  */
         val = FloatSwap(data[index]);
         fwrite((void *)&val, sizeof(float), 1, fp);
         index++;
      }

fclose(fp);

Finally, here's the visualization which I created in ParaView. The dataset is the probability density function corresponding to the so-called 3d_z² atomic orbital, and it's visualized using volume rendering and an isosurface, which highlights its characteristic shape: two lobes and a torus (I once heard the profile of rather overweight chemistry teacher being compared unflatteringly to this).

Using ParaView to visualize a 3D dataset

Half-way through HECToR

2010-12-16T10:55:00.005Z

HECToR, the UK national supercomputer service, had its third birthday in October. NAG, as many readers will know, provides the Computational Science and Engineering (CSE) support for the service, helping users with application problems, and with porting and tuning their codes to make them as efficient as possible. We also provide an extensive programme of training courses throughout the UK covering both basic and advanced topics. To date we have had over 900 attendees on these courses and delivered them in 16 locations. We are currently putting together our programme for next year, upgrading our course materials in response to the latest hardware upgrade (from a Cray XT6 to an XE6), and developing material for some new topics that we haven't addressed before.

A novel aspect of the HECToR service is the Distributed CSE Service (DCSE) which funds dedicated resources to work on specific codes. Those resources can come from within the research community itself, or from specialist teams (including our own HPC team), and to date we have funded more than 46 years of effort for 48 projects.

DCSE projects have addressed a wide range of issues, but two themes recur frequently:

Adopting shared-memory techniques. This is necessary to share data between multiple cores on the same socket where there is not enough memory for each core to have its own copy, and also to make efficient use of the increased number of cores per node.
More efficient I/O. Reading and writing data is often a major bottleneck in applications and parallelising the I/O and using libraries that compress data efficiently can deliver impressive performance improvements.

In the last three years HECToR has gone from 2 cores per node to 24, and from a total of 11,328 cores (on the original XT4) to 44,544 (on the new XE6). However memory per core has dropped from 3GB to 1.33G, and clock speed has dropped from 2.8GHz to 2.1GHz. These changes mirror those that we are seeing more generally in our industry, namely more cores running more slowly and with less memory. Our experiences on HECToR demonstrate that software needs to be adapted to make efficient use of newer hardware. The good news is that this sort of software engineering work is often very cost-effective - many of the DCSE projects have saved more money in HECToR resources than they cost in labour. Hopefully this will lead to more and better scientific outputs in the future.

Christmas Comes Early

2010-11-25T14:01:00.001Z

Christmas came early for one little girl at NAG Ltd. Being in charge of office services has its ups and downs (hands down a blocked loo or physically being up a ladder changing light bulbs) are part of my normal day. But a call from the local fire brigade certainly cheered me up. Since the fire service has stopped its annual building inspections we are now reliant on paying somebody to come in every year to carry out a risk assessment and for some reason these people are always rather bland and lacking in making their work (although serious) fun. Anyone who was here two years ago and was lucky enough to have training on the use of fire extinguishers would know this has to have been the best training course ever; it was a laugh from start to finish. The request was to familiarise themselves with the layout of the building so if there was ever a real emergency and they would have some background knowledge. They had promised to arrive in a full rig and I had expected one or two officers, but four was hitting the jackpot! We had a nice chat over a quick cuppa about fire exits (staff should use the designated fire exits and not the main entrance), number of hydrants on site, lighting etc and then we did a tour of the building. Being only 5ft 1 ½” (5ft 6” in heels) I felt really small as they were all over 6ft but they were a really nice bunch who do a fabulous job and hopefully will only ever be here is a more social capacity rather than a fully fledged work one.

Now you see it...

2010-11-22T17:06:00.041Z

A few weeks ago, I presented the HECToR scientific visualisation training course to a group of researchers at the University of Reading. This course looks at the use of visualisation as a tool for the improved understanding of numerical data (as produced, for example, by calculations or simulations run on HECToR). We started by exploring the different types of data - characterizing each according to its structure, dependencies and dimensionality - before reviewing the different visualisation techniques (such as contouring, particle tracing, volume rendering) that are applicable to each type. Some examples of the techniques, used to display a variety of data types from several application areas, are given in the figures below.

Figure 1. Fluid flow around a cylinder

Figure 2. Visualisation of MRI data, showing an aneurism

Figure 3. Potential energy in the space between three atoms

Next, we briefly reviewed a few applications for performing visualisation, of which there are a plethora in both the public domain and commercial arenas. Some of these have comparable functionalities - thus, for example, Figure 1 was produced using IRIS Explorer, ADVISE was used for Figure 2, and Figure 3 was created with ParaView. The plot in Figure 4 (which illustrates the use of a routine from the NAG Toolbox for MATLAB®) was made using MATLAB.

Figure 4. Solution of a Sturm–Liouville equation

Figure 4 can be used to demonstrate another aspect of the training course - namely, the highlighting of good and bad practice in visualisation. Depending on the resolution of the display medium for the figure, it is usually apparent that the blue solid curve is easier to see than the dashed green one. This may not necessarily be a significant effect, although you might be reminded of other instances you've seen in presentations where - for example - a yellow curve has been rendered more or less invisible by being displayed against a white background, or a blue curve against black. In general, it can be surprising how little attention is apparently paid to issues such as clarity and reproducibility in visualisations by users who've spent a lot of time generating and checking their results, only to see their impact lost or impeded because of a poorly-designed display.

As alluded to above, one of the things that should be taken into account when designing a visualisation is the use of colour. Indeed, Figure 3 also illustrates how care should be taken when selecting colours: many viewers will find its juxtaposition of the red surface against a green background difficult to look at (try reading green text on a red field - or vice versa - for a further illustration of this). The way in which we perceive colour is a complex subject, but one salient point is that, whilst our eyes see absolute colours, our brains perceive differences in colours - i.e., the appearance of a colour depends on its surroundings. This apparently surprising fact underpins the workings of a whole range of optical illusions, of which my current favourite is the so-called same colour illusion:

Figure 5. Are the two orange circles the same colour?

Looking at Figure 5, it can be difficult to see that the two orange circles are the same colour - although, bearing in mind the discussion above, most people would probably accept it. But what's much harder to believe (although it's also true) is that the squares which surround them are also the same colour as each other. The reason we perceive them as different is that they have different surroundings (and also because of the way in which our brains try and compensate for the effect of shadows).

Demonstrations of illusions such as this (and the examples of good and - mostly - bad practice in visualisation) made the training session quite a lively one. In fact, I suspect that - notwithstanding the undoubted fascination that can be exerted by the study of visualisation techniques and applications - this was the part of the presentation that lingered longest in the minds of many of the attendees. I could be wrong, though.

Fortran Fame

2010-11-17T16:19:00.007Z

I was standing on the NAG booth here at SC10 yesterday, muttering to myself about the difficulty of maintaining an internet connection in the exhibition hall. The organisers claim that there's more bandwidth available at SC10 than is used by any but the five largest countries in the world - I think they're exaggerating because I can never get any of it.

Anyway, along came a dapper fellow wearing a baseball cap who started a converstaion about the way standards were being implemented in the NAG Fortran compiler. I asked whether he was interested in using Fortran. "Well, I used to be", he said. Belatedly I looked down at his name badge. It turned out to be Tom Lahey, founding father of the Lahey Fortran compiler. What could I do but shake his hand and bow.

It's a hard life.

2010-11-17T06:50:00.013Z

Every year in November NAG staff attend SC, formerly known as SuperComputing, an annual showcase where computer hardware and software vendors can display their goods and services, along with a huge technical agenda covering all aspects of HPC and scientific computing.

The show moves around different cities in the United States. This year SC10 is in New Orleans, Louisiana, and running from 14th to 21st November, with around 10,000 registered attendees. It's the 23rd outing since the first one which took place in Orlando, Florida in 1988. NAG has attended every single one of these events, and is one of just a handful of outfits to do so.

New Orleans was founded in 1718 by French speculators interested in making money from trade and finance, and named after the then French Regent, Philippe, Duke of Orléans. In the mid-eighteenth century it was ceded to the Spanish. They're the ones who actually created most of the characteristic architecture in the heart of the city now known as the French Quarter. In 1801 it returned to French control, only to be sold to the United States a couple of years later as part of the "Louisiana Purchase". In more recent times New Orleans was badly hit by Hurricane Katrina, and there's a nice museum detailing all the sad events of that calamity next door to St. Louis' Church in downtown New Orleans.

So, SC is always a great place to find out about the latest developments in computing technology, talk to NAG customers, meet old friends and make new ones, and this year is no exception. By the way, if it's not too late and you want to come to talk to us, the NAG booth is number 3131 in the exhibition hall.

The trouble is the parties. They start on the Sunday evening, with the "Exhibitors' Party". This is open to anyone with an X on their name badge, meaning that they're working one of the booths in the exhibition hall. The food is piled high, and there's enough drink to float the Queen Mary.

Monday evening is the "Opening Night Party", which is at the nearby Metropolitan Night Club. A jazz band plays the hits of Fats Domino, Louis Armstrong and other New Orleans greats. The food is piled high, and there's enough drink to float the Queen Mary.

Tuesday evening. Two big parties going on, one hosted by Cray and AMD which takes place in a ballroom at the Marriott Hotel, and one hosted by SGI and Intel, which is at Pat O'Brien's Bar on Bourbon Street in the French Quarter. In the interests of vendor neutrality, of course NAG staff must attend both parties. O'Brien's has another great jazz quartet playing, and is famous for "The Hurricane", a drink containing rum and a secret mix of fruits and spices. It's lethal, and guess what, there's enough of it to float the Queen Mary.

There are more parties to come. So, although we have the compensation of being able to learn more about what NAG users want in the way of fast numerical software, and finding out what our contemporaries are working on, the fact is, we've all got blisters on our feet from walking between the different parties.

It's a hard life. But someone's got to do it.

Looking forward to meeting people with new ideas at Supercomputing.

2010-11-11T10:59:00.003Z

I’m someone who enjoys talking with people at exhibition booths and the International High Performance Computing event, which is being staged next week (http://sc10.supercomputing.org/ ), will provide good opportunities.

Here is a chance to spend time with many current and future users of NAG products and services as well as with NAG collaborators. Attendees at this conference and exhibition work in such a variety of different areas that the conversations are always fascinating and informative.

Meeting with people face-to-face inevitably triggers more new ideas, and results in innovative solutions being created faster, than any electronic communication method.

I look forward to meeting some of you there (NAG have Booth #3131 – right at the heart of the main hall).

Creative Programming

2010-11-05T09:44:00.007Z

As I was being asked to write this blog post my brain awoke from it's number crunching stupor and I realised that I couldn't think of anything to write. The past few hours of sitting, programming in front of a computer seemed to have wiped out the artistic half of my mind.

I know at least for myself that staring at a monitor seems to obliterate all of my thoughts until I've finished the click clacking and click clicking of my current task and I can finally put blessed ink to physical paper to plan my next step.

So in order to re-enliven my thinking and my creativity, I sought a nice picture of a stunning scene. There's nothing like the beauty of creation to awaken something deep in your heart.

As I was quickly sketching the above picture to get the other half of my brain warmed up to it's full creative potential, I started wondering about creativity in programming.

Firstly there is the beauty of the written code itself, the ideas encapsulated within it and the intelligence of the design. Secondly there is the created product of what the non-technical viewer might see, such as a website, 3D artwork, mathematical art, 3D movies etc.

There is creativity in all that we do, and no less in programming. It is very satisfying to have finally cracked a problem successfully, and even more so if you can look at the flow and structure of your program and know that it is good and pleasing.

But despite the beauty that can be found in computer generated output or even in code (if you are able to comprehend it) I think I will always have to turn to a beautiful, fiery, ruby red sunset, or the playful sound of a nearby splashing stream as a way of awakening my heart.

(photo kindly contributed by Paul Wilson from his flikr account)

Comparing HPC across China, USA and Europe

2010-10-30T00:25:00.000+01:00

In my earlier blog post today on China announcing the world's faster supercomputer, I said I'd be back with more later on the comparisons with the USA, Europe and others. In this morning's blog, I made the point that the world's fastest supercomputer, in itself, is not world changing. But leading supercomputers, critically matched with appropriate expertise in programming and using them, togther with the vision to ensure use across basic research, industry and defence applications can indeed be strategically beneficial to a nation - including real economic impact.

There are plenty of reports and studies describing the strategic impact of HPC within a given organisation or at national levels (some are catalogued by IDC here), so let's take it as a premise for the following thoughts.

With this in mind, there are some comparisons to be made between the approaches to supercomputing across the USA, Europe and China.

The USA has long enjoyed near total dominance of the hardware technology underlying the leading supercomputers. The USA has invested repeatedly in ensuring that American supercomputer manufacturers have the technology to deploy the world's largest supercomputers. The last time the USA lost the public leadership crown of fastest supercomputer, a huge investment in Amercian supercomputer technology followed. As well as national support for the development, the strong implicit requirement for USA organisations to "buy American" ensures a continued USA supercomputer manufacturing industry. As a result of the sheer size of the country, the USA has a large HPC R&D community. There is also significant usage of HPC in American industry and additional support for this through government initiatives like INCITE.

Despite the recent growth of Bull, and their highly-rated supercomputers and in-house HPC expertise, Europe has not seen it essential to have a home-grown supercomputer manufacturing industry (beyond the component level R&D which is strong in Europe). Europe has always highlighted its expertise in applied HPC software development and in software-related HPC R&D as its distinguishing strength on the international stage. European organisations have always bought supercomputers from around the world - and since they are likely to continue to enjoy access to products from around the world, have no need to develop independent supply from wholly European sources.

But the news of China's investment in hardware, software and people (and stated ambition of independence of supply) should make a clear message to USA, Europe and others that they cannot rely on their continued leadership of HPC. And thus the economic benefits of HPC might soon be driving Chinese growth rather than the European or American economies.

There is one other potentially killer advantage that China might have. All the predictions of the path to the next levels of supercomputer performance (Exascale) involve major changes in the technology - much greater levels of parallelism (seeing this now with GPUs), memory performance challenegs, resilience, etc. As Dave Turek of IBM has said before, China is mostly not hindered by legacy code - they can start anew with the best HPC ideas and methods of today, looking to the future. In USA and Europe, our obsession with "protecting our investment" in established applications, means we first have to figure out how to get from yesterday's software technology to current methods, then to the future. Is "protecting our investment" actually constraining our future?

[More on legacy code and revolutions vs evolution in a blog coming here soon ...]

Why does the China supercomputer matter to western governments?

2010-10-29T12:16:00.002+01:00

There is a lot of fuss in the mainstream media (BBC, FT, CNET, even the Daily Mail!) the last few days about the world's fastest supercomputer being in China for the first time. And much ado on Twitter (me too - @hpcnotes).

But much of the mainstream reporting, twitter-fest, and blogging is missing the point I think. China deploying the world's fastest supercomputer is news (the fastest supercomputer has almost always been American for decades, with the occasional Japanese crown). But the machine alone is not the big news.

Imagine that China announced a new prototype passenger aircraft, half the cost of the latest Boeing or Airbus. It has 50% greater fuel efficiency too. And an order of magnitude greater predicted reliability statistics. That would be major news. Sure it uses a lot of US designed components too.

But what if China announced this new aircraft wasn't just a prototype. It was a commercially available product now. And they have the capacity to make lots of them - faster than Boeing or Airbus. And they have a plan to train huge numbers of future aircraft maintenance engineers, aerodynamic designers, structural engineers, etc. In other words, China can not only build a world-beating passenger aircraft, but it is building the capability to do so without US designed components in the future. And it is building the expertise capacity to be a world leader in aircraft maintenance.

That would be very important.

And, while we are not quite there yet, that is where this China supercomputer news is going. American scientists and HPC professionals have been calling for rounded investment in people and software not just hardware for years. Europe has been proud of it's relative HPC software expertise, but the recent IDC led EU HPC strategic recommendations report shows that much more investment is needed.

The mass commentary talks about the Chinese hardware milestone. But public material from Chinese experts also talks about a plan to deploy several top supercomputers, to train huge numbers of HPC programmers, to invest in applications development and commercial use of HPC and to develop end-to-end nationally independent supercomputing technology.

If that happens, then China will have the ability to develop that super aircraft industry. And automobile. And the many household products that are designed with supercomputers. And materials science. And ...

Get the point? It's not the world's fastest supercomputer that matters most. It's not just national pride. It's the ambition and comphrensive plan behind the world's leading supercomputer that matters.

[More later today on the comparison with USA, EU and others]

Carbon Footprint

2010-10-22T08:58:00.001+01:00

In these austere times I am always on the lookout for ways to reduce our expenditure without reducing our service to customers. NAG have just been awarded a matched funding grant from http://www.sustainableroutes.co.uk/ for £1000 to help us reduce our carbon footprint and we intend to use the grant (and some) to improve our video conferencing and cut down on some travelling costs and CO2 for meetings and training. Currently we are looking at some state of the art video conference equipment which we hope to start trialling with our Manchester office. There are some other useful tips and grants available on the Sustainable Routes website so I would recommend go taking a look.

Source-level debugging of Python in Emacs

2010-10-19T08:00:00.000+01:00

To help me track problems in Python code (yes, even sometimes in my own...) I usually rely on good-old print/trace debugging. Owing to Python's speed of interpretation and execution, this is a pretty convenient approach. Python does have its own interactive debugger though— pdb—for those odd occasions when it's desirable to poke about in a program while it's running. The debugging mode in Emacs even supports pdb by default, but there's a little snag: you probably don't have a pdb in your path, so M-x pdb will just fail. Solution? Add a pdb script to your path and make it executable

#!/bin/sh
pdb_path=`python -c "import pdb, sys; sys.stdout.write(pdb.__file__ + '\n')"`
exec python ${pdb_path} "$@"

2010: A Retail Store Odyssey

2010-10-07T04:24:00.001+01:00

Stanley Kubrick’s film “2001: A Space Odyssey” was the epic 1968 science fiction film that explored human evolution, technology and artificial intelligence with both realism (and surrealism) and remains one of the top films of all time. In it, two astronauts battle the computer HAL for control of their spaceship and for their lives while investigating a series of strange monoliths left from an earlier civilization. For many years and for many people, the film has been symbolic of our struggle to master, and not be mastered by, computers. Kubrick and his co-author Sir Arthur C. Clark were both brilliant and far ahead of their time. In many respects, they still are.

In 2010 we are certainly wrestling with computers that occasionally seem to get the better of us. In our world, computers are ubiquitous and the software behind them has a pervasive impact on our daily lives but hardly in the way Kubrick and Clark envisioned. Consider this absolutely mundane sequence of events, at least in the developed world. We hop into our car to go shopping for groceries. As we turn the key one or more microprocessors start, employing sophisticated software to optimize efficiency, performance and emissions. Our cars talk to us, connecting phone calls routed through a cellular network. They entertain us with hundreds of channels from a satellite radio connection and give us visual and voice directions to where we want to go. All of this amazing hardware comes to life through the software that makes it work. And, of course the software makes considerable use of mathematics to accomplish what it does for us. So, you may ask, what’s the reference to a “retail odyssey”? We haven’t even gotten to the grocery store yet.

In my view, the most amazing thing at the Sainsbury I frequent in the UK or the Dominick’s near home isn’t the check out where they let me scan my items and coupons and pay with my credit card, all with a few taps of the touch-screen. It’s the realization that I’ve just walked through a store with literally thousands of unique items to meet my needs, each residing in a database linking the bar code on the package with a price, an inventory level, cost, supplier and even a “loyalty card” database that permits analysis of which shoppers bought which products and in which combinations. While you are pondering this miracle of modern technology, ponder this question: who set the price of the 2-pound package of Folgers coffee at the end of Aisle 2 and how did they do it?

The answer, if it’s not already obvious, is sophisticated software from companies like NAG partner DemandTec (NASDAQ: DMAN) whose demand management software is helping retailers and manufacturers worldwide optimize revenues, prices and inventories. We’ve worked with DemandTec since 2004, providing them sophisticated mathematical and statistical software to enable their application to help retailers manage demand.

One of the benefits of working with cutting edge companies like DemandTec is that we get to participate in events such as I have been involved with lately. NAG has partnered with DemandTec to sponsor the Chicago Regional scholarship competition in a national event called the DemandTec Retail Challenge. In it, high school seniors in teams of two get to apply their problem-solving and mathematical skills as pricing analysts managing an assortment of products in a grocery store. They set the price, order inventory and make decisions about promotions in a simulation-based competition with other teams in the Chicago area. The eventual winners will have maximized profitability and successfully communicated their approach to the problem to experts in the field. The winners earn a college scholarship and the right to compete in the national version of the contest at NASDAQ in New York City in January 2011. The national champions get a significant additional scholarship and the right to ring the closing bell for the trading day. For those of us at NAG it’s both a way of giving back to the community and helping the next generation apply their academic skills to real-world problems. From the conversations I’ve had with them thus far I suspect that the computer HAL would be no match for them.

The Shoulders of Giants

2010-09-30T12:02:00.000+01:00

Last week NAG held its 34th AGM and prior to the business part of the meeting we had the pleasure of hearing a lecture by Professor Nick Higham of the University of Manchester on the subject How and How Not to Compute the Exponential of a Matrix. The lecture was filmed and can be viewed here. Matrix functions are useful in a number of areas, in particular for describing the solutions of certain types of differential equations. NAG is working with Nick to include his algorithms in future marks of the Library.

One of the things I particularly enjoyed about Nick's talk was the way that he set the subject in its proper historical context, starting with the work of 19th Century mathematicians such as Sylvester and Cayley. Working at the cutting edge of computing technology its easy to forget how much of what we do depends on work that goes back hundreds of years. It wasn't always this way: Fermat, having claimed to have found a proof of what became known as his last theorem, added the waspish comment "and perhaps, posterity will thank me for having shown it that the ancients did not know everything".

A few years ago SIAM interviewed a number of prominent people who had been active in the early days of numerical analysis and scientific computing, including NAG's own Brian Ford. You can find transcripts of the interviews and related material on SIAM's website. Reading these transcripts is illuminating, particularly to a (relative) youngster like me, especially given the number of references to NAG's early history that they contain. I hope that SIAM keeps up this activity, and continues to collect personal reminiscences from prominent people in the field.

HTML5: Possible implications for technical documentation

2010-09-23T10:06:00.002+01:00

A New Version of HTML

As you may have heard, HTML5, a new version of HTML, the main markup language for the web is under development. This is the first version for quite some time, HTML4.01, the current version, was released in December, 1999. There was an XML version, XHTML 1.0, XHTML (Released in January 2000 and revised August 2002), but it did not contain any new features.

Much of the buzz around HTML5 concentrates on new features such as video (as an alternative to using flash plugins etc) and canvas (a JavaScript API for 2 dimensional graphics), but in this article I want to look at the features that impact on the kind of technical, mathematical documents used by NAG.

SVG and MathML

SVG is an XML format for describing scalar graphics, using a similar graphics model to PostScript or PDF, but using XML syntax and CSS for styling, just as for HTML. It has been around for some years, but HTML5, for the first time, specifies how it works in the context of HTML (rather than XML). In a related move, Internet Explorer 9 will support SVG. (Other common browsers have supported it for some years.)

Many of the function plots in the NAG documentation are generated from gnuplot. Gnuplot can save to SVG as well as the EPS and PNG outputs that we currently use. Also the remaining diagrams are currently stored as EPS (for use in PDF versions of the document) and EPS can usually be converted without any loss of information to SVG.

There are many advantages to using a scalable format for plots, especially documents that are read on smaller devices such as tablets (or phones), but also printing at higher resolution; being able to scale and zoom in as required without the loss of quality inherent in scaling bitmap images is a big win.

The other XML-derived language supported natively by the HTML5 parser is MathML. MathML should be familiar to most NAG users as it has been used in the XHTML version of our documentation for years. However being able to use MathML in HTML should hopefully simplify the installation of this documentation, and also the prominence given to MathML by its inclusion in HTML5 will also hopefully encourage the remaining browser vendors, who do not natively support it, to implement mathematics rendering in their browsers.

Current State of (Future) Browsers

All the common desktop browsers are moving towards supporting HTML5 in the near future. The following sections list the state of SVG and MathML rendering in some common cases.

Internet Explorer 9 Beta

IE 9 finally adds SVG support to Microsoft's Internet Explorer browser. Design Science's MathPlayer may still be used to get very high quality rendering of MathML. There are currently some problems in interfacing MathPlayer and IE 9, and some of the markup that worked in IE 8 no longer works, however this may be fixed before the full product comes out, or a few lines of IE-specific javascript, as used in the example below, may be used to work around the main problems.

Warning: IE 9 Beta, unlike all the other test releases of browsers mentioned here will replace your existing Internet Explorer installation. If you wish to test IE 9 rendering without removing your IE 8 (or earlier) installation make sure that you install the Platform Preview rather than the Beta. The platform preview releases use the IE 9 rendering engine but without all the normal browser menus and facilities, and it does not replace your existing browser.

IE 9 Beta 1 was tested.

Firefox 4.0 Beta

Firefox 4 uses the new HTML5 parser which automatically places MathML elements in the MathML namespace and SVG elements in the SVG namespace. The rendering of MathML and SVG in HTML 5 is then essentially exactly as has been available in earlier Firefox releases in XHTML documents, but now also available in HTML.

Firefox 4 Beta 6 was tested.

WebKit Nightlies

A Nightly release is a public binary release that is essentially just a snapshot of the current development build. As such it is expected to have bugs and unimplemented features. However these releases, which show the development state of the core rendering engine used by both Chrome and Safari browsers show that it has a reasonably well developed SVG renderer, and that MathML support, while very new, is improving. (These WebKit browsers have very good CSS support so until the native MathML is ready for use it will probably be possible to use a CSS rendering of MathML in these cases, as is done for the current NAG XHTML documentation).

WebKit-r67637, from 16th September was tested.

An Example

We do not yet have any definite plans or dates to move towards using HTML5, it depends a lot on the timing of the above test releases becoming mainstream (which is likely to be very soon) and how long it takes for use of older browsers to die out. (This can take a long time, many people only change their browser when they get a new machine, Internet Explorer 6 (released in 2001 and replaced by IE 7 in 2006 and IE 8 in 2009) is still in widespread use for example.)

As an experiment we have converted one example routine document into HTML5, using inline MathML and SVG. Note in the example all links to external documents have been made to point back to the test document so that it is self-standing. The document has 207 MathML fragments, which are essentially unchanged from the existing XHTML+MathML version of the document, the two plots shown in the example section, which were previously only available as monochrome scalable images in the PDF version of the document, or a coloured PNG bitmaps in the XHTML version, are provided as scalable SVG which is directly included into the HTML file, not referenced using <img src=....

The document is reported as valid HTML5+SVGMathML by the validator.nu online validator.

The document appears to be perfectly usable in Firefox 4 and IE 9. The MathML support in WebKit is still very new, and there also appear to be some problems with the SVG which get rendered with a black background in the Windows nightly build of 16th September. Note however that, as their name suggests nightly builds are unstable test releases and should not be taken as indicative of the final product. The fact that the SVG was fully rendered apart from some colour problems, and that the MathML was recognized, albeit with some rendering problems, is I think a good sign that WebKit based browsers (Safari, Chrome, and many mobile phone browsers) will support these documents in the near future.

Fortran Library Mark 22 Routine Document (D02AGF)
XHML+PNG+MathML	Mark 22 Documentation	d02agf.xml
PDF	Mark 22 Documentation	d02agf.pdf
HTML5+SVG+MathML	Experimental	d02agf.html

Conclusions

HTML5 will be mainstream in currently popular browsers very soon, it will take a few years before its use may be assumed, but it promises to be a big improvement for technical documents.

Both MathML and SVG allow the use of scalable formats freeing the web from the use of inappropriate bitmap images. This simplifies document distribution and should improve the quality both for printing and for use on a wider range of devices, including greatly increased accessibility.

Do you want ice with your supercomputer?

2010-09-13T13:41:00.000+01:00

“Would you like ice with your drink?” It’s a common question of course. One that divides people – few will think “I don’t mind” – most have a firm preference one way or the other. There are people who hate ice with their drink and those who freak if there is none. National stereotypes have a role to play – in the USA the question is not always asked – it’s assumed you want ice with everything. In the UK, you often have to ask specifically to get ice.

Yet the role of ice in making our drinks chilled is misleading. I once had a discussion with a leading American member of the international HPC community about this. “No ice”, he was complaining as we headed out of a European country, “they had no ice for the drink”.

“I don’t get this obsession with ice”, I chipped in. “What?!” He looked at me as if I were mad. “Why do you like your coke warm?”

“Ah, but that’s just it”, I replied. “I hate warm drinks – I really like my coke chilled. But surely, in this modern world over a century after the invention of the refrigerator, it’s not unreasonable to expect the fluid to be chilled – without the need to drop lumps of solid water into it?”

“Ah, fair point”, he conceded.

What has this got to do with supercomputing? Perhaps the common thread is that usually we just accept the habitual choices of ways to do things – and don’t often step back to think – “are those the only choices?”

Maybe we should step back a little more often and ask ourselves what we are trying to achieve with HPC – and are the usual choices the only ways forward? Or are there different ways to approach the problem that will deliver simpler, better or cheaper performance?

Perhaps your business/research goals mean you need to conduct more complex modelling or you need faster performance. Maybe the drive of computing technology towards many-core processors rather than faster processors is limiting your ability to achieve this. (I have had several conversations recently, where companies are buying older technology because their software won’t run on multicore).

The “ice or no ice” question might be whether or not to upgrade your HPC with the latest multicore processors. But what about the “just chill the fluid” option? Well, how about upgrading the software instead, or as well?

NAG has plenty of case studies to show where enhancements to software have achieved huge gains in performance or capability (e.g., www.hector.ac.uk/cse/reports).

Sometimes buying more compute power is the right answer. Sometimes, extracting more efficient performance from what you have is the answer. Bringing them together - a balance of hardware upgrades and software innovations might well give you the best chance of optimising cost efficiency, performance and sustainability of performance.

Working on the ADVISE project

2010-09-08T16:37:00.034+01:00

For the past three and a half years, my colleagues and I have been working on ADVISE, a TSB-funded collaborative research project which has been developing a new toolkit for visualization and analysis. Besides NAG, the partners in the project were VSNi and the University of Leeds. VSNi have expertise in statistics, as implemented in their GenStat product, while Leeds have an international reputation for their work in visualization research. As for NAG, we've had some success with IRIS Explorer, a popular visualization toolkit which allows users to construct applications by connecting modules together via a visual programming interface.

Figure 1. The ADVISE pipeline builder

We retained that interface in ADVISE (see Figure 1) because it has proved to be a rather intuitive way to create, modify and interact with applications. Thus, in this figure, the user selects modules from the repository on the right and connects them together in the area on the left. The widgets for controlling one of the modules are in the pane on the right at the bottom, whilst messages from the system are displayed in the area at top right.

In a similar spirit of re-use, the visualization and analysis functionality encapsulated within the ADVISE modules has come from porting just about all of the modules from IRIS Explorer into the new environment, and creating new modules that generate and process scripts of GenStat commands. We've used ADVISE to visualize and analyze a variety of data - see, for example, Figure 2, which is a display of some of the results from Christoper Goodyer's simulations of diffusion through skin.

Figure 2. Visualization produced by the final module in the pipeline of Figure 1

So much for the re-use of old technology, but what's new in ADVISE? Well, its architecture makes use of recent technology developments in web services and distributed computing. This has several advantages, including the fact that it's easier to integrate applications built using ADVISE with the web (for example, running inside a browser), and that it's possible to connect ADVISE applications to other services (which could, for example, act as data sources).

Figure 3. Using ADVISE to create a browser-based visualization of air quality data

To illustrate this ease of integration, Figure 3 shows a web-based application that's been created in ADVISE for the visualization of air-quality data. The window at the back shows the interface for selecting the location and duration of the data to be visualized, the next window shows that data displayed as a coloured elevated surface and the window in front shows the same data displayed as a 3D histogram. Widgets in the web page (linked to ADVISE modules) give the user a simple interface to the application - for example, allowing control over the type of display (surface or histogram), other parameters associated with the visualization, and selection of the next dataset to be displayed.

If you want to know more about the ADVISE project, or the system we produced, head over to our web page, which contains more pictures of visualizations created with the toolkit together with papers, posters and talk slides from throughout the life of the project. One picture you won't find over there, however, is the one below, which shows the whole ADVISE team (with the exception of Jungwook Seo and Colin Myers from Leeds) in all their glory at the end of the final project review meeting last month. There were (a lot of) other pictures taken at the same event, but this is the only one in which the project members aren't holding glasses.

NAG at Quant Congress, New York

2010-08-27T15:10:00.011+01:00

Last month I was in New York with my sales colleagues Mike Modica and Rick Guido from NAG's US office, attending Quant Congress USA. This is a meeting which is devoted to the latest developments in financial derivatives, risk management and the associated use of numerical techniques; NAG has been associated with it for the past five years because our software has found extensive use in the quantitative analysis community. Each day of the conference program opened with a plenary session of two presentations, followed by around twenty talks divided across two streams. The meeting had nearly a hundred attendees, not all of whom attended both days (perhaps it's worth noting that, whilst the technical program undoubtedly had its own attractions and merits, the attendance might have been enhanced still further by the fact that the two days of the meeting saw heavy rain fall on the city - in the middle of an otherwise-sweltering week). NAG was one of five exhibitors at the meeting, and we found that the relatively small size of the event meant that everyone who wanted to view our stand - almost exclusively during the breaks from the talks - was able to.

We were kept busy by a steady stream of stand visitors throughout the event, which made for an interesting couple of days. Most of the queries and comments we received were focussed on technical queries about the functionality of the NAG libraries, with a particular interest, of course, in those routines which have proved valuable in the finance industry. There was also a lot of interest in the work that NAG has been doing on numerical routines for GPUs, partly because of the speedup that has been observed by customers when running Monte Carlo simulations on GPUs - as compared to running on a single CPU - using NAG routines (see this post for more details).

We also had a lot of questions about the NAG Toolbox for MATLAB^®, while some users - or prospective users - were interested in being able to call NAG routines from within Python. I had an interesting discussion on this latter topic with Professor Dennis Allison (who is famous for, amongst other things, being one of the founders of Dr Dobb's Journal). I found that providing answers to a few of the most technical queries about numerics in quantitative finance was somewhat beyond my abilities, notwithstanding the recent education I'd received in some aspects of this field from my colleague Marcin Krzysztofik (and our collaborations involving the NAG option pricing routines), but we were able to deal with all outstanding issues by email within a few days of the conference's end.

Whilst in New York, we also made a series of site visits, mainly in the financial sector. During one of these, I was told by the customer that they thought the technical support which they received from NAG was excellent - in fact, they said it was often better than what they received from their own in-house team. Since the support service is one of the things that NAG prides itself on (and which the customer pays for), this was good to hear, particularly as I felt it made up for some less gratifying - but still memorable - experiences I'd had earlier in the trip. These included celebrating my arrival in the city by falling off a barstool - apparently without any provocation at all (apart from the fact that the barroom singer was playing a cover of Tom Petty's "Free Fallin'" at the time). Any damage that my self-esteem may have suffered was, however, assuaged by the (re)discovery of this bar. Always nice to be where everybody knows your name, I think.

Why quality has always been essential for NAG’s own internal process (and a holiday).

2010-08-26T14:50:00.003+01:00

Quality is one of the things that matters a lot to our users.

That is certainly true, but another reality is that the NAG Library could not exist at all without the checks and tests that result in the best quality software. The NAG Library is particularly complex and involves many different people doing very different types of task and activities. The Library is created from intricate actions including:

- Investigating lots of numerical problems and potential approaches and solutions
- Designing hundreds of methods for the various chapters
- Producing example code, data and plots
- Creating installer and user notes
- Writing the detailed documentation and descriptions of routines
- Building the multiple different processor and compiler implementations of the Library
- etc…

The best approach to quality is essential, in making the NAG Library, precisely because it is this multifaceted combination of numerical mathematics, software engineering and evolved process. If it were not for the quality checks, at all points through the product organization, it would be impossible to get all the bits to mesh and work together.
NAG learnt this fact many, many years ago!

NAG quality methods are continuously scrutinised because of this internal need, as much as for users. Levels of peer review of documentation are increased, as new expert knowledge come to the fore. The frequency of automated implementation testing carried out during the development process is escalated, as the availability of test platforms increases. And so on.

The results are self-evident. The quality and reliability of the Library never falters, even as new content is added and further software environments are supported. The NAG quality process benefits all parties.
On a much lighter note I’ve just returned from our family holiday. I, my wife and our three young children flew to Athens, stayed for 3 days to revise the history, and then moved, about 100 nautical miles by ferry, to an island to explore and relax on beautiful beaches...

…in summary the vacation was made up of different sections, and was enjoyed by all five of us, because of the quality and detail that went into the planning and booking. Perhaps it was also made easier since we are all hardened travellers.

This leaves me with a question, both for people and organisations
– ‘Can we all manage just as well with a top quality process but only limited experience?’

I’d be interested to hear your view…

What's in a Name?

2010-08-20T11:07:00.000+01:00

Newcomers to NAG often ask about the names of NAG routines. Certainly the names appear strange and unfamiliar at first, but to many NAG users the name encapsulates the heritage of quality and precision offered by NAG.

Originally the NAG Libraries were created by contributors from U.K. academia. Some of these wrote in Algol 60 and others in Fortran. It was agreed to have two libraries of identical content, so that every contribution in one language had a counterpart in the other. Thus both communities were appeased. Early printed documentation included both language versions.

It was clearly desirable to be able to distinguish one language version of an algorithm from the same algorithm implemented in the alternate language and in general this is still desirable if a single program needs to call both instances of the algorithm.

Equally the early developers saw that, as new algorithms were developed, replacement routines would be required with perhaps different calling sequences, so that any naming scheme chosen needed to accommodate the ongoing development of the libraries. Another desirable feature would be to naturally group by name routines that lay in the same field of numerical computation.

With these goals in mind the founders of NAG considered two lines of approach.

The first idea looked at the concept of 'meaningful' names, names which might convey the algorithm being used or the problem addressed. Even today this is difficult, but in the early 1970s standard Fortran did not permit more than 6 characters for any routine name. And NAG certainly did not want to have Fortran and Algol names so radically different from each other that the algorithmic connection was obscured. The problem of updating also needed to be solved.

The alternative looked much more appealing. At that time a modified SHARE classification system had been published and NAG decided to base its routine names on that. This gave rise to the names we now see. Each name is 6 characters long, the final letter denoting the language of implementation. In those days 'F' denoted Fortran and 'A' denoted Algol 60. Subsequently this would be extended as NAG produced both single and double precision libraries in Fortran and as NAG produced other language versions such as C. The first characters, generally the first 3, indicated the broad subject ( or chapter) area to which the algorithm belonged. (Characters in positions 2 and 3 are always digits.) Thus E04WDF is a routine in the FORTRAN Library that lies in the E04 chapter. This chapter deals with 'Minimizing or Maximizing a Function'. The letters 'WD' uniquely identify the algorithm.

Of course modern Fortran no longer has a restriction of only 6 characters for routine names and so it is perhaps worthwhile to consider whether the NAG naming scheme can usefully take advantage of this relaxation.

The debate rages within NAG. Certainly the old problems remain if we wish to have 'meaningful' names, since routines will change with time. Do we describe the problem to be tackled or the algorithm used? If the former then NAG may wish to offer several techniques. In the C06 chapter we offer both Weeks' and Crump's method for Inverse Laplace transforms for example. Problem addressed is evidently not enough to uniquely identify an algorithm. Description of the routine by algorithm alone falls down heavily in the optimisation chapter where 'Sequential Quadratic Programming' and 'Quasi-Newton' and 'Modified-Newton' account for the vast majority of routines. We are thus led towards a combination of the two.

Looking in detail at the effect of combining problem and method into a meaningful name we rapidly encounter the very powerful and flexible routines of the NAG optimisation chapter. Here one routine often addresses several related, but different, problems and has both sparse and dense counterparts. Long names might indeed become very long if they are to convey any useful meaning and not deceive.

We see an early attempt by NAG to provide long, meaningful, names in its C Library. Users will be familiar with nag_opt_lin_lsq to denote the C routine e04ncc. This routine also solves convex qp problems, something the reader would not deduce from this name, yet this is probably its most important role.

Proponents of meaningful long names argue that the solution to this is to have separate routines with only a single purpose. Their opponents argue that it is easy for users to effectively choose their own names for routines, given the variety of options open to them: pre-processors, wrappers and aliasing for example.

What do you think? How many characters in a long name would you tolerate and how would you name NAG routines?

A day in the life of...

2010-08-18T09:48:00.005+01:00

... a Numerical Analyst:

Wow, this could be fairly boring, but here goes...

I arrive at my desk, slightly dishevelled from my daily cycle to work, whip my laptop out of the draw and unhibernate my old faithful friend. He and I have had some ups and downs over the past year but I am happy to say that we are still on good terms.

Having gone through the usual email checking routine, I ssh (remotely logon) onto one of the local linux work servers, to pick up from where I finished the day before.

And then it's on to slogging at the current routine, improving it's accuracy or speed, checking it works properly, documenting the code, screaming (under my breath) in frustration at the latest problem, or quickly writing a blog post to give my mind a break from the overflows and NaN's (think they are grandmas? Think again) of the numererical analyst's life.

But it's not all slog. One of the highlights for me, is in the cracking of problems. The simple answers, and the not so simple answers that when found, bring a smile of victory to the lips. Another highlight is in automating repetetive tasks, and creating/finding useful tools to speed up the work cycle. As I read in a magazine a while back, that if you do something a few times, it's worth making a quick tool to do the job for you.

At the end of the day, it's back home to my lovely wife, and to feed and bath my young son, only a Numerical Analyst in the back of my mind.

Hey, you sass that hoopy Fortran?

2010-07-30T09:25:00.000+01:00

We have to work with lots of different Fortran compilers at NAG. So far at Mark 22 of our Fortran Library we have built using products from six distinct sources (GNU, IBM, Intel®, NAG, PGI and Sun). The list of what we build with naturally changes from Mark to Mark as different vendors fall by the wayside or as commercial considerations preclude us from creating a particular Library implementation.

With the six vendors listed above we are blessed that the base Fortran coverage is Fortran 95. At previous Marks we still supported compilers covering only Fortran 77 (g77 and pgf77 for example). This forced us to jump through all manner of hoops so that we could do some of the really useful stuff from Fortran 90—like, woah, dynamically allocating memory—in a Fortran 77 way.

Some hoops still exist though, even with Fortran 90 code.

For example, we recently ran into the following:

    PROGRAM str_read
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: wp = KIND(0.0D0)
!      .. Local Scalars ..
       REAL (kind=wp)                  :: rval
       INTEGER                         :: ioerr
       CHARACTER (200)                 :: my_str
!      .. Executable Statements ..
       my_str = '1.0D400'
       READ (my_str,'(E16.0)',iostat=ioerr) rval
       PRINT *, 'IOERR = ', ioerr
    END PROGRAM str_read

Ignoring the fact that the code fragment is a little contextless, just imagine that we're trying to trap overflowing real values being read from a CHARACTER. On five of the six compilers given above, two return a nonzero ioerr from the READ, two carry on happily and set rval to Infinity, and one core dumps! (At the time of writing, I didn't have access to all six compilers.)

Conferring with the Fortran standard (e.g., the draft Fortran 95) we see the dreaded phrase The set of input/output error conditions is processor dependent. (my italics), so we can't even complain to the compiler vendors about this behaviour! We have good relationships with the vendors, so when our building process reveals a compiler bug we'll report it. But as you can imagine, sometimes it can be a battle to achieve what we want across many Fortran platforms in as simple and maintainable a way as possible.

The next level of complication comes from ensuring that what we do in Fortran is callable from non-Fortran environments. We usually take Microsoft Excel as a sufficiently-far-removed example. I won't go into details of the usual issues raised by cross-language programming; however, now we're at a base level of Fortran 95 we've been looking at pushing our envelope more and trying out some Fortran 2003 features, specifically C Interoperability.

M'colleague Nicolas has been working on a suite of image processing routines. In pure Fortran a NAG_IMAGE TYPE exists to hold the pixel values and other metadata for the image. To make this passable from outside Fortran we're using Fortran 2003. It's understandable with new features in the language that there should be a settling-in period, but it's no exaggeration to say that every compiler we tried did, at some point, fall over on our new code. Thanks to the hard work of the compiler developers on the various teams we're now at a stable state, but there are still some vendors who are lagging behind with Fortran 2003 and for which we have to exclude the image-processing code from our testing.

Looking forward (!) to Fortran 2008, I'm interested to see how the new special-function INTRINSICs (BESSEL_J*, BESSEL_Y*, ERF* and *GAMMA) will behave across different compilers. In theory these new INTRINSICs mean we'd be able to withdraw a sizeable chunk of our S Chapter. But I wonder how easy that will be?

(The images in this post are from now voyager, but are probably copyright PolyGram Filmed Entertainment.)

A Life Well-Lived: Erwin Ruppenthal (1958 – 2010)

2010-07-23T04:37:00.005+01:00

There’s probably some unwritten rule that blogs are off-limits to memorials. If so, I’m happy to break it for our colleague Erwin (or Erv as he was known at home). He died early this week after a long struggle with brain cancer, at home, surrounded by his wife Heidi and sons Alex and Kurt. If we take a closer look at his life we can learn a lot about the value of work and about having a purpose in life greater than earning a paycheck.

Erwin taught me a number of lessons about work and I’ll tell you about some of them so that things he valued might live on in each of us, but first let me tell you a little about him. He came to NAG in early August 1990 to work in IT at our US office. What I’m told was that he was an earnest, hardworking German citizen in his early 30’s who most would describe as shy and self-effacing. It was nine years later that I encountered him on my first day as his boss.

He was a man of two countries and two cultures who embraced both. Following sports was a passion: American baseball, American football and European football featured prominently. In his office you could find emblems of the Chicago Bears and Chicago Cubs along with a Bayern Munich FC flag. I read in his obituary that he almost never missed a Naperville Central high school football game over many years. If you knew him, none of this loyalty and dedication would surprise you because he brought the same qualities to work every day,

I’ve worked with few people who could equal his dedication to his work, our mission as a company and the customers we serve. His loyalty to serving our customers showed in the countless problems he solved late at night and early in the morning from home while monitoring our helpdesk inbox or remote web site. He had no hesitation in coming to me if he felt we weren’t being fair to a customer or living up to our promises. He confirmed in me that our first principle had to be to do the right thing by a customer and trust that they would return the same for us.

He taught me some powerful management lessons though I doubt he ever spent a day in business school. Early on, I would see something we were doing that needed change, blindingly obvious (at least to me). Erwin would sit patiently as I explained, asking questions and listening. At first I thought he was just humoring his boss so I’d go back a week later and repeat the process and come away with the same impression. At some point later I would notice that a change had been made with no fanfare. Erwin would take my idea, however half-baked, improve it and implement it. Therein I learned the value of patience, trust and that overused cliché, empowerment. With Erwin, I mostly needed to plant a seed and get out of the way.

Erwin didn’t spend money frivolously in his personal life, as far as I could tell, but when he did advocate spending it at work, he pursued quality, durability and things which made for better customer service. We have (still) a truly ancient UNIX server that serves as an archive of customer data. Through two office moves and occasional disk failures Erwin kept it going reliably. His care for such things could be summed up by a short story of the first of those moves. The move was from one building to another in an office park in suburban Chicago separated by perhaps two football fields of lawn and street. It was late on a cold, moonlit Friday night in early January 2002. The movers had finished earlier in the evening and all of the staff had gone home. All that was left to do was to power down the servers, move them across the street and re-start them. They were already ancient and Erwin wasn’t about to trust them to the movers so he and I methodically shut everything down and loaded it all on a 4-wheel cart. A few minutes later we were pushing the cart gently down the middle of the deserted street when I jokingly wondered what would happen if a police car came by. He had prepared for everything but that. Even in the small things, he planned carefully and followed through. That server is still running.

I could tell you more but these things tell us the enduring values Erwin brought to work: loyalty to people and a purpose bigger than himself, patience, persistence and reliability. He was more than a co-worker for everyone he touched. He was a true colleague and a trusted friend. Erwin may be gone but his memory and his values live on in each of us.

Time Machines and Supercomputers

2010-07-19T15:17:00.003+01:00

I found a Linpack App for the iPhone last week. Nothing special, just a bit of five minute fun. It seems a 3G model achieves about 20 MFLOPS. [Note 1]

What's that got to do with time machines? Well it got me thinking "I wonder when 20 MFLOPS was the performance of a leading edge supercomputer?" Actually, it was before the start of the Top500 list (1993), so finding out was beyond the research I was prepared to do for this blog.

So I thought instead about the first supercomputer I used in anger. As soon as I name it, if anyone is still reading this waffle, you will immediately fall into two camps - those who think I'm too young to be nostalgic about old supercomputers yet - and those who think I'm too old to be talking about modern supercomoputers :-).

It was a Cray T3D.

You're still waiting for the time machine bit ... hang on in there.

My application on that T3D sustained about 25 GFLOPS. Which is about the same as a high end PC of recent years. What this means to me is that anyone who cares to apply the effort today with a high end PC, could get comparable results to that work of 15-20 years ago that needed the supercomputer.

Or, in other words, that supercomputer gave us a 15-20 years time advantage over everyone who didn't have supercomputers - or a few years over others with smaller supercomputers. [Note 2]

That is one of the key benefits of High Performance Computing - the ability to get a result before a competitor - you could say HPC is a time machine for simulation and modelling.

Now for the [Notes] - which actually contain the real story!

Note 1 : It's not really true to say the iPhone 3G can do 20 MFLOPs - all we can say is that particular App achieved 20 MFLOPs on that iPhone 3G. The result is a factor of both the software and the hardware. Better performance can come from optimising the application as much as from buying a more powerful phone.

Note 2 : If fact, even with the same supercomputer, it would be hard for most people to replicate the results - simply because there was as much value in the software (physics, algorithms, performance engineering, implementation, etc) and the associated validation and verification program as there was in the supercomputer.

The supercomputer offered us a time machine. But the attention to performance and scalability in the application enabled us to actually use that time machine to get results faster than others - even if those others used the same supercomputer. And the validation and verification effort meant that we could trust what our time machine was telling us.

Fantasy Football – a classic Portfolio Optimisation problem

2010-07-08T09:24:00.007+01:00

England out the World Cup, German colleagues, customers, collaborators (actually not just Germans, Americans, Scots,…) and so called friends all sending me mocking e-mails and texts.

“I hear OXO are making a new product. The packaging is white with a red cross and they're calling it the laughing stock.”

How could I channel my frustration? NAG Blog to the rescue.

In my last blog I promised to reflect on my early career at NAG which isn’t that long ago compared with some of colleagues.

My commute to the NAG office is longer than I’d like, but there is a positive side…

· Valuable thinking time

· Audio books - I recently enjoyed “No One Would Listen” by Harry Markopolos. The exclusive story of the Harry Markopolos–lead investigation into Bernie Madoff and his $65 billion Ponzi scheme.

Back to “thinking time,” on one commute into Oxford I’d been pondering the previous day’s internal training course. I was recalling my senior technical colleague’s wise words. That wise colleague was David (some of you will have met him). He had given an internal training course to the Sales and Marketing group. His presentation was a simple introduction to the topic of Optimisation and he had touched on modern portfolio theory even explaining the efficient frontier.

Some of you may be lucky enough to have heard this talk before….

Imagine an English ice cream manufacturer. He might run a very profitable business, but recognise that his business is very dependant on long sunny, dry summers. Wishing to diversify and minimise his risk he chooses to start another business thus protecting himself against a cold, rainy summer with an umbrella company.

Subsequently he went on to talk about stocks and shares and the benefits of diversification. This struck a chord with me. I learnt early on in my adult life the importance of diversification. Remember GEC that became Marconi. Where did Marconi go wrong? As a result of shares options I had a portfolio dominated by Marconi and consequently suffered a very large paper loss. Oh, why didn’t I meet David earlier? He went on to speak about how diversification can be achieved by holding a selection of stocks from a spread of geographic regions, industry sectors etc.

Of course this whole training course was aimed at helping us understand NAG’s Optimisation Routines . For those of you wishing to learn more about NAG’s optimisation routines you should refer to the NAG Library Manual and the appropriate chapter introduction.

Anyway, as I drove back into work that morning I had that light bulb moment. Fantasy football is like portfolio optimisation! Well, you can imagine how I sprinted into the office from the work car park eager to share my discovery with David.

Me: “Remember your lecture about optimisation and portfolio optimisation? Fantasy football is like this, isn’t it?”

David: “What is Fantasy Football?”

Me: “Well

·You have a set of rules.

· You have a limited pot of money that you can spend to pick a squad of 15 football players.

·You then have to pick X Goal Keepers, Y defenders, …

· You can only pick Z players from any one team

· Players win points by

o Keeping a clean sheet

o Scoring a goal

oAssisting a goal

· Players lose points by

oConceding goals
o Getting red or yellow cards.”

Pause for thought.

Me: “These are constraints aren’t they?”

Patronising nod from David

Me: “And this is an Integer problem, isn’t it?”

Laughter. David: “Well, at least you were listening last week. Now, go and get on with your work.”

Me: “No, you don’t get it. You’re going to help me program this up and we’ll see how the NAG Optimisers perform! I bet they won’t beat me. I’ve won the office fantasy football competition two years in a row. I would be interested to see how it performs though.”

Two or three days pass. I find a way of getting all the historic data from “The Official Fantasy Game Of The Premier League” into an Excel spreadsheet.

A week later I manage to persuade David to code up one of NAG’s optimisation routines to optimise my portfolio of football players. What does this mean? Well one is supplied with a list of players. Each player has a value and their previous season’s points total is listed. So one might choose a strategy which optimises a squad of 15 players based on a maximum spend of £100 million. A more risky approach might be to assume you will always have 11 fit players so you pick the 4 cheapest players and then look to maximise your return from an optimal 11.

Fortunately NAG supplies Visual Basic Declaration Statements and C Types with its DLLs. This makes these libraries easy to use for VB programmers who could easily code up one of NAG’s Optimisers to pick a Fantasy Football side. Yes, NAG’s Libraries are easy to link to Excel and we include a simple Portfolio Optimisation example on our website.

Let me explain this screen shot.

The two teams that were picked entirely by NAG’s Solvers were David’s “Too Hot” and Sven’s “Old NAG’s Best 15.”

· “Too Hot” was based on picking the four cheapest players (with the highest return) and then the “Optimised 11”

· “Old NAG’s Best 15” was based on optimising the entire 15 i.e. taking no account that only 11 players can be picked and used in each match.

In my next blog I’ll share with you how “Too Hot” and “Old NAG’s Best 15” finished the season.

The power within graphics cards

2010-07-02T08:44:00.012+01:00

NAG is a company that likes to be present at scientific events all around the world. One of the main reasons for this is we get face-to-face time with NAG users and this gives us the chance to listen to their needs and reiterate that 40 years on, NAG is still a not-for-profit company for which collaboration is essential.

NAG at the GPGPU Developer Briefing in New York City

Last week I had the pleasure of representing NAG at two events in North America: GPGPU Developer Briefing in New York City and the 6th World Congress of Bachelier Finance Society in Toronto. While manning the NAG stand together with my colleagues we had many interesting conversations with those attending and made them aware of one of NAG’s most recent developments – the NAG Library for GPUs.

6th World Congress of Bachelier Finance Society in Toronto

We expected lots of interest in GPUs at the briefing in NYC, and were all happy to describe our successes with Monte Carlo simulations on GPUs. Currently the NAG Library for GPUs consists of two random number generators: pseudorandom L’Ecuyer MRG and quasirandom Sobol, a generator of sample paths for a Brownian bridge, and statistical distributions. Other features, such as PDE methods or optimization will follow soon. It was really good to see the high level of interest in the Library and how to call it from Excel and other environments.

The Bachelier Society Congress in Toronto attracted an audience primarily from financial academia. We spoke to various NAG users who talked about their use of NAG routines for optimization, wavelet transforms, and RNG to name a few, but this audience also expressed keen interest in our GPU work, in particular the work with a major bank where the NAG Library for GPUs performed a Monte Carlo simulation over 200 faster than on a single-core CPU! Of course, this level of speed up cannot always be achieved, but it is an indicator of the power of GPUs combined with the robust, efficient and stable implementation of random number generators in the NAG Library for GPUs.

My first experiences calling the NAG Library for SMP and Multicore

2010-06-23T14:00:00.010+01:00

A quick introduction: I work for NAG as a Technical Support Engineer in Germany. I look after clients in Germany, Switzerland and Austria and am interested in approximation and interpolation methods as well as .NET languages. Naturally, I’m big fan of the NAG library for .NET as well as the NAG Library for SMP and Multicore.

As a non FORTRAN programmer to be able to run some tests with the new NAG Library for SMP and Multicore I decided to choose a few simple tests that I could run using both CPUs of my laptop. I created some .NET console applications, wrappers around the FORTRAN routines that were able to pass all the arguments to the routines in a .NET fashion and implemented a time counting function. There are many parts of the Library parallelized, including quadrature, partial differential equations, interpolation, curve and surface fitting, linear algebra, correlation and regression analysis, multivariate methods, random number generators, time series analysis, sorting and special functions. I chose routines from the correlation, curve and surface fitting and random number generation chapters.

The results were calculated on a Intel® Core™2 Duo Prozessor P8700 (2.53GHz,1066MHz,3MB) machine with Windows Vista as operating system and Microsoft Visual Studio 2008 (32‐bit project) as compiler. Each function was used to solve a large enough problem to allow for parallelism. Each test case was run with one or two threads. The number of threads was set using the system variable ‘OMP_NUM_THREADS’ in the Visual Studio console window.

A ‘QueryPerformanceCounter’ function, was implemented in all the examples to calculate the time which the routines needed for calculating the result. Below is the class implemented for each of the examples which can be used by first creating the constructor HiPerfTimer and then using the methods Start, Stop and Duration:

HiPerfTimer pt = new HiPerfTimer();
pt.Start();
//routine
pt.Stop();
pt.Duration();

internal class HiPerfTimer
    {
        [DllImport("Kernel32.dll")]
        private static extern bool QueryPerformanceCounter(
            out long lpPerformanceCount);

        [DllImport("Kernel32.dll")]
        private static extern bool QueryPerformanceFrequency(
            out long lpFrequency);

        private long startTime, stopTime;
        private long freq;

        // Constructor

        public HiPerfTimer()
        {
            startTime = 0;
            stopTime = 0;

            if (QueryPerformanceFrequency(out freq) == false)
            {
                throw new Win32Exception();
            }
        }


        // Start the timer

        public void Start()
        {

            Thread.Sleep(0);

            QueryPerformanceCounter(out startTime);
        }


        // Stop the timer

        public void Stop()
        {
            QueryPerformanceCounter(out stopTime);
        }


        // Returns the duration of the timer (in seconds)

        public double Duration
        {
            get
            {
                return (double)(stopTime - startTime) / (double)freq;
            }
        }
    }

After compiling the examples with the C# compiler option ‘csc’, the number of threads is first set to the value 1 and then set to the value 2. The only result which is printed is the time which the routine needs for calculating the result.

Correlation

NAG routine C06PKF calculates the circular convolution or correlation of two complex vectors of period n.
In this example to complex vectors are build up with a period of n (n=5000000) and the correlation and the circular convolution are calculated.

C: \C06PKF\C06PKF>csc c06pkf.cs

C:\ C06PKF\C06PKF> set OMP_NUM_THREADS=1
C:\C06PKF\C06PKF>c06pkf
Duration: 5.74 sec

C: \C06PKF\C06PKF>set OMP_NUM_THREADS=2
C: \C06PKF\C06PKF>c06pkf
Duration: 2.98 sec

Random Number Generation

NAG routine G05YKF generates a quasi-random sequence from a log-normal distribution. It must be preceded by a call to one of the initialization routines G05YLF or G05YNF. The number N (N = 10000000) of quasi-random numbers of dimension IDIM (IDIM = 4) are passed and also the mean (IDIM) of the underlying Normal distribution for each dimension and the std (IDIM) standard deviation of the underlying Normal distribution.

C: \ G05YKF\G05YKF>csc c06pkf.cs

C: \G05YKF\G05YKF>set OMP_NUM_THREADS=1
C: \G05YKF\G05YKF>g05ykf
Duration: 3.52 sec

C: \G05YKF\G05YKF>set OMP_NUM_THREADS=2
C:\ G05YKF\G05YKF>g05ykf
Duration: 2.06 sec

Approximation

NAG Routine E02CAF forms an approximation to the weighted, least-squares Chebyshev series surface fit to data arbitrarily distributed on lines parallel to one independent coordinate axis. It determines a bivariate polynomial approximation of degree k in x and l in y to the set of data points , with weights , for and . That is, the data points are on lines , but the x values may be different on each line. The polynomial is represented in double Chebyshev series form.
N, the number of lines on which data points are given, was set to 1000
K, the required degree of x was set to 100 and
L, the required degree of y was set also to 100.
Also the fitting polynomial was evaluated at the data points using the routine E02CBF.

C:\E02CAF\E02CAF >csc e02caf.cs

C:\ E02CAF\E02CAF >set OMP_NUM_THREADS=1
C:\ E02CAF\E02CAF >e02caf
Duration: 4.77 sec

C:\ E02CAF\E02CAF >set OMP_NUM_THREADS=2
C:\ E02CAF\E02CAF >e02caf
Duration: 2.48 sec

Fig 1: Timings for the SMP enabled NAG routines with one or two threads
My experience with presenting these examples at different clients has been very positive, but also interesting! “Warum verwenden Sie diese Bibliothek in .NET? Das ergibt doch keinen Sinn, da wir nur in Fortran programmieren!” “Why are you using the Library in .NET and not directly in Fortran? This doesn’t make any sense, we are coding only in Fortran!”
Others were excited to learn and understand they might use their second processor!
My conclusions are as follows

Never show a C# program to HPC Computer Programmers who love Fortran and think languages such as C are radical!
NAG Library for SMP and Multicore really is easy to use and for those who are current NAG Fortran Library users the transition really is painless. Users just need to take care of the number of threads used.

I should admit a bias as a NAG employee, but perhaps more importantly I’d like to use this blog to talk to users / prospective users and offer help. As I said in my introduction I am interesting in this, but also hearing from those who have interests in approximation and interpolation. I plan to work with some of my developer colleagues to expand some of NAG’s interpolation solvers.

Technical computing futures part 2: GPU and manycore success

2010-06-22T17:31:00.001+01:00

In my previous blog, I suggested that the HPC revolution towards GPUs (or similar many-core technologies) as the primary processor has a lot in common with the move from RISC to commodity x86 processors a few years ago. A new technology appears to offer cheaper (or better) performance than the incumbent, for some porting and tuning pain. Of course, I’m not the first HPC blogger to have made this observation, but I hope to follow it a little further.

In particular, my previous blog suggested the outcome might be: “at first the uptake is tentative ... but in a few years time, we might well look back with nostalgia to when GPU’s were not the dominant processor for HPC systems” – in other words, hard going initially, but GPU/many-core will “win” eventually. I even ended up with an ambitious promise for my next blog (i.e. this one): “an idea of what/who will emerge as the dominant solution ...”

Continuing the basis of using the past to guess the future, my prediction is that the next steady state of HPC processors will be GPU-like/manycore technologies (for most of the FLOPS at least) and, just like the current steady state (x86), those few companies with the strongest financial muscle will eventually own the dominant market share. However, other companies will have pioneered many of the technologies that make that dominant market share possible, enjoying good market share surges in the process.

I can even have a go at predicting some of the path that might get us to the next steady state of HPC architecture. NVIDIA has already shown us that GPUs for HPC are sometimes a good solution – and importantly, that a good programming ecosystem (CUDA) really helps adoption. Over the last year or so, I’d say the HPC community has moved from “if GPUs can work in this case ...” to “how do I make GPUs work across my workload?”

As Intel’s Knights processors bring us many-core but with a familiar x86 instruction set, we might learn that getting good performance across a broad range of applications is possible, but critically dependent on software tools and hard work by skilled parallel programmers. AMD’s Fusion with tighter links between CPU & GPU, could show that the nature of the integration between the many-core/GPU unit and the rest of the system (be it CPU, network, main memory etc) will affect not only maximum performance on specific applications, but maybe more importantly the ease of getting “good enough” performance across a range of applications.

I don't know of any GPU/many-core/accelerator announcements from IBM, but it’s always possible IBM will throw in another useful contribution before the dust settles. They were one of the first into many-core processors for HPC acceleration with Cell and they cannot be easily counted out of top end HPC solutions - e.g. the forthcoming Blue Waters (POWER7) and Sequoia (BG/Q) chart-toppers.
But back to my “winner” prediction. When the revolution settles into a new steady state of mostly GPU/many-core for HPC processors, there won’t be (can’t be) critical distinctions between the various products anymore for most applications. Whichever product we consider (whether GPU or x86-based or whatever), many-core is sufficiently different from few-core (e.g. 1-8 cores) to mean that the early winners have been those users who are easily able to move their key applications across to get step changes in cost and performance.

The big winners in the next stages of the GPU/manycore emergence will be those users who can move the bulk of their high-value-generating HPC usage to many-core processors with the most attractive transition (economy and speed) compared to their competitors.

So what about the dominant solution I promised? For the technology to be pervasive, first there must be greater commonality between offerings (I stop short of standardization) so that programmers have at least a hope of portability. Secondly, users need to be able to extract the available performance. Ideally these would mean a software method that makes many-core programming “good enough easily enough” is discovered – and if so, that software method will be the dominant solution, across all hardware.

Or, if the magic bullet is still not market ready, skilled parallel programmers will be the dominant solution for achieving competitive performance and cost benefits - just like it is for HPC using commodity x86 processors today.

The meaning of recognition

2010-06-16T13:41:00.003+01:00

A few years ago, NAG decided to brush up its public image and issued its staff with The Company Shirt (actually, it turned out that there was so much money in the marketing budget that each of us was able to have our own shirt, as opposed to being obliged to take it in turns to sport a single item of clothing). It's a rather splendid garment (you can see my colleague Mike Dewar elegantly modelling his below) that proudly but discreetly displays the company logo and, for good measure, the full name of the company (lest the abbreviation be misinterpreted as an exhortation to complain endlessly). The idea - which, I'd imagine, is common to just about every organization in the world - is that the shirt can be worn on exhibition stands, when giving commercial presentations or making customer visits so that a (somewhat loosely) unified image of the company is presented to the outside world. Members of staff have acceded to this idea with varying degrees of alacrity; speaking for myself - following an initial period of uneasiness where I suspected (quite without foundation) that the next step on the road to an improved image would be The Company Song - I've been happy to wear The Shirt on every appropriate occasion.

And on some less appropriate ones as well. For example, last night I was taking part in a choir rehearsal as part of the preparation for a well-known religious leader's visit to the UK later this year. Not having had time to get changed after work, I was wearing The Shirt as we collectively negotiated the joys of counting bars, leaping fourths, subdividing triplets and other more applied forms of numerical analysis. Approaching the conductor - who'd been brought in from another parish in order to adeptly marshal our enthusiastic but slightly unfocussed efforts - with a technical question at the end of the rehearsal, I was a little surprised when he asked if I worked for NAG. Wondering if he was about to quiz me about - say - our optimization routines, I replied - somewhat cautiously - in the affirmative. "Great stuff," he responded. "I used the NAG Library all the time at university when I was programming in Fortran - it was really, really good." Owing to the context, the generous and unlooked-for compliment was so surprising that I forgot to say that these days, the Library wasn't only available to Fortran programmers (on reflection, perhaps that was just as well on this occasion) but I was also reminded that this kind of encounter isn't at all uncommon. Given the remarkable age of the company, perhaps it's only to be expected that you frequently bump into users - or ex-users - of your products, but it's still gratifying when they're able to share positive experiences - or happy memories - of it. Maybe I should start work on that Company Song after all.

Revealing the future of technical computing: part 1

2010-06-08T16:47:00.011+01:00

I recall some years ago porting an application code I worked with, which was developed and used almost exclusively on a high end supercomputer, to my PC. Naively (I was young), I was shocked to find that, per-processor, the code ran (much) faster on my PC than on the supercomputer. With very little optimization effort.

How could this be – this desktop machine costing only a few hundred pounds was matching the performance of a four processor HPC node costing many times that? Since I was also starting to get involved in HPC procurements, I naturally asked why we spend millions on special supercomputers, when for a twentieth of the price, we’d get the same throughput from a bunch of high-spec PCs?

The answer then (and now) was that I was extrapolating from only one application, and that application could be run as lots of separate test cases with no reduction in capability (i.e. we didn’t need large memory etc, just lots of parameter space). However, the other major workload (which I also ported and also ran fast on the PC) would not have been able to do the size of problem we wanted on a PC – we needed the larger memory and extra grunt from parallel processing. (We did look at the newfangled Network Of Workstations emerging at the time but decided it might be a wolf in sheep’s clothing. Sorry.)

In the end, we had to find a balance between (a) speed at lowest cost for the one application; (b) the best capability for the other application (i.e. fastest solution time for the largest problems); (c) ease of programming – to get a good enough (fast-enough) code developed with the limited developer effort and funding we had; and (d) whole life affordability.

Why do I foist this reminiscence on you? Because the current GPU crisis (maybe “crisis” is a bit strong – "PR storm" perhaps?) looks very much the same to me. The desktop HPC surprise of my youth has evolved into the dominant HPC processor and so for some years now, we have been developing and running our applications on clusters of general purpose processors – and a new upstart is trying to muscle in with the same tactic – “look how fast and how cheap” – the GPU (or similar technologies – e.g. Larrabee, sorry Knights-thingy).

The issues are the same: (a) for some applications, GPUs offer substantial performance improvements for considerably less cost than a “normal” HPC processor; (b) for other applications, the limits such as off-card bandwidth etc mean that GPU’s cannot deliver the required capability; (c) the underlying concern is ease of programming for GPUs; (d) affordability – sure GPU’s are cheap to buy, but what about power costs when in bulk, or code porting costs, etc?

Maybe the result will be the same as when commodity processors and clusters eventually exploded to leave custom supercomputer hardware as the minority solution. At first the uptake (now) is tentative - and painful. Some will have great success stories, many will get burnt. But in a few years time, we might well look back with nostalgia to when GPU’s were not the dominant processor for HPC systems.

I’ll continue on the future of HPC in my next blog in a few days, including an idea of what/who will emerge as the dominant solution ...

Why is writing good numerical software so hard?

2010-06-07T17:01:00.011+01:00

People sometimes say to us "Why does NAG continue to exist? Surely all the good numerical algorithms have already been devised and implemented?". Well, this question ignores the fact that people keep finding new kinds of problem to solve. It also forgets that we need to make our current software continue to work robustly and efficiently as new hardware and software infrastructures develop - we need our code to run on all the platforms and in all the environments that our customers use.

Making the NAG library work with a new compiler is always a non-trivial task (and sometimes it can seem more like a nightmare - just ask some of my colleagues!) It may be because we've got some code that's not yet been tested in a wide enough variety of environments. Or it may be because compiler optimization switches can reorder our code in ways that make the floating-point arithmetic behave in a way that we hadn't anticipated (which could be the compiler's fault or could be our fault, depending). But whoever is to blame, we have to fix it.

When you are dealing with complicated numerical codes, tracking down a problem might take a lot of hard work with a debugger (or with print statements if you're unlucky enough to be in an environment where a debugger just won't work). And, what look like even the simplest operations might turn out to be a lot more complicated than you'd expect.

As an example, take the method of dividing one complex number by another. Languages like Fortran have in-built support for complex numbers and complex arithmetic, but not all do, and so the NAG Library does contain some routines to help.

Performing a complex division is simple in principle. If the numbers are called X and Y, in general they are composed of real and imaginary parts, let's call them X.re and X.im, Y.re and Y.im. Then to compute the value Z = X / Y, just multiply both the top and bottom halves of the quotient by the complex conjugate of Y (i.e. the complex number you get when you negate the imaginary part of Y). Since we're multiplying the top and bottom by the same thing, there's no net effect on the answer, but the clever part is that when you multiply a complex number by its own conjugate, the answer is real. That means that Y * conjugate(Y) is real, and the complex division operation is reduced to an (easier) complex multiplication followed by a few real divisions. The complex multiplication is done by cross multiplying (in real arithmetic) the real and imaginary parts of X and Y, then adding and subtracting appropriate pieces to form the real and imaginary parts of the result. Counting all the floating-point operations up, this turns out to be 6 multiplies, 3 additions or subtractions, and 2 divides - a total of 11 operations - quite significant compared to a single real division, but still not a lot of work.

So - it can't be all that hard to get it right can it?

Wrong. We need to take some care here. For example, if the complex numbers X and Y are very large (and I'm talking about larger than you might ordinarily meet - say about the size of or bigger than the square root of the largest number you have on your computer), there's a chance that although the number Z = X / Y is a perfectly reasonable number, some of those intermediate 11 computations could cause arithmetic overflow. You'll get garbage as a result - not good. People are not happy if they divide a large number by itself and get something not equal to 1!

The implementation of complex division that had been in the NAG Library since the early 1970s took care to avoid unnecessary overflow and underflow by doing some cunning scaling and rearranging into 4 multiplies, 4 adds or subtracts, and 3 divides. The algorithm is documented here in the Mark 22 NAG Library Manual

nagdoc_fl22/pdf/A02/a02acf.pdf

and I'd never given it much thought - I always assumed this algorithm was fine.

Well, we recently discovered that the method used by routine A02AC was not as perfect as I'd thought! The reason? Old age, really. The algorithm was based on a method by Robert Smith, published in Communications of the ACM way back in August 1962, and this had been developed in the dark days before the IEEE standard for floating-point arithmetic was designed in the 1980s (Standard for Binary Floating-Point Arithmetic (IEEE 754-1985)).

The IEEE standard introduced the concept of denormalized numbers - tiny but non-zero numbers which have less precision than regular numbers. It turned out that if you fed these into the NAG routine, you might still end up with garbage. For example, dividing the number X which has real and imaginary parts equal to the same denormalized number would return not the expected value 1, but a complex number composed of real part infinity and imaginary part NaN (NaN stands for "Not a Number"). Yuk.

Does it matter? Probably not much - not many people would call a NAG routine for something as "simple" as complex division (they'd rely on their compiler) and even if they did the chances are they would not be using denormalized numbers. But, a curse of working at NAG is that as soon as you find out about a problem like this you feel obliged to do something about it, and that's what my colleague Mat did. He located newer published methods which did take account of denormalized numbers, by introducing extra scaling operations, and based a revised version of the NAG routine on one of those. I won't bore you with the details of it here - suffice it to say that the newer method is a bit more complicated, but deals with the problem nicely, and it will debut in the next version of the NAG Library.

This simple case illustrates that even for the most basic numerical software things can go wrong if you don't take care. Imagine what it's like for something complicated like the quadratic formula! (I'm only half joking - safely finding the roots of a quadratic equation takes a frighteningly large number of lines of code).

To people who provide technical support

2010-06-03T13:43:00.012+01:00

After feedback from very satisfied NAG users, it seems to be high time to publicly appreciate the folks who work in front line technical support teams. To point out what great value they provide to customers, not just for NAG but across many software organisations worldwide. It is a very important job that involves much more than answering questions. It is about asking the right return questions and about careful communications.

The idea of a TV detective springs to mind; working the crime scene at an upmarket golf club near Hollywood, with a smart party in the background… image Detective Lt Colombo, in the clubhouse, picking his words with oh-so-much care, and asking questions in the most un-intimidating way. At the same time, he has four suspects lined up at the back of his mind and is waiting to discount each one as the facts of the case (or support problem) are drawn out into the light.

Who at NAG?
------------ What! Not a doctor?

The support job sounds easy to start with: When someone gets in touch with a support desk it is because they need help. But it needs special skills to uncover all the important details that provide the clues as to what is going wrong. The support role needs an ability to understand what the user is trying to achieve and the mechanisms he or she wants to use. It needs enough knowledge to be able to challenge the user’s ideas of how their own system or application should behave and it needs the balance of authority and diplomacy to work through the possible solutions.

A little like the down at heel TV detective investigating the crime at the exclusive golf club, it is common for the person handling a support request to have no pre-knowledge about the user’s application or environment. So it can be a steep learning curve in support. But the rapid self-education in the user’s field can be critical to make for clear technical conversations. This is particularly so when communicating with NAG Library users, since they work in very many different industries and subject specialties.

The final twist in the detective plot is often the most important. Here, unlike the TV detective who has to expose the evil perpetrator, the very best support people don’t really say that the user may have made any mistakes. Instead they leave the user with new and valuable knowledge and often a favoured contact that they are happy to turn to if they need real help again in the years to come.

To end I need to be clear that this is not meant to suggest that support people are dishevelled, or that they dress badly (I am in no position to comment these subjects). And support teams actually do have an in-depth knowledge of many fields - at NAG these are the same type of people who develop the NAG library and provide all the detail support information http://www.nag.co.uk/Forms/support.asp and the diverse technical papers http://www.nag.co.uk/doc/techrep/index.asp that help users make the best use of the NAG library every day.

The analogy is intended to show how, like good detective work, the job is skilled, subtle and is critical. Many organisations are built on great support teams who provide astoundingly good help. Just listen to the users of NAG software saying – 'keep up the good work'.

NAG at ISC 2010

2010-06-01T13:59:00.007+01:00

NAG is attending International Supercomputing 2010 in Hamburg, Europe's largest Supercomputing conference. Attendence is apparently at its highest ever and many of the sessions have been standing-room only. During Thomas Sterling's talk yesterday, on "Parallel Computing in the Years to Come", there were people sitting in the aisles and listening from outside in the corridor. The conference starts with the announcement of the latest Top 500 List, and the big news was that a Chinese system, Nebulae, was now number 2, behind Jaguar at Oak Ridge National Laboratory in the US. This is the second Chinese system to get into the top 10 and, like Tianhe-1, currently at number 7 in the list, it gets its performance from GPUs. In fact its theoretical peak performance is higher than Jaguar's, but it is hard to achieve that performance with existing software development tools. GPUs are very much the flavour of the moment. There are some interesting new products for GPU programming, in particular both Allinea and TotalView have versions of their parallel debuggers for CUDA code. The NAG booth has been very busy and it seems that at least half the people who drop by are interested in our own software for GPUs, currently available for beta testing. The feedbacks and suggestions for additional functionality will be a great help to us in bringing this product to market. The other interesting development in this area is that we're seeing increasing interest in software for OpenCL, largely due to the adoption of ATI cards. We'd love to hear from people who are using OpenCL and would like to make use of numerical components such as ours.

NAG Toolbox for MATLAB Documentation: Features in Development

2010-05-27T13:54:00.000+01:00

The documentation for the NAG Toolbox for MATLAB aims to integrate as far as possible with the host MATAB environment and provide the same integrated documentation experience as found with built in MATLAB functions.

MATLAB is itself in continuous development with a 6-month release cycle. New features appear in the MATLAB GUI, and as far as possible we try to use these features if it is possible to keep the documentation working on all the versions of MATLAB that we support.

The Function Browser

One relatively new feature, added at MATLAB release 2008b is the function browser. This is available in the help menu, or by pressing Shift-F1 (or equivalent keyboard shortcut, depending on platform and local customisation). It is essentially a modified table of contents to (just) the function reference pages, with built in search box and popup windows showing abbreviated summary documentation.

Currently the functions of the NAG Toolbox do not show up in this function browser, however this has now been enabled in the development builds and should, subject to final testing, be appearing in newer releases of the toolbox scheduled for later in the year. the following screenshot shows how it should look:

The screenshot shows the drop down function menu showing the hierarchical table of contents for the NAG Toolbox functions, and a popup help box showing a summary for c05aj that appears if the mouse hovers over that entry. These tooltips may be “torn off” the function browser in which case they stay in the GUI and can be positioned as desired, this is shown by the summary for c05ad that is still visible.

In order to make this work the function browser (apparently) looks for certain named sections on the function's reference page, and so I have had to change the section titles from being numbered to unnumbered, Purpose, Syntax, Description rather than 1 Purpose, 2 Syntax, 3 Description as in the current documentation. The tooltip essentially shows the Description section. I should also note that while versions of the function browser appeared in MATLAB 2008b, the internal details were apparently different and at the time of writing, with the development build, the NAG functions appear in Function Browsers in MATLAB versions starting from 2009b. (The screenshot is of MATLAB 2010a).

Improving the doc function

One way of accessing the function reference page (html documentation) is to open the help browser and then navigate the table of contents or search to obtain the documentation for the function required. This of course works for the NAG Toolbox as well as for built in MATLAB functions. However when working on a function in the command window, an alternative method of accessing documentation is to use the doc function. Unfortunately (for several MATLAB releases) doc is now documented as only working for MATLAB functions built into MATLAB, or supplied in MathWorks Toolboxes. A command such as doc s18ad does not bring up the full html reference page, but opens the help browser with a summary information derived from the ASCII comments in the s18ad.m file, which is the information produced by the command

help
s18ad

. Fortunately the docsearch function does work, and one useful by-product of NAG's function naming scheme is that searching for the function name reliably returns the reference page as the most relevant result, and the reference page is displayed in the help browser. To make the doc and help information more useful, an HTML link invoking MATLAB's docsearch function has been placed into the comments and so, as shown in the following screenshot, in future, while doc s18ad does not bring up the full reference page, there is a hypertext link on the summary page that is displayed that will invoke docsearch s18ad and so the full reference page will be accessible with a single click.

Note that the phrase “NAG Toolbox Help Files” is now a hyperlink, with the status bar showing that clicking on it will activate the MATLAB command docsearch.

Future Improvements

If you have any further ideas for improving the documentation or integrating with currently unused features in the MATLAB GUI, please do let us know!

Reflecting on 40 years of software development

2010-05-18T21:14:00.006+01:00

As regular readers of this blog will know, NAG celebrated its fortieth anniversary on May 13th with a series of informal talks during the day and a dinner party in the evening. Over 200 people - members, collaborators, current and ex-staff - took part in what was a thoroughly enjoyable day. Katie, our Marketing Communications Manager and Queen of the Blogosphere, put together an amazing scrapbook with photographs going back to the early days. Somebody commented to me that a healthy proportion of the pictures involve food and drink - I guess thats just the sort of Company that we are ☺.

Talking to some of the older participants it was clear that - while computing has changed completely in the last four decades - many of the issues which faced the early NAG pioneers are still with us today. Portability remains a key concern, particularly since after a long period of relative stability we are seeing a great divergence in hardware and a smorgasbord of proprietary, non-standard mechanisms for parallel computing. How to satisfy users of different language has always been another issue. Its a little-known fact that the original aim of the NAG project was to produce an Algol 60 library and that the Fortran library that persists to this day was a secondary activity. Over the years we've produced libraries in Ada, Pascal, Fortran 90 and C, with only the latter still under active development. Alongside that work we've produced many bindings for other languages - the first work that I did for NAG was to produce C header files for the Fortran Library while a PhD student at the University of Bath.

Another area where NAG has always been at the forefront is computer-generated documentation. When I first visited NAG more than twenty years ago most of the routine documents were produced on a golfball typewriter, although there were computer-readable versions maintained by Mike O'Donohoe at Cambridge University. His ingenuity at laying out complex mathematical formulae using only ASCII characters was quite astonishing! From there we adopted TSSD (a UK-developed contemporary of TeX), LaTeX, SGML and now XML from which we produce both pdf and HTML variants.

I have a sneaking suspicion that when NAG hosts its fiftieth birthday in ten years time, some of these issues will still be causing us headaches. For every problem that is solved, another will arise. Last Thursday we reminisced about the days when RAM was measured in kilobytes and IEEE arithmetic had yet to be invented. In ten years we might find the idea of having only a single CPU on a computer quite quaint, and a new generation of functional programmers will struggle to understand what object-oriented languages are all about. Its hard to predict the future but, whatever the issues, I'm sure that NAG will be working hard to address them and to continue to deliver high-quality software to our customers.

NAG's 40th Anniversary

2010-05-13T10:51:00.005+01:00

NAG might be 40 years old today but the office feels like we are having a 21st birthday party with a real buzz about the place. Absolutely amazed that Jack Dongarra is in the building and about to give a seminar on the Impact of Architecture and Technology for Extreme Scale on Software and Algorithm Design having travelled from the States just to attend our Anniversary celebrations. Hopefully the weather will hold till later when we have afternoon tea and Pimms after hearing from our founding Directors Brian Ford and Steve Hague on how they created NAG and made it the success it is today.

Jack Dongarra Presenting at NAG's 40th Anniversary Open Day

Congratulations NAG Prize 2010 Winners!

2010-05-03T10:05:00.035+01:00

Last week some of my colleagues attended the first Manchester SIAM Student Chapter conference at the University of Manchester. Going by the feedback it was a really good day with some great presentations by Professor Sven Hammarling (NAG), Professor Bill Lionheart and Dr Oliver Dorn (both of the School of Mathematics, University of Manchester) as well as many student presentations - so well done to those who organised the event.

John Holden (centre) and Craig Lucas (right) chat to an attendee of at the Manchester SIAM Student Conference

The NAG Student Prize was established in 2006 to recognise the long standing relationship NAG has with Manchester and reward the highest achieving students in maths and computation related courses. Mike Dewar, NAG Chief Technology Officer awarded this years prizes at a ceremony later in the day.We’re delighted that this years prizes go to Xiaolu Han and Michael Crabb for their outstanding work in the Mathematical Finance and Mathematics and Computational Science MSc courses respectively.

NAG Prize 2010 Winners Xiaolu Han and Michael Crabb with Mike Dewar, NAG Cheif Technical Officer (centre)

Congratulations and well done to both Xiaolu and Michael.

Waiting for software and school children

2010-04-29T13:32:00.006+01:00

I spent a few hours with colleagues at a major national laboratory earlier this week talking about national security, the role of software and school children. On the surface, they would seem to have little connection but as I drove back to the office the connections started to form. I invite you to follow along. In North America, Europe and other major economies throughout the world we spend huge sums of money in an attempt to detect and prevent national border intrusions and terrorist activities. In the process we've created a massive industry to x-ray, scan and monitor everything that moves in airports and elsewhere. In the process, we've put hundreds of thousand of people into jobs that apparently weren't needed 40 years ago, spent billions in taxes and worst of all, thrown sand into the gears of both business productivity and leisure time by a process ostensibly designed to keep us safe. We've created a huge "tax" in time and money when the real problem is to predict, detect and prevent a few "outlier" events from occurring. My question: where's the software? Isn't the promise of software to be able to predict and detect things on which we have mind-numbing amounts of data? In some respects, the situation in security is much like the situation in HPC: the hardware is constantly evolving and relatively easy to buy and software is difficult and time-consuming to write, validate and document. You were wondering where the children came in? Ah yes. While I had lunch in the lab's common dining area, two bus loads of primary school children on a field trip to the lab arrived to eat their lunch. No doubt they had spent the morning oohing and ahhing over particle accelerators and superconducting magnets. I couldn't help but think how we could get them interested in working on the software that is part of the fabric of modern society. Every time I go through the "remove coat, shoes, laptop, liquids, mobile phone and walk slowly through the scanner dance" I wonder if we in software will find a solution or if we'll have to wait for them to do it for us.

Plotting Lambert's W Function

2010-04-22T12:21:00.007+01:00

For this post I wanted to come up with something a little easier on the eye than my previous offerings (Loading DTDs using DOM in Python, Colorised svn diffs and How do you like your Fortran 95 interfaces?).

Here you go:

I've been developing code for computing complex values of Lambert's W function, and I thought I'd share the Riemann-surface plots above, which are, from top to bottom, Re(W), Im(W) and abs(W) for the branches -1 (Blues), 0 (Greens) and 1 (Oranges). The route I've taken to get the plots has been a good learning experience for me.

First approach

My initial idea was to write a mex wrapper to my algorithmic Fortran code, along the lines of what happens for the NAG Toolbox for MATLAB. I also wanted to see if I could provide a vectorized interface, as native MATLAB functions have.

Whew, it's not much fun debugging mex files, is it?

The main thing I learned from this experience was to use mxCreateNumericMatrix to create my function's return values: in this case the W value, the residual error and the fail flag.

/* prhs[2] is the array of input zs */
dims = mxGetDimensions(prhs[2]);
plhs[0] = mxCreateNumericMatrix(dims[0], dims[1], 
                                mxDOUBLE_CLASS,
                                mxCOMPLEX);

if (nlhs >= 2)
{
  plhs[1] = mxCreateNumericMatrix(dims[0], dims[1],
                                  mxDOUBLE_CLASS,
                                  mxREAL);
  
  if (nlhs >= 3)
    {
      plhs[2] = mxCreateNumericMatrix(dims[0], dims[1],
                                      mxINT_CLASS,
                                      mxREAL);
    }

}

Other approaches I had tried, such as populating local mxArray *s which were only assigned back to the plhs arrays at the end (and if required) caused all sorts of upset to MATLAB when it exited my mex routine.

With the working code, I managed to get some good plots in a MATLABy way. I felt a little discontent though, and then I remembered the advice of my old (and sadly departed) A-level statistics teacher, about how it's better to be able to solve a problem two ways than one.

Second approach

I thought back to Mike Croucher's PyNAG tutorials on calling the NAG C Library from Python, and decided it might be fun to call my algorithmic Fortran code from Python and use matplotlib to do the plotting. While I was at it I decided I'd see if I could get Fortran 2003's C Interoperability to help me make the call to Fortran from Python more straightforward; in the 'unassisted' case, calling Fortran from Python is very similar to calling Fortran from C, with the usual headaches that entails.

The primary hoop I had to jump through to allow C interoperability was to map Fortran character strings to arrays:

...
   USE, INTRINSIC :: iso_c_binding, ONLY : c_char
   ...
   CHARACTER (kind=c_char), DIMENSION (6), &
      INTENT (IN) :: in_str
   CHARACTER (kind=c_char), DIMENSION (10), &
      INTENT (OUT) :: out_str
   CHARACTER (len=SIZE(in_str)) :: local_in_str
   CHARACTER (len=SIZE(out_str)) :: local_out_str
   INTEGER :: i
   ...
   DO i = 1, SIZE(in_str)
      local_in_str(i:i) = in_str(i)
   END DO
   ...
!  Here's where all the real stuff happens. Operate on
!  local_in_str and local_out_str.
   ...
   DO i = 1, SIZE(out_str)
      out_str(i) = local_out_str(i:i)
   END DO

END

Then with a similarly-vectorized and Pythonic wrapper around my Fortran call I generated the plots above using

from mpl_toolkits.mplot3d import \
     Axes3D
from matplotlib import \
     cm
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = Axes3D(fig)
n = 100
lower = -1
upper = -lower
threshold = 0.001
X1, Y1 = np.meshgrid(np.linspace(lower, upper, n),
                     np.linspace(lower, -threshold, n))
X2, Y2 = np.meshgrid(np.linspace(lower, upper, n),
                     np.linspace(threshold, upper, n))
cmaps = {-1: cm.Blues,
         0: cm.Greens,
         1: cm.Oranges}

# Loop over some branches
for b in [-1, 0, 1]:
    w = nag_LambertW(X1 + 1j*Y1, b=b, fail=1)[0]
    ax.plot_surface(X1, Y1, w.imag,
                    rstride=1, cstride=1, cmap=cmaps[b])
    w = nag_LambertW(X2 + 1j*Y2, b=b, fail=1)[0]
    ax.plot_surface(X2, Y2, w.imag,
                    rstride=1, cstride=1, cmap=cmaps[b])

plt.show()

I had to do two separate subplots for each branch to avoid matplotlib interpolating across the discontinuity at the branch cuts. Being a bit of a matplotlib newb, I'm now wondering

is there a better way of telling matplotlib to avoid discontinuities?
is it possible to apply a colour map continuously across two subplots?
doesn't plt.title('Title') work with surface plots?

New product launch = busy, but happy marketing department!

2010-04-15T10:44:00.000+01:00

Today is an exciting day for NAG. The launch of the latest release of the NAG Library for SMP & multicore, Mark 22.

A new product launch (or release) is always a busy time for any marketing department. Months of planning and preparing new literature, presentations, collating benchmark data, new advertising campaigns, exhibition stands etc, and then comes the fun part - announcing to the outside world.

As with all marketing activities, there's a plethora of best practice guides in achieving a successful product launch. Our goal in the marketing department at NAG is to ensure that our news reaches every single user, potential user and specific media.

How do we achieve this goal?

The relatively recent emergence of new communication methods, blogging, microblogging and social media have transformed how we communicate with our audience and has made this launch different from those previous to it. What's so great about these new channels to market is that we can communicate and engage directly and instantaneously with those most important to us as an organisation. This is priceless.

Using these new methods in conjunction with more traditional marketing such as press releases and email blasts will serve to create some of the "buzz" that's needed to achieve early mindshare. However, even with the slickest of launch operations, if the product itself falls short of expectation and fails to impress then all the marketing effort is wasted.

Thankfully at NAG, our products reflect our reputation for quality and correctness. All new numerical routines have to undergo stringent quality checks. This vital testing process results in a portfolio of world-class products that simply wouldn't reach launch if they weren't something we could be proud of.

The NAG Library for SMP & multicore launched on 15th April 2010.

Scratch

2010-04-09T17:20:00.004+01:00

How early can you be taught to program? A thought that has struck me since my wife gave birth to our first born son a month ago. There is no denying that the thought processes that have to be learnt in order to program are invaluable and transferrable abilities. Skills such as logical thinking, reasoning, planning, modularity, that are applicable in all forms of jobs not just ones involving programming.

Having only learnt to program three years ago, I wish I had had more of a chance to program in my early years. So if ever my son is interested in programming I might start him off with a kids programming language called scratch that I heard about in a magazine.

Scratch is a pictorial drag and drop language that feels a bit like lego in style. Different types of variables are different shapes and only fit in the right shaped holes. It's not a very powerful language, but a good educational tool for the youngsters as it is easy to use, understand, and they can create a lot of fun interactive programs very easily.

I’m all for encouraging young ones to learn useful skills, and who knows, these new prodigy scratchers may well be updating and maintaining the NAG library in years to come.

Easter egg or not?

2010-04-01T17:47:00.014+01:00

Hi, I'm Marcin, technical support engineer here at NAG and I've recently encountered an interesting problem. I coauthored a recent paper on option pricing with the NAG Library (Option pricing in Excel). My coauthor Jeremy subsequently blogged about it (Extending, imitating and collaborating). Now it’s time for the next step.

While working on a demonstration of NAG random number generators in Excel and their application in option pricing using the Monte Carlo method I noticed a considerable discrepancy in results while using unscrambled Sobol quasirandom numbers. Well, to start with, how do you price a European option when the underlying asset follows geometric Brownian motion? Basically for each random number you simulate what the price of the asset will be in, say T=1 year, moving in one step from the starting time point t=0 to the end point t=T, calculate the payoff, and the average value of all payoffs will be our desired option price (Monte Carlo methods for option pricing). I’m not going to get into the details of Monte Carlo method, that will be covered in a paper I'm currently working on. The idea is that this simulation works very well and gives good convergence to the result obtained using the close-form formula, no surprises here. The stock prices have log-normal distribution, as expected. This is illustrated in the figure below. This histogram groups the stock prices from 20000 simulations.

One could however want to price Asian options, which are path-dependent. The payoff of such options depends not on the final price of the asset (as in European options), but on the average price of the asset throughout the life of the option. We therefore need to generate a path of stock prices, so we’ll move from t=0 to t=T not in one step as previously, but in many steps, say 252 (this is the average number of trading days per year). In practice we generate 252 random numbers, generate a path of the asset price and for the final step we calculate the payoff. This is done for every Monte Carlo simulation, so as you’ve probably noticed we need 252 times more random numbers than in case of the European option. The quantity is not an issue, the Sobol sequence is long enough. What happens is that in case of such pathwise Monte Carlo simulation the distribution is not log-normal, it doesn’t actually resemble any distribution, as can be seen below.

Why is that? Why do the quasirandom numbers fail? The simulation using the same data, but pseudorandom numbers instead of quasirandom, does return stock prices that are log-normally distributed (below).

What is going on here? Well, I know the answer, it took me some time to realise, and it would have taken more time if not the help of my colleagues Jacques and Robert. Do you, dear reader of this blog, know why it is like this? Is this just an easter egg, or an interesting property of quasirandom numbers? Happy Easter!

My early career at NAG…

2010-03-24T10:41:00.003Z

“Don’t worry you’ve done a mathematics and finance degree and have an engineering masters and you’ve programmed in Pascal and can write VB macros you’ll be just fine.”

“You’ve far more technical experience than any previous sales persons (the last manager used to sell used cars!). If you’re not careful you’ll end up in the technical division.”

Well with such resounding support and confidence from my manager how could I fail? Buoyed by a good track record at previous companies I now felt less daunted by this technology shift… industrial and electronic weighing scales to telecoms to numerical software… what could be more natural?

Thinking about my days selling weighing scales, I can recall a previous manager, John Ruskin Utley’s wise words (that isn’t quite his real name, in fact I suspect his middle name didn’t even begin with R).

“Listen to your customers, listen to your distributors/resellers and ensure those in technical development get the market feedback. Rest assured laddy if you don’t hit your sales target and you blame the product those in development will be asking why you didn’t communicate what you wanted when you had the chance.”

What’s so complicated about weighing scales I hear you ask? Let me give you a simple example… in 1998 £1 = 380,000 Turkish lira… clearly an electronic supermarket scale in Turkey is going to have the capability to display a few more digits than an electronic weighing scale in England!

How does this all relate to NAG? Well it’s all about listening to customers and distributors and feeding back market needs or not as you will see!

One of my first overseas trips for NAG was to one of our European distributors (who have now changed ownership and sales staff). Giving the distributor’s sales manager a hard time about lack of sales we chatted through the NAG product portfolio. “Give me your SWOT analysis on our products” I prompted. Distributor said “Well I can’t sell the compiler as the users complain that it throws up too many error messages and won’t allow them to compile their code.”

“I’ll get to the bottom of this”, I said firmly. “Let me feed this back to technical management. I’m sure they’ll want to fix this especially if it is impacting the reputation of our product and limiting sales.”

Tee hee, ha ha. Well you can imagine the reaction of my technical colleagues. Oh how they must have laughed.

The NAG Fortran Compiler is the world’s best checking compiler. Clearly my distributor’s client was an inexperienced programmer probably used to programming with a compiler that had limited error checking capabilities. Of course the distributor (and enthusiastic salesmen = me) were subsequently given training about compilers generally and in particular the NAG Fortran Compiler.

“How fast do you want the wrong answer?” Well I get it now… fast code is only any good if the answer is correct. NAG's Numerical Libraries are very much based on the philosophy… Accuracy first Performance second Of course many of you will know and be grateful that NAG generally achieve both.

A few more anecdotes about compilers. Ian Chivers, Rhymney Consulting and Fortranplus (known Fortran specialist) once told me “of course any decent programmer has at least two compilers. The NAG Fortran Compiler for checking and perhaps another for comparison or perhaps speed.”

Finally some background behind the NAG Fortran Compiler, NAG developed this in its early days of Numerical Library development not because NAG wanted to become a compiler company, but because the compilers available were inadequate to build accurate, robust numerical libraries and by producing our own we were able to advance the Fortan language. In fact NAG produced the world’s first Fortran 90 compiler and still keeps pace with the latest Fortran Standards with lead developer Malcolm Cohen playing an active role in the Fortran Standards committee.

Well that’s quite enough from me in my first blog. More stories from an enthusiastic salesmen next time ;-)