I Need to to solve an generalized Eigenvalue Problem and to compare two Methods (QR and QZ) concerning their convergence rate and execution time.
I started with the Basic QR-Algorithm, implemented in Matlab, tested the results and in the last step i included double-shifts and Deflation strategies.
Now it is time to improve my QZ-Algorithm: I think ist way more complicated, as i did not find any QZ-Algorithm implementations in the Internet. So i grabed the paper of Moler and Stewart and wrote my own implementation. I've been very glad that i did manage to write an algorithm producing the correct results. To have a fair comparison i would like to enhance this algorithm as to include single/double shifts so i can compare the Performance between the QR- and QZ-Algorithm. I think the implicit shift should be more stable, but i want the fastest Solution (less iterations).
You can find the paper here: http://www.dtic.mil/dtic/tr/fulltext/u2/746896.pdf And an extensional paper here: http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19730018792.pdf
This is my algorithm so far (works correctly, but without shifting, so it is quite slow):
function [lambda, iteration] = QZAlgo(A, B, precision, maxIterations, useHouseholder) len = length(A); % 1. Step: Preparation of Matrices % Transform A to Hessenberg form and B to upper triangular form if (useHouseholder ~= 0) [Q, B] = TriangularFormByHouseholder(B); else [Q, B] = TriangularFormByGivens(B); Q = Q'; end A = Q*A; [A, B] = HessenbergFormByGivensWhileMaintainingB(A, B); % 2. Step: QZ-Iteration for iteration=1:maxIterations p = 1; S = A - p*B; if (useHouseholder ~= 0) [Q, R] = TriangularFormByHouseholder(S); else [Q, R] = TriangularFormByGivens(S); Q = Q'; end B1 = Q*B; [B2, Z] = RepairTriangularFormByGivensFromRight(B1); B = Q*B*Z; A = Q*A*Z; lambda = diag(A) ./ diag(B); if (CheckForRequiredPrecision(A, B, precision)) break; end end end
The Triangular-Functions are basically QR-Decompositions. The HessenbergFormByGivensWhileMaintainingB-Function creates an QR-Decomposition of A while repairing changes of B with Givens-Transformations of the right side.
The Algorithm works so far, but is, as already mentioned, quite slow.
Do anybody know how i can implement some shifting-strategies to the algorithm?
ANY help, no matter of what Kind and no matter about which shifting strategy, is GREATLY appreciated. I will spend you a crate of beer of your choice, if we can manage to get this working, if you like ;-)
Thanks in Advance!
EDIT: If you Need to have a look in some of the subfuctions, i can provide them here, too.
EDIT2: I tried to have a look at CLAPACK (http://www.netlib.org/clapack/) to figure it out, but the Code is so confusing to me that i had to give it up...
EDIT3: Some of my thoughts about the explicit shift
- Choose shift p: I've Chosen the element in the right bottom Corner of S=A-B as i thought it must be similiar to the procedure of the QR-Algorithm
- Calculate QR-Decomposition of A-p*B
- Calculate B'=QBZ just the same way i did above
- Calculate A=SZ + pB = Q(A-pB) + p*B' But this Procedure did not work. I am clueless at this Point :-(