From help-request at octave dot org Wed Oct 19 09:33:01 2005
Subject: calling a FORTRAN subroutine that accepts a function returning a vector/matrix
From: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
To: "John Weatherwax" <wax0530 at hotmail dot com>
Cc: help at octave dot org
Date: Wed, 19 Oct 2005 10:32:39 -0400


--q1dxSwy21d
Content-Type: text/plain; charset=us-ascii
Content-Description: message body text
Content-Transfer-Encoding: 7bit

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



--q1dxSwy21d
Content-Type: text/plain
Content-Disposition: inline;
	filename="foo.cc"
Content-Transfer-Encoding: 7bit

#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;
}

--q1dxSwy21d
Content-Type: text/plain
Content-Disposition: inline;
	filename="bar.f"
Content-Transfer-Encoding: 7bit

      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

--q1dxSwy21d--



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------