Prev

Next

Index-> contents reference index search external

Up-> CppAD Appendix numeric_ad LuRatio

CppAD-> Install Introduction AD ADFun preprocessor multi_thread utility ipopt_solve Example speed Appendix

Appendix-> Faq directory Theory glossary Bib wish_list whats_new deprecated compare_c numeric_ad addon License

numeric_ad-> BenderQuad opt_val_hes LuRatio

LuRatio-> lu_ratio.cpp

Headings-> Syntax Description Include Matrix Storage sign ip jp LU ---..A ---..P ---..L ---..U ---..Factor ---..Determinant ratio ---..Purpose SizeVector ADvector Example

@(@\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 and Stability Calculation
Syntax
# include <cppad/cppad.hpp>
sign = LuRatio(ip, jp, LU, ratio)

Description
Computes an LU factorization of the matrix A where A is a square matrix. A measure of the numerical stability called ratio is calculated. This ratio is useful when the results of LuRatio are used as part of an ADFun object.

Include
This routine is designed to be used with AD objects and requires the cppad/cppad.hpp file to be included.

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
     ADvector &LU
and the size of LU must equal @(@ n * n @)@ (see description of ADvector 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]]

ratio
The argument ratio has prototype
        AD<Base> &ratio
On input, the value of ratio does not matter. On output it is a measure of how good the choice of pivots is. For @(@ p = 0 , \ldots , n-1 @)@, the p-th pivot element is the element of maximum absolute value of a @(@ (n-p) \times (n-p) @)@ sub-matrix. The ratio of each element of sub-matrix divided by the pivot element is computed. The return value of ratio is the maximum absolute value of such ratios over with respect to all elements and all the pivots.

Purpose
Suppose that the execution of a call to LuRatio is recorded in the ADFun<Base> object F . Then a call to Forward of the form
     F.Forward(k, xk)
with k equal to zero will revaluate this Lu factorization with the same pivots and a new value for A . In this case, the resulting ratio may not be one. If ratio is too large (the meaning of too large is up to you), the current pivots do not yield a stable LU factorization of A . A better choice for the pivots (for this value of A ) will be made if you recreate the ADFun object starting with the Independent variable values that correspond to the vector xk .

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.

ADvector
The type ADvector must be a simple vector class with elements of type AD<Base> . The routine CheckSimpleVector will generate an error message if this is not the case.

Example
The file lu_ratio.cpp contains an example and test of using LuRatio. It returns true if it succeeds and false otherwise.

---

Input File: cppad/core/lu_ratio.hpp