Function hmatrixtd

src/ortfac.cpp:2613–2768

Reduction of a Hermitian matrix which is given by its higher or lower triangular part to a real tridiagonal matrix using unitary similarity transformation: Q'*A*Q = T. Input parameters: A - matrix to be transformed array with elements [0..N-1, 0..N-1]. N - size of matrix A. IsUpper - storage format. If IsUpper = True, then matrix A is given

     October 31, 1992
*************************************************************************/
void hmatrixtd(ap::complex_2d_array& a,
     int n,
     bool isupper,
     ap::complex_1d_array& tau,
     ap::real_1d_array& d,
     ap::real_1d_array& e)
{
    int i;
    ap::complex alpha;
    ap::complex taui;
    ap::complex v;
    ap::complex_1d_array t;
    ap::complex_1d_array t2;
    ap::complex_1d_array t3;

    if( n<=0 )
    {
        return;
    }
    for(i = 0; i <= n-1; i++)
    {
        ap::ap_error::make_assertion(ap::fp_eq(a(i,i).y,0), "");
    }
    if( n>1 )
    {
        tau.setbounds(0, n-2);
        e.setbounds(0, n-2);
    }
    d.setbounds(0, n-1);
    t.setbounds(0, n-1);
    t2.setbounds(0, n-1);
    t3.setbounds(0, n-1);
    if( isupper )
    {
        
        //
        // Reduce the upper triangle of A
        //
        a(n-1,n-1) = a(n-1,n-1).x;
        for(i = n-2; i >= 0; i--)
        {
            
            //
            // Generate elementary reflector H = I+1 - tau * v * v'
            //
            alpha = a(i,i+1);
            t(1) = alpha;
            if( i>=1 )
            {
                ap::vmove(&t(2), 1, &a(0, i+1), a.getstride(), "N", ap::vlen(2,i+1));
            }
            complexgeneratereflection(t, i+1, taui);
            if( i>=1 )
            {
                ap::vmove(&a(0, i+1), a.getstride(), &t(2), 1, "N", ap::vlen(0,i-1));
            }
            alpha = t(1);
            e(i) = alpha.x;