MF_LUdecompose MD_LUdecompose ME_LUdecompose
MCF_LUdecompose MCD_LUdecompose MCE_LUdecompose
MF_LUdecomposewEdit MD_LUdecomposewEdit ME_LUdecomposewEdit
MCF_LUdecomposewEdit MCD_LUdecomposewEdit MCE_LUdecomposewEdit
MFb_LUdecompose MDb_LUdecompose MEb_LUdecompose
MCFb_LUdecompose MCDb_LUdecompose MCEb_LUdecompose
MFb_LUdecomposewEdit MDb_LUdecomposewEdit MEb_LUdecomposewEdit
MCFb_LUdecomposewEdit MCDb_LUdecomposewEdit MCEb_LUdecomposewEdit
MFb_LUdecompose_sizeBuf MDb_LUdecompose_sizeBuf MEb_LUdecompose_sizeBuf
MCFb_LUdecompose_sizeBuf MCDb_LUdecompose_sizeBuf MCEb_LUdecompose_sizeBuf
FunctionLU decomposition
Syntax C/C++#include <MFstd.h>
int MF_LUdecompose( fMatrix MLU, uiVector Ind, fMatrix MA, ui len );
int MF_LUdecomposewEdit( fMatrix MLU, uiVector Ind, fMatrix MA, ui len, float thresh );
int MCF_LUdecomposewEdit( cfMatrix MLU, uiVector Ind, cfMatrix MA, ui len, float thresh );
int MFb_LUdecompose( fMatrix MLU, uiVector Ind, fMatrix MA, ui len, fVector Buf );
int MFb_LUdecomposewEdit( fMatrix MLU, uiVector Ind, fMatrix MA, ui len, float thresh, fVector Buf );
int MCFb_LUdecomposewEdit( cfMatrix MLU, uiVector Ind, cfMatrix MA, ui len, float thresh, cfVector Buf );
ui MFb_LUdecompose_sizeBuf( ui len );
C++ MatObj#include <OptiVec.h>
void matrix<T>::LUdecompose( vector<ui> Ind, const matrix<T>& MA );
void matrix<T>::LUdecompose( vector<ui>* Ind, const matrix<T>& MA );
void matrix<T>::LUdecomposewEdit( vector<ui> Ind, const matrix<T>& MA, const T& thresh );
void matrix<T>::LUdecomposewEdit( vector<ui>* Ind, const matrix<T>& MA, const T& thresh );
void matrix<complex<T>>::LUdecomposewEdit( vector<ui>* Ind, const matrix<complex<T>>& MA, const T& thresh );
void matrix<T>::b_LUdecompose( vector<ui> Ind, const matrix<T>& MA, vector<T> Buf );
void matrix<T>::b_LUdecompose( vector<ui>* Ind, const matrix<T>& MA, vector<T> Buf );
void matrix<T>::b_LUdecomposewEdit( vector<ui> Ind, const matrix<T>& MA, const T& thresh, vector<T> Buf );
void matrix<T>::b_LUdecomposewEdit( vector<ui>* Ind, const matrix<T>& MA, const T& thresh, vector<T> Buf );
void matrix<complex<T>>::b_LUdecomposewEdit( vector<ui>* Ind, const matrix<complex<T>>& MA, const T& thresh, vector<complex<T>> Buf );
ui matrix<T>::b_LUdecompose_sizeBuf();
Pascal/Delphiuses MFstd;
function MF_LUdecompose( MLU:fMatrix; Ind:uiVector; MA:fMatrix; len:UIntSize ):Integer;
function MF_LUdecomposewEdit( MLU:fMatrix; Ind:uiVector; MA:fMatrix; len:UIntSize; thresh:Single ):Integer;
function MCF_LUdecomposewEdit( MLU:cfMatrix; Ind:uiVector; MA:cfMatrix; len:UIntSize; thresh:Single ):Integer;
function MFb_LUdecompose( MLU:fMatrix; Ind:uiVector; MA:fMatrix; len:UIntSize; Buf:fVector ):Integer;
function MFb_LUdecomposewEdit( MLU:fMatrix; Ind:uiVector; MA:fMatrix; len:UIntSize; thresh:Single; Buf:fVector ):Integer;
function MCFb_LUdecomposewEdit( MLU:cfMatrix; Ind:uiVector; MA:cfMatrix; len:UIntSize; thresh:Single; Buf:cfVector ):Integer;
function MFb_LUdecompose_sizeBuf( len:UIntSize ):UIntSize;
DescriptionMA is decomposed into a product MA = L * U, where L is lower-triangular with the diagonal elements equal to 1, and U is upper-triangular. As the combined number of non-trivial elements of L and U just fit into a matrix of the same dimensions as MA, the result is stored in a single matrix MLU rather than in two distinct matrices L and U. Actually, it is not the "true" matrices L and U which are combined into MLU, but rather a row-wise permutation, whose indices are output in the vector Ind.

MA may or may not be overwritten by MLU.

There are applications where it makes sense to make near-singular matrices decomposable by "pivot editing": If, in the partial pivoting process employed in the decomposition and inversion, no diagonal element larger than the threshold is available, the scaling of the matrix is done with the threshold value rather than with the tiny diagonal element. That way, divisions by (near-)zero are avoided.
By default, this pivot editing is switched off. If you wish to switch it on for all calls to MF_LUdecompose and to the functions based upon it, namely MF_inv and MF_solve, you can do so by calling MF_LUDsetEdit. However, as this method is not thread-safe, you cannot use it in order to set different thresholds for different calls to the functions mentioned. Instead, use their "wEdit" variants, i.e. MF_LUdecomposewEdit, MF_invwEdit and MF_solvewEdit. They take the desired threshold as the additional argument thresh. Note that thresh is always real, also in the complex versions.

The return value indicates if the factorization was successful and if the number of permuations was even or uneven:
Return valueMeaning
0Matrix is singular; result cannot be used
+1Factorization successful; even number of permutations
-1Factorization successful; uneven number of permutations
+2Matrix (nearly) singular; could be factorized by pivot editing; result partially usable; even number of permutations
-2Matrix (nearly) singular; could be factorized by pivot editing; result partially usable; uneven number of permutations

To check if MF_LUdecompose was successful, in single-thread programs, you may also call MF_LUDresult, whose return value will be FALSE (0), if the MA could be decomposed, and TRUE (1) for singular MA. In multi-thread programs, on the other hand, it would not be clear wich instance of MF_LUdecompose the call to MF_LUDresult would refer to. So, here, inspection of the return value of MF_LUdecompose is the only option.

The LU decomposition functions need buffer memory. The "normal" versions (prefixes MF_, MCF_ etc.) allocate it themselves, whereas the version with the prefixes MFb_, MCFb_ etc. take a vector Buf as an additional argument. The required size of Buf can be obtained by calling MFb_LUdecompose_sizeBuf() etc., whereby the size is given as the number of elements of Buf in the respective data type (rather than in bytes).

Error handlingIn the case of a singular matrix, MLU will be undefined, and failure is indicated in the return value. Additionally, an internal flag is set, which can be read by calling MF_LUDresult.
Return valueCode 0, +-1 or +-2, indicating success, even or uneven number of row permutations (see above)
See alsochapter 10

MatrixLib Table of Contents  OptiVec home