\(\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 representation 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