00001 #ifndef vnl_adjugate_txx_ 00002 #define vnl_adjugate_txx_ 00003 //: 00004 // \file 00005 // \author fsm 00006 00007 #include "vnl_adjugate.h" 00008 #include <vnl/vnl_matrix.h> 00009 #include <vnl/algo/vnl_determinant.h> 00010 00011 // This is a rudimentary implementation. It could be improved by noting 00012 // that adj(A B) = adj(B) adj(A) for all matrices A, B (invertible or 00013 // not) and then using a matrix decomposition for larger matrices. 00014 // 00015 // E.g. using a singular value decomposition A = U D V^* gives 00016 // adj(A) = V adj(D) U^*. 00017 // 00018 // On the other hand, SVD decomposition makes no sense for e.g. integer matrices 00019 // and we want to keep T as general as possible. 00020 00021 template <class T> 00022 void vnl_adjugate(vnl_matrix<T> const &A, vnl_matrix<T> *out) 00023 { 00024 int n = A.rows(); 00025 A.assert_size(n, n); 00026 out->assert_size(n, n); 00027 00028 vnl_matrix<T> sub(n-1, n-1); 00029 for (int i=0; i<n; ++i) 00030 for (int j=0; j<n; ++j) 00031 { 00032 for (int u=0; u<n-1; ++u) 00033 for (int v=0; v<n-1; ++v) 00034 sub[u][v] = A[v+(v<i?0:1)][u+(u<j?0:1)]; 00035 (*out)[i][j] = vnl_determinant(sub, false); 00036 } 00037 } 00038 00039 template <class T> 00040 vnl_matrix<T> vnl_adjugate(vnl_matrix<T> const &A) 00041 { 00042 vnl_matrix<T> adj(A.rows(), A.cols()); 00043 vnl_adjugate(A, &adj); 00044 return adj; 00045 } 00046 00047 //-------------------------------------------------------------------------------- 00048 00049 #undef VNL_ADJUGATE_INSTANTIATE 00050 #define VNL_ADJUGATE_INSTANTIATE(T) \ 00051 template void vnl_adjugate(vnl_matrix<T > const &, vnl_matrix<T > *); \ 00052 template vnl_matrix<T > vnl_adjugate(vnl_matrix<T > const &) 00053 00054 #endif