00001 # ifndef CPPAD_LU_RATIO_INCLUDED 00002 # define CPPAD_LU_RATIO_INCLUDED 00003 00004 /* -------------------------------------------------------------------------- 00005 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 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 LuRatio$$ 00017 $spell 00018 cppad.hpp 00019 xk 00020 Cpp 00021 Lu 00022 bool 00023 const 00024 ip 00025 jp 00026 std 00027 ADvector 00028 $$ 00029 00030 $index LuRatio$$ 00031 $index linear, Lu factor equation$$ 00032 $index equation, Lu factor$$ 00033 $index determinant, Lu factor$$ 00034 $index solve, Lu factor$$ 00035 00036 $section LU Factorization of A Square Matrix and Stability Calculation$$ 00037 00038 $head Syntax$$ 00039 $code# include <cppad/cppad.hpp>$$ 00040 $pre 00041 $$ 00042 $syntax%%sign% = LuRatio(%ip%, %jp%, %LU%, %ratio%)%$$ 00043 00044 00045 $head Description$$ 00046 Computes an LU factorization of the matrix $italic A$$ 00047 where $italic A$$ is a square matrix. 00048 A measure of the numerical stability called $italic ratio$$ is calculated. 00049 This ratio is useful when the results of $code LuRatio$$ are 00050 used as part of an $xref/ADFun/$$ object. 00051 00052 $head Include$$ 00053 This routine is designed to be used with AD objects and 00054 requires the $code cppad/cppad.hpp$$ file to be included. 00055 00056 $head Matrix Storage$$ 00057 All matrices are stored in row major order. 00058 To be specific, if $latex Y$$ is a vector 00059 that contains a $latex p$$ by $latex q$$ matrix, 00060 the size of $latex Y$$ must be equal to $latex p * q $$ and for 00061 $latex i = 0 , \ldots , p-1$$, 00062 $latex j = 0 , \ldots , q-1$$, 00063 $latex \[ 00064 Y_{i,j} = Y[ i * q + j ] 00065 \] $$ 00066 00067 $head sign$$ 00068 The return value $italic sign$$ has prototype 00069 $syntax% 00070 int %sign% 00071 %$$ 00072 If $italic A$$ is invertible, $italic sign$$ is plus or minus one 00073 and is the sign of the permutation corresponding to the row ordering 00074 $italic ip$$ and column ordering $italic jp$$. 00075 If $italic A$$ is not invertible, $italic sign$$ is zero. 00076 00077 $head ip$$ 00078 The argument $italic ip$$ has prototype 00079 $syntax% 00080 %SizeVector% &%ip% 00081 %$$ 00082 (see description of $xref/LuFactor/SizeVector/SizeVector/$$ below). 00083 The size of $italic ip$$ is referred to as $italic n$$ in the 00084 specifications below. 00085 The input value of the elements of $italic ip$$ does not matter. 00086 The output value of the elements of $italic ip$$ determine 00087 the order of the rows in the permuted matrix. 00088 00089 $head jp$$ 00090 The argument $italic jp$$ has prototype 00091 $syntax% 00092 %SizeVector% &%jp% 00093 %$$ 00094 (see description of $xref/LuFactor/SizeVector/SizeVector/$$ below). 00095 The size of $italic jp$$ must be equal to $italic n$$. 00096 The input value of the elements of $italic jp$$ does not matter. 00097 The output value of the elements of $italic jp$$ determine 00098 the order of the columns in the permuted matrix. 00099 00100 $head LU$$ 00101 The argument $italic LU$$ has the prototype 00102 $syntax% 00103 %ADvector% &%LU% 00104 %$$ 00105 and the size of $italic LU$$ must equal $latex n * n$$ 00106 (see description of $xref/LuRatio/ADvector/ADvector/$$ below). 00107 00108 $subhead A$$ 00109 We define $italic A$$ as the matrix corresponding to the input 00110 value of $italic LU$$. 00111 00112 $subhead P$$ 00113 We define the permuted matrix $italic P$$ in terms of $italic A$$ by 00114 $syntax% 00115 %P%(%i%, %j%) = %A%[ %ip%[%i%] * %n% + %jp%[%j%] ] 00116 %$$ 00117 00118 $subhead L$$ 00119 We define the lower triangular matrix $italic L$$ in terms of the 00120 output value of $italic LU$$. 00121 The matrix $italic L$$ is zero above the diagonal 00122 and the rest of the elements are defined by 00123 $syntax% 00124 %L%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ] 00125 %$$ 00126 for $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , i$$. 00127 00128 $subhead U$$ 00129 We define the upper triangular matrix $italic U$$ in terms of the 00130 output value of $italic LU$$. 00131 The matrix $italic U$$ is zero below the diagonal, 00132 one on the diagonal, 00133 and the rest of the elements are defined by 00134 $syntax% 00135 %U%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ] 00136 %$$ 00137 for $latex i = 0 , \ldots , n-2$$ and $latex j = i+1 , \ldots , n-1$$. 00138 00139 $subhead Factor$$ 00140 If the return value $italic sign$$ is non-zero, 00141 $syntax% 00142 %L% * %U% = %P% 00143 %$$ 00144 If the return value of $italic sign$$ is zero, 00145 the contents of $italic L$$ and $italic U$$ are not defined. 00146 00147 $subhead Determinant$$ 00148 $index determinant$$ 00149 If the return value $italic sign$$ is zero, 00150 the determinant of $italic A$$ is zero. 00151 If $italic sign$$ is non-zero, 00152 using the output value of $italic LU$$ 00153 the determinant of the matrix $italic A$$ is equal to 00154 $syntax% 00155 %sign% * %LU%[%ip%[0], %jp%[0]] * %...% * %LU%[%ip%[%n%-1], %jp%[%n%-1]] 00156 %$$ 00157 00158 $head ratio$$ 00159 The argument $italic ratio$$ has prototype 00160 $syntax% 00161 AD<%Base%> &%ratio% 00162 %$$ 00163 On input, the value of $italic ratio$$ does not matter. 00164 On output it is a measure of how good the choice of pivots is. 00165 For $latex p = 0 , \ldots , n-1$$, 00166 the $th p$$ pivot element is the element of maximum absolute value of a 00167 $latex (n-p) \times (n-p)$$ sub-matrix. 00168 The ratio of each element of sub-matrix divided by the pivot element 00169 is computed. 00170 The return value of $italic ratio$$ is the maximum absolute value of 00171 such ratios over with respect to all elements and all the pivots. 00172 00173 $subhead Purpose$$ 00174 Suppose that the execution of a call to $code LuRatio$$ 00175 is recorded in the $syntax%ADFun<%Base%>%$$ object $italic F$$. 00176 Then a call to $xref/Forward/$$ of the form 00177 $syntax% 00178 %F%.Forward(%k%, %xk%) 00179 %$$ 00180 with $italic k$$ equal to zero will revaluate this Lu factorization 00181 with the same pivots and a new value for $italic A$$. 00182 In this case, the resulting $italic ratio$$ may not be one. 00183 If $italic ratio$$ is too large (the meaning of too large is up to you), 00184 the current pivots do not yield a stable LU factorization of $italic A$$. 00185 A better choice for the pivots (for this value of $italic A$$) 00186 will be made if you recreate the $code ADFun$$ object 00187 starting with the $xref/Independent/$$ variable values 00188 that correspond to the vector $italic xk$$. 00189 00190 $head SizeVector$$ 00191 The type $italic SizeVector$$ must be a $xref/SimpleVector/$$ class with 00192 $xref/SimpleVector/Elements of Specified Type/elements of type size_t/$$. 00193 The routine $xref/CheckSimpleVector/$$ will generate an error message 00194 if this is not the case. 00195 00196 $head Base$$ 00197 This type is determined by the argument $italic ratio$$ 00198 which has type $syntax%AD<%Base%>%$$. 00199 The type $italic Base$$ must be 00200 $code float$$, $code double$$, or any type in the 00201 $cref/AD levels above/glossary/AD Levels Above Base/$$ 00202 above $code float$$ or $code double$$. 00203 00204 $head ADvector$$ 00205 The type $italic ADvector$$ must be a 00206 $xref/SimpleVector//simple vector class/$$ with elements of type 00207 $syntax%AD<%Base%>%$$. 00208 The routine $xref/CheckSimpleVector/$$ will generate an error message 00209 if this is not the case. 00210 00211 00212 $head Example$$ 00213 $children% 00214 example/lu_ratio.cpp 00215 %$$ 00216 The file $xref/LuRatio.cpp/$$ 00217 contains an example and test of using $code LuRatio$$. 00218 It returns true if it succeeds and false otherwise. 00219 00220 $end 00221 -------------------------------------------------------------------------- 00222 */ 00223 namespace CppAD { // BEGIN CppAD namespace 00224 00225 // Lines different from the code in cppad/lu_factor.hpp end with // 00226 template <class SizeVector, class ADvector, class Base> // 00227 int LuRatio(SizeVector &ip, SizeVector &jp, ADvector &LU, AD<Base> &ratio) // 00228 { 00229 typedef ADvector FloatVector; // 00230 typedef AD<Base> Float; // 00231 00232 // check numeric type specifications 00233 CheckNumericType<Float>(); 00234 00235 // check simple vector class specifications 00236 CheckSimpleVector<Float, FloatVector>(); 00237 CheckSimpleVector<size_t, SizeVector>(); 00238 00239 size_t i, j; // some temporary indices 00240 const Float zero( 0 ); // the value zero as a Float object 00241 size_t imax; // row index of maximum element 00242 size_t jmax; // column indx of maximum element 00243 Float emax; // maximum absolute value 00244 size_t p; // count pivots 00245 int sign; // sign of the permutation 00246 Float etmp; // temporary element 00247 Float pivot; // pivot element 00248 00249 // ------------------------------------------------------- 00250 size_t n = ip.size(); 00251 CPPAD_ASSERT_KNOWN( 00252 jp.size() == n, 00253 "Error in LuFactor: jp must have size equal to n" 00254 ); 00255 CPPAD_ASSERT_KNOWN( 00256 LU.size() == n * n, 00257 "Error in LuFactor: LU must have size equal to n * m" 00258 ); 00259 // ------------------------------------------------------- 00260 00261 // initialize row and column order in matrix not yet pivoted 00262 for(i = 0; i < n; i++) 00263 { ip[i] = i; 00264 jp[i] = i; 00265 } 00266 // initialize the sign of the permutation 00267 sign = 1; 00268 // initialize the ratio // 00269 ratio = Float(1); // 00270 // --------------------------------------------------------- 00271 00272 // Reduce the matrix P to L * U using n pivots 00273 for(p = 0; p < n; p++) 00274 { // determine row and column corresponding to element of 00275 // maximum absolute value in remaining part of P 00276 imax = jmax = n; 00277 emax = zero; 00278 for(i = p; i < n; i++) 00279 { for(j = p; j < n; j++) 00280 { CPPAD_ASSERT_UNKNOWN( 00281 (ip[i] < n) & (jp[j] < n) 00282 ); 00283 etmp = LU[ ip[i] * n + jp[j] ]; 00284 00285 // check if maximum absolute value so far 00286 if( AbsGeq (etmp, emax) ) 00287 { imax = i; 00288 jmax = j; 00289 emax = etmp; 00290 } 00291 } 00292 } 00293 for(i = p; i < n; i++) // 00294 { for(j = p; j < n; j++) // 00295 { etmp = abs(LU[ ip[i] * n + jp[j] ] / emax); // 00296 ratio = // 00297 CondExpGt(etmp, ratio, etmp, ratio); // 00298 } // 00299 } // 00300 CPPAD_ASSERT_KNOWN( 00301 (imax < n) & (jmax < n) , 00302 "AbsGeq must return true when second argument is zero" 00303 ); 00304 if( imax != p ) 00305 { // switch rows so max absolute element is in row p 00306 i = ip[p]; 00307 ip[p] = ip[imax]; 00308 ip[imax] = i; 00309 sign = -sign; 00310 } 00311 if( jmax != p ) 00312 { // switch columns so max absolute element is in column p 00313 j = jp[p]; 00314 jp[p] = jp[jmax]; 00315 jp[jmax] = j; 00316 sign = -sign; 00317 } 00318 // pivot using the max absolute element 00319 pivot = LU[ ip[p] * n + jp[p] ]; 00320 00321 // check for determinant equal to zero 00322 if( pivot == zero ) 00323 { // abort the mission 00324 return 0; 00325 } 00326 00327 // Reduce U by the elementary transformations that maps 00328 // LU( ip[p], jp[p] ) to one. Only need transform elements 00329 // above the diagonal in U and LU( ip[p] , jp[p] ) is 00330 // corresponding value below diagonal in L. 00331 for(j = p+1; j < n; j++) 00332 LU[ ip[p] * n + jp[j] ] /= pivot; 00333 00334 // Reduce U by the elementary transformations that maps 00335 // LU( ip[i], jp[p] ) to zero. Only need transform elements 00336 // above the diagonal in U and LU( ip[i], jp[p] ) is 00337 // corresponding value below diagonal in L. 00338 for(i = p+1; i < n; i++ ) 00339 { etmp = LU[ ip[i] * n + jp[p] ]; 00340 for(j = p+1; j < n; j++) 00341 { LU[ ip[i] * n + jp[j] ] -= 00342 etmp * LU[ ip[p] * n + jp[j] ]; 00343 } 00344 } 00345 } 00346 return sign; 00347 } 00348 } // END CppAD namespace 00349 00350 # endif
1.4.7