The LA_LUDC procedure computes the LU decomposition of an *n*-column by *m*-row array as:

*A* = *P L U*

where *P* is a permutation matrix, *L *is lower trapezoidal with unit diagonal elements (lower triangular if *n* = *m*), and *U* is upper trapezoidal (upper triangular if *n* = *m*).

LA_LUDC is based on the following LAPACK routines:

Output Type |
LAPACK Routine |

Float |
sgetrf |

Double |
dgetrf |

Complex |
cgetrf |

Double complex |
zgetrf |

## Examples

The following example uses the LU decomposition on a given array, then determines the residual error of using the resulting lower and upper arrays to recompute the original array:

`PRO ExLA_LUDC`

`; Create a random array:`

`n = 20`

`seed = 12321`

`array = RANDOMN(seed, n, n, /RAN1)`

`; Compute LU decomposition.`

`aludc = array ; make a copy`

`LA_LUDC, aludc, index`

`; Extract the lower and upper triangular arrays.`

`l = IDENTITY(n)`

`u = FLTARR(n, n)`

FOR j = 1,n - 1 DO l[0:j-1,j] = aludc[0:j-1,j]

FOR j=0,n - 1 DO u[j:*,j] = aludc[j:*,j]

`; Reconstruct array, but with rows permuted.`

arecon = l ## u

`; Adjust from LAPACK back to IDL indexing.`

`Index = Index - 1`

`; Permute the array rows back into correct order.`

`; Note that we need to loop in reverse order.`

FOR i = n - 1,0,-1 DO BEGIN & $

temp = arecon[*,i]

arecon[*, i] = arecon[*,index[i]]

arecon[*, index[i]] = temp

`ENDFOR`

PRINT, 'LA_LUDC Error:', MAX(ABS(arecon - array))

`END`

ExLA_LUDC

When this program is compiled and run, IDL prints:

LA_LUDC error: 4.76837e-007

## Syntax

LA_LUDC, *Array*, *Index* [, /DOUBLE] [, STATUS=*variable*]

## Arguments

### Array

A named variable containing the real or complex array to decompose. This procedure returns *Array* as its LU decomposition.

### Index

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

**Note: **Row numbers within *Index* start at one rather than zero.

## 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 *Array* 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:

- STATUS = 0: The computation was successful.
- STATUS > 0: One of the diagonal elements of
*U*is zero. The STATUS value specifies which value along the diagonal (starting at one) is zero.

**Note: **If STATUS is not specified, any error messages will output to the screen.

## Version History

5.6 |
Introduced |

## Resources and References

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