$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$
# ifndef CPPAD_DET_OF_MINOR_HPP # define CPPAD_DET_OF_MINOR_HPP # include <vector> # include <cstddef> namespace CppAD { // BEGIN CppAD namespace template <class Scalar> 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 ) { const size_t R0 = r[m]; // R(0) size_t Cj = c[m]; // C(j) (case j = 0) size_t Cj1 = m; // C(j-1) (case j = 0) // check for 1 by 1 case if( n == 1 ) return a[ R0 * m + Cj ]; // initialize determinant of the minor M Scalar detM = Scalar(0); // initialize sign of factor for next sub-minor int s = 1; // remove row with index 0 in M from all the sub-minors of M r[m] = r[R0]; // for each column of M for(size_t j = 0; j < n; j++) { // element with index (0,j) in the minor M Scalar M0j = a[ R0 * m + Cj ]; // remove column with index j in M to form next sub-minor S of M c[Cj1] = c[Cj]; // compute determinant of the current sub-minor S Scalar detS = det_of_minor(a, m, n - 1, r, c); // restore column Cj to represenation of M as a minor of A c[Cj1] = Cj; // include this sub-minor term in the summation if( s > 0 ) detM = detM + M0j * detS; else detM = detM - M0j * detS; // advance to next column of M Cj1 = Cj; Cj = c[Cj]; s = - s; } // restore row zero to the minor representation for M r[m] = R0; // return the determinant of the minor M return detM; } } // END CppAD namespace # endif