Bugfix in MatrixLogarithm.h

This commit is contained in:
jdh8 2012-08-18 21:28:05 +08:00
commit 15dabd4db7
5 changed files with 32 additions and 29 deletions

View File

@ -490,7 +490,7 @@ Read-write access to sub-vectors:
<tr><td>\code vec1.head(n)\endcode</td><td>\code vec1.head<n>()\endcode</td><td>the first \c n coeffs </td></tr> <tr><td>\code vec1.head(n)\endcode</td><td>\code vec1.head<n>()\endcode</td><td>the first \c n coeffs </td></tr>
<tr><td>\code vec1.tail(n)\endcode</td><td>\code vec1.tail<n>()\endcode</td><td>the last \c n coeffs </td></tr> <tr><td>\code vec1.tail(n)\endcode</td><td>\code vec1.tail<n>()\endcode</td><td>the last \c n coeffs </td></tr>
<tr><td>\code vec1.segment(pos,n)\endcode</td><td>\code vec1.segment<n>(pos)\endcode</td> <tr><td>\code vec1.segment(pos,n)\endcode</td><td>\code vec1.segment<n>(pos)\endcode</td>
<td>the \c n coeffs in \n the range [\c pos : \c pos + \c n [</td></tr> <td>the \c n coeffs in the \n range [\c pos : \c pos + \c n - 1]</td></tr>
<tr class="alt"><td colspan="3"> <tr class="alt"><td colspan="3">
Read-write access to sub-matrices:</td></tr> Read-write access to sub-matrices:</td></tr>

View File

@ -6,10 +6,10 @@ using namespace std;
int main() int main()
{ {
MatrixXf m = MatrixXf::Random(3,3); MatrixXd m = MatrixXd::Random(3,3);
m = (m + MatrixXf::Constant(3,3,1.2)) * 50; m = (m + MatrixXd::Constant(3,3,1.2)) * 50;
cout << "m =" << endl << m << endl; cout << "m =" << endl << m << endl;
VectorXf v(3); VectorXd v(3);
v << 1, 2, 3; v << 1, 2, 3;
cout << "m * v =" << endl << m * v << endl; cout << "m * v =" << endl << m * v << endl;
} }

View File

@ -6,10 +6,10 @@ using namespace std;
int main() int main()
{ {
Matrix3f m = Matrix3f::Random(); Matrix3d m = Matrix3d::Random();
m = (m + Matrix3f::Constant(1.2)) * 50; m = (m + Matrix3d::Constant(1.2)) * 50;
cout << "m =" << endl << m << endl; cout << "m =" << endl << m << endl;
Vector3f v(1,2,3); Vector3d v(1,2,3);
cout << "m * v =" << endl << m * v << endl; cout << "m * v =" << endl << m * v << endl;
} }

View File

@ -173,10 +173,11 @@ int MatrixLogarithmAtomic<MatrixType>::getPadeDegree(float normTminusI)
{ {
const float maxNormForPade[] = { 2.5111573934555054e-1 /* degree = 3 */ , 4.0535837411880493e-1, const float maxNormForPade[] = { 2.5111573934555054e-1 /* degree = 3 */ , 4.0535837411880493e-1,
5.3149729967117310e-1 }; 5.3149729967117310e-1 };
for (int degree = 3; degree <= maxPadeDegree; ++degree) int degree = 3;
for (; degree <= maxPadeDegree; ++degree)
if (normTminusI <= maxNormForPade[degree - minPadeDegree]) if (normTminusI <= maxNormForPade[degree - minPadeDegree])
break;
return degree; return degree;
assert(false); // this line should never be reached
} }
/* \brief Get suitable degree for Pade approximation. (specialized for RealScalar = double) */ /* \brief Get suitable degree for Pade approximation. (specialized for RealScalar = double) */
@ -185,10 +186,11 @@ int MatrixLogarithmAtomic<MatrixType>::getPadeDegree(double normTminusI)
{ {
const double maxNormForPade[] = { 1.6206284795015624e-2 /* degree = 3 */ , 5.3873532631381171e-2, const double maxNormForPade[] = { 1.6206284795015624e-2 /* degree = 3 */ , 5.3873532631381171e-2,
1.1352802267628681e-1, 1.8662860613541288e-1, 2.642960831111435e-1 }; 1.1352802267628681e-1, 1.8662860613541288e-1, 2.642960831111435e-1 };
for (int degree = 3; degree <= maxPadeDegree; ++degree) int degree = 3;
for (; degree <= maxPadeDegree; ++degree)
if (normTminusI <= maxNormForPade[degree - minPadeDegree]) if (normTminusI <= maxNormForPade[degree - minPadeDegree])
break;
return degree; return degree;
assert(false); // this line should never be reached
} }
/* \brief Get suitable degree for Pade approximation. (specialized for RealScalar = long double) */ /* \brief Get suitable degree for Pade approximation. (specialized for RealScalar = long double) */
@ -215,10 +217,11 @@ int MatrixLogarithmAtomic<MatrixType>::getPadeDegree(long double normTminusI)
3.6688019729653446926585242192447447e-2L, 5.9290962294020186998954055264528393e-2L, 3.6688019729653446926585242192447447e-2L, 5.9290962294020186998954055264528393e-2L,
8.6998436081634343903250580992127677e-2L, 1.1880960220216759245467951592883642e-1L }; 8.6998436081634343903250580992127677e-2L, 1.1880960220216759245467951592883642e-1L };
#endif #endif
for (int degree = 3; degree <= maxPadeDegree; ++degree) int degree = 3
for (; degree <= maxPadeDegree; ++degree)
if (normTminusI <= maxNormForPade[degree - minPadeDegree]) if (normTminusI <= maxNormForPade[degree - minPadeDegree])
break;
return degree; return degree;
assert(false); // this line should never be reached
} }
/* \brief Compute Pade approximation to matrix logarithm */ /* \brief Compute Pade approximation to matrix logarithm */
@ -424,7 +427,7 @@ void MatrixLogarithmAtomic<MatrixType>::computePade11(MatrixType& result, const
* This class holds the argument to the matrix function until it is * This class holds the argument to the matrix function until it is
* assigned or evaluated for some other reason (so the argument * assigned or evaluated for some other reason (so the argument
* should not be changed in the meantime). It is the return type of * should not be changed in the meantime). It is the return type of
* matrixBase::log() and most of the time this is the only way it * MatrixBase::log() and most of the time this is the only way it
* is used. * is used.
*/ */
template<typename Derived> class MatrixLogarithmReturnValue template<typename Derived> class MatrixLogarithmReturnValue

View File

@ -462,11 +462,11 @@ template <typename MatrixType, typename RealScalar, typename PlainObject, int Is
inline int MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::getPadeDegree(float normIminusT) inline int MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::getPadeDegree(float normIminusT)
{ {
const float maxNormForPade[] = { 2.7996156e-1f /* degree = 3 */ , 4.3268868e-1f }; const float maxNormForPade[] = { 2.7996156e-1f /* degree = 3 */ , 4.3268868e-1f };
int degree = 3;
for (int degree = 3; degree <= 4; degree++) for (; degree <= 4; degree++)
if (normIminusT <= maxNormForPade[degree - 3]) if (normIminusT <= maxNormForPade[degree - 3])
break;
return degree; return degree;
assert(false); // this line should never be reached
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger>
@ -474,11 +474,11 @@ inline int MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::getPadeDegr
{ {
const double maxNormForPade[] = { 1.882832775783710e-2 /* degree = 3 */ , 6.036100693089536e-2, const double maxNormForPade[] = { 1.882832775783710e-2 /* degree = 3 */ , 6.036100693089536e-2,
1.239372725584857e-1, 1.998030690604104e-1, 2.787629930861592e-1 }; 1.239372725584857e-1, 1.998030690604104e-1, 2.787629930861592e-1 };
int degree = 3;
for (int degree = 3; degree <= 7; degree++) for (; degree <= 7; degree++)
if (normIminusT <= maxNormForPade[degree - 3]) if (normIminusT <= maxNormForPade[degree - 3])
break;
return degree; return degree;
assert(false); // this line should never be reached
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger>
@ -508,11 +508,11 @@ inline int MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::getPadeDegr
3.907876732697568523164749432441966e-2L, 6.266303975524852476985111609267074e-2L, 3.907876732697568523164749432441966e-2L, 6.266303975524852476985111609267074e-2L,
9.133823549851655878933476070874651e-2L }; 9.133823549851655878933476070874651e-2L };
#endif #endif
int degree = 3;
for (int degree = 3; degree <= maxPadeDegree; degree++) for (; degree <= maxPadeDegree; degree++)
if (normIminusT <= maxNormForPade[degree - 3]) if (normIminusT <= maxNormForPade[degree - 3])
break;
return degree; return degree;
assert(false); // this line should never be reached
} }
template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger> template <typename MatrixType, typename RealScalar, typename PlainObject, int IsInteger>
void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computePade(const int& degree, const ComplexMatrix& IminusT) void MatrixPower<MatrixType,RealScalar,PlainObject,IsInteger>::computePade(const int& degree, const ComplexMatrix& IminusT)