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 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. 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) ];
  }
}