core/vnl/algo/vnl_adjugate.txx
Go to the documentation of this file.
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