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].
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:
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