CppAD: A C++ Algorithmic Differentiation Package  20171217
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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 
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 $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)