New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
Bug 430 - Porting from Eigen 2 to 3 - runtime size assertion with Matrix product
Summary: Porting from Eigen 2 to 3 - runtime size assertion with Matrix product
Status: NEW
Alias: None
Product: Eigen
Classification: Unclassified
Component: Eigenvalues (show other bugs)
Version: 3.0
Hardware: All All
: Normal Unknown
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks:
 
Reported: 2012-03-07 23:27 UTC by Marcus D. Hanwell
Modified: 2014-04-10 22:08 UTC (History)
4 users (show)



Attachments

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.

Note You need to log in before you can comment on or make changes to this bug.