bidiagonal.cpp

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

/*************************************************************************
Unpacking matrix P which reduces matrix A to bidiagonal form.
The subroutine returns transposed matrix P.

Input parameters:
    QP      -   matrices Q and P in compact form.
                Output of ToBidiagonal subroutine.
    M       -   number of rows in matrix A.
    N       -   number of columns in matrix A.
    TAUP    -   scalar factors which are used to form P.
                Output of ToBidiagonal subroutine.
    PTRows  -   required number of rows of matrix P^T. N >= PTRows >= 0.

Output parameters:
    PT      -   first PTRows columns of matrix P^T
                Array[0..PTRows-1, 0..N-1]
                If PTRows=0, the array is not modified.

  -- ALGLIB --
     Copyright 2005-2007 by Bochkanov Sergey
*************************************************************************/
void rmatrixbdunpackpt(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;

    ap::ap_error::make_assertion(ptrows<=n, "RMatrixBDUnpackPT: PTRows>N!");
    ap::ap_error::make_assertion(ptrows>=0, "RMatrixBDUnpackPT: PTRows<0!");
    if( m==0||n==0||ptrows==0 )
    {
        return;
    }
    
    //
    // prepare PT
    //
    pt.setbounds(0, ptrows-1, 0, n-1);
    for(i = 0; i <= ptrows-1; i++)
    {
        for(j = 0; j <= n-1; j++)
        {
            if( i==j )
            {
                pt(i,j) = 1;
            }
            else
            {
                pt(i,j) = 0;
            }
        }
    }
    
    //
    // Calculate
    //
    rmatrixbdmultiplybyp(qp, m, n, taup, pt, ptrows, n, true, true);
}


/*************************************************************************
Multiplication by matrix P which reduces matrix A to  bidiagonal form.

The algorithm allows pre- or post-multiply by P or P'.

Input parameters:
    QP          -   matrices Q and P in compact form.
                    Output of RMatrixBD subroutine.
    M           -   number of rows in matrix A.
    N           -   number of columns in matrix A.
    TAUP        -   scalar factors which are used to form P.
                    Output of RMatrixBD subroutine.
    Z           -   multiplied matrix.
                    Array whose indexes range within [0..ZRows-1,0..ZColumns-1].
    ZRows       -   number of rows in matrix Z. If FromTheRight=False,
                    ZRows=N, otherwise ZRows can be arbitrary.
    ZColumns    -   number of columns in matrix Z. If FromTheRight=True,
                    ZColumns=N, otherwise ZColumns can be arbitrary.
    FromTheRight -  pre- or post-multiply.
    DoTranspose -   multiply by P or P'.

Output parameters:
    Z - product of Z and P.
                Array whose indexes range within [0..ZRows-1,0..ZColumns-1].
                If ZRows=0 or ZColumns=0, the array is not modified.

  -- ALGLIB --
     Copyright 2005-2007 by Bochkanov Sergey
*************************************************************************/
void rmatrixbdmultiplybyp(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;
    ap::real_1d_array v;
    ap::real_1d_array work;
    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, "RMatrixBDMultiplyByP: incorrect Z size!");
    
    //
    // init
    //
    mx = ap::maxint(m, n);
    mx = ap::maxint(mx, zrows);
    mx = ap::maxint(mx, zcolumns);
    v.setbounds(0, mx);
    work.setbounds(0, mx);
    v.setbounds(0, mx);
    work.setbounds(0, mx);
    if( m>=n )
    {
        
        //
        // setup
        //
        if( fromtheright )
        {
            i1 = n-2;
            i2 = 0;
            istep = -1;
        }
        else
        {
            i1 = 0;
            i2 = n-2;
            istep = +1;
        }
        if( !dotranspose )
        {
            i = i1;
            i1 = i2;
            i2 = i;
            istep = -istep;
        }
        
        //
        // Process
        //
        if( n-1>0 )
        {
            i = i1;
            do
            {
                ap::vmove(&v(1), &qp(i, i+1), ap::vlen(1,n-1-i));
                v(1) = 1;
                if( fromtheright )
                {
                    applyreflectionfromtheright(z, taup(i), v, 0, zrows-1, i+1, n-1, work);
                }
                else
                {
                    applyreflectionfromtheleft(z, taup(i), v, i+1, n-1, 0, zcolumns-1, work);
                }
                i = i+istep;
            }
            while(i!=i2+istep);
        }
    }
    else
    {
        
        //
        // setup
        //
        if( fromtheright )
        {
            i1 = m-1;
            i2 = 0;
            istep = -1;
        }
        else
        {
            i1 = 0;
            i2 = m-1;
            istep = +1;
        }
        if( !dotranspose )
        {
            i = i1;
            i1 = i2;
            i2 = i;
            istep = -istep;
        }
        
        //
        // Process
        //
        i = i1;
        do
        {
            ap::vmove(&v(1), &qp(i, i), ap::vlen(1,n-i));
            v(1) = 1;
            if( fromtheright )
            {
                applyreflectionfromtheright(z, taup(i), v, 0, zrows-1, i, n-1, work);
            }
            else
            {
                applyreflectionfromtheleft(z, taup(i), v, i, n-1, 0, zcolumns-1, work);
            }
            i = i+istep;
        }
        while(i!=i2+istep);
    }
}


/*************************************************************************
Unpacking of the main and secondary diagonals of bidiagonal decomposition
of matrix A.

Input parameters:
    B   -   output of RMatrixBD subroutine.
    M   -   number of rows in matrix B.
    N   -   number of columns in matrix B.

Output parameters:
    IsUpper -   True, if the matrix is upper bidiagonal.
                otherwise IsUpper is False.
    D       -   the main diagonal.
                Array whose index ranges within [0..Min(M,N)-1].
    E       -   the secondary diagonal (upper or lower, depending on
                the value of IsUpper).
                Array index ranges within [0..Min(M,N)-1], the last
                element is not used.

  -- ALGLIB --
     Copyright 2005-2007 by Bochkanov Sergey
*************************************************************************/
void rmatrixbdunpackdiagonals(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(0, n-1);
        e.setbounds(0, n-1);
        for(i = 0; i <= n-2; i++)
        {
            d(i) = b(i,i);
            e(i) = b(i,i+1);
        }
        d(n-1) = b(n-1,n-1);
    }
    else
    {
        d.setbounds(0, m-1);
        e.setbounds(0, m-1);
        for(i = 0; i <= m-2; i++)
        {
            d(i) = b(i,i);
            e(i) = b(i+1,i);
        }
        d(m-1) = b(m-1,m-1);
    }
}


/*************************************************************************
Obsolete 1-based subroutine.
See RMatrixBD for 0-based replacement.
*************************************************************************/
void tobidiagonal(ap::real_2d_array& a,
     int m,
     int n,
     ap::real_1d_array& tauq,
     ap::real_1d_array& taup)
{
    ap::real_1d_array work;
    ap::real_1d_array t;
    int minmn;
    int maxmn;
    int i;
    double ltau;
    int mmip1;
    int nmi;
    int ip1;
    int nmip1;
    int mmi;

    minmn = ap::minint(m, n);
    maxmn = ap::maxint(m, n);
    work.setbounds(1, maxmn);
    t.setbounds(1, maxmn);
    taup.setbounds(1, minmn);
    tauq.setbounds(1, minmn);
    if( m>=n )
    {
        
        //
        // Reduce to upper bidiagonal form
        //
        for(i = 1; i <= n; i++)
        {
            
            //
            // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
            //
            mmip1 = m-i+1;
            ap::vmove(t.getvector(1, mmip1), a.getcolumn(i, i, m));
            generatereflection(t, mmip1, ltau);
            tauq(i) = ltau;
            ap::vmove(a.getcolumn(i, i, m), t.getvector(1, mmip1));
            t(1) = 1;
            
            //
            // Apply H(i) to A(i:m,i+1:n) from the left
            //
            applyreflectionfromtheleft(a, ltau, t, i, m, i+1, n, work);
            if( i<n )
            {
                
                //
                // Generate elementary reflector G(i) to annihilate
                // A(i,i+2:n)
                //
                nmi = n-i;
                ip1 = i+1;
                ap::vmove(&t(1), &a(i, ip1), ap::vlen(1,nmi));
                generatereflection(t, nmi, ltau);
                taup(i) = ltau;
                ap::vmove(&a(i, ip1), &t(1), ap::vlen(ip1,n));
                t(1) = 1;
                
                //
                // Apply G(i) to A(i+1:m,i+1:n) from the right
                //
                applyreflectionfromtheright(a, ltau, t, i+1, m, i+1, n, work);
            }
            else
            {
                taup(i) = 0;
            }
        }
    }
    else
    {
        
        //
        // Reduce to lower bidiagonal form
        //
        for(i = 1; i <= m; i++)
        {
            
            //
            // Generate elementary reflector G(i) to annihilate A(i,i+1:n)
            //
            nmip1 = n-i+1;
            ap::vmove(&t(1), &a(i, i), ap::vlen(1,nmip1));
            generatereflection(t, nmip1, ltau);
            taup(i) = ltau;
            ap::vmove(&a(i, i), &t(1), ap::vlen(i,n));
            t(1) = 1;
            
            //
            // Apply G(i) to A(i+1:m,i:n) from the right
            //
            applyreflectionfromtheright(a, ltau, t, i+1, m, i, n, work);
            if( i<m )
            {
                
                //
                // Generate elementary reflector H(i) to annihilate
                // A(i+2:m,i)
                //
                mmi = m-i;
                ip1 = i+1;
                ap::vmove(t.getvector(1, mmi), a.getcolumn(i, ip1, m));
                generatereflection(t, mmi, ltau);
                tauq(i) = ltau;
                ap::vmove(a.getcolumn(i, ip1, m), t.getvector(1, mmi));
                t(1) = 1;
                
                //
                // Apply H(i) to A(i+1:m,i+1:n) from the left
                //
                applyreflectionfromtheleft(a, ltau, t, i+1, m, i+1, n, work);
            }
            else
            {
                tauq(i) = 0;
            }
        }
    }
}


/*************************************************************************
Obsolete 1-based subroutine.
See RMatrixBDUnpackQ for 0-based replacement.
*************************************************************************/
void unpackqfrombidiagonal(const ap::real_2d_array& qp,
     int m,
     int n,
     const ap::real_1d_array& tauq,
     int qcolumns,
     ap::real_2d_array& q)
{
    int i;
    int j;
    int ip1;
    ap::real_1d_array v;
    ap::real_1d_array work;
    int vm;

    ap::ap_error::make_assertion(qcolumns<=m, "UnpackQFromBidiagonal: QColumns>M!");
    if( m==0||n==0||qcolumns==0 )
    {
        return;
    }
    
    //
    // init