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