\(\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}} }\)

LuSolve

View page source

Compute Determinant and Solve Linear Equations

Syntax

# include <cppad/utility/lu_solve.hpp>

signdet = LuSolve ( n , m , A , B , X , logdet )

Description

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.

Include

The file cppad/utility/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 ]\]

signdet

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.

n

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 .

m

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.

A

The argument A has the prototype

const 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.

B

The argument B has the prototype

const 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.

X

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.

logdet

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.

Float

The type Float must satisfy the conditions for a NumericType . 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

FloatVector

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.

LeqZero

Including the file lu_solve.hpp defines the template function

      template <class Float >

      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> & x )

      bool LeqZero< std::complex<double> > ( const std::complex<double> & x )

are defined by including lu_solve.hpp . These return true if x is zero and false otherwise.

AbsGeq

Including the file lu_solve.hpp defines the template function

      template <class Float >

      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 .

Example

The file lu_solve.cpp contains an example and test of using this routine.

Source

The file lu_solve.hpp contains the current source code that implements these specifications.