IDL Analyst Reference Guide: Linear Systems

IMSL_SP_BDFAC

The IMSL_SP_BDFAC procedure computes the LU factorization of a matrix stored in band storage mode.

Syntax

IMSL_SP_BDFAC, nlca, nuca, n_rows, a, pivot, factor [, BLK_FACTOR=value] [, CONDITION=variable] [, /DOUBLE]

Arguments

a

Array of size (nlca + nuca + 1) x n containing the n x n banded coefficient matrix in band storage mode A(i,j). See Band Storage Format for a description of band storage mode.

factor

A named variable that will contain an array of size (2*nlca + nuca + 1) x n_rows containing the LU factorization of A with column pivoting. The keywords FACTOR and CONDITION cannot be used together.

n_rows

Number of rows in a.

nlca

Number of lower codiagonals in a.

nuca

Number of upper codiagonals in a.

pivot

A named variable that will contain a one-dimensional array containing the pivot sequence. The keywords PIVOT and CONDITION cannot be used together.

Keywords

BLK_FACTOR

The blocking factor. This keyword must be set no larger than 32. Default: BLK_FACTOR = 1.

CONDITION

Named variable into which an estimate of the L₁ condition number is stored. The keyword CONDITION cannot be used with arguments pivot or factor.

DOUBLE

If present and nonzero, double precision is used.

Discussion

The IMSL_SP_BDFAC function computes the LU factorization of A with based on the blocked LU factorization algorithm given in Du Croz, et al, (1990). Level-3 BLAS invocations were replaced by in-line loops. The blocking factor BLK_FACTOR has the default value of 1, but can be reset to any positive value not exceeding 32.

An estimate of the L₁ 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/e (where e 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.

Example

Consider the 1000 x 1000 banded matrix below:

This example computes the solution to Ax₁ = b₁ and Ax₂ = b₂, where b₁ and b₂ are random vectors. The factorization is computed just once, using IMSL_SP_BDFAC, and the solutions are computed using IMSL_SP_BDSOL.

n_rows = 1000L  
nlca = 1L  
nuca = 1L  
a = DBLARR(n_rows*(nlca+nuca+1))  
a(1:n_rows-1) = 4  
a(n_rows:2*n_rows-1) = -1  
a(2*n_rows:*) = 4  
; Fill A with the values of the bands.  
seed = 123L  
b1 = RANDOMU(seed, n_rows)  
b2 = RANDOMU(seed, n_rows)  
; Fill random vectors  
IMSL_SP_BDFAC, nlca, nuca, n_rows, a, pivot, factor  
; Compute the factorization using IMSL_SP_BDFAC.  
x1 = IMSL_SP_BDSOL(b1, nlca, nuca, Factor = factor, Pivot = pivot)  
; Compute solution of Ax1 = b1 above, and output residual below.  
PM, TOTAL(ABS(IMSL_SP_MVMUL(n_rows, n_rows, nlca, nuca, $  
   a, x1)-b1))  
  
   1.2367884e-13  
  
x2 = IMSL_SP_BDSOL(b2, nlca, nuca, Factor = factor, Pivot = pivot)  
; Compute the solution of Ax2 = b2 above, and output residual.  
PM, TOTAL(ABS(IMSL_SP_MVMUL(n_rows, n_rows, nlca, nuca, $  
   a, x2)-b2))  
  
   9.1537888e-14