JBLopen.com

I was surprised by the amount of people who asked for the code used to render the fractal in the screen shot of a previous blog entry. Since I used an external function written in C to generate the fractal the code needs a bit of explanation in order to get it working. As such, I've decided to post this entry for everyone who asked for the source (All right, all right all three of them). At the same time this is an excellent demonstration of how easy and powerful Scilab external routines can be.

This discussion is mostly aimed at Unix/Linux users, while the process is similar for Windows it is much easier to work with *NIX shared objects than Windows DLL.

The full Mandelbrot fractal. (Click for full resolution.)

Scilab extensions

Currently there is a number of ways to extend Scilab with custom functions. In this entry I will only discuss of the simplest method but for sake of completeness here is a brief list of the extension methods available.

External function call (Long format)

This is the simplest way of creating an external routine for use within Scilab. In this scenario a dynamic shared object is loaded and functions are registered in Scilab by specifying the functions names and arguments. In this blog entry we are going to explore this method since it's the quickest albeit being somewhat more costly and less powerful than other methods. Functions loaded this way can also be used as external function used in some ODE and optimization functions.

External function call using Scilab gateways / Toolbox

The second method is very similar to the first one except that a specially written gateway function takes care of passing the arguments greatly simplifying the calling of external functions. This method is also faster since it induce less overhead. Writing a whole toolbox is not much more complicated than this when using a loader and the full Scilab API.

Scilab module

A Scilab module is similar to a toolbox except it is built alongside Scilab becoming part of the main Scilab code tree. There is little advantages of this method over a toolbox unless one needs a deeply modified version of Scilab for specific uses. Modifying or updating a module also requires a complete rebuild of Scilab.

MEX

Scilab also supports the MEX API commonly used with Matlab external modules. Using this, someone can write eternal routines compatible with both Scilab and Matlab. I would not however recommend this method for libraries intended only for Scilab.

Writing a fractal function for Scilab

Writing an external function is rather simple and can be written like a plain C function which is the beauty of this method. No need to learn the Scilab API or master the complex toolbox loading mechanism. Not all Scilab data types are available though, so we must constraint ourselves to double , int and strings. Complex data must be passed as two arrays of double precision values. All arguments are passed by reference and matrices are passed as standard one-dimensional C arrays. Just a note of warning, arrays in Scilab follow the Fortran convention and are stored in memory column by column known as the column-major matrix ordering. So when interfacing Scilab with other C matrix algebra libraries care must be taken to perform the correct conversion.

Without further delay, the fractal function looks like this :

void GenerateFractal( double xvect[], double yvect[], int *xsize, int *ysize, int *nmax, double Output[] )
{
    int i, j, k;
    double x, xtemp, y, x0, y0;

    for( i = 0; i < *xsize; i++ )
    {
        for( j = 0; j < *ysize; j++ )
        {
            x = xvect[i];
            y = yvect[j];
            x0 = x;
            y0 = y;
            k = 0;

            while( (x*x + y*y) < 4 && k < *nmax )
            {
                xtemp = x*x - y*y + x0;
                y = 2*x*y + y0;
                x = xtemp;
                k++;
            }

            Output[ j + i * (*xsize) ] = k + 1 - (log(log((sqrt(x*x + y*y)))))/(log(2));
        }
    }

}

Most of the code is self explanatory, I elected to pass coordinate vectors for each points to be evaluated, xvect and yvect in the code alongside with their respective size xsize and ysize. The maximum number of iterations is passed as nmax followed by the output vector. The fractal code itself is the simplest implementation of evaluating the Mandlebrot set to generate a fractal. The only slightly curious part is the last line. Since we are only evaluating the Mandelbrot set for graphing we can use many different functions to give a better output, the function used give a better colour gradient in the sections that escape quickly. I do not remember where I found it however.

Before calling it, the C file must be compiled into a shared library, one can compile a shared library simply by invoking this command :

gcc -fPIC -O2 -shared -o libfract.so fract.c

When building for a 64 bit distribution of Scilab we must be careful to add the -m64 option.

Calling a custom C function from within Scilab

Before actually calling an external function it must be linked in order to load it into the Scilab address space. The link command loads a single function from a specific shared object like this :

so1 = link('fractal.so', 'GenerateFractal', 'c' );

The third argument specifies the language of the linked object, either 'c' or 'fortran'. Afterwards the function is available to be called using the call command. Without using a gateway function the call command is a little bit complicated since it must specify exactly the arguments of the called function. It is also relatively unsafe since most errors will generate a segfault as Scilab has no easy way of verifying the validity of the call sequence. Input arguments are listed first by specifying for each the passed Scilab variable, the position of the argument and it's target type. Output arguments are listed after the “out” keyword and must be specified by giving the size, position and type of each. The correct call line for the GenerateFractal function is as follow :

res = call( 'GenerateFractal', xvect, 1, "d", yvect, 2, "d", xsize, 3, "i", ysize, 4, "i", nmax, 5, "i", "out", [xsize,ysize],  6, "d" );

The abridged script goes like this, a complete version is available for download at the end of this entry :

so1 = link('libfract.so', 'GenerateFractal', 'c' );

nmax = 1000;
xmin = 0.2675;
xmax = 0.2685;
ymin = 0.591;
ymax = 0.592;
xsize = 1000;
ysize = 1000;

map = hotcolormap(512);
cmap = [map(:,3), map(:,2), map(:,1)];

f=get("current_figure");
f.color_map = cmap;
f.auto_resize = "off";
f.axes_size = [1000,1000];
f.figure_size = [800,800];

xvect = linspace( xmin, xmax, xsize );
yvect = linspace( ymin, ymax, ysize );

res = call( 'GenerateFractal', xvect, 1, "d", yvect, 2, "d", xsize, 3, "i", ysize, 4, "i", nmax, 5, "i", "out", [xsize,ysize],  6, "d" );

fact = (max(res)) / 512
res = res./fact;
clf();
a=get("current_axes");
a.margins = [0,0,0,0];
Matplot(res, "080");

The source code in the archive contains more comments and is released in the public domain.

Higly zoomed region of the Mandelbrot fractal. (Click for full resolution.)

Going further

Now it is certainly out of question to show C code for generating a fractal without also giving a parallel version of said code. With OpenMP it's as easy as one additional line of code :

#pragma omp parallel for shared(Output,xvect,yvect,xsize,ysize,nmax) private(i,j,x,y,x0,y0,k,xtemp) schedule(dynamic,100)
    for( i = 0; i < *xsize; i++ )
    {

This is another advantages of external functions is that Scilab for now is mostly single threaded, but it's possible to write parallel code in C or Fortran without too much pain. To compile this version one needs to add the -fopenmp option when invoking GCC or else the omp pragma will be ignored.

Conclusion

In this entry I shown the simplest way to create custom functions in C for use within Scilab. This method is best suited when relatively few functions are needed and they are not called very often. For example in this case the function is called once to process a bunch of data which would have been excessively slow to do purely in Scilab. In general, algorithms requiring to iterate over a lot of data are well suited for external functions. It's also a wonderful framework for algorithms development where Scilab becomes a test harness for the code under study.

Source code : scilabfract.tar.gz