Prev Next Index-> contents reference index search external Up-> CppAD Appendix compare_c det_of_minor_c CppAD-> Install Introduction AD ADFun preprocessor multi_thread utility ipopt_solve Example speed Appendix Appendix-> Faq directory Theory glossary Bib wish_list whats_new deprecated compare_c numeric_ad addon License compare_c-> det_of_minor_c det_by_minor_c uniform_01_c correct_det_by_minor_c repeat_det_by_minor_c elapsed_seconds_c time_det_by_minor_c main_compare_c det_of_minor_c Headings-> Syntax Purpose Determinant of A a m n r c d Source Code

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