1. Alexey Korepanov
  2. eigen-QZ

Commits

Alexey Korepanov  committed 29561ed

Added random shifts, changed the way 0 on sub-diagonal of A is found, minor changes.

  • Participants
  • Parent commits bb3bfec
  • Branches default

Comments (0)

Files changed (1)

File Eigen/src/Eigenvalues/RealQZ.h

View file
  • Ignore whitespace
     * check for complex numbers, as we don't do complex here
     * what if we use FullPivHouseholderQr? Will result improve?
     * or may be there's some sort of balancing for bad matrices?
+    * chase precision loss
 */
 
 #include <algorithm>
 #include <vector>
 
+//#define CHECK_A_ERR { Scalar errA=(Q*A*Z - A_).norm(); std::cout << __LINE__ << ", errA: " << errA << "\n"; }
+#define CHECK_A_ERR
+
 #include <Eigen/Dense>
 using namespace Eigen;
 
             Q(size,size),
             Z(size,size),
             workspace(size),
-            m_computed(false) { }
+            m_computed(false),
+            m_converged(false) { }
 
         // constructive constructor
         RealQZ(const MatrixType& A_, const MatrixType& B_, bool computeQZ = true) :
             Z(A_.rows(),A_.cols()),
             workspace(A_.rows()),
             m_computed(false),
+            m_converged(false),
             m_computeQZ(computeQZ) { 
                 compute(A_, B_, computeQZ);
             }
         // actual computation routine
         void compute(const MatrixType& A_, const MatrixType& B_, bool computeQZ = true);
 
+        // return true if converged, false else
+        bool converged() const { 
+            eigen_assert(m_computed && "no computation done yet");
+            return m_converged; 
+        }
+
         // get matrix Q
         const MatrixType & matrixQ() const {
             eigen_assert(m_computed && m_computeQZ && "matrix Q was not computed");
         MatrixType A, B, Q, Z;
         Matrix<Scalar,Dynamic,1> workspace;
         bool m_computed;
+        bool m_converged;
         bool m_computeQZ;
 
         unsigned long iter;
-        bool converged;
 
         struct {
             bool operator() (Scalar a, Scalar b) { return internal::abs(a) < internal::abs(b); }
     // entrance point to Francis-like QR iteration - Hessenberg-Triangular form
     hessenbergTriangular();
 
-    converged = true;
+    m_converged = false;
 
     // Big Loop
     for (;;) {
         //     | 0   0  A33|
         // then A33 is upper quazi-triangular, and A22 is unreduced upper Hessenberg;
         // B is partitioned accondingly
+CHECK_A_ERR
         pqr(p,q,r);
+CHECK_A_ERR
+
+//std::cout << "new iteration; p: " << p << ", q: " << q << ", r: " << r << ", A:\n"<< A << "\n";
+//std::cout << "A B^{-1}:\n" << A*B.inverse() << "\n\n";
+//sleep(1);
 
         if (q==dim) {
             // if q==dim, we are done
+            m_converged = true;
             break;
         } else {
             // otherwise we work on A22, B22
                 }
                 
             } else {
-                // B_22 is not singular
-                // compute (x,y,z,0..0)* - the first column of A22 B22^{-1}
+                // B22 is not singular, A22 is unreduced upper-Hessenberg
+                // Let M = A22 B22^{-1}, and compute first column (x y z 0...0) of
+                // (M - l1 I) (M-l2 I) where l1 and l2 are eigenvalues of
+                // bottom-right 2x2 block of M; l1 and l2 may be complex,
+                // but result is always real
                 Scalar x=0.0, y=0.0, z;
                 if (r==3) {
+                    // special case when B22 and A22 are 3x3 matrices, 
+                    // computations are somewhat easier here
                     std::vector<Scalar> xx(9), yy(3);
                     const Scalar
-                        a11=A(p+0,p+0),
-                        a12=A(p+0,p+1),
-                        a21=A(p+1,p+0),
-                        a22=A(p+1,p+1),
-                        a23=A(p+1,p+2),
-                        a32=A(p+2,p+1),
-                        a33=A(p+2,p+2),
-                        b12=B(p+0,p+1),
-                        b13=B(p+0,p+2),
+                        a11=A(p+0,p+0), a12=A(p+0,p+1), 
+                        a21=A(p+1,p+0), a22=A(p+1,p+1), a23=A(p+1,p+2),
+                        a32=A(p+2,p+1), a33=A(p+2,p+2),
+                        b12=B(p+0,p+1), b13=B(p+0,p+2),
                         b23=B(p+1,p+2),
                         b11i=1.0/B(p+0,p+0),
                         b22i=1.0/B(p+1,p+1),
                     // r>=4
                     std::vector<Scalar> xx(11), yy(7);
                     const Scalar
-                        a11=A(p+0,p+0),
-                        a12=A(p+0,p+1),
-                        a21=A(p+1,p+0),
-                        a22=A(p+1,p+1),
-                        a32=A(p+2,p+1),
+                        a11=A(p+0,p+0), a12=A(p+0,p+1),
+                        a21=A(p+1,p+0), a22=A(p+1,p+1), a32=A(p+2,p+1),
                         b12=B(p+0,p+1),
                         b11i=1.0/B(p+0,p+0),
                         b22i=1.0/B(p+1,p+1),
+                        a87=A(l22-1,l22-2), a88=A(l22-1,l22-1), a89=A(l22-1,l22-0),
+                        a98=A(l22-0,l22-1), a99=A(l22-0,l22-0), 
+                        b78=B(l22-2,l22-1), b79=B(l22-2,l22-0),
+                        b89=B(l22-1,l22-0),
                         b77i=1.0/B(l22-2,l22-2),
                         b88i=1.0/B(l22-1,l22-1),
-                        b99i=1.0/B(l22-0,l22-0),
-                        a99=A(l22-0,l22-0),
-                        a88=A(l22-1,l22-1),
-                        a89=A(l22-1,l22-0),
-                        a98=A(l22-0,l22-1),
-                        a87=A(l22-1,l22-2),
-                        b89=B(l22-1,l22-0),
-                        b78=B(l22-2,l22-1),
-                        b79=B(l22-2,l22-0);
+                        b99i=1.0/B(l22-0,l22-0);
                     xx[0] = a11*a11*b11i*b11i;
                     xx[1] = a12*a21*b11i*b22i;
                     xx[2] = -a11*a21*b12*b11i*b11i*b22i;
                     z = a21*a32*b11i*b22i;
                 }
 
+                // on iterations 15, 31, 47 and so modify x,y,z randomly
+                if (((iter-1)%16)==1) {
+                    Scalar xr = internal::random<Scalar>(0.5,2.0),
+                           yr = internal::random<Scalar>(0.5,2.0),
+                           zr = internal::random<Scalar>(0.5,2.0);
+                    x *= xr;
+                    y *= yr;
+                    z *= zr;
+                }
+
                 // variables for Householder rotations
                 Vector2s essential2;
                 Vector1s essential1;
                 // update Z
                 applyH_l<1>(Z,l22,l22-1,essential1,tau);
 
-                std::cout << "\n\n************************************\n"
-                    << "loop over k completed (iter=" << iter << ", p=" << p << ", l22=" << l22 << ")\n"
-                    << "*************************************\n";
-                std::cout << "A:\n" << A << "\nB:\n" << B << "\nA B^{-1}:\n" << A*B.inverse() << "\n\n";
+                //std::cout << "\n\n************************************\n"
+                //    << "loop over k completed (iter=" << iter << ", p=" << p << ", l22=" << l22 << ")\n"
+                //    << "*************************************\n";
+                //std::cout << "A:\n" << A << "\nB:\n" << B << "\nA B^{-1}:\n" << A*B.inverse() << "\n\n";
                 //sleep(2);
                 
-                if (iter++>20) {
-                    std::cout << "**********************************\nLeaving on iteration " 
-                        << iter << "\n***********************************\n\n";
-                    converged = false;
+                if (iter++>2000) {
+                    //std::cout << "**********************************\nLeaving on iteration " 
+                    //    << iter << "\n***********************************\n\n";
                     break;
                 }
             } // B22 is not singular
         } // q!=dim
     } // Big Loop
 
-    if (converged) {
+    if (m_converged) {
         // work on 2x2 blocks with real eigenvalues
         for (Index i=dim-2; i>=0; i--) {
             if (A(i+1,i) != 0.0 && B(i,i)!=0.0 && B(i+1,i+1)!=0.0) {
 
             }
         }
-
     }
 
     m_computed = true;
 template<typename MatrixType>
 void RealQZ<MatrixType>::hessenbergTriangular() {
     
-    std::cout << "Before Hessenberg-Triangular:\nA:\n" << A << "\nB:\n" << B << "\n\n";
+    //std::cout << "Before Hessenberg-Triangular:\nA:\n" << A << "\nB:\n" << B << "\n\n";
     
     // perform QR decomposition of B, overwrite B with R, save Q
     HouseholderQR<MatrixType> qrB(B);
         }
     }
     
-    std::cout << "After Hessenberg-Triangular:\nA:\n" << A << "\nB:\n" << B << "\nA B^{-1}:\n" << A*B.inverse() << "\n\n";
+    //std::cout << "After Hessenberg-Triangular:\nA:\n" << A << "\nB:\n" << B << "\nA B^{-1}:\n" << A*B.inverse() << "\n\n";
 }
 
 
         char nz=0;
         // compute p,q
         for (Index i=dim-1; i>0; i--) {
-            if ( internal::isMuchSmallerThan(A(i,i-1), internal::abs(A(i-1,i-1))+internal::abs(A(i,i)) ) ) {
+            if (internal::abs(A(i,i-1)) <= NumTraits<Scalar>::epsilon() * (internal::abs(A(i-1,i-1))+internal::abs(A(i,i))) ) { 
+//std::cout << "Zeroing A("<<i<<","<<i-1<<") in \n" << A/*.block(i-1,i-1,2,2)*/ << "\n";
                 A(i,i-1) = 0.0;
                 if (nz==1) {
                     nz = 0;