text: // $Id: det_of_minor.hpp 3757 2015-11-30 12:03:07Z bradbell $
  # ifndef CPPAD_DET_OF_MINOR_HPP
  # define CPPAD_DET_OF_MINOR_HPP
  
  /* --------------------------------------------------------------------------
  CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-15 Bradley M. Bell
  
  CppAD is distributed under multiple licenses. This distribution is under
  the terms of the
                      Eclipse Public License Version 1.0.
  
  A copy of this license is included in the COPYING file of this distribution.
  Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
  -------------------------------------------------------------------------- */
  /*
  $begin det_of_minor$$
  $spell
  	CppAD
  	cppad.hpp
  	hpp
  	std
  	Det
  	const
  	namespace
  $$
  
  
  $section Determinant of a Minor$$
  $mindex det_of_minor matrix$$
  
  $head Syntax$$
  $codei%# include <cppad/speed/det_of_minor.hpp>
  %$$
  $icode%d% = det_of_minor(%a%, %m%, %n%, %r%, %c%)%$$
  
  
  $head Inclusion$$
  The template function $code det_of_minor$$ is defined in the $code CppAD$$
  namespace by including
  the file $code cppad/speed/det_of_minor.hpp$$
  (relative to the CppAD distribution directory).
  It is only intended for example and testing purposes,
  so it is not automatically included by
  $cref/cppad.hpp/cppad/$$.
  
  $head Purpose$$
  This template function
  returns the determinant of a minor of the matrix $latex A$$
  using expansion by minors.
  The elements of the $latex n \times n$$ minor $latex M$$
  of the matrix $latex A$$ are defined,
  for $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , n-1$$, by
  $latex \[
  	M_{i,j} = A_{R(i), C(j)}
  \]$$
  where the functions
  $latex R(i)$$ is defined by the $cref/argument r/det_of_minor/r/$$ and
  $latex C(j)$$ is defined by the $cref/argument c/det_of_minor/c/$$.
  $pre
  
  $$
  This template function
  is for example and testing purposes only.
  Expansion by minors is chosen as an example because it uses
  a lot of floating point operations yet does not require much source code
  (on the order of $icode m$$ factorial floating point operations and
  about 70 lines of source code including comments).
  This is not an efficient method for computing a determinant;
  for example, using an LU factorization would be better.
  
  $head Determinant of A$$
  If the following conditions hold, the minor is the
  entire matrix $latex A$$ and hence $code det_of_minor$$
  will return the determinant of $latex A$$:
  
  $list number$$
  $latex n = m$$.
  $lnext
  for $latex i = 0 , \ldots , m-1$$, $latex r[i] = i+1$$,
  and $latex r[m] = 0$$.
  $lnext
  for $latex j = 0 , \ldots , m-1$$, $latex c[j] = j+1$$,
  and $latex c[m] = 0$$.
  $lend
  
  $head a$$
  The argument $icode a$$ has prototype
  $codei%
  	const std::vector<%Scalar%>& %a%
  %$$
  and is a vector with size $latex m * m$$
  (see description of $cref/Scalar/det_of_minor/Scalar/$$ below).
  The elements of the $latex m \times m$$ matrix $latex A$$ are defined,
  for $latex i = 0 , \ldots , m-1$$ and $latex j = 0 , \ldots , m-1$$, by
  $latex \[
  	A_{i,j} = a[ i * m + j]
  \] $$
  
  $head m$$
  The argument $icode m$$ has prototype
  $codei%
  	size_t %m%
  %$$
  and is the number of rows (and columns) in the square matrix $latex A$$.
  
  $head n$$
  The argument $icode n$$ has prototype
  $codei%
  	size_t %n%
  %$$
  and is the number of rows (and columns) in the square minor $latex M$$.
  
  $head r$$
  The argument $icode r$$ has prototype
  $codei%
  	std::vector<size_t>& %r%
  %$$
  and is a vector with $latex m + 1$$ elements.
  This vector defines the function $latex R(i)$$
  which specifies the rows of the minor $latex M$$.
  To be specific, the function $latex R(i)$$
  for $latex i = 0, \ldots , n-1$$ is defined by
  $latex \[
  \begin{array}{rcl}
  	R(0)   & = & r[m]
  	\\
  	R(i+1) & = & r[ R(i) ]
  \end{array}
  \] $$
  All the elements of $icode r$$ must have value
  less than or equal $icode m$$.
  The elements of vector $icode r$$ are modified during the computation,
  and restored to their original value before the return from
  $code det_of_minor$$.
  
  $head c$$
  The argument $icode c$$ has prototype
  $codei%
  	std::vector<size_t>& %c%
  %$$
  and is a vector with $latex m + 1$$ elements
  This vector defines the function $latex C(i)$$
  which specifies the rows of the minor $latex M$$.
  To be specific, the function $latex C(i)$$
  for $latex j = 0, \ldots , n-1$$ is defined by
  $latex \[
  \begin{array}{rcl}
  	C(0)   & = & c[m]
  	\\
  	C(j+1) & = & c[ C(j) ]
  \end{array}
  \] $$
  All the elements of $icode c$$ must have value
  less than or equal $icode m$$.
  The elements of vector $icode c$$ are modified during the computation,
  and restored to their original value before the return from
  $code det_of_minor$$.
  
  $head d$$
  The result $icode d$$ has prototype
  $codei%
  	%Scalar% %d%
  %$$
  and is equal to the determinant of the minor $latex M$$.
  
  $head Scalar$$
  If $icode x$$ and $icode y$$ are objects of type $icode Scalar$$
  and $icode i$$ is an object of type $code int$$,
  the $icode Scalar$$ must support the following operations:
  $table
  $bold Syntax$$
  	$cnext $bold Description$$
  	$cnext $bold Result Type$$
  $rnext
  $icode%Scalar% %x%$$
  	$cnext default constructor for $icode Scalar$$ object.
  $rnext
  $icode%x% = %i%$$
  	$cnext set value of $icode x$$ to current value of $icode i$$
  $rnext
  $icode%x% = %y%$$
  	$cnext set value of $icode x$$ to current value of $icode y$$
  $rnext
  $icode%x% + %y%$$
  	$cnext value of $icode x$$ plus $icode y$$
  	$cnext $icode Scalar$$
  $rnext
  $icode%x% - %y%$$
  	$cnext value of $icode x$$ minus $icode y$$
  	$cnext $icode Scalar$$
  $rnext
  $icode%x% * %y%$$
  	$cnext value of $icode x$$ times value of $icode y$$
  	$cnext $icode Scalar$$
  $tend
  
  $children%
  	speed/example/det_of_minor.cpp%
  	omh/det_of_minor_hpp.omh
  %$$
  
  $head Example$$
  The file
  $cref det_of_minor.cpp$$
  contains an example and test of $code det_of_minor.hpp$$.
  It returns true if it succeeds and false otherwise.
  
  $head Source Code$$
  The file
  $cref det_of_minor.hpp$$
  contains the source for this template function.
  
  
  $end
  ---------------------------------------------------------------------------
  */
  // BEGIN C++
  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 representaion 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
  // END C++
  # endif
  