@c DO NOT EDIT!  Generated automatically by munge-texi.pl.

@c Copyright (C) 2007-2013 John W. Eaton and David Bateman
@c Copyright (C) 2007 Paul Thomas and Christoph Spiel
@c
@c This file is part of Octave.
@c
@c Octave is free software; you can redistribute it and/or modify it
@c under the terms of the GNU General Public License as published by the
@c Free Software Foundation; either version 3 of the License, or (at
@c your option) any later version.
@c 
@c Octave is distributed in the hope that it will be useful, but WITHOUT
@c ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
@c FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
@c for more details.
@c 
@c You should have received a copy of the GNU General Public License
@c along with Octave; see the file COPYING.  If not, see
@c <http://www.gnu.org/licenses/>.

@node External Code Interface
@appendix External Code Interface
@cindex dynamic-linking
@cindex Dynamically Linked Functions
@cindex Octave API

"The sum of human wisdom is not contained in any one language"
  ---Ezra Pound

Octave is a fantastic language for solving many problems in science and
engineering.  However, it is not the only computer language and there
are times when you may want to use code written in other languages.
Good reasons for doing so include: 1) not re-inventing the wheel; existing
function libraries which have been thoroughly tested and debugged or
large scale simulation codebases are a good example, 2) accessing unique
capabilities of a different language; for example the well-known regular
expression functions of Perl (but don't do that because @code{regexp}
already exists in Octave).

Performance should generally @strong{not} be a reason for using compiled
extensions.  Although compiled extensions can run faster, particularly
if they replace a loop in Octave code, this is almost never the best path
to take.  First, there are many techniques to speed up Octave performance while
remaining within the language.  Second, Octave is a high-level language that
makes it easy to perform common mathematical tasks.  Giving that up means
shifting the focus from solving the real problem to solving a computer
programming problem.  It means returning to low-level constructs such as
pointers, memory management, mathematical overflow/underflow, etc.  Because
of the low level nature, and the fact that the compiled code is executed outside
of Octave, there is the very real possibility of crashing the interpreter and
losing work.

Before going further, you should first determine if you really need to bother
writing code outside of Octave.

@itemize @bullet
@item
Can I get the same functionality using the Octave scripting language alone?

Even when a function already exists outside the language, it may be
better to simply reproduce the behavior in an m-file rather than attempt to
interface to the outside code.

@item
Is the code thoroughly optimized for Octave?

If performance is an issue you should always start with the in-language
techniques for getting better performance.  Chief among these is vectorization
(@pxref{Vectorization and Faster Code Execution}) which not only makes the
code concise and more understandable but improves performance (10X-100X).
If loops must be used, make sure that the allocation of space for variables
takes place outside the loops using an assignment to a matrix of the right
size, or zeros.

@item
Does the code make as much use as possible of existing built-in library
routines?

These routines are highly optimized and many do not carry the overhead
of being interpreted.

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

It will take time to learn Octave's interface for external code and
there will inevitably be issues with tools such as compilers.
@end itemize

With that said, Octave offers a versatile interface for including chunks
of compiled code as dynamically linked extensions.  These dynamically linked
functions can be called from the interpreter in the same manner as any
ordinary function.  The interface is bi-directional and external code can
call Octave functions (like @code{plot}) which otherwise might be very
difficult to develop.

The interface is centered around supporting the languages C++, C, and Fortran.
Octave itself is written in C++ and can call external C++/C code through its
native oct-file interface.  The C language is also supported through the
mex-file interface for compatibility with @sc{matlab}.  Fortran code is easiest
to reach through the oct-file interface.

Because many other languages provide C or C++ APIs it is relatively simple
to build bridges between Octave and other languages.  This is also a way to
bridge to hardware resources which often have device drivers written in C.

@menu
* Oct-Files::
* Mex-Files::
* Standalone Programs::
@end menu

@node Oct-Files
@section Oct-Files
@cindex oct-files
@cindex mkoctfile
@cindex oct

@menu
* Getting Started with Oct-Files::
* Matrices and Arrays in Oct-Files::
* Character Strings in Oct-Files::
* Cell Arrays in Oct-Files::
* Structures in Oct-Files::
* Sparse Matrices in Oct-Files::
* Accessing Global Variables in Oct-Files::
* Calling Octave Functions from Oct-Files::
* Calling External Code from Oct-Files::
* Allocating Local Memory in Oct-Files::
* Input Parameter Checking in Oct-Files::
* Exception and Error Handling in Oct-Files::
* Documentation and Test of Oct-Files::
@c * Application Programming Interface for Oct-Files::  
@end menu

@node Getting Started with Oct-Files
@subsection Getting Started with Oct-Files

Oct-files are pieces of C++ code that have been compiled with the Octave
API into a dynamically loadable object.  They take their name from the file
which contains the object which has the extension @file{.oct}.

Finding a C++ compiler, using the correct switches, adding the right include
paths for header files, etc. is a difficult task.  Octave automates this by
providing the @code{mkoctfile} command with which to build oct-files.  The
command is available from within Octave or at the shell command line.

@c mkoctfile scripts/miscellaneous/mkoctfile.m
@anchor{XREFmkoctfile}
@deftypefn  {Command} {} mkoctfile [-options] file @dots{}
@deftypefnx {Function File} {[@var{output}, @var{status}] =} mkoctfile (@dots{})

The @code{mkoctfile} function compiles source code written in C,
C++, or Fortran.  Depending on the options used with @code{mkoctfile}, the
compiled code can be called within Octave or can be used as a stand-alone
application.

@code{mkoctfile} can be called from the shell prompt or from the Octave
prompt.  Calling it from the Octave prompt simply delegates the
call to the shell prompt.  The output is stored in the @var{output}
variable and the exit status in the @var{status} variable.

@code{mkoctfile} accepts the following options, all of which are optional
except for the file name of the code you wish to compile:

@table @samp
@item -I DIR
Add the include directory DIR to compile commands.

@item -D DEF
Add the definition DEF to the compiler call.

@item -l LIB
Add the library LIB to the link command.

@item -L DIR
Add the library directory DIR to the link command.

@item  -M
@itemx --depend
Generate dependency files (.d) for C and C++ source files.

@item -R DIR
Add the run-time path to the link command.

@item @nospell{-Wl,@dots{}}
Pass flags though the linker like @nospell{"-Wl,-rpath=@dots{}"}.
The quotes are needed since commas are interpreted as command
separators.

@item -W@dots{}
Pass flags though the compiler like @nospell{"-Wa,OPTION"}.

@item -c
Compile but do not link.

@item -g
Enable debugging options for compilers.

@item  -o FILE
@itemx --output FILE
Output file name.  Default extension is .oct
(or .mex if @samp{--mex} is specified) unless linking
a stand-alone executable.

@item  -p VAR
@itemx --print VAR
Print the configuration variable VAR@.  Recognized variables are:

@example
   ALL_CFLAGS                  INCFLAGS
   ALL_CXXFLAGS                INCLUDEDIR
   ALL_FFLAGS                  LAPACK_LIBS
   ALL_LDFLAGS                 LD_CXX
   AR                          LDFLAGS
   BLAS_LIBS                   LD_STATIC_FLAG
   CC                          LFLAGS
   CFLAGS                      LIBDIR
   CPICFLAG                    LIBOCTAVE
   CPPFLAGS                    LIBOCTINTERP
   CXX                         LIBS
   CXXFLAGS                    OCTAVE_HOME
   CXXPICFLAG                  OCTAVE_LIBS
   DEPEND_EXTRA_SED_PATTERN    OCTAVE_LINK_DEPS
   DEPEND_FLAGS                OCTAVE_LINK_OPTS
   DL_LD                       OCTAVE_PREFIX
   DL_LDFLAGS                  OCTINCLUDEDIR
   F77                         OCTLIBDIR
   F77_INTEGER8_FLAG           OCT_LINK_DEPS
   FFLAGS                      OCT_LINK_OPTS
   FFTW3F_LDFLAGS              RANLIB
   FFTW3F_LIBS                 RDYNAMIC_FLAG
   FFTW3_LDFLAGS               READLINE_LIBS
   FFTW3_LIBS                  SED
   FFTW_LIBS                   SPECIAL_MATH_LIB
   FLIBS                       XTRA_CFLAGS
   FPICFLAG                    XTRA_CXXFLAGS
@end example

@item --link-stand-alone
Link a stand-alone executable file.

@item --mex
Assume we are creating a MEX file.  Set the default output extension
to ".mex".

@item  -s
@itemx --strip
Strip the output file.

@item  -v
@itemx --verbose
Echo commands as they are executed.

@item file
The file to compile or link.  Recognized file types are

@example
@group
   .c    C source
   .cc   C++ source
   .C    C++ source
   .cpp  C++ source
   .f    Fortran source (fixed form)
   .F    Fortran source (fixed form)
   .f90  Fortran source (free form)
   .F90  Fortran source (free form)
   .o    object file
   .a    library file
@end group
@end example

@end table
@end deftypefn


Consider the following short example which introduces the basics of
writing a C++ function that can be linked to Octave.

@example
@group
@verbatim
#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 ();
}
@end verbatim

@end group
@end example

The first critical line is @code{#include <octave/oct.h>} which 
makes available most of the definitions necessary for a C++ oct-file.
Note that @file{octave/oct.h} is a C++ header and cannot be directly
@code{#include}'ed in a C source file, nor any other language.

Included by @file{oct.h} is a definition for the macro
@w{@code{DEFUN_DLD}} which creates a dynamically loaded function.  This
macro takes four arguments:

@enumerate 1
@item The function name as it will be seen in Octave,

@item The list of arguments to the function of type @code{octave_value_list},

@item The number of output arguments, which can and often is omitted if
not used, and

@item The string to use for the help text of the function.
@end enumerate

The return type of functions defined with @w{@code{DEFUN_DLD}} is always
@code{octave_value_list}.

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

@example
mkoctfile helloworld.cc
@end example

This will create a file called @file{helloworld.oct} that is the compiled
version of the function.  It should be noted that it is perfectly
acceptable to have more than one @w{@code{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 @w{@code{DEFUN_DLD}}
macro or the @code{autoload} (@ref{Function Files}) function should be used.

The rest of the function 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

@example
@group
helloworld (1, 2, 3)
@print{} Hello World has 3 input arguments and 0 output arguments.
@end group
@end example

Subsequent sections show how to use specific classes from Octave's core
internals.  Base classes like dMatrix (a matrix of double values) are
found in the directory @file{liboctave/array}.  The definitive reference for
how to use a particular class is the header file itself.  However, it is
often enough just to study the examples in the manual in order to be able
to use the class.

@node Matrices and Arrays in Oct-Files
@subsection 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

@table @code
@item Matrix
A double precision matrix class defined in @file{dMatrix.h},

@item ComplexMatrix
A complex matrix class defined in @file{CMatrix.h}, and

@item BoolMatrix
A boolean matrix class defined in @file{boolMatrix.h}.
@end table

These are the basic two-dimensional matrix types of Octave.  In
addition there are a number of multi-dimensional array types including

@table @code
@item NDArray
A double precision array class defined in @file{dNDArray.h}

@item ComplexNDarray
A complex array class defined in @file{CNDArray.h}

@item boolNDArray
A boolean array class defined in @file{boolNDArray.h}

@item  int8NDArray
@itemx int16NDArray
@itemx int32NDArray
@itemx int64NDArray
8, 16, 32, and 64-bit signed array classes defined in
@file{int8NDArray.h}, @file{int16NDArray.h}, etc.

@item  uint8NDArray
@itemx uint16NDArray
@itemx uint32NDArray
@itemx uint64NDArray
8, 16, 32, and 64-bit unsigned array classes defined in
@file{uint8NDArray.h}, @file{uint16NDArray.h}, etc.
@end table

There are several basic ways of constructing matrices or
multi-dimensional arrays.  Using the class @code{Matrix} as an example
one can

@itemize @bullet
@item
Create an empty matrix or array with the empty constructor.  For example:

@example
Matrix a;
@end example

This can be used for all matrix and array types.

@item
Define the dimensions of the matrix or array with a dim_vector which has
the same characteristics as the vector returned from @code{size}.  For example:

@example
@group
dim_vector dv (2);
dv(0) = 2; dv(1) = 3;  // 2 rows, 3 columns
Matrix a (dv);
@end group
@end example

This can be used on all matrix and array types.

@item
Define the number of rows and columns in the matrix.  For example:

@example
Matrix a (2, 2)
@end example

However, this constructor can only be used with matrix types.
@end itemize

These types all share a number of basic methods and operators.  Many bear
a resemblance to functions that exist in the interpreter.  A selection of
useful methods include

@deftypefn  {Method} {T&} operator () (octave_idx_type)
@deftypefnx {Method} {T&} elem (octave_idx_type)
The @code{()} operator or @code{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 @code{octave_idx_type}, that is the index into the matrix or
array.  Additionally, the matrix type allows two argument versions of the
@code{()} operator and elem method, giving the row and column index of the
value to obtain or set.
@end deftypefn

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 @nospell{fortran_vec} method discussed below.

@deftypefn {Method} {} octave_idx_type numel (void) const
The total number of elements in the matrix or array.
@end deftypefn

@deftypefn {Method} {size_t} byte_size (void) const
The number of bytes used to store the matrix or array.
@end deftypefn

@deftypefn {Method} {dim_vector} dims (void) const
The dimensions of the matrix or array in value of type dim_vector.
@end deftypefn

@deftypefn {Method} {int} ndims (void) const
The number of dimensions of the matrix or array.  Matrices are 2-D,
but arrays can be N-dimensional.
@end deftypefn

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

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

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

The typical way to extract a matrix or array from the input arguments of
@w{@code{DEFUN_DLD}} function is as follows

@example
@verbatim
#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 ();
}
@end verbatim

@end example

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 @code{NDArray} and adds these together.  Note that the array_value
method is called without using the @code{is_matrix_type} type, and instead the
error_state is checked before returning @code{A + B}.  The reason to
prefer this is that the arguments might be a type that is not an
@code{NDArray}, but it would make sense to convert it to one.  The
@code{array_value} method allows this conversion to be performed
transparently if possible, and sets @code{error_state} if it is not.

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

@example
@group
addtwomatrices (ones (2, 2), eye (2, 2))
      @result{}  2  1
          1  2
@end group
@end example

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

@multitable @columnfractions .3 .4 .3
@headitem Type @tab Function @tab Source Code
@item @code{RowVector} @tab @code{row_vector_value} @tab @file{dRowVector.h}
@item @code{ComplexRowVector} @tab @code{complex_row_vector_value} @tab @file{CRowVector.h}
@item @code{ColumnVector} @tab @code{column_vector_value} @tab @file{dColVector.h}
@item @code{ComplexColumnVector} @tab @code{complex_column_vector_value} @tab @file{CColVector.h}
@item @code{Matrix} @tab @code{matrix_value} @tab @file{dMatrix.h}
@item @code{ComplexMatrix} @tab @code{complex_matrix_value} @tab @file{CMatrix.h}
@item @code{boolMatrix} @tab @code{bool_matrix_value} @tab @file{boolMatrix.h}
@item @code{charMatrix} @tab @code{char_matrix_value} @tab @file{chMatrix.h}
@item @code{NDArray} @tab @code{array_value} @tab @file{dNDArray.h}
@item @code{ComplexNDArray} @tab @code{complex_array_value} @tab @file{CNDArray.h}
@item @code{boolNDArray} @tab @code{bool_array_value} @tab @file{boolNDArray.h}
@item @code{charNDArray} @tab @code{char_array_value} @tab @file{charNDArray.h}
@item @code{int8NDArray} @tab @code{int8_array_value} @tab @file{int8NDArray.h}
@item @code{int16NDArray} @tab @code{int16_array_value} @tab @file{int16NDArray.h}
@item @code{int32NDArray} @tab @code{int32_array_value} @tab @file{int32NDArray.h}
@item @code{int64NDArray} @tab @code{int64_array_value} @tab @file{int64NDArray.h}
@item @code{uint8NDArray} @tab @code{uint8_array_value} @tab @file{uint8NDArray.h}
@item @code{uint16NDArray} @tab @code{uint16_array_value} @tab @file{uint16NDArray.h}
@item @code{uint32NDArray} @tab @code{uint32_array_value} @tab @file{uint32NDArray.h}
@item @code{uint64NDArray} @tab @code{uint64_array_value} @tab @file{uint64NDArray.h}
@end multitable

@node Character Strings in Oct-Files
@subsection Character Strings in Oct-Files

A character string in Octave is just a special @code{Array} class.
Consider the example:

@example
@verbatim
#include <octave/oct.h>

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

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

      if (! error_state)
        {
          retval(1) = octave_value (ch, '\'');  // Single Quote String

          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, '"');  // Double Quote String
        }
    }
  return retval;
}
@end verbatim

@end example

An example of the use of this function is

@example
@group
s0 = ["First String"; "Second String"];
[s1,s2] = stringdemo (s0)
@result{} s1 = Second String
        First String

@result{} s2 = First String
        Second String

typeinfo (s2)
@result{} sq_string
typeinfo (s1)
@result{} string
@end group
@end example

One additional complication of strings in Octave is the difference
between single quoted and double quoted strings.  To find out if an
@code{octave_value} contains a single or double quoted string use
one of the predicate tests shown below.

@example
@group
if (args(0).is_sq_string ())
  octave_stdout << "First argument is a single quoted string\n";
else if (args(0).is_dq_string ())
  octave_stdout << "First argument is a double quoted string\n";
@end group
@end example

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

@example
@group
octave_value_list retval;
charNDArray ch;
@dots{}
// Create single quoted string
retval(1) = octave_value (ch);        // default constructor is sq_string
           OR
retval(1) = octave_value (ch, '\'');  // explicitly create sq_string

// Create a double quoted string
retval(0) = octave_value (ch, '"');
@end group
@end example

@node Cell Arrays in Oct-Files
@subsection Cell Arrays in Oct-Files

Octave's cell type is also available from within oct-files.  A cell
array is just an array of @code{octave_value}s, and thus each element of the
cell array can be treated just like any other @code{octave_value}.  A simple
example is

@example
@verbatim
#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.numel (); i++)
          {
            retval(i) = c(i);          // using operator syntax
            //retval(i) = c.elem (i);  // using method syntax
          }
    }

  return retval;
}
@end verbatim

@end example

Note that cell arrays are used less often in standard oct-files and so
the @file{Cell.h} header file must be explicitly included.  The rest of the
example extracts the @code{octave_value}s one by one from the cell array and
returns them as individual return arguments.  For example:

@example
@group
[b1, b2, b3] = celldemo (@{1, [1, 2], "test"@})
@result{}
b1 =  1
b2 =

   1   2

b3 = test
@end group
@end example

@node Structures in Oct-Files
@subsection Structures in Oct-Files

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

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

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

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

  if (args.length () == 2)
    {
      octave_scalar_map arg0 = args(0).scalar_map_value ();
      //octave_map arg0 = args(0).map_value ();

      if (! error_state)
        {
          std::string arg1 = args(1).string_value ();

          if (! error_state)
            {
              octave_value tmp = arg0.contents (arg1);
              //octave_value tmp = arg0.contents (arg1)(0);

              if (tmp.is_defined ())
                {
                  octave_scalar_map st;

                  st.assign ("selected", tmp);

                  retval = octave_value (st);
                }
              else
                error ("structdemo: struct does not have a field named '%s'\n",
                       arg1.c_str ());
            }
          else
            error ("structdemo: ARG2 must be a character string");
        }
      else
        error ("structdemo: ARG1 must be a struct");
    }
  else
    print_usage ();

  return retval;
}
@end verbatim

@end example

An example of its use is

@example
@group
x.a = 1; x.b = "test"; x.c = [1, 2];
structdemo (x, "b")
@result{} selected = test
@end group
@end example

The example above specifically uses the @code{octave_scalar_map} class which
is for representing a single struct.  For structure arrays the
@code{octave_map} class is used instead.  The commented code shows how the
demo could be modified to handle a structure array.  In that case the
@code{contents} method returns a @code{Cell} which may have more than one
element.  Therefore, to obtain the underlying @code{octave_value} in
this single-struct example we write

@example
octave_value tmp = arg0.contents (arg1)(0);
@end example

@noindent
where the trailing (0) is the () operator on the @code{Cell} object.  If
this were a true structure array with multiple elements we could iterate
over the elements using the () operator.

Structures are a relatively complex data container and there are more
functions available in @file{oct-map.h} which make coding with them easier
than relying on just @code{contents}.

@node Sparse Matrices in Oct-Files
@subsection Sparse Matrices in Oct-Files

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

@table @code
@item SparseMatrix
A double precision sparse matrix class

@item SparseComplexMatrix
A complex sparse matrix class

@item SparseBoolMatrix
A boolean sparse matrix class
@end table

All of these classes inherit from the @code{Sparse<T>} template class,
and so all have similar capabilities and usage.  The @code{Sparse<T>}
class was based on Octave's @code{Array<T>} class, and so users familiar
with Octave's @code{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 @code{Array} classes.  However,
there are a few differences due the different nature of sparse objects,
and these will be described.  First, although it is fundamentally
possible to have N-dimensional sparse objects, the Octave sparse classes do
not allow them at this time; All instances of the sparse classes
must be 2-dimensional.  This means that @code{SparseMatrix} is actually
more similar to Octave's @code{Matrix} class than its @code{NDArray} class.

@menu
* Array and Sparse Class Differences::
* Creating Sparse Matrices in Oct-Files::
* Using Sparse Matrices in Oct-Files::
@end menu

@node Array and Sparse Class Differences
@subsubsection Array and Sparse Class Differences

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

@example
@group
SparseMatrix sm;
@dots{}
int nel = sm.nelem ();
@end group
@end example

@noindent
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 @code{numel} rather than @code{nelem}.  Note that for very
large matrices, where the product of the two dimensions is larger than
the representation of an unsigned int, then @code{numel} can overflow.
An example is @code{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 @code{numel} should therefore be avoided useless it is known
it won't overflow.

Extreme care must be take with the elem method and the @qcode{"()"} 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

@example
@group
SparseMatrix sm;
@dots{}
for (int j = 0; j < nc; j++)
  for (int i = 0; i < nr; i++)
    std::cerr << " (" << i << "," << j << "): " << sm(i,j) << std::endl;
@end group
@end example

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

@example
@group
SparseMatrix sm;
@dots{}
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;
@end group
@end example

Finally, as the sparse types aren't represented by a contiguous
block of memory, the @nospell{@code{fortran_vec}} method of the @code{Array<T>}
is not available.  It is, however, replaced by three separate methods
@code{ridx}, @code{cidx} and @code{data}, that access the raw compressed
column format that Octave sparse matrices are stored in.  These methods can be
used in a manner similar to @code{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.

@node Creating Sparse Matrices in Oct-Files
@subsubsection Creating Sparse Matrices in Oct-Files

There are several useful alternatives for creating a sparse matrix.
The first is to create three vectors representing the row index, column index,
and data values, and from these create the matrix.
The second alternative is to create a sparse matrix with the appropriate
amount of space and then fill in the values.  Both techniques have their
advantages and disadvantages.

Below is an example of creating a small sparse matrix using the first
technique

@example
@group
int nz, nr, nc;
nz = 4, nr = 3, nc = 4;

ColumnVector ridx (nz);
ColumnVector cidx (nz);
ColumnVector data (nz);

ridx(0) = 1; cidx(0) = 1; data(0) = 1;
ridx(1) = 2; cidx(1) = 2; data(1) = 2;
ridx(2) = 2; cidx(2) = 4; data(2) = 3;
ridx(3) = 3; cidx(3) = 4; data(3) = 4;
SparseMatrix sm (data, ridx, cidx, nr, nc);
@end group
@end example

@noindent
which creates the matrix given in section
@ref{Storage of Sparse Matrices}.  Note that the compressed matrix
format is not used at the time of the creation of the matrix itself,
but is used internally.

As discussed in the chapter on Sparse Matrices, the values of the sparse
matrix are stored in increasing column-major ordering.  Although the data
passed by the user need not respect this requirement, pre-sorting the
data will significantly speed up creation of the sparse matrix.

The disadvantage of this technique for creating a sparse matrix is
that there is a brief time when two copies of the data exist.  For
extremely memory constrained problems this may not be the best
technique for creating a sparse matrix.

The alternative is to first create a sparse matrix with the desired
number of non-zero elements and then later fill those elements in.
Sample code:

@example
@group
int nz, nr, nc;
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;
@end group
@end example

This creates the same matrix as previously.  Again, although not
strictly necessary, it is significantly faster if the sparse matrix is
created and the elements are added in column-major ordering.  The reason
for this is that when 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 method of creating
a sparse matrix.  First, it is possible to create a sparse matrix
with fewer elements than are actually inserted in the matrix.  Therefore,

@example
@group
int nr, nc;
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;
@end group
@end example

@noindent 
is perfectly valid.  However, it is a very bad idea because as each new
element is added to the sparse matrix the matrix needs to request more
space and reallocate 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
@var{nz} greater than 4 is also valid.  The disadvantage is that the matrix
occupies more memory than strictly needed.

It is not always possible to know the number of non-zero elements prior
to filling a matrix.  For this reason the additional unused storage of 
a sparse matrix can be removed after its creation with the
@code{maybe_compress} function.  In addition, @code{maybe_compress} can
deallocate the unused storage, but it can also remove zero elements
from the matrix.  The removal of zero elements from the matrix is
controlled by setting the argument of the @code{maybe_compress} function
to be @code{true}.  However, the cost of removing the zeros is high because it
implies re-sorting the elements.  If possible, it is better
if the user does not add the unnecessary zeros in the first place.
An example of the use of @code{maybe_compress} is

@example
@group
int nz, nr, nc;
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
@end group
@end example

The use of the @code{maybe_compress} function should be avoided if
possible as it will slow the creation of the matrix.

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

@example
octave_value arg;
@dots{}
int nz, nr, nc;
nz = 6, nr = 3, nc = 4;   // Assume we know the max # 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 # of nz.
@end example

@noindent
which is probably the most efficient means of creating a 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 @code{change_capacity} exists to reallocate the sparse memory.
The above example would then be modified as

@example
octave_value arg;
@dots{}
int nz, nr, nc;
nz = 6, nr = 3, nc = 4;   // Assume we know the max # 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 # of nz.
@end example

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 old and new copies
at the same time, additional memory is needed.  Therefore, if possible this
should be avoided.

@node Using Sparse Matrices in Oct-Files
@subsubsection Using Sparse Matrices in Oct-Files

Most of the same operators and functions on sparse matrices that are
available from the Octave command line are also available within oct-files.
The basic means of extracting a sparse matrix from an @code{octave_value}
and returning it as an @code{octave_value}, can be seen in the
following example.

@example
@group
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 ();
@dots{}
retval(2) = sbm;
retval(1) = scm;
retval(0) = sm;
@end group
@end example

The conversion to an @code{octave_value} is handled by the sparse
@code{octave_value} constructors, and so no special care is needed.

@node Accessing Global Variables in Oct-Files
@subsection Accessing Global Variables in Oct-Files

Global variables allow variables in the global scope to be
accessed.  Global variables can be accessed within oct-files by using
the support functions @code{get_global_value} and @code{set_global_value}.
@code{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 if no global variable of the desired name is found.
An example of the use of these two functions is

@example
@verbatim
#include <octave/oct.h>

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

  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;
}
@end verbatim

@end example

An example of its use is

@example
@group
global a b
b = 10;
globaldemo ("b")
@result{} 10
globaldemo ("c")
@result{} "Global variable not found"
num2str (a)
@result{} 42
@end group
@end example

@node Calling Octave Functions from Oct-Files
@subsection 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 @code{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

@enumerate 1
@item Function Handle

@item Anonymous Function Handle

@item Inline Function

@item String
@end enumerate

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

@example
@verbatim
#include <octave/oct.h>
#include <octave/parse.h>

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

  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: INPUT must be string, inline, or function handle");
    }
  return retval;
}
@end verbatim

@end example

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

@example
@group
funcdemo (@@sin,1)
@result{} 0.84147
funcdemo (@@(x) sin (x), 1)
@result{} 0.84147
funcdemo (inline ("sin (x)"), 1)
@result{} 0.84147
funcdemo ("sin",1)
@result{} 0.84147
funcdemo (@@atan2, 1, 1)
@result{} 0.78540
@end group
@end example

When the user function is passed as a string the treatment of the
function is different.  In some cases it is necessary to have the
user supplied function as an @code{octave_function} object.  In that
case the string argument can be used to create a temporary function
as demonstrated below.

@example
@group
std::octave fcn_name = unique_symbol_name ("__fcn__");
std::string fcode = "function y = ";
fcode.append (fcn_name);
fcode.append ("(x) y = ");
fcn = extract_function (args(0), "funcdemo", fcn_name,
                        fcode, "; endfunction");
@dots{}
if (fcn_name.length ())
  clear_function (fcn_name);
@end group
@end example

There are two important things to know in this case.  First, the number of
input arguments to the user function is fixed, and in the above example is
a single argument.  Second, to avoid leaving the temporary function in the
Octave symbol table it should be cleared after use.  Also, by convention
internal function names begin and end with the character sequence @samp{__}.

@node Calling External Code from Oct-Files
@subsection 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
that the declarations of the external C functions may need to be explicitly
defined as C functions to the compiler.  If the declarations of the
external C functions are in the header @file{foo.h}, then the tactic to
ensure that the C++ compiler treats these declarations as C code is

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

Calling Fortran code, however, can pose more difficulties.  This is due to
differences in the manner in which 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 @code{XSTOPX} function to
replace the Fortran @code{STOP} function.  @code{XSTOPX} uses the Octave
exception handler to treat failing cases in the Fortran code
explicitly.  Note that Octave supplies its own replacement @sc{blas}
@code{XERBLA} function, which uses @code{XSTOPX}.

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

There is no great harm in using @w{@code{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 @w{@code{F77_FCN}}
is used with code that calls @code{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

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

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

DEFUN_DLD (fortrandemo, 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.numel ();
          OCTAVE_LOCAL_BUFFER (char, ctmp, 128);

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

          retval(1) = std::string (ctmp);
          retval(0) = a;
        }
    }
  return retval;
}
@end verbatim

@end example

@noindent
and the Fortran function is

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

@end verbatim

@end example

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.  Both the Fortran and C++ files need to be compiled
in order for the example to work.

@example
@group
mkoctfile fortrandemo.cc fortransub.f
[b, s] = fortrandemo (1:3)
@result{}
  b = 1.00000   0.50000   0.33333
  s = There are   3 values in the input vector
[b, s] = fortrandemo (0:3)
error: fortrandemo: fortransub: divide by zero
@end group
@end example

@node Allocating Local Memory in Oct-Files
@subsection 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 great care must be
taken to avoid memory leaks.  The preferred manner in which to allocate
memory for use locally is to use the @w{@code{OCTAVE_LOCAL_BUFFER}} macro.
An example of its use is

@example
OCTAVE_LOCAL_BUFFER (double, tmp, len)
@end example

@noindent
that returns a pointer @code{tmp} of type @code{double *} of length
@code{len}.

In this case Octave itself will worry about reference counting and variable
scope and will properly free memory without programmer intervention.

@node Input Parameter Checking in Oct-Files
@subsection Input Parameter Checking in Oct-Files

As oct-files are compiled functions they open up the possibility of
crashing Octave through careless function calls or memory faults.
It is quite important that each and every function have a sufficient level
of parameter checking 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, in some cases 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, or if the input is complex
and a real value was expected.  Therefore, oct-files often need additional
input parameter checking.

There are several functions within Octave that can be useful for the
purposes of parameter checking.  These include the methods of the
octave_value class like @code{is_real_matrix}, @code{is_numeric_type}, etc.
Often, with a knowledge of the Octave m-file language, you can guess at what
the corresponding C++ routine will.  In addition there are some more
specialized input validation functions of which a few are demonstrated below.

@example
@verbatim
#include <octave/oct.h>

DEFUN_DLD (paramdemo, args, nargout, "Parameter Check Demo")
{
  octave_value retval;
  int nargin = args.length ();

  if (nargin != 1)
    print_usage ();
  else if (nargout != 0)
    error ("paramdemo: OUTPUT argument required");
  else
    {
      NDArray m = args(0).array_value ();
      double min_val = -10.0;
      double max_val = 10.0;
      octave_stdout << "Properties of input array:\n";
      if (m.any_element_is_negative ())
        octave_stdout << "  includes negative values\n";
      if (m.any_element_is_inf_or_nan ())
        octave_stdout << "  includes Inf or NaN values\n";
      if (m.any_element_not_one_or_zero ())
        octave_stdout << "  includes other values than 1 and 0\n";
      if (m.all_elements_are_int_or_inf_or_nan ())
        octave_stdout << "  includes only int, Inf or NaN values\n";
      if (m.all_integers (min_val, max_val))
        octave_stdout << "  includes only integers in [-10,10]\n";
    }
  return retval;
}
@end verbatim

@end example

@noindent
An example of its use is:

@example
@group
paramdemo ([1, 2, NaN, Inf])
@result{} Properties of input array:
     includes Inf or NaN values
     includes other values than 1 and 0
     includes only int, Inf or NaN values
@end group
@end example

@node Exception and Error Handling in Oct-Files
@subsection Exception and Error Handling in Oct-Files

Another important feature of Octave is its ability to react to the user
typing @key{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 @key{Control-C},
the @w{@code{OCTAVE_QUIT}} macro is supplied.  For example:

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

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

The @w{@code{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
@w{@code{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 @w{@code{OCTAVE_QUIT}} macro
in this case.  The alternative in this case is

@example
@group
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
@dots{}  some code that calls a "foreign" function @dots{}
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
@end group
@end example

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
@nospell{fortran_vec} method or the @w{@code{OCTAVE_LOCAL_BUFFER}} macro.

The Octave unwind_protect mechanism (@ref{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

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

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

DEFUN_DLD (unwinddemo, args, nargout, "Unwind Demo")
{
  octave_value retval;
  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)
        {
          // Declare unwind_protect frame which lasts as long as
          // the variable frame has scope.
          unwind_protect frame;
          frame.protect_var (current_liboctave_warning_handler);

          set_liboctave_warning_handler (my_err_handler);
          retval = octave_value (quotient (a, b));
        }
    }
  return retval;
}
@end verbatim

@end example

As can be seen in the example:

@example
@group
unwinddemo (1, 0)
@result{} Inf
1 / 0
@result{} warning: division by zero
   Inf
@end group
@end example

The warning for division by zero (and in fact all warnings) are disabled in the
@code{unwinddemo} function.

@node Documentation and Test of Oct-Files
@subsection Documentation and Test of Oct-Files

The documentation of an oct-file is the fourth string parameter of the
@w{@code{DEFUN_DLD}} macro.  This string can be formatted in the same manner
as the help strings for user functions (@pxref{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 needs to
be taken into account.  There are several possible solutions, but the most
common is illustrated in the following example,

@example
@group
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")
@{
@dots{}
@}
@end group
@end example

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

Octave also includes the ability to embed test and demonstration
code for a function within the code itself (@pxref{Test and Demo Functions}).
This can be used from within oct-files (or in fact any file) with
certain provisos.  First, the test and demo functions of Octave look
for @code{%!} as the first two characters of a line to identify test
and demonstration code.  This is a requirement for oct-files as well.
In addition, the test and demonstration code must be wrapped in a comment
block to avoid it being interpreted by the compiler.  Finally, the Octave
test and demonstration code must have access to the original 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

@example
@group
/*
%!assert (sin ([1,2]), [sin(1),sin(2)])
%!error (sin ())
%!error (sin (1,1))
*/
@end group
@end example

@c @node Application Programming Interface for Oct-Files
@c @subsection Application Programming Interface for Oct-Files
@c 
@c WRITE ME, using Coda section 1.3 as a starting point.

@node Mex-Files
@section Mex-Files
@cindex mex-files
@cindex mex

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 @sc{matlab} users.  However, as mex-files expose
@sc{matlab}'s internal API, and the internal structure of Octave is
different, a mex-file can never have the same performance in Octave as
the equivalent oct-file.  In particular, to support the manner in which
variables are passed to mex functions there are a significant number of
additional copies of memory blocks when calling or returning from a
mex-file function.  For this reason, it is recommended that any new code
be written with the oct-file interface previously discussed.

@menu
* Getting Started with Mex-Files::
* Working with Matrices and Arrays in Mex-Files::
* Character Strings in Mex-Files::
* Cell Arrays with Mex-Files::
* Structures with Mex-Files::
* Sparse Matrices with Mex-Files::
* Calling Other Functions in Mex-Files::
@c * Application Programming Interface for Mex-Files::  
@end menu

@node Getting Started with Mex-Files
@subsection Getting Started with Mex-Files

The basic command to build a mex-file is either @code{mkoctfile --mex}
or @code{mex}.  The first command can be used either from within Octave or from
the command line.  However, to avoid issues with @sc{matlab}'s own @code{mex}
command, the use of the command @code{mex} is limited to within Octave.
Compiled mex-files have the extension @file{.mex}.

@c mex scripts/miscellaneous/mex.m
@anchor{XREFmex}
@deftypefn {Command} {} mex [options] file @dots{}
Compile source code written in C, C++, or Fortran, to a MEX file.
This is equivalent to @code{mkoctfile --mex [options] file}.
@seealso{@ref{XREFmkoctfile,,mkoctfile}}
@end deftypefn


@c mexext scripts/miscellaneous/mexext.m
@anchor{XREFmexext}
@deftypefn {Function File} {} mexext ()
Return the filename extension used for MEX files.
@seealso{@ref{XREFmex,,mex}}
@end deftypefn


Consider the following short example:

@example
@group
@verbatim
#include "mex.h"

void
mexFunction (int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[])
{
  mexPrintf ("Hello, World!\n");

  mexPrintf ("I have %d inputs and %d outputs\n", nrhs, nlhs);
}
@end verbatim

@end group
@end example

The first line @code{#include "mex.h"} makes available all of the definitions
necessary for a mex-file.  One important difference between Octave and
@sc{matlab} is that the header file @qcode{"matrix.h"} is implicitly included
through the inclusion of @qcode{"mex.h"}.  This is necessary to avoid a conflict
with the Octave file @qcode{"Matrix.h"} for operating systems and compilers that
don't distinguish between filenames in upper and lower case.

The entry point into the mex-file is defined by @code{mexFunction}.  The
function takes four arguments:

@enumerate 1
@item The number of return arguments (# of left-hand side args).

@item An array of pointers to return arguments.

@item The number of input arguments (# of right-hand side args).

@item An array of pointers to input arguments.
@end enumerate

Note that the function name definition is not explicitly included in
@code{mexFunction} and so there can only be a single @code{mexFunction}
entry point per file.  Instead, the name of the function as seen in Octave is
determined by the name of the mex-file itself minus the extension.  Therefore,
if the above function is in the file @file{myhello.c}, it can be compiled with

@example
mkoctfile --mex myhello.c
@end example

@noindent
which creates a file @file{myhello.mex}.  The function can then be run from
Octave as

@example
@group
myhello (1,2,3)
@result{} Hello, World!
@result{} I have 3 inputs and 0 outputs
@end group
@end example

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 @file{myhello.m} that might
contain the text

@example
%MYHELLO Simple test of the functionality of a mex-file.
@end example

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 there cannot be multiple entry points in a single mex-file,
one can use the @code{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

@example
@group
@verbatim
#include "mex.h"

void
mexFunction (int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[])
{
  const char *nm;

  nm = mexFunctionName ();
  mexPrintf ("You called function: %s\n", nm);
  if (strcmp (nm, "myfunc") == 0)
    mexPrintf ("This is the principal function\n", nm);

  return;
}
@end verbatim

@end group
@end example

@noindent
is in file @file{myfunc.c}, and it is compiled with

@example
@group
mkoctfile --mex myfunc.c
ln -s myfunc.mex myfunc2.mex
@end group
@end example

@noindent
then as can be seen by

@example
@group
myfunc ()
@result{} You called function: myfunc
    This is the principal function
myfunc2 ()
@result{} You called function: myfunc2
@end group
@end example

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

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

@node Working with Matrices and Arrays in Mex-Files
@subsection Working with Matrices and Arrays in Mex-Files

The basic mex type of all variables is @code{mxArray}.  Any object,
such as a matrix, cell array, or structure is stored in this basic
type.  As such, @code{mxArray} serves basically the same purpose as the
octave_value class in oct-files in that it acts as a container for the
more specialized types.

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

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

@example
@group
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);
@end group
@end example

There are also the functions @code{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 @code{mxSetPr}.

Note the type @code{mwSize} used above, and also @code{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.  @code{mwSize} is used to define array dimensions
and the maximum number or elements, while @code{mwIndex} is used to define
indexing into arrays.

An example that demonstrates how to work with arbitrary real or complex
double precision arrays is given by the file @file{mypow2.c} shown below.

@example
@verbatim
#include "mex.h"

void
mexFunction (int nlhs, mxArray* plhs[],
             int nrhs, const mxArray* prhs[])
{
  mwSize n;
  mwIndex i;
  double *vri, *vro;

  if (nrhs != 1 || ! mxIsNumeric (prhs[0]))
    mexErrMsgTxt ("ARG1 must be a matrix");

  n = mxGetNumberOfElements (prhs[0]);
  plhs[0] = mxCreateNumericArray (mxGetNumberOfDimensions (prhs[0]),
                                  mxGetDimensions (prhs[0]),
                                  mxGetClassID (prhs[0]),
                                  mxIsComplex (prhs[0]));
  vri = mxGetPr (prhs[0]);
  vro = mxGetPr (plhs[0]);

  if (mxIsComplex (prhs[0]))
    {
      double *vii, *vio;
      vii = mxGetPi (prhs[0]);
      vio = mxGetPi (plhs[0]);

      for (i = 0; i < n; i++)
        {
          vro[i] = vri[i] * vri[i] - vii[i] * vii[i];
          vio[i] = 2 * vri[i] * vii[i];
        }
    }
  else
    {
      for (i = 0; i < n; i++)
        vro[i] = vri[i] * vri[i];
    }
}
@end verbatim

@end example

@noindent
with an example of its use

@example
@group
b = randn (4,1) + 1i * randn (4,1);
all (b.^2 == mypow2 (b))
@result{} 1
@end group
@end example

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

@node Character Strings in Mex-Files
@subsection 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 @file{stringdemo.cc} is given in the
file @file{mystring.c}, as shown below.

@example
@verbatim
#include <string.h>
#include "mex.h"

void
mexFunction (int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[])
{
  mwSize m, n;
  mwIndex i, j;
  mxChar *pi, *po;

  if (nrhs != 1 || ! mxIsChar (prhs[0])
      || mxGetNumberOfDimensions (prhs[0]) > 2)
    mexErrMsgTxt ("ARG1 must be a char matrix");

  m = mxGetM (prhs[0]);
  n = mxGetN (prhs[0]);
  pi = mxGetChars (prhs[0]);
  plhs[0] = mxCreateNumericMatrix (m, n, mxCHAR_CLASS, mxREAL);
  po = mxGetChars (plhs[0]);

  for (j = 0; j < n; j++)
    for (i = 0; i < m; i++)
      po[j*m + m - 1 - i] = pi[j*m + i];
}
@end verbatim

@end example

@noindent
An example of its expected output is

@example
@group
mystring (["First String"; "Second String"])
@result{} Second String
   First String
@end group
@end example

Other functions in the mex interface for handling character strings are
@code{mxCreateString}, @code{mxArrayToString}, and
@code{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.

@node Cell Arrays with Mex-Files
@subsection Cell Arrays with Mex-Files

One can perform exactly the same operations on Cell arrays in mex-files
as in oct-files.  An example that reduplicates the function of
the @file{celldemo.cc} oct-file in a mex-file is given by @file{mycell.c}
as shown below.

@example
@verbatim
#include "mex.h"

void
mexFunction (int nlhs, mxArray* plhs[],
             int nrhs, const mxArray* prhs[])
{
  mwSize n;
  mwIndex i;

  if (nrhs != 1 || ! mxIsCell (prhs[0]))
    mexErrMsgTxt ("ARG1 must be a cell");

  n = mxGetNumberOfElements (prhs[0]);
  n = (n > nlhs ? nlhs : n);

  for (i = 0; i < n; i++)
    plhs[i] = mxDuplicateArray (mxGetCell (prhs[0], i));
}
@end verbatim

@end example

@noindent
The output is identical to the oct-file version as well.

@example
@group
[b1, b2, b3] = mycell (@{1, [1, 2], "test"@})
@result{}
b1 =  1
b2 =

   1   2

b3 = test
@end group
@end example

Note in the example the use of the @code{mxDuplicateArray} function.  This
is needed as the @code{mxArray} pointer returned by @code{mxGetCell}
might be deallocated.  The inverse function to @code{mxGetCell}, used for
setting Cell values, is @code{mxSetCell} and is defined as

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

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

@example
@group
mxArray *mxCreateCellArray (int ndims, const int *dims);
mxArray *mxCreateCellMatrix (int m, int n);
@end group
@end example

@node Structures with Mex-Files
@subsection Structures with Mex-Files

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

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

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

@example
@group
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);
@end group
@end example

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 @code{index} over the elements of the
arrays of elements per @code{field}; Whereas, the oct-file structure
includes a Cell Array per field of the structure.

An example that demonstrates the use of structures in a mex-file can be
found in the file @file{mystruct.c} shown below.

@example
@verbatim
#include "mex.h"

void
mexFunction (int nlhs, mxArray* plhs[],
             int nrhs, const mxArray* prhs[])
{
  int i;
  mwIndex j;
  mxArray *v;
  const char *keys[] = { "this", "that" };

  if (nrhs != 1 || ! mxIsStruct (prhs[0]))
    mexErrMsgTxt ("expects struct");

  for (i = 0; i < mxGetNumberOfFields (prhs[0]); i++)
    for (j = 0; j < mxGetNumberOfElements (prhs[0]); j++)
      {
        mexPrintf ("field %s(%d) = ", mxGetFieldNameByNumber (prhs[0], i), j);
        v = mxGetFieldByNumber (prhs[0], j, i);
        mexCallMATLAB (0, NULL, 1, &v, "disp");
      }

  v = mxCreateStructMatrix (2, 2, 2, keys);

  mxSetFieldByNumber (v, 0, 0, mxCreateString ("this1"));
  mxSetFieldByNumber (v, 0, 1, mxCreateString ("that1"));
  mxSetFieldByNumber (v, 1, 0, mxCreateString ("this2"));
  mxSetFieldByNumber (v, 1, 1, mxCreateString ("that2"));
  mxSetFieldByNumber (v, 2, 0, mxCreateString ("this3"));
  mxSetFieldByNumber (v, 2, 1, mxCreateString ("that3"));
  mxSetFieldByNumber (v, 3, 0, mxCreateString ("this4"));
  mxSetFieldByNumber (v, 3, 1, mxCreateString ("that4"));

  if (nlhs)
    plhs[0] = v;
}
@end verbatim

@end example

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

@example
a(1).f1 = "f11"; a(1).f2 = "f12"; 
a(2).f1 = "f21"; a(2).f2 = "f22";
b = mystruct (a);
@result{}  field f1(0) = f11
    field f1(1) = f21
    field f2(0) = f12
    field f2(1) = f22
b
@result{} 2x2 struct array containing the fields:

     this
     that

b(3)
@result{} scalar structure containing the fields:

     this = this3
     that = that3
@end example

@node Sparse Matrices with Mex-Files
@subsection 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, in addition to using @code{mxGetM}, @code{mxGetN},
@code{mxSetM}, @code{mxSetN}, @code{mxGetPr}, @code{mxGetPi},
@code{mxSetPr}, and @code{mxSetPi}, also supplies the following functions.

@example
@group
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);
@end group
@end example

@noindent
@code{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.  @code{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
@code{mxGetJc} define the number of non-zero elements in each column of
the sparse matrix.  Therefore,

@example
@group
mwSize nz, n;
mwIndex *Jc;
mxArray *m;
@dots{}
n = mxGetN (m);
Jc = mxGetJc (m);
nz = Jc[n];
@end group
@end example

@noindent
returns the actual number of non-zero elements stored in the matrix in
@code{nz}.  As the arrays returned by @code{mxGetPr} and @code{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
@code{mxGetIr}.  A complete example of the use of sparse matrices in
mex-files is given by the file @file{mysparse.c} shown below.

@example
@verbatim
#include "mex.h"

void
mexFunction (int nlhs, mxArray *plhs[],
             int nrhs, const mxArray *prhs[])
{
  mwSize m, n, nz;
  mxArray *v;
  mwIndex i;
  double *pr, *pi;
  double *pr2, *pi2;
  mwIndex *ir, *jc;
  mwIndex *ir2, *jc2;

  if (nrhs != 1 || ! mxIsSparse (prhs[0]))
    mexErrMsgTxt ("ARG1 must be a sparse matrix");

  m = mxGetM (prhs[0]);
  n = mxGetN (prhs[0]);
  nz = mxGetNzmax (prhs[0]);

  if (mxIsComplex (prhs[0]))
    {
      mexPrintf ("Matrix is %d-by-%d complex sparse matrix", m, n);
      mexPrintf (" with %d elements\n", nz);

      pr = mxGetPr (prhs[0]);
      pi = mxGetPi (prhs[0]);
      ir = mxGetIr (prhs[0]);
      jc = mxGetJc (prhs[0]);

      i = n;
      while (jc[i] == jc[i-1] && i != 0) i--;

      mexPrintf ("last non-zero element (%d, %d) = (%g, %g)\n",
                 ir[nz-1]+ 1, i, pr[nz-1], pi[nz-1]);

      v = mxCreateSparse (m, n, nz, mxCOMPLEX);
      pr2 = mxGetPr (v);
      pi2 = mxGetPi (v);
      ir2 = mxGetIr (v);
      jc2 = mxGetJc (v);

      for (i = 0; i < nz; i++)
        {
          pr2[i] = 2 * pr[i];
          pi2[i] = 2 * pi[i];
          ir2[i] = ir[i];
        }
      for (i = 0; i < n + 1; i++)
        jc2[i] = jc[i];

      if (nlhs > 0)
        plhs[0] = v;
    }
  else if (mxIsLogical (prhs[0]))
    {
      mxLogical *pbr, *pbr2;
      mexPrintf ("Matrix is %d-by-%d logical sparse matrix", m, n);
      mexPrintf (" with %d elements\n", nz);

      pbr = mxGetLogicals (prhs[0]);
      ir = mxGetIr (prhs[0]);
      jc = mxGetJc (prhs[0]);

      i = n;
      while (jc[i] == jc[i-1] && i != 0) i--;
      mexPrintf ("last non-zero element (%d, %d) = %d\n",
                 ir[nz-1]+ 1, i, pbr[nz-1]);

      v = mxCreateSparseLogicalMatrix (m, n, nz);
      pbr2 = mxGetLogicals (v);
      ir2 = mxGetIr (v);
      jc2 = mxGetJc (v);

      for (i = 0; i < nz; i++)
        {
          pbr2[i] = pbr[i];
          ir2[i] = ir[i];
        }
      for (i = 0; i < n + 1; i++)
        jc2[i] = jc[i];

      if (nlhs > 0)
        plhs[0] = v;
    }
  else
    {
      mexPrintf ("Matrix is %d-by-%d real sparse matrix", m, n);
      mexPrintf (" with %d elements\n", nz);

      pr = mxGetPr (prhs[0]);
      ir = mxGetIr (prhs[0]);
      jc = mxGetJc (prhs[0]);

      i = n;
      while (jc[i] == jc[i-1] && i != 0) i--;
      mexPrintf ("last non-zero element (%d, %d) = %g\n",
                 ir[nz-1]+ 1, i, pr[nz-1]);

      v = mxCreateSparse (m, n, nz, mxREAL);
      pr2 = mxGetPr (v);
      ir2 = mxGetIr (v);
      jc2 = mxGetJc (v);

      for (i = 0; i < nz; i++)
        {
          pr2[i] = 2 * pr[i];
          ir2[i] = ir[i];
        }
      for (i = 0; i < n + 1; i++)
        jc2[i] = jc[i];

      if (nlhs > 0)
        plhs[0] = v;
    }
}
@end verbatim

@end example

A sample usage of @code{mysparse} is

@example
@group
sm = sparse ([1, 0; 0, pi]);
mysparse (sm)
@result{}
Matrix is 2-by-2 real sparse matrix with 2 elements
last non-zero element (2, 2) = 3.14159
@end group
@end example

@node Calling Other Functions in Mex-Files
@subsection Calling Other Functions in Mex-Files

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

@example
@verbatim
#include "mex.h"

void
mexFunction (int nlhs, mxArray* plhs[],
             int nrhs, const mxArray* prhs[])
{
  char *str;

  mexPrintf ("Starting file myfeval.mex\n");

  mexPrintf ("I have %d inputs and %d outputs\n", nrhs, nlhs);

  if (nrhs < 1 || ! mxIsString (prhs[0]))
    mexErrMsgTxt ("ARG1 must be a function name");

  str = mxArrayToString (prhs[0]);

  mexPrintf ("I'm going to call the function %s\n", str);

  if (nlhs == 0)
    nlhs = 1;  // Octave's automatic 'ans' variable

  /* Cast prhs just to get rid of 'const' qualifier and stop compile warning */
  mexCallMATLAB (nlhs, plhs, nrhs-1, (mxArray**)prhs+1, str);

  mxFree (str);
}
@end verbatim

@end example

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

@example
@group
a = myfeval ("sin", 1)
@result{} Starting file myfeval.mex
   I have 2 inputs and 1 outputs
   I'm going to call the interpreter function sin
   a =  0.84147
@end group
@end example

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

@c @node Application Programming Interface for Mex-Files
@c @subsection Application Programming Interface for Mex-Files
@c 
@c WRITE ME, refer to mex.h and mexproto.h

@node Standalone Programs
@section 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 of the Octave algorithms.  The
following C++ program, uses class Matrix from @file{liboctave.a} or
@file{liboctave.so}.

@example
@verbatim
#include <iostream>
#include <octave/oct.h>

int
main (void)
{
  std::cout << "Hello Octave world!\n";

  int n = 2;
  Matrix a_matrix = Matrix (n, n);

  for (octave_idx_type i = 0; i < n; i++)
    for (octave_idx_type j = 0; j < n; j++)
      a_matrix(i,j) = (i + 1) * 10 + (j + 1);

  std::cout << a_matrix;

  return 0;
}
@end verbatim

@end example

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

@example
@group
$ mkoctfile --link-stand-alone standalone.cc -o standalone
$ ./standalone
Hello Octave world!
  11 12
  21 22
$
@end group
@end example

Note that the application @code{standalone} 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 built-in 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

@example
@verbatim
#include <iostream>
#include <octave/oct.h>
#include <octave/octave.h>
#include <octave/parse.h>
#include <octave/toplev.h>

int
main (void)
{
  string_vector argv (2);
  argv(0) = "embedded";
  argv(1) = "-q";

  octave_main (2, argv.c_str_vec (), 1);

  octave_idx_type n = 2;
  octave_value_list in;

  for (octave_idx_type i = 0; i < n; i++)
    in(i) = octave_value (5 * (i + 2));

  octave_value_list out = feval ("gcd", in, 1);

  if (! error_state && out.length () > 0)
    std::cout << "GCD of ["
              << in(0).int_value ()
              << ", "
              << in(1).int_value ()
              << "] is " << out(0).int_value ()
              << std::endl;
  else
    std::cout << "invalid\n";

  clean_up_and_exit (0);
}
@end verbatim

@end example

@noindent
which, as before, is compiled and run as a standalone application with

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

It is worth noting that, if only built-in functions are to be called from 
a C++ standalone program, then it does not need to initialize the 
interpreter to do so.  The general rule is that, for a built-in 
function named @code{function_name} in the interpreter, there will be
a C++ function named @code{Ffunction_name} (note the prepended capital 
@code{F}) accessible in the C++ API@.  The declarations for all built-in 
functions are collected in the header file @code{builtin-defun-decls.h}.
This feature should be used with care as the list of built-in functions can
change.  No guarantees can be made that a function that is currently built in
won't be implemented as a .m file or as a dynamically linked function in the
future.  An example of how to call built-in functions from C++ can be seen in the
code

@example
@verbatim
#include <iostream>
#include <octave/oct.h>
#include <octave/builtin-defun-decls.h>

int
main (void)
{
  int n = 2;
  Matrix a_matrix = Matrix (n, n);

  for (octave_idx_type i = 0; i < n; i++)
    for (octave_idx_type j = 0; j < n; j++)
      a_matrix(i,j) = (i + 1) * 10 + (j + 1);

  std::cout << "This is a matrix:" << std::endl 
            << a_matrix            << std::endl;

  octave_value_list in;
  in(0) = a_matrix;

  octave_value_list out = Fnorm (in, 1);
  double norm_of_the_matrix = out(0).double_value ();

  std::cout << "This is the norm of the matrix:" << std::endl 
            << norm_of_the_matrix                << std::endl;
  
  return 0;
}
@end verbatim

@end example

@noindent
which, again, is compiled and run as a standalone application with

@example
@group
$ mkoctfile --link-stand-alone standalonebuiltin.cc -o standalonebuiltin
$ ./standalonebuiltin 
This is a matrix:
 11 12
 21 22

This is the norm of the matrix:
34.4952
$
@end group
@end example