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

[John Weatherwax]

Hello,

I have been trying to get some legacy FORTRAN partial differential
equations software running in octave.  Much like the ODE solvers that
already exist, this requires writing an octave wrapper to extract a
pointer to a user defined function (a function that returns the
residual derivatives given the inputs).  This function would then be
used inside a octave subroutine that the FORTRAN code would repeatedly
call.  John Eaton was kind enough to provide the code below which
works very well for scalars, i.e. when the function FORTRAN expects is
something like

SUBROUTINE F(T,X,DXDT)

with everything a scalar.  In my PDE problem however, the FORTRAN
routine I need to call has the following form (I was perhaps unclear
on this in an earlier post)

SUBROUTINE F( T, X, U, UX, UXX, FVAL, NPDE )
DIMENSION U(NPDE), UX(NPDE), UXX(NPDE), FVAL(NPDE)

Here the INPUTS: T and X are double precision scalars, U,UX,UXX are
double precision vectors of length NPDE, and the output FVAL is vector
of length NPDE.  For the past couple of days I have tried to modify
the given code to return a vector type argument but have been unable
to get the pointers to work out correctly with the octave C++ API.
The idea would be to have the user supply a PDE in the octave language
like follows:

function [ res ] = F(t,x,u,ux,uxx)

res = [ u(2)*u(2)*uxx(1) - u(1)*u(2) - u(1)^2 + 10     + 2*u(2)*ux(2)*ux(1);
        u(1)*u(1)*uxx(2) + u(1)*u(2) - u(2)^2 + uxx(1) + 2*u(1)*ux(1)*ux(2); ];

and have the FORTRAN code repeatedly call this function (as "F" above)
returning the vector "res".  I would appreciate it if some kind person
with more expertise than I do would offer suggestions as to how to
modify the code below to accomplish what I would like to do.  I could
send the code I attempted but since there are a great number of other
details it might be better to hold off on flooding this group with code.

A great thanks to any who respond,

John Weatherwax

PS.  FYI: at other places in the legacy FORTRAN we must call a
function that returns a MATRIX like

SUBROUTINE BNDRY( T, X, U, UX, DBDU, DBDUX, DZDT, NPDE )
DIMENSION U(NPDE), UX(NPDE), DZDT(NPDE)
DIMENSION DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE)

Here the scalar inputs are T and X, the vector inputs are U and UX and
the MATRIX outputs are DBDU, DBDUX.  So ideally any solutions could
cover this case also ... (I am hopefull...)







Return-Path: <address@hidden>
Received: from ll.mit.edu (llpost [155.34.234.40])
        by llinfo.llan.ll.mit.edu (8.12.10+Sun/8.12.6) with ESMTP id 
j4NI8GiH020932
        for <address@hidden>; Mon, 23 May 2005 14:08:17 -0400 (EDT)
Received: (from address@hidden)
        by ll.mit.edu (8.12.10/8.8.8) id j4NI8GC7025394
        for <address@hidden>; Mon, 23 May 2005 14:08:16 -0400 (EDT)
Received: from UNKNOWN(144.92.13.24), claiming to be "mail.cae.wisc.edu"
via SMTP by llmail, id smtpdAAAt2aiOV; Mon May 23 14:07:42 2005
Received: from portkey.cae.wisc.edu (portkey.cae.wisc.edu [144.92.13.118])
        by mail.cae.wisc.edu (8.12.11/8.12.11) with ESMTP id j4NI72Gg013170;
        Mon, 23 May 2005 13:07:03 -0500 (CDT)
Received: from devzero.bogus.domain (68-232-188-117.pittpa.adelphia.net 
[68.232.188.117])
        (authenticated bits=0)
	by portkey.cae.wisc.edu (8.13.4/8.13.4/Debian-1) with ESMTP id 
j4NI70GQ032169
        (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NOT);
        Mon, 23 May 2005 13:07:02 -0500
From: "John W. Eaton" <address@hidden>
MIME-Version: 1.0
Content-Type: multipart/mixed; boundary="m/VBRpPSbI"
Content-Transfer-Encoding: 7bit
Message-ID: <address@hidden>
Date: Mon, 23 May 2005 14:02:26 -0400
To: John Weatherwax <address@hidden>
Cc: address@hidden
Subject: Calling a octave function from a fortran function that accepts a
function as an argument
In-Reply-To: <address@hidden>
References: <address@hidden>
X-Mailer: VM 7.19 under Emacs 21.4.1
X-CAE-MailScanner-Information: Please contact address@hidden if this 
message contains a virus or has been corrupted in delivery.
X-CAE-MailScanner: Found to be clean (starburst)


--m/VBRpPSbI
Content-Type: text/plain; charset=us-ascii
Content-Description: message body and .signature
Content-Transfer-Encoding: 7bit

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

--



--m/VBRpPSbI
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&, 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;
}

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

     subroutine bar (fcn, x, r)
     external fcn
     double precision x, r
     call fcn (x, r)
     return
     end

--m/VBRpPSbI--

_________________________________________________________________
Is your PC infected? Get a FREE online computer virus scan from McAfee® 
Security. http://clinic.mcafee.com/clinic/ibuy/campaign.asp?cid=3963



-------------------------------------------------------------
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
-------------------------------------------------------------