CppAD: A C++ Algorithmic Differentiation Package  20171217
lu_invert.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 LuInvert\$\$
17 \$escape #\$\$
18 \$spell
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
38 %\$\$
39 \$codei%LuInvert(%ip%, %jp%, %LU%, %X%)%\$\$
40
41
43 Solves the matrix equation \$icode%A% * %X% = %B%\$\$
44 using an LU factorization computed by \$cref LuFactor\$\$.
45
48 but it can also be included separately with out the rest of
50
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
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
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
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
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
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
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
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
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 %\$\$
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
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++
164
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