help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Calling a octave function from a fortran function that accepts a functio


From: John W. Eaton
Subject: Calling a octave function from a fortran function that accepts a function as an argument
Date: Mon, 23 May 2005 14:02:26 -0400

On 23-May-2005, John Weatherwax wrote:

|     I am having trouble with "Calling a octave function from a fortran 
| function that accepts a function as an argument"
| 
|     Many FORTRAN routines (think quadrature or ODE's) require a function 
| to be given at the same time as the  computational routine .

Yes, of course.  Octave already contains examples of such functions.

| Does 
| anyone have a SIMPLE example of how to create an *.oct file that will 
| allow an octave function to be called from the FORTRAN computational 
| routine? 
| 
| For instance in octave (testfun.m):
| 
| function [ res ] = testfun( x )
|   res = sin( x ) ;
| endfunction
| 
| In FORTRAN (fortsub2.f):
| 
|       subroutine fortsub2(f,xin,nxin,xout)
| C     Call f(.) on every element of xin:
| C     return in xout:
|       external f
|       real*8 xin(*),xout(*)
|       integer*4 nxin
|       do i=1,nin
|          xout(i)=f(xin(i))
|       enddo
|       return
|       end
| 
| Now in C++ I have been unsuccessfull in creating the *.cc file to make 
| an *.oct file with mkoctfile.  I want to call the routine with something 
| like:
| 
|  > xin = rand( 4,1 );
|  > [ xout ] = exOct2FortFnCall( "testfun", xin );
| 
| The result in xout should be xout = sin( xin );
| 
| My feable attempt at this procedure is included below.  Any help 
| pointers would be greatly appreciated.  I know there are source code 
| examples for things like LSODE but I was hoping for a smaller example 
| that I could play with.

OK, I try the code below.  It is about as simple as I can make it
without throwing out useful things like checking arguments to avoid
crashing if the functions are called incorrectly.  You should be able
to save the attachments and run

  mkoctfile foo.cc bar.f

then run Octave and you should be able to do

  octave:1> foo (@sin, pi)
  ans =  1.2246e-16
  octave:2> foo (@cos, pi)
  ans = -1
  octave:3> foo (inline ("x^2"), 2)
  ans = 4
  octave:4> foo (inline ("sin(x)"), pi/4)
  ans = 0.70711

For this simplest of examples, you must pass the function argument
using a function handle or as an inline function.  It will take a bit
more work to also allow name of the function to be passed as a
character string.

jwe

-- 


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

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

extern "C"
{
  F77_RET_T F77_FCN (bar, BAR) (bar_fcn_ptr, const double&, double&);
}

static octave_function *user_fcn;

static F77_RET_T
bar_fcn (const double& x, double& r)
{
  octave_value_list args;
  args(0) = x;
  int nargout = 1;

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

  if (! error_state && retval.length () == 1)
    {
      r = retval(0).double_value ();

      if (error_state)
        error ("failure evaluating user-supplied function");
    }
  else
    error ("failure evaluating user-supplied function");

  F77_RETURN (0);
}

DEFUN_DLD (foo, args, ,
  "r = foo (fcn, x)")
{
  octave_value retval;

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

      if (! error_state)
        {
          double r;

          double x = args(1).double_value ();

          if (! error_state)
            {
              F77_FCN (bar, BAR) (bar_fcn, x, r);

              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, x, r)
      external fcn
      double precision x, r
      call fcn (x, r)
      return
      end

reply via email to

[Prev in Thread] Current Thread [Next in Thread]