CppAD: A C++ Algorithmic Differentiation Package  20171217
det_of_minor.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_SPEED_DET_OF_MINOR_HPP
2 # define CPPAD_SPEED_DET_OF_MINOR_HPP
3
4 /* --------------------------------------------------------------------------
5 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
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 $begin det_of_minor$$16 spell 17 CppAD 18 hpp 19 std 20 Det 21 const 22 namespace 23 cppad 24$$ 25 26 27$section Determinant of a Minor$$28 mindex det_of_minor matrix$$
29
30 $head Syntax$$31 codei%# include <cppad/speed/det_of_minor.hpp> 32 %$$ 33$icode%d% = det_of_minor(%a%, %m%, %n%, %r%, %c%)%$$34 35 36 head Inclusion$$
37 The template function $code det_of_minor$$is defined in the code CppAD$$ 38 namespace by including 39 the file$code cppad/speed/det_of_minor.hpp$$40 (relative to the CppAD distribution directory). 41 42 head Purpose$$
43 This template function
44 returns the determinant of a minor of the matrix $latex A$$45 using expansion by minors. 46 The elements of the latex n \times n$$ minor$latex M$$47 of the matrix latex A$$ are defined,
48 for $latex i = 0 , \ldots , n-1$$and latex j = 0 , \ldots , n-1$$, by 49$latex $50 M_{i,j} = A_{R(i), C(j)} 51$$$52 where the functions 53 latex R(i)$$ is defined by the $cref/argument r/det_of_minor/r/$$and 54 latex C(j)$$ is defined by the$cref/argument c/det_of_minor/c/$$. 55 pre 56 57$$
58 This template function
59 is for example and testing purposes only.
60 Expansion by minors is chosen as an example because it uses
61 a lot of floating point operations yet does not require much source code
62 (on the order of $icode m$$factorial floating point operations and 63 about 70 lines of source code including comments). 64 This is not an efficient method for computing a determinant; 65 for example, using an LU factorization would be better. 66 67 head Determinant of A$$ 68 If the following conditions hold, the minor is the 69 entire matrix$latex A$$and hence code det_of_minor$$
70 will return the determinant of $latex A$$: 71 72 list number$$ 73$latex n = m$$. 74 lnext 75 for latex i = 0 , \ldots , m-1$$, $latex r[i] = i+1$$, 76 and latex r[m] = 0$$. 77$lnext
78 for $latex j = 0 , \ldots , m-1$$, latex c[j] = j+1$$, 79 and$latex c[m] = 0$$. 80 lend 81 82 head a$$
83 The argument $icode a$$has prototype 84 codei% 85 const std::vector<%Scalar%>& %a% 86 %$$ 87 and is a vector with size$latex m * m$$88 (see description of cref/Scalar/det_of_minor/Scalar/$$ below).
89 The elements of the $latex m \times m$$matrix latex A$$ are defined, 90 for$latex i = 0 , \ldots , m-1$$and latex j = 0 , \ldots , m-1$$, by
91 $latex $92 A_{i,j} = a[ i * m + j] 93$ $$94 95 head m$$ 96 The argument$icode m$$has prototype 97 codei% 98 size_t %m% 99 %$$
100 and is the number of rows (and columns) in the square matrix $latex A$$. 101 102 head n$$ 103 The argument$icode n$$has prototype 104 codei% 105 size_t %n% 106 %$$
107 and is the number of rows (and columns) in the square minor $latex M$$. 108 109 head r$$ 110 The argument$icode r$$has prototype 111 codei% 112 std::vector<size_t>& %r% 113 %$$
114 and is a vector with $latex m + 1$$elements. 115 This vector defines the function latex R(i)$$ 116 which specifies the rows of the minor$latex M$$. 117 To be specific, the function latex R(i)$$
118 for $latex i = 0, \ldots , n-1$$is defined by 119 latex $120 \begin{array}{rcl} 121 R(0) & = & r[m] 122 \\ 123 R(i+1) & = & r[ R(i) ] 124 \end{array} 125$$$ 126 All the elements of$icode r$$must have value 127 less than or equal icode m$$.
128 The elements of vector $icode r$$are modified during the computation, 129 and restored to their original value before the return from 130 code det_of_minor$$. 131 132$head c$$133 The argument icode c$$ has prototype
134 $codei% 135 std::vector<size_t>& %c% 136 %$$137 and is a vector with latex m + 1$$ elements 138 This vector defines the function$latex C(i)$$139 which specifies the rows of the minor latex M$$.
140 To be specific, the function $latex C(i)$$141 for latex j = 0, \ldots , n-1$$ is defined by 142$latex $143 \begin{array}{rcl} 144 C(0) & = & c[m] 145 \\ 146 C(j+1) & = & c[ C(j) ] 147 \end{array} 148$ $$149 All the elements of icode c$$ must have value
150 less than or equal $icode m$$. 151 The elements of vector icode c$$ are modified during the computation, 152 and restored to their original value before the return from 153$code det_of_minor$$. 154 155 head d$$
156 The result $icode d$$has prototype 157 codei% 158 %Scalar% %d% 159 %$$ 160 and is equal to the determinant of the minor$latex M$$. 161 162 head Scalar$$
163 If $icode x$$and icode y$$ are objects of type$icode Scalar$$164 and icode i$$ is an object of type $code int$$, 165 the icode Scalar$$ must support the following operations: 166$table
167 $bold Syntax$$168 cnext bold Description$$ 169$cnext $bold Result Type$$170 rnext 171 icode%Scalar% %x%$$ 172$cnext default constructor for $icode Scalar$$object. 173 rnext 174 icode%x% = %i%$$ 175$cnext set value of $icode x$$to current value of icode i$$ 176$rnext
177 $icode%x% = %y%$$178 cnext set value of icode x$$ to current value of$icode y$$179 rnext 180 icode%x% + %y%$$
181  $cnext value of$icode x$$plus icode y$$
182  $cnext$icode Scalar$$183 rnext 184 icode%x% - %y%$$
185  $cnext value of$icode x$$minus icode y$$
186  $cnext$icode Scalar$$187 rnext 188 icode%x% * %y%$$
189  $cnext value of$icode x$$times value of icode y$$
190  $cnext$icode Scalar$$191 tend 192 193 children% 194 speed/example/det_of_minor.cpp% 195 omh/det_of_minor_hpp.omh 196 %$$
197
198 $head Example$$199 The file 200 cref det_of_minor.cpp$$ 201 contains an example and test of$code det_of_minor.hpp$$. 202 It returns true if it succeeds and false otherwise. 203 204 head Source Code$$
205 The file
206 $cref det_of_minor.hpp$$207 contains the source for this template function. 208 209 210$end
211 ---------------------------------------------------------------------------
212 */
213 // BEGIN C++
214 # include <vector>
215 # include <cstddef>
216
217 namespace CppAD { // BEGIN CppAD namespace
218 template <class Scalar>
220  const std::vector<Scalar>& a ,
221  size_t m ,
222  size_t n ,
223  std::vector<size_t>& r ,
224  std::vector<size_t>& c )
225 {
226  const size_t R0 = r[m]; // R(0)
227  size_t Cj = c[m]; // C(j) (case j = 0)
228  size_t Cj1 = m; // C(j-1) (case j = 0)
229
230  // check for 1 by 1 case
231  if( n == 1 ) return a[ R0 * m + Cj ];
232
233  // initialize determinant of the minor M
234  Scalar detM = Scalar(0);
235
236  // initialize sign of factor for next sub-minor
237  int s = 1;
238
239  // remove row with index 0 in M from all the sub-minors of M
240  r[m] = r[R0];
241
242  // for each column of M
243  for(size_t j = 0; j < n; j++)
244  { // element with index (0,j) in the minor M
245  Scalar M0j = a[ R0 * m + Cj ];
246
247  // remove column with index j in M to form next sub-minor S of M
248  c[Cj1] = c[Cj];
249
250  // compute determinant of the current sub-minor S
251  Scalar detS = det_of_minor(a, m, n - 1, r, c);
252
253  // restore column Cj to represenation of M as a minor of A
254  c[Cj1] = Cj;
255
256  // include this sub-minor term in the summation
257  if( s > 0 )
258  detM = detM + M0j * detS;
259  else detM = detM - M0j * detS;
260
261  // advance to next column of M
262  Cj1 = Cj;
263  Cj = c[Cj];
264  s = - s;
265  }
266
267  // restore row zero to the minor representation for M
268  r[m] = R0;
269
270  // return the determinant of the minor M
271  return detM;
272 }
273 } // END CppAD namespace
274 // END C++
275 # endif
Scalar det_of_minor(const std::vector< Scalar > &a, size_t m, size_t n, std::vector< size_t > &r, std::vector< size_t > &c)