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.
LA_TRIDC, AL, A, AU, U2, Index [, /DOUBLE] [, STATUS=variable]
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 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.
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.
An output vector that contains the (n - 2) elements of the second superdiagonal of the upper array.
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].
| Note |
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.
Set this keyword to a named variable that will contain the status of the computation. Possible values are:
| Note |
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
Introduced 5.6