Iterative reduction to Hessenberg upper triangular form

[thijssteel]

Description

This PR adds an implementation of an iterative algorithm for the reduction of a pencil to Hessenberg-upper triangular form. It leverages the connection between generalized and standard eigenvalue problems. This introduces accuracy problems, which are then solved using iterative refinement. For a pencil (A,B), the pseudocode is:

• X = AB^{-1}
• Q^*XQ = H (Hessenberg reduction of X)
• A <- QA; B <- QB
• B = RZ (RQ decomposition of B)
• A <- AZ*; B <- BZ*

This process is repeated untill (A,B) is in Hessenberg-upper triangular form.

Because the reduction of a matrix to Hessenberg form is highly optimized, this algorithm can result in significant speedups despite needing multiple iterations.

I wrote a paper on this algorithm with prof. Vandebril, currently in review. The journal is not clear on preprint policy, otherwise i'd share a link.

Hybrid

The iterative refinement scheme is not guaranteed to work, especially not for pencils that have many infinite eigenvalues. To dampen this effect, my implementation performs some of the work using XGGHD3. As more iterations fail, more and more of the work is done using XGGHD3 so the scheme is guaranteed to converge and will never be more than a constant factor slower than XGGHD3.

Timings

Wall time of the algorithm performed on a compute server with two intel Xeon E5-2650 v2 cpu's. Pencil consist of two matrices of size n with randomly drawn entries. Programs were linked with MKL.

This test is a best case, where not many refinement steps were needed. For pencils with many infinite eigenvalues, DGGHD4 can be up to twice as slow as DGGHD3. In that case, DGGHD4 essentially performs all of the work relying on DGGHD3, but introduces extra overhead.

n	DGGHD4	DGGHD3
500	0.11433215600000000	0.23248758600000000
707	0.24008591900000001	0.57098310500000005
1000	0.52512665599999997	1.5334543110000001
1414	1.1062537070000000	3.9855276819999998
2000	2.6581395319999999	11.229642351000001
2828	7.2890489220000001	29.942265634000002
4000	20.061862538000000	83.943255199999996
5657	54.722602287999997	217.99881286100000
8000	156.28234842000001	625.75760643599995

Draft

This PR is a draft and only contains the double precision implementation. I will add the other types and tests later if the feedback on this implementation is positive.

[codecov[bot]]

Codecov Report

Merging #645 (f859e30) into master (1782d90) will not change coverage. The diff coverage is 0.00%.

:exclamation: Current head f859e30 differs from pull request most recent head 1bd7d49. Consider uploading reports for the commit 1bd7d49 to get more accurate results

@@           Coverage Diff           @@
##           master     #645   +/-   ##
=======================================
  Coverage    0.00%    0.00%           
=======================================
  Files        1894     1894           
  Lines      184034   184047   +13     
=======================================
- Misses     184034   184047   +13     

Impacted Files	Coverage Δ
SRC/dgghd3.f	0.00% <0.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 1782d90...1bd7d49. Read the comment docs.