Hi,

I am currently dealing with two problems:

1. I have to calculate a LU decomposition for somewhat ill defined matrices. Python and Mathematica however to do not have any problems dealing with them. If you run the following the code, L for example is not unit diagonal. What am I doing wrong here?

Code: Select all

#include <iostream>
#include <Eigen/Core>
#include <Eigen/LU>

using namespace std;
using namespace Eigen;

int main() {
    MatrixXd to_det(12,12);
    int sites = 12;

    to_det <<
        1.2431e+16, -3.94892e+15,  1.30443e+15, -4.57433e+14,  1.29723e+15, -1.84442e+15,  1.22654e+16, -9.04374e+15,            0,            0,            0,      0,
        -2.32088e+16,  7.37345e+15, -2.43579e+15,  8.54133e+14, -2.42194e+15,  3.44331e+15,  -2.2898e+16,  1.68838e+16,            0,            0,            0,    0,
        5.121e+15, -1.62693e+15,  5.37446e+14, -1.88462e+14,  5.34398e+14, -7.59769e+14,  5.05245e+15, -3.72542e+15,            0,            0,            0,      0,
        -2.45812e+15,  7.80795e+14, -2.57906e+14,  9.04445e+13, -2.56516e+14,  3.64736e+14, -2.42549e+15,  1.78839e+15,            0,            0,            0,      0,
        6.84946e+15, -2.17507e+15,   7.1835e+14, -2.51945e+14,  7.14782e+14, -1.01649e+15,   6.7597e+15, -4.98399e+15,            0,            0,            0,        0,
        -2.0641e+15,  6.55336e+14, -2.16413e+14,  7.59079e+13, -2.15402e+14,  3.06359e+14, -2.03729e+15,  1.50208e+15,            0,            0,            0,       0,
        2.3618e+15, -7.49896e+14,  2.47647e+14, -8.68613e+13,  2.46469e+14, -3.50533e+14,  2.33105e+15, -1.71868e+15,            0,            0,            0,       0,
        -1.23279e+16,  3.91476e+15, -1.29291e+15,  4.53457e+14, -1.28649e+15,  1.82952e+15, -1.21663e+16,  8.97036e+15,            0,            0,            0,      0,
        0,            0,            0,            0,            0,            0,            0,            0,            2,            0,            0,            0,
        0,            0,            0,            0,            0,            0,            0,            0,            0,            2,            0,            0,
        0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            2,            0,
        0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            0,            2;

    FullPivLU<MatrixXd> lu(sites, sites);
    lu = lu.compute(to_det);
    cout << "Matrix LU" << endl;
    cout << lu.matrixLU() << endl;
    double det = lu.determinant();
    cout << det << endl;
   MatrixXd u = lu.matrixLU().triangularView<Upper>();
    cout << "U" << endl;
    cout << u << endl;
    MatrixXd l(sites, sites);
    l.Identity(sites, sites);
    l.block(0,0,sites,sites).triangularView<StrictlyLower>() = lu.matrixLU();
    cout << "L" << endl;
    cout << l << endl;
    return 0;
}


2. When doing the QR decomposition, is there any way to gain access to the pivot vector that is available in the equivalent matlab routine?

Thanks a lot!

1- The Idendity() function is a static member of the class Matrix that returns an identity matrix. to set the matrix l to the identity, use l.setIdentity(sites,sites);

2- Use qr.colsPermutation() to get the respective PermutationMatrix.

Now what's the numerical issue you are hitting? I tried your matrix and with FullPivLU, the reconstructed matrix just match the input, and solving with a random right hand side, works well too.

Hello,

I am having a similar problem with JacobiSVD where when i am trying to compute using the FullPivHouseholderQRPreconditioner and since it does not run with Thin U|V , i compute Full U|V to finally use solve . This leads me to an error stating that the assertion for U & V has failed . Can you please point me in the direction as how to use Full PivQR with Jacobi SVD to get the decomposition .

Hi,

sorry for the delay and thanks for the answers! There were two issues: 1. I expected the LU decomposition to give me the LDU decomposition and 2. I missed the last digit when copying the file, giving errorneous results in my code compared to numpy and mathematica. It's all good now!