[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
calling a FORTRAN subroutine that accepts a function returning a vector/
From: 
John Weatherwax 
Subject: 
calling a FORTRAN subroutine that accepts a function returning a vector/matrix 
Date: 
Tue, 18 Oct 2005 20:49:59 0500 
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...)
ReturnPath: <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 (68232188117.pittpa.adelphia.net
[68.232.188.117])
(authenticated bits=0)
by portkey.cae.wisc.edu (8.13.4/8.13.4/Debian1) with ESMTP id
j4NI70GQ032169
(version=TLSv1/SSLv3 cipher=DHERSAAES256SHA bits=256 verify=NOT);
Mon, 23 May 2005 13:07:02 0500
From: "John W. Eaton" <address@hidden>
MIMEVersion: 1.0
ContentType: multipart/mixed; boundary="m/VBRpPSbI"
ContentTransferEncoding: 7bit
MessageID: <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
InReplyTo: <address@hidden>
References: <address@hidden>
XMailer: VM 7.19 under Emacs 21.4.1
XCAEMailScannerInformation: Please contact address@hidden if this
message contains a virus or has been corrupted in delivery.
XCAEMailScanner: Found to be clean (starburst)
m/VBRpPSbI
ContentType: text/plain; charset=usascii
ContentDescription: message body and .signature
ContentTransferEncoding: 7bit
On 23May2005, 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.2246e16
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
ContentType: text/plain
ContentDisposition: inline;
filename="foo.cc"
ContentTransferEncoding: 7bit
#include <octave/oct.h>
#include <octave/parse.h>
#include <octave/f77fcn.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 usersupplied function");
}
else
error ("failure evaluating usersupplied 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
ContentType: text/plain
ContentDisposition: inline;
filename="bar.f"
ContentTransferEncoding: 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

 calling a FORTRAN subroutine that accepts a function returning a vector/matrix,
John Weatherwax <=