laed5.hpp

Go to the documentation of this file.

 
//
// Copyright (c) 2025, 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_LAED5_HH
#define TLAPACK_LAED5_HH
 
#include "tlapack/base/utils.hpp"
 
namespace tlapack {
 
template <class d_t, class z_t, class delta_t, class real_t, class idx_t>
void laed5(idx_t i, d_t& d, z_t& z, delta_t& delta, real_t rho, real_t& dlam)
{
    real_t w, b, c, tau, temp;
    real_t del = d[1] - d[0];
    if (i == 0) {
        w = real_t(1.0) + real_t(2.0) * rho * (z[1] * z[1] - z[0] * z[0]) / del;
 
        if (w > 0) {
            b = del + rho * (z[0] * z[0] + z[1] * z[1]);
            c = rho * z[0] * z[0] * del;
 
            // B > 0, always
 
            tau = real_t(2.0) * c / (b + sqrt(abs(b * b - real_t(4.0) * c)));
            dlam = d[0] + tau;
            delta[0] = -z[0] / tau;
            delta[1] = z[1] / (del - tau);
        }
        else {
            b = -del + rho * (z[0] * z[0] + z[1] * z[1]);
            c = rho * z[1] * z[1] * del;
 
            if (b > 0) {
                tau = real_t(-2.0) * c / (b + sqrt(b * b + real_t(4.0) * c));
            }
            else {
                tau = (b - sqrt(b * b + real_t(4.0) * c)) / real_t(2.0);
            }
 
            dlam = d[1] + tau;
            delta[0] = -z[0] / (del + tau);
            delta[1] = -z[1] / tau;
        }
 
        temp = sqrt(delta[0] * delta[0] + delta[1] * delta[1]);
        delta[0] = delta[0] / temp;
        delta[1] = delta[1] / temp;
    }
    else {
        // Now I = 2
        b = -del + rho * (z[0] * z[0] + z[1] * z[1]);
        c = rho * z[1] * z[1] * del;
 
        if (b > 0) {
            tau = (b + sqrt(b * b + real_t(4.0) * c)) / real_t(2.0);
        }
        else {
            tau = real_t(2.0) * c / (-b + sqrt(b * b + real_t(4.0) * c));
        }
 
        dlam = d[1] + tau;
        delta[0] = -z[0] / (del + tau);
        delta[1] = -z[1] / tau;
        temp = sqrt(delta[0] * delta[0] + delta[1] * delta[1]);
        delta[0] = delta[0] / temp;
        delta[1] = delta[1] / temp;
    }
 
    return;
}
}  // namespace tlapack
 
#endif  // TLAPACK_LAED5_HH