core/vnl/vnl_power.h
Go to the documentation of this file.
00001 // This is core/vnl/vnl_power.h
00002 #ifndef vnl_power_h_
00003 #define vnl_power_h_
00004 //:
00005 // \file
00006 // \brief Calculates nth power of a small vnl_matrix_fixed (not using svd)
00007 // \author Peter Vanroose
00008 // \date   21 July 2009
00009 //
00010 // \verbatim
00011 //  Modifications
00012 //   <none yet>
00013 // \endverbatim
00014 
00015 #include <vnl/vnl_matrix_fixed.h>
00016 #include <vnl/vnl_matrix.h>
00017 #include <vnl/vnl_inverse.h> // used for negative powers
00018 #include <vcl_cassert.h>
00019 
00020 //: Calculates nth power of a vnl_matrix_fixed (not using svd)
00021 //  This allows you to write e.g.
00022 //
00023 //  x = vnl_power(A,7) * vnl_power(B,-4) * b;
00024 //
00025 // Note that this function is inlined (except for the call to vnl_inverse()),
00026 // which makes it much faster than a full-fledged square matrix power
00027 // implementation using svd, which belongs in vnl/algo.
00028 //
00029 //  \relatesalso vnl_matrix_fixed
00030 
00031 template <class T, unsigned int d>
00032 vnl_matrix_fixed<T,d,d> vnl_power(vnl_matrix_fixed<T,d,d> const& m, int n)
00033 {
00034   assert(n >= 0 || d <= 4); // to allow the use of vnl_inverse()
00035   if (n == 0)
00036     return vnl_matrix_fixed<T,d,d>().set_identity();
00037   else if (n == 1 || m.is_identity())
00038     return m;
00039   else if (n < 0)
00040     return vnl_inverse(vnl_power(m, -n));
00041   else {
00042     vnl_matrix_fixed<T,d,d> r = vnl_power(m, n/2);
00043     return n%2==0 ? r * r : r * r * m;
00044   }
00045 }
00046 
00047 //: Calculates nth power of a square vnl_matrix (not using svd)
00048 //  This allows you to write e.g.
00049 //
00050 //  x = vnl_power(A,7) * vnl_power(B,-4) * b;
00051 //
00052 // Note that this function is inlined (except for the call to vnl_inverse()),
00053 // which makes it much faster than a full-fledged square matrix power
00054 // implementation using svd, which belongs in vnl/algo.
00055 //
00056 //  \relatesalso vnl_matrix
00057 
00058 template <class T>
00059 vnl_matrix<T> vnl_power(vnl_matrix<T> const& m, int n)
00060 {
00061   assert(m.rows() == m.columns());
00062   assert(n >= 0 || m.rows() <= 4);
00063   if (n == 0)
00064     return vnl_matrix<T>(m.rows(),m.columns()).set_identity();
00065   else if (n == 1 || m.is_identity())
00066     return m;
00067   else if (n < 0 && m.rows() == 1)
00068     return vnl_power(vnl_matrix_fixed<T,1,1>(m),n).as_ref();
00069   else if (n < 0 && m.rows() == 2)
00070     return vnl_power(vnl_matrix_fixed<T,2,2>(m),n).as_ref();
00071   else if (n < 0 && m.rows() == 3)
00072     return vnl_power(vnl_matrix_fixed<T,3,3>(m),n).as_ref();
00073   else if (n < 0 && m.rows() == 4)
00074     return vnl_power(vnl_matrix_fixed<T,4,4>(m),n).as_ref();
00075   else {
00076     vnl_matrix<T> r = vnl_power(m, n/2);
00077     return n%2==0 ? r * r : r * r * m;
00078   }
00079 }
00080 
00081 #endif // vnl_power_h_