Using Mathlink with the cernlib FORTRAN package MINUIT

Question

I could successfully use Mathlink to call a simple fortran subroutine. However, when I went to a realistic problem, I failed. To be specific, I want to fit a set of data with some function using the cernlib package MINUIT (for function minimization and error analysis). Then within FORTRAN, I can get the best parameter sets and their errors, and many more.

Now, I want to plot some figures with the resulting best parameters in Mathematica. Of course, I can read the parameters by eyes... But it would be nice to do this with Mathlink. However, I failed. I created an executable, but it cannot be recoganized by Mathlink. Install simply hang until I aborted it. This is a bit strange, since the fortran code was successfully called in C when Mathlink was not using. Here is the shortened FORTRAN code, to be called myprog.f, (still somewhat lengthy) for a MINUIT example here:

  ! n is not used, just to give the Mathlink function an argument
  subroutine myprog(n, fmin) 
  !program myprog  
  implicit none
  integer SYSIN/1/, SYSRD/5/, SYSWR/6/, SYSPU/7/
  integer n, ndat, npari, nparx, istat
  real(8) fmin, fedm, errdef
  external myfcn

  !     initialize unit numbers
  call mintio(SYSRD, SYSWR, SYSPU)
  !driven card running the fit
  open(unit=SYSRD, file='myprog.mincards', status='OLD') 
  !output, if commented out, the output will be written to terminal
  open(unit=SYSWR, file='myprog.out', status='UNKNOWN')      

  !     invoke minuit
  call minuit(myfcn, 0)      
  ! current status of minimization
  CALL mnstat(fmin, fedm, errdef, npari, nparx, istat)
  end


  subroutine myfcn(npar, g, f, x, iflag, futil)
  implicit none
  integer npar, iflag, i, ndat, n
  integer SYSIN/1/, SYSRD/5/, SYSWR/6/, SYSPU/7/
  parameter (ndat=5)
  real*8 f, g(*), x(*), xdat(ndat), ydat(ndat), edat(ndat), func,        
 &       answer/15.7/, rms, avg, diff, futil
  !     Save all data that must be preserved between calls!
  save xdat, ydat, edat         ! MUST

  if (iflag .eq. 1) then
  !     Initialization mode
     do i=1,ndat
        xdat(i) = i
        ydat(i) = answer + i - float(ndat + 1) / 2.
        edat(i) = 1.
     end do
  end if
  !     compute LSQ
  f = 0
  do i = 1, ndat
     f = f + (ydat(i) - func(xdat(i), x(1)))**2 / edat(i)
  end do
  end

  function func(x, par)
  implicit none
  real*8 x, par(*), func
  func = par(1)
  end

The driven card myprog.mincards, where the initial values of the parameters are set, contains

set title
Minuit data cards for myprog example program
parameters
1 'average' 0.0 1.0

set print 1
mig
mino
return
stop    

The subroutine myprog(n,fmin) (fmin gives the minimized value of chi-square) can be successfully called in C with the following code

#include "stdio.h"
#ifdef __cplusplus
extern "C"{
#endif

double myprog(int n) 
{
  double fmin;  
  myprog_(&n, &fmin);
  return fmin;
}

#ifdef __cplusplus
}
#endif

int main(void)
{
printf("%f\n", myprog());
return 0;
}

The result is 10. I thought as long as the FORTRAN subroutine can be called in C, one should be able to call it in Mathematica with Mathlink. Here is the C wrapper (I paste the complete code since it might be helpful to others):

#include "mathlink.h"

#ifdef __cplusplus
extern "C"{
#endif

double myprog(int n) 
{
  double fmin;

  myprog_(&n, &fmin);
  MLPutDouble(stdlink, fmin);
  return fmin;
}

#ifdef __cplusplus
}
#endif

#if WINDOWS_MATHLINK

#if __BORLANDC__
#pragma argsused
#endif

int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
{
char  buff[512];
char FAR * buff_start = buff;
char FAR * argv[32];
char FAR * FAR * argv_end = argv + 32;

hinstPrevious = hinstPrevious; /* suppress warning */

if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;
MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
return MLMain( (int)(argv_end - argv), argv);
}

#else

int main(int argc, char* argv[])
{
return MLMain(argc, argv);
}

#endif

The Mathlink template file myprog.tm is

:Begin:
:Function:       myprog
:Pattern:        Minuit[n_Integer]
:Arguments:      { n }
:ArgumentTypes:  { Integer }
:ReturnType:     Manual
:End:

:Evaluate: Minuit::usage = "Minuit[n] does fit using MINUIT. n can be an arbitary integer. This function returns the value of the least chi-square."

I use Windows 7 with MinGW. I wrote math.bat to create an executable for Mathlink

gfortran -c myprog.f
gcc -c callfmath.c -o myprogc.o
mprep myprog.tm -o myprogtm.c
gcc -c myprogtm.c -o myprogtm.o
gcc *.o libminuit.a libML32i3.a -lm -lpthread -mwindows -lstdc++ -lgfortran -o myprog_math
rm *.o

where libminuit.a is the library for MINUIT.

But the generated myprog_math.exe does not do what I want. It does not create a link to Mathematica, and is not Mathlink-aware...

If anybody can tell me what was wrong with the procedure, I would be appreciate it very much!!!

EDIT: The MINUIT source files can be downloaded here, and libminuit.a made by gfortran.

Answer 1

Here I describe the steps I took to make it work on my Macbook Pro using Mac OS X 10.7.4. It shouldn't be very different on Windows but the idea is the same. I'll try to make this post as simple as possible explaining all the details that I can.

Fortran Libraries

Lets begin by compiling the fortran libraries. To do this you will need to download the MINUIT routines. I will be working on the directory /Users/jmlopez/Desktop/MathLink. In this directory I made a directory called fortran. The MINUIT download will be placed in this fortran folder.

Once you have extracted the MINUIT folder from the zip file you will see a lot of fortran files. One in particular is called NEWd506cm.FOR. You will need to change the name of this file to d506cm.F

Now, with your terminal (Mac and Linux) navigate to this folder and enter the command:

gfortran -c *.F

This will now compile all the fortran files. The next step is to put the .o files into a library. We can do this by entering

ar rvs ../libminuit.a *.o

This command puts all the .o files together in the file libminuit.a which is located in one directory above the current one. Now we need your fortran file that makes use of this library. I will post it again here hoping that if someone copies and pastes the code into a file it will work with no problems, I also got rid of the argument n since you clearly do not need it there. Here is the file myprog.f

      ! FILE: myprog.f
      subroutine myprog(fmin) 
      !program myprog  
      implicit none
      integer SYSIN/1/, SYSRD/5/, SYSWR/6/, SYSPU/7/
      integer ndat, npari, nparx, istat
      real(8) fmin, fedm, errdef
      external myfcn

      !     initialize unit numbers
      call mintio(SYSRD, SYSWR, SYSPU)
      !driven card running the fit
      open(unit=SYSRD, file='myprog.mincards', status='OLD')
      !output, if commented out, the output will be written to terminal
      open(unit=SYSWR, file='myprog.out', status='UNKNOWN')

      !     invoke minuit
      call minuit(myfcn, 0)
      ! current status of minimization
      CALL mnstat(fmin, fedm, errdef, npari, nparx, istat)
      end

      subroutine myfcn(npar, g, f, x, iflag, futil)
      implicit none
      integer npar, iflag, i, ndat, n
      integer SYSIN/1/, SYSRD/5/, SYSWR/6/, SYSPU/7/
      parameter (ndat=5)
      real*8 f, g(*), x(*), xdat(ndat), ydat(ndat), edat(ndat), func
      real*8 answer/15.7/, rms, avg, diff, futil
      !     Save all data that must be preserved between calls!
      save xdat, ydat, edat         ! MUST

      if (iflag .eq. 1) then
      !     Initialization mode
      do i=1,ndat
            xdat(i) = i
            ydat(i) = answer + i - float(ndat + 1) / 2.
            edat(i) = 1.
      end do
      end if
      !     compute LSQ
      f = 0
      do i = 1, ndat
            f = f + (ydat(i) - func(xdat(i), x(1)))**2 / edat(i)
      end do
      end

      function func(x, par)
      implicit none
      real*8 x, par(*), func
      func = par(1)
      end

Now lets compile it this file to obtain the object myprog.o

gfortran -c myprog.f

We should be set with our fortran libraries. At this point I have the files myprog.o and libminuit.a in the directory /Users/jmlopez/Desktop/MathLink/fortran.

C program

Just to make sure that we can work with these libraries let us make a simple C program that does this. Now we will work in the directory /Users/jmlopez/Desktop/MathLink/cxx. Under this directory I will now make the file fortranCall.c:

// FILE: fortranCall.c
#include <stdio.h>

double myprog_(double*);

double myprog(void) {
  double fmin;  
  myprog_(&fmin);
  return fmin;
}

int main(void) {
    printf("%f\n", myprog());
    return 0;
}

Notice that I have declared a prototype for the fortran function. In this way the c compiler will not complain that it doesn't know what myprog_ is. Now we can compile this file by doing:

 gcc -lgfortran ../fortran/libminuit.a ../fortran/myprog.o fortranCall.c -o fortranCall

Now we have our file fortranCall. If we run it we will obtain an error.

Manuel-Lopezs-MacBook-Pro:cxx jmlopez$ ./fortranCall
At line 13 of file myprog.f (unit = 5, file = 'stdin')
Fortran runtime error: File 'myprog.mincards' does not exist

I should note here that if for some reason your program crashes the MathLink connection with Mathematica will die. When I first did it the connection was made but as soon I as I ran it it the connection broke and I had no idea what was happening until I remembered that I did not have the myprog.mincards file. Once we have this file in the same directory we should obtain this:

Manuel-Lopezs-MacBook-Pro:cxx jmlopez$ ./fortranCall
10.000000

One thing to notice is that the program generates the file myprog.out. At the end of the file it says:

CALL TO USER FUNCTION WITH IFLAG = 3

..........MINUIT TERMINATED AND RETURNS TO USER PROGRAM. 

I'm not sure if this is good or bad but it does bring problems to the MathLink program.

Mathematica

I will now be working in /Users/jmlopez/Desktop/MathLink/mma. In this folder I will have 5 files: mlprog.m, MLFile.nb, mlprog.tm, myprog.mincards and mlprog.cpp.

mlprog.cpp

// FILE: mlprog.cpp
#include "mathlink.h"

extern "C" {
    void myprog_(double*);
}

void MyProg(void) {
    double fmin = 200.0; //myprog_(&fmin);
    MLPutReal64(stdlink, fmin);
}

int main(int argc, char* argv[]) {
    return MLMain(argc, argv);
}

Notice how I commented the part where we call the fortran routine. We are doing this to verify that we can get a working MathLink program. This file declares the function MyProg, this function takes no arguments and returns only one value: 200.

mlprog.tm

:Evaluate: BeginPackage["mlprog`"]

:Evaluate: MyProg::usage = "MyProg[], returns a real number. "
:Evaluate: MyProg::argerr = "This function takes no arguments.";
:Evaluate: MyProg[args___] := (Message[MyProg::argerr]; $Failed)

:Evaluate: EndPackage[]

:Evaluate: Begin["mlprog`Private`"]

:Begin:
:Function: MyProg
:Pattern: MyProg[]
:Arguments: {}
:ArgumentTypes: {Manual}
:ReturnType: Manual
:End:

:Evaluate: End[]

In this template file I declared some mma commands so that we can take care of some input errors. Since you don't need any input, don't write anything there.

Finally, the file mlprog.m only needs one line:

Install["mlprog"];

Now we move on to the terminal. Assuming that we are in the directory containing all these files I now do the following:

export MMA=/Applications/Mathematica.app/SystemFiles/Links/MathLink/DeveloperKit/CompilerAdditions

I'm doing this since that particular directory is not in my path and I'm currently too lazy to include it there. This directory contains the mprep binary and other files necessary by mathematica in MAC OS X.

${MMA}/mprep mlprog.tm -o mlprogtm.cpp
g++ -c -I${MMA} mlprogtm.cpp
g++ -c -I${MMA} mlprog.cpp
g++ -I${MMA} mlprogtm.o mlprog.o ../fortran/myprog.o ../fortran/libminuit.a -L${MMA} -lMLi3 -lstdc++ -lgfortran -framework Foundation -o mlprog

Now we have the executable mlprog. Here is a screenshot of my mathematica notebook using that executable:

So far so good. So now let us uncomment that section in the mlprog.cpp file to make use of the fortran routine. We repeat all the steps above to compile the program again and this is how it goes down for me:

What just happened here? The function worked the first time but after that, the connection died. We saw that before using the fortran routine everything was fine. And you have mentioned that you managed to make it work with another simple fortran routine. I guess all this means is that somewhere in the fortran routine mainly from the library, there's an instruction to exit the program which terminates the link automatically. This is what the end of the file it generated says:

 CALL TO USER FUNCTION WITH IFLAG = 3

 ..........MINUIT TERMINATED AND RETURNS TO USER PROGRAM.            
 ***************************************************************************
  MINUIT RELEASE 96.03  INITIALIZED.   DIMENSIONS 100/ 50  EPSMAC=  0.89E-15
 ***************************************************************************
                          MINUIT DATA BLOCK NO.   1
 ***************************************************************************
 stop                                              
 ******************************************************************************
 END OF DATA ON UNIT NO.  5
 TWO CONSECUTIVE EOFs ON PRIMARY INPUT FILE WILL TERMINATE EXECUTION.

 ..........MINUIT TERMINATED BY UNRECOVERABLE READ ERROR ON INPUT. 

It seems that the first time it calls it works fine, but on the second time there is an error. I don't know how to fix this since I'm not familiar with the fortran routines but there you have it. An example of using fortran from Mathematica. I hope that you won't have much trouble adjusting this in windows.

UPDATE

I looked again in the c program. After I ran the program twice I get the same myprog.out file. Now, on the mathematica side, after I ran it once the first time I get the result 10. But myprog.out has nothing in it. After I ran it again, the mathematica tells me that there's been an error and now the file has all the output it should have with the last error thing appended to it. This somehow makes me think that, in the C program, maybe the file gets closed after the program has been finished but it is never closed.

Ha, this is indeed what is happening, in your fortran code there is no closing statement. One thing to keep in mind is that the Mathlink program doesn't exit after a function has been called. It is always there waiting for instructions. So, since you open a file, it is still there in memory and there is no way for the buffer to get released. Ok, I will now attempt to fix the fortran file:

      close(SYSRD)
      close(SYSWR)

Add those two lines at the end of the myprog subroutine. Compile that file and lets try again ...

There you have it. Always remember to put things in their place after you are done using them. If you open a file, remember to close it. If you use some memory, remember to give it back. Now everytime you call the function it will write to file and you will always have the output.