core/vnl/vnl_matrix_exp.txx
Go to the documentation of this file.
00001 // This is core/vnl/vnl_matrix_exp.txx
00002 #ifndef vnl_matrix_exp_txx_
00003 #define vnl_matrix_exp_txx_
00004 //:
00005 // \file
00006 // \author fsm
00007 #include "vnl_matrix_exp.h"
00008 #include <vcl_cassert.h>
00009 #ifdef DEBUG
00010 #include <vcl_iostream.h>
00011 #endif
00012 
00013 template <class Matrix>
00014 bool vnl_matrix_exp(Matrix const &X, Matrix &expX, double max_err)
00015 {
00016   assert(X.rows() == X.cols());
00017   assert(X.rows() == expX.rows());
00018   assert(X.cols() == expX.cols());
00019   assert(max_err > 0);
00020 
00021   double norm_X = X.operator_inf_norm();
00022 #ifdef DEBUG
00023   vcl_cerr << "norm_X = " << norm_X << vcl_endl;
00024 #endif
00025 
00026   // exponential series
00027   expX.set_identity();
00028   Matrix acc(X);
00029   double norm_acc_bound = norm_X;
00030   for (unsigned n=1; true; ++n) {
00031     expX += acc;
00032 #ifdef DEBUG
00033     vcl_cerr << "n=" << n << vcl_endl;
00034 #endif
00035 
00036     if (norm_X < n) {
00037       double err_bound = norm_acc_bound / (1 - norm_X/n);
00038 #ifdef DEBUG
00039       vcl_cerr << "err_bound = " << err_bound << vcl_endl;
00040 #endif
00041       if (err_bound < max_err)
00042         break;
00043     }
00044 
00045     acc = acc * X;
00046     acc /= n+1;
00047 
00048     norm_acc_bound *= norm_X/(n+1);
00049   }
00050 
00051   return true;
00052 }
00053 
00054 
00055 template <class Matrix>
00056 Matrix vnl_matrix_exp(Matrix const &X)
00057 {
00058   Matrix expX(X.rows(), X.cols());
00059 #ifndef NDEBUG
00060   bool retval =
00061 #endif
00062   vnl_matrix_exp(X, expX, 1e-10);
00063 
00064   assert(retval);
00065   return expX;
00066 }
00067 
00068 //------------------------------------------------------------------------------
00069 
00070 #undef VNL_MATRIX_EXP_INSTANTIATE
00071 #define VNL_MATRIX_EXP_INSTANTIATE(Matrix) \
00072 template bool vnl_matrix_exp(Matrix const&, Matrix &, double); \
00073 template Matrix vnl_matrix_exp(Matrix const&)
00074 
00075 #endif // vnl_matrix_exp_txx_