CppAD: A C++ Algorithmic Differentiation Package  20171217
lu_solve.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 LuSolve\$\$
17 \$escape #\$\$
18 \$spell
20  det
21  exp
22  Leq
23  typename
24  bool
25  const
26  namespace
27  std
28  Geq
29  Lu
31  signdet
32  logdet
33 \$\$
34
35
36 \$section Compute Determinant and Solve Linear Equations\$\$
37 \$mindex LuSolve Lu\$\$
38
39 \$pre
40 \$\$
41
43 %\$\$
44 \$icode%signdet% = LuSolve(%n%, %m%, %A%, %B%, %X%, %logdet%)%\$\$
45
46
48 Use an LU factorization of the matrix \$icode A\$\$ to
49 compute its determinant
50 and solve for \$icode X\$\$ in the linear of equation
51 \$latex \[
52  A * X = B
53 \] \$\$
54 where \$icode A\$\$ is an
55 \$icode n\$\$ by \$icode n\$\$ matrix,
56 \$icode X\$\$ is an
57 \$icode n\$\$ by \$icode m\$\$ matrix, and
58 \$icode B\$\$ is an \$latex n x m\$\$ matrix.
59
62 but it can also be included separately with out the rest of
64
66 This routine is an easy to user interface to
67 \$cref LuFactor\$\$ and \$cref LuInvert\$\$ for computing determinants and
68 solutions of linear equations.
69 These separate routines should be used if
70 one right hand side \$icode B\$\$
71 depends on the solution corresponding to another
72 right hand side (with the same value of \$icode A\$\$).
73 In this case only one call to \$code LuFactor\$\$ is required
74 but there will be multiple calls to \$code LuInvert\$\$.
75
76
78 All matrices are stored in row major order.
79 To be specific, if \$latex Y\$\$ is a vector
80 that contains a \$latex p\$\$ by \$latex q\$\$ matrix,
81 the size of \$latex Y\$\$ must be equal to \$latex p * q \$\$ and for
82 \$latex i = 0 , \ldots , p-1\$\$,
83 \$latex j = 0 , \ldots , q-1\$\$,
84 \$latex \[
85  Y_{i,j} = Y[ i * q + j ]
86 \] \$\$
87
89 The return value \$icode signdet\$\$ is a \$code int\$\$ value
90 that specifies the sign factor for the determinant of \$icode A\$\$.
91 This determinant of \$icode A\$\$ is zero if and only if \$icode signdet\$\$
92 is zero.
93
95 The argument \$icode n\$\$ has type \$code size_t\$\$
96 and specifies the number of rows in the matrices
97 \$icode A\$\$,
98 \$icode X\$\$,
99 and \$icode B\$\$.
100 The number of columns in \$icode A\$\$ is also equal to \$icode n\$\$.
101
103 The argument \$icode m\$\$ has type \$code size_t\$\$
104 and specifies the number of columns in the matrices
105 \$icode X\$\$
106 and \$icode B\$\$.
107 If \$icode m\$\$ is zero,
108 only the determinant of \$icode A\$\$ is computed and
109 the matrices \$icode X\$\$ and \$icode B\$\$ are not used.
110
112 The argument \$icode A\$\$ has the prototype
113 \$codei%
114  const %FloatVector% &%A%
115 %\$\$
116 and the size of \$icode A\$\$ must equal \$latex n * n\$\$
117 (see description of \$cref/FloatVector/LuSolve/FloatVector/\$\$ below).
118 This is the \$latex n\$\$ by \$icode n\$\$ matrix that
119 we are computing the determinant of
120 and that defines the linear equation.
121
123 The argument \$icode B\$\$ has the prototype
124 \$codei%
125  const %FloatVector% &%B%
126 %\$\$
127 and the size of \$icode B\$\$ must equal \$latex n * m\$\$
128 (see description of \$cref/FloatVector/LuSolve/FloatVector/\$\$ below).
129 This is the \$latex n\$\$ by \$icode m\$\$ matrix that
130 defines the right hand side of the linear equations.
131 If \$icode m\$\$ is zero, \$icode B\$\$ is not used.
132
134 The argument \$icode X\$\$ has the prototype
135 \$codei%
136  %FloatVector% &%X%
137 %\$\$
138 and the size of \$icode X\$\$ must equal \$latex n * m\$\$
139 (see description of \$cref/FloatVector/LuSolve/FloatVector/\$\$ below).
140 The input value of \$icode X\$\$ does not matter.
141 On output, the elements of \$icode X\$\$ contain the solution
142 of the equation we wish to solve
143 (unless \$icode signdet\$\$ is equal to zero).
144 If \$icode m\$\$ is zero, \$icode X\$\$ is not used.
145
147 The argument \$icode logdet\$\$ has prototype
148 \$codei%
149  %Float% &%logdet%
150 %\$\$
151 On input, the value of \$icode logdet\$\$ does not matter.
152 On output, it has been set to the
153 log of the determinant of \$icode A\$\$
154 (but not quite).
155 To be more specific,
156 the determinant of \$icode A\$\$ is given by the formula
157 \$codei%
158  %det% = %signdet% * exp( %logdet% )
159 %\$\$
160 This enables \$code LuSolve\$\$ to use logs of absolute values
161 in the case where \$icode Float\$\$ corresponds to a real number.
162
164 The type \$icode Float\$\$ must satisfy the conditions
165 for a \$cref NumericType\$\$ type.
166 The routine \$cref CheckNumericType\$\$ will generate an error message
167 if this is not the case.
168 In addition, the following operations must be defined for any pair
169 of \$icode Float\$\$ objects \$icode x\$\$ and \$icode y\$\$:
170
171 \$table
172 \$bold Operation\$\$ \$cnext \$bold Description\$\$ \$rnext
173 \$codei%log(%x%)%\$\$ \$cnext
174  returns the logarithm of \$icode x\$\$ as a \$icode Float\$\$ object
175 \$tend
176
178 The type \$icode FloatVector\$\$ must be a \$cref SimpleVector\$\$ class with
179 \$cref/elements of type Float/SimpleVector/Elements of Specified Type/\$\$.
180 The routine \$cref CheckSimpleVector\$\$ will generate an error message
181 if this is not the case.
182
184 Including the file \$code lu_solve.hpp\$\$ defines the template function
185 \$codei%
186  template <typename %Float%>
187  bool LeqZero<%Float%>(const %Float% &%x%)
188 %\$\$
189 in the \$code CppAD\$\$ namespace.
190 This function returns true if \$icode x\$\$ is less than or equal to zero
191 and false otherwise.
192 It is used by \$code LuSolve\$\$ to avoid taking the log of
193 zero (or a negative number if \$icode Float\$\$ corresponds to real numbers).
194 This template function definition assumes that the operator
195 \$code <=\$\$ is defined for \$icode Float\$\$ objects.
196 If this operator is not defined for your use of \$icode Float\$\$,
197 you will need to specialize this template so that it works for your
198 use of \$code LuSolve\$\$.
199 \$pre
200
201 \$\$
202 Complex numbers do not have the operation or \$code <=\$\$ defined.
203 In addition, in the complex case,
204 one can take the log of a negative number.
205 The specializations
206 \$codei%
207  bool LeqZero< std::complex<float> > (const std::complex<float> &%x%)
208  bool LeqZero< std::complex<double> >(const std::complex<double> &%x%)
209 %\$\$
210 are defined by including \$code lu_solve.hpp\$\$.
211 These return true if \$icode x\$\$ is zero and false otherwise.
212
214 Including the file \$code lu_solve.hpp\$\$ defines the template function
215 \$codei%
216  template <typename %Float%>
217  bool AbsGeq<%Float%>(const %Float% &%x%, const %Float% &%y%)
218 %\$\$
219 If the type \$icode Float\$\$ does not support the \$code <=\$\$ operation
220 and it is not \$code std::complex<float>\$\$ or \$code std::complex<double>\$\$,
221 see the documentation for \$code AbsGeq\$\$ in \$cref/LuFactor/LuFactor/AbsGeq/\$\$.
222
223 \$children%
224  example/utility/lu_solve.cpp%
225  omh/lu_solve_hpp.omh
226 %\$\$
228 The file
229 \$cref lu_solve.cpp\$\$
230 contains an example and test of using this routine.
231 It returns true if it succeeds and false otherwise.
232
234 The file \$cref lu_solve.hpp\$\$ contains the
235 current source code that implements these specifications.
236
237 \$end
238 --------------------------------------------------------------------------
239 */
240 // BEGIN C++
241 # include <complex>
242 # include <vector>
243
244 // link exp for float and double cases
246
252
254
255 // LeqZero
256 template <typename Float>
257 inline bool LeqZero(const Float &x)
258 { return x <= Float(0); }
259 inline bool LeqZero( const std::complex<double> &x )
260 { return x == std::complex<double>(0); }
261 inline bool LeqZero( const std::complex<float> &x )
262 { return x == std::complex<float>(0); }
263
264 // LuSolve
265 template <typename Float, typename FloatVector>
267  size_t n ,
268  size_t m ,
269  const FloatVector &A ,
270  const FloatVector &B ,
271  FloatVector &X ,
272  Float &logdet )
273 {
274  // check numeric type specifications
275  CheckNumericType<Float>();
276
277  // check simple vector class specifications
278  CheckSimpleVector<Float, FloatVector>();
279
280  size_t p; // index of pivot element (diagonal of L)
281  int signdet; // sign of the determinant
282  Float pivot; // pivot element
283
284  // the value zero
285  const Float zero(0);
286
287  // pivot row and column order in the matrix
288  std::vector<size_t> ip(n);
289  std::vector<size_t> jp(n);
290
291  // -------------------------------------------------------
293  size_t(A.size()) == n * n,
294  "Error in LuSolve: A must have size equal to n * n"
295  );
297  size_t(B.size()) == n * m,
298  "Error in LuSolve: B must have size equal to n * m"
299  );
301  size_t(X.size()) == n * m,
302  "Error in LuSolve: X must have size equal to n * m"
303  );
304  // -------------------------------------------------------
305
306  // copy A so that it does not change
307  FloatVector Lu(A);
308
309  // copy B so that it does not change
310  X = B;
311
312  // Lu factor the matrix A
313  signdet = LuFactor(ip, jp, Lu);
314
315  // compute the log of the determinant
316  logdet = Float(0);
317  for(p = 0; p < n; p++)
318  { // pivot using the max absolute element
319  pivot = Lu[ ip[p] * n + jp[p] ];
320
321  // check for determinant equal to zero
322  if( pivot == zero )
323  { // abort the mission
324  logdet = Float(0);
325  return 0;
326  }
327
328  // update the determinant
329  if( LeqZero ( pivot ) )
330  { logdet += log( - pivot );
331  signdet = - signdet;
332  }
333  else logdet += log( pivot );
334
335  }
336
337  // solve the linear equations
338  LuInvert(ip, jp, Lu, X);
339
340  // return the sign factor for the determinant
341  return signdet;
342 }
343 } // END CppAD namespace
344 // END C++
345 # endif