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