bidiagonal.cpp

The document describes an AP library adapted for C++. The AP library for C++ contains a basic set of

    //
    q.setbounds(1, m, 1, qcolumns);
    v.setbounds(1, m);
    work.setbounds(1, qcolumns);
    
    //
    // prepare Q
    //
    for(i = 1; i <= m; i++)
    {
        for(j = 1; j <= qcolumns; j++)
        {
            if( i==j )
            {
                q(i,j) = 1;
            }
            else
            {
                q(i,j) = 0;
            }
        }
    }
    if( m>=n )
    {
        for(i = ap::minint(n, qcolumns); i >= 1; i--)
        {
            vm = m-i+1;
            ap::vmove(v.getvector(1, vm), qp.getcolumn(i, i, m));
            v(1) = 1;
            applyreflectionfromtheleft(q, tauq(i), v, i, m, 1, qcolumns, work);
        }
    }
    else
    {
        for(i = ap::minint(m-1, qcolumns-1); i >= 1; i--)
        {
            vm = m-i;
            ip1 = i+1;
            ap::vmove(v.getvector(1, vm), qp.getcolumn(i, ip1, m));
            v(1) = 1;
            applyreflectionfromtheleft(q, tauq(i), v, i+1, m, 1, qcolumns, work);
        }
    }
}


/*************************************************************************
Obsolete 1-based subroutine.
See RMatrixBDMultiplyByQ for 0-based replacement.
*************************************************************************/
void multiplybyqfrombidiagonal(const ap::real_2d_array& qp,
     int m,
     int n,
     const ap::real_1d_array& tauq,
     ap::real_2d_array& z,
     int zrows,
     int zcolumns,
     bool fromtheright,
     bool dotranspose)
{
    int i;
    int ip1;
    int i1;
    int i2;
    int istep;
    ap::real_1d_array v;
    ap::real_1d_array work;
    int vm;
    int mx;

    if( m<=0||n<=0||zrows<=0||zcolumns<=0 )
    {
        return;
    }
    ap::ap_error::make_assertion(fromtheright&&zcolumns==m||!fromtheright&&zrows==m, "MultiplyByQFromBidiagonal: incorrect Z size!");
    
    //
    // init
    //
    mx = ap::maxint(m, n);
    mx = ap::maxint(mx, zrows);
    mx = ap::maxint(mx, zcolumns);
    v.setbounds(1, mx);
    work.setbounds(1, mx);
    if( m>=n )
    {
        
        //
        // setup
        //
        if( fromtheright )
        {
            i1 = 1;
            i2 = n;
            istep = +1;
        }
        else
        {
            i1 = n;
            i2 = 1;
            istep = -1;
        }
        if( dotranspose )
        {
            i = i1;
            i1 = i2;
            i2 = i;
            istep = -istep;
        }
        
        //
        // Process
        //
        i = i1;
        do
        {
            vm = m-i+1;
            ap::vmove(v.getvector(1, vm), qp.getcolumn(i, i, m));
            v(1) = 1;
            if( fromtheright )
            {
                applyreflectionfromtheright(z, tauq(i), v, 1, zrows, i, m, work);
            }
            else
            {
                applyreflectionfromtheleft(z, tauq(i), v, i, m, 1, zcolumns, work);
            }
            i = i+istep;
        }
        while(i!=i2+istep);
    }
    else
    {
        
        //
        // setup
        //
        if( fromtheright )
        {
            i1 = 1;
            i2 = m-1;
            istep = +1;
        }
        else
        {
            i1 = m-1;
            i2 = 1;
            istep = -1;
        }
        if( dotranspose )
        {
            i = i1;
            i1 = i2;
            i2 = i;
            istep = -istep;
        }
        
        //
        // Process
        //
        if( m-1>0 )
        {
            i = i1;
            do
            {
                vm = m-i;
                ip1 = i+1;
                ap::vmove(v.getvector(1, vm), qp.getcolumn(i, ip1, m));
                v(1) = 1;
                if( fromtheright )
                {
                    applyreflectionfromtheright(z, tauq(i), v, 1, zrows, i+1, m, work);
                }
                else
                {
                    applyreflectionfromtheleft(z, tauq(i), v, i+1, m, 1, zcolumns, work);
                }
                i = i+istep;
            }
            while(i!=i2+istep);
        }
    }
}


/*************************************************************************
Obsolete 1-based subroutine.
See RMatrixBDUnpackPT for 0-based replacement.
*************************************************************************/
void unpackptfrombidiagonal(const ap::real_2d_array& qp,
     int m,
     int n,
     const ap::real_1d_array& taup,
     int ptrows,
     ap::real_2d_array& pt)
{
    int i;
    int j;
    int ip1;
    ap::real_1d_array v;
    ap::real_1d_array work;
    int vm;

    ap::ap_error::make_assertion(ptrows<=n, "UnpackPTFromBidiagonal: PTRows>N!");
    if( m==0||n==0||ptrows==0 )
    {
        return;
    }
    
    //
    // init
    //
    pt.setbounds(1, ptrows, 1, n);
    v.setbounds(1, n);
    work.setbounds(1, ptrows);
    
    //
    // prepare PT
    //
    for(i = 1; i <= ptrows; i++)
    {
        for(j = 1; j <= n; j++)
        {
            if( i==j )
            {
                pt(i,j) = 1;
            }
            else
            {
                pt(i,j) = 0;
            }
        }
    }
    if( m>=n )
    {
        for(i = ap::minint(n-1, ptrows-1); i >= 1; i--)
        {
            vm = n-i;
            ip1 = i+1;
            ap::vmove(&v(1), &qp(i, ip1), ap::vlen(1,vm));
            v(1) = 1;
            applyreflectionfromtheright(pt, taup(i), v, 1, ptrows, i+1, n, work);
        }
    }
    else
    {
        for(i = ap::minint(m, ptrows); i >= 1; i--)
        {
            vm = n-i+1;
            ap::vmove(&v(1), &qp(i, i), ap::vlen(1,vm));
            v(1) = 1;
            applyreflectionfromtheright(pt, taup(i), v, 1, ptrows, i, n, work);
        }
    }
}


/*************************************************************************
Obsolete 1-based subroutine.
See RMatrixBDMultiplyByP for 0-based replacement.
*************************************************************************/
void multiplybypfrombidiagonal(const ap::real_2d_array& qp,
     int m,
     int n,
     const ap::real_1d_array& taup,
     ap::real_2d_array& z,
     int zrows,
     int zcolumns,
     bool fromtheright,
     bool dotranspose)
{
    int i;
    int ip1;
    ap::real_1d_array v;
    ap::real_1d_array work;
    int vm;
    int mx;
    int i1;
    int i2;
    int istep;

    if( m<=0||n<=0||zrows<=0||zcolumns<=0 )
    {
        return;
    }
    ap::ap_error::make_assertion(fromtheright&&zcolumns==n||!fromtheright&&zrows==n, "MultiplyByQFromBidiagonal: incorrect Z size!");
    
    //
    // init
    //
    mx = ap::maxint(m, n);
    mx = ap::maxint(mx, zrows);
    mx = ap::maxint(mx, zcolumns);
    v.setbounds(1, mx);
    work.setbounds(1, mx);
    v.setbounds(1, mx);
    work.setbounds(1, mx);
    if( m>=n )
    {
        
        //
        // setup
        //
        if( fromtheright )
        {
            i1 = n-1;
            i2 = 1;
            istep = -1;
        }
        else
        {
            i1 = 1;
            i2 = n-1;
            istep = +1;
        }
        if( !dotranspose )
        {
            i = i1;
            i1 = i2;
            i2 = i;
            istep = -istep;
        }
        
        //
        // Process
        //
        if( n-1>0 )
        {
            i = i1;
            do
            {
                vm = n-i;
                ip1 = i+1;
                ap::vmove(&v(1), &qp(i, ip1), ap::vlen(1,vm));
                v(1) = 1;
                if( fromtheright )
                {
                    applyreflectionfromtheright(z, taup(i), v, 1, zrows, i+1, n, work);
                }
                else
                {
                    applyreflectionfromtheleft(z, taup(i), v, i+1, n, 1, zcolumns, work);
                }
                i = i+istep;
            }
            while(i!=i2+istep);
        }
    }
    else
    {
        
        //
        // setup
        //
        if( fromtheright )
        {
            i1 = m;
            i2 = 1;
            istep = -1;
        }
        else
        {
            i1 = 1;
            i2 = m;
            istep = +1;
        }
        if( !dotranspose )
        {
            i = i1;
            i1 = i2;
            i2 = i;
            istep = -istep;
        }
        
        //
        // Process
        //
        i = i1;
        do
        {
            vm = n-i+1;
            ap::vmove(&v(1), &qp(i, i), ap::vlen(1,vm));
            v(1) = 1;
            if( fromtheright )
            {
                applyreflectionfromtheright(z, taup(i), v, 1, zrows, i, n, work);
            }
            else
            {
                applyreflectionfromtheleft(z, taup(i), v, i, n, 1, zcolumns, work);
            }
            i = i+istep;
        }
        while(i!=i2+istep);
    }
}


/*************************************************************************
Obsolete 1-based subroutine.
See RMatrixBDUnpackDiagonals for 0-based replacement.
*************************************************************************/
void unpackdiagonalsfrombidiagonal(const ap::real_2d_array& b,
     int m,
     int n,
     bool& isupper,
     ap::real_1d_array& d,
     ap::real_1d_array& e)
{
    int i;

    isupper = m>=n;
    if( m==0||n==0 )
    {
        return;
    }
    if( isupper )
    {
        d.setbounds(1, n);
        e.setbounds(1, n);
        for(i = 1; i <= n-1; i++)
        {
            d(i) = b(i,i);
            e(i) = b(i,i+1);
        }
        d(n) = b(n,n);
    }
    else
    {
        d.setbounds(1, m);
        e.setbounds(1, m);
        for(i = 1; i <= m-1; i++)
        {
            d(i) = b(i,i);
            e(i) = b(i+1,i);
        }
        d(m) = b(m,m);
    }
}