Friday, 29 January 2010

Reverse Communication

We all know that computers are becoming increasingly parallel. The latest generation of x86 processors have 12 cores and Intel recently announced an experimental 48 core processor. A lot of attention has focussed on the associated reductions in clock speeds which makes software run slower, but less attention has focussed on other issues such as changes in the amounts of cache available and more contention for memory bandwidth.

The fact is that access to memory is fast becoming the scarce resource, rather than CPU cycles. Most mathematical algorithms in the NAG Library (and other similar software) are designed for ease-of-use with the user passing in a collection of data structures describing the problem to be solved and the routine only exiting when it has calculated the answer or an error is encountered. However we do have a number of routines which take a different approach, called reverse communication. Here, the user provides some initial data and when the routine wants more it stops. State is preserved between calls but the routine never "sees" all the input data at once. An example is f11be which solves a sparse system of equations Ax=b. Instead of supplying A and b explicitly, the routine stops when it needs information like the value of Au where u is an intermediate estimate of the solution.

There are two main advantages to this approach. First, the user can carry out other computations with the data while iterating over it, reducing memory access. Secondly, the user can represent thne data any way that they want, rather than being forced to use the layout mandated by the routine being called. Of course there are downsides too - the user may unwittingly choose a representation for their data which is very inefficient, and of course their program will be more complicated. Nevertheless this may be a much more efficient way of using libraries like ours in the future, and its possible that many of our new routines will follow this approach in years to come.

Monday, 18 January 2010

How do you like your Fortran 95 interfaces?

Fortran 95 has been fully implemented for a couple of years now by all the compilers we support, so we'd like to be able to provide a Fortran 95 way of enjoying our algorithmic content. What kinds of features would you expect or like to see from such an endeavour?

At the very least, by the Mark 23 release of our standard Fortran Library all example programs will have been 'F90-ed' to have module procedures for user-supplied routines, allocation of local arrays, and a range of other cosmetic improvements. Here's one of the simpler instances of this, the example program for the Nelder–Mead routine E04CBF: e04cbfe.f90

So that's just calling the regular Fortran (77) Library interface in a more modern context. If this isn't satisfactory then what you're after are more-modern interfaces too. An uber-example of this would be for the sparse NLP solver E04VHF
 SUBROUTINE nag_e04vhf(start,objadd,objrow,prob,usrfun,a,g,xlow,xupp, &
    xnames,flow,fupp,fnames,x,xstate,xmul,f,fstate,fmul,ns,ninf,sinf, &
!         .. Implicit None Statement ..
!         .. Non-Generic Interface Blocks ..
       SUBROUTINE usrfun(status,x,f,g,cuser,iuser,ruser)
!               .. Use Statements ..
          USE nag_precisions, ONLY : wp
!               .. Scalar Arguments ..
          INTEGER, INTENT (INOUT) :: status
!               .. Array Arguments ..
          REAL (kind=wp), OPTIONAL, INTENT (INOUT) :: f(:), g(:)
          REAL (kind=wp), INTENT (INOUT) :: ruser(:)
          REAL (kind=wp), INTENT (IN) :: x(:)
          INTEGER, INTENT (INOUT) :: iuser(:)
          CHARACTER (8), INTENT (INOUT) :: cuser(:)
       END SUBROUTINE usrfun
!         .. Structure Arguments ..
    TYPE (nag_sparse_matrix), OPTIONAL, INTENT (IN) :: a, g
    TYPE (nag_comm), INTENT (INOUT) :: comm
    TYPE (nag_error), OPTIONAL, INTENT (INOUT) :: error
    TYPE (nag_user), OPTIONAL, INTENT (INOUT) :: user
!         .. Scalar Arguments ..
    REAL (kind=wp), OPTIONAL, INTENT (IN) :: objadd
    REAL (kind=wp), OPTIONAL, INTENT (OUT) :: sinf
    INTEGER, INTENT (INOUT)      :: objrow
    INTEGER, INTENT (IN) :: start
!         .. Array Arguments ..
    REAL (kind=wp), INTENT (INOUT) :: x(:)
    REAL (kind=wp), OPTIONAL, INTENT (INOUT) :: f(:), fmul(:)
    REAL (kind=wp), OPTIONAL, INTENT (IN)  :: flow(:), fupp(:), xlow(:), &
    REAL (kind=wp), OPTIONAL, INTENT (OUT) :: xmul(:)
    INTEGER, OPTIONAL, INTENT (INOUT) :: fstate(:), xstate(:)
    CHARACTER (8), OPTIONAL, INTENT (IN) :: fnames(:), xnames(:)
 END SUBROUTINE nag_e04vhf
How can we get to that point? The desirable approach for us from an engineering perspective is to use our XML technologies to generate Fortran 95 wrappers to our algorithmic base. As an example, an XML-generated interface for G01HBF could be
 FUNCTION nag_multi_normal_f90(tail,a,b,xmu,sig,tol,lwk,ifail)
!         .. Use Statements ..
    USE nag_precisions, ONLY : wp
!         .. Function Return Value ..
    REAL (kind=wp)            :: nag_multi_normal_f90
!         .. Scalar Arguments ..
    REAL (kind=wp), OPTIONAL, INTENT (IN) :: tol
    CHARACTER (1), INTENT (IN) :: tail
!         .. Array Arguments ..
    REAL (kind=wp), INTENT (IN) :: a(:), b(:), sig(:,:), xmu(:)
 END FUNCTION nag_multi_normal_f90
That looks a bit like the interface for g01hb in the NAG Toolbox for MATLAB, because we'd be using similar XML technologies to generate this kind of Fortran interface as is used to generate the MATLAB interfaces.

What level of Fortran 95 interfaces would you like to see?
Are there any particular routines you feel would benefit from detailed human-driven, case-by-case design decisions (e.g., optimization, above)?
Alternatively, what do you hate about our current Fortran 77 interfaces?

Friday, 15 January 2010

Lawrence intro / Painleve equations

Introduction Let me introduce myself first. I'm Lawrence Mulholland and I'm Group Leader for a Group called "Numerical Libraries" which includes a large proportion of our internal developers. I look after project development as we move from one major release of the Library based products to the next; so I am currently managing a FL23 project which aims to produce its first implementations of FL23 at the latter part of this year.

Painleve Equations

I have recently been invited to a workshop hosted by ICMS in Edinburgh on Painleve equations:

This seems like an interesting area of research; and since we are looking to introduce a routine for the Hypergeometric function 2_F_1 in FL23, it seems like a logical extension to be investigating the even more challenging Painleve equations which seem to have a fair range of applications.

What I'm not sure about is: how useful would routines that provide particular solutions to the Painleve equations be to our customer base and to the research community in particular ?

Does anyone out there have any interest in this?