calling a FORTRAN subroutine that accepts a function returning a vector/matrix

[John W. Eaton]

On 18-Oct-2005, John Weatherwax wrote:

| [question about writing functions to call Fortran code that in turn
| calls user-defined functions that accept matrices and vectors]

I'm attaching a modified version of the original example.  Beyond
this, I think you should look at the other examples of code for
interfacing with Fortran that are in the Octave sources.

jwe


#include <octave/oct.h>
#include <octave/parse.h>
#include <octave/f77-fcn.h>

typedef F77_RET_T (*bar_fcn_ptr) (const double*, const int&, const int&,
                                  const int&, double*, const int&);

extern "C"
{
  // Can declare args as const if we know the Fortran subroutine does
  // not change the values.  Can pass scalars by reference.

  F77_RET_T F77_FCN (bar, BAR) (bar_fcn_ptr, const double*, const int&,
                                const int&, const int&, double*);
}

static octave_function *user_fcn;

static F77_RET_T
bar_fcn (const double *a, const int& lda, const int& nr, const int& nc,
         double *r, const int& ldr)
{
  Matrix a_mat (nr, nc);
  // LDA is irrelevant here, because we know that NR == LDA.
  for (int i = 0; i < nr*nc; i++)
    a_mat(i) = a[i];

  octave_value_list args;
  args(0) = a_mat;
  int nargout = 1;

  octave_value_list retval = feval (user_fcn, args, nargout);

  if (! error_state && retval.length () == 1)
    {
      Matrix r_mat = retval(0).matrix_value ();

      if (! error_state)
        {
          if (r_mat.rows () == nr && r_mat.columns () == nc)
            {
              // LDR is irrelevant here, becase we know that NR == LDR.
              for (int i = 0; i < nr*nc; i++)
                r[i] = r_mat(i);
            }
          else
            error ("expecting result to be same size as input matrix");
        }
      else
        error ("failure evaluating user-supplied function");
    }
  else
    error ("failure evaluating user-supplied function");

  F77_RETURN (0);
}

DEFUN_DLD (foo, args, ,
  "r = foo (fcn, a)\n\
\n\
FCN is also of the form R = FCN (A)")
{
  octave_value retval;

  if (args.length () == 2)
    {
      user_fcn = args(0).function_value ();

      if (! error_state)
        {
          Matrix a = args(1).matrix_value ();

          if (! error_state)
            {
              int nr = a.rows ();
              int nc = a.columns ();
              int lda = nr;

              Matrix r (nr, nc);

              F77_FCN (bar, BAR) (bar_fcn, a.data (), lda, nr, nc,
                                  r.fortran_vec ());

              if (! error_state)
                retval = r;
            }
          else
            error ("foo: expecting scalar for second argument");
        }
      else
        error ("foo: expecting function handle for first argument");
    }
  else
    print_usage ("foo");

  return retval;
}

      subroutine bar (fcn, a, lda, nr, nc, r, ldr)
      external fcn
      integer lda, nr, nc, ldr
      double precision a(lda,*), r(ldr,*)
      call fcn (a, lda, nr, nc, r)
      return
      end