octave c++函数中调用fortran77子程序

2019独角兽企业重金招聘Python工程师标准>>>

代码都放在工作目录~/octave_workplace/f77下

一、fortran子程序tnine.f

      SUBROUTINE TNINE (IOPT, PARMOD, PS, X, Y, Z, BX, BY, BZ) 
      INTEGER IOPT 
      DOUBLE PRECISION PARMOD(10), PS, X, Y, Z, BX, BY, BZ 

C     This is just a test subroutine body, to check connexions. 
C     Put the sum of PARMOD in PS, and X, Y, Z into BX, BY, BZ 

      INTEGER I 
      
      PS = 0D0 
      DO 1 I=1, 10 
         PS = PS + PARMOD (I) 
 1    CONTINUE 

      BX = X 
      BY = Y 
      BZ = Z 

      END 

这里参考了octave自带例子程序examples/fortransub.f

二、c++程序t96.cc

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

extern "C" 
{ 
  int F77_FUNC (tnine, TNINE) (const int& IOPT, const double* PARMOD, 
                                 double& PS, 
                                 const double& X, const double& Y, 
                                 const double &Z, 
                                 double& BX, double& BY, double& BZ ); 
} 

DEFUN_DLD (t96, args, , 
           "- Loadable Function: [PS, BX, BY, BZ] = t96 (PM, X, Y, Z) Returns the sum of PM in PS and X, Y, and Z in BX, BY, and BZ.") 
{ 
  octave_value_list retval; 

  const int dummy_integer = 0; 
  Matrix pm; 
  const double x = args(1).double_value(), y = args(2).double_value(), 
    z = args(3).double_value(); 
  double ps, bx, by, bz; 

  pm = args(0).matrix_value (); 

  F77_XFCN (tnine, TNINE, 
            (dummy_integer, pm.fortran_vec(), ps, x, y, z, bx, by, bz) ); 

  if (f77_exception_encountered) 
    { 
      error ("unrecoverable error in t96"); 
      return retval; 
    } 

  retval(0) = octave_value (ps); 
  retval(1) = octave_value (bx); 
  retval(2) = octave_value (by); 
  retval(3) = octave_value (bz); 
  return retval; 
} 

三、Compile this  in the Bourne Again Shell（can also in Octave） and run it in Octave like: 

>> mkoctfile t96.cc tnine.f
>>  [p, x, y, z] = t96 (1:10, sqrt (2), pi, e)
p =  55
x =  1.4142
y =  3.1416
z =  2.7183
>>

四、官方例子：

fortransub.f

      subroutine fortransub (n, a, s)
      implicit none
      character*(*) s
      real*8 a(*)
      integer*4 i, n, ioerr
      do i = 1, n
        if (a(i) .eq. 0d0) then
          call xstopx ('fortransub: divide by zero')
        else
          a(i) = 1d0 / a(i)
        endif
      enddo
      write (unit = s, fmt = '(a,i3,a,a)', iostat = ioerr)
     $       'There are ', n,
     $       ' values in the input vector', char(0)
      if (ioerr .ne. 0) then
        call xstopx ('fortransub: error writing string')
      endif
      return
      end

fortrandemo.cc

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

extern "C"
{
  F77_RET_T
  F77_FUNC (fortransub, FORTSUB)
    (const F77_INT&, F77_DBLE*, F77_CHAR_ARG_DECL F77_CHAR_ARG_LEN_DECL);
}

DEFUN_DLD (fortrandemo, args,  /* nargout */, "Fortran Demo")
{
  if (args.length () != 1)
    print_usage ();

  NDArray a = args(0).array_value ();

  double *av = a.fortran_vec ();
  octave_idx_type na = a.numel ();

  OCTAVE_LOCAL_BUFFER (char, ctmp, 128);

  F77_XFCN (fortransub, FORTSUB,
            (na, av, ctmp F77_CHAR_ARG_LEN (128)));

  return ovl (a, std::string (ctmp));
}

编译：

>> mkoctfile fortrandemo.cc fortransub.f

转载于:https://my.oschina.net/u/2245781/blog/1815474