move_bulge.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_MOVE_BULGE_HH
#define TLAPACK_MOVE_BULGE_HH
 
#include "tlapack/base/utils.hpp"
#include "tlapack/lapack/lahqr_shiftcolumn.hpp"
#include "tlapack/lapack/larfg.hpp"
 
namespace tlapack {
 
template <TLAPACK_SMATRIX matrix_t, TLAPACK_CVECTOR vector_t>
void move_bulge(matrix_t& H,
                vector_t& v,
                complex_type<type_t<matrix_t>> s1,
                complex_type<type_t<matrix_t>> s2)
{
    using T = type_t<matrix_t>;
    using real_t = real_type<T>;
 
    using idx_t = size_type<matrix_t>;
    using range = pair<idx_t, idx_t>;
    const real_t zero(0);
    const real_t eps = ulp<real_t>();
 
    CreateStatic<vector_t, 3> new_3_vector;
 
    // Perform delayed update of row below the bulge
    // Assumes the first two elements of the row are zero
    scalar_type<type_t<vector_t>, T> refsum = v[0] * v[2] * H(3, 2);
    H(3, 0) = -refsum;
    H(3, 1) = -refsum * conj(v[1]);
    H(3, 2) = H(3, 2) - refsum * conj(v[2]);
 
    // Generate reflector to move bulge down
    T tau, beta;
    v[0] = H(1, 0);
    v[1] = H(2, 0);
    v[2] = H(3, 0);
    larfg(FORWARD, COLUMNWISE_STORAGE, v, tau);
    beta = v[0];
    v[0] = tau;
 
    // Check for bulge collapse
    if (H(3, 0) != zero or H(3, 1) != zero or H(3, 2) != zero) {
        // The bulge hasn't collapsed, typical case
        H(1, 0) = beta;
        H(2, 0) = zero;
        H(3, 0) = zero;
    }
    else {
        // The bulge has collapsed, attempt to reintroduce using
        // 2-small-subdiagonals trick
        T vt_[3];
        auto vt = new_3_vector(vt_);
        auto H2 = slice(H, range{1, 4}, range{1, 4});
        lahqr_shiftcolumn(H2, vt, s1, s2);
        larfg(FORWARD, COLUMNWISE_STORAGE, vt, tau);
        vt[0] = tau;
 
        refsum = conj(vt[0]) * H(1, 0) + conj(vt[1]) * H(2, 0);
        if (abs1(H(2, 0) - refsum * vt[1]) + abs1(refsum * vt[2]) >
            eps * (abs1(H(0, 0)) + abs1(H(1, 1)) + abs1(H(2, 2)))) {
            // Starting a new bulge here would create non-negligible fill. Use
            // the old one.
            H(1, 0) = beta;
            H(2, 0) = zero;
            H(3, 0) = zero;
        }
        else {
            // Fill-in is negligible, use the new reflector.
            H(1, 0) = H(1, 0) - refsum;
            H(2, 0) = zero;
            H(3, 0) = zero;
            v[0] = vt[0];
            v[1] = vt[1];
            v[2] = vt[2];
        }
    }
}
 
}  // namespace tlapack
 
#endif  // TLAPACK_MOVE_BULGE_HH