1 #ifndef vnl_algo_determinant_hxx_
2 #define vnl_algo_determinant_hxx_
3 /*
4 fsm
5 */
6 #include "vnl_determinant.h"
7
8 #include <cassert>
9 #ifdef _MSC_VER
10 # include <vcl_msvc_warnings.h>
11 #endif
12 #include <vnl/algo/vnl_qr.h>
13
14
15 template <class T>
vnl_determinant(T const * row0,T const * row1)16 T vnl_determinant(T const *row0, T const *row1)
17 {
18 return row0[0]*row1[1] - row0[1]*row1[0];
19 }
20
21 template <class T>
vnl_determinant(T const * row0,T const * row1,T const * row2)22 T vnl_determinant(T const *row0, T const *row1, T const *row2)
23 {
24 return // the extra '+' makes it work nicely with emacs indentation.
25 + row0[0]*row1[1]*row2[2]
26 - row0[0]*row2[1]*row1[2]
27 - row1[0]*row0[1]*row2[2]
28 + row1[0]*row2[1]*row0[2]
29 + row2[0]*row0[1]*row1[2]
30 - row2[0]*row1[1]*row0[2];
31 }
32
33 template <class T>
vnl_determinant(T const * row0,T const * row1,T const * row2,T const * row3)34 T vnl_determinant(T const *row0, T const *row1, T const *row2, T const *row3)
35 {
36 return
37 + row0[0]*row1[1]*row2[2]*row3[3]
38 - row0[0]*row1[1]*row3[2]*row2[3]
39 - row0[0]*row2[1]*row1[2]*row3[3]
40 + row0[0]*row2[1]*row3[2]*row1[3]
41 + row0[0]*row3[1]*row1[2]*row2[3]
42 - row0[0]*row3[1]*row2[2]*row1[3]
43 - row1[0]*row0[1]*row2[2]*row3[3]
44 + row1[0]*row0[1]*row3[2]*row2[3]
45 + row1[0]*row2[1]*row0[2]*row3[3]
46 - row1[0]*row2[1]*row3[2]*row0[3]
47 - row1[0]*row3[1]*row0[2]*row2[3]
48 + row1[0]*row3[1]*row2[2]*row0[3]
49 + row2[0]*row0[1]*row1[2]*row3[3]
50 - row2[0]*row0[1]*row3[2]*row1[3]
51 - row2[0]*row1[1]*row0[2]*row3[3]
52 + row2[0]*row1[1]*row3[2]*row0[3]
53 + row2[0]*row3[1]*row0[2]*row1[3]
54 - row2[0]*row3[1]*row1[2]*row0[3]
55 - row3[0]*row0[1]*row1[2]*row2[3]
56 + row3[0]*row0[1]*row2[2]*row1[3]
57 + row3[0]*row1[1]*row0[2]*row2[3]
58 - row3[0]*row1[1]*row2[2]*row0[3]
59 - row3[0]*row2[1]*row0[2]*row1[3]
60 + row3[0]*row2[1]*row1[2]*row0[3];
61 }
62
63 //--------------------------------------------------------------------------------
64
65 template <class T>
vnl_determinant(vnl_matrix<T> const & M,bool balance)66 T vnl_determinant(vnl_matrix<T> const &M, bool balance)
67 {
68 unsigned n = M.rows();
69 assert(M.cols() == n);
70
71 switch (n)
72 {
73 case 1: return M[0][0];
74 case 2: return vnl_determinant(M[0], M[1]);
75 case 3: return vnl_determinant(M[0], M[1], M[2]);
76 case 4: return vnl_determinant(M[0], M[1], M[2], M[3]);
77 default:
78 if (balance)
79 {
80 vnl_matrix<T> tmp(M);
81 typedef typename vnl_numeric_traits<T>::abs_t abs_t;
82 abs_t scalings(1);
83 for (int t=0; t<5; ++t)
84 {
85 // normalize rows.
86 for (unsigned int i=0; i<n; ++i) {
87 abs_t rn = tmp.get_row(i).rms();
88 if (rn > 0) {
89 scalings *= rn;
90 tmp.scale_row(i, abs_t(1)/rn);
91 }
92 }
93 // normalize columns.
94 for (unsigned int i=0; i<n; ++i) {
95 abs_t rn = tmp.get_column(i).rms();
96 if (rn > 0) {
97 scalings *= rn;
98 tmp.scale_column(i, abs_t(1)/rn);
99 }
100 }
101 }
102 T balanced_det = vnl_qr<T>(tmp).determinant();
103 //std::clog << __FILE__ ": scalings, balanced_det = " << scalings << ", " << balanced_det << std::endl;
104 return T(scalings) * balanced_det;
105 }
106 else
107 return vnl_qr<T>(M).determinant();
108 }
109 }
110
111
112 //--------------------------------------------------------------------------------
113
114 #define VNL_DETERMINANT_INSTANTIATE_1(T) \
115 template VNL_ALGO_EXPORT T vnl_determinant(T const *, T const *); \
116 template VNL_ALGO_EXPORT T vnl_determinant(T const *, T const *, T const *); \
117 template VNL_ALGO_EXPORT T vnl_determinant(T const *, T const *, T const *, T const *)
118
119 #define VNL_DETERMINANT_INSTANTIATE_2(T) \
120 template VNL_ALGO_EXPORT T vnl_determinant(vnl_matrix<T > const &, bool)
121
122 #undef VNL_DETERMINANT_INSTANTIATE
123 #define VNL_DETERMINANT_INSTANTIATE(T) \
124 VNL_DETERMINANT_INSTANTIATE_1(T); \
125 VNL_DETERMINANT_INSTANTIATE_2(T)
126
127 #endif // vnl_algo_determinant_hxx_
128