The Lapack routine dspsv solves the linear system of equations Ax=b, where A is a symmetric matrix in packed storage format. However, there appear to be no Lapack functions that compute the determinant of such a matrix. We need to compute the determinant, for instance, in order to compute the multivariate normal density function. 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 1×1 or 2×2. 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 1×1, then the determinant of the block is simply the value of the single element in the block. If the diagonal block is of dimension 2×2 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. The following C code snip demonstrates the procedure.

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, D is placed in A
** where A represents a qxq matrix, b a 1xq vector
*/
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 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) ];
  }
}