Introducción informal a Matlab y Octave - Temas avanzados (II)

7.4  Extender Octave con otros lenguajes.

7.4.1  Una advertencia antes de empezar a programar.

1. Seguir con gcc3, escribir código en Fortran 77 y compilar con g77.
2. Pasarnos a la rama inestable de Octave, compilarla con gcc4 y gfortran y rezar para que no rompa nada.

7.4.2  Extender Octave con C++

#include <octave/oct.h>
 DEFUN_DLD (eqlorentz,args, ,
     "Ecuacion de Lorentz en C++")
     {
     ColumnVector xdot (3);
     ColumnVector x (args(0).vector_value());
     int a=10;
     int b=28;
     float c=8./3;
     xdot(0) = a*(x(1)-x(0));
     xdot(1) = x(0)*(b-x(2))-x(1);
     xdot(2) = x(0)*x(1)-c*x(2);
     
     return octave_value (xdot);
     }
 

bash ~/> mkoctfile eqlorentz.cc
 

>> help eqlorentz
 eqlorentz is the dynamically-linked function from the file
 /home/guillem/sync/CursoScripting/ejercicios/eqlorentz.oct
 Ecuacion de Lorentz en C++
 >> eqlorentz([1,1,1],1)
 ans =
     0.00000
    26.00000
    -1.66667
 

7.4.2.1  Llamar funciones desde C++

/*testfh.cc*/
 # include <octave/oct.h>
 
 DEFUN_DLD(testfh,args, ,
           "testing how C++ can call an octave function")
 {
   octave_value_list tmp;
   octave_value_list inval;
   octave_function *input_fcn=0;
   if (args(0).is_function_handle() || args(0).is_inline_function())
     input_fcn = args(0).function_value();
   else
     {
       error("this is not a function handle nor an inline function");
       return octave_value(-1);
     }
   double x = args(1).double_value();
   inval.append(octave_value(x));
   tmp = input_fcn->do_multi_index_op(1,inval);
   return tmp;
 }
 

bash~/> mkoctfile testfh.cc
 

>> testfh('sin',.123)
 error: this is not a function handle nor an inline function
 >> testfh(@sin,pi/2)
 ans = 1
 >> testfh(@(x) sin(x)*exp(x/2),pi/4)
 ans = 1.0472
 >> testfh(inline('sin(x)*exp(x/2)'),pi/4)
 ans = 1.0472
 

7.4.3  Extender Octave con Fortran

7.4.3.1  ¿Por qué Fortran?

7.4.3.2  La difícil comunicación entre C y Fortran

a=2
 

#include <stdio.h>
 int a;
 int *ptr;
 int main(int argc, char *argv[])
 {
   a=2;
   ptr =  &a;
   printf("a es %i \n",a);
   printf("la direccion de a es %p \n",ptr);
   printf("el puntero ptr contiene %i \n",*ptr);
   return 0;
 }
 

a es 2
 la direccion de a es 0x403020
 el puntero ptr contiene 2
 

#include <stdio.h>
 
 int a,b;
 int sum(int a,int b)
 { 
   int result;
   result=a+b;
   return result;
 }
 
 int main(int argc, char *argv[])
 {
   a=2;
   b=3;
   int resultado=sum(a,b);
   printf("la suma de %i y %i es %i \n",a,b,resultado);
   return 0;
 
 }
 

program test
 integer :: a,b,c
 a=2
 b=3
 call sum(a,b,c)
 write(*,*) c
 end program test
 
 
 subroutine sum(arg1,arg2,resultado)
   integer, intent(in) :: arg1,arg2
   integer, intent(out) :: resultado
   resultado=arg1+arg2
 end subroutine sum
 

7.4.3.3  Punteros y arrays

array[15];
 array[3][5];
 

7.4.3.4  Escritura de wrappers para funciones en Fortran.

subroutine eqlorentzfsub(x,xdot)
   real(8),dimension(:),intent(in) :: x
   real(8),dimension(:),intent(out) :: xdot
   real(8) :: a=10d0,b=28d0
   real(8) :: c
   c=8d0/3d0
   xdot(1)= a*(x(2)-x(1))
   xdot(2)= x(1)*(b-x(3)-x(2))
   xdot(3)= x(1)*x(2)-c*x(3)
 end subroutine eqlorentzfsub
 

#include <octave/oct.h>
 #include "f77-fcn.h"
 extern "C" int F77_FUNC (eqlorentzfsub,
                            EQLORENTZFSUB)(double*,double*);
 DEFUN_DLD(eqlorentzf,args, ,"xdot=eqlorentz(x)")
 {
   octave_value_list retval;
   ColumnVector wargin(args(0).vector_value());
   ColumnVector wargout(3);
   F77_FUNC(eqlorentzfsub,EQLORENTZFSUB)(wargin.fortran_vec(),
                                         wargout.fortran_vec());
   retval.append(wargout);
   return retval;
 }
 

1. Recogeremos la variable de entrada de la función que Octave dejará en la variable args. Si es un escalar podemos utilizar una variable propia de C, si queremos un vector o una matriz tendremos que utilizar un tipo propio de Octave. En el caso del wrapper anterior declaramos un ColumnVector llamado wargin inicializado con la variable args(0) en forma de vector.
2. Declararemos la variable de salida, en este caso otro ColumnVector, al que también asignaremos un puntero. También declararemos la variable de retorno de la función de Octave, siempre como octave_value_list.
3. Llamaremos la función previamente compilada en Fortran mediante F77_FUNC a la que pasaremos como argumentos los punteros que hemos creado para ella, en nuestro caso argin y argout. Previamente y fuera de la función DEFUN_DLD declararemos la función externa donde definiremos todos los argumentos como punteros a las variables de entrada y salida.
4. Pasaremos la variable de salida wargout a la variable retval de retorno de la función en Fortran.

bash$~/> gfortran -c eqlorentzfsub.f
 
 bash$~/> mkoctfile eqlorentzf.cc eqlorentzfsub.o
 

>> eqlorentzf([1,2,3])
 ans =
    10
    23
    -6
 

      subroutine eqlorentzfsub(x,xdot)
       real*8 x(*),xdot(*)
       real*8 a,b,c
 
       a=10d0
       b=28d0
       c=8d0/3d0
 
       xdot(1)= a*(x(2)-x(1))
       xdot(2)= x(1)*(b-x(3)-x(2))
       xdot(3)= x(1)*x(2)-c*x(3)
 
       end subroutine eqlorentzfsub
 

      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     Este es un ejemplo para probar un wrapper.
 C     Asigna la suma de PARMOD a PS, y X, Y, Z a BX, BY, BZ
       INTEGER I
       PS = 0D0
       DO I=1, 10
          PS = PS + PARMOD (I)
       END DO
       BX = X
       BY = Y
       BZ = Z
       END
 

#include <octave/oct.h>
 #include "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) \n\
  \n\
 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;
 }
 

bash $~/> mkoctfile t96.cc tnine.f
 

>> help t96
 t96 is the dynamically-linked function from the file
 /home/guillem/Desktop/sandbox/t96.oct
 - Funcion: [PS, BX, BY, BZ] = t96 (PM, X, Y, Z)
 Asigna la suma de PM a PS y  X, Y, Z a BX, BY, BZ.
 >> [uno,dos,tres,cuatro]=t96(1:10,pi,e,4)
 uno = 55
 dos = 3.1416
 tres = 2.7183
 cuatro = 4
 

subroutine f95sub2(din,dout)
 type              ::        mytype
   real*8          ::        mydouble
   integer*4       ::        myint
 end type mytype
 type (mytype)     ::        din , dout
 dout%mydouble = din%mydouble ** float( din%myint )
 dout%myint = din%myint * din%myint
 end subroutine f95sub2
 

bash $ ~/> gfortran -c f95sub.f90
 

#include <octave/oct.h>
 #include <octave/f77-fcn.h>
 struct mystruct {
   double mydouble;
   int myint;
   };
 extern "C" int F77_FUNC (f95sub,f95SUB) ( mystruct &, mystruct &); 
 DEFUN_DLD (f95demo, args , ,"[w,z] = f95demo( x , y ) \
 returns w  = x ^y and z = y * y for integer y") {
   octave_value_list retval;
   mystruct dinptr , doutptr;
   dinptr.mydouble = args(0).scalar_value();
   dinptr.myint = int( args(1).scalar_value() );
   F77_FUNC (f95sub,f95SUB) (dinptr,doutptr );
   retval.append( octave_value( doutptr.mydouble ) );
   retval.append( octave_value( double( doutptr.myint  ) ) );
   return retval;
 }
 

bash$~/>mkoctfile f95demo.cc f95sub.o