[Rd] Lapack, determinant, multivariate normal density, solution to linear system, C language

[shotwelm]

r-devel list,

I have recently written an R package that solves a linear least squares
problem, and computes the multivariate normal density function. The bulk
of the code is written in C, with interfacing code to the BLAS and
Lapack libraries. The motivation here is speed. I ran into a problem
computing the determinant of a symmetric matrix in packed storage.
Apparently, there are no explicit routines for this as part of Lapack.
While there IS an explicit routine for this in Linpack, I did not want
to use the older library. Also, right before I needed the determinant of
the matrix A, I had used the Lapack routine dspsv to solve the linear
system Ax=b, which does much of the work of computing a determinant
also. In fact, the solution I came up with involves the output of this
routine (which might be obvious to Lapack designers, but not me)

My modest Googleing turned up very little unique material (as is typical
with BLAS/Lapack/Linpack queries). Hence, I am writing the r-devel list
partly to document the solution I've come up with, but mainly to elicit
additional wisdom from seasoned R programmers.

My solution to the problem is illustrated in the appended discussion and
C code. Thanks for your input.

-Matt Shotwell

--------------

The Lapack routine dspsv solves the linear system of equations Ax=b,
where A is a symmetric matrix in packed storage format. The dspsv
function performs the factorization A=UDU', where U is a unitriangular
matrix and D is a block diagonal matrix where the blocks are of
dimension 1x1 or 2x2. In addition to the solution for x, the dspsv
function also returns the matrices U and D. The matrix D may then be
used to compute the determinant of A. Recall from linear algebra that
det(A) = det(UDU') = det(U)det(D)det(U'). Since U is unitriangular,
det(U) = 1. The determinant of D is the product of the determinants of
the diagonal blocks. If a diagonal block is of dimension 1x1, then the
determinant of the block is simply the value of the single element in
the block. If the diagonal block is of dimension 2x2 then the
determinant of the block may be computed according to the well-known
formula b11*b22-b12*b21, where bij is the value in the i'th row and j'th
column of the block.

  int i, q, info, *ipiv, one = 1;
  double *b, *A, *D, det;

  /* 
  ** A and D are upper triangular matrices in packed storage
  ** A[] = a00, a01, a11, a02, a12, a22, a03, a13, a23, a33, ...
  ** use the following macro to address the element in the 
  ** i'th row and j'th column for i <= j
  */
  #define UMAT(i, j) (i + j * ( j + 1 ) / 2)

  /* 
  ** additional code should be here
  ** - set q
  ** - allocate ipiv...
  ** - allocate and fill A and b... 
  */

  /* 
  ** solve Ax=b using A=UDU' factorization where 
  ** A represents a qxq matrix, b a 1xq vector.
  ** dspsv outputs the elements of the matrix D
  ** is upper triangular packed storage
  ** in the memory addressed by A. That is, A is
  ** replaced by D when dspsv returns.
  */
  F77_CALL(dspsv)("U", &q, &one, A, ipiv, b, &q, &info); 
  if( info > 0 ) { /*issue warning, system is singular*/ }
  if( info < 0 ) { /*issue error, invalid argument*/ }

  /* 
  ** compute the determinant det = det(A)
  ** if ipiv[i] > 0, then D(i,i) is a 1x1 block diagonal
  ** if ipiv[i] = ipiv[i-1] < 0, then D(i-1,i-1), 
  ** D(i-1,i), and D(i,i) form the upper triangle 
  ** of a 2x2 block diagonal
  */
  D = A;
  det = 1.0;
  for( i = 0; i < q; i++ ) { 
    if( ipiv[ i ] > 0 ) {
      det *= D[ UMAT(i,i) ];
    } else if( i > 0 && ipiv[ i ] < 0 && ipiv[ i-1 ] == ipiv[ i ] ) { 
      det *= D[ UMAT(i,i) ] * D[ UMAT(i-1,i-1) ] -\
             D[ UMAT(i-1,i) ] * D[ UMAT(i-1,i) ];
    }
  }