det_of_minor_c

View page source

Determinant of a Minor

Syntax

d = det_of_minor ( a , m , n , r , c )

Purpose

returns the determinant of a minor of the matrix \(A\) using expansion by minors. The elements of the \(n \times n\) minor \(M\) of the matrix \(A\) are defined, for \(i = 0 , \ldots , n-1\) and \(j = 0 , \ldots , n-1\), by

\[M_{i,j} = A_{R(i), C(j)}\]

where the functions \(R(i)\) is defined by the argument r and \(C(j)\) is defined by the argument c .

This 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 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.

Determinant of A

If the following conditions hold, the minor is the entire matrix \(A\) and hence det_of_minor will return the determinant of \(A\):

  1. \(n = m\).

  2. for \(i = 0 , \ldots , m-1\), \(r[i] = i+1\), and \(r[m] = 0\).

  3. for \(j = 0 , \ldots , m-1\), \(c[j] = j+1\), and \(c[m] = 0\).

a

The argument a has prototype

const double * a

and is a vector with size \(m * m\). The elements of the \(m \times m\) matrix \(A\) are defined, for \(i = 0 , \ldots , m-1\) and \(j = 0 , \ldots , m-1\), by

\[A_{i,j} = a[ i * m + j]\]

m

The argument m has prototype

size_t m

and is the size of the square matrix \(A\).

n

The argument n has prototype

size_t n

and is the size of the square minor \(M\).

r

The argument r has prototype

size_t * r

and is a vector with \(m + 1\) elements. This vector defines the function \(R(i)\) which specifies the rows of the minor \(M\). To be specific, the function \(R(i)\) for \(i = 0, \ldots , n-1\) is defined by

\begin{eqnarray} R(0) & = & r[m] \\ R(i+1) & = & r[ R(i) ] \end{eqnarray}

All the elements of r must have value less than or equal m . The elements of vector r are modified during the computation, and restored to their original value before the return from det_of_minor .

c

The argument c has prototype

size_t * c

and is a vector with \(m + 1\) elements This vector defines the function \(C(i)\) which specifies the rows of the minor \(M\). To be specific, the function \(C(i)\) for \(j = 0, \ldots , n-1\) is defined by

\begin{eqnarray} C(0) & = & c[m] \\ C(j+1) & = & c[ C(j) ] \end{eqnarray}

All the elements of c must have value less than or equal m . The elements of vector c are modified during the computation, and restored to their original value before the return from det_of_minor .

d

The result d has prototype

double d

and is equal to the determinant of the minor \(M\).

Source Code

double det_of_minor(
   const double*        a  ,
   size_t               m  ,
   size_t               n  ,
   size_t*              r  ,
   size_t*              c  )
{  size_t R0, Cj, Cj1, j;
   double detM, M0j, detS;
   int s;

   R0 = r[m]; /* R(0) */
   Cj = c[m]; /* C(j)    (case j = 0) */
   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 */
   detM = 0.;

   /* initialize sign of factor for neat sub-minor */
   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(j = 0; j < n; j++)
   {  /* element with index (0,j) in the minor M */
      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 */
      detS = det_of_minor(a, m, n - 1, r, c);

      /* restore column Cj to representation 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 neat 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;
}