Combining MATLAB with NAG and MAPLE

Marko Laine

University of Helsinki, Computing Centre

`Marko.Laine@Helsinki.FI`

This is mostly the text that that appeared in Proceedings of the Workshop on Symbolic and Numeric Computing, Helsinki 1993 with some minor corrections.

GENMEX is a set of tools to automatically generate MEX-files for use with MATLAB. MATLAB is an interactive program for matrix manipulation, numerical computation and scientific visualization. MEX-files in MATLAB provide a way of dynamically loading an executable to MATLAB which can the be called as it were a MATLAB-function. MEX-files make it possible to efficiently and interactively use numerical Fortran libraries from MATLAB. MEX-files and their use is described in MATLAB External Interface Guide

MATLAB-NAG-MAPLE-link

GENMEX provides an automatic generation of MEX-files for routines in the NAG library with default values assigned for control variables (eg. accuracy, IFAIL, ...) and calculated values for workspace variables and dimensional parameters (leading dimensions in Fortran arrays). GENMEX generates a MEX-file which is linked with corresponding NAG-routine and an M-file in which the default are assigned. An M-file is text file written in MATLAB script language.

Argument subprograms for NAG routines are entered as MATLAB-functions which can therefore be internal MATLAB functions, M-files or maybe MEX-files calling other NAG- or Fortran written routines.

The role of Maple is for automatic code generation for derivatives needed with many NAG routines. MATLAB functions are provided with which the user can write equations in Maple format and then Maple is called to generate code for wanted functionals. The generated Fortran function is compiled into another MEX-function that can directly be called from MATLAB and that can be given as an argument to GENMEX generated NAG-function. See examples below.

GENMEX function

The MATLAB function GENMEX does the actual generation. Because MEX-files, at least on the SUN, are linked statically and are usually quite large, it is economic to generate the MEX-files to some temporary directory which is cleaned up periodically such as /tmp or /usr/tmp.

GENMEX --- Generate MEX-file for NAG-routine

Purpose

Generates and compiles a MEX-file written in C and a cover M-file to specified NAG-Fortran-library routine.

Syntax

genmex('routine') or genmex routine

Example

Calculate (NAG-example)

$\int_0^{2\pi}{x\sin(30x)\over\sqrt{1-({x\over 2\pi})^2}}\,dx$

using routine D01AJF.

>> genmex d01ajf
Writing d01ajf_.c
Writing d01ajf.m
[result,abserr]=d01ajf(f,a,b,optargs)
compiling ...
>> type fun.m      % we have written an m-file for argument subprogram

function y=fun(x)
y=x*sin(30*x)/sqrt(1-x^2/(4*pi^2));

>> [result,abserr]=d01ajf('fun',0,2*pi);
result =
   -2.5433

abserr =
   1.2752e-05

Bugs

Not all NAG routines are available at the present.

Parameters

The generated MATLAB cover function (M-function) for library routine has syntax

[out1, out2, ... ]=routine(in1, in2, ... ,optargs)

The input arguments in1, in2, ... are those which are declared of type IN or INOUT in routine defaults file. Output parameters out1, ... are those declared of type OUT or INOUT. If an argument has a default value it can be altered in optional string argument optargs. This style of giving parameters emulates the so-called named arguments that can be used for example in Mathematica or S-PLUS. optargs is a MATLAB string and it has the form

'var1=value1,var2=value2,...'

where the names of variables correspond to those given in specification file and which also should match the names in NAG-manuals and naghelp-online manuals. Example

[result,abserr]=d01ajf('myfun',0,1,'epsabs=0.01,ifail=1');

A parameter can also be declared as WORK and it is then allocated automatically in M-file but not given as input or output.

Defaults file

GENMEX gets its information about library routines from a defaults file. This is a textual database which contains a record for each routine. Here is an example for routine d01akf, one-dimensional quadrature

          d01akf:S:12
                f:RF(REAL):IN:
                a::IN:
                b::IN:
                epsabs::IN:0.0
                epsrel::IN:0.00001
                result::OUT:
                abserr::OUT:
                w:REAL(lw):WORK:
                lw:INT:IN:2000
                iw:INT(liw):WORK:
                liw:INT:IN:lw/8+2
                ifail:INT:IN:-1

The four fields for each 12 arguments contain

    name:type(dimension):in/out:default value

Here the input arguments with no default values are f, a, and b. Output arguments are result and abserr. The fist line of generated M-file which defines the function arguments for MATLAB function becomes:

  [result,abserr]=d01akf(f,a,b,optargs)

The defaults are assigned in generated M-file, which is ordinary text file and can therefore easily be edited to supply more suitable defaults. Defaults can also be changed when calling the function in MATLAB using the last argument optargs. See the examples below.

An simple editor is also provided for editing or creating default files. At this moment defaults are available for most NAG routines thanks to Dr. Terry Robb who allowed us to use the defaults in his Mathematica-NAG-IMSL-package InterCall (See article of P.C. Abbott in these proceedings).

MAPLE-link

Many NAG-routines need symbolic derivatives. The routines must be provided with a Fortran subroutine which calculates derivatives, gradients, Hessians or Jacobian of the user function we are working with.

MAPLE can be used to automatically generate Fortran code for a function once it is given in MAPLE-form. This link contains MATLAB functions that call MAPLE to generate Fortran-MEX-code for wanted functionals. This work is done by Antti Nevanlinna at the University of Helsinki Computing Centre. Pirkka Peltola describes in these proceedings his work on automatic code generation with Maple, which can also be integrated with GENMEX.

We have following MATLAB/MAPLE functions:

deriv (partial) derivative

hss Hessian matrix

jacob Jacobian matrix

grad Gradient

Example

>> jacob('sin(x),cos(y)','x,y','jac1');
>> help jac1

  Jacobian matrix of [sin(x),cos(y)]
  variables: [x,y]

  Usage: result = jac1(x);  % x is a vector of length 2

>> jac1([(2*pi) (pi/2)])

ans =
     1     0
     0    -1

GENMEX examples

The use of GENMEX is illustrated with two examples. They are both from NAG manual and further information on these routines and their parameters can be acquired there.

E04ABF/BBF Minimum of function of one variable

Find the minimum of the function sin(x)/x

Routine e04abf wants the user function in the form

SUBROUTINE FUNC(X,Y)

so we must write the M-function:

function [x,y]=sinxx(x,y)
y=sin(x)./x;

Plotting the function suggests that the minimum is somewhere between [4,5] so we use

>> [x,y]=sinxx(0:0.1:10); plot(x,y)
>> [x0,f]=e04abf('sinxx',4,5)
x0 =
    4.4934
f =
   -0.2172

For the routine e04bbf we must supply the derivative of the function:

>> deriv('sin(x)/x','x','x','sind')
>> type sinxxd.m   % this the argument subprogram for e04bbf
function [x,y,d]=sinxxd(x,y,d);
y=sin(x)/x;
d=sind(x);

>> [x,f]=e04bbf('sinxxd',4,5)
x =
    4.4934
f =
   -0.2172

Solving initial-value problems

We solve NAG example initial value problem using d02ebf. First generate MEX and M-files for d02ebf.

>> genmex d02ebf
Routine name: d02ebf_ is of type: S and has 13 parameters
     x: type:        R, i/o:  IN
  xend: type:        R, i/o:  IN
     n: type:        I, i/o:  IN, default: rows(y)
     y: type:     R(n), i/o: INOUT
   tol: type:        R, i/o:  IN, default: 10^(-6)
irelab: type:        I, i/o:  IN, default: 0
   fcn: type:        S(R,R(*n),R(*n))
  mped: type:        I, i/o:  IN, default: 1
pederv: type:        S(R,R(*n),R(*n,*n))
output: type:        S(R,R(*n))
     w: type:  R(n,iw), i/o: WORK
    iw: type:        I, i/o:  IN, default: n*(12+n)+50
 ifail: type:        I, i/o:  IN, default: -1

The argument subprograms fcn and fcnder calculate right hand sides of the system and its Jacobian for given x and y.

function [t,y,f]=fcn(t,y,f)
% nag d02ebf example
f(1) = -0.04*y(1) + 1.0e4*y(2)*y(3);
f(2) = 0.04*y(1) - 1.0e4*y(2)*y(3) - 3.0e7*y(2)*y(2);
f(3) = 3.0e7*y(2)*y(2);

function [x,y,pw]=fcnder(x,y,pw)
% nag d02ebf example derivative
pw(1,1) = -0.04d0;
pw(1,2) = 1.0d4*y(3);
pw(1,3) = 1.0d4*y(2);
pw(2,1) = 0.04d0;
pw(2,2) = -1.0d4*y(3) - 6.0d7*y(2);
pw(2,3) = -1.0d4*y(2);
pw(3,1) = 0.0d0;
pw(3,2) = 6.0d7*y(2);
pw(3,3) = 0.0d0;

The argument subprogram output decides what points of the solution are computed (xnodes) and saves them to matrix named solution (This is adapted from InterCall manual example).

function [x,y]=output(x,y)
% output routine for d02ebf
global solution xnodes
[r,c]=size(solution);
solution(r+1,:)=[x,y'];
x=xnodes(1);
xnodes=xnodes(2:length(xnodes));

In MATLAB we write:

>> x0=0; y0=[1 0 0]; xend=10;
>> xnodes=[0.0001 0.001 0.01 .1 .2 1 3 5 7 10 11];
>> solution=[];
>> global solution xnodes
>> d02ebf(x0,xend,y0,'fcn2','fcn2der','output','mped=0,irelab=2,tol=0.0001');

Now we have the solution in the matrix solution and can for example plot it

>> for i=1:3;subplot(2,2,i);plot(solution(:,1),solution(:,i+1));end

Firure of the solution of a stiff system

Solutions of a stiff system

To do

One thing I would like to add to GENMEX is automatic generation of loops on specified routine arguments. While MATLAB is essentially vector oriented most NAG routines, eg. those for special functions operate on single scalar. It is easy to add loop on argument to cover M-file by hand afterwards if one needs it but it would be much more efficient to do the looping in the generated C-code.

Generating MEX-file for users own Fortran code could also be made easier. Users with information of Fortran and MEX programming can easily add their own modules to the system but this could be automated. Portability to other environments where these products are available should be checked.

We are currently integrating this system to the ESC-environment (See Pirkka Peltola's article). The idea is to replace the use of APL with perhaps more user friendly MATLAB.

Availability

The GENMEX package is available from the author. It seems to be in an quite usable stage at this point of development though it has not been used very extensively. We have SUN SPARCStation with SUNOS 4.4.1, MATLAB 4.1, NAG-Fortran-library MARK 15 and Maple V release 1. GENMEX is written in ANSI-C and compiles at least under GNU-C compiler.

References

Heikki Apiola, Marko Laine, and Pirkka Peltola: Unix ESC introductory guide, environment for scientific computing. Guides of computing centre, Computing Centre, University of Helsinki, 1993. No 19.

The MathWorks, Inc: MATLAB Reference Guide, User's Guide, External Interface Guide, New Features Guide, 1992. MATLAB version 4.1.

The Numerical Algorithms Group, Ltd: NAG Fortran Library Manual, Mark 15, 1991.

Terry D. Robb: InterCall manual. Analytica, Nedlands, Australia, 1993.

Marko Laine <Marko.Laine@Helsinki.FI>

`deriv`	(partial) derivative
`hss`	Hessian matrix
`jacob`	Jacobian matrix
`grad`	Gradient