Hi I am trying to generate an arbitrary Gauss quadrature rule by using the Golub-Welsh algorithm.
I need to code this on C++ for my personal project. This algorithm involves the eigendecomposition of a matrix in which the only non-zero elements are the subdiagonal and superdiagonal. Here, $n$ could go up to $512$.
To illustrate in Matlab code:
n = 16;
beta = .5./sqrt(1-(2*(1:n-1)).-2);
T = diag(beta,1) + diag(beta,-1);
[V,D] = eig(T);
I want to implement the eigenvalue decomposition in C++ and not use Matlab routines for this since I want to parallelize it. What is the best way to find the eigendecomposition of this type of matrix? Is bisection method acceptable for my use case? How about divide and conquer or QR method or Lanczos?
(I was more inclined towards Lanczos since it seems that it is sparse friendly.)