Prev Next

@(@\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 Factor and Solve with Recorded Pivoting

Syntax
int LuVecAD(
     size_t 
n,
     size_t 
m,
     VecAD<
double> &Matrix,
     VecAD<
double> &Rhs,
     VecAD<
double> &Result,
     AD<
double> &logdet)


Purpose
Solves the linear equation @[@ Matrix * Result = Rhs @]@ where Matrix is an @(@ n \times n @)@ matrix, Rhs is an @(@ n x m @)@ matrix, and Result is an @(@ n x m @)@ matrix.

The routine LuSolve uses an arbitrary vector type, instead of VecAD , to hold its elements. The pivoting operations for a ADFun object corresponding to an LuVecAD solution will change to be optimal for the matrix being factored.

It is often the case that LuSolve is faster than LuVecAD when LuSolve uses a simple vector class with elements of type double , but the corresponding ADFun objects have a fixed set of pivoting operations.

Storage Convention
The matrices stored in row major order. To be specific, if @(@ A @)@ contains the vector storage for an @(@ n x m @)@ matrix, @(@ i @)@ is between zero and @(@ n-1 @)@, and @(@ j @)@ is between zero and @(@ m-1 @)@, @[@ A_{i,j} = A[ i * m + j ] @]@ (The length of @(@ A @)@ must be equal to @(@ n * m @)@.)

n
is the number of rows in Matrix , Rhs , and Result .

m
is the number of columns in Rhs and Result . It is ok for m to be zero which is reasonable when you are only interested in the determinant of Matrix .

Matrix
On input, this is an @(@ n \times n @)@ matrix containing the variable coefficients for the equation we wish to solve. On output, the elements of Matrix have been overwritten and are not specified.

Rhs
On input, this is an @(@ n \times m @)@ matrix containing the right hand side for the equation we wish to solve. On output, the elements of Rhs have been overwritten and are not specified. If m is zero, Rhs is not used.

Result
On input, this is an @(@ n \times m @)@ matrix and the value of its elements do not matter. On output, the elements of Rhs contain the solution of the equation we wish to solve (unless the value returned by LuVecAD is equal to zero). If m is zero, Result is not used.

logdet
On input, the value of logdet does not matter. On output, it has been set to the log of the determinant of Matrix (but not quite). To be more specific, if signdet is the value returned by LuVecAD, the determinant of Matrix is given by the formula @[@ det = signdet \exp( logdet ) @]@ This enables LuVecAD to use logs of absolute values.

Example
The file lu_vec_ad_ok.cpp contains an example and test of LuVecAD. It returns true if it succeeds and false otherwise.
Input File: example/general/lu_vec_ad.cpp