CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_LU_SOLVE_INCLUDED 00003 # define CPPAD_LU_SOLVE_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell 00007 00008 CppAD is distributed under multiple licenses. This distribution is under 00009 the terms of the 00010 Eclipse 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 $icode%signdet% = LuSolve(%n%, %m%, %A%, %B%, %X%, %logdet%)%$$ 00051 00052 00053 $head Description$$ 00054 Use an LU factorization of the matrix $icode A$$ to 00055 compute its determinant 00056 and solve for $icode X$$ in the linear of equation 00057 $latex \[ 00058 A * X = B 00059 \] $$ 00060 where $icode A$$ is an 00061 $icode n$$ by $icode n$$ matrix, 00062 $icode X$$ is an 00063 $icode n$$ by $icode m$$ matrix, and 00064 $icode 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 $cref LuFactor$$ and $cref LuInvert$$ for computing determinants and 00074 solutions of linear equations. 00075 These separate routines should be used if 00076 one right hand side $icode B$$ 00077 depends on the solution corresponding to another 00078 right hand side (with the same value of $icode 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 $icode signdet$$ is a $code int$$ value 00096 that specifies the sign factor for the determinant of $icode A$$. 00097 This determinant of $icode A$$ is zero if and only if $icode signdet$$ 00098 is zero. 00099 00100 $head n$$ 00101 The argument $icode n$$ has type $code size_t$$ 00102 and specifies the number of rows in the matrices 00103 $icode A$$, 00104 $icode X$$, 00105 and $icode B$$. 00106 The number of columns in $icode A$$ is also equal to $icode n$$. 00107 00108 $head m$$ 00109 The argument $icode m$$ has type $code size_t$$ 00110 and specifies the number of columns in the matrices 00111 $icode X$$ 00112 and $icode B$$. 00113 If $icode m$$ is zero, 00114 only the determinant of $icode A$$ is computed and 00115 the matrices $icode X$$ and $icode B$$ are not used. 00116 00117 $head A$$ 00118 The argument $icode A$$ has the prototype 00119 $codei% 00120 const %FloatVector% &%A% 00121 %$$ 00122 and the size of $icode A$$ must equal $latex n * n$$ 00123 (see description of $cref/FloatVector/LuSolve/FloatVector/$$ below). 00124 This is the $latex n$$ by $icode n$$ matrix that 00125 we are computing the determinant of 00126 and that defines the linear equation. 00127 00128 $head B$$ 00129 The argument $icode B$$ has the prototype 00130 $codei% 00131 const %FloatVector% &%B% 00132 %$$ 00133 and the size of $icode B$$ must equal $latex n * m$$ 00134 (see description of $cref/FloatVector/LuSolve/FloatVector/$$ below). 00135 This is the $latex n$$ by $icode m$$ matrix that 00136 defines the right hand side of the linear equations. 00137 If $icode m$$ is zero, $icode B$$ is not used. 00138 00139 $head X$$ 00140 The argument $icode X$$ has the prototype 00141 $codei% 00142 %FloatVector% &%X% 00143 %$$ 00144 and the size of $icode X$$ must equal $latex n * m$$ 00145 (see description of $cref/FloatVector/LuSolve/FloatVector/$$ below). 00146 The input value of $icode X$$ does not matter. 00147 On output, the elements of $icode X$$ contain the solution 00148 of the equation we wish to solve 00149 (unless $icode signdet$$ is equal to zero). 00150 If $icode m$$ is zero, $icode X$$ is not used. 00151 00152 $head logdet$$ 00153 The argument $icode logdet$$ has prototype 00154 $codei% 00155 %Float% &%logdet% 00156 %$$ 00157 On input, the value of $icode logdet$$ does not matter. 00158 On output, it has been set to the 00159 log of the determinant of $icode A$$ 00160 (but not quite). 00161 To be more specific, 00162 the determinant of $icode A$$ is given by the formula 00163 $codei% 00164 %det% = %signdet% * exp( %logdet% ) 00165 %$$ 00166 This enables $code LuSolve$$ to use logs of absolute values 00167 in the case where $icode Float$$ corresponds to a real number. 00168 00169 $head Float$$ 00170 The type $icode Float$$ must satisfy the conditions 00171 for a $cref NumericType$$ type. 00172 The routine $cref 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 $icode Float$$ objects $icode x$$ and $icode y$$: 00176 00177 $table 00178 $bold Operation$$ $cnext $bold Description$$ $rnext 00179 $codei%log(%x%)%$$ $cnext 00180 returns the logarithm of $icode x$$ as a $icode Float$$ object 00181 $tend 00182 00183 $head FloatVector$$ 00184 The type $icode FloatVector$$ must be a $cref SimpleVector$$ class with 00185 $cref/elements of type Float/SimpleVector/Elements of Specified Type/$$. 00186 The routine $cref 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 $codei% 00192 template <typename %Float%> 00193 bool LeqZero<%Float%>(const %Float% &%x%) 00194 %$$ 00195 in the $code CppAD$$ namespace. 00196 This function returns true if $icode 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 $icode Float$$ corresponds to real numbers). 00200 This template function definition assumes that the operator 00201 $code <=$$ is defined for $icode Float$$ objects. 00202 If this operator is not defined for your use of $icode 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 $codei% 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 $icode x$$ is zero and false otherwise. 00218 00219 $head AbsGeq$$ 00220 Including the file $code lu_solve.hpp$$ defines the template function 00221 $codei% 00222 template <typename %Float%> 00223 bool AbsGeq<%Float%>(const %Float% &%x%, const %Float% &%y%) 00224 %$$ 00225 If the type $icode 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 $cref/LuFactor/LuFactor/AbsGeq/$$. 00228 00229 $children% 00230 example/lu_solve.cpp% 00231 omh/lu_solve_hpp.omh 00232 %$$ 00233 $head Example$$ 00234 The file 00235 $cref lu_solve.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 C++ 00247 # include <complex> 00248 # include <vector> 00249 00250 // link exp for float and double cases 00251 # include <cppad/base_require.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 size_t(A.size()) == n * n, 00300 "Error in LuSolve: A must have size equal to n * n" 00301 ); 00302 CPPAD_ASSERT_KNOWN( 00303 size_t(B.size()) == n * m, 00304 "Error in LuSolve: B must have size equal to n * m" 00305 ); 00306 CPPAD_ASSERT_KNOWN( 00307 size_t(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 C++ 00351 # endif