Attachment #352 for bug #614

This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen




template<typename MatrixType>
void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const
{
  res.resizeLike(m_A);
  switch (m_A.rows()) {
    case 0:
      break;

  res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
}

// this function assumes that res has the correct size (see bug 614)
template<typename MatrixType>
void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
{

						   9.134603732914548552537150753385375e-2L; // quadruple precision
  MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
  RealScalar normIminusT;
  int degree, degree2, numberOfSquareRoots;
  bool hasExtraSquareRoot = false;

  // Sqrt at most 64 times to avoid infinite loop for ill conditioned (I - T)
  for (numberOfSquareRoots = 0; numberOfSquareRoots < 64; ++numberOfSquareRoots) {
    IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
    normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
    if (normIminusT < maxNormForPade) {

    }
    MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
    T = sqrtT.template triangularView<Upper>();
    ++numberOfSquareRoots;
  }
  computePade(degree, IminusT, res, p);


Lines 181-186 Link Here

(-)a/unsupported/test/matrix_power.cpp (-1 / +1 lines)
181	CALL_SUBTEST_1(testMatrixVector(Matrix2f(), Vector2f(), 1e-4));	181	CALL_SUBTEST_1(testMatrixVector(Matrix2f(), Vector2f(), 1e-4));
182	CALL_SUBTEST_5(testMatrixVector(Matrix3cf(), Vector3cf(), 1e-4));	182	CALL_SUBTEST_5(testMatrixVector(Matrix3cf(), Vector3cf(), 1e-4));
183	CALL_SUBTEST_8(testMatrixVector(Matrix4f(), Vector4f(), 1e-4));	183	CALL_SUBTEST_8(testMatrixVector(Matrix4f(), Vector4f(), 1e-4));
184	CALL_SUBTEST_6(testMatrixVector(MatrixXf(8,8), VectorXf(8), 1e-3));	184	CALL_SUBTEST_6(testMatrixVector(MatrixXf(2,2), VectorXf(2), 1e-3)); // see bug 614
185	CALL_SUBTEST_9(testMatrixVector(MatrixXe(7,7), VectorXe(7), 1e-13));	185	CALL_SUBTEST_9(testMatrixVector(MatrixXe(7,7), VectorXe(7), 1e-13));
186	}	186	}

Return to bug 614

Lines 294-299 Link Here

(-)a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h (-3 / +5 lines)
294	template<typename MatrixType>	294	template<typename MatrixType>
295	void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const	295	void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const
296	{	296	{
		297	res.resizeLike(m_A);
297	switch (m_A.rows()) {	298	switch (m_A.rows()) {
298	case 0:	299	case 0:
299	break;	300	break;
Lines 320-325 Link Here
320	res += MatrixType::Identity(IminusT.rows(), IminusT.cols());	321	res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
321	}	322	}
322		323
		324	// this function assumes that res has the correct size (see bug 614)
323	template<typename MatrixType>	325	template<typename MatrixType>
324	void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const	326	void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealScalar p) const
325	{	327	{
Lines 357-366 Link Here
357	9.134603732914548552537150753385375e-2L; // quadruple precision	359	9.134603732914548552537150753385375e-2L; // quadruple precision
358	MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();	360	MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>();
359	RealScalar normIminusT;	361	RealScalar normIminusT;
360	int degree, degree2, numberOfSquareRoots = 0;	362	int degree, degree2, numberOfSquareRoots;
361	bool hasExtraSquareRoot = false;	363	bool hasExtraSquareRoot = false;
362		364
363	while (true) {	365	// Sqrt at most 64 times to avoid infinite loop for ill conditioned (I - T)
		366	for (numberOfSquareRoots = 0; numberOfSquareRoots < 64; ++numberOfSquareRoots) {
364	IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;	367	IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
365	normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();	368	normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
366	if (normIminusT < maxNormForPade) {	369	if (normIminusT < maxNormForPade) {
Lines 372-378 Link Here
372	}	375	}
373	MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);	376	MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT);
374	T = sqrtT.template triangularView<Upper>();	377	T = sqrtT.template triangularView<Upper>();
375	++numberOfSquareRoots;
376	}	378	}
377	computePade(degree, IminusT, res, p);	379	computePade(degree, IminusT, res, p);
378		380