lacaml / README.md
LACAML  Linear Algebra for OCaml
What is LACAML?
This OCamllibrary interfaces two widely used mathematical FORTRANlibraries:
This allows developers to write highperformance numerical code for applications that require linear algebra.
Features

The BLAS and LAPACKlibraries have evolved over about two decades of time and are therefore extremely mature both in terms of stability and performance.

LACAML interfaces most of the functions in BLAS and LAPACK (many hundreds!). It supports among other things linear equations, least squares problems, eigenvalue problems, singular value decomposition (SVD), Cholesky and QRfactorization, etc.

Many convenience functions for creating and manipulating matrices.

Powerful printing functions for large vectors and matrices and supplemental information (e.g. row and column headers). Users can specify easily how much context to print. For example, it is usually sufficient to print small blocks of the four corners of a large result matrix to manually verify the correctness of an algorithm. LACAML uses this approach to limit the output to humanmanageable size.

Integration into the OCamltoplevel allows for easy experimentation for students and researchers as well as demonstration for lecturers. Values of vector and matrix type will be printed automatically without cluttering the screen.

The OCamlinterface was designed in a way to combine both the possibility of gaining optimum efficiency (e.g. by allowing the creation of work arrays outside of loops) with simplicity (thanks to labels and default arguments).

The code is precisionindependent and supports both real and complex transforms in a consistent way. There are four modules that implement the same interface modulo the precision type and specialized real/complex functions. If you refer to elements in this interface only, your code becomes precision and (if meaningful) real/complex independent, too: you can choose at anytime whether you want to use singleprecision or doubleprecision simply by referring to the required module.

You can fully exploit the library within multithreaded programs. Many numerical routines are likely to run for a long time, but they will never block other threads. This also means that you can execute several routines at the same time on several processors if you use POSIXthreads in OCaml.

To make things easy for developers used to the "real" implementation in FORTRAN but also for beginners who need detailed documentation, both function and argument names have been kept compatible to the ones used in the BLAS and LAPACKdocumentation. Only exception: you need not prefix functions with
s
,d
,c
orz
to indicate the precision and type of numbers, because the OCaml module system provides us with a more convenient means of choosing them. 
(Almost) all errors are handled within OCaml. Typical mistakes like passing nonconforming matrices, parameters that are out of range, etc., will be caught before calling Fortran code and will raise exceptions. These exceptions will explain the error in detail, for example the received illegal parameter and the range of expected legal values.
The only exception to the above is illegal contents of vectors and matrices. This can happen, for example, when freshly allocated matrices are used without initialization. Some LAPACKalgorithms may not be able to deal with floats that correspond to NaNs, infinities, or are subnormal. Checking matrices on every call would seem excessive. Some functions also expect matrices with certain properties, e.g. positivedefiniteness, which would be way too costly to verify beforehand. It is the task of the user to make sure that data contained in matrices is suitable for the application of the intended functions.
Using LACAML
You can make use of this library by referring to the corresponding module for the required precision and number type. E.g.:
open Lacaml.S (* Singleprecision real numbers *) open Lacaml.D (* Doubleprecision real numbers *) open Lacaml.C (* Singleprecision complex numbers *) open Lacaml.Z (* Doubleprecision complex numbers *)
These modules become available if you link the lacaml
library with your
application. The widely used OCamltool findlib
will take care of linking
lacaml
correctly. If you do not use this tool, you will also have to link
in the bigarray
library provided by the OCamldistribution.
The Lacaml.?
modules implement the BLAS/LAPACKinterface. Their
corresponding submodules Vec
and Mat
provide for vector and matrix
operations that relate to the given precision and number type.
Most functions were implemented using optional arguments (= default arguments). If you do not provide them at the callsite, sane defaults will be used instead. Here is an example of a function call:
let rank = gelss in_mat out_mat in (* ... *)
This example computes the solution to a general least squares problem (=
linear regression) using the SVDalgorithm with in_mat
as the matrix
containing the predictor variables and out_mat
as the matrix containing
(possibly many) response variables (this function can handle several response
variables at once). The result is the rank of the matrix. The matrices
provided in the arguments will be overwritten with further results (here:
the singular vectors and the solution matrix).
If the above happened in a loop, this would be slightly inefficient, because
a workarray would have to be allocated (and later deallocated) at each call.
You can hoist the creation of this work array out of the loop, e.g. (m
,
n
, nrhs
are problem dependent parameters):
let work = gelss_min_work ~m ~n ~nrhs in for i = 1 to 1000 do (* ... *) let rank = gelss in_mat ~work out_mat in (* ... *) done
All matrices can be accessed in a restricted way, i.e. you can specify
submatrices for all matrix parameters. For example, if some matrix is called
a
in the interface documentation, you can specify the left upper corner of
the wanted submatrix for the operation by setting ar
for the row and ac
for the column (1 by default). A vector y
would have an extra optional
parameter ofsy
(also 1 by default). Parameters like m
or n
typically
specify the numbers of rows or columns.
Printing vectors and matrices
Here is a toplevel example of printing a large random matrix:
# #require "lacaml";; # open Lacaml.D;; # let mat = Mat.random 100 200;; val mat : Lacaml.D.mat = C1 C2 C3 C198 C199 C200 R1 0.314362 0.530711 0.309887 ... 0.519965 0.230156 0.0479154 R2 0.835658 0.581404 0.161607 ... 0.749358 0.630019 0.858998 R3 0.403421 0.458116 0.497516 ... 0.210811 0.422094 0.589661 ... ... ... ... ... ... ... R98 0.352474 0.878897 0.357842 ... 0.150786 0.74011 0.353253 R99 0.104805 0.984924 0.319127 ... 0.143679 0.858269 0.859059 R100 0.419968 0.333358 0.237761 ... 0.483535 0.0224016 0.513944
Only the corner sections of the matrix, which would otherwise be too large
to display readably, are being printed, and ellipses (...
) are used in
place of the removed parts of the matrix.
If the user required even less context, the Lacaml.Io.Toplevel.lsc
function,
which is also available in each precision module for convenience (here:
Lacaml.D
), could be used to indicate how much. In the following example
only twobytwo blocks are requested in each corner of the matrix:
# lsc 2;;  : unit = () # mat;;  : Lacaml.D.mat = C1 C2 C199 C200 R1 0.314362 0.530711 ... 0.230156 0.0479154 R2 0.835658 0.581404 ... 0.630019 0.858998 ... ... ... ... ... R99 0.104805 0.984924 ... 0.858269 0.859059 R100 0.419968 0.333358 ... 0.0224016 0.513944
Applications can use the standard Format
module in the OCamldistribution
together with LACAML printing functions to output vectors and matrices.
Here is an example using labels and showing the high customizability of the
printing functions:
open Lacaml.D open Lacaml.Io let () = let rows, cols = 200, 100 in let a = Mat.random rows cols in Format.printf "@[<2>This is an indented random matrix:@\n@\n%a@]@." (Lacaml.Io.pp_lfmat ~row_labels: (Array.init rows (fun i > Printf.sprintf "Row %d" (i + 1))) ~col_labels: (Array.init cols (fun i > Printf.sprintf "Col %d" (i + 1))) ~vertical_context:(Some (Context.create 2)) ~horizontal_context:(Some (Context.create 3)) ~ellipsis:"*" ~print_right:false ~print_foot:false ()) a
The above code might print:
This is an indented random matrix: Col 1 Col 2 Col 3 Col 98 Col 99 Col 100 Row 1 0.852078 0.316723 0.195646 * 0.513697 0.656419 0.545189 Row 2 0.606197 0.411059 0.158064 * 0.368989 0.2174 0.9001 * * * * * * * Row 199 0.684374 0.939027 0.000699582 * 0.117598 0.285587 0.654935 Row 200 0.929341 0.823264 0.895798 * 0.198334 0.725029 0.621723
Many other options, e.g. for different padding, printing numbers in other formats or with different precision, etc., are available for output customization.
Other sources of usage information
API documentation
Please refer to the APIdocumentation that will be generated as HTMLfiles
if requested (make doc
). It will give more details on the numerous
functions and supported arguments in the library. It can also be found
online.
BLAS/LAPACK man pages
BLAS and LAPACK binary packages for Unix operating systems usually come with appropriate manpages. E.g. to quickly look up how to factorize a positivedefinite, complex, single precision matrix, you might enter:
man cpotrf
The corresponding function in Lacaml would be Lacaml.C.potrf
. The naming
conventions and additional documentation for BLAS and LAPACK can be found
at their respective websites.
Examples
The examples
directory contains several demonstrations of how to use this
library for various linear algebra problems.
Improving Performance
It is highly recommended that users install a variant of BLAS (or even LAPACK)
that has been optimized for their system. Processor vendors (e.g. Intel)
usually sell the most optimized implementations for their CPUarchitectures.
Some computer and OSvendors like Apple distribute their own implementations
with their products (e.g. vecLib
).
There is also ATLAS, a very efficient and compatible substitute for BLAS. It specializes code for the architecture it is compiled on. Binary packages (e.g. RPMs) for Linux should be available from your distribution vendor's site (you must recompile the package to make sure it is suited to your distribution, see the package documentation for more details.).
Another alternative for BLAS is OpenBLAS.
Contact Information and Contributing
In the case of bugs, feature requests, contributions and similar, please contact the maintainers:
 Markus Mottl markus.mottl@gmail.com
 Christophe Troestler Christophe.Troestler@umons.ac.be
Uptodate information concerning this tool should be available at: https://bitbucket.org/mmottl/lacaml
Enjoy!
Markus Mottl in Rutherford, NJ on July 9, 2012