IDL Reference Guide: Procedures and Functions

LA_TRIDC

The LA_TRIDC procedure computes the LU decomposition of a tridiagonal (n x n) array as Array = L U, where L is a product of permutation and unit lower bidiagonal arrays, and U is upper triangular with nonzero elements only in the main diagonal and the first two superdiagonals.

LA_TRIDC is based on the following LAPACK routines:

For more details, see Anderson et al., LAPACK Users' Guide, 3rd ed., SIAM, 1999.

Syntax

LA_TRIDC, AL, A, AU, U2, Index [, /DOUBLE] [, STATUS=variable]

Arguments

AL

A named vector of length (n - 1) containing the subdiagonal elements of an array. This procedure returns AL as the (n - 1) elements of the lower bidiagonal array from the LU decomposition.

A

A named vector of length n containing the main diagonal elements of an array. This procedure returns A as the n diagonal elements of the upper array from the LU decomposition.

AU

A named vector of length (n - 1) containing the superdiagonal elements of an array. This procedure returns AU as the (n - 1) superdiagonal elements of the upper array.

U2

An output vector that contains the (n - 2) elements of the second superdiagonal of the upper array.

Index

An output vector that records the row permutations which occurred as a result of partial pivoting. For 1 < j < n, row j of the matrix was interchanged with row Index[j].

Keywords

DOUBLE

Set this keyword to use double-precision for computations and to return a double-precision (real or complex) result. Set DOUBLE = 0 to use single-precision for computations and to return a single-precision (real or complex) result. The default is /DOUBLE if AL is double precision, otherwise the default is DOUBLE = 0.

STATUS

Set this keyword to a named variable that will contain the status of the computation. Possible values are:

Examples

Create a test program to compute the LU decomposition of a tridiagonal array:

PRO EX_LA_TRIDC  
    ; Create a random tridiagonal array.  
    n = 9  
    seed = 12321  
    AL = RANDOMN(seed, n-1)  
    A  = RANDOMN(seed, n)  
    AU = RANDOMN(seed, n-1)  
  
    ; Construct tridiagonal array.  
    Array = DIAG_MATRIX(AL, -1) + DIAG_MATRIX(A) + $  
        DIAG_MATRIX(AU, 1)  
  
    ; Compute the LU decomposition.  
    LA_TRIDC, AL, A, AU, U2, Index  
  
    ; Adjust from LAPACK back to IDL indexing.  
    Index = Index - 1  
  
    ; Create upper and lower arrays.  
    Upper = DIAG_MATRIX(A) + $  
        DIAG_MATRIX(AU, 1) + DIAG_MATRIX(U2, 2)  
    Lower = DIAG_MATRIX(AL, -1) + IDENTITY(n)  
  
    ; To conserve storage, LA_TRIDC keeps all lower diagonal  
    ; elements in AL, regardless of row. The Index array  
    ; tells which subdiagonals need to be shifted down.  
    ; Loop starts at 1 since there aren't any subdiagonals  
    ; to the left of the first diagonal element.  
    FOR i = 1,n-2 DO BEGIN  
        IF (Index[i] NE i) THEN $  
            Lower[0:i-1,[i,i+1]] = Lower[0:i-1,[i+1,i]]  
    ENDFOR  
  
    ; Permute the row order.  
    FOR i = n-2, 0, -1 DO BEGIN  
        IF (Index[i] NE i) THEN $  
            Lower[*,[i,i+1]] = Lower[*,[i+1,i]]  
    ENDFOR  
  
    ; Reconstruct the array and check the difference:  
    Arecon = Lower ## Upper  
    PRINT, 'LA_TRIDC error:', MAX(ABS(Arecon - Array))  
END  

When this program is compiled and run, IDL prints:

LA_TRIDC error: 1.50427e-008