**BioStatMatt » R**, and kindly contributed to R-bloggers)

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

**leave a comment**for the author, please follow the link and comment on his blog:

**BioStatMatt » R**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...