Prev Next det_of_minor.hpp Headings

@(@\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}} }@)@
Source: det_of_minor
# 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

Input File: omh/det_of_minor_hpp.omh