KKTDirect 0.5
Robert Bridson (rbridson@cs.ubc.ca)

This project is released into the public domain.

This is an alpha release of KKTDirect, a direct solver package for
sparse symmetric indefinite matrices of so-called KKT structure. I
make no guarantees of correctness or the ability to fail gracefully.
The interface is not well documented; see kktsolve.cpp, fulltest.cpp or
kktsolve_mex.cpp (which does implement a simple MATLAB interface to
the solver) for an outline of how one might call it. In particular,
kktsolve.cpp provides a single function for factoring and solving a linear
system which hides the complexity of the intermediate calls to the
supernodal version of the solver. Use at your own risk.

The code should be able to solve problems with matrices of the form:
   [ A   B ]
   [ B^T -C ]
where A is symmetric positive definite, B has full column rank, and C
is symmetric positive semi-definite (typically C=0). If the LDL^T
factorization of the unpermuted system works (apart from memory problems
from fill-in) then this code should work; it has succeeded on more
challenging problems but I haven't derived guarantees. (Note that
in particular the factorization may be unstable if A is much smaller
than B, even though the matrix as a whole may be well-conditioned.)

A report describing how the algorithms work is included as a
PDF preprint (kktdirect.pdf) for a document-in-progress
describing the underlying ideas and some numerical experiments.
Note that the numerical experiments in the document were run with
version 0.3 of the code.

To make the software, edit the Makefile for your platform
(determining where BLAS/LAPACK libraries are etc. or setting USE_NAIVE_* if
unavailable) then run make.
The following should be produced:
  libkktdirect.a: a library containing the KKTDirect functions
  fulltest: a commandline executable which takes an ASCII description of a
            matrix from standard input and runs several of the KKTDirect
            routines to test the library is working.
  stokes: a commandline executable which takes a grid size as argument and
          outputs an ASCII description of the matrix arising from a standard
          staggered grid discretization of the 2D Stokes problem (suitable
          for testing with fulltest).
  stokes3d: the same thing, but for the 3D problem. Note that unlike
            stokes, this produces a singular matrix (with null-space
            corresponding to constant pressures) and thus factorization
            may be reported to fail. However, only the very last pivot
            will be zero, so a sensible solution will be returned with
            the residual close to optimal.
  kktsolve_mex: a MATLAB MEX file for solving saddlepoint problems (see
                kktsolve.m for more details)

-----------------------------------------------------------------------------
New since version 0.4:
- The BLAS calls are cleaned up a bit.
- There is now an option to make the library without an external BLAS and
  LAPACK (#define USE_NAIVE_BLAS and #define USE_NAIVE_LAPACK). It is
  suggested to use this only as a last resort if you cannot get a working
  BLAS and LAPACK linked: the performance and in some cases reliability
  are not going to be as good.
- fixed an uninitialized memory bug in symbolic factorization

-----------------------------------------------------------------------------
New since version 0.3:
- The symbolic factorization for left-looking supernodal runs uses modern
  algorithms to run much faster.
- A bug was found and corrected in the supernodal signed Cholesky code,
  where occasionally a supernode that begins with a constrained node is
  erroneously reordered to have a different (unconstrained) first node.
  This can occasionally reduce fill-in, but as a side effect invalidates
  the elimination tree and postordering, which can cause bugs later on
  in factorization. This occurs rarely enough it seems unlikely to be
  worth pursuing as a strategy to reduce fill, and thus version 0.4 instead
  makes sure never to change the first node in a supernode when modifying
  orderings based on constraints.
- A few small API changes to make the order of function parameters a little
  more consistent; expect to see further changes in the future to make the
  library easier to use.

Edit the Makefile as needed for your platform.
Three things are created by "make all":
- fulltest, which reads the structure of a compressed-sparse-column matrix
            from stdin and writes a KKT-minimum-degree ordering back to
            stdout.
- stokes, which given a grid size as a command-line argument (default 30
          if omitted), write to stdout an example compressed-sparse-column
          matrix to stdout. (In particular, a 2D Stokes problem with
          solid wall boundary conditions around a square domain -- note
          that it has a null-space consisting of pressure values all one.)
- kktmd_mex, which is a MEX extension to MATLAB (see kktmd.m)

New in version 0.5 is the ability to use "naive" versions of the BLAS and
LAPACK routines for the case where you cannot link to working external
libraries. The file dense.h in this case has plain C++ versions of the
required routines (dgemm, dpotrf, dpotrs, dsytrf, dsytrs) specialized as
desired; these are not expected to be nearly as good as proper BLAS and
LAPACK, but may help if you have nothing else.

Also new in version 0.5 is "kktsolve.cpp", an example of how to call KKTDirect
from C++. This provides a simple single function call which factors and solves
a linear system, hiding all the complexity of the intermediate steps.