--------------------------------------------------
lines 7-132 of file: example/general/lu_vec_ad.cpp
--------------------------------------------------

{xrst_begin lu_vec_ad.cpp}
{xrst_spell
   logdet
   rhs
   signdet
}

Lu Factor and Solve with Recorded Pivoting
##########################################

Syntax
******
| ``int lu_vec_ad`` (
| |tab| ``size_t`` *n* ,
| |tab| ``size_t`` *m* ,
| |tab| ``VecAD`` < *double* > & *Matrix* ,
| |tab| ``VecAD`` < *double* > & *Rhs* ,
| |tab| ``VecAD`` < *double* > & *Result* ,
| |tab| *AD* < ``double`` > & ``logdet`` )

Purpose
*******
Solves the linear equation

.. math::

   Matrix * Result = Rhs

where *Matrix* is an :math:`n \times n` matrix,
*Rhs* is an :math:`n x m` matrix, and
*Result* is an :math:`n x m` matrix.

The routine :ref:`LuSolve-name` uses an arbitrary vector type,
instead of :ref:`VecAD-name` ,
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
:ref:`elements of type double<SimpleVector@Elements of Specified Type>` ,
but the corresponding :ref:`ADFun-name` objects have a fixed
set of pivoting operations.

Storage Convention
******************
The matrices stored in row major order.
To be specific, if :math:`A` contains the vector storage for an
:math:`n x m` matrix,
:math:`i` is between zero and :math:`n-1`,
and :math:`j` is between zero and :math:`m-1`,

.. math::

   A_{i,j} = A[ i * m + j ]

(The length of :math:`A` must be equal to :math:`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
:math:`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
:math:`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
:math:`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

.. math::

   det = signdet \exp( logdet )

This enables ``lu_vec_ad`` to use logs of absolute values.

Example
*******
{xrst_toc_hidden
   example/general/lu_vec_ad_ok.cpp
}
The file
:ref:`lu_vec_ad_ok.cpp-name`
contains an example and test of ``lu_vec_ad`` .

{xrst_end lu_vec_ad.cpp}