Hi.

I added explicit conversion to scalar, employed coeff and coeffRef, took care of updating blocks in splitOffTwoRows and pushDownZero, rewrote block operations you mentioned, changed indentation to two spaces.

There's slight performance decrease (like from 30.5 seconds to 31.5 seconds for solution of one thousand 100x100 problems on my laptop).
I'll take a look at unit tests.

Best

I don't see how HessenbergDecomposition would help... There's a way to speed up Hessenberg-Triangular reduction working with blocks instead of individual matrix elements, but I don't understand it well yet.

I added unit test.

hi,

I've pulled your changes into the main repository and granted you write access:

https://bitbucket.org/eigen/eigen/chang ... 8aaaba572/

I also took the liberty to update the license (to MPL2) and optimized a bit the hessenberg to not apply Givens on zero entries.

Regarding the performance, I get similar ones than Lapack. Not surprising since Lapack implements the same algorithm from the same sources.

Btw, in your repo there was a changeset hanging (29561ed775b7 Added random shifts, changed the way 0 on sub-diagonal of A is found, minor ). You can see what I mean there:

https://bitbucket.org/Khumarahn/eigen-qz/changesets

Since it is just before the "rewrite", I guess that was on purpose?

Hello,

there seems to be an issue when the input matrices are SPD. For instance, when adding:

Code: Select all

// test a SPD case
  {
    MatrixType a = MatrixType::Random(dim,dim);
    MatrixType b = MatrixType::Random(dim,dim);
    MatrixType a1 = MatrixType::Random(dim,dim);
    MatrixType b1 = MatrixType::Random(dim,dim);
    MatrixType spdA =  a.adjoint() * a + a1.adjoint() * a1;
    MatrixType spdB =  b.adjoint() * b + b1.adjoint() * b1;

    qz.compute(spdA, spdB);
    VERIFY_IS_APPROX(qz.matrixQ()*qz.matrixS()*qz.matrixZ(), A);
    VERIFY_IS_APPROX(qz.matrixQ()*qz.matrixT()*qz.matrixZ(), B);
  }


to the unit test, I get errors. The problem seems to be that complex eigenvalues are found while they are clearly all real. Here is a sequence of the matrix m_S during at each iteration:

Code: Select all

-2.262   -1.924   -1.376     3.13
  1.754    2.792  -0.5502   -1.909
      0   0.6373    1.499   -1.241
      0        0   0.1911   0.8552

   4.189    -1.371   -0.5542    -3.665
  0.0526     1.389    0.2854    -1.396
       0    0.2885   -0.6363   -0.2507
       0         0  -0.06798     1.614

   -4.172      1.205      3.254     -1.956
-0.001908     -1.339       1.39     0.0665
       -0  0.0004965     -1.647     0.2759
        0          0     0.1171     0.6983

   -4.171      1.206      2.492      2.865
-8.244e-05      -1.34     0.1751       1.38
       -0  1.712e-09    -0.8657    -0.9574
        0          0    -0.0857       1.27

<found a pair of complex eigenvalues>

 -4.171    1.206    2.492    2.865
-8.244e-05    -1.34   0.1751     1.38
     -0        0  -0.8657  -0.9574
      0        0  -0.0857     1.27


Since there is not much comments in the code, it's pretty hard to investigate more deeply.

Hello.

I will investigate the problem, it looks very interesting. I was going to add more comments, but had no time lately.

For commit 29561ed775b7, it was rather a mistake, and can be ignored. I would remove it if I knew how.

one more thing: I don't really understand how the shift vector (x, y, z) is computed: it seems you simply used the diagonal of the 2x2 block as its eigenvalues. That does not look correct, is there any hidden trick?

No, I compute x,y,z honestly. I take 2x2 trailing blocks SS and TT of S and T respectively, and compute eigenvalues of SS TT^{-1}. I found a good way to compute x,y,z in "An algorithm for generalized eigenvalue problems" by C.B.Moler and G.W.Stewart.

x,y,z are the non-zero elements in (M- lambda1 I)(M-lambda2 I) where M = S T^{-1} (or more precisely diagonal blocks of S and T where we currently work), and lambda1 and lambda2 are the eigenvalues of SS TT^{-1}

I found the problem. I appreciate your attention, Gael
https://bitbucket.org/eigen/eigen/chang ... d5d4c038f7
https://bitbucket.org/eigen/eigen/chang ... 70d48db6d6

Still there are symmetric matrices A and B

Code: Select all

A << 3.8125, -0.25, -2.375, -1.5625,
-0.25, 1.4375, -0.25, -0.5,
-2.375, -0.25, 3.1875, 0.0625,
-1.5625, -0.5, 0.0625, 1.625;
B << 2.875, -0.25, -1.1875, -1.8125,
-0.25, 1.5625, -0.125, -1.4375,
-1.1875, -0.125, 2.5625, -0.125,
-1.8125, -1.4375, -0.125, 2.9375;


for which algorithm finds complex eigenvalues:

Code: Select all

After 8 iterations:
S:
   -4.0343   -3.07279   -3.31122   0.770352
         0 0.00912231    2.04283  -0.417684
         0 0.00606289 -0.0219617   -1.78274
         0          0         -0    1.23393
T:
-2.92056 -1.88709 -2.37729  1.07864
       0 -0.71892  2.48593 -1.38954
       0        0  1.09379 -2.75667
       0        0       -0   1.8123
Q:
 -0.771923   -0.07831 -0.0382368  -0.629715
  0.166838   0.905745  -0.245869  -0.302223
  0.589391   -0.26662   0.286468  -0.706731
  0.170054   -0.32001  -0.925215  -0.112482
Z:
  0.243743  -0.659492  -0.608881  -0.367319
  0.865047 -0.0172427  0.0650458   0.497157
  0.215839   0.750663  -0.560001  -0.276255
 -0.381701 -0.0357335  -0.558061   0.735931
ERR: 7.0138e-19, 6.5009e-19, 7.28519e-19, 1.23276e-18


For errors on the bottom are |QSZ - A|/|A|, |QTZ-B|/|B| , |QQ*-I|, |ZZ*-I|