Posts Tagged ‘fortran’

Installing and running fmmlib from MATLAB on Mac OS X

Tuesday, December 26th, 2017

No better way to finish 2017 than trying to compile some old fortran code. This time it’s the fast multipole method library from nyu, fmmlib.

Here’s how I did it on a MacBook Pro running OS X 10.12 with MATLAB r2018a pre-release.

First install gfortran using homebrew via:

brew install gcc

Now we need to convince to use find and use gfortran.

Download gfortran.xml

and

Download gfortran.xml

In MATLAB issue:

copyfile('~/Downloads/xcode7_mexopts/gfortran.xml',fullfile( matlabroot, 'bin', 'maci64', 'mexopts' ))

On the command line, install gfortrani8 in the same directory as gfortran. I did this via:

cp ~/Downloads/gfortrani8 /usr/local/Cellar/gcc/7.2.0/bin/gfortrani8

Now, the following should succeed in MATLAB to setup the gfortran compiler:

mex -setup -v FORTRAN

For posterity: If you get “sdk” related issues, you may need to copy all the 10.13 lines in gfortran.xml and create analogous lines for 10.14 etc.

It seems like the '-compatibleArrayDims' flag should be used with fmmlib, but honestly I’m not sure. Something to revisit if it crashes on giant input I guess.

Finally, I had trouble locating the libgfortran.a library and convincing MATLAB’s Mex to use it. Therefore, I added:

[~,libgfortran] = system('gfortran -print-file-name=libgfortran.dylib');
libgfortran = libgfortran(1:end-1);

to my Mex script so that I can issue lines like this:

mex('fmm3d_r2012a.c',fortran_object_files{:},libgfortran);

This is enough to compile all of the fortran source files one by one and then link them during compilation and linking of the fmm3d_r2012a.c file.

I’ve added a matlab/compile.m script to my GitHub mirror/fork of the fmmlib library. So if you do the steps above you should be to just issue:

cd ~/Downloads
git clone 
cd ~/Downloads/fmmlib3d-1.2/matlab/
compile

Source

ARPACK driver translated from fortran to c++

Wednesday, December 1st, 2010

In an attempt to understand how to make a C++ interface to ARPACK to solve eigenvalue problems. I’ve translated the “real nonsymmetric generalized eigenvalue problem” example (EXAMPLES/NONSYM/dndrv3.f) from fortran into a cpp program.

Here’s my translation, save it in a file called dndrv3.cpp:


/* dndrv3.f translated by Alec Jacobson
*/
#include <assert.h>
#include <cstdio>
#include "f2c.h"

/*     %-----------------------------% */
/*     | Define leading dimensions   | */
/*     | for all arrays.             | */
/*     | MAXN:   Maximum dimension   | */
/*     |         of the A allowed.   | */
/*     | MAXNEV: Maximum NEV allowed | */
/*     | MAXNCV: Maximum NCV allowed | */
/*     %-----------------------------% */
#define maxn 256
#define maxnev 10
#define maxncv  25
#define ldv maxn
// need pointer to ldv constant value
integer LDV = ldv;


/*     %--------------% */
/*     | Local Arrays | */
/*     %--------------% */
integer iparam[11], ipntr[14];
logical select[maxncv];
// Fortran is column-major order!
doublereal ax[maxn],
           mx[maxn],
           d[3*maxncv],
           resid[maxn],
           v[maxncv*ldv],
           workd[3*maxn],
           workev[3*maxncv],
           workl[3*maxncv*maxncv+6*maxncv],
           md[maxn],
           me[maxn-1];

/*     %---------------% */
/*     | Local Scalars | */
/*     %---------------% */
char * bmat;
char * which;
integer ido,
        n,
        nev,
        ncv,
        lworkl,
        info,
        ierr,
        j,
        nconv,
        maxitr,
        ishfts,
        mode;
doublereal tol,
           sigmar,
           sigmai,
           h;
logical first,
        rvec;

/*     %------------% */
/*     | Parameters | */
/*     %------------% */
doublereal zero = 0.0,
           one = 1.0;

/*     %-----------------------------% */
/*     | BLAS & LAPACK routines used | */
/*     %-----------------------------% */
extern "C"{
  extern doublereal dnrm2_(integer *, doublereal *, integer *);
  extern doublereal dlapy2_(doublereal *, doublereal *);
  extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
  integer *);
  extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 
    integer *, doublereal *, integer *);
  extern /* Subroutine */ int dpttrf_(integer *, doublereal *, doublereal *,
    integer *);
  extern /* Subroutine */ int dpttrs_(integer *, integer *, doublereal *, 
    doublereal *, doublereal *, integer *, integer *);

/*     %-----------------------------% */
/*     | other routines used         | */
/*     %-----------------------------% */
/* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
extern /* Subroutine */ int dmout_(integer *, integer *, integer *, 
  doublereal *, integer *, integer *, char *, ftnlen);

extern /* Subroutine */ int av_(integer *, doublereal *, doublereal *);
extern /* Subroutine */ int mv_(integer *, doublereal *, doublereal *);
extern /* Subroutine */ int dnaupd_(integer *, char *, integer *, char *, 
  integer *, doublereal *, doublereal *, integer *, doublereal *, 
  integer *, integer *, integer *, doublereal *, doublereal *, 
  integer *, integer *, ftnlen, ftnlen);
extern /* Subroutine */ int dneupd_(logical *, char *, logical *, 
  doublereal *, doublereal *, doublereal *, integer *, doublereal *,
  doublereal *, doublereal *, char *, integer *, char *, integer *,
  doublereal *, doublereal *, integer *, doublereal *, integer *, 
  integer *, integer *, doublereal *, doublereal *, integer *, 
  integer *, ftnlen, ftnlen, ftnlen);
}

/*     matrix vector multiplication subroutine */
int av_(integer *n, doublereal *v, doublereal *w)
{
  integer j;
  doublereal one, two, dd, dl, du, s, h, rho;
  rho = 10.;
  one = 1.0;
  two = 2.0;

  /*     Compute the matrix vector multiplication y<---A*x */
  /*     where A is stiffness matrix obtained from the finite element */
  /*     discretization of the 1-dimensional convection diffusion operator */
  /*                           d^2u/dx^2 + rho*(du/dx) */
  /*     on the interval [0,1] with zero Dirichlet boundary condition using */
  /*     linear elements. */

  h = one / (doublereal)(*n+1);
  s = rho/two;
  dd = two/h;
  dl = -one/h -s;
  du = -one/h + s;

  w[0] = dd*v[0] + du * v[1];
  for(int j = 1;j<*n-1;j++)
  {
    w[j] = dl*v[j-1] + dd*v[j] + du*v[j+1];
  }
  w[*n-1] = dl*v[*n-2] + dd*v[*n-1];

}

int mv_(integer *n, doublereal *v, doublereal *w)
{
  integer j;
  doublereal one, four, h;
  one = 1.0;
  four = 4.0;

/*     Compute the matrix vector multiplication y<---M*x */
/*     where M is the mass matrix formed by using piecewise linear */
/*     elements on [0,1]. */

  w[0] = four*v[0] + one*v[1];
  for(int j = 1;j<*n-1;j++)
  {
    w[j] = one*v[j-1] + four*v[j] + one*v[j+1];
  }
  w[*n-1] = one*v[*n-2] + four*v[*n-1];

  h = one / (doublereal)(*n + 1);
  // increment of array, need as pointer by dscal_
  integer incx = 1;
  dscal_(n,&h,w,&incx);
}

int main(void){
  /*     %----------------------------------------------------% */
  /*     | The number N is the dimension of the matrix.  A    | */
  /*     | generalized eigenvalue problem is solved (BMAT =   | */
  /*     | 'G').  NEV is the number of eigenvalues to be      | */
  /*     | approximated.  The user can modify NEV, NCV, WHICH | */
  /*     | to solve problems of different sizes, and to get   | */
  /*     | different parts of the spectrum.  However, The     | */
  /*     | following conditions must be satisfied:            | */
  /*     |                    N <= MAXN,                      | */
  /*     |                  NEV <= MAXNEV,                    | */
  /*     |              NEV + 2 <= NCV <= MAXNCV              | */
  /*     %----------------------------------------------------% */
  // Set problem size
  n = 100;
  // Set number of eigen values desired
  nev = 4;
  //  NCV is the largest number of basis vectors that will be used in 
  //  the Implicitly Restarted Arnoldi Process.  Work per major iteration
  //  is proportional to N*NCV*NCV.
  ncv = 20;
  assert(n<=maxn);
  assert(nev<=maxnev);
  assert(ncv<=maxncv);
  bmat = "G";
  which = "LM";

/*     %------------------------------------------------% */
/*     | M is the mass matrix formed by using piecewise | */
/*     | linear elements on [0,1].                      | */
/*     %------------------------------------------------% */

  h = one / (doublereal)(n+1);
  for(int j = 0;j<n-1;j++)
  {
    md[j] = 4.0 * h;
    me[j] = one*h;
  }
  md[n-1] = 4*h;

  // Factor md using LAPACK's symmetric positive definite factor routine
  dpttrf_(&n,md,me,&ierr);
  if( ierr != 0 )
  {
    printf("ERROR with _pttrf\n");
    exit(1);
  }

/*     %-----------------------------------------------------% */
/*     | The work array WORKL is used in DNAUPD as           | */
/*     | workspace.  Its dimension LWORKL is set as          | */
/*     | illustrated below.  The parameter TOL determines    | */
/*     | the stopping criterion. If TOL<=0, machine          | */
/*     | precision is used.  The variable IDO is used for    | */
/*     | reverse communication, and is initially set to 0.   | */
/*     | Setting INFO=0 indicates that a random vector is    | */
/*     | generated in DNAUPD to start the Arnoldi iteration. | */
/*     %-----------------------------------------------------% */
  lworkl = 3*ncv*ncv+6*ncv;
  tol    = 0.0;
  ido    = 0;
  info   = 0;

/*     %---------------------------------------------------% */
/*     | This program uses exact shifts with respect to    | */
/*     | the current Hessenberg matrix (IPARAM[0] = 1).    | */
/*     | IPARAM[2] specifies the maximum number of Arnoldi | */
/*     | iterations allowed.  Mode 2 of DNAUPD is used     | */
/*     | (IPARAM[6] = 2).  All these options can be        | */
/*     | changed by the user. For details, see the         | */
/*     | documentation in DNAUPD.                          | */
/*     %---------------------------------------------------% */

  // use the exact shift strategy 
  // ISHIFT = 1: exact shift with respect to the current
  //   Hessenberg matrix H.  This is equivalent to
  //   restarting the iteration from the beginning
  //   after updating the starting vector with a linear 
  //   combination of Ritz vectors associated with the  
  //   "wanted" eigenvalues. 
  ishfts = 1;
  maxitr = 300;
  mode = 2;

  iparam[0] = ishfts;
  iparam[2] = maxitr;
  iparam[6] = mode;

/*     %-------------------------------------------% */
/*     | M A I N   L O O P (Reverse communication) | */
/*     %-------------------------------------------% */
  while(true)
  {

/*        %---------------------------------------------% */
/*        | Repeatedly call the routine DNAUPD and take | */
/*        | actions indicated by parameter IDO until    | */
/*        | either convergence is indicated or maxitr   | */
/*        | has been exceeded.                          | */
/*        %---------------------------------------------% */
    dnaupd_(
      &ido,
      bmat,
      &n,
      which,
      &nev,
      &tol,
      resid,
      &ncv,
      v,
      &LDV,
      iparam,
      ipntr,
      workd,
      workl,
      &lworkl,
      &info,
      (ftnlen)1,
      (ftnlen)2);
    if (ido == -1 || ido == 1)
    {
/*           %----------------------------------------% */
/*           | Perform  y <--- OP*x = inv[M]*A*x      | */
/*           | The user should supply his/her own     | */
/*           | matrix vector routine and a linear     | */
/*           | system solver.  The matrix-vector      | */
/*           | subroutine should take workd(ipntr(1)) | */
/*           | as input, and the final result should  | */
/*           | be returned to workd(ipntr(2)).        | */
/*           %----------------------------------------% */
      // Perform matrix vector multiply: A*word[ipntr[0]-1] into
      // workd[ipntr[1]-1]
      av_(&n, &workd[ipntr[0] - 1], &workd[ipntr[1] - 1]);
      // dpttrs_ needs pointer to NRHS, number of right hand sides
      integer num_rhs= 1;
      // DPTTRS solves a tridiagonal system of the form
      //    M * Y = A*X
      //    you can read this as computing y:
      //    Y = M^-1 * A * X
      //    workd[ipntr[1]-1] is both the input (A*X) and the output (Y)
      // using the L*D*L' factorization of A computed by DPTTRF.  D is a
      // diagonal matrix specified in the vector D, L is a unit bidiagonal
      // matrix whose subdiagonal is specified in the vector E, and X and B
      // are N by NRHS matrices.
      dpttrs_(&n, &num_rhs, md, me, &workd[ipntr[1] - 1], &n, &ierr);
      if( ierr != 0 )
      {
        printf("ERROR with _pttrs\n");
        exit(1);
      }
    }else if ( ido == 2 ) 
    {

/*           %-------------------------------------% */
/*           |        Perform  y <--- M*x          | */
/*           | The matrix vector multiplication    | */
/*           | routine should take workd(ipntr(1)) | */
/*           | as input and return the result to   | */
/*           | workd(ipntr(2)).                    | */
/*           %-------------------------------------% */
      mv_(&n, &workd[ipntr[0] - 1], &workd[ipntr[1] - 1]);
    }
/*     %-----------------------------------------% */
/*     | Either we have convergence, or there is | */
/*     | an error.                               | */
/*     %-----------------------------------------% */

    else if (info < 0) 
    {

/*        %--------------------------% */
/*        | Error message. Check the | */
/*        | documentation in DNAUPD. | */
/*        %--------------------------% */
      printf("Error with _naupd, info = %d\n"
        "Check the documentation of _naupd.\n",
        info);
      exit(1);
    }else{
/*        %-------------------------------------------% */
/*        | No fatal errors occurred.                 | */
/*        %-------------------------------------------% */

        // Break out of reverse communication loop
      break;
    }
  }

/*        %-------------------------------------------% */
/*        | Post-Process using DNEUPD.                | */
/*        |                                           | */
/*        | Computed eigenvalues may be extracted.    | */
/*        |                                           | */
/*        | Eigenvectors may also be computed now if  | */
/*        | desired.  (indicated by rvec = .true.)    | */
/*        %-------------------------------------------% */
  rvec = true;
  dneupd_(
    &rvec, 
    "A", 
    select, 
    d, 
    &d[0+1*maxncv], 
    v, 
    &LDV, 
    &sigmar, 
    &sigmai,
    workev, 
    bmat, 
    &n,
    which, 
    &nev, 
    &tol, 
    resid, 
    &ncv, 
    v, 
    &LDV,
    iparam, 
    ipntr,
    workd,
    workl,
    &lworkl,
    &ierr,
    (ftnlen) 1,
    (ftnlen)1, 
    (ftnlen)2);

/*        %-----------------------------------------------% */
/*        | The real part of the eigenvalue is returned   | */
/*        | in the first column of the two dimensional    | */
/*        | array D, and the IMAGINARY part is returned   | */
/*        | in the second column of D.  The corresponding | */
/*        | eigenvectors are returned in the first NEV    | */
/*        | columns of the two dimensional array V if     | */
/*        | requested.  Otherwise, an orthogonal basis    | */
/*        | for the invariant subspace corresponding to   | */
/*        | the eigenvalues in D is returned in V.        | */
/*        %-----------------------------------------------% */
  if( ierr != 0 )
  {

/*           %------------------------------------% */
/*           | Error condition:                   | */
/*           | Check the documentation of DNEUPD. | */
/*           %------------------------------------% */
    printf("Error with _neupd, info = %d\n"
      "Check the doumentation of _neupd\n",
      ierr);
  }else
  {
    first = true;
    nconv = iparam[4];
    for(int j = 0; j<iparam[4];j++)
    {
/*              %---------------------------% */
/*              | Compute the residual norm | */
/*              |                           | */
/*              |  ||  A*x - lambda*M*x ||  | */
/*              |                           | */
/*              | for the NCONV accurately  | */
/*              | computed eigenvalues and  | */
/*              | eigenvectors.  (iparam(5) | */
/*              | indicates how many are    | */
/*              | accurate to the requested | */
/*              | tolerance)                | */
/*              %---------------------------% */
      // need point to x,y increment values for daxpy
      integer incx = 1;
      integer incy = 1;
      if( d[1*maxncv+j] == zero)
      {

/*                 %--------------------% */
/*                 | Ritz value is real | */
/*                 %--------------------% */
        // Multiply A * v(:,j) -> ax
        av_(&n, &v[j*ldv + 0], ax);
        // Multiply M * v(:,j) -> mx
        mv_(&n, &v[j*ldv + 0], mx);
        // need pointer to minus d[0,j]
        doublereal minus_d_0_j = -d[0*maxncv+j];
        // ax <-- -d[0,j] * mx + ax
        daxpy_(&n, &minus_d_0_j, mx, &incx, ax, &incy);
        // d[2,j] <-- ||ax||
        d[2*maxncv+j] = dnrm2_(&n,ax,&incx);
        d[2*maxncv+j] /= dabs(d[0*maxncv + j]);
      }else if (first) 
      {

/*                 %------------------------% */
/*                 | Ritz value is complex  | */
/*                 | Residual of one Ritz   | */
/*                 | value of the conjugate | */
/*                 | pair is computed.      | */
/*                 %------------------------% */

        // Multiply A * v(:,j) -> ax
        av_(&n, &v[j*ldv + 0], ax);
        // Multiply M * v(:,j) -> mx
        mv_(&n, &v[j*ldv + 0], mx);
        // need pointer to minus d[0,j]
        doublereal minus_d_0_j = -d[0*maxncv+j];
        // ax <-- -d[0,j] * mx + ax
        daxpy_(&n, &minus_d_0_j, mx, &incx, ax, &incy);
        // Multiply M * v(:,j+1) -> mx
        mv_(&n, &v[(j+1)*ldv + 0], mx);
        // ax <-- d[1,j] * mx + ax
        daxpy_(&n, &d[1*maxncv +j], mx, &incx, ax, &incy);
        // d[2,j] <-- ||ax||
        d[2*maxncv+j] = dnrm2_(&n,ax,&incx);
        // d[2,j] <-- d[2,j]^2
        d[2*maxncv+j] *= d[2*maxncv+j];
        // Multiply A * v(:,j+1) -> ax
        av_(&n, &v[(j+1)*ldv + 0], ax);
        // Multiply M * v(:,j+1) -> mx
        mv_(&n, &v[(j+1)*ldv + 0], mx);
        // need pointer to minus d[0,j]
        minus_d_0_j = -d[0*maxncv+j];
        // ax <-- -d[0,j] * mx + ax
        daxpy_(&n, &minus_d_0_j, mx, &incx, ax, &incy);
        // Multiply M * v(:,j) -> mx
        mv_(&n, &v[(j)*ldv + 0], mx);
        // need pointer to minus d[1,j]
        doublereal minus_d_1_j = -d[1*maxncv+j];
        // ax <-- d[1,j] * mx + ax
        daxpy_(&n, &minus_d_1_j, mx, &incx, ax, &incy);
        // d[2,j] <-- sqrt(d[2,j]^2 + ||ax||^2)
        doublereal ax_sqr = dnrm2_(&n,ax,&incx);
        d[2*maxncv+j] = dlapy2_(&d[2*maxncv+j],&ax_sqr);
        // d[2,j] <-- d[2,j] / sqrt(d[0,j]^2 + d[1,j]^2)
        d[2*maxncv+j] = d[2*maxncv+j] / dlapy2_(&d[0*maxncv+j],&d[1*maxncv+j]);
        d[2*maxncv+j+1] = d[2*maxncv+j];
        first = false;
      }else
      {
        first = true;
      }
    }

/*           %-----------------------------% */
/*           | Display computed residuals. | */
/*           %-----------------------------% */
    // need pointer to maxncv
    integer MAXNCV = maxncv;
    integer m = 6, lda = 3, idigit = -6;
    dmout_(&m, &nconv, &lda, d, &MAXNCV, &idigit, 
        "Ritz values (Real,Imag) and relative residuals", 
        (ftnlen)46);

  }

/*        %------------------------------------------% */
/*        | Print additional convergence information | */
/*        %------------------------------------------% */
  if( info == 1 )
  {
    printf("Maximum number of iterations reached.\n");
  }else if(info == 3)
  {
    printf("No shifts could be applied during implicit"
        " Arnoldi update, try increasing NCV.\n");
  }
  printf("_NRDV3\n"
   " ====== \n"
   " \n"
   " Size of the matrix is %d\n"
   " The number of Ritz values requested is %d\n"
   " The number of Arnoldi vectors generated"
   " (NCV) is %d\n"
   " What portion of the spectrum: %s\n"
   " The number of converged Ritz values is %d\n"
   " The number of Implicit Arnoldi update,"
   " iterations taken is %d\n"
   " The number of OP*x is %d\n"
   " The convergence criterion is %g\n"
   " \n",
   n,
   nev,
   ncv,
   which,
   nconv,
   iparam[2],
   iparam[8],
   tol);

/*     %---------------------------% */
/*     | Done with program dndrv3. | */
/*     %---------------------------% */
  return 0;
}

This assumes that you used f2c to compile ARPACK, so you should be able to compile this on a mac with something like:


g++ dndrv3.cpp -o dndrv3_cpp -framework Accelerate -L[path to libarpack.a] -larpack -lf2c

Then if you issue:


./dndrv3_cpp

You should see:


 Ritz values (Real,Imag) and relative residuals
 ----------------------------------------------
               Col   1       Col   2       Col   3
  Row   1:    4.96055E+01   0.00000E+00   9.79801E-01
  Row   2:    4.10548E+02   0.00000E+00   2.73244E-01
  Row   3:    2.00256E+03   0.00000E+00   1.38874E-01
  Row   4:    3.01419E+03   0.00000E+00   7.42380E-02
  
_NRDV3
 ====== 
 
 Size of the matrix is 100
 The number of Ritz values requested is 4
 The number of Arnoldi vectors generated (NCV) is 20
 What portion of the spectrum: LM
 The number of converged Ritz values is 4
 The number of Implicit Arnoldi update, iterations taken is 14
 The number of OP*x is 218
 The convergence criterion is 1.11022e-16

Which is the same as what you should see if you ran the original fortran version (up to double precision print accuracy).

Compiling and using ARPACK++ on Mac OS X

Saturday, November 27th, 2010

Hilariously arpack++ is supposed to make arpack more easily accessible. In the end, once it’s finally working this is probably true. But getting it working is not at all easy.

Arpack++ is a bunch of header files that provide an interface to the amazingly complicated fortran code of arpack.So the first step is to compile arpack itself.

Next follow these instructions to build a superlu static or dynamic library using xcode. If you want a universal library (32-bit and 64-bit) be sure to enable both i386 and x86_64 when building the final release.

Note: I skipped UMFPACK because I knew that I wouldn’t be getting that far into arpack++’s fancy features.

Now follow instructions 9 through 11 of these instructions in order to patch the outdated arpack++ code. I repeat them here:

9. So now download Arpack++ from (http://www.ime.unicamp.br/~chico/arpack++/), and its patch from (http://reuter.mit.edu/index.php/software/arpackpatch/).
10. extract Arpack++ and copy the patch file to the arpack++ folder, and type
patch -p1 < arpack++1.2.patch.diff

I don’t mess with the make files for the examples.

Instead travel to examples/nonsym and build an example with something like:


g++ simple.cc -I../../../include -I ../../matrices/nonsym/ -I../ -o simple ~/ARPACK/libarpack.a -framework Accelerate ~/Downloads/SuperLU_4.1/lib/superlu/build/Release/libsuperlu.a -lf2c

You’ll have to change the paths in the linked libraries to match yours.

Update: Upon trying to compile x86_64 architecture 64-bit executables I notice that I get linker warnings, like:


ld warning: for symbol _debug_ tentative definition of size 192 from /Users/ajx/ARPACK/libarpack.a(dneigh.o) is is smaller than the real definition of size 96 from /var/folders/XJ/XJuT0FQMG4a+QhMSnFcAwU+++TI/-Tmp-//ccoufDAR.o

Then when I execute the programs I get runtime errors like:


Arpack error in Eupp.
-> Error in ARPACK Eupd fortran code.
Arpack error in FindEigenvectors.
-> Could not find any eigenvector.
Eigenvalues:

I guess for now I am bound to the 32-bit i386 version…

Compiling L-BFGS-B nonlinear optimizer for mac os x

Monday, August 16th, 2010

Just running:


make

on the fortran source code of L-BFGS-B gave me linker errors on my mac OS X 10.6 machine:


invalid parameter driver1.o
Undefined symbols:
  "_s_wsle", referenced from:
      _subsm_ in routines.o
    ....

I fixed this by installing g77, following these instructions.

Then changed the value of F77 in the makefile of L-BFGS-B:


#FC = f77
FC = g77

Then running


make

produced i386 executables and object files.

Note: I’m not sure what you would need to change to produce x86_64 files, luckily for me I wanted i386 all along.

Update: So the examples that g77 creates above work fine. But when I tried to call the fortran subroutines from c or cpp code I got weird errors. To solve this I switched gears and opted to use f2c/f77 and universal binaries, this so far is working much better.