\(\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_vec_ad.cpp¶
View page sourceLu Factor and Solve with Recorded Pivoting¶
Syntax¶
int lu_vec_ad
(size_t
n ,size_t
m ,VecAD
< double > & Matrix ,VecAD
< double > & Rhs ,VecAD
< double > & Result ,double
> & logdet
)Purpose¶
Solves the linear equation
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 lu_vec_ad
solution
will change to be optimal for the matrix being factored.
It is often the case that
LuSolve
is faster than lu_vec_ad
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\),
(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 lu_vec_ad
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 lu_vec_ad
,
the determinant of Matrix is given by the formula
This enables lu_vec_ad
to use logs of absolute values.
Example¶
The file
lu_vec_ad_ok.cpp
contains an example and test of lu_vec_ad
.