Prev Next det_of_minor_c

@(@\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}} }@)@
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{array}{rcl} R(0) & = & r[m] \\ R(i+1) & = & r[ R(i) ] \end{array} @]@ 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{array}{rcl} C(0) & = & c[m] \\ C(j+1) & = & c[ C(j) ] \end{array} @]@ 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;
}
Input File: test_more/compare_c/det_by_minor.c