00001 # ifndef CPPAD_DET_OF_MINOR_INCLUDED 00002 # define CPPAD_DET_OF_MINOR_INCLUDED 00003 00004 /* -------------------------------------------------------------------------- 00005 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell 00006 00007 CppAD is distributed under multiple licenses. This distribution is under 00008 the terms of the 00009 Common Public License Version 1.0. 00010 00011 A copy of this license is included in the COPYING file of this distribution. 00012 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00013 -------------------------------------------------------------------------- */ 00014 /* 00015 $begin det_of_minor$$ 00016 $spell 00017 cppad.hpp 00018 hpp 00019 std 00020 Det 00021 const 00022 CppADvector 00023 namespace 00024 $$ 00025 00026 $index det_of_minor$$ 00027 $index determinant, matrix minor$$ 00028 $index matrix, minor determinant$$ 00029 $index minor, matrix determinant$$ 00030 00031 $section Determinant of a Minor$$ 00032 00033 $head Syntax$$ 00034 $syntax%# include <cppad/speed/det_of_minor.hpp> 00035 %$$ 00036 $syntax%%d% = det_of_minor(%a%, %m%, %n%, %r%, %c%)%$$ 00037 00038 00039 $head Inclusion$$ 00040 The template function $code det_of_minor$$ is defined in the $code CppAD$$ 00041 namespace by including 00042 the file $code cppad/speed/det_of_minor.hpp$$ 00043 (relative to the CppAD distribution directory). 00044 It is only intended for example and testing purposes, 00045 so it is not automatically included by 00046 $cref/cppad.hpp/cppad/$$. 00047 00048 $head Purpose$$ 00049 This template function 00050 returns the determinant of a minor of the matrix $latex A$$ 00051 using expansion by minors. 00052 The elements of the $latex n \times n$$ minor $latex M$$ 00053 of the matrix $latex A$$ are defined, 00054 for $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , n-1$$, by 00055 $latex \[ 00056 M_{i,j} = A_{R(i), C(j)} 00057 \]$$ 00058 where the functions 00059 $latex R(i)$$ is defined by the $cref/argument r/det_of_minor/r/$$ and 00060 $latex C(j)$$ is defined by the $cref/argument c/det_of_minor/c/$$. 00061 $pre 00062 00063 $$ 00064 This template function 00065 is for example and testing purposes only. 00066 Expansion by minors is chosen as an example because it uses 00067 a lot of floating point operations yet does not require much source code 00068 (on the order of $italic m$$ factorial floating point operations and 00069 about 70 lines of source code including comments). 00070 This is not an efficient method for computing a determinant; 00071 for example, using an LU factorization would be better. 00072 00073 $head Determinant of A$$ 00074 If the following conditions hold, the minor is the 00075 entire matrix $latex A$$ and hence $code det_of_minor$$ 00076 will return the determinant of $latex A$$: 00077 00078 $list number$$ 00079 $latex n = m$$. 00080 $lnext 00081 for $latex i = 0 , \ldots , m-1$$, $latex r[i] = i+1$$, 00082 and $latex r[m] = 0$$. 00083 $lnext 00084 for $latex j = 0 , \ldots , m-1$$, $latex c[j] = j+1$$, 00085 and $latex c[m] = 0$$. 00086 $lend 00087 00088 $head a$$ 00089 The argument $italic a$$ has prototype 00090 $syntax% 00091 const std::vector<%Scalar%>& %a% 00092 %$$ 00093 and is a vector with size $latex m * m$$ 00094 (see description of $cref/Scalar/det_of_minor/Scalar/$$ below). 00095 The elements of the $latex m \times m$$ matrix $latex A$$ are defined, 00096 for $latex i = 0 , \ldots , m-1$$ and $latex j = 0 , \ldots , m-1$$, by 00097 $latex \[ 00098 A_{i,j} = a[ i + j * m ] 00099 \] $$ 00100 00101 $head m$$ 00102 The argument $italic m$$ has prototype 00103 $syntax% 00104 size_t %m% 00105 %$$ 00106 and is the size of the square matrix $latex A$$. 00107 00108 $head n$$ 00109 The argument $italic n$$ has prototype 00110 $syntax% 00111 size_t %n% 00112 %$$ 00113 and is the size of the square minor $latex M$$. 00114 00115 $head r$$ 00116 The argument $italic r$$ has prototype 00117 $syntax% 00118 std::vector<size_t>& %r% 00119 %$$ 00120 and is a vector with $latex m + 1$$ elements. 00121 This vector defines the function $latex R(i)$$ 00122 which specifies the rows of the minor $latex M$$. 00123 To be specific, the function $latex R(i)$$ 00124 for $latex i = 0, \ldots , n-1$$ is defined by 00125 $latex \[ 00126 \begin{array}{rcl} 00127 R(0) & = & r[m] 00128 \\ 00129 R(i+1) & = & r[ R(i) ] 00130 \end{array} 00131 \] $$ 00132 All the elements of $italic r$$ must have value 00133 less than or equal $italic m$$. 00134 The elements of vector $italic r$$ are modified during the computation, 00135 and restored to their original value before the return from 00136 $code det_of_minor$$. 00137 00138 $head c$$ 00139 The argument $italic c$$ has prototype 00140 $syntax% 00141 std::vector<size_t>& %c% 00142 %$$ 00143 and is a vector with $latex m + 1$$ elements 00144 This vector defines the function $latex C(i)$$ 00145 which specifies the rows of the minor $latex M$$. 00146 To be specific, the function $latex C(i)$$ 00147 for $latex j = 0, \ldots , n-1$$ is defined by 00148 $latex \[ 00149 \begin{array}{rcl} 00150 C(0) & = & c[m] 00151 \\ 00152 C(j+1) & = & c[ C(j) ] 00153 \end{array} 00154 \] $$ 00155 All the elements of $italic c$$ must have value 00156 less than or equal $italic m$$. 00157 The elements of vector $italic c$$ are modified during the computation, 00158 and restored to their original value before the return from 00159 $code det_of_minor$$. 00160 00161 $head d$$ 00162 The result $italic d$$ has prototype 00163 $syntax% 00164 %Scalar% %d% 00165 %$$ 00166 and is equal to the determinant of the minor $latex M$$. 00167 00168 $head Scalar$$ 00169 If $italic x$$ and $italic y$$ are objects of type $italic Scalar$$ 00170 and $italic i$$ is an object of type $code int$$, 00171 the $italic Scalar$$ must support the following operations: 00172 $table 00173 $bold Syntax$$ 00174 $cnext $bold Description$$ 00175 $cnext $bold Result Type$$ 00176 $rnext 00177 $syntax%%Scalar% %x%$$ 00178 $cnext default constructor for $italic Scalar$$ object. 00179 $rnext 00180 $syntax%%Scalar% %x%(%i%)%$$ 00181 $cnext construct $italic x$$ with value corresponding to $italic i$$ 00182 $rnext 00183 $syntax%%x% = %y%$$ 00184 $cnext set value of $italic x$$ to current value of $italic y$$ 00185 $rnext 00186 $syntax%%x% + %y%$$ 00187 $cnext value of $italic x$$ plus of $italic y$$ 00188 $cnext $italic Scalar$$ 00189 $rnext 00190 $syntax%%x% - %y%$$ 00191 $cnext value of $italic x$$ minus of $italic y$$ 00192 $cnext $italic Scalar$$ 00193 $rnext 00194 $syntax%%x% * %y%$$ 00195 $cnext value of $italic x$$ times value of $italic y$$ 00196 $cnext $italic Scalar$$ 00197 $tend 00198 00199 $children% 00200 speed/example/det_of_minor.cpp% 00201 omh/det_of_minor_hpp.omh 00202 %$$ 00203 00204 $head Example$$ 00205 The file 00206 $cref/det_of_minor.cpp/$$ 00207 contains an example and test of $code det_of_minor.hpp$$. 00208 It returns true if it succeeds and false otherwise. 00209 00210 $head Source Code$$ 00211 The file 00212 $xref/det_of_minor.hpp/$$ 00213 contains the source for this template function. 00214 00215 00216 $end 00217 --------------------------------------------------------------------------- 00218 */ 00219 // BEGIN PROGRAM 00220 namespace CppAD { // BEGIN CppAD namespace 00221 template <class Scalar> 00222 Scalar det_of_minor( 00223 const std::vector<Scalar>& a , 00224 size_t m , 00225 size_t n , 00226 std::vector<size_t>& r , 00227 std::vector<size_t>& c ) 00228 { size_t R0; // R(0) 00229 size_t Cj; // C(j) 00230 size_t Cj1; // C(j-1) 00231 size_t j; // column index in M 00232 Scalar detS; // determinant of sub-minor of M 00233 Scalar M0j; // an element of the first row of M 00234 int s; // sign of the current sub-minor in summation 00235 00236 // value of R(0) and C(0) 00237 R0 = r[m]; 00238 Cj = c[m]; 00239 00240 // check for 1 by 1 case 00241 if( n == 1 ) 00242 return a[ R0 * m + Cj ]; 00243 00244 // initialize determinant of the minor M 00245 Scalar detM(0); 00246 00247 // initialize sign of permutation 00248 s = 1; 00249 00250 // remove row with index R0 from sub-minors of M 00251 r[m] = r[R0]; 00252 00253 // initialize previous column index; Cj = c[ Cj1 ] 00254 Cj1 = m; 00255 00256 // for each column of M 00257 for(j = 0; j < n; j++) 00258 { // M(0, j) 00259 M0j = a[ R0 * m + Cj ]; 00260 00261 // remove column index j from sub-minor of M 00262 c[Cj1] = c[Cj]; 00263 00264 // compute determinant of sub-minor of M 00265 // where row R0 and column Cj are removed 00266 detS = det_of_minor(a, m, n - 1, r, c); 00267 00268 // restore column Cj to sub-minor of M 00269 c[Cj1] = Cj; 00270 00271 // include in summation 00272 if( s > 0 ) 00273 detM = detM + M0j * detS; 00274 else detM = detM - M0j * detS; 00275 00276 // advance to next column of M 00277 Cj1 = Cj; 00278 Cj = c[Cj]; 00279 s = -s; 00280 } 00281 00282 // restore row index zero to the minor M 00283 r[m] = R0; 00284 00285 // return the determinant of the minor M 00286 return detM; 00287 } 00288 } // END CppAD namespace 00289 // END PROGRAM 00290 # endif