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}} }@)@
Compute Determinant and Solve Linear Equations

 include <cppad/utility/lu_solve.hpp>
signdet = LuSolve(nmABXlogdet)

Use an LU factorization of the matrix A to compute its determinant and solve for X in the linear of equation @[@ A * X = B @]@ where A is an n by n matrix, X is an n by m matrix, and B is an @(@ n x m @)@ matrix.

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

Factor and Invert
This routine is an easy to user interface to LuFactor and LuInvert for computing determinants and solutions of linear equations. These separate routines should be used if one right hand side B depends on the solution corresponding to another right hand side (with the same value of A ). In this case only one call to LuFactor is required but there will be multiple calls to LuInvert.

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 ] @]@

The return value signdet is a int value that specifies the sign factor for the determinant of A . This determinant of A is zero if and only if signdet is zero.

The argument n has type size_t and specifies the number of rows in the matrices A , X , and B . The number of columns in A is also equal to n .

The argument m has type size_t and specifies the number of columns in the matrices X and B . If m is zero, only the determinant of A is computed and the matrices X and B are not used.

The argument A has the prototype
FloatVector &A
and the size of A must equal @(@ n * n @)@ (see description of FloatVector below). This is the @(@ n @)@ by n matrix that we are computing the determinant of and that defines the linear equation.

The argument B has the prototype
FloatVector &B
and the size of B must equal @(@ n * m @)@ (see description of FloatVector below). This is the @(@ n @)@ by m matrix that defines the right hand side of the linear equations. If m is zero, B is not used.

The argument X has the prototype
FloatVector &X
and the size of X must equal @(@ n * m @)@ (see description of FloatVector below). The input value of X does not matter. On output, the elements of X contain the solution of the equation we wish to solve (unless signdet is equal to zero). If m is zero, X is not used.

The argument logdet has prototype
Float &logdet
On input, the value of logdet does not matter. On output, it has been set to the log of the determinant of A (but not quite). To be more specific, the determinant of A is given by the formula
det = signdet * exp( logdet )
This enables LuSolve to use logs of absolute values in the case where Float corresponds to a real number.

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

The type FloatVector must be a SimpleVector class with elements of type Float . The routine CheckSimpleVector will generate an error message if this is not the case.

Including the file lu_solve.hpp defines the template function
     template <typename 
     bool LeqZero<
Float>(const Float &x)
in the CppAD namespace. This function returns true if x is less than or equal to zero and false otherwise. It is used by LuSolve to avoid taking the log of zero (or a negative number if Float corresponds to real numbers). This template function definition assumes that the operator <= is defined 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 LuSolve.

Complex numbers do not have the operation or <= defined. In addition, in the complex case, one can take the log of a negative number. The specializations
     bool LeqZero< std::complex<float> > (const std::complex<float> &
     bool LeqZero< std::complex<double> >(const std::complex<double> &
are defined by including lu_solve.hpp. These return true if x is zero and false otherwise.

Including the file lu_solve.hpp defines the template function
     template <typename 
     bool AbsGeq<
Float>(const Float &x, const Float &y)
If the type Float does not support the <= operation and it is not std::complex<float> or std::complex<double>, see the documentation for AbsGeq in LuFactor .

The file lu_solve.cpp contains an example and test of using this routine. It returns true if it succeeds and false otherwise.

The file lu_solve.hpp contains the current source code that implements these specifications.
Input File: cppad/utility/lu_solve.hpp