Attachment #355 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
{

  for (Index i=1; i < m_A.cols(); ++i) {
    res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
    if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) {
      res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1);
    }
    else if (2*abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) || 2*abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1))) {
      res.coeffRef(i-1,i) = (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1));
    }
    else {
      int unwindingNumber = std::ceil((numext::imag(logTdiag[i]-logTdiag[i-1]) - M_PI) / (2*M_PI));
      Scalar w = internal::matrix_power_unwinder<Scalar>::run(m_A.coeff(i,i), m_A.coeff(i-1,i-1), unwindingNumber);
      res.coeffRef(i-1,i) = RealScalar(2) * std::exp(RealScalar(0.5)*p*(logTdiag[i]+logTdiag[i-1])) * std::sinh(p * w)
          / (m_A.coeff(i,i) - m_A.coeff(i-1,i-1));
    }
    res.coeffRef(i-1,i) *= m_A.coeff(i-1,i);
  }
}


  int degree, degree2, numberOfSquareRoots = 0;
  bool hasExtraSquareRoot = false;

  /* FIXME
   * For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite
   * loop.  We should move 0 eigenvalues to bottom right corner.  We need not
   * worry about tiny values (e.g. 1e-300) because they will reach 1 if
   * repetitively sqrt'ed.
   *
   * [ T  A ]^p   [ T^p  (T^-1 T^p A) ]
   * [      ]   = [                   ]
   * [ 0  0 ]     [  0         0      ]
   */
  for (Index i=0; i < m_A.cols(); ++i)
    eigen_assert(m_A(i,i) != RealScalar(0))

  while (true) {
    IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
    normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();


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 (-4 / +20 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 332-348 Link Here
332	for (Index i=1; i < m_A.cols(); ++i) {	334	for (Index i=1; i < m_A.cols(); ++i) {
333	res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);	335	res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
334	if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) {	336	if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) {
335	res.coeffRef(i-1,i) = p * pow(m_A.coeff(i-1,i), p-1);	337	res.coeffRef(i-1,i) = p * pow(m_A.coeff(i,i), p-1);
336	}	338	}
337	else if (2abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) \|\| 2abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1))) {	339	else if (2abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) \|\| 2abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1))) {
338	res.coeffRef(i-1,i) = m_A.coeff(i-1,i) * (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1));	340	res.coeffRef(i-1,i) = (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.coeff(i-1,i-1));
339	}	341	}
340	else {	342	else {
341	int unwindingNumber = std::ceil((numext::imag(logTdiag[i]-logTdiag[i-1]) - M_PI) / (2*M_PI));	343	int unwindingNumber = std::ceil((numext::imag(logTdiag[i]-logTdiag[i-1]) - M_PI) / (2*M_PI));
342	Scalar w = internal::matrix_power_unwinder<Scalar>::run(m_A.coeff(i,i), m_A.coeff(i-1,i-1), unwindingNumber);	344	Scalar w = internal::matrix_power_unwinder<Scalar>::run(m_A.coeff(i,i), m_A.coeff(i-1,i-1), unwindingNumber);
343	res.coeffRef(i-1,i) = m_A.coeff(i-1,i) * RealScalar(2) * std::exp(RealScalar(0.5)p(logTdiag[i]+logTdiag[i-1])) *	345	res.coeffRef(i-1,i) = RealScalar(2) * std::exp(RealScalar(0.5)p(logTdiag[i]+logTdiag[i-1])) * std::sinh(p * w)
344	std::sinh(p * w) / (m_A.coeff(i,i) - m_A.coeff(i-1,i-1));	346	/ (m_A.coeff(i,i) - m_A.coeff(i-1,i-1));
345	}	347	}
		348	res.coeffRef(i-1,i) *= m_A.coeff(i-1,i);
346	}	349	}
347	}	350	}
348		351
Lines 360-365 Link Here
360	int degree, degree2, numberOfSquareRoots = 0;	363	int degree, degree2, numberOfSquareRoots = 0;
361	bool hasExtraSquareRoot = false;	364	bool hasExtraSquareRoot = false;
362		365
		366	/* FIXME
		367	* For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite
		368	* loop. We should move 0 eigenvalues to bottom right corner. We need not
		369	* worry about tiny values (e.g. 1e-300) because they will reach 1 if
		370	* repetitively sqrt'ed.
		371	*
		372	* [ T A ]^p [ T^p (T^-1 T^p A) ]
		373	* [ ] = [ ]
		374	* [ 0 0 ] [ 0 0 ]
		375	*/
		376	for (Index i=0; i < m_A.cols(); ++i)
		377	eigen_assert(m_A(i,i) != RealScalar(0))
		378
363	while (true) {	379	while (true) {
364	IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;	380	IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
365	normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();	381	normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
366		382