CppAD: A C++ Algorithmic Differentiation Package 20110419
lu_solve.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_LU_SOLVE_INCLUDED
00003 # define CPPAD_LU_SOLVE_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 /*
00017 $begin LuSolve$$
00018 $escape #$$
00019 $spell
00020         cppad.hpp
00021         det
00022         exp
00023         Leq
00024         typename
00025         bool
00026         const
00027         namespace
00028         std
00029         Geq
00030         Lu
00031         CppAD
00032         signdet
00033         logdet
00034 $$
00035 
00036 $index LuSolve$$
00037 $index linear, equation$$
00038 $index equation, linear$$
00039 $index determinant, Lu$$
00040 $index solve, linear equation$$
00041 
00042 $section Compute Determinant and Solve Linear Equations$$
00043 
00044 $pre
00045 $$
00046 
00047 $head Syntax$$ $code# include <cppad/lu_solve.hpp>$$
00048 $pre
00049 $$
00050 $syntax%%signdet% = LuSolve(%n%, %m%, %A%, %B%, %X%, %logdet%)%$$
00051 
00052 
00053 $head Description$$
00054 Use an LU factorization of the matrix $italic A$$ to
00055 compute its determinant 
00056 and solve for $italic X$$ in the linear of equation
00057 $latex \[
00058         A * X = B
00059 \] $$
00060 where $italic A$$ is an 
00061 $italic n$$ by $italic n$$ matrix,
00062 $italic X$$ is an 
00063 $italic n$$ by $italic m$$ matrix, and
00064 $italic B$$ is an $latex n x m$$ matrix.
00065 
00066 $head Include$$
00067 The file $code cppad/lu_solve.hpp$$ is included by $code cppad/cppad.hpp$$
00068 but it can also be included separately with out the rest of 
00069 the $code CppAD$$ routines.
00070 
00071 $head Factor and Invert$$
00072 This routine is an easy to user interface to
00073 $xref/LuFactor/$$ and $xref/LuInvert/$$ for computing determinants and
00074 solutions of linear equations.
00075 These separate routines should be used if
00076 one right hand side $italic B$$
00077 depends on the solution corresponding to another
00078 right hand side (with the same value of $italic A$$).
00079 In this case only one call to $code LuFactor$$ is required
00080 but there will be multiple calls to $code LuInvert$$.
00081 
00082 
00083 $head Matrix Storage$$
00084 All matrices are stored in row major order.
00085 To be specific, if $latex Y$$ is a vector
00086 that contains a $latex p$$ by $latex q$$ matrix,
00087 the size of $latex Y$$ must be equal to $latex  p * q $$ and for
00088 $latex i = 0 , \ldots , p-1$$,
00089 $latex j = 0 , \ldots , q-1$$,
00090 $latex \[
00091         Y_{i,j} = Y[ i * q + j ]
00092 \] $$
00093 
00094 $head signdet$$
00095 The return value $italic signdet$$ is a $code int$$ value
00096 that specifies the sign factor for the determinant of $italic A$$.
00097 This determinant of $italic A$$ is zero if and only if $italic signdet$$
00098 is zero. 
00099 
00100 $head n$$
00101 The argument $italic n$$ has type $code size_t$$ 
00102 and specifies the number of rows in the matrices
00103 $italic A$$,
00104 $italic X$$,
00105 and $italic B$$.
00106 The number of columns in $italic A$$ is also equal to $italic n$$.
00107 
00108 $head m$$
00109 The argument $italic m$$ has type $code size_t$$ 
00110 and specifies the number of columns in the matrices
00111 $italic X$$
00112 and $italic B$$.
00113 If $italic m$$ is zero,
00114 only the determinant of $italic A$$ is computed and
00115 the matrices $italic X$$ and $italic B$$ are not used.
00116 
00117 $head A$$
00118 The argument $italic A$$ has the prototype
00119 $syntax%
00120         const %FloatVector% &%A%
00121 %$$
00122 and the size of $italic A$$ must equal $latex n * n$$
00123 (see description of $xref/LuSolve/FloatVector/FloatVector/$$ below).
00124 This is the $latex n$$ by $italic n$$ matrix that 
00125 we are computing the determinant of 
00126 and that defines the linear equation.
00127 
00128 $head B$$
00129 The argument $italic B$$ has the prototype
00130 $syntax%
00131         const %FloatVector% &%B%
00132 %$$
00133 and the size of $italic B$$ must equal $latex n * m$$
00134 (see description of $xref/LuSolve/FloatVector/FloatVector/$$ below).
00135 This is the $latex n$$ by $italic m$$ matrix that 
00136 defines the right hand side of the linear equations.
00137 If $italic m$$ is zero, $italic B$$ is not used.
00138 
00139 $head X$$
00140 The argument $italic X$$ has the prototype
00141 $syntax%
00142         %FloatVector% &%X%
00143 %$$
00144 and the size of $italic X$$ must equal $latex n * m$$
00145 (see description of $xref/LuSolve/FloatVector/FloatVector/$$ below).
00146 The input value of $italic X$$ does not matter.
00147 On output, the elements of $italic X$$ contain the solution
00148 of the equation we wish to solve
00149 (unless $italic signdet$$ is equal to zero).
00150 If $italic m$$ is zero, $italic X$$ is not used.
00151 
00152 $head logdet$$
00153 The argument $italic logdet$$ has prototype
00154 $syntax%
00155         %Float% &%logdet%
00156 %$$
00157 On input, the value of $italic logdet$$ does not matter.
00158 On output, it has been set to the 
00159 log of the determinant of $italic A$$ 
00160 (but not quite).
00161 To be more specific,
00162 the determinant of $italic A$$ is given by the formula
00163 $syntax%
00164         %det% = %signdet% * exp( %logdet% )
00165 %$$
00166 This enables $code LuSolve$$ to use logs of absolute values
00167 in the case where $italic Float$$ corresponds to a real number.
00168 
00169 $head Float$$
00170 The type $italic Float$$ must satisfy the conditions
00171 for a $xref/NumericType/$$ type.
00172 The routine $xref/CheckNumericType/$$ will generate an error message
00173 if this is not the case.
00174 In addition, the following operations must be defined for any pair
00175 of $italic Float$$ objects $italic x$$ and $italic y$$:
00176 
00177 $table
00178 $bold Operation$$ $cnext $bold Description$$  $rnext
00179 $syntax%log(%x%)%$$ $cnext
00180         returns the logarithm of $italic x$$ as a $italic Float$$ object
00181 $tend
00182 
00183 $head FloatVector$$
00184 The type $italic FloatVector$$ must be a $xref/SimpleVector/$$ class with
00185 $xref/SimpleVector/Elements of Specified Type/elements of type Float/$$.
00186 The routine $xref/CheckSimpleVector/$$ will generate an error message
00187 if this is not the case.
00188 
00189 $head LeqZero$$
00190 Including the file $code lu_solve.hpp$$ defines the template function 
00191 $syntax%
00192         template <typename %Float%>
00193         bool LeqZero<%Float%>(const %Float% &%x%)
00194 %$$
00195 in the $code CppAD$$ namespace.
00196 This function returns true if $italic x$$ is less than or equal to zero
00197 and false otherwise.
00198 It is used by $code LuSolve$$ to avoid taking the log of
00199 zero (or a negative number if $italic Float$$ corresponds to real numbers).
00200 This template function definition assumes that the operator 
00201 $code <=$$ is defined for $italic Float$$ objects. 
00202 If this operator is not defined for your use of $italic Float$$,
00203 you will need to specialize this template so that it works for your
00204 use of $code LuSolve$$.
00205 $pre
00206 
00207 $$
00208 Complex numbers do not have the operation or $code <=$$ defined.
00209 In addition, in the complex case, 
00210 one can take the log of a negative number.
00211 The specializations
00212 $syntax%
00213         bool LeqZero< std::complex<float> > (const std::complex<float> &%x%)
00214         bool LeqZero< std::complex<double> >(const std::complex<double> &%x%)
00215 %$$ 
00216 are defined by including $code lu_solve.hpp$$.
00217 These return true if $italic x$$ is zero and false otherwise.
00218 
00219 $head AbsGeq$$
00220 Including the file $code lu_solve.hpp$$ defines the template function 
00221 $syntax%
00222         template <typename %Float%>
00223         bool AbsGeq<%Float%>(const %Float% &%x%, const %Float% &%y%)
00224 %$$
00225 If the type $italic Float$$ does not support the $code <=$$ operation
00226 and it is not $code std::complex<float>$$ or $code std::complex<double>$$,
00227 see the documentation for $code AbsGeq$$ in $xref/LuFactor/AbsGeq/LuFactor/$$. 
00228 
00229 $children%
00230         example/lu_solve.cpp%
00231         omh/lu_solve_hpp.omh
00232 %$$
00233 $head Example$$
00234 The file 
00235 $xref/LuSolve.cpp/$$
00236 contains an example and test of using this routine.
00237 It returns true if it succeeds and false otherwise.
00238 
00239 $head Source$$
00240 The file $cref/lu_solve.hpp/$$ contains the
00241 current source code that implements these specifications.
00242 
00243 $end
00244 --------------------------------------------------------------------------
00245 */
00246 // BEGIN PROGRAM
00247 # include <complex>
00248 # include <vector>
00249 
00250 // link exp for float and double cases
00251 # include <cppad/std_math_unary.hpp>
00252 
00253 # include <cppad/local/cppad_assert.hpp>
00254 # include <cppad/check_simple_vector.hpp>
00255 # include <cppad/check_numeric_type.hpp>
00256 # include <cppad/lu_factor.hpp>
00257 # include <cppad/lu_invert.hpp>
00258 
00259 namespace CppAD { // BEGIN CppAD namespace
00260 
00261 // LeqZero
00262 template <typename Float>
00263 inline bool LeqZero(const Float &x)
00264 {       return x <= Float(0); }
00265 inline bool LeqZero( const std::complex<double> &x )
00266 {       return x == std::complex<double>(0); }
00267 inline bool LeqZero( const std::complex<float> &x )
00268 {       return x == std::complex<float>(0); }
00269 
00270 // LuSolve
00271 template <typename Float, typename FloatVector>
00272 int LuSolve(
00273         size_t             n      ,
00274         size_t             m      , 
00275         const FloatVector &A      , 
00276         const FloatVector &B      , 
00277         FloatVector       &X      , 
00278         Float        &logdet      )
00279 {       
00280         // check numeric type specifications
00281         CheckNumericType<Float>();
00282 
00283         // check simple vector class specifications
00284         CheckSimpleVector<Float, FloatVector>();
00285 
00286         size_t        p;       // index of pivot element (diagonal of L)
00287         int     signdet;       // sign of the determinant
00288         Float     pivot;       // pivot element
00289 
00290         // the value zero
00291         const Float zero(0);
00292 
00293         // pivot row and column order in the matrix
00294         std::vector<size_t> ip(n);
00295         std::vector<size_t> jp(n);
00296 
00297         // -------------------------------------------------------
00298         CPPAD_ASSERT_KNOWN(
00299                 A.size() == n * n,
00300                 "Error in LuSolve: A must have size equal to n * n"
00301         );
00302         CPPAD_ASSERT_KNOWN(
00303                 B.size() == n * m,
00304                 "Error in LuSolve: B must have size equal to n * m"
00305         );
00306         CPPAD_ASSERT_KNOWN(
00307                 X.size() == n * m,
00308                 "Error in LuSolve: X must have size equal to n * m"
00309         );
00310         // -------------------------------------------------------
00311 
00312         // copy A so that it does not change
00313         FloatVector Lu(A);
00314 
00315         // copy B so that it does not change
00316         X = B;
00317 
00318         // Lu factor the matrix A
00319         signdet = LuFactor(ip, jp, Lu);
00320 
00321         // compute the log of the determinant
00322         logdet  = Float(0);
00323         for(p = 0; p < n; p++)
00324         {       // pivot using the max absolute element
00325                 pivot   = Lu[ ip[p] * n + jp[p] ];
00326 
00327                 // check for determinant equal to zero
00328                 if( pivot == zero )
00329                 {       // abort the mission
00330                         logdet = Float(0);
00331                         return   0;
00332                 }
00333 
00334                 // update the determinant
00335                 if( LeqZero ( pivot ) )
00336                 {       logdet += log( - pivot );
00337                         signdet = - signdet;
00338                 }
00339                 else    logdet += log( pivot );
00340 
00341         }
00342 
00343         // solve the linear equations
00344         LuInvert(ip, jp, Lu, X);
00345 
00346         // return the sign factor for the determinant
00347         return signdet;
00348 }
00349 } // END CppAD namespace 
00350 // END PROGRAM
00351 # endif