SP_BDPDFAC Function
Computes the RTR Cholesky factorization of symmetric positive definite matrix, A, in band symmetric storage mode.
Usage
result = SP_BDPDFAC(a, n, ncoda)
Input Parameters
a—Array of size (ncoda + 1) × n containing the n × n banded coefficient matrix in band symmetric storage mode A(i,j). See the chapter introduction for a description of band symmetric storage mode.
n—Number rows in a.
ncoda—Number of upper codiagonals in a.
Returned Value
result—An array of size (ncoda + 1) × n containing the RTR factorization of A in band symmetric storage mode.
Input Keywords
Double—If present and nonzero, double precision is used.
Output Keywords
Condition—Specifies a named variable into which an estimate of the L1 condition number is stored.
Discussion
The function SP_BDPDFAC computes the RTR Cholesky factorization of A. R is an upper triangular band matrix.
The L1 condition number of A is computed using Higham’s modifications to Hager’s method, as given in Higham (1988). If the estimated condition number is greater than 1/ε (where ε is the machine precision), a warning message is issued. This indicates that very small changes in A may produce large changes in the solution x.
The function SP_BDPDFAC fails if any submatrix of R is not positive definite or if R has a zero diagonal element. These errors occur only if A is very close to a singular matrix or to a matrix which is not positive definite.
The function SP_BDPDFAC is partially based on the LINPACK subroutines CPBFA and SPBSL; see Dongarra et al. (1979).
Example 1
Solve a system of linear equations Ax = b, using both SP_BDPDFAC and SP_BDPDSOL, where:
, 
n = 4L
ncoda = 2L
; Define A in band symmetric storage mode.
a = DBLARR((ncoda+1)*n)
a(0:n-1) = [0, 0, -1, 1]
a(n:2L*n-1) = [0, 0, 2, -1]
a(2L*n:*) = [2, 4, 7, 3]
b = [6, -11, -11, 19]
; Use SP_BDPDFAC to compute the factorization.
factor = SP_BDPDFAC(a, n, ncoda)
; Compute the solution
x = SP_BDPDSOL(b, ncoda, Factor=factor)
PM, x ; PV-WAVE prints the following:
; 4.0000000
; -6.0000000
; 2.0000000
; 9.0000000
