1 #ifndef vnl_adjugate_hxx_
2 #define vnl_adjugate_hxx_
3 //:
4 // \file
5 // \author fsm
6
7 #include "vnl_adjugate.h"
8 #include <vnl/vnl_matrix.h>
9 #include <vnl/algo/vnl_determinant.h>
10
11 // This is a rudimentary implementation. It could be improved by noting
12 // that adj(A B) = adj(B) adj(A) for all matrices A, B (invertible or
13 // not) and then using a matrix decomposition for larger matrices.
14 //
15 // E.g. using a singular value decomposition A = U D V^* gives
16 // adj(A) = V adj(D) U^*.
17 //
18 // On the other hand, SVD decomposition makes no sense for e.g. integer matrices
19 // and we want to keep T as general as possible.
20
21 template <class T>
vnl_adjugate(vnl_matrix<T> const & A,vnl_matrix<T> * out)22 void vnl_adjugate(vnl_matrix<T> const &A, vnl_matrix<T> *out)
23 {
24 int n = A.rows();
25 A.assert_size(n, n);
26 out->assert_size(n, n);
27
28 vnl_matrix<T> sub(n-1, n-1);
29 for (int i=0; i<n; ++i)
30 for (int j=0; j<n; ++j)
31 {
32 for (int u=0; u<n-1; ++u)
33 for (int v=0; v<n-1; ++v)
34 sub[u][v] = A[v+(v<i?0:1)][u+(u<j?0:1)];
35 (*out)[i][j] = vnl_determinant(sub, false);
36 }
37 }
38
39 template <class T>
vnl_adjugate(vnl_matrix<T> const & A)40 vnl_matrix<T> vnl_adjugate(vnl_matrix<T> const &A)
41 {
42 vnl_matrix<T> adj(A.rows(), A.cols());
43 vnl_adjugate(A, &adj);
44 return adj;
45 }
46
47 //--------------------------------------------------------------------------------
48
49 #undef VNL_ADJUGATE_INSTANTIATE
50 #define VNL_ADJUGATE_INSTANTIATE(T) \
51 template VNL_ALGO_EXPORT void vnl_adjugate(vnl_matrix<T > const &, vnl_matrix<T > *); \
52 template VNL_ALGO_EXPORT vnl_matrix<T > vnl_adjugate(vnl_matrix<T > const &)
53
54 #endif
55