Prev Next Index-> contents reference index search external Up-> CppAD utility LuDetAndSolve LuFactor CppAD-> Install Introduction AD ADFun preprocessor multi_thread utility ipopt_solve Example speed Appendix utility-> ErrorHandler NearEqual speed_test SpeedTest time_test test_boolofvoid NumericType CheckNumericType SimpleVector CheckSimpleVector nan pow_int Poly LuDetAndSolve RombergOne RombergMul Runge45 Rosen34 OdeErrControl OdeGear OdeGearControl CppAD_vector thread_alloc index_sort to_string set_union sparse_rc sparse_rcv LuDetAndSolve-> LuSolve LuFactor LuInvert LuFactor-> lu_factor.cpp lu_factor.hpp Headings-> Syntax Description Include Matrix Storage sign ip jp LU ---..A ---..P ---..L ---..U ---..Factor ---..Determinant SizeVector FloatVector Float AbsGeq Example Source

$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$
LU Factorization of A Square Matrix
 Syntax
 include <cppad/utility/lu_factor.hpp>  sign = LuFactor(ip, jp, LU)

Description
Computes an LU factorization of the matrix A where A is a square matrix.

Include
The file cppad/lu_factor.hpp is included by cppad/cppad.hpp but it can also be included separately with out the rest of the CppAD routines.

Matrix Storage
All matrices are stored in row major order. To be specific, if $Y$ is a vector that contains a $p$ by $q$ matrix, the size of $Y$ must be equal to $p * q$ and for $i = 0 , \ldots , p-1$, $j = 0 , \ldots , q-1$, $$Y_{i,j} = Y[ i * q + j ]$$

sign
The return value sign has prototype       int sign  If A is invertible, sign is plus or minus one and is the sign of the permutation corresponding to the row ordering ip and column ordering jp . If A is not invertible, sign is zero.

ip
The argument ip has prototype       SizeVector &ip  (see description of SizeVector below). The size of ip is referred to as n in the specifications below. The input value of the elements of ip does not matter. The output value of the elements of ip determine the order of the rows in the permuted matrix.

jp
The argument jp has prototype       SizeVector &jp  (see description of SizeVector below). The size of jp must be equal to n . The input value of the elements of jp does not matter. The output value of the elements of jp determine the order of the columns in the permuted matrix.

LU
The argument LU has the prototype       FloatVector &LU  and the size of LU must equal $n * n$ (see description of FloatVector below).

A
We define A as the matrix corresponding to the input value of LU .

P
We define the permuted matrix P in terms of A by       P(i, j) = A[ ip[i] * n + jp[j] ] 
L
We define the lower triangular matrix L in terms of the output value of LU . The matrix L is zero above the diagonal and the rest of the elements are defined by       L(i, j) = LU[ ip[i] * n + jp[j] ]  for $i = 0 , \ldots , n-1$ and $j = 0 , \ldots , i$.

U
We define the upper triangular matrix U in terms of the output value of LU . The matrix U is zero below the diagonal, one on the diagonal, and the rest of the elements are defined by       U(i, j) = LU[ ip[i] * n + jp[j] ]  for $i = 0 , \ldots , n-2$ and $j = i+1 , \ldots , n-1$.

Factor
If the return value sign is non-zero,       L * U = P  If the return value of sign is zero, the contents of L and U are not defined.

Determinant
If the return value sign is zero, the determinant of A is zero. If sign is non-zero, using the output value of LU the determinant of the matrix A is equal to  sign * LU[ip[0], jp[0]] * ... * LU[ip[n-1], jp[n-1]] 
SizeVector
The type SizeVector must be a SimpleVector class with elements of type size_t . The routine CheckSimpleVector will generate an error message if this is not the case.

FloatVector
The type FloatVector must be a simple vector class . The routine CheckSimpleVector will generate an error message if this is not the case.

Float
This notation is used to denote the type corresponding to the elements of a FloatVector . The type Float must satisfy the conditions for a NumericType type. The routine CheckNumericType will generate an error message if this is not the case. In addition, the following operations must be defined for any pair of Float objects x and y :
 Operation Description log(x) returns the logarithm of x as a Float object

AbsGeq
Including the file lu_factor.hpp defines the template function       template <typename Float>      bool AbsGeq<Float>(const Float &x, const Float &y)  in the CppAD namespace. This function returns true if the absolute value of x is greater than or equal the absolute value of y . It is used by LuFactor to choose the pivot elements. This template function definition uses the operator <= to obtain the absolute value for Float objects. If this operator is not defined for your use of Float , you will need to specialize this template so that it works for your use of LuFactor.  Complex numbers do not have the operation <= defined. The specializations  bool AbsGeq< std::complex<float> >      (const std::complex<float> &x, const std::complex<float> &y) bool AbsGeq< std::complex<double> >      (const std::complex<double> &x, const std::complex<double> &y)  are define by including lu_factor.hpp These return true if the sum of the square of the real and imaginary parts of x is greater than or equal the sum of the square of the real and imaginary parts of y .

Example
The file lu_factor.cpp contains an example and test of using LuFactor by itself. It returns true if it succeeds and false otherwise.  The file lu_solve.hpp provides a useful example usage of LuFactor with LuInvert.

Source
The file lu_factor.hpp contains the current source code that implements these specifications.