\(\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}} }\)

det_of_minor.hpp

View page source

Source: det_of_minor

# ifndef CPPAD_DET_OF_MINOR_HPP

# define CPPAD_DET_OF_MINOR_HPP

# include <vector>
# include <cassert>
# include <cstddef>

namespace CppAD { // BEGIN CppAD namespace

// BEGIN_DET_OF_MINOR
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  )
{  assert( a.size() == m * m );
   assert( r.size() == m + 1 );
   assert( c.size() == m + 1 );
   // END_DET_OF_MINOR
   //
   // R0 = R(0)
   size_t R0 = r[m];
   assert( R0 < m );
   //
   // Cj = C(0)
   size_t Cj = c[m];
   assert( Cj < m );
   //
   //
   // check if this is a 1 by 1 minor
   if( n == 1 ) return a[ R0 * m + Cj ];
   //
   // detM
   // initialize determinant of the minor M
   Scalar detM(0);
   //
   // sign
   // initialize sign of factor for next sub-minor
   int sign = 1;
   //
   // r
   // remove row with index 0 in M from all the sub-minors of M
   r[m] = r[R0];
   //
   // C(j-1)
   // initial index in c for previous column of the minor M
   size_t Cj1 = m;
   //
   // for each column of M
   for(size_t j = 0; j < n; j++)
   {
      // M[0,j] = A[ R0, Cj ]
      // element with index (0, j) in the minor M
      assert( Cj < 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];
      //
      // detS
      // compute determinant of S, the sub-minor of M with
      // row R(0) and column C(j) removed.
      Scalar detS = det_of_minor(a, m, n - 1, r, c);
      //
      // restore column with index j in represenation of M as a minor of A
      c[Cj1] = Cj;
      //
      // detM
      // include this sub-minor term in the summation
      if( sign > 0 )
         detM = detM + M0j * detS;
      else
         detM = detM - M0j * detS;
      //
      // advance to next column of M
      Cj1  = Cj;
      Cj   = c[Cj];
      sign = - sign;
   }
   //
   // r
   // 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