CppAD: A C++ Algorithmic Differentiation Package  20171217
lu_ratio.hpp
Go to the documentation of this file.
3
4 /* --------------------------------------------------------------------------
6
8 the terms of the
9  Eclipse Public License Version 1.0.
10
11 A copy of this license is included in the COPYING file of this distribution.
13 -------------------------------------------------------------------------- */
14
15 /*
16 $begin LuRatio$$17 spell 18 cppad.hpp 19 xk 20 Cpp 21 Lu 22 bool 23 const 24 ip 25 jp 26 std 27 ADvector 28$$ 29 30 31$section LU Factorization of A Square Matrix and Stability Calculation$$32 mindex LuRatio linear equation solve$$
33
34 $head Syntax$$35 code# include <cppad/cppad.hpp>$$ 36$pre
37 $$38 icode%sign% = LuRatio(%ip%, %jp%, %LU%, %ratio%)%$$
39
40
41 $head Description$$42 Computes an LU factorization of the matrix icode A$$ 43 where$icode A$$is a square matrix. 44 A measure of the numerical stability called icode ratio$$ is calculated.
45 This ratio is useful when the results of $code LuRatio$$are 46 used as part of an cref ADFun$$ object. 47 48$head Include$$49 This routine is designed to be used with AD objects and 50 requires the code cppad/cppad.hpp$$ file to be included.
51
52 $head Matrix Storage$$53 All matrices are stored in row major order. 54 To be specific, if latex Y$$ is a vector 55 that contains a$latex p$$by latex q$$ matrix,
56 the size of $latex Y$$must be equal to latex p * q$$ and for 57$latex i = 0 , \ldots , p-1$$, 58 latex j = 0 , \ldots , q-1$$,
59 $latex $60 Y_{i,j} = Y[ i * q + j ] 61$ $$62 63 head sign$$ 64 The return value$icode sign$$has prototype 65 codei% 66 int %sign% 67 %$$
68 If $icode A$$is invertible, icode sign$$ is plus or minus one 69 and is the sign of the permutation corresponding to the row ordering 70$icode ip$$and column ordering icode jp$$.
71 If $icode A$$is not invertible, icode sign$$ is zero. 72 73$head ip$$74 The argument icode ip$$ has prototype
75 $codei% 76 %SizeVector% &%ip% 77 %$$78 (see description of cref/SizeVector/LuFactor/SizeVector/$$ below). 79 The size of$icode ip$$is referred to as icode n$$ in the
80 specifications below.
81 The input value of the elements of $icode ip$$does not matter. 82 The output value of the elements of icode ip$$ determine 83 the order of the rows in the permuted matrix. 84 85$head jp$$86 The argument icode jp$$ has prototype
87 $codei% 88 %SizeVector% &%jp% 89 %$$90 (see description of cref/SizeVector/LuFactor/SizeVector/$$ below). 91 The size of$icode jp$$must be equal to icode n$$.
92 The input value of the elements of $icode jp$$does not matter. 93 The output value of the elements of icode jp$$ determine 94 the order of the columns in the permuted matrix. 95 96$head LU$$97 The argument icode LU$$ has the prototype
98 $codei% 99 %ADvector% &%LU% 100 %$$101 and the size of icode LU$$ must equal$latex n * n$$102 (see description of cref/ADvector/LuRatio/ADvector/$$ below).
103
104 $subhead A$$105 We define icode A$$ as the matrix corresponding to the input 106 value of$icode LU$$. 107 108 subhead P$$
109 We define the permuted matrix $icode P$$in terms of icode A$$ by 110$codei%
111  %P%(%i%, %j%) = %A%[ %ip%[%i%] * %n% + %jp%[%j%] ]
112 %$$113 114 subhead L$$
115 We define the lower triangular matrix $icode L$$in terms of the 116 output value of icode LU$$. 117 The matrix$icode L$$is zero above the diagonal 118 and the rest of the elements are defined by 119 codei% 120 %L%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ] 121 %$$
122 for $latex i = 0 , \ldots , n-1$$and latex j = 0 , \ldots , i$$. 123 124$subhead U$$125 We define the upper triangular matrix icode U$$ in terms of the
126 output value of $icode LU$$. 127 The matrix icode U$$ is zero below the diagonal, 128 one on the diagonal, 129 and the rest of the elements are defined by 130$codei%
131  %U%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ]
132 %$$133 for latex i = 0 , \ldots , n-2$$ and $latex j = i+1 , \ldots , n-1$$. 134 135 subhead Factor$$ 136 If the return value$icode sign$$is non-zero, 137 codei% 138 %L% * %U% = %P% 139 %$$
140 If the return value of $icode sign$$is zero, 141 the contents of icode L$$ and$icode U$$are not defined. 142 143 subhead Determinant$$
144 If the return value $icode sign$$is zero, 145 the determinant of icode A$$ is zero. 146 If$icode sign$$is non-zero, 147 using the output value of icode LU$$
148 the determinant of the matrix $icode A$$is equal to 149 codei% 150 %sign% * %LU%[%ip%[0], %jp%[0]] * %...% * %LU%[%ip%[%n%-1], %jp%[%n%-1]] 151 %$$ 152 153$head ratio$$154 The argument icode ratio$$ has prototype
155 $codei% 156 AD<%Base%> &%ratio% 157 %$$158 On input, the value of icode ratio$$ does not matter. 159 On output it is a measure of how good the choice of pivots is. 160 For$latex p = 0 , \ldots , n-1$$, 161 the th p$$ pivot element is the element of maximum absolute value of a
162 $latex (n-p) \times (n-p)$$sub-matrix. 163 The ratio of each element of sub-matrix divided by the pivot element 164 is computed. 165 The return value of icode ratio$$ is the maximum absolute value of 166 such ratios over with respect to all elements and all the pivots. 167 168$subhead Purpose$$169 Suppose that the execution of a call to code LuRatio$$
170 is recorded in the $codei%ADFun<%Base%>%$$object icode F$$. 171 Then a call to$cref Forward$$of the form 172 codei% 173 %F%.Forward(%k%, %xk%) 174 %$$
175 with $icode k$$equal to zero will revaluate this Lu factorization 176 with the same pivots and a new value for icode A$$. 177 In this case, the resulting$icode ratio$$may not be one. 178 If icode ratio$$ is too large (the meaning of too large is up to you),
179 the current pivots do not yield a stable LU factorization of $icode A$$. 180 A better choice for the pivots (for this value of icode A$$) 181 will be made if you recreate the$code ADFun$$object 182 starting with the cref Independent$$ variable values
183 that correspond to the vector $icode xk$$. 184 185 head SizeVector$$ 186 The type$icode SizeVector$$must be a cref SimpleVector$$ class with
187 $cref/elements of type size_t/SimpleVector/Elements of Specified Type/$$. 188 The routine cref CheckSimpleVector$$ will generate an error message 189 if this is not the case. 190 191$head ADvector$$192 The type icode ADvector$$ must be a
193 $cref/simple vector class/SimpleVector/$$with elements of type 194 codei%AD<%Base%>%$$. 195 The routine$cref CheckSimpleVector$$will generate an error message 196 if this is not the case. 197 198 199 head Example$$
200 $children% 201 example/general/lu_ratio.cpp 202 %$$203 The file cref lu_ratio.cpp$$ 204 contains an example and test of using$code LuRatio.
205 It returns true if it succeeds and false otherwise.
206
207 \$end
208 --------------------------------------------------------------------------
209 */
211
212 // Lines different from the code in cppad/lu_factor.hpp end with //
213 template <class SizeVector, class ADvector, class Base> //
215 {
218
219  // check numeric type specifications
220  CheckNumericType<Float>();
221
222  // check simple vector class specifications
223  CheckSimpleVector<Float, FloatVector>();
224  CheckSimpleVector<size_t, SizeVector>();
225
226  size_t i, j; // some temporary indices
227  const Float zero( 0 ); // the value zero as a Float object
228  size_t imax; // row index of maximum element
229  size_t jmax; // column indx of maximum element
230  Float emax; // maximum absolute value
231  size_t p; // count pivots
232  int sign; // sign of the permutation
233  Float etmp; // temporary element
234  Float pivot; // pivot element
235
236  // -------------------------------------------------------
237  size_t n = size_t(ip.size());
239  size_t(jp.size()) == n,
240  "Error in LuFactor: jp must have size equal to n"
241  );
243  size_t(LU.size()) == n * n,
244  "Error in LuFactor: LU must have size equal to n * m"
245  );
246  // -------------------------------------------------------
247
248  // initialize row and column order in matrix not yet pivoted
249  for(i = 0; i < n; i++)
250  { ip[i] = i;
251  jp[i] = i;
252  }
253  // initialize the sign of the permutation
254  sign = 1;
255  // initialize the ratio //
256  ratio = Float(1); //
257  // ---------------------------------------------------------
258
259  // Reduce the matrix P to L * U using n pivots
260  for(p = 0; p < n; p++)
261  { // determine row and column corresponding to element of
262  // maximum absolute value in remaining part of P
263  imax = jmax = n;
264  emax = zero;
265  for(i = p; i < n; i++)
266  { for(j = p; j < n; j++)
268  (ip[i] < n) & (jp[j] < n)
269  );
270  etmp = LU[ ip[i] * n + jp[j] ];
271
272  // check if maximum absolute value so far
273  if( AbsGeq (etmp, emax) )
274  { imax = i;
275  jmax = j;
276  emax = etmp;
277  }
278  }
279  }
280  for(i = p; i < n; i++) //
281  { for(j = p; j < n; j++) //
282  { etmp = fabs(LU[ ip[i] * n + jp[j] ] / emax); //
283  ratio = //
284  CondExpGt(etmp, ratio, etmp, ratio); //
285  } //
286  } //
288  (imax < n) & (jmax < n) ,
289  "AbsGeq must return true when second argument is zero"
290  );
291  if( imax != p )
292  { // switch rows so max absolute element is in row p
293  i = ip[p];
294  ip[p] = ip[imax];
295  ip[imax] = i;
296  sign = -sign;
297  }
298  if( jmax != p )
299  { // switch columns so max absolute element is in column p
300  j = jp[p];
301  jp[p] = jp[jmax];
302  jp[jmax] = j;
303  sign = -sign;
304  }
305  // pivot using the max absolute element
306  pivot = LU[ ip[p] * n + jp[p] ];
307
308  // check for determinant equal to zero
309  if( pivot == zero )
310  { // abort the mission
311  return 0;
312  }
313
314  // Reduce U by the elementary transformations that maps
315  // LU( ip[p], jp[p] ) to one. Only need transform elements
316  // above the diagonal in U and LU( ip[p] , jp[p] ) is
317  // corresponding value below diagonal in L.
318  for(j = p+1; j < n; j++)
319  LU[ ip[p] * n + jp[j] ] /= pivot;
320
321  // Reduce U by the elementary transformations that maps
322  // LU( ip[i], jp[p] ) to zero. Only need transform elements
323  // above the diagonal in U and LU( ip[i], jp[p] ) is
324  // corresponding value below diagonal in L.
325  for(i = p+1; i < n; i++ )
326  { etmp = LU[ ip[i] * n + jp[p] ];
327  for(j = p+1; j < n; j++)
328  { LU[ ip[i] * n + jp[j] ] -=
329  etmp * LU[ ip[p] * n + jp[j] ];
330  }
331  }
332  }
333  return sign;
334 }
335 } // END CppAD namespace
336
337 # endif