Solve and Inverse with Symmetric Indefinite Decomposition
DESCRIPTION:
Given the symmetric indefinite decomposition of a real symmetric or complex Hermitian matrix, either solves a system of linear equations with that matrix as coefficient matrix or else computes the inverse of the matrix.
USAGE: 
solve.lu.Hermitian(a, b, tol=0, norma)
REQUIRED ARGUMENTS:
a:
An object of class lu.Hermitian, representing the symmetric indefinite decomposition of a real symmetric or complex Hermitian matrix.
OPTIONAL ARGUMENTS:
b:
A real or complex matrix or vector. The number of rows in b must equal the dimension of the matrix underlying a.
tol:
Tolerance for reciprocal condition estimation. If tol is negative, no condition estimation is done. Otherwise, the reciprocal one/infinity norm condition estimate is computed and the solve or inverse computation is done only if the condition estimate is greater than tol. By default, tol = 0.
norma:
The one/infinity norm of the matrix for use in condition estimation. By default, it is assumed that the norm is available as an attribute of a, since the default for lu.Hermitian is to return one/infinity norm of the matrix as an attribute.
VALUE:
If A is the matrix whose symmetric-indefinite decomposition is represented by a, an object of class "Matrix" is returned that is the solution x to the system of equations A %*% x = b If b is not supplied, the inverse of A is returned. Attributes include a copy of the call to solve, and the one/infinity norm reciprocal condition estimate if tol is nonnegative.
DETAILS:
Based on the functions dsycon, dsytri, zhecon, zhetri from LAPACK (Anderson et al. 1994).
REFERENCES:
Anderson, E., et al. (1994). LAPACK User's Guide, 2nd edition, SIAM, Philadelphia.
SEE ALSO:
rcond.Hermitian , solve , solve.Hermitian .
EXAMPLES: 
 n <- 5
 a <- Matrix( rnorm(n*n), nrow = n, ncol = n)
 a[row(a) > col(a)] <- t(a)[row(a) > col(a)]  # construct symmetric matrix
 class(a) <- Matrix.class(x)
 b <- rnorm(n)
 z <- lu(a)                                   # symmetric-indefinite decomp
 a %*% solve(a,b) - b                         # residual
(solve(a) %*% b) - solve(a,b)