# Using Mathlink with the cernlib FORTRAN package MINUIT

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);
return fmin;
}

#ifdef __cplusplus
}
#endif

#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!!!

-

## migrated from stackoverflow.comMay 17 '12 at 14:05

This question came from our site for professional and enthusiast programmers.

Could you give a link to the fortran file? I have a hard time figuring out what fortran likes (6 space at beginning of line rule). Also, on your C file you called the function myprog without an argument (printf("%f\n", myprog()); ), does C really work like that? I thought myprog required an integer as an input so that it could return a double. – jmlopez May 17 '12 at 22:31
When I tried to make the library libminuit.a I get an error from the file minuit.F. This file "Can't open module file 'd506cm.mod' for reading". Do you know where I can get this? – jmlopez May 17 '12 at 22:42
@jmlopez thanks for your effort. the source files can be downloaded from the link I provided at the end of the post. here is d506cm.mod – unstable May 17 '12 at 23:52
Ok, I got it to work on my Mac, I cant write an answer right now, but, just check that mathematica has access to the .mincards file. Mine didnt because I forgot to include it and I got a linkd error. Anyway, let us know if you've made progress. – jmlopez May 18 '12 at 5:17
I don't have any clue to the mathlink solution yet. However, I called my fortran code successfully using NETLink. I post my memo here. Any suggestion making the mathlink version work would still be appreciated. – unstable May 18 '12 at 8:29

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

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

void MyProg(void) {
double fmin = 200.0; //myprog_(&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["mlprogPrivate"] :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.

-
We don't often get posts with other languages, but for future reference, see my edit to your answer to for how to add the appropriate syntax highlighting for other languages (not all languages are supported though... see this answer for a list of supported ones and their short forms). When a decent highlighter is not available, I generally use lang-none – R. M. May 18 '12 at 18:10
@jmlopez thanks a lot for your effort! I'll have a second look at my code tonight to see wheter it would work. – unstable May 18 '12 at 18:52
Compared with your implementation, I finally identified the problem, which still leaves me puzzled: If I delete the part from #if WINDOWS_MATHLINK to #else, and of course the last #endif in my C wrapper shown in my question, then it works. However, keeping the same code, I can work out the classic addtwo example perfectly. Does anyone have an idea about the reason? – unstable May 18 '12 at 23:19
@unstable, I have never tried making MathLink programs in Windows, but sometime ago when I was reading about it said you need that code you had in the WINDOWS_MATHLINK`. Weird. Anyway, you should try to make it work with the addtwo example, make a call to the fortran code and see if it still works. Just don't forget to add the libraries when you compile. – jmlopez May 18 '12 at 23:29