C89 Cholesky decomposition using LAPACK

lapack_cholesky.c

/**
  @file lapack_cholesky.c
  @author Pedro Rittner
  @date November 25, 2014
  @brief Small example for how to perform a Cholesky decomposition using LAPACK

  A small example performing a Cholesky decomposition using LAPACK's C API.
  This is meant to demystify the API somewhat and to make it clearer what
  is happening in terms of the function call/arguments/etc.

  Obviously this depends on GCC, GFortran, and LAPACK being installed.

  If you're using Ubuntu 14.04 or later, install the following packages:
  - gfortran
  - liblapack-dev
  - liblapacke-dev
  - libblas-dev
  - libblas3gf
  - libopenblas-dev

  To build using GCC from the shell:

  $ gcc lapack_cholesky.c -llapack
  
  Add flags for flavoring, of course. Generally I would prefer to write this
  using C99 conventions but have written it so that ANSI-C/ISO-C89 works.

  A word on LDA - Leading Dimension of Array:

  Suppose that you have a matrix A of size 100x100 which is stored in an array
  100x100. In this case LDA is the same as N. Now suppose that you want to work
  only on the submatrix A(91:100 , 1:100); in this case the # of rows is 10
  but LDA=100. Assuming the fortran column-major ordering (which is the case in
  LAPACK), the LDA is used to define the distance in memory between elements of
  two consecutive columns which have the same row index. If you call
  B = A(91:100 , 1:100) then B(1,1) and B(1,2) are 100 memory locations
  far from each other. 

  From: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=217

 */

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

/* If the LAPACK library and/or headers are not in your compiler's include path,
   make sure to use the compiler's -I and -L options to specify where your
   compiler should look for them.*/
#include <lapacke.h>

/* Prints matrix in column-major format. */
static void show_matrix(const double* A, const int n)
{
    int i = 0, j = 0;
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            printf("%2.5f ", A[i * n + j]);
        }
        printf("\n");
    }
}

/* Zeros out the lower diagonal portion of a matrix. */
static void zero_ldiag(double* m, const int n)
{
    int i = 0, j = 0;
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            if (i < j)
            {
                m[j * n + i] = 0.0;
            }
        }
    }
}

int main(void)
{
    int n = 3;       /* Matrix dimensions (since this is a square matrix) */
    char uplo = 'L'; /* We ask LAPACK for the lower diagonal matrix L */
    int info = 0;    /* "Info" return value, used for error-checking */
    int lda = n;     /* Leading dimension of array - see above */

    /* Column-major format; doesn't matter since it's positive-definitive */
    double m1[] = {25, 15, -5,
                   15, 18,  0,
                   -5,  0, 11};
    
    /* LAPACK API uses updating of input variables instead of returning */
    dpotrf_(&uplo, &n, m1, &lda, &info);
    
    /* Will be 0 IFF call succeeds */
    assert(info == 0);
    
    /* Note that because m1 is over-written byy LAPACK, you need to 0 out the 
       lower diagonal entries yourself, since LAPACK will not do it for you. */
    zero_ldiag(m1, n);
    
    /* Should print this:
       5.00000 3.00000 -1.00000 
       0.00000 3.00000 1.00000
       0.00000 0.00000 3.00000 
    */
    show_matrix(m1, n);
    
    return 0;
}

Very nice code! Just one question. If you "ask LAPACK for the lower diagonal matrix L", why matrix showed by zero_ldiag(m1, n) is an upper diagonal matrix?

On my system (Mac os) dpotrf_() does not work. Instead:
LAPACKE_dpotrf() works but you need to change arguments. Instead of addresses you give values. Like that:
LAPACKE_dpotrf(LAPACK_COL_MAJOR, uplo, n, m1, lda);
It now produces the correct output.