Go to the documentation of this file.00001
00002 #ifndef vnl_matrix_exp_txx_
00003 #define vnl_matrix_exp_txx_
00004
00005
00006
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
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_