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}} }@)@

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

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.
If the following conditions hold, the minor is the entire matrix @(@ A @)@ and hence

`det_of_minor`

will return the determinant of @(@
A
@)@:
- @(@ n = m @)@.
- for @(@ i = 0 , \ldots , m-1 @)@, @(@ r[i] = i+1 @)@, and @(@ r[m] = 0 @)@.
- for @(@ j = 0 , \ldots , m-1 @)@, @(@ c[j] = j+1 @)@, and @(@ c[m] = 0 @)@.

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]
@]@
The argument

*m*

has prototype

size_t *m*

and is the size of the square matrix @(@
A
@)@.
The argument

*n*

has prototype

size_t *n*

and is the size of the square minor @(@
M
@)@.
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`

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

.
The result

*d*

has prototype

double *d*

and is equal to the determinant of the minor @(@
M
@)@.
doubledet_of_minor(constdouble* 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 )returna[ 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;elsedetM = 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 */returndetM; }

Input File: test_more/compare_c/det_by_minor.c