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 the possibility of Octave crashing due to errors in the user code.  For example a segmentation violation in the user's 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

          [file helloworld.cc not included]

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, int16NDArray.h, etc.

uint8NDArray
uint16NDArray
uint32NDArray
uint64NDArray
     8, 16, 32 and 64-bit unsigned array classes defined in uint8NDArray.h, uint16NDArray.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

 -- Method: T& operator () (octave_idx_type)
 -- Method: 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 functions do significant error checking and so in some circumstances the user might prefer to access the data of the array or matrix directly through the fortran_vec method discussed below.

 -- Method: octave_idx_type nelem (void) const
     The total number of elements in the matrix or array.

 -- Method: size_t byte_size (void) const
     The number of bytes used to store the matrix or array.

 -- Method: dim_vector dims (void) const
     The dimensions of the matrix or array in value of type dim_vector.

 -- Method: 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.

 -- Method: T* fortran_vec (void)
     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

          [file addtwomatrices.cc not included]

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

Character Strings in Oct-Files
------------------------------

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

          [file stringdemo.cc not included]

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;
     ...
     // Create single quoted string
     retval(1) = octave_value (ch, true, '\'');

     // Create a double quoted string
     retval(0) = octave_value (ch, true);

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

Octave's cell type is equally accessible within 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

          [file celldemo.cc not included]

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

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

          [file structdemo.cc not included]

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 iterate 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 object.  We can equally iterate of the elements of the Cell array to address the elements of the structure array.

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 similarity 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 larger than the representation of 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 Sparse Matrices in Oct-Files
.....................................

You have several alternatives for creating a sparse matrix. You 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, you can create a sparse matrix with the appropriate amount of space and then fill in the values.  Both techniques have their advantages and disadvantages.

Here is an example of how to create a small sparse matrix with the first technique

     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 of Sparse Matrices.  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 possible 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 valid.  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 possible to create a sparse matrix with too much storage, so having NZ above equaling 6 is also valid. The disadvantage is that the matrix occupies more memory than strictly needed.

It is not always easy to know 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.

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

          [file globaldemo.cc not included]

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.

          [file funcdemo.cc not included]

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 object.  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 XSTOPX, 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 XSTOPX 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

          [file fortdemo.cc not included]

and the fortran function is

          [file fortsub.f not included]

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
-------------------------------------

As oct-files are compiled functions they have the possibility of causing Octave to abort abnormally.  It is therefore important that each and every function has the minimum of parameter checking needed to ensure that Octave behaves well.

The minimum requirement, as previously discussed, is to check the number of input arguments before using them to avoid referencing a non existent argument.  However, it some case this might not be sufficient as the underlying code imposes further constraints.  For example an external function call might be undefined if the input arguments are not integers, or if one of the arguments is zero.  Therefore, oct-files often need additional input parameter checking.

There are several functions within Octave that might be useful for the purposes of parameter checking.  These include the methods of the octave_value class like is_real_matrix, etc, but equally include more specialized functions.  Some of the more common ones are demonstrated in the following example

          [file paramdemo.cc not included]

and an example of its use is

     paramdemo ([1, 2, NaN, Inf])
     => Properties of input array:
           includes Inf or NaN values
           includes other values than 1 and 0
           includes only int, Inf or NaN values

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

          [file unwinddemo.cc not included]

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
-----------------------------------

The documentation of an oct-file is the fourth string parameter of the DEFUN_DLD macro.  This string can be formatted in the same manner as the help strings for user functions (Documentation Tips), however there are some issue that are particular to the formatting of help strings within oct-files.

The major issue is that the help string will typically be longer than a single line of text, and so the formatting of long help strings need to be taken into account.  There are several manners in which to treat this issue, but the most common is illustrated in the following example

     DEFUN_DLD (do_what_i_want, args, nargout,
       "-*- texinfo -*-\n\
     @deftypefn {Function File} {} do_what_i_say (@var{n})\n\
     A function that does what the user actually wants rather\n\
     than what they requested.\n\
     @end deftypefn")
     {
     ...
     }

where, as can be seen, end line of text within the help string is terminated by \n\ which is an embedded new-line in the string together with a C++ string continuation character.  Note that the final \ must be the last character on the line.

Octave also includes the ability to embed the test and demonstration code for a function within the code itself (Test and Demo Functions). This can be used from within oct-files (or in fact any file) with certain provisos.  Firstly, the test and demo functions of Octave look for a %! as the first characters on a new-line to identify test and demonstration code.  This is equally a requirement for oct-files.  Furthermore the test and demonstration code must be included in a comment block of the compiled code to avoid it being interpreted by the compiler.  Finally, the Octave test and demonstration code must have access to the source code of the oct-file and not just the compiled code as the tests are stripped from the compiled code.  An example in an oct-file might be

     /*

     %!error (sin())
     %!error (sin(1,1))
     %!assert (sin([1,2]),[sin(1),sin(2)])

     */

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 internal API of an alternative product 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 command line.  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

          [file firstmexdemo.c not included]

This simple example demonstrates the basics of writing a mex-file.  The entry point into the mex-file is defined by mexFunction.  Note that the function name is not explicitly included in the mexFunction and so there can only be a single mexFunction entry point per-file.  Also the name of the function is determined by the name of the mex-file itself.  Therefore if the above function is in the file firstmexdemo.c, it can be compiled with

     mkoctfile --mex firstmexdemo.c

which creates a file firstmexdemo.mex.  The function can then be run from Octave as

     firstmexdemo()
     => 1.2346

It should be noted that the mex-file contains no help string for the functions it contains.  To document mex-files, there should exist an m-file in the same directory as the mex-file itself.  Taking the above as an example, we would therefore have a file firstmexdemo.m that might contain the text

     %FIRSTMEXDEMO Simple test of the functionality of a mex-file.

In this case, the function that will be executed within Octave will be given by the mex-file, while the help string will come from the m-file.  This can also be useful to allow a sample implementation of the mex-file within the Octave language itself for testing purposes.

Although we cannot have multiple entry points into a single mex-file, we can use the mexFunctionName function to determine what name the mex-file was called with.  This can be used to alter the behavior of the mex-file based on the function name.  For example if

          [file myfunc.c not included]

is in file myfunc.c, and it is compiled with

     mkoctfile --mex myfunc.c
     ln -s myfunc.mex myfunc2.mex

Then as can be seen by

     myfunc()
     => You called function: myfunc
         This is the principal function
     myfunc2()
     => You called function: myfunc2

the behavior of the mex-file can be altered depending on the functions name.

Allow the user should only include mex.h in their code, Octave declares additional functions, typedefs, etc, available to the user to write mex-files in the headers mexproto.h and mxarray.h.

Working with Matrices and Arrays in Mex-Files
---------------------------------------------

The basic mex type of all variables is mxArray.  All variables, such as matrices, cell arrays or structures are all stored in this basic type, and this type serves basically the same purpose as the octave_value class in oct-files.  That is it acts as a container for the more specialized types.

The mxArray structure contains at a minimum, the variable it represents name, its dimensions, its type and whether the variable is real or complex.  It can however contain a number of additional fields depending on the type of the mxArray.  There are a number of functions to create mxArray structures, including mxCreateCellArray, mxCreateSparse and the generic mxCreateNumericArray.

The basic functions to access the data contained in an array is mxGetPr.  As the mex interface assumes that the real and imaginary parts of a complex array are stored separately, there is an equivalent function mxGetPi that get the imaginary part.  Both of these functions are for use only with double precision matrices.  There also exists the generic function mxGetData and mxGetImagData that perform the same operation on all matrix types.  For example

     mxArray *m;
     mwSize *dims;
     UINT32_T *pr;

     dims = (mwSize *) mxMalloc (2 * sizeof(mwSize));
     dims[0] = 2;
     dims[1] = 2;
     m = mxCreateNumericArray (2, dims, mxUINT32_CLASS, mxREAL);
     pr =  = (UINT32_T *) mxGetData (m);

There are also the functions mxSetPr, etc, that perform the inverse, and set the data of an Array to use the block of memory pointed to by the argument of mxSetPr.

Note the type mwSize used above, and mwIndex are defined as the native precision of the indexing in Octave on the platform on which the mex-file is built. This allows both 32- and 64-bit platforms to support mex-files. mwSize is used to define array dimension and maximum number or elements, while mwIndex is used to define indexing into arrays.

An example that demonstration how to work with arbitrary real or complex double precision arrays is given by the file mypow2.c as given below.

          [file mypow2.c not included]

with an example of its use

     b = randn(4,1) + 1i * randn(4,1);
     all(b.^2 == mypow2(b))
     => 1

The example above uses the functions mxGetDimensions, mxGetNumberOfElements, and mxGetNumberOfDimensions to work with the dimensions of multi-dimensional arrays.  The functions mxGetM, and mxGetN are also available to find the number of rows and columns in a matrix.

Character Strings in Mex-Files
------------------------------

As mex-files do not make the distinction between single and double quoted strings within Octave, there is perhaps less complexity in the use of strings and character matrices in mex-files.  An example of their use, that parallels the demo in stringdemo.cc, is given in the file mystring.c, as seen below.

          [file mystring.c not included]

An example of its expected output is

     mystring(["First String"; "Second String"])
     => s1 = Second String
             First String

Other functions in the mex interface for handling character strings are mxCreateString, mxArrayToString, and mxCreateCharMatrixFromStrings.  In a mex-file, a character string is considered to be a vector rather than a matrix.  This is perhaps an arbitrary distinction as the data in the mxArray for the matrix is consecutive in any case.

Cell Arrays with Mex-Files
--------------------------

We can perform exactly the same operations in Cell arrays in mex-files as we can in oct-files.  An example that reduplicates the functional of the celldemo.cc oct-file in a mex-file is given by mycell.c as below

          [file mycell.c not included]

which as can be seen below has exactly the same behavior as the oct-file version.

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

        1   2

     b3 = test

Note in the example the use of the mxDuplicateArry function.  This is needed as the mxArray pointer returned by mxGetCell might be deallocated.  The inverse function to mxGetCell is mcSetCell and is defined as

     void mxSetCell (mxArray *ptr, int idx, mxArray *val);

Finally, to create a cell array or matrix, the appropriate functions are

     mxArray *mxCreateCellArray (int ndims, const int *dims);
     mxArray *mxCreateCellMatrix (int m, int n);

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

The basic function to create a structure in a mex-file is mxCreateStructMatrix, which creates a structure array with a two dimensional matrix, or mxCreateStructArray.

     mxArray *mxCreateStructArray (int ndims, int *dims,
                                   int num_keys,
                                   const char **keys);
     mxArray *mxCreateStructMatrix (int rows, int cols,
                                    int num_keys,
                                    const char **keys);

Accessing the fields of the structure can then be performed with the mxGetField and mxSetField or alternatively with the mxGetFieldByNumber and mxSetFieldByNumber functions.

     mxArray *mxGetField (const mxArray *ptr, mwIndex index,
                          const char *key);
     mxArray *mxGetFieldByNumber (const mxArray *ptr,
                                  mwIndex index, int key_num);
     void mxSetField (mxArray *ptr, mwIndex index,
                      const char *key, mxArray *val);
     void mxSetFieldByNumber (mxArray *ptr, mwIndex index,
                              int key_num, mxArray *val);

A difference between the oct-file interface to structures and the mex-file version is that the functions to operate on structures in mex-files directly include an index over the elements of the arrays of elements per field.  Whereas the oct-file structure includes a Cell Array per field of the structure.

An example that demonstrates the use of structures in mex-file can be found in the file mystruct.c, as seen below

          [file mystruct.c not included]

An example of the behavior of this function within Octave is then

     a(1).f1 = "f11"; a(1).f2 = "f12";
     a(2).f1 = "f21"; a(2).f2 = "f22";
     b = mystruct(a)
     => field f1(0) = f11
         field f1(1) = f21
         field f2(0) = f12
         field f2(1) = f22
         b =
         {
           this =

           (,
             [1] = this1
             [2] = this2
             [3] = this3
             [4] = this4
           ,)

           that =

           (,
             [1] = that1
             [2] = that2
             [3] = that3
             [4] = that4
           ,)

         }

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

The Octave format for sparse matrices is identical to the mex format in that it is a compressed column sparse format.  Also in both, sparse matrices are required to be two dimensional.  The only difference is that the real and imaginary parts of the matrix are stored separately.

The mex-file interface, as well as using mxGetM, mxGetN, mxSetM, mxSetN, mxGetPr, mxGetPi, mxSetPr and mxSetPi, the mex-file interface supplies the functions

     mwIndex *mxGetIr (const mxArray *ptr);
     mwIndex *mxGetJc (const mxArray *ptr);
     mwSize mxGetNzmax (const mxArray *ptr);

     void mxSetIr (mxArray *ptr, mwIndex *ir);
     void mxSetJc (mxArray *ptr, mwIndex *jc);
     void mxSetNzmax (mxArray *ptr, mwSize nzmax);

mxGetNzmax gets the maximum number of elements that can be stored in the sparse matrix.  This is not necessarily the number of non-zero elements in the sparse matrix.  mxGetJc returns an array with one additional value than the number of columns in the sparse matrix.  The difference between consecutive values of the array returned by mxGetJc define the number of non-zero elements in each column of the sparse matrix.  Therefore

     mwSize nz, n;
     mwIndex *Jc;
     mxArray *m;
     ...
     n = mxGetN (m);
     Jc = mxGetJc (m);
     nz = Jc[n];

returns the actual number of non-zero elements stored in the matrix in nz.  As the arrays returned by mxGetPr and mxGetPi only contain the non-zero values of the matrix, we also need a pointer to the rows of the non-zero elements, and this is given by mxGetIr.  A complete example of the use of sparse matrices in mex-files is given by the file mysparse.c as seen below

          [file mysparse.c not included]

Calling Other Functions in Mex-Files
------------------------------------

It is also possible call other Octave functions from within a mex-file using mexCallMATLAB.  An example of the use of mexCallMATLAB can be see in the example below

          [file myfeval.c not included]

If this code is in the file myfeval.c, and is compiled to myfeval.mex, then an example of its use is

     myfeval("sin", 1)
     a = myfeval("sin", 1)
     => Hello, World!
         I have 2 inputs and 1 outputs
         I'm going to call the interpreter function sin
         a =  0.84147

Note that it is not possible to use function handles or inline functions within a mex-file.

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.

          [file standalone.cc not included]

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

     $ mkoctfile --link-stand-alone standalone.cc -o standalone
     $ ./standalone
     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. The above allows the Octave math libraries to be used by an application. It does not however allow the script files, oct-files or builtin functions of Octave to be used by the application. To do that the Octave interpreter needs to be initialized first. An example of how to do this can then be seen in the code

          [file embedded.cc not included]

which is compiled and run as before as a standalone application with

     $ mkoctfile --link-stand-alone embedded.cc -o embedded
     $ ./embedded
     GCD of [10, 15] is 5
     $