mlimberger
Registered Member

Hello.
So far, I used SimplicialLDLT to solve for linear equations system with a sparse, symmetric, positive definite lefthand side (LHS). Now, the LHS has to be extended and becomes sparse, symmetric but not positive definite anymore. The size is around 15000 x15000. The structure of LHS has the characteristic, that all diagonal elements but one have nonzero entries: diag(LHS) = [x,x,x,x,x,x,x,x,...,x,x,0] The RHS then is [x,x,x,x,x,...,x,x,0]. Therefore I switched from SimplicialLDLT (which is only valid for symmetric positive definite problems, right?) to SparseLu from which I thought it is mybe the best way for solving. Can you agree or do you think there are better alternatives? Based on x = LHS^(1)*RHS, my implementation now looks like:  ... #include "Eigen/SparseLU" #include "Eigen/OrderingMethods" ... x = Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int> >(LHS).solve(RHS);  Unfortunetaly I receive the following error message: Eigen::DenseCoeffsBase<Derived, 1>::Scalar& Eigen::DenseCoeffsBase<Derived, 1>::operator[](Eigen::DenseCoeffsBase<Derived, 1>::Index) [with Derived = Eigen::Ref<Eigen::Matrix<int, 0x00000000000000001, 1>, 0, Eigen::InnerStride<1> >, Eigen::DenseCoeffsBase<Derived, 1>::Scalar = int, Eigen::DenseCoeffsBase<Derived, 1>::Index = long int]: Assertion `index >= 0 && index < size()' failed. Sounds for me like the reson for the error is based on the COLAMDOrdering. I also testet NaturalOrdering and AMDOrdering, but the same error message appeared. I would greatly appreciate it if you could give me some help, support or comments on how to solve this kind of least squares in an efficient way. Thank you very much in advance! Marco 
ggael
Moderator

Do you use the devel branch? or the 3.2beta1? There has been a recent fix in SparseLU. If that still fails with the latest revision, could you send the matrix?

mlimberger
Registered Member

Hi ggael.
Thanks for your fast reply! I use the Eigen 3.2beta1 version. I tested SparseLu with a small 3x3 LHS test matrix where the same problem appears: Eigen::MatrixXd lhs(3,3); lhs.coeffRef(0,0)=3; lhs.coeffRef(0,1)=11; lhs.coeffRef(0,2)=2; lhs.coeffRef(1,0)=11; lhs.coeffRef(1,1)=5; lhs.coeffRef(1,2)=2; lhs.coeffRef(2,0)=2; lhs.coeffRef(2,1)=2; lhs.coeffRef(2,2)=0; and 3x1 RHS: Eigen::VectorXd rhs(3,1); rhs<<1,1,1; In case of using a dense matrix, the result yields: x = lhs.inverse()*rhs = [0.266667, 0.233333, 1.38333]  But when I use SparseLu with Eigen::SparseMatrix<double>lhs(3,3); ... same entries... and x = Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int> >(lhs).solve(rhs); the error message appears. Any idea why SparseLu cannot solve it? Or do you anyway think that SparseLu is not the the best way for such a linear equations system? Regards, Marco 
ggael
Moderator

please, try with the devel branch: http://bitbucket.org/eigen/eigen/get/default.tar.bz2. This issue is likely already fixed.

mlimberger
Registered Member

Same Error with the devel branch Hmmm...Do you have any suggestion what other solver I could try?
The important thing is that it works properly also for relatively large systems. 
ggael
Moderator

I can reproduce the issue with the beta1 but it does work fine with the devel branch so make sure you properly configured your compiler to use the devel branch and not another copy of Eigen.
Here is my test program:

mlimberger
Registered Member

The compiler was configured correctly. But you are right, your Testprogram is working fine. My testprogram was:
In this case, the error message appears. But I realized that if I add an additional lhs.makeCompressed(); is works:
For the real case with the large matrix, I initialize the matrix also with the expected size and add single elements by LHS.coeffRef(i,j). Can this cause an error when using SparseLu? At least, the error still remains here, even if I compress the matrix. 
ggael
Moderator

hm, right SparseLU is expecting a compressed matrix.

mlimberger
Registered Member

Unfortunately, it still does not work for the larger matrix. The error message that arises is:
const Eigen::internal::solve_retval<Eigen::SparseLU<_MatrixType, _OrderingType>, Rhs> Eigen::SparseLU<_MatrixType, _OrderingType>::solve(const Eigen::MatrixBase<OtherDerived>&) const [with Rhs = Eigen::Matrix<double, 0x00000000000000001, 1>, _MatrixType = Eigen::SparseMatrix<double>, _OrderingType = Eigen::COLAMDOrdering<int>]: Assertion `m_factorizationIsOk && "SparseLU is not initialized."' failed. Using dense matrices with x = LHS.jacobiSvd(Eigen::ComputeThinUEigen::ComputeThinV).solve(RHS); works. Is there a convenient way that I send you the LHS matrix (733x733) and RHS vector (733x1)? I am really stucked with this problem. 
Dee33
Registered Member

Please use the functions saveMarket() and saveMarketVector() defined in unsupported/Eigen/src/SparseExtra/MarketIO.h
Then file a bug here and attach the required files : http://eigen.tuxfamily.org/bz/ 
mlimberger
Registered Member

Thank you so much for your support! I opened a Bug session. So let's see...

Dee33
Registered Member

SparseLU now accepts non compressed matrices : https://bitbucket.org/eigen/eigen/commits/02d91f29a8ff/ Changeset: 02d91f29a8ff User: dnuentsa Date: 20130716 15:15:53 Summary: Fix Sparse LU for matrices in non compressed mode Affected #: 1 file For following posts : visit http://eigen.tuxfamily.org/bz/show_bug.cgi?id=630 
tonyhenz
Registered Member

Hi Dee33 and all,
Now i have met exact the same problem. According to the previous discussions, i have used the newest version of Eigen, for example 3.2rc1. However, this issue still exists. But it is a little strange, i have tested my problem, if the sparse matrix dimension is small, to say 70000 * 70000, the SparseLU works fine, but when the dimension increases to 200000 * 200000, the error presented above will show again. I have tested my problem by using BiCGSTAB solver, it works fine but quite slowly. Any suggestions will be appreciated. 
ggael
Moderator

You mean you get this exact same error message:
? You can also check you are really using 3.2rc1 by printing cout << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION; You should get 1.92. 
tonyhenz
Registered Member

Hi ggael, Thanks for your prompt reply. I have printed the output, which is 1.92. And i found the error message is not exactly the same, but slightly different, please see below. The strange thing is, for a different matrix dimension it works quite nice, so i am confused. ../eigen3/Eigen/src/Core/PlainObjectBase.h:265: void Eigen::PlainObjectBase<Derived>::resize(typename Eigen::internal::traits<T>:$$ts<T>::Index) [with Derived = Eigen::Matrix<int, 0x00000000000000001, 1, 0, 0x00000000000000001, 1>]: Assertion `((SizeAtCompileTime == D$$e == Dynamic && (MaxSizeAtCompileTime==Dynamic  size<=MaxSizeAtCompileTime))  SizeAtCompileTime == size) && size>=0' failed. 
Registered users: andreas_k, Baidu [Spider], Bing [Bot], elvisangelaccio, frankragr, Google [Bot], google01103, hugo.pereira@free.fr, jakedarcy, kdeuserk, koriun, Kver, lruffini, Majestic12 [Bot], mmolch, nuttshell, paviluf, pedrorodriguez, predatux, ramskulls, rv8ter, Saabhero, Steve Guilford, wojtask, Xiceph, Yahoo [Bot]