lae2.hpp

Go to the documentation of this file.

 
//
// Copyright (c) 2021-2023, University of Colorado Denver. All rights reserved.
//
// This file is part of <T>LAPACK.
// <T>LAPACK is free software: you can redistribute it and/or modify it under
// the terms of the BSD 3-Clause license. See the accompanying LICENSE file.
 
#ifndef TLAPACK_LAE2_HH
#define TLAPACK_LAE2_HH
 
#include "tlapack/base/utils.hpp"
 
namespace tlapack {
 
template <TLAPACK_SCALAR T>
void lae2(T a, T b, T c, T& s1, T& s2)
{
    // Constants
    const T zero(0);
    const T one(1);
    const T two(2);
    const T half(0.5);
 
    // Compute the eigenvalues
    T sm = a + c;
    T df = a - c;
    T adf = abs(df);
    T tb = b + b;
    T ab = abs(tb);
    T acmx, acmn;
    if (abs(a) > abs(c)) {
        acmx = a;
        acmn = c;
    }
    else {
        acmx = c;
        acmn = a;
    }
 
    T rt;
    if (adf > ab) {
        rt = adf * sqrt(one + square(ab / adf));
    }
    else if (adf < ab) {
        rt = ab * sqrt(one + square(adf / ab));
    }
    else {
        // This case includes case AB=ADF=0
        rt = ab * sqrt(two);
    }
 
    if (sm < zero) {
        s1 = half * (sm - rt);
        // Order of execution important.
        // To get fully accurate smaller eigenvalue,
        // next line needs to be executed in higher precision.
        s2 = (acmx / s1) * acmn - (b / s1) * b;
    }
    else if (sm > zero) {
        s1 = half * (sm + rt);
        // Order of execution important.
        // To get fully accurate smaller eigenvalue,
        // next line needs to be executed in higher precision.
        s2 = (acmx / s1) * acmn - (b / s1) * b;
    }
    else {
        // Includes case s1 = s2 = 0
        s1 = half * rt;
        s2 = -half * rt;
    }
}
}  // namespace tlapack
 
#endif  // TLAPACK_LAE2_HH