CppAD: A C++ Algorithmic Differentiation Package  20171217
lu_factor.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 LuFactor\$\$
17 \$escape #\$\$
18 \$spell
20  Cpp
21  Geq
22  Lu
23  bool
24  const
25  ip
26  jp
27  namespace
28  std
29  typename
30 \$\$
31
32
33 \$section LU Factorization of A Square Matrix\$\$
34 \$mindex LuFactor linear equation solve\$\$
35
36 \$pre
37 \$\$
38
40 %\$\$
41 \$icode%sign% = LuFactor(%ip%, %jp%, %LU%)%\$\$
42
43
45 Computes an LU factorization of the matrix \$icode A\$\$
46 where \$icode A\$\$ is a square matrix.
47
50 but it can also be included separately with out the rest of
52
54 All matrices are stored in row major order.
55 To be specific, if \$latex Y\$\$ is a vector
56 that contains a \$latex p\$\$ by \$latex q\$\$ matrix,
57 the size of \$latex Y\$\$ must be equal to \$latex p * q \$\$ and for
58 \$latex i = 0 , \ldots , p-1\$\$,
59 \$latex j = 0 , \ldots , q-1\$\$,
60 \$latex \[
61  Y_{i,j} = Y[ i * q + j ]
62 \] \$\$
63
65 The return value \$icode sign\$\$ has prototype
66 \$codei%
67  int %sign%
68 %\$\$
69 If \$icode A\$\$ is invertible, \$icode sign\$\$ is plus or minus one
70 and is the sign of the permutation corresponding to the row ordering
71 \$icode ip\$\$ and column ordering \$icode jp\$\$.
72 If \$icode A\$\$ is not invertible, \$icode sign\$\$ is zero.
73
75 The argument \$icode ip\$\$ has prototype
76 \$codei%
77  %SizeVector% &%ip%
78 %\$\$
79 (see description of \$cref/SizeVector/LuFactor/SizeVector/\$\$ below).
80 The size of \$icode ip\$\$ is referred to as \$icode n\$\$ in the
81 specifications below.
82 The input value of the elements of \$icode ip\$\$ does not matter.
83 The output value of the elements of \$icode ip\$\$ determine
84 the order of the rows in the permuted matrix.
85
87 The argument \$icode jp\$\$ has prototype
88 \$codei%
89  %SizeVector% &%jp%
90 %\$\$
91 (see description of \$cref/SizeVector/LuFactor/SizeVector/\$\$ below).
92 The size of \$icode jp\$\$ must be equal to \$icode n\$\$.
93 The input value of the elements of \$icode jp\$\$ does not matter.
94 The output value of the elements of \$icode jp\$\$ determine
95 the order of the columns in the permuted matrix.
96
98 The argument \$icode LU\$\$ has the prototype
99 \$codei%
100  %FloatVector% &%LU%
101 %\$\$
102 and the size of \$icode LU\$\$ must equal \$latex n * n\$\$
103 (see description of \$cref/FloatVector/LuFactor/FloatVector/\$\$ below).
104
106 We define \$icode A\$\$ as the matrix corresponding to the input
107 value of \$icode LU\$\$.
108
110 We define the permuted matrix \$icode P\$\$ in terms of \$icode A\$\$ by
111 \$codei%
112  %P%(%i%, %j%) = %A%[ %ip%[%i%] * %n% + %jp%[%j%] ]
113 %\$\$
114
116 We define the lower triangular matrix \$icode L\$\$ in terms of the
117 output value of \$icode LU\$\$.
118 The matrix \$icode L\$\$ is zero above the diagonal
119 and the rest of the elements are defined by
120 \$codei%
121  %L%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ]
122 %\$\$
123 for \$latex i = 0 , \ldots , n-1\$\$ and \$latex j = 0 , \ldots , i\$\$.
124
126 We define the upper triangular matrix \$icode U\$\$ in terms of the
127 output value of \$icode LU\$\$.
128 The matrix \$icode U\$\$ is zero below the diagonal,
129 one on the diagonal,
130 and the rest of the elements are defined by
131 \$codei%
132  %U%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ]
133 %\$\$
134 for \$latex i = 0 , \ldots , n-2\$\$ and \$latex j = i+1 , \ldots , n-1\$\$.
135
137 If the return value \$icode sign\$\$ is non-zero,
138 \$codei%
139  %L% * %U% = %P%
140 %\$\$
141 If the return value of \$icode sign\$\$ is zero,
142 the contents of \$icode L\$\$ and \$icode U\$\$ are not defined.
143
145 If the return value \$icode sign\$\$ is zero,
146 the determinant of \$icode A\$\$ is zero.
147 If \$icode sign\$\$ is non-zero,
148 using the output value of \$icode LU\$\$
149 the determinant of the matrix \$icode A\$\$ is equal to
150 \$codei%
151 %sign% * %LU%[%ip%[0], %jp%[0]] * %...% * %LU%[%ip%[%n%-1], %jp%[%n%-1]]
152 %\$\$
153
155 The type \$icode SizeVector\$\$ must be a \$cref SimpleVector\$\$ class with
156 \$cref/elements of type size_t/SimpleVector/Elements of Specified Type/\$\$.
157 The routine \$cref CheckSimpleVector\$\$ will generate an error message
158 if this is not the case.
159
161 The type \$icode FloatVector\$\$ must be a
162 \$cref/simple vector class/SimpleVector/\$\$.
163 The routine \$cref CheckSimpleVector\$\$ will generate an error message
164 if this is not the case.
165
167 This notation is used to denote the type corresponding
168 to the elements of a \$icode FloatVector\$\$.
169 The type \$icode Float\$\$ must satisfy the conditions
170 for a \$cref NumericType\$\$ type.
171 The routine \$cref CheckNumericType\$\$ will generate an error message
172 if this is not the case.
173 In addition, the following operations must be defined for any pair
174 of \$icode Float\$\$ objects \$icode x\$\$ and \$icode y\$\$:
175
176 \$table
177 \$bold Operation\$\$ \$cnext \$bold Description\$\$ \$rnext
178 \$codei%log(%x%)%\$\$ \$cnext
179  returns the logarithm of \$icode x\$\$ as a \$icode Float\$\$ object
180 \$tend
181
183 Including the file \$code lu_factor.hpp\$\$ defines the template function
184 \$codei%
185  template <typename %Float%>
186  bool AbsGeq<%Float%>(const %Float% &%x%, const %Float% &%y%)
187 %\$\$
188 in the \$code CppAD\$\$ namespace.
189 This function returns true if the absolute value of
190 \$icode x\$\$ is greater than or equal the absolute value of \$icode y\$\$.
191 It is used by \$code LuFactor\$\$ to choose the pivot elements.
192 This template function definition uses the operator
193 \$code <=\$\$ to obtain the absolute value for \$icode Float\$\$ objects.
194 If this operator is not defined for your use of \$icode Float\$\$,
195 you will need to specialize this template so that it works for your
196 use of \$code LuFactor\$\$.
197 \$pre
198
199 \$\$
200 Complex numbers do not have the operation \$code <=\$\$ defined.
201 The specializations
202 \$codei%
203 bool AbsGeq< std::complex<float> >
204  (const std::complex<float> &%x%, const std::complex<float> &%y%)
205 bool AbsGeq< std::complex<double> >
206  (const std::complex<double> &%x%, const std::complex<double> &%y%)
207 %\$\$
208 are define by including \$code lu_factor.hpp\$\$
209 These return true if the sum of the square of the real and imaginary parts
210 of \$icode x\$\$ is greater than or equal the
211 sum of the square of the real and imaginary parts of \$icode y\$\$.
212
213 \$children%
214  example/utility/lu_factor.cpp%
215  omh/lu_factor_hpp.omh
216 %\$\$
218 The file
219 \$cref lu_factor.cpp\$\$
220 contains an example and test of using \$code LuFactor\$\$ by itself.
221 It returns true if it succeeds and false otherwise.
222 \$pre
223
224 \$\$
225 The file \$cref lu_solve.hpp\$\$ provides a useful example usage of
226 \$code LuFactor\$\$ with \$code LuInvert\$\$.
227
229 The file \$cref lu_factor.hpp\$\$ contains the
230 current source code that implements these specifications.
231
232 \$end
233 --------------------------------------------------------------------------
234 */
235 // BEGIN C++
236
237 # include <complex>
238 # include <vector>
239
243
245
246 // AbsGeq
247 template <typename Float>
248 inline bool AbsGeq(const Float &x, const Float &y)
249 { Float xabs = x;
250  if( xabs <= Float(0) )
251  xabs = - xabs;
252  Float yabs = y;
253  if( yabs <= Float(0) )
254  yabs = - yabs;
255  return xabs >= yabs;
256 }
257 inline bool AbsGeq(
258  const std::complex<double> &x,
259  const std::complex<double> &y)
260 { double xsq = x.real() * x.real() + x.imag() * x.imag();
261  double ysq = y.real() * y.real() + y.imag() * y.imag();
262
263  return xsq >= ysq;
264 }
265 inline bool AbsGeq(
266  const std::complex<float> &x,
267  const std::complex<float> &y)
268 { float xsq = x.real() * x.real() + x.imag() * x.imag();
269  float ysq = y.real() * y.real() + y.imag() * y.imag();
270
271  return xsq >= ysq;
272 }
273
274 // Lines that are different from code in cppad/core/lu_ratio.hpp end with //
275 template <class SizeVector, class FloatVector> //
276 int LuFactor(SizeVector &ip, SizeVector &jp, FloatVector &LU) //
277 {
278  // type of the elements of LU //
279  typedef typename FloatVector::value_type Float; //
280
281  // check numeric type specifications
282  CheckNumericType<Float>();
283
284  // check simple vector class specifications
285  CheckSimpleVector<Float, FloatVector>();
286  CheckSimpleVector<size_t, SizeVector>();
287
288  size_t i, j; // some temporary indices
289  const Float zero( 0 ); // the value zero as a Float object
290  size_t imax; // row index of maximum element
291  size_t jmax; // column indx of maximum element
292  Float emax; // maximum absolute value
293  size_t p; // count pivots
294  int sign; // sign of the permutation
295  Float etmp; // temporary element
296  Float pivot; // pivot element
297
298  // -------------------------------------------------------
299  size_t n = ip.size();
301  size_t(jp.size()) == n,
302  "Error in LuFactor: jp must have size equal to n"
303  );
305  size_t(LU.size()) == n * n,
306  "Error in LuFactor: LU must have size equal to n * m"
307  );
308  // -------------------------------------------------------
309
310  // initialize row and column order in matrix not yet pivoted
311  for(i = 0; i < n; i++)
312  { ip[i] = i;
313  jp[i] = i;
314  }
315  // initialize the sign of the permutation
316  sign = 1;
317  // ---------------------------------------------------------
318
319  // Reduce the matrix P to L * U using n pivots
320  for(p = 0; p < n; p++)
321  { // determine row and column corresponding to element of
322  // maximum absolute value in remaining part of P
323  imax = jmax = n;
324  emax = zero;
325  for(i = p; i < n; i++)
326  { for(j = p; j < n; j++)
328  (ip[i] < n) & (jp[j] < n)
329  );
330  etmp = LU[ ip[i] * n + jp[j] ];
331
332  // check if maximum absolute value so far
333  if( AbsGeq (etmp, emax) )
334  { imax = i;
335  jmax = j;
336  emax = etmp;
337  }
338  }
339  }
341  (imax < n) & (jmax < n) ,
342  "LuFactor can't determine an element with "
343  "maximum absolute value.\n"
344  "Perhaps original matrix contains not a number or infinity.\n"
345  "Perhaps your specialization of AbsGeq is not correct."
346  );
347  if( imax != p )
348  { // switch rows so max absolute element is in row p
349  i = ip[p];
350  ip[p] = ip[imax];
351  ip[imax] = i;
352  sign = -sign;
353  }
354  if( jmax != p )
355  { // switch columns so max absolute element is in column p
356  j = jp[p];
357  jp[p] = jp[jmax];
358  jp[jmax] = j;
359  sign = -sign;
360  }
361  // pivot using the max absolute element
362  pivot = LU[ ip[p] * n + jp[p] ];
363
364  // check for determinant equal to zero
365  if( pivot == zero )
366  { // abort the mission
367  return 0;
368  }
369
370  // Reduce U by the elementary transformations that maps
371  // LU( ip[p], jp[p] ) to one. Only need transform elements
372  // above the diagonal in U and LU( ip[p] , jp[p] ) is
373  // corresponding value below diagonal in L.
374  for(j = p+1; j < n; j++)
375  LU[ ip[p] * n + jp[j] ] /= pivot;
376
377  // Reduce U by the elementary transformations that maps
378  // LU( ip[i], jp[p] ) to zero. Only need transform elements
379  // above the diagonal in U and LU( ip[i], jp[p] ) is
380  // corresponding value below diagonal in L.
381  for(i = p+1; i < n; i++ )
382  { etmp = LU[ ip[i] * n + jp[p] ];
383  for(j = p+1; j < n; j++)
384  { LU[ ip[i] * n + jp[j] ] -=
385  etmp * LU[ ip[p] * n + jp[j] ];
386  }
387  }
388  }
389  return sign;
390 }
391 } // END CppAD namespace
392 // END C++
393 # endif
Check that exp is true, if not print msg and terminate execution.
bool AbsGeq(const Float &x, const Float &y)
Definition: lu_factor.hpp:248
Define the CppAD error checking macros (all of which begin with CPPAD_ASSERT_)
int LuFactor(SizeVector &ip, SizeVector &jp, FloatVector &LU)
Definition: lu_factor.hpp:276