
/* Calling GNU Fortran subroutines and functions from Pure works pretty much
   like calling C routines. You just have to remember the following:

   * Subroutine identifiers have an extra trailing underscore.

   * By default, parameters are passed per reference. You can use the 'value'
     attribute in the Fortran code to override this.

   * Matrices are traversed in column-major order, so the Fortran subroutine
     will actually work on the transpose of your Pure matrices.

   The Fortran source code for the following examples can be found in
   fortran.f90. This first needs to be compiled to a shared library as
   follows:

     gfortran -shared -o fortran.so fortran.f90

   Add -fPIC on 64 bit systems. Or just run 'make' in this directory. */

using "lib:fortran";

/* A simple scalar function. */

extern int fact_(int) = fact;

fact 10;
map fact (1..10);

/* Matrix example (Gaussian elimination). Pure matrices can be passed as void*
   parameters (meaning that the raw data pointer will be passed), enabling the
   Fortran routine to modify the matrix in-place. */

extern void gauss_(int m, int n, void *A, void *index);

/* Here we use a little wrapper routine which takes care of allocating the
   result matrix. */

gauss A::matrix = index, A when
  // Get the dimensions and pack the matrix into contiguous storage. This also
  // allocates a new matrix which can be modified in-place.
  n,m = dim A; A = pack A;
  // Create the index vector to be passed to the Fortran routine.
  index = imatrix n;
  // Call the Fortran routine to do all the hard work. Note that in this
  // example the Fortran routine is prepared to work on the columns of the
  // input matrix, so that we don't have to transpose inputs and outputs here.
  gauss_ m n A index;
  // Translate to zero-based indices.
  index = map pred index;
  // Apply the row permutation to the result matrix, in order to bring the
  // matrix into triangle form. (You may want to comment this out to just
  // return the result as is.)
  A = rowcat [row A (index!i) | i = 0..n-1];
end if dmatrixp A;

/* For convenience, print a double matrix in "short" format a la Octave. */

using system;
__show__ x::matrix
= strcat [printd j (x!(i,j))|i=0..n-1; j=0..m-1] + "\n"
with printd 0 = sprintf "\n%10.5f"; printd _ = sprintf "%10.5f" end
when n,m = dim x end if dmatrixp x;

/* Example: */

let x = dmatrix {2,1,-1,8; -3,-1,2,-11; -2,1,2,-3};
x; gauss x;
