C, MEX, BLAS, LAPACK, Fortran and Java

Calling C from MATLAB using the MEX interface

C is presently the lingua franca of computing. If you are a C programmer, you should learn how to call the BLAS from C and how to call C routines from MATLAB using MEX. If you are not a C programmer, this is a good time to learn a little C.

MEX stands for MATLAB Executable. MEX-files are dynamically linked subroutines produced from C or Fortran source code that, when compiled, can be run from within MATLAB in the same way as MATLAB M-files or built-in functions. You can find complete documentation on how to use MEX in the MATLAB help facility, but the best way to learn is to look at the examples here and change them carefully. I strongly recommend that you use a Sun Solaris account when you use MEX. Take a look at

clgs.c (C program for classical Gram-Schmidt (unstable), together with MEX interface to BLAS in the Sun Performance library)

Download this file and save it in whatever working directory you are using. We can compile this program from inside MATLAB, but we need to use the right MEX interface options, set in the file

mexopts.sh (mex options file)

Download this file and save it as a file .matlab/R14/mexopts.sh in your home directory. (If you use your own computer, you are on your own. You would have to figure out how to set up the mexopts.sh options to call your C compiler (assuming you have one), and you would have to install BLAS if they are not installed.) Then, from the MATLAB prompt, type mex clgs.c and then [Q,R] = clgs(A) as usual from the MATLAB prompt. What if you also have a clgs.m file? Which one will it use? You can find out by typing which clgs. Another way to make sure that MATLAB is really calling the C code is to call it with only one output argument: Q = clgs(A): this should generate an error message saying clgs.c requires two output arguments.

When you write your own C program to be called from MATLAB, use clgs.c as a template. If you are new at C, don't be alarmed when you find you have tons of error messages. The first thing to check is that your comments are written correctly. There are two styles of comments: /*.....*/ often running over several lines (those extra *'s on the lines in between are purely stylistic), and // for one-liners, usually on the same line as actual code, but following it. Another common error is forgetting semicolons which are required after every statement, unlike in Matlab. The tough bugs, however, usually concern allocation of space and parameter passing. Be very careful to follow the template in the example file to avoid trouble.

You can also look (NOT READY YET...) at

lapacktest.c (a test program for LAPACK: it calls the QR factorization routine but doesn't do anything useful and returns nothing)

This includes a sample mex interface for calling LAPACK routines. Use it as a template for calling LAPACK. You have to change it appropriately to do anything useful.

The BLAS

The BLAS are Basic Linear Algebra Subroutines. They come in 3 kinds:

Level 1 for basic vector operations (e.g. DAXPY: adding a scalar times a vector to another vector, and DDOT: the vector inner product)
Level 2 for basic matrix-vector operations (e.g. matrix-vector multiply)
Level 3 for basic matrix-matrix operations such as matrix-matrix multiply

They are callable from C and Fortran. Because vector operations are built into Fortran 90, it is not necessary to explicitly call the Level 1 BLAS when using Fortran 90. Highly efficient machine-specific implementations of the BLAS are available for most high-performance computers, and following the example shown in clgs.c allows you to access the precompiled BLAS in the Sun performance library. The BLAS enable LAPACK (and other) routines to achieve high performance with transportable software. You can easily find information about the BLAS on the web. For example, the Sun documentation on DAXPY is here. Notice the difference between the C and Fortran interfaces! This is discussed more below and in the comments in clgs.c. I have put a local copy of the Sun performance library user's guide here .

LAPACK

LAPACK is a free package of linear algebra subroutines, based on the BLAS. These were written in Fortran, and the compiled routines are callable from Fortran and C. You can also obtain "f2c" versions, but these are not recommended. When doing numerical computing in C or Fortran, it is always an excellent idea to use LAPACK as much as possible. As of Version 6, MATLAB actually calls LAPACK routines too. For details on individual routines, go to the LAPACK site. You can also find this information and much more information on numerical software at NETLIB, as well as the LAPACK user's guide.

Tips on the Differences Between Fortran and C Programming

When invoking a compiled Fortran library subroutine, e.g. a BLAS or an LAPACK routine, from C, it is necessary to put the name in lower case followed by an underscore (_), and declare it using ``extern", e.g. ``extern void daxpy_(int *, ...)". The reason is that the Fortran compiler adds the underscore to the routine name when it generates the compiled code. The underscore is not used when calling the routine from a Fortran program. However, it's better to call the C versions of the BLAS and LAPACK routines in the Sun Performance library directly; these are wrappers for the Fortran routines and therefore take care of the interface for you. This way, you don't need to write either the underscore or the extern declarations, and you use standard C parameter conventions, passing scalars that don't return values by value, not by reference.
Fortran uses only call by reference. To call a Fortran routine from C, it's necessary to pass pointers to the relevant variables, not the variables themselves; hence the &'s in the calling statements and the *'s in the formal parameter lists. Again, this is avoided when you call the C wrapper routines in the Sun performance library.
Two-dimensional arrays in Fortran are stored by columns, i.e. as if the matrix is a one-dimensional array [A_11 A_21 ... A_m1 A_12 ... A_m2 ...]. Thus, to pass the j-th column of a matrix A in Fortran 77, you pass A(1,j), which is the address in memory where A(1,j) is stored. In Fortran 90, you pass either A(1,j) or A(:,j). MATLAB also stores matrices by columns. In C, the easiest way to implement two-dimensional arrays is as a pointer to an array of pointers. Thus A[j] points to the beginning of the j-th column, and A[j][i] is the i-th row element of the j-th column of A. Thus the floating point data in each column must be stored consecutively in memory.
Fortran arrays start indexing at 1 by default (this is changeable in Fortran 90); C arrays start at 0.
Dynamic allocation, commonly used in C, is part of the Fortran 90 standard but not the Fortran 77 standard, so Fortran LAPACK routines typically require the passing of work arrays. Again this is avoided when you use the C interface in the Sun performance library.
The first 6 columns in Fortran source code normally must be left blank. Columns after column 72 are normally ignored.
Fortran 77 will not accept semicolons; C insists on them.
Both C and Fortran can be "mexed" so routines are callable from MATLAB.
Fortran 77 uses single quotes and C uses double quotes for strings.
In C, a function may or may not return a value. In Fortran, it must; the alternative is subroutine, which does not.

Fortran Programming

Fortran is traditionally the language of numerical computing and is no longer the dinosaur it used to be (Fortran 77). Fortan is probably still more widely used than C in scientific applications, though this is shifting. Fortran 90 looks a lot like MATLAB: vector operations are built in, and so the Level 1 BLAS are not needed. If you insist on using Fortran instead of C, you will have to figure out the mex interface yourself.

Java

Java emerged in the 1990's as an important new language. At present, Java is not a language of choice for scientific computing because it is relatively slow compared to C and Fortran.....though not as slow as it was when it was first announced. See the Java Numerics web page to learn about efforts to promote the use of Java for numerical computing. Java methods can be called directly from MATLAB, without using the MEX interface. Type help java at the MATLAB prompt for more information. You can call the Fortran BLAS using the Java Native Interface, but this should only be attempted by experienced Java programmers. It's better to get Java versions from the web.

Back to course home page