CppAD: A C++ Algorithmic Differentiation Package  20171217
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
lu_invert.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_UTILITY_LU_INVERT_HPP
2 # define CPPAD_UTILITY_LU_INVERT_HPP
3 
4 /* --------------------------------------------------------------------------
5 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
6 
7 CppAD is distributed under multiple licenses. This distribution is under
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.
12 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
13 -------------------------------------------------------------------------- */
14 
15 /*
16 $begin LuInvert$$
17 $escape #$$
18 $spell
19  cppad.hpp
20  Lu
21  Cpp
22  jp
23  ip
24  const
25  namespace
26  typename
27  etmp
28 $$
29 
30 
31 $section Invert an LU Factored Equation$$
32 $mindex LuInvert linear$$
33 
34 $pre
35 $$
36 
37 $head Syntax$$ $codei%# include <cppad/utility/lu_invert.hpp>
38 %$$
39 $codei%LuInvert(%ip%, %jp%, %LU%, %X%)%$$
40 
41 
42 $head Description$$
43 Solves the matrix equation $icode%A% * %X% = %B%$$
44 using an LU factorization computed by $cref LuFactor$$.
45 
46 $head Include$$
47 The file $code cppad/lu_invert.hpp$$ is included by $code cppad/cppad.hpp$$
48 but it can also be included separately with out the rest of
49 the $code CppAD$$ routines.
50 
51 $head Matrix Storage$$
52 All matrices are stored in row major order.
53 To be specific, if $latex Y$$ is a vector
54 that contains a $latex p$$ by $latex q$$ matrix,
55 the size of $latex Y$$ must be equal to $latex p * q $$ and for
56 $latex i = 0 , \ldots , p-1$$,
57 $latex j = 0 , \ldots , q-1$$,
58 $latex \[
59  Y_{i,j} = Y[ i * q + j ]
60 \] $$
61 
62 $head ip$$
63 The argument $icode ip$$ has prototype
64 $codei%
65  const %SizeVector% &%ip%
66 %$$
67 (see description for $icode SizeVector$$ in
68 $cref/LuFactor/LuFactor/SizeVector/$$ specifications).
69 The size of $icode ip$$ is referred to as $icode n$$ in the
70 specifications below.
71 The elements of $icode ip$$ determine
72 the order of the rows in the permuted matrix.
73 
74 $head jp$$
75 The argument $icode jp$$ has prototype
76 $codei%
77  const %SizeVector% &%jp%
78 %$$
79 (see description for $icode SizeVector$$ in
80 $cref/LuFactor/LuFactor/SizeVector/$$ specifications).
81 The size of $icode jp$$ must be equal to $icode n$$.
82 The elements of $icode jp$$ determine
83 the order of the columns in the permuted matrix.
84 
85 $head LU$$
86 The argument $icode LU$$ has the prototype
87 $codei%
88  const %FloatVector% &%LU%
89 %$$
90 and the size of $icode LU$$ must equal $latex n * n$$
91 (see description for $icode FloatVector$$ in
92 $cref/LuFactor/LuFactor/FloatVector/$$ specifications).
93 
94 $subhead L$$
95 We define the lower triangular matrix $icode L$$ in terms of $icode LU$$.
96 The matrix $icode L$$ is zero above the diagonal
97 and the rest of the elements are defined by
98 $codei%
99  %L%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ]
100 %$$
101 for $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , i$$.
102 
103 $subhead U$$
104 We define the upper triangular matrix $icode U$$ in terms of $icode LU$$.
105 The matrix $icode U$$ is zero below the diagonal,
106 one on the diagonal,
107 and the rest of the elements are defined by
108 $codei%
109  %U%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ]
110 %$$
111 for $latex i = 0 , \ldots , n-2$$ and $latex j = i+1 , \ldots , n-1$$.
112 
113 $subhead P$$
114 We define the permuted matrix $icode P$$ in terms of
115 the matrix $icode L$$ and the matrix $icode U$$
116 by $icode%P% = %L% * %U%$$.
117 
118 $subhead A$$
119 The matrix $icode A$$,
120 which defines the linear equations that we are solving, is given by
121 $codei%
122  %P%(%i%, %j%) = %A%[ %ip%[%i%] * %n% + %jp%[%j%] ]
123 %$$
124 (Hence
125 $icode LU$$ contains a permuted factorization of the matrix $icode A$$.)
126 
127 
128 $head X$$
129 The argument $icode X$$ has prototype
130 $codei%
131  %FloatVector% &%X%
132 %$$
133 (see description for $icode FloatVector$$ in
134 $cref/LuFactor/LuFactor/FloatVector/$$ specifications).
135 The matrix $icode X$$
136 must have the same number of rows as the matrix $icode A$$.
137 The input value of $icode X$$ is the matrix $icode B$$ and the
138 output value solves the matrix equation $icode%A% * %X% = %B%$$.
139 
140 
141 $children%
142  example/utility/lu_invert.cpp%
143  omh/lu_invert_hpp.omh
144 %$$
145 $head Example$$
146 The file $cref lu_solve.hpp$$ is a good example usage of
147 $code LuFactor$$ with $code LuInvert$$.
148 The file
149 $cref lu_invert.cpp$$
150 contains an example and test of using $code LuInvert$$ by itself.
151 It returns true if it succeeds and false otherwise.
152 
153 $head Source$$
154 The file $cref lu_invert.hpp$$ contains the
155 current source code that implements these specifications.
156 
157 $end
158 --------------------------------------------------------------------------
159 */
160 // BEGIN C++
161 # include <cppad/core/cppad_assert.hpp>
164 
165 namespace CppAD { // BEGIN CppAD namespace
166 
167 // LuInvert
168 template <typename SizeVector, typename FloatVector>
169 void LuInvert(
170  const SizeVector &ip,
171  const SizeVector &jp,
172  const FloatVector &LU,
173  FloatVector &B )
174 { size_t k; // column index in X
175  size_t p; // index along diagonal in LU
176  size_t i; // row index in LU and X
177 
178  typedef typename FloatVector::value_type Float;
179 
180  // check numeric type specifications
181  CheckNumericType<Float>();
182 
183  // check simple vector class specifications
184  CheckSimpleVector<Float, FloatVector>();
185  CheckSimpleVector<size_t, SizeVector>();
186 
187  Float etmp;
188 
189  size_t n = ip.size();
191  size_t(jp.size()) == n,
192  "Error in LuInvert: jp must have size equal to n * n"
193  );
195  size_t(LU.size()) == n * n,
196  "Error in LuInvert: Lu must have size equal to n * m"
197  );
198  size_t m = size_t(B.size()) / n;
200  size_t(B.size()) == n * m,
201  "Error in LuSolve: B must have size equal to a multiple of n"
202  );
203 
204  // temporary storage for reordered solution
205  FloatVector x(n);
206 
207  // loop over equations
208  for(k = 0; k < m; k++)
209  { // invert the equation c = L * b
210  for(p = 0; p < n; p++)
211  { // solve for c[p]
212  etmp = B[ ip[p] * m + k ] / LU[ ip[p] * n + jp[p] ];
213  B[ ip[p] * m + k ] = etmp;
214  // subtract off effect on other variables
215  for(i = p+1; i < n; i++)
216  B[ ip[i] * m + k ] -=
217  etmp * LU[ ip[i] * n + jp[p] ];
218  }
219 
220  // invert the equation x = U * c
221  p = n;
222  while( p > 0 )
223  { --p;
224  etmp = B[ ip[p] * m + k ];
225  x[ jp[p] ] = etmp;
226  for(i = 0; i < p; i++ )
227  B[ ip[i] * m + k ] -=
228  etmp * LU[ ip[i] * n + jp[p] ];
229  }
230 
231  // copy reordered solution into B
232  for(i = 0; i < n; i++)
233  B[i * m + k] = x[i];
234  }
235  return;
236 }
237 } // END CppAD namespace
238 // END C++
239 # endif
#define CPPAD_ASSERT_KNOWN(exp, msg)
Check that exp is true, if not print msg and terminate execution.
Define the CppAD error checking macros (all of which begin with CPPAD_ASSERT_)
void LuInvert(const SizeVector &ip, const SizeVector &jp, const FloatVector &LU, FloatVector &B)
Definition: lu_invert.hpp:169
Scalar value_type