CppAD: A C++ Algorithmic Differentiation Package 20110419
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_LU_FACTOR_INCLUDED 00003 # define CPPAD_LU_FACTOR_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-08 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 LuFactor$$ 00018 $escape #$$ 00019 $spell 00020 cppad.hpp 00021 Cpp 00022 Geq 00023 Lu 00024 bool 00025 const 00026 ip 00027 jp 00028 namespace 00029 std 00030 typename 00031 $$ 00032 00033 $index LuFactor$$ 00034 $index linear, Lu factor equation$$ 00035 $index equation, Lu factor$$ 00036 $index determinant, Lu factor$$ 00037 $index solve, Lu factor$$ 00038 00039 $section LU Factorization of A Square Matrix$$ 00040 00041 $pre 00042 $$ 00043 00044 $head Syntax$$ $code# include <cppad/lu_factor.hpp>$$ 00045 $pre 00046 $$ 00047 $syntax%%sign% = LuFactor(%ip%, %jp%, %LU%)%$$ 00048 00049 00050 $head Description$$ 00051 Computes an LU factorization of the matrix $italic A$$ 00052 where $italic A$$ is a square matrix. 00053 00054 $head Include$$ 00055 The file $code cppad/lu_factor.hpp$$ is included by $code cppad/cppad.hpp$$ 00056 but it can also be included separately with out the rest of 00057 the $code CppAD$$ routines. 00058 00059 $head Matrix Storage$$ 00060 All matrices are stored in row major order. 00061 To be specific, if $latex Y$$ is a vector 00062 that contains a $latex p$$ by $latex q$$ matrix, 00063 the size of $latex Y$$ must be equal to $latex p * q $$ and for 00064 $latex i = 0 , \ldots , p-1$$, 00065 $latex j = 0 , \ldots , q-1$$, 00066 $latex \[ 00067 Y_{i,j} = Y[ i * q + j ] 00068 \] $$ 00069 00070 $head sign$$ 00071 The return value $italic sign$$ has prototype 00072 $syntax% 00073 int %sign% 00074 %$$ 00075 If $italic A$$ is invertible, $italic sign$$ is plus or minus one 00076 and is the sign of the permutation corresponding to the row ordering 00077 $italic ip$$ and column ordering $italic jp$$. 00078 If $italic A$$ is not invertible, $italic sign$$ is zero. 00079 00080 $head ip$$ 00081 The argument $italic ip$$ has prototype 00082 $syntax% 00083 %SizeVector% &%ip% 00084 %$$ 00085 (see description of $xref/LuFactor/SizeVector/SizeVector/$$ below). 00086 The size of $italic ip$$ is referred to as $italic n$$ in the 00087 specifications below. 00088 The input value of the elements of $italic ip$$ does not matter. 00089 The output value of the elements of $italic ip$$ determine 00090 the order of the rows in the permuted matrix. 00091 00092 $head jp$$ 00093 The argument $italic jp$$ has prototype 00094 $syntax% 00095 %SizeVector% &%jp% 00096 %$$ 00097 (see description of $xref/LuFactor/SizeVector/SizeVector/$$ below). 00098 The size of $italic jp$$ must be equal to $italic n$$. 00099 The input value of the elements of $italic jp$$ does not matter. 00100 The output value of the elements of $italic jp$$ determine 00101 the order of the columns in the permuted matrix. 00102 00103 $head LU$$ 00104 The argument $italic LU$$ has the prototype 00105 $syntax% 00106 %FloatVector% &%LU% 00107 %$$ 00108 and the size of $italic LU$$ must equal $latex n * n$$ 00109 (see description of $xref/LuFactor/FloatVector/FloatVector/$$ below). 00110 00111 $subhead A$$ 00112 We define $italic A$$ as the matrix corresponding to the input 00113 value of $italic LU$$. 00114 00115 $subhead P$$ 00116 We define the permuted matrix $italic P$$ in terms of $italic A$$ by 00117 $syntax% 00118 %P%(%i%, %j%) = %A%[ %ip%[%i%] * %n% + %jp%[%j%] ] 00119 %$$ 00120 00121 $subhead L$$ 00122 We define the lower triangular matrix $italic L$$ in terms of the 00123 output value of $italic LU$$. 00124 The matrix $italic L$$ is zero above the diagonal 00125 and the rest of the elements are defined by 00126 $syntax% 00127 %L%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ] 00128 %$$ 00129 for $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , i$$. 00130 00131 $subhead U$$ 00132 We define the upper triangular matrix $italic U$$ in terms of the 00133 output value of $italic LU$$. 00134 The matrix $italic U$$ is zero below the diagonal, 00135 one on the diagonal, 00136 and the rest of the elements are defined by 00137 $syntax% 00138 %U%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ] 00139 %$$ 00140 for $latex i = 0 , \ldots , n-2$$ and $latex j = i+1 , \ldots , n-1$$. 00141 00142 $subhead Factor$$ 00143 If the return value $italic sign$$ is non-zero, 00144 $syntax% 00145 %L% * %U% = %P% 00146 %$$ 00147 If the return value of $italic sign$$ is zero, 00148 the contents of $italic L$$ and $italic U$$ are not defined. 00149 00150 $subhead Determinant$$ 00151 $index determinant$$ 00152 If the return value $italic sign$$ is zero, 00153 the determinant of $italic A$$ is zero. 00154 If $italic sign$$ is non-zero, 00155 using the output value of $italic LU$$ 00156 the determinant of the matrix $italic A$$ is equal to 00157 $syntax% 00158 %sign% * %LU%[%ip%[0], %jp%[0]] * %...% * %LU%[%ip%[%n%-1], %jp%[%n%-1]] 00159 %$$ 00160 00161 $head SizeVector$$ 00162 The type $italic SizeVector$$ must be a $xref/SimpleVector/$$ class with 00163 $xref/SimpleVector/Elements of Specified Type/elements of type size_t/$$. 00164 The routine $xref/CheckSimpleVector/$$ will generate an error message 00165 if this is not the case. 00166 00167 $head FloatVector$$ 00168 The type $italic FloatVector$$ must be a 00169 $xref/SimpleVector//simple vector class/$$. 00170 The routine $xref/CheckSimpleVector/$$ will generate an error message 00171 if this is not the case. 00172 00173 $head Float$$ 00174 This notation is used to denote the type corresponding 00175 to the elements of a $italic FloatVector$$. 00176 The type $italic Float$$ must satisfy the conditions 00177 for a $xref/NumericType/$$ type. 00178 The routine $xref/CheckNumericType/$$ will generate an error message 00179 if this is not the case. 00180 In addition, the following operations must be defined for any pair 00181 of $italic Float$$ objects $italic x$$ and $italic y$$: 00182 00183 $table 00184 $bold Operation$$ $cnext $bold Description$$ $rnext 00185 $syntax%log(%x%)%$$ $cnext 00186 returns the logarithm of $italic x$$ as a $italic Float$$ object 00187 $tend 00188 00189 $head AbsGeq$$ 00190 Including the file $code lu_factor.hpp$$ defines the template function 00191 $syntax% 00192 template <typename %Float%> 00193 bool AbsGeq<%Float%>(const %Float% &%x%, const %Float% &%y%) 00194 %$$ 00195 in the $code CppAD$$ namespace. 00196 This function returns true if the absolute value of 00197 $italic x$$ is greater than or equal the absolute value of $italic y$$. 00198 It is used by $code LuFactor$$ to choose the pivot elements. 00199 This template function definition uses the operator 00200 $code <=$$ to obtain the absolute value 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 LuFactor$$. 00204 $pre 00205 00206 $$ 00207 Complex numbers do not have the operation $code <=$$ defined. 00208 The specializations 00209 $syntax% 00210 bool AbsGeq< std::complex<float> > 00211 (const std::complex<float> &%x%, const std::complex<float> &%y%) 00212 bool AbsGeq< std::complex<double> > 00213 (const std::complex<double> &%x%, const std::complex<double> &%y%) 00214 %$$ 00215 are define by including $code lu_factor.hpp$$ 00216 These return true if the sum of the square of the real and imaginary parts 00217 of $italic x$$ is greater than or equal the 00218 sum of the square of the real and imaginary parts of $italic y$$. 00219 00220 $children% 00221 example/lu_factor.cpp% 00222 omh/lu_factor_hpp.omh 00223 %$$ 00224 $head Example$$ 00225 The file 00226 $xref/LuFactor.cpp/$$ 00227 contains an example and test of using $code LuFactor$$ by itself. 00228 It returns true if it succeeds and false otherwise. 00229 $pre 00230 00231 $$ 00232 The file $xref/lu_solve.hpp/$$ provides a useful example usage of 00233 $code LuFactor$$ with $code LuInvert$$. 00234 00235 $head Source$$ 00236 The file $cref/lu_factor.hpp/$$ contains the 00237 current source code that implements these specifications. 00238 00239 $end 00240 -------------------------------------------------------------------------- 00241 */ 00242 // BEGIN PROGRAM 00243 00244 # include <complex> 00245 # include <vector> 00246 00247 # include <cppad/local/cppad_assert.hpp> 00248 # include <cppad/check_simple_vector.hpp> 00249 # include <cppad/check_numeric_type.hpp> 00250 00251 namespace CppAD { // BEGIN CppAD namespace 00252 00253 // AbsGeq 00254 template <typename Float> 00255 inline bool AbsGeq(const Float &x, const Float &y) 00256 { Float xabs = x; 00257 if( xabs <= Float(0) ) 00258 xabs = - xabs; 00259 Float yabs = y; 00260 if( yabs <= Float(0) ) 00261 yabs = - yabs; 00262 return xabs >= yabs; 00263 } 00264 inline bool AbsGeq( 00265 const std::complex<double> &x, 00266 const std::complex<double> &y) 00267 { double xsq = x.real() * x.real() + x.imag() * x.imag(); 00268 double ysq = y.real() * y.real() + y.imag() * y.imag(); 00269 00270 return xsq >= ysq; 00271 } 00272 inline bool AbsGeq( 00273 const std::complex<float> &x, 00274 const std::complex<float> &y) 00275 { float xsq = x.real() * x.real() + x.imag() * x.imag(); 00276 float ysq = y.real() * y.real() + y.imag() * y.imag(); 00277 00278 return xsq >= ysq; 00279 } 00280 00281 // Lines that are different from code in cppad/local/lu_ratio.hpp end with // 00282 template <class SizeVector, class FloatVector> // 00283 int LuFactor(SizeVector &ip, SizeVector &jp, FloatVector &LU) // 00284 { 00285 // type of the elements of LU // 00286 typedef typename FloatVector::value_type Float; // 00287 00288 // check numeric type specifications 00289 CheckNumericType<Float>(); 00290 00291 // check simple vector class specifications 00292 CheckSimpleVector<Float, FloatVector>(); 00293 CheckSimpleVector<size_t, SizeVector>(); 00294 00295 size_t i, j; // some temporary indices 00296 const Float zero( 0 ); // the value zero as a Float object 00297 size_t imax; // row index of maximum element 00298 size_t jmax; // column indx of maximum element 00299 Float emax; // maximum absolute value 00300 size_t p; // count pivots 00301 int sign; // sign of the permutation 00302 Float etmp; // temporary element 00303 Float pivot; // pivot element 00304 00305 // ------------------------------------------------------- 00306 size_t n = ip.size(); 00307 CPPAD_ASSERT_KNOWN( 00308 jp.size() == n, 00309 "Error in LuFactor: jp must have size equal to n" 00310 ); 00311 CPPAD_ASSERT_KNOWN( 00312 LU.size() == n * n, 00313 "Error in LuFactor: LU must have size equal to n * m" 00314 ); 00315 // ------------------------------------------------------- 00316 00317 // initialize row and column order in matrix not yet pivoted 00318 for(i = 0; i < n; i++) 00319 { ip[i] = i; 00320 jp[i] = i; 00321 } 00322 // initialize the sign of the permutation 00323 sign = 1; 00324 // --------------------------------------------------------- 00325 00326 // Reduce the matrix P to L * U using n pivots 00327 for(p = 0; p < n; p++) 00328 { // determine row and column corresponding to element of 00329 // maximum absolute value in remaining part of P 00330 imax = jmax = n; 00331 emax = zero; 00332 for(i = p; i < n; i++) 00333 { for(j = p; j < n; j++) 00334 { CPPAD_ASSERT_UNKNOWN( 00335 (ip[i] < n) & (jp[j] < n) 00336 ); 00337 etmp = LU[ ip[i] * n + jp[j] ]; 00338 00339 // check if maximum absolute value so far 00340 if( AbsGeq (etmp, emax) ) 00341 { imax = i; 00342 jmax = j; 00343 emax = etmp; 00344 } 00345 } 00346 } 00347 CPPAD_ASSERT_KNOWN( 00348 (imax < n) & (jmax < n) , 00349 "LuFactor can't determine an element with " 00350 "maximum absolute value.\n" 00351 "Perhaps original matrix contains not a number or infinity.\n" 00352 "Perhaps your specialization of AbsGeq is not correct." 00353 ); 00354 if( imax != p ) 00355 { // switch rows so max absolute element is in row p 00356 i = ip[p]; 00357 ip[p] = ip[imax]; 00358 ip[imax] = i; 00359 sign = -sign; 00360 } 00361 if( jmax != p ) 00362 { // switch columns so max absolute element is in column p 00363 j = jp[p]; 00364 jp[p] = jp[jmax]; 00365 jp[jmax] = j; 00366 sign = -sign; 00367 } 00368 // pivot using the max absolute element 00369 pivot = LU[ ip[p] * n + jp[p] ]; 00370 00371 // check for determinant equal to zero 00372 if( pivot == zero ) 00373 { // abort the mission 00374 return 0; 00375 } 00376 00377 // Reduce U by the elementary transformations that maps 00378 // LU( ip[p], jp[p] ) to one. Only need transform elements 00379 // above the diagonal in U and LU( ip[p] , jp[p] ) is 00380 // corresponding value below diagonal in L. 00381 for(j = p+1; j < n; j++) 00382 LU[ ip[p] * n + jp[j] ] /= pivot; 00383 00384 // Reduce U by the elementary transformations that maps 00385 // LU( ip[i], jp[p] ) to zero. Only need transform elements 00386 // above the diagonal in U and LU( ip[i], jp[p] ) is 00387 // corresponding value below diagonal in L. 00388 for(i = p+1; i < n; i++ ) 00389 { etmp = LU[ ip[i] * n + jp[p] ]; 00390 for(j = p+1; j < n; j++) 00391 { LU[ ip[i] * n + jp[j] ] -= 00392 etmp * LU[ ip[p] * n + jp[j] ]; 00393 } 00394 } 00395 } 00396 return sign; 00397 } 00398 } // END CppAD namespace 00399 // END PROGRAM 00400 # endif