Re: Re: Re: recursion and fsolve

[Paul THOMAS]

David,

There was a break in the experiment last night so I wrote a recursive fortran routine to calculate the factorial of a number. I did not have time to explore further. I tried both the gfortran binary from http://f77.linksysnet.com and the gfortran in the gcc-3.5.0 I built myself, with identical results.

#include <octave/oct.h>

#include <octave/f77-fcn.h>

extern "C" int F77_FUNC (g95sub,G95SUB) ( double&, double&);

DEFUN_DLD (g95demo1, args , , "y=g95demo1(x)\

returns the factorial of x for integer values")

{

double dinptr = args(0).scalar_value();

double doutptr;

F77_FUNC(g95sub,G95SUB)(dinptr,doutptr );

return octave_value(doutptr);

}

recursive subroutine g95sub(din,dout)

real*8 :: din ,dout

if ( din > 1.0 ) then

call g95sub( din - 1.0 , dout )

dout = din * dout

else

dout = 1.0

end if

return

end subroutine g95sub

I had problems with the library libgfortran, in trying to build octave, so I separated the fortran compile from mkoctfile:

> /gcc35/gfortran -c g95sub.f95 -o g95sub.o

> mkoctfile g95demo1.cc g95sub.o -L/gcc35/lib -lgfortran

et voila, g95demo1(x) gives you the factorial of x! Note that I have not tried recursion thru' the octave interface but I cannot see any reason why it should not work. The compiler objected to using assumed length arrays in the arguments. However, I have not checked if the fortran95 standard permits that in recursive routines. You should be able to render the fortran in fsolve re-entrant by declaring all the routines to be recursive, putting "end subroutine mysubroutinenname" and changing .f to .f95

Best regards

Paul

> Message du 06/08/04 09:58
> De : "David Bateman"
> A : "Paul THOMAS"
> Copie à : "David Bateman" , "John W. Eaton" , address@hidden, "Fred Metcalf"
> Objet : Re: Re: recursion and fsolve
> Daprès Paul THOMAS (le 05/08/2004):
> > Dear All,
> >
> > I have been testing gfortran, the fortran compiler for gcc3.5. Unfortunately,
> > gfortran and g95 forked about a year and a half ago and it is the former that
> > is associated with gcc.
> >
> > I have mainly been examining the fortran95 aspects but have also taken a
> > look-see at what it does with fortran77. I have only found one broken bit that
> > affects octave - multiple entry points are not supported. I rewrote the two
> > subroutines in ranlib that contain ENTRY and compiled every last bit of fortran
> > in octave.
>
> If multiple entry points only occur in randlib my preference would be to
> get the mersenne twister and zigurat code from octave-forge in to octave
> and eliminate randlib entirely. Making changes to code in libcruft is a
> bit problematic due to the fact that we might want to update versions
> of teh code used. Therefore if the changes you make aren't fed back to
> the original package maintainers and acceptted, we'll enter into a cycle
> where we have to modify the code anytime the libcruft libraries are
> updated.
>
> > The object modules are OK; those that I have tested work fine.
> > Unfortunately, the distribution between the libraries is such that an octave
> > build fails with gcc-3.5.0 because of multiple definitions all over the place.
> > --accept-multiple-definition allows the build to run to the end but octave then
> > falls over in kpathsea, for some reason. I suspect that the snapshot of g++ I
> > downloaded has regressed some place. I do not see any interest in
> > investigating further for the present.
>
> Ok.
>
> > However, dynamically loaded functions can be written with gfortran
> > functions and subroutines and they work fine. I haven't tried a recursive
> > function yet but will give it a try and will report back with the results.
>
> This, I'd be interested in. If it works a test can be added to configure.in
> to test if the fortran compiler supports recursion and we can enable it in
> functions like quad and fsolve.
>
> Of course that means we'd break the rule above about not modifying libcruft
> code as ddaspk.f would probably need to have the line "SAVE LID, LENID, NONNEG"
> eliminated somehow.
>
> In any case I'm interested to hear your results..
>
> Cheers
> David
>
> --
> David Bateman address@hidden
> Motorola CRM +33 1 69 35 48 04 (Ph)
> Parc Les Algorithmes, Commune de St Aubin +33 1 69 35 77 01 (Fax)
> 91193 Gif-Sur-Yvette FRANCE
>
> The information contained in this communication has been classified as:
>
> [x] General Business Information
> [ ] Motorola Internal Use Only
> [ ] Motorola Confidential Proprietary
>
>
>
> -------------------------------------------------------------
> 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
> -------------------------------------------------------------
>
>