Диагонализация самосопряженной (эрмитовой) матрицы 2×2

Диагонализация эрмитовой матрицы 2х2 проста, это можно сделать аналитически. Однако когда дело доходит до вычисления собственных значений и собственных векторов более чем в 10-6 раз, важно сделать это как можно более эффективным. Особенно, если недиагональные элементы могут исчезать, невозможно использовать одну формулу для собственных векторов: необходимо выражение if, что, конечно, замедляет код. Таким образом, я подумал, что использование Eigen, где указано, что диагонализация матриц 2×2 и 3×3 оптимизирована, все равно будет хорошим выбором:

с помощью

const std::complex<double> I ( 0.,1. );
inline double block_distr ( double W )
{
return (-W/2. + rand() * W/RAND_MAX);
}

тестовая петля будет

...
SelfAdjointEigenSolver<Matrix<complex< double >, 2, 2> > ces;
Matrix<complex< double >, 2, 2> X;

for (int i = 0 ; i <iter_MAX; ++i) {
a00=block_distr(100.);
a11=block_distr(100.);
re_a01=block_distr(100.);
im_a01=block_distr(100.);

X(0,0)=a00;
X(1,0)=re_a01-I*im_a01;
//only the lower triangular part is referenced! X(0,1)=0.; <--- not necessary
X(1,1)=a11;
ces.compute(X,ComputeEigenvectors);
}

Запись цикла без собственных значений с использованием непосредственно формул для собственных значений и собственных векторов эрмитовой матрицы и оператора if для проверки того, равна ли диагональ выключения нулю, в 5 раз быстрее. Я не использую Eigen должным образом или такие накладные расходы нормально? Существуют ли другие библиотеки lib.s, которые оптимизированы для небольших самосопряженных матриц?