Created attachment 347 [details] Matrix power segfaults for trivial 2x2 matrix. The attached file demonstrates pow() segfaulting on trivial matrix powers. I am using g++ (Gentoo 4.7.3 p1.0, pie-0.5.5) 4.7.3. Not all powers segfault, for example A.pow(1.0) does not. When I run the debugger (after compiling with -ggdb3), and step through until it segfaults, I get the following backtrace: ~/tmp/c++/eigen$ gdb ./matrix_power GNU gdb (Gentoo 7.6 p1) 7.6 Copyright (C) 2013 Free Software Foundation, Inc. License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html> This is free software: you are free to change and redistribute it. There is NO WARRANTY, to the extent permitted by law. Type "show copying" and "show warranty" for details. This GDB was configured as "x86_64-pc-linux-gnu". For bug reporting instructions, please see: <http://bugs.gentoo.org/>... Reading symbols from /home/hazelnusse/tmp/c++/eigen/matrix_power...done. (gdb) (gdb) list 1 // Compiled with: 2 // 3 // $ g++ -Wall matrix_power.cc -O2 -I/home/hazelnusse/usr/include/eigen3 -o matrix_power 4 // 5 #include <iostream> 6 #include <Eigen/Dense> 7 #include <unsupported/Eigen/MatrixFunctions> 8 int main() 9 { 10 Eigen::MatrixXcd A(2,2); (gdb) b main Breakpoint 1 at 0x401cd0: file matrix_power.cc, line 9. (gdb) run Starting program: /home/hazelnusse/tmp/c++/eigen/matrix_power warning: Could not load shared library symbols for linux-vdso.so.1. Do you need "set solib-search-path" or "set sysroot"? Breakpoint 1, main () at matrix_power.cc:9 9 { (gdb) n 10 Eigen::MatrixXcd A(2,2); (gdb) 9 { (gdb) 10 Eigen::MatrixXcd A(2,2); (gdb) 11 A << 1, 1, 4, -2; (gdb) 20 std::cout << A.pow(3.1) << "\n\n"; (gdb) print A $1 = Eigen::Matrix<std::complex<double>,2,2,ColMajor> (data ptr: 0x632010) = {[0,0] = {_M_value = 1 + 0 * I}, [1,0] = {_M_value = 4 + 0 * I}, [0,1] = {_M_value = 1 + 0 * I}, [1,1] = {_M_value = -2 + 0 * I}} (gdb) n Program received signal SIGSEGV, Segmentation fault. 0x00000000004068b7 in Eigen::MatrixPowerTriangularAtomic<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::compute2x2 (this=0x7fffffffd5d0, res=Eigen::Matrix<std::complex<double>,0,0,ColMajor> (data ptr: 0x0), p=0.10000000000000009) at /home/hazelnusse/usr/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h:330 330 res.coeffRef(0,0) = pow(m_A.coeff(0,0), p); (gdb) bt #0 0x00000000004068b7 in Eigen::MatrixPowerTriangularAtomic<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::compute2x2 (this=0x7fffffffd5d0, res=Eigen::Matrix<std::complex<double>,0,0,ColMajor> (data ptr: 0x0), p=0.10000000000000009) at /home/hazelnusse/usr/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h:330 #1 0x00000000004165db in Eigen::MatrixPowerTriangularAtomic<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::compute (this=this@entry=0x7fffffffd5d0, res=Eigen::Matrix<std::complex<double>,0,0,ColMajor> (data ptr: 0x0), p=0.10000000000000009) at /home/hazelnusse/usr/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h:304 #2 0x0000000000417432 in Eigen::MatrixPower<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::computeFracPower<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> > (this=this@entry=0x7fffffffd6e0, res=Eigen::Matrix<std::complex<double>,2,2,ColMajor> (data ptr: 0x6320c0) = {...}, p=<optimized out>) at /home/hazelnusse/usr/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h:459 #3 0x00000000004187ed in Eigen::MatrixPower<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::compute (this=this@entry=0x7fffffffd6e0, res=Eigen::Matrix<std::complex<double>,2,2,ColMajor> (data ptr: 0x6320c0) = {...}, p=1.0717734625362931) at /home/hazelnusse/usr/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h:335 #4 0x0000000000418919 in Eigen::MatrixPowerReturnValue<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> >::evalTo<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> > ( this=this@entry=0x7fffffffd8f0, res=Eigen::Matrix<std::complex<double>,2,2,ColMajor> (data ptr: 0x6320c0) = {...}) at /home/hazelnusse/usr/include/eigen3/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h:511 #5 0x00000000004189a6 in evalTo<Eigen::Matrix<std::complex<double>, -1, -1> > (dst=Eigen::Matrix<std::complex<double>,2,2,ColMajor> (data ptr: 0x6320c0) = {...}, this=0x7fffffffd8f0) at /home/hazelnusse/usr/include/eigen3/Eigen/src/Core/ReturnByValue.h:61 #6 Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1>::Matrix<Eigen::MatrixPowerReturnValue<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> > > (this=0x7fffffffd830, other=...) at /home/hazelnusse/usr/include/eigen3/Eigen/src/Core/Matrix.h:296 #7 0x0000000000418aaa in eval (this=0x7fffffffd8f0) at /home/hazelnusse/usr/include/eigen3/Eigen/src/Core/DenseBase.h:357 #8 Eigen::operator<< <Eigen::ReturnByValue<Eigen::MatrixPowerReturnValue<Eigen::Matrix<std::complex<double>, -1, -1, 0, -1, -1> > > > (s=..., m=...) at /home/hazelnusse/usr/include/eigen3/Eigen/src/Core/IO.h:244 #9 0x0000000000401de9 in main () at matrix_power.cc:20 (gdb) I'm not sure how to proceed next...
With your example given in bug 465 I observed it works fine if you use fixed size matrices, like Matrix2d instead of MatrixXd(2,2), and with MatrixXd it also works fine with integer powers.
Created attachment 348 [details] Proposed patch, including regression test
Comment on attachment 348 [details] Proposed patch, including regression test Review of attachment 348 [details]: ----------------------------------------------------------------- The patch solves the issue, however I think it should be specified better which method resizes its output parameter and which does not. E.g., computeBig() called in the default case automatically resizes, but case 1 does not (actually, case 1 would fail, but apparently is never called, because that case is handled elsewhere). I would accept the patch as a fix, but still leave the entire MatrixPower-module in a "needs review/cleanup"-status.
The function compute2x2 is frequently used in computeBig, so calling resize in compute2x2 is inefficient. I suggest calling res.resize(m_A.rows(), m_A.cols()); before switch. Sorry for my overlooking that res was not allocated.
(In reply to comment #4) > res.resize(m_A.rows(), m_A.cols()); This is equivalent to: res.resizeLike(m_A);
Created attachment 352 [details] Proposed patch, including regression test This also fixes the infinite loop in computeBig (see mailing list).
just one question: why limiting the iteration count to 64? I mean why 64? Is there any good justification for this choice? Is it reasonable for any matrix orders? Cannot we easily detect the non convergence?
Hmm... maybe just test if there is 0 on T's diag. T(i,i) == 0 seems okay because nonzero (even 1e-307) entries converge to 1 if continually sqrt'ed. I shouldn't have worried about "==" too much. The hard limit 64 is given in http://www.mathworks.com/matlabcentral/fileexchange/30527-arbitrary-real-power-of-a-matrix-by-schur-pade-algorithm/content/powerm_pade.m but not in the article. The MATLAB code is for double precision matrices.
Created attachment 355 [details] Proposed patch, including regression test Sorry that compute2x2 produced wrong result for adjacent equal eigenvalues.
Comment on attachment 355 [details] Proposed patch, including regression test Missing semicolon at eigen_assert(m_A(i,i) != RealScalar(0))
Patch looks fine, I pushed it to the repository. changeset: 5380:01867724d43d tag: tip parent: 5378:1e6d2684d0b5 user: Chen-Pang He <jdh8@ms63.hinet.net> date: Mon Jun 24 13:58:51 2013 +0100 summary: Fix segfault and bug with equal eivals in matrix power (bug 614). Regarding your FIXME, the documentation says that the matrix power is defined via the matrix logarithm, and the logarithm is not defined for singular matrices, so one might argue that the matrix power does not need to accept singular matrices. Nevertheless, this would of course be good to implement. If you want, you can open a new bug/enhancement request for this. I also want to repeat what Christoph said: (In reply to comment #3) > The patch solves the issue, however I think it should be specified better which > method resizes its output parameter and which does not. E.g., computeBig() > called in the default case automatically resizes, but case 1 does not > (actually, case 1 would fail, but apparently is never called, because that case > is handled elsewhere). > I would accept the patch as a fix, but still leave the entire > MatrixPower-module in a "needs review/cleanup"-status. Thus I intend to make a cleaning pass through the whole MatrixFunctions module (not just the matrix power).
-- GitLab Migration Automatic Message -- This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity. You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/614.