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

Bug 430

Summary: Porting from Eigen 2 to 3 - runtime size assertion with Matrix product
Product: Eigen Reporter: Marcus D. Hanwell <marcus.hanwell>
Component: EigenvaluesAssignee: Nobody <eigen.nobody>
Status: NEW ---    
Severity: Unknown CC: chtz, gael.guennebaud, jacob.benoit.1, jitseniesen
Priority: Normal    
Version: 3.0   
Hardware: All   
OS: All   
Whiteboard:

Description Marcus D. Hanwell 2012-03-07 23:27:42 UTC
I ported OpenQube from Eigen 2 to 3 recently, and commit https://github.com/OpenChemistry/openqube/commit/1511ea03e020ca4e21653e1c570d9de8740ee5b4 shows the changes made, I am now seeing a runtime failure when compiling against Eigen 3.0.5 about size assertions failing. I added a little extra debug to try to get an idea of what is happening. The diff from master is,

diff --git a/src/slaterset.cpp b/src/slaterset.cpp
index 527b3d2..d95547d 100644
--- a/src/slaterset.cpp
+++ b/src/slaterset.cpp
@@ -29,6 +29,7 @@
 #include <Eigen/QR>
 
 #include <cmath>
+#include <iostream>
 
 #include <QtCore/QtConcurrentMap>
 #include <QtCore/QFuture>
@@ -249,16 +250,19 @@ void SlaterSet::calculationComplete()
   m_slaterShells[0].cube->lock()->unlock();
 }
 
+
 bool SlaterSet::initialize()
 {
   m_normalized.resize(m_overlap.cols(), m_overlap.rows());
 
   SelfAdjointEigenSolver<MatrixXd> s(m_overlap);
   MatrixXd p = s.eigenvectors();
-  MatrixXd m = p * s.eigenvalues().array().inverse().array().sqrt()
+  MatrixXd m = p * s.eigenvalues().array().inverse().sqrt()
                     .matrix().asDiagonal() * p.inverse();
   m_normalized = m * m_eigenVectors;
 
+  std::cout << "overlap: " << m_overlap << "\nm: " << m << std::endl;
+
   if (!(m_overlap*m*m).isIdentity())
     qDebug() << "Identity test FAILED - do you need a newer version of Eigen?";
   //    std::cout << m_normalized << std::endl << std::endl;

The output I see at runtime (along with the backtrace) is here,
File "/home/marcus/ssd/src/VTKData/Data/2h2o.aux" opened. 
Number of atoms = 6 
Number of atomic orbitals = 12 
Number of atomic orbital types = 12 
Number of zeta values = 12 
Number of PQN values = 12 
Number of electrons = 16 
Number of atomic coordinates = 18 
Size of lower half triangle of overlap matrix = 78 
Size of eigen vectors matrix = 144 
Size of lower half triangle of density matrix = 78 
No bonds found. Running simple bond perception...
Bonds found: 4
Num electrons: 16
overlap:       1       0       0       0       0  0.0002 -0.0002  0.0003    0.16  0.1596  0.0007  0.0186
      0       1       0       0 -0.0002  -0.002  0.0028 -0.0038  0.3877 -0.0974 -0.0037 -0.0494
      0       0       1       0  0.0002  0.0028 -0.0014  0.0034  0.0383 -0.3577  0.0029  0.0585
      0       0       0       1 -0.0003 -0.0038  0.0034 -0.0036 -0.0474 -0.1273 -0.0024 -0.0599
      0 -0.0002  0.0002 -0.0003       1       0       0       0  0.0006  0.0012  0.1625   0.155
 0.0002  -0.002  0.0028 -0.0038       0       1       0       0  0.0037  0.0038 -0.3188  0.2312
-0.0002  0.0028 -0.0014  0.0034       0       0       1       0 -0.0019 -0.0066  0.1721 -0.0987
 0.0003 -0.0038  0.0034 -0.0036       0       0       0       1  0.0027  0.0047  0.1574  0.2944
   0.16  0.3877  0.0383 -0.0474  0.0006  0.0037 -0.0019  0.0027       1  0.2314  0.0032  0.0551
 0.1596 -0.0974 -0.3577 -0.1273  0.0012  0.0038 -0.0066  0.0047  0.2314       1  0.0063  0.0669
 0.0007 -0.0037  0.0029 -0.0024  0.1625 -0.3188  0.1721  0.1574  0.0032  0.0063       1  0.2307
 0.0186 -0.0494  0.0585 -0.0599   0.155  0.2312 -0.0987  0.2944  0.0551  0.0669  0.2307       1
m:       1.0178    0.0161734   -0.0176657  -0.00976272  0.000172855   0.00074683  -0.00058752  0.000539823   -0.0745298   -0.0741942  0.000655697  -0.00217328
   0.0161734      1.08429    0.0371897   0.00369942  -0.00341627   -0.0065675   0.00224412  -0.00527425    -0.246753     0.108133  -0.00537548    0.0355131
  -0.0176657    0.0371897      1.06558    0.0171757   0.00407951   0.00802903  -0.00278621    0.0068215   -0.0686538     0.221512   0.00680877   -0.0437271
 -0.00976272   0.00369942    0.0171757      1.00825  -0.00267665    -0.004915   0.00161314  -0.00437983    0.0101869    0.0684497  -0.00503297    0.0296277
 0.000172855  -0.00341627   0.00407951  -0.00267665      1.01768  -0.00598006   0.00465166    0.0248341   0.00301178   0.00335374   -0.0768208   -0.0718596
  0.00074683   -0.0065675   0.00802903    -0.004915  -0.00598006      1.08885   -0.0441569    0.0149336   0.00555181   0.00624344     0.222118    -0.181809
 -0.00058752   0.00224412  -0.00278621   0.00161314   0.00465166   -0.0441569      1.02212  -0.00491714  -0.00265811 -9.33026e-05    -0.116389    0.0821662
 0.000539823  -0.00527425    0.0068215  -0.00437983    0.0248341    0.0149336  -0.00491714      1.04116   0.00536323   0.00566414   -0.0552142     -0.15631
  -0.0745298    -0.246753   -0.0686538    0.0101869   0.00301178   0.00555181  -0.00265811   0.00536323      1.11767    -0.166618    0.0052314   -0.0332231
  -0.0741942     0.108133     0.221512    0.0684497   0.00335374   0.00624344 -9.33026e-05   0.00566414    -0.166618      1.11792   0.00411964   -0.0385886
 0.000655697  -0.00537548   0.00680877  -0.00503297   -0.0768208     0.222118    -0.116389   -0.0552142    0.0052314   0.00411964      1.11753    -0.168632
 -0.00217328    0.0355131   -0.0437271    0.0296277   -0.0718596    -0.181809    0.0821662     -0.15631   -0.0332231   -0.0385886    -0.168632      1.12221
ChemistryCxxTests: /home/marcus/ssd/build/openchemistry/prefix/include/eigen3/Eigen/src/Core/ProductBase.h:154: typename Eigen::ProductBase<Derived, Lhs, Rhs>::Base::CoeffReturnType Eigen::ProductBase<Derived, Lhs, Rhs>::coeff(Eigen::ProductBase<Derived, Lhs, Rhs>::Index, Eigen::ProductBase<Derived, Lhs, Rhs>::Index) const [with Derived = Eigen::GeneralProduct<Eigen::GeneralProduct<Eigen::Matrix<double, -0x00000000000000001, -0x00000000000000001>, Eigen::Matrix<double, -0x00000000000000001, -0x00000000000000001>, 5>, Eigen::Matrix<double, -0x00000000000000001, -0x00000000000000001>, 5>, Lhs = Eigen::GeneralProduct<Eigen::Matrix<double, -0x00000000000000001, -0x00000000000000001>, Eigen::Matrix<double, -0x00000000000000001, -0x00000000000000001>, 5>, Rhs = Eigen::Matrix<double, -0x00000000000000001, -0x00000000000000001>, typename Eigen::ProductBase<Derived, Lhs, Rhs>::Base::CoeffReturnType = double, Eigen::ProductBase<Derived, Lhs, Rhs>::Index = long int]: Assertion `this->rows() == 1 && this->cols() == 1' failed.

Program received signal SIGABRT, Aborted.
0x00007ffff361b975 in raise () from /lib/libc.so.6
(gdb) bt
#0  0x00007ffff361b975 in raise () from /lib/libc.so.6
#1  0x00007ffff361cdeb in abort () from /lib/libc.so.6
#2  0x00007ffff36149ce in __assert_fail_base () from /lib/libc.so.6
#3  0x00007ffff3614a72 in __assert_fail () from /lib/libc.so.6
#4  0x00007ffff336002f in Eigen::ProductBase<Eigen::GeneralProduct<Eigen::GeneralProduct<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 5>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 5>, Eigen::GeneralProduct<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 5>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >::coeff (
    this=0x7fffffffcf70, row=0, col=0)
    at /home/marcus/ssd/build/openchemistry/prefix/include/eigen3/Eigen/src/Core/ProductBase.h:154
#5  0x00007ffff335ef91 in Eigen::DenseCoeffsBase<Eigen::GeneralProduct<Eigen::GeneralProduct<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 5>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 5>, 0>::coeff (
    this=0x7fffffffcf70, row=0, col=0)
    at /home/marcus/ssd/build/openchemistry/prefix/include/eigen3/Eigen/src/Core/DenseCoeffsBase.h:111
#6  0x00007ffff335d796 in Eigen::MatrixBase<Eigen::GeneralProduct<Eigen::GeneralProduct<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 5>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 5> >::isIdentity (
    this=0x7fffffffcf70, prec=9.9999999999999998e-13)
    at /home/marcus/ssd/build/openchemistry/prefix/include/eigen3/Eigen/src/Core/CwiseNullaryOp.h:714
#7  0x00007ffff3358897 in OpenQube::SlaterSet::initialize (this=0x45a3c0)
    at /home/marcus/ssd/src/openqube/src/slaterset.cpp:266
#8  0x00007ffff3357e2f in OpenQube::SlaterSet::calculateCubeMO (this=0x45a3c0, 
    cube=0x4b43d0, state=4)
    at /home/marcus/ssd/src/openqube/src/slaterset.cpp:154
#9  0x00007ffff3325a49 in OpenQube::BasisSet::blockingCalculateCubeMO (
    this=0x45a3c0, cube=0x4b43d0, mo=4)
    at /home/marcus/ssd/src/openqube/src/basisset.cpp:28
#10 0x00007ffff7fd0588 in vtkOpenQubeElectronicData::CalculateMO (
    this=0x45a0f0, orbitalNumber=4)
    at /home/marcus/ssd/src/VTK/Chemistry/vtkOpenQubeElectronicData.cxx:307
#11 0x00007ffff7fcfb4f in vtkOpenQubeElectronicData::GetMO (this=0x45a0f0, 
    orbitalNumber=4)
    at /home/marcus/ssd/src/VTK/Chemistry/vtkOpenQubeElectronicData.cxx:214
#12 0x0000000000422a78 in TestOpenQubeMOPACOrbital (argc=3, 
    argv=0x7fffffffe040)
    at /home/marcus/ssd/src/VTK/Chemistry/Testing/Cxx/TestOpenQubeMOPACOrbital.cxx:87
#13 0x000000000040ade3 in main (ac=3, av=0x7fffffffe040)
    at /home/marcus/ssd/build/VTK/Chemistry/Testing/Cxx/ChemistryCxxTests.cxx:283
Comment 1 Marcus D. Hanwell 2012-03-07 23:32:34 UTC
To add, if I move to the preceding commit (before porting) the code works as expected. That is against 2.0.17.
Comment 2 Christoph Hertzberg 2012-03-08 02:05:44 UTC
I guess the problem is that  in
>   if (!(m_overlap*m*m).isIdentity())

isIdentity checks the expression element-wise, which in Eigen3 is prohibited for Product-expressions (due to bad efficiency).

A work-around for you is to add an .eval() before .isIdentity()

The very-least fixing this bug needs is that a meaningful assert-message is added. Furthermore, isIdentity could call .eval() itself if the expression's complexity requires that.


Some unrelated side-notes: 
p.inverse() is the same as p.adjoint(), because p are orthonormal (for SelfAdjoint matrices). Furthermore you don't need to copy p, just declare p as const MatrixXd &.
Also, what you do looks exactly like operatorInverseSqrt():
http://eigen.tuxfamily.org/dox-devel/classEigen_1_1SelfAdjointEigenSolver.html#a811ad0873e06be5404fc91f64f0f658d
Comment 3 Gael Guennebaud 2012-03-08 08:27:40 UTC
if thta's really the issue, this should be fixed in isIdentity by honoring the prescribed nesting type
Comment 4 Marcus D. Hanwell 2012-03-08 23:18:02 UTC
Adding the eval call works, and fixes the immediate issue. I appreciate the side notes too - I will integrate your suggestions once things have settled down. I will leave it to you to decide on the most correct fix for the Eigen API - at a minimum I would like to see the assert message you suggested so that the source of the error when porting is a little more obvious.
Comment 5 Nobody 2019-12-04 11:31:26 UTC
-- 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/430.