det_of_minor.hpp

View page source

det_of_minor: Source Code

# include <cassert>
# include <cstddef>
# include <cmpad/vector.hpp>

namespace cmpad { // BEGIN cmpad namespace

// BEGIN PROTOTYPE
template <class Vector>
typename Vector::value_type det_of_minor(
   const Vector&                   a  ,
   size_t                          n  ,
   size_t                          m  ,
   cmpad::vector<size_t>&          r  ,
   cmpad::vector<size_t>&          c  )
{  assert( a.size() == n * n );
   assert( r.size() == n + 1 );
   assert( c.size() == n + 1 );
   // END PROTOTYPE
   //
   // scalar_type
   typedef typename Vector::value_type scalar_type;
   //
   // R0 = R(0)
   size_t R0 = r[n];
   assert( R0 < n );
   //
   // Cj = C(0)
   size_t Cj = c[n];
   assert( Cj < n );
   //
   //
   // check if this is a 1 by 1 minor
   if( m == 1 ) return a[ R0 * n + Cj ];
   //
   // detM
   // initialize determinant of the minor M
   scalar_type 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[n] = r[R0];
   //
   // C(j-1)
   // initial index in c for previous column of the minor M
   size_t Cj1 = n;
   //
   // for each column of M
   for(size_t j = 0; j < m; j++)
   {
      // M[0,j] = A[ R0, Cj ]
      // element with index (0, j) in the minor M
      assert( Cj < n );
      scalar_type M0j = a[ R0 * n + 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_type detS = det_of_minor(a, n, m - 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[n] = R0;
   //
   // return the determinant of the minor M
   return detM;
}

} // END cmpad namespace