Dynamically Linked Functions
****************************

Octave has the possibility of including compiled code as dynamically linked extensions and then using these extensions as if they were part of Octave itself.  Octave has the option of directly calling C++ code through its native oct-file interface or C code through its mex interface.  It can also indirectly call functions written in any other language through a simple wrapper.  The reasons to write code in a compiled language might be either to link to an existing piece of code and allow it to be used within Octave, or to allow improved performance for key pieces of code.

Before going further, you should first determine if you really need to use dynamically linked functions at all.  Before proceeding with writing any dynamically linked function to improve performance you should address ask yourself

  * Can I get the same functionality using the Octave scripting language only.

  * Is it thoroughly optimized Octave code? Vectorization of Octave code, doesn't just make it concise, it generally significantly improves its performance.  Above all, if loops must be used, make sure that the allocation of space for variables takes place outside the loops using an assignment to a like matrix or zeros.

  * Does it make as much use as possible of existing built-in library routines? These are highly optimized and many do not carry the overhead of being interpreted.

  * Does writing a dynamically linked function represent useful investment of your time, relative to staying in Octave?

Also, as oct- and mex-files are dynamically linked to octave, they introduce to possibility of having Octave abort due to coding errors in the user code.  For example a segmentation violation in the users code will cause Octave to abort.

Oct-Files
=========

Getting Started with Oct-Files
------------------------------

The basic command to build oct-files is mkoctfile and it can be call from within octave or from the command line.

(mkoctfile)

Consider the short example

     #include <octave/oct.h>

     DEFUN_DLD (helloworld, args, nargout,
       "Hello World Help String")
     {
       int nargin = args.length ();
       octave_stdout << "Hello World has " << nargin
             << " input arguments and "
             << nargout << " output arguments.\n";
       return octave_value_list ();
     }

This example although short introduces the basics of writing a C++ function that can be dynamically linked to Octave.  The easiest way to make available most of the definitions that might be necessary for an oct-file in Octave is to use the #include <octave/oct.h> header.

The macro that defines the entry point into the dynamically loaded function is DEFUN_DLD.  This macro takes four arguments, these being

 1. The function name as it will be seen in Octave,

 2. The list of arguments to the function of type octave_value_list,

 3. The number of output arguments, which can and often is omitted if not used, and

 4. The string that will be seen as the help text of the function.

The return type of functions defined with DEFUN_DLD is always octave_value_list.

There are a couple of important considerations in the choice of function name.  Firstly, it must be a valid Octave function name and so must be a sequence of letters, digits and underscores, not starting with a digit.  Secondly, as Octave uses the function name to define the filename it attempts to find the function in, the function name in the DEFUN_DLD macro must match the filename of the oct-file.  Therefore, the above function should be in a file helloworld.cc, and it should be compiled to an oct-file using the command

     mkoctfile helloworld.cc

This will create a file call helloworld.oct, that is the compiled version of the function.  It should be noted that it is perfectly acceptable to have more than one DEFUN_DLD function in a source file.  However, there must either be a symbolic link to the oct-file for each of the functions defined in the source code with the DEFUN_DLD macro or the autoload (Function Files) function should be used.

The rest of this function then shows how to find the number of input arguments, how to print through the octave pager, and return from the function.  After compiling this function as above, an example of its use is

     helloworld(1,2,3)
     => Hello World has 3 input arguments and 0 output arguments.

Matrices and Arrays in Oct-Files
--------------------------------

Octave supports a number of different array and matrix classes, the majority of which are based on the Array class.  The exception is the sparse matrix types discussed separately below.  There are three basic matrix types

Matrix
     A double precision matrix class defined in dMatrix.h,

ComplexMatrix
     A complex matrix class defined in CMatrix.h, and

BoolMatrix
     A boolean matrix class defined in boolMatrix.h.

These are the basic two-dimensional matrix types of octave.  In additional there are a number of multi-dimensional array types, these being

NDArray
     A double precision array class defined in dNDArray.h,

ComplexNDarray
     A complex array class defined in CNDArray.h,

boolNDArray
     A boolean array class defined in boolNDArray.h,

int8NDArray, int16NDArray, int32NDArray, int64NDArray
     8, 16, 32 and 64-bit signed array classes defined in int8NDArray.h, etc, and

uint8NDArray, uint16NDArray, uint32NDArray, uint64NDArray
     8, 16, 32 and 64-bit unsigned array classes defined in uint8NDArray.h, etc.

There are several basic means of constructing matrices of multi-dimensional arrays.  Considering the Matrix type as an example

  * We can create an empty matrix or array with the empty constructor.  For example

          Matrix a;

     This can be used on all matrix and array types

  * Define the dimensions of the matrix or array with a dim_vector.  For example

          dim_vector dv(2);
          dv(0) = 2; dv(1) = 2;
          Matrix a(dv);

     This can be used on all matrix and array types

  * Define the number of rows and columns in the Matrix.  For example

          Matrix a(2,2)

     However, this constructor can only be used with the matrix types.

These types all share a number of basic methods and operators, a selection of which include

T& operator (octave_idx_type), T& elem(octave_idx_type)
     The () operator or elem method allow the values of the matrix or array to be read or set.  These can take a single argument, which is of type octave_idx_type, that is the index into the matrix or array.  Additionally, the matrix type allows two argument versions of the () operator and elem method, giving the row and column index of the value to obtain or set.

     Note that these function do significant error checking and so in some circumstances the user might prefer the access the data of the array or matrix directly through the fortran_vec method discussed below.

octave_idx_type nelem ()
     The total number of elements in the matrix or array.

size_t byte_size ()
     The number of bytes used to store the matrix or array.

dim_vector dims()
     The dimensions of the matrix or array in value of type dim_vector.

void resize (const dim_vector&)
     A method taking either an argument of type dim_vector, or in the case of a matrix two arguments of type octave_idx_type defining the number of rows and columns in the matrix.

T* fortran_vec ()
     This method returns a pointer to the underlying data of the matrix or a array so that it can be manipulated directly, either within Octave or by an external library.

Operators such an +, -, or * can be used on the majority of the above types.  In addition there are a number of methods that are of interest only for matrices such as transpose, hermitian, solve, etc.

The typical way to extract a matrix or array from the input arguments of DEFUN_DLD function is as follows

     #include <octave/oct.h>

     DEFUN_DLD (addtwomatrices, args, , "Add A to B")
     {
       int nargin = args.length ();
       if (nargin != 2)
         print_usage ();
       else
         {
           NDArray A = args(0).array_value();
           NDArray B = args(1).array_value();
           if (! error_state)
             return octave_value(A + B);
         }
       return octave_value_list ();
     }

To avoid segmentation faults causing Octave to abort, this function explicitly checks that there are sufficient arguments available before accessing these arguments.  It then obtains two multi-dimensional arrays of type NDArray and adds these together.  Note that the array_value method is called without using the is_matrix_type type, and instead the error_state is checked before returning A + B.  The reason to prefer this is that the arguments might be a type that is not an NDArray, but it would make sense to convert it to one.  The array_value method allows this conversion to be performed transparently if possible, and sets error_state if it is not.

A + B, operating on two NDArray's returns an NDArray, which is cast to an octave_value on the return from the function.  An example of the use of this demonstration function is

     addtwomatrices(ones(2,2),ones(2,2))
           =>  2  2
               2  2

A list of the basic Matrix and Array types, the methods to extract these from an octave_value and the associated header is listed below.

RowVector   row_vector_value   dRowVector.h
ComplexRowVector   complex_row_vector_value   CRowVector.h
ColumnVector   column_vector_value   dColVector.h
ComplexColumnVector   complex_column_vector_value   CColVector.h
Matrix   matrix_value   dMatrix.h
ComplexMatrix   complex_matrix_value   CMatrix.h
boolMatrix   bool_matrix_value   boolMatrix.h
charMatrix   char_matrix_value   chMatrix.h
NDArray   array_value   dNDArray.h
ComplexNDArray   complex_array_value   CNDArray.h
boolNDArray   bool_array_value   boolNDArray.h
charNDArray   char_array_value   charNDArray.h
int8NDArray   int8_array_value   int8NDArray.h
int16NDArray   int16_array_value   int16NDArray.h
int32NDArray   int32_array_value   int32NDArray.h
int64NDArray   int64_array_value   int64NDArray.h
uint8NDArray   uint8_array_value   uint8NDArray.h
uint16NDArray   uint16_array_value   uint16NDArray.h
uint32NDArray   uint32_array_value   uint32NDArray.h
uint64NDArray   uint64_array_value   uint64NDArray.h

Using Sparse Matrices in Oct-Files
----------------------------------

There are three classes of sparse objects that are of interest to the user.

SparseMatrix
     A double precision sparse matrix class

SparseComplexMatrix
     A complex sparse matrix class

SparseBoolMatrix
     A boolean sparse matrix class

All of these classes inherit from the Sparse<T> template class, and so all have similar capabilities and usage.  The Sparse<T> class was based on Octave Array<T> class, and so users familiar with Octave's Array classes will be comfortable with the use of the sparse classes.

The sparse classes will not be entirely described in this section, due to their similar with the existing Array classes.  However, there are a few differences due the different nature of sparse objects, and these will be described.  Firstly, although it is fundamentally possible to have N-dimensional sparse objects, the Octave sparse classes do not allow them at this time.  So all operations of the sparse classes must be 2-dimensional.  This means that in fact SparseMatrix is similar to Octave's Matrix class rather than its NDArray class.

The Differences between the Array and Sparse Classes
....................................................

The number of elements in a sparse matrix is considered to be the number of non-zero elements rather than the product of the dimensions.  Therefore

       SparseMatrix sm;
       ...
       int nel = sm.nelem ();

returns the number of non-zero elements.  If the user really requires the number of elements in the matrix, including the non-zero elements, they should use numel rather than nelem.  Note that for very large matrices, where the product of the two dimensions is large that the representation of the an unsigned int, then numel can overflow. An example is speye(1e6) which will create a matrix with a million rows and columns, but only a million non-zero elements.  Therefore the number of rows by the number of columns in this case is more than two hundred times the maximum value that can be represented by an unsigned int. The use of numel should therefore be avoided useless it is known it won't overflow.

Extreme care must be take with the elem method and the "()" operator, which perform basically the same function.  The reason is that if a sparse object is non-const, then Octave will assume that a request for a zero element in a sparse matrix is in fact a request to create this element so it can be filled.  Therefore a piece of code like

       SparseMatrix sm;
       ...
       for (int j = 0; j < nc; j++)
         for (int i = 0; i < nr; i++)
           std::cerr << " (" << i << "," << j << "): " << sm(i,j)
                     << std::endl;

is a great way of turning the sparse matrix into a dense one, and a very slow way at that since it reallocates the sparse object at each zero element in the matrix.

An easy way of preventing the above from happening is to create a temporary constant version of the sparse matrix.  Note that only the container for the sparse matrix will be copied, while the actual representation of the data will be shared between the two versions of the sparse matrix.  So this is not a costly operation.  For example, the above would become

       SparseMatrix sm;
       ...
       const SparseMatrix tmp (sm);
       for (int j = 0; j < nc; j++)
         for (int i = 0; i < nr; i++)
           std::cerr << " (" << i << "," << j << "): " << tmp(i,j)
                     << std::endl;

Finally, as the sparse types aren't just represented as a contiguous block of memory, the fortran_vec method of the Array<T> is not available.  It is however replaced by three separate methods ridx, cidx and data, that access the raw compressed column format that the Octave sparse matrices are stored in. Additionally, these methods can be used in a manner similar to elem, to allow the matrix to be accessed or filled.  However, in that case it is up to the user to respect the sparse matrix compressed column format discussed previous.

Creating Spare Matrices in Oct-Files
....................................

The user has several alternatives in how to create a sparse matrix. They can first create the data as three vectors representing the row and column indexes and the data, and from those create the matrix. Or alternatively, they can create a sparse matrix with the appropriate amount of space and then fill in the values.  Both techniques have their advantages and disadvantages.

An example of how to create a small sparse matrix with the first technique might be seen the example

       int nz = 4, nr = 3, nc = 4;
       ColumnVector ridx (nz);
       ColumnVector cidx (nz);
       ColumnVector data (nz);

       ridx(0) = 0; ridx(1) = 0; ridx(2) = 1; ridx(3) = 2;
       cidx(0) = 0; cidx(1) = 1; cidx(2) = 3; cidx(3) = 3;
       data(0) = 1; data(1) = 2; data(2) = 3; data(3) = 4;

       SparseMatrix sm (data, ridx, cidx, nr, nc);

which creates the matrix given in section Storage.  Note that the compressed matrix format is not used at the time of the creation of the matrix itself, however it is used internally.

As previously mentioned, the values of the sparse matrix are stored in increasing column-major ordering.  Although the data passed by the user does not need to respect this requirement, the pre-sorting the data significantly speeds up the creation of the sparse matrix.

The disadvantage of this technique of creating a sparse matrix is that there is a brief time where two copies of the data exists.  Therefore for extremely memory constrained problems this might not be the right technique to create the sparse matrix.

The alternative is to first create the sparse matrix with the desired number of non-zero elements and then later fill those elements in.  The easiest way to do this is

       int nz = 4, nr = 3, nc = 4;
       SparseMatrix sm (nr, nc, nz);
       sm(0,0) = 1; sm(0,1) = 2; sm(1,3) = 3; sm(2,3) = 4;

That creates the same matrix as previously.  Again, although it is not strictly necessary, it is significantly faster if the sparse matrix is created in this manner that the elements are added in column-major ordering.  The reason for this is that if the elements are inserted at the end of the current list of known elements then no element in the matrix needs to be moved to allow the new element to be inserted.  Only the column indexes need to be updated.

There are a few further points to note about this technique of creating a sparse matrix.  Firstly, it is not illegal to create a sparse matrix with fewer elements than are actually inserted in the matrix.  Therefore

       int nz = 4, nr = 3, nc = 4;
       SparseMatrix sm (nr, nc, 0);
       sm(0,0) = 1; sm(0,1) = 2; sm(1,3) = 3; sm(2,3) = 4;

is perfectly legal.  However it is a very bad idea.  The reason is that as each new element is added to the sparse matrix the space allocated to it is increased by reallocating the memory.  This is an expensive operation, that will significantly slow this means of creating a sparse matrix.  Furthermore, it is not illegal to create a sparse matrix with too much storage, so having NZ above equaling 6 is also legal. The disadvantage is that the matrix occupies more memory than strictly needed.

It is not always easy to known the number of non-zero elements prior to filling a matrix.  For this reason the additional storage for the sparse matrix can be removed after its creation with the maybe_compress function.  Furthermore, the maybe_compress can deallocate the unused storage, but it can equally remove zero elements from the matrix.  The removal of zero elements from the matrix is controlled by setting the argument of the maybe_compress function to be 'true'.  However, the cost of removing the zeros is high because it implies resorting the elements.  Therefore, if possible it is better is the user doesn't add the zeros in the first place.  An example of the use of maybe_compress is

       int nz = 6, nr = 3, nc = 4;
       SparseMatrix sm1 (nr, nc, nz);
       sm1(0,0) = 1; sm1(0,1) = 2; sm1(1,3) = 3; sm1(2,3) = 4;
       sm1.maybe_compress ();  // No zero elements were added

       SparseMatrix sm2 (nr, nc, nz);
       sm2(0,0) = 1; sm2(0,1) = 2; sm(0,2) = 0; sm(1,2) = 0;
       sm1(1,3) = 3; sm1(2,3) = 4;
       sm2.maybe_compress (true);  // Zero elements were added

The use of the maybe_compress function should be avoided if possible, as it will slow the creation of the matrices.

A third means of creating a sparse matrix is to work directly with the data in compressed row format.  An example of this technique might be

       octave_value arg;

       ...

       int nz = 6, nr = 3, nc = 4;   // Assume we know the max no nz
       SparseMatrix sm (nr, nc, nz);
       Matrix m = arg.matrix_value ();

       int ii = 0;
       sm.cidx (0) = 0;
       for (int j = 1; j < nc; j++)
         {
           for (int i = 0; i < nr; i++)
             {
               double tmp = foo (m(i,j));
               if (tmp != 0.)
                 {
                   sm.data(ii) = tmp;
                   sm.ridx(ii) = i;
                   ii++;
                 }
             }
           sm.cidx(j+1) = ii;
        }
       sm.maybe_compress ();  // If don't know a-priori the final no of nz.

which is probably the most efficient means of creating the sparse matrix.

Finally, it might sometimes arise that the amount of storage initially created is insufficient to completely store the sparse matrix.  Therefore, the method change_capacity exists to reallocate the sparse memory. The above example would then be modified as

       octave_value arg;

       ...

       int nz = 6, nr = 3, nc = 4;   // Assume we know the max no nz
       SparseMatrix sm (nr, nc, nz);
       Matrix m = arg.matrix_value ();

       int ii = 0;
       sm.cidx (0) = 0;
       for (int j = 1; j < nc; j++)
         {
           for (int i = 0; i < nr; i++)
             {
               double tmp = foo (m(i,j));
               if (tmp != 0.)
                 {
                   if (ii == nz)
                     {
                       nz += 2;   // Add 2 more elements
                       sm.change_capacity (nz);
                     }
                   sm.data(ii) = tmp;
                   sm.ridx(ii) = i;
                   ii++;
                 }
             }
           sm.cidx(j+1) = ii;
        }
       sm.maybe_mutate ();  // If don't know a-priori the final no of nz.

Note that both increasing and decreasing the number of non-zero elements in a sparse matrix is expensive, as it involves memory reallocation.  Also as parts of the matrix, though not its entirety, exist as the old and new copy at the same time, additional memory is needed.  Therefore if possible this should be avoided.

Using Sparse Matrices in Oct-Files
..................................

Most of the same operators and functions on sparse matrices that are available from the Octave are equally available with oct-files. The basic means of extracting a sparse matrix from an octave_value and returning them as an octave_value, can be seen in the following example

        octave_value_list retval;

        SparseMatrix sm = args(0).sparse_matrix_value ();
        SparseComplexMatrix scm = args(1).sparse_complex_matrix_value ();
        SparseBoolMatrix sbm = args(2).sparse_bool_matrix_value ();

        ...

        retval(2) = sbm;
        retval(1) = scm;
        retval(0) = sm;

The conversion to an octave-value is handled by the sparse octave_value constructors, and so no special care is needed.

Using Strings in Oct-Files
--------------------------

In Octave a string is just a special Array class.  Consider the example

     #include <octave/oct.h>

     DEFUN_DLD (stringdemo, args, , "String Demo")
     {
       int nargin = args.length();
       octave_value_list retval;

       if (nargin != 1)
         print_usage ();
       else
         {
           charMatrix ch = args(0).char_matrix_value ();

           if (! error_state)
             {
               if (args(0).is_sq_string ())
                 retval(1) = octave_value (ch, true);
               else
                 retval(1) = octave_value (ch, true, '\'');

               octave_idx_type nr = ch.rows();
               for (octave_idx_type i = 0; i < nr / 2; i++)
                 {
                   std::string tmp = ch.row_as_string (i);
                   ch.insert (ch.row_as_string(nr - i - 1).c_str(), i, 0);
                   ch.insert (tmp.c_str(), nr - i - 1, 0);
                 }
               retval(0) = octave_value (ch, true);
             }
         }
       return retval;
     }

An example of the of the use of this function is

     s0 = ["First String";"Second String"];
     [s1,s2] = stringdemo (s0)
     =>s1 = Second String
             First String

     =>s2 = First String
             Second String

     typeinfo (s2)
     => sq_string
     typeinfo (s1)
     => string

One additional complication of strings in Octave is the difference between single quoted and double quoted strings.  To find out if an octave_value contains a single or double quoted string an example is

         if (args(0).is_sq_string())
           octave_stdout << "First argument is a singularly quoted string\n";
         else if (args(0).is_dq_string())
           octave_stdout << "First argument is a doubly quoted string\n";

Note however, that both types of strings are represented by the charNDArray type, and so when assigning to an octave_value, the type of string should be specified.  For example

     octave_value_list retval;
     charNDArray c;
     ...
     retval(0) = octave_value (ch, true); // Create a double quoted string
     retval(1) = octave_value (ch, true, "'"); // Create single quoted string

Cell Arrays in Oct-Files
------------------------

Octave's cell type is equally accessible within an oct-files.  A cell array is just an array of octave_values, and so each element of the cell array can then be treated just like any other octave_value.  A simple example is

     #include <octave/oct.h>
     #include <octave/Cell.h>

     DEFUN_DLD (celldemo, args, , "Cell Demo")
     {
       octave_value_list retval;
       int nargin = args.length();

       if (nargin != 1)
         print_usage ();
       else
         {
           Cell c = args (0).cell_value ();
           if (! error_state)
             for (octave_idx_type i = 0; i < c.nelem (); i++)
               retval(i) = c.elem (i);
         }
       return retval;
     }

Note that cell arrays are used less often in standard oct-files and so the Cell.h header file must be explicitly included.  The rest of this example extracts the octave_values one by one from the cell array and returns be as individual return arguments.  For example consider

     [b1,b2,b3] = celldemo({1, [1,2], "test"})
     =>
     b1 =  1
     b2 =

        1   2

     b3 = test

Using Structures in Oct-Files
-----------------------------

A structure in Octave is map between a number of fields represented and their values.  The "Standard Template Library" map class is used, with the pair consisting of a std::string and an octave Cell variable.

A simple example demonstrating the use of structures within oct-files is

     #include <octave/oct.h>
     #include <octave/ov-struct.h>

     DEFUN_DLD (structdemo, args, , "Struct demo.")
     {
       int nargin = args.length ();
       octave_value retval;

       if (nargin != 2)
         print_usage ();
       else
         {
           Octave_map arg0 = args(0).map_value ();
           std::string arg1 = args(1).string_value ();

           if (! error_state && arg0.contains (arg1))
             {
               // The following two lines might be written as
               //    octave_value tmp;
               //    for (Octave_map::iterator p0 = arg0.begin() ;
               //        p0 != arg0.end(); p0++ )
               //      if (arg0.key (p0) == arg1)
               //        {
               //          tmp = arg0.contents (p0) (0);
               //          break;
               //        }
               // though using seek is more concise.
               Octave_map::const_iterator p1 = arg0.seek (arg1);
               octave_value tmp =  arg0.contents( p1 ) (0);
               Octave_map st;
               st.assign ("selected", tmp);
               retval = octave_value (st);
             }
         }
       return retval;
     }

An example of its use is

     x.a = 1; x.b = "test"; x.c = [1,2];
     structdemo(x,"b")
     => selected = test

The commented code above demonstrates how to iterated over all of the fields of the structure, where as the following code demonstrates finding a particular field in a more concise manner.

As can be seen the contents method of the Octave_map class returns a Cell which allows Structure Arrays to be represented.  Therefore, to obtain the underlying octave_value we write

     octave_value tmp = arg0.contents (p1) (0);

where the trailing (0) is the () operator on the Cell array.

Accessing Global Variables in Oct-Files
---------------------------------------

Global variables allow variables in the global scope to be accessed.  Global variables can easily be accessed with oct-files using the support functions get_global_value and set_global_value.  get_global_value takes two arguments, the first is a string representing the variable name to obtain.  The second argument is a boolean argument specifying what to do in the case that no global variable of the desired name is found.  An example of the use of these two functions is

     #include <octave/oct.h>

     DEFUN_DLD (globaldemo, args, , "Global demo.")
     {
       int nargin = args.length();
       octave_value retval;

       if (nargin != 1)
         print_usage ();
       else
         {
           std::string s = args(0).string_value ();
           if (! error_state)
             {
               octave_value tmp = get_global_value (s, true);
               if (tmp.is_defined ())
                 retval = tmp;
               else
                 retval = "Global variable not found";

               set_global_value ("a", 42.0);
             }
         }
       return retval;
     }

An example of its use is

     global a b
     b = 10;
     globaldemo ("b")
     => 10
     globaldemo ("c")
     => "Global variable not found"
     num2str (a)
     => 42

Calling Octave Functions from Oct-Files
---------------------------------------

There is often a need to be able to call another octave function from within an oct-file, and there are many examples of such within octave itself.  For example the quad function is an oct-file that calculates the definite integral by quadrature over a user supplied function.

There are also many ways in which a function might be passed.  It might be passed as one of

 1. Function Handle

 2. Anonymous Function Handle

 3. Inline Function

 4. String

The example below demonstrates an example that accepts all four means of passing a function to an oct-file.

     #include <octave/oct.h>
     #include <octave/parse.h>

     DEFUN_DLD (funcdemo, args, nargout, "Function Demo")
     {
       int nargin = args.length();
       octave_value_list retval;

       if (nargin < 2)
         print_usage ();
       else
         {
           octave_value_list newargs;
           for (octave_idx_type i = nargin - 1; i > 0; i--)
             newargs (i - 1) = args(i);
           if (args(0).is_function_handle () ||
               args(0).is_inline_function ())
             {
               octave_function *fcn = args(0).function_value ();
               if (! error_state)
                 retval = feval (fcn, newargs, nargout);
             }
           else if (args(0).is_string ())
             {
               std::string fcn = args (0).string_value ();
               if (! error_state)
                 retval = feval (fcn, newargs, nargout);
             }
           else
             error ("funcdemo: expected string, inline or function handle");
         }
       return retval;
     }

The first argument to this demonstration is the user supplied function and the following arguments are all passed to the user function.

     funcdemo(@sin,1)
     => 0.84147
     funcdemo(@(x) sin(x), 1)
     => 0.84147
     funcdemo (inline("sin(x)"), 1)
     => 0.84147
     funcdemo("sin",1)
     => 0.84147
     funcdemo (@atan2, 1, 1)
     => 0.78540

When the user function is passed as a string, the treatment of the function is different.  In some cases it is necessary to always have the user supplied function as an octave_function.  In that case the string argument can be used to create a temporary function like

       std::octave fcn_name = unique_symbol_name ("__fcn__");
       std::string fname = "function y = ";
       fname.append (fcn_name);
       fname.append ("(x) y = ");
       fcn = extract_function (args(0), "funcdemo", fcn_name,
                               fname, "; endfunction");
       ...
       if (fcn_name.length())
         clear_function (fcn_name);

There are two important things to know in this case.  The number of input arguments to the user function is fixed, and in the above is a single argument, and secondly to avoid leaving the temporary function in the Octave symbol table it should be cleared after use.

Calling External Code from Oct-Files
------------------------------------

Linking external C code to Octave is relatively simple, as the C functions can easily be called directly from C++.  One possible issue is the declarations of the external C functions might need to be explicitly defined as C functions to the compiler.  If the declarations of the external C functions are in the header foo.h, then the manner in which to ensure that the C++ compiler treats these declarations as C code is

     #ifdef __cplusplus
     extern "C"
     {
     #endif
     #include "foo.h"
     #ifdef __cplusplus
     }  /* end extern "C" */
     #endif

Calling Fortran code however can pose some difficulties.  This is due to differences in the manner in compilers treat the linking of Fortran code with C or C++ code.  Octave supplies a number of macros that allow consistent behavior across a number of compilers.

The underlying Fortran code should use the XSTOPX function to replace the Fortran STOP function.  XSTOPX uses the Octave exception handler to treat failing cases in the fortran code explicitly.  Note that Octave supplies its own replacement blas XERBLA function, which uses XSTOPX.

If the underlying code calls XSTOP, then the F77_XFCN macro should be used to call the underlying fortran function.  The Fortran exception state can then be checked with the global variable f77_exception_encountered.  If XSTOP will not be called, then the F77_FCN macro should be used instead to call the Fortran code.

There is no harm in using F77_XFCN in all cases, except that for Fortran code that is short running and executes a large number of times, there is potentially an overhead in doing so.  However, if F77_FCN is used with code that calls XSTOP, Octave can generate a segmentation fault.

An example of the inclusion of a Fortran function in an oct-file is given in the following example, where the C++ wrapper is

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

     extern "C"
     {
       F77_RET_T
       F77_FUNC (fortsub, FORTSUB)
             (const int&, double*, F77_CHAR_ARG_DECL
              F77_CHAR_ARG_LEN_DECL);
     }

     DEFUN_DLD (fortdemo , args , , "Fortran Demo.")
     {
       octave_value_list retval;
       int nargin = args.length();
       if (nargin != 1)
         print_usage ();
       else
         {
           NDArray a = args(0).array_value ();
           if (! error_state)
             {
               double *av = a.fortran_vec ();
               octave_idx_type na = a.nelem ();
               OCTAVE_LOCAL_BUFFER (char, ctmp, 128);

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

               if ( f77_exception_encountered )
                 error ("fortdemo: error in fortran");
               else
                 {
                   retval(1) = std::string (ctmp);
                   retval(0) = a;
                 }
             }
         }
       return retval;
     }

and the fortran function is

           subroutine fortsub (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('fortsub:divide by zero')
              else
                 a (i) = 1d0 / a (i)
              end if
           end do
           write(unit = s, fmt = '(a,i3,a,a)', iostat = ioerr)
          1     'There are ', n, ' values in the input vector', char(0)
           if (ioerr .ne. 0) then
              call xstopx('fortsub: error writing string')
           end if
           return
           end

This example demonstrates most of the features needed to link to an external Fortran function, including passing arrays and strings, as well as exception handling.  An example of the behavior of this function is

     [b, s] = fortdemo(1:3)
     =>
       b = 1.00000   0.50000   0.33333
       s = There are   3 values in the input vector
     [b, s] = fortdemo(0:3)
     error: fortsub:divide by zero
     error: exception encountered in Fortran subroutine fortsub_
     error: fortdemo: error in fortran

Allocating Local Memory in Oct-Files
------------------------------------

Allocating memory within an oct-file might seem easy as the C++ new/delete operators can be used.  However, in that case care must be taken to avoid memory leaks.  The preferred manner in which to allocate memory for use locally is to use the OCTAVE_LOCAL_BUFFER macro.  An example of its use is

     OCTAVE_LOCAL_BUFFER (double, tmp, len)

that returns a pointer tmp of type double * of length len.

Input Parameter Checking in Oct-Files
-------------------------------------

WRITE ME

Exception and Error Handling in Oct-Files
-----------------------------------------

Another important feature of Octave is its ability to react to the user typing Control-C even during calculations.  This ability is based on the C++ exception handler, where memory allocated by the C++ new/delete methods are automatically released when the exception is treated.  When writing an oct-file, to allow Octave to treat the user typing Control-C, the OCTAVE_QUIT macro is supplied.  For example

       for (octave_idx_type i = 0; i < a.nelem (); i++)
         {
           OCTAVE_QUIT;
           b.elem(i) = 2. * a.elem(i);
         }

The presence of the OCTAVE_QUIT macro in the inner loop allows Octave to treat the user request with the Control-C.  Without this macro, the user must either wait for the function to return before the interrupt is processed, or press Control-C three times to force Octave to exit.

The OCTAVE_QUIT macro does impose a very small speed penalty, and so for loops that are known to be small it might not make sense to include OCTAVE_QUIT.

When creating an oct-file that uses an external libraries, the function might spend a significant portion of its time in the external library.  It is not generally possible to use the OCTAVE_QUIT macro in this case.  The alternative in this case is

       BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
       ...  some code that calls a "foreign" function ...
       END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;

The disadvantage of this is that if the foreign code allocates any memory internally, then this memory might be lost during an interrupt, without being deallocated.  Therefore, ideally Octave itself should allocate any memory that is needed by the foreign code, with either the fortran_vec method or the OCTAVE_LOCAL_BUFFER macro.

The Octave unwind_protect mechanism (The unwind_protect Statement) can also be used in oct-files.  In conjunction with the exception handling of Octave, it is important to enforce that certain code is run to allow variables, etc to be restored even if an exception occurs.  An example of the use of this mechanism is

     #include <octave/oct.h>
     #include <octave/unwind-prot.h>

     void
     err_hand (const char *fmt, ...)
     {
       // Do nothing!!
     }

     DEFUN_DLD (unwinddemo, args, nargout, "Unwind Demo")
     {
       int nargin = args.length();
       octave_value retval;
       if (nargin < 2)
         print_usage ();
       else
         {
           NDArray a = args(0).array_value ();
           NDArray b = args(1).array_value ();

           if (! error_state)
             {
               unwind_protect::begin_frame("Funwinddemo");
               unwind_protect_ptr(current_liboctave_warning_handler);
               set_liboctave_warning_handler(err_hand);
               retval = octave_value ( quotient (a, b));
               unwind_protect::run_frame("Funwinddemo");
             }
         }
       return retval;
     }

As can be seen in the example

     unwinddemo(1,0)
     => Inf
     1 / 0
     => warning: division by zero
         Inf

The division by zero (and in fact all warnings) is disabled in the unwinddemo function.

Documentation and Test of Oct-Files
-----------------------------------

WRITE ME, reference how to use texinfo in oct-file and embed test code.

Application Programming Interface for Oct-Files
-----------------------------------------------

WRITE ME, using Coda section 1.3 as a starting point.

Mex-Files
=========

Octave includes an interface to allow legacy mex-files to be compiled and used with Octave.  This interface can also be used to share code between Octave and non Octave users.  However, as mex-files expose the intern API of a product alternative to Octave, and the internal structure of Octave is different to this product, a mex-file can never have the same performance in Octave as the equivalent oct-file.  In particular to support the manner in which mex-files access the variables passed to mex functions, there are a significant number of additional copies of memory when calling or returning from a mex function.  For this reason, new code should be written using the oct-file interface discussed above if possible.

Getting Started with Mex-Files
------------------------------

The basic command to build a mex-file is either mkoctfile --mex or mex.  The first can either be used from within Octave or from the commandline.  However, to avoid issues with the installation of other products, the use of the command mex is limited to within Octave.

(mex)

(mexext)

One important difference between the use of mex with other products and with Octave is that the header file "matrix.h" is implicitly included through the inclusion of "mex.h".  This is to avoid a conflict with the Octave file "Matrix.h" with operating systems and compilers that don't distinguish between filenames in upper and lower case

Consider the short example

     #include "mex.h"

     void
     mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     {
       mxArray *v = mxCreateDoubleMatrix (1, 1, mxREAL);
       double *data = mxGetPr (v);
       *data = 1.23456789;
       plhs[0] = v;
     }

This simple example demonstrates the basics of writing a mex-file.

WRITE ME

Using Structures with Mex-Files
-------------------------------

WRITE ME

Sparse Matrices with Mex-Files
------------------------------

WRITE ME

Calling External Functions in Mex-Files
---------------------------------------

WRITE ME

Standalone Programs
===================

The libraries Octave itself uses, can be utilized in standalone applications.  These applications then have access, for example, to the array and matrix classes as well as to all the Octave algorithms.  The following C++ program, uses class Matrix from liboctave.a or liboctave.so.

     #include <iostream>
     #include <octave/oct.h>
     int
     main (void)
     {
       std::cout << "Hello Octave world!\n";
       const int size = 2;
       Matrix a_matrix = Matrix(size, size);
       for (octave_idx_type row = 0; row < size; ++row)
         {
           for (octave_idx_type column = 0; column < size; ++column)
             {
               a_matrix(row, column) = (row + 1)*10 + (column + 1);
             }
         }
       std::cout << a_matrix;
       return 0;
     }

mkoctfile can then be used to build a standalone application with a command like

     $ mkoctfile --link-stand-alone hello.cc -o hello
     $ ./hello
     Hello Octave world!
       11 12
       21 22
     $

Note that the application hello will be dynamically linked against the octave libraries and any octave support libraries.