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