C++ real_2d_array::getstride方法代码示例

/*************************************************************************
Example of usage of an IterativeEstimateNorm subroutine

Input parameters:
    A   -   matrix.
            Array whose indexes range within [1..N, 1..N].

Return:
    Matrix norm estimated by the subroutine.

  -- ALGLIB --
     Copyright 2005 by Bochkanov Sergey
*************************************************************************/
double demoiterativeestimate1norm(const ap::real_2d_array& a, int n)
{
    double result;
    int i;
    double s;
    ap::real_1d_array x;
    ap::real_1d_array t;
    ap::real_1d_array v;
    ap::integer_1d_array iv;
    int kase;

    kase = 0;
    t.setbounds(1, n);
    iterativeestimate1norm(n, v, x, iv, result, kase);
    while(kase!=0)
    {
        if( kase==1 )
        {
            for(i = 1; i <= n; i++)
            {
                s = ap::vdotproduct(&a(i, 1), 1, &x(1), 1, ap::vlen(1,n));
                t(i) = s;
            }
        }
        else
        {
            for(i = 1; i <= n; i++)
            {
                s = ap::vdotproduct(&a(1, i), a.getstride(), &x(1), 1, ap::vlen(1,n));
                t(i) = s;
            }
        }
        ap::vmove(&x(1), 1, &t(1), 1, ap::vlen(1,n));
        iterativeestimate1norm(n, v, x, iv, result, kase);
    }
    return result;
}

/*************************************************************************
LU inverse
*************************************************************************/
static bool rmatrixinvmatlu(ap::real_2d_array& a,
     const ap::integer_1d_array& pivots,
     int n)
{
    bool result;
    ap::real_1d_array work;
    int i;
    int iws;
    int j;
    int jb;
    int jj;
    int jp;
    double v;

    result = true;
    
    //
    // Quick return if possible
    //
    if( n==0 )
    {
        return result;
    }
    work.setbounds(0, n-1);
    
    //
    // Form inv(U)
    //
    if( !rmatrixinvmattr(a, n, true, false) )
    {
        result = false;
        return result;
    }
    
    //
    // Solve the equation inv(A)*L = inv(U) for inv(A).
    //
    for(j = n-1; j >= 0; j--)
    {
        
        //
        // Copy current column of L to WORK and replace with zeros.
        //
        for(i = j+1; i <= n-1; i++)
        {
            work(i) = a(i,j);
            a(i,j) = 0;
        }
        
        //
        // Compute current column of inv(A).
        //
        if( j<n-1 )
        {
            for(i = 0; i <= n-1; i++)
            {
                v = ap::vdotproduct(&a(i, j+1), 1, &work(j+1), 1, ap::vlen(j+1,n-1));
                a(i,j) = a(i,j)-v;
            }
        }
    }
    
    //
    // Apply column interchanges.
    //
    for(j = n-2; j >= 0; j--)
    {
        jp = pivots(j);
        if( jp!=j )
        {
            ap::vmove(&work(0), 1, &a(0, j), a.getstride(), ap::vlen(0,n-1));
            ap::vmove(&a(0, j), a.getstride(), &a(0, jp), a.getstride(), ap::vlen(0,n-1));
            ap::vmove(&a(0, jp), a.getstride(), &work(0), 1, ap::vlen(0,n-1));
        }
    }
    return result;
}

static void getbdsvderror(const ap::real_1d_array& d,
     const ap::real_1d_array& e,
     int n,
     bool isupper,
     const ap::real_2d_array& u,
     const ap::real_2d_array& c,
     const ap::real_1d_array& w,
     const ap::real_2d_array& vt,
     double& materr,
     double& orterr,
     bool& wsorted)
{
    int i;
    int j;
    int k;
    double locerr;
    double sm;

    
    //
    // decomposition error
    //
    locerr = 0;
    for(i = 0; i <= n-1; i++)
    {
        for(j = 0; j <= n-1; j++)
        {
            sm = 0;
            for(k = 0; k <= n-1; k++)
            {
                sm = sm+w(k)*u(i,k)*vt(k,j);
            }
            if( isupper )
            {
                if( i==j )
                {
                    locerr = ap::maxreal(locerr, fabs(d(i)-sm));
                }
                else
                {
                    if( i==j-1 )
                    {
                        locerr = ap::maxreal(locerr, fabs(e(i)-sm));
                    }
                    else
                    {
                        locerr = ap::maxreal(locerr, fabs(sm));
                    }
                }
            }
            else
            {
                if( i==j )
                {
                    locerr = ap::maxreal(locerr, fabs(d(i)-sm));
                }
                else
                {
                    if( i-1==j )
                    {
                        locerr = ap::maxreal(locerr, fabs(e(j)-sm));
                    }
                    else
                    {
                        locerr = ap::maxreal(locerr, fabs(sm));
                    }
                }
            }
        }
    }
    materr = ap::maxreal(materr, locerr);
    
    //
    // check for C = U'
    // we consider it as decomposition error
    //
    locerr = 0;
    for(i = 0; i <= n-1; i++)
    {
        for(j = 0; j <= n-1; j++)
        {
            locerr = ap::maxreal(locerr, fabs(u(i,j)-c(j,i)));
        }
    }
    materr = ap::maxreal(materr, locerr);
    
    //
    // orthogonality error
    //
    locerr = 0;
    for(i = 0; i <= n-1; i++)
    {
        for(j = i; j <= n-1; j++)
        {
            sm = ap::vdotproduct(&u(0, i), u.getstride(), &u(0, j), u.getstride(), ap::vlen(0,n-1));
            if( i!=j )
            {
                locerr = ap::maxreal(locerr, fabs(sm));
            }
            else
//.........这里部分代码省略.........

/********************************************************************
real TRSM kernel
********************************************************************/
bool ialglib::_i_rmatrixrighttrsmf(int m,
     int n,
     const ap::real_2d_array& a,
     int i1,
     int j1,
     bool isupper,
     bool isunit,
     int optype,
     ap::real_2d_array& x,
     int i2,
     int j2)
{
    if( m>alglib_r_block || n>alglib_r_block )
        return false;

    //
    // local buffers
    //
    double *pdiag;
    int i;
    double __abuf[alglib_r_block*alglib_r_block+alglib_simd_alignment];
    double __xbuf[alglib_r_block*alglib_r_block+alglib_simd_alignment];
    double __tmpbuf[alglib_r_block+alglib_simd_alignment];
    double * const abuf   = (double * const) alglib_align(__abuf,  alglib_simd_alignment);
    double * const xbuf   = (double * const) alglib_align(__xbuf,  alglib_simd_alignment);
    double * const tmpbuf = (double * const) alglib_align(__tmpbuf,alglib_simd_alignment);

    //
    // Prepare
    //
    bool uppera;
    mcopyblock(n, n, &a(i1,j1), optype, a.getstride(), abuf);
    mcopyblock(m, n, &x(i2,j2), 0, x.getstride(), xbuf);
    if( isunit )
        for(i=0,pdiag=abuf; i<n; i++,pdiag+=alglib_r_block+1)
            *pdiag = 1.0;
    if( optype==0 )
        uppera = isupper;
    else
        uppera = !isupper;

    //
    // Solve Y*A^-1=X where A is upper or lower triangular
    //
    if( uppera )
    {
        for(i=0,pdiag=abuf; i<n; i++,pdiag+=alglib_r_block+1)
        {
            double beta  = 1.0/(*pdiag);
            double alpha = -beta;
            vcopy(i, abuf+i, alglib_r_block, tmpbuf, 1);
            mv(m, i, xbuf, tmpbuf, xbuf+i, alglib_r_block, alpha, beta);
        }
        mcopyunblock(m, n, xbuf, 0, &x(i2,j2), x.getstride());
    }
    else
    {
        for(i=n-1,pdiag=abuf+(n-1)*alglib_r_block+(n-1); i>=0; i--,pdiag-=alglib_r_block+1)
        {
            double beta = 1.0/(*pdiag);
            double alpha = -beta;
            vcopy(n-1-i, pdiag+alglib_r_block, alglib_r_block, tmpbuf, 1);
            mv(m, n-1-i, xbuf+i+1, tmpbuf, xbuf+i, alglib_r_block, alpha, beta);
        }
        mcopyunblock(m, n, xbuf, 0, &x(i2,j2), x.getstride());
    }
    return true;
}

/*************************************************************************
triangular inverse
*************************************************************************/
static bool rmatrixinvmattr(ap::real_2d_array& a,
     int n,
     bool isupper,
     bool isunittriangular)
{
    bool result;
    bool nounit;
    int i;
    int j;
    double v;
    double ajj;
    ap::real_1d_array t;

    result = true;
    t.setbounds(0, n-1);
    
    //
    // Test the input parameters.
    //
    nounit = !isunittriangular;
    if( isupper )
    {
        
        //
        // Compute inverse of upper triangular matrix.
        //
        for(j = 0; j <= n-1; j++)
        {
            if( nounit )
            {
                if( ap::fp_eq(a(j,j),0) )
                {
                    result = false;
                    return result;
                }
                a(j,j) = 1/a(j,j);
                ajj = -a(j,j);
            }
            else
            {
                ajj = -1;
            }
            
            //
            // Compute elements 1:j-1 of j-th column.
            //
            if( j>0 )
            {
                ap::vmove(&t(0), 1, &a(0, j), a.getstride(), ap::vlen(0,j-1));
                for(i = 0; i <= j-1; i++)
                {
                    if( i<j-1 )
                    {
                        v = ap::vdotproduct(&a(i, i+1), 1, &t(i+1), 1, ap::vlen(i+1,j-1));
                    }
                    else
                    {
                        v = 0;
                    }
                    if( nounit )
                    {
                        a(i,j) = v+a(i,i)*t(i);
                    }
                    else
                    {
                        a(i,j) = v+t(i);
                    }
                }
                ap::vmul(&a(0, j), a.getstride(), ap::vlen(0,j-1), ajj);
            }
        }
    }
    else
    {
        
        //
        // Compute inverse of lower triangular matrix.
        //
        for(j = n-1; j >= 0; j--)
        {
            if( nounit )
            {
                if( ap::fp_eq(a(j,j),0) )
                {
                    result = false;
                    return result;
                }
                a(j,j) = 1/a(j,j);
                ajj = -a(j,j);
            }
            else
            {
                ajj = -1;
            }
            if( j<n-1 )
            {
                
//.........这里部分代码省略.........

/********************************************************************
real rank-1 kernel
********************************************************************/
bool ialglib::_i_rmatrixrank1f(int m,
     int n,
     ap::real_2d_array& a,
     int ia,
     int ja,
     ap::real_1d_array& u,
     int uoffs,
     ap::real_1d_array& v,
     int voffs)
{
    double *arow0, *arow1, *pu, *pv, *vtmp, *dst0, *dst1;
    int m2 = m/2;
    int n2 = n/2;
    int stride  = a.getstride();
    int stride2 = 2*a.getstride();
    int i, j;

    //
    // update pairs of rows
    //
    arow0 = &a(ia,ja);
    arow1 = arow0+stride;
    pu    = &u(uoffs);
    vtmp  = &v(voffs);
    for(i=0; i<m2; i++,arow0+=stride2,arow1+=stride2,pu+=2)
    {
        //
        // update by two
        //
        for(j=0,pv=vtmp, dst0=arow0, dst1=arow1; j<n2; j++, dst0+=2, dst1+=2, pv+=2)
        {
            dst0[0] += pu[0]*pv[0];
            dst0[1] += pu[0]*pv[1];
            dst1[0] += pu[1]*pv[0];
            dst1[1] += pu[1]*pv[1];
        }

        //
        // final update
        //
        if( n%2!=0 )
        {
            dst0[0] += pu[0]*pv[0];
            dst1[0] += pu[1]*pv[0];
        }
    }

    //
    // update last row
    //
    if( m%2!=0 )
    {
        //
        // update by two
        //
        for(j=0,pv=vtmp, dst0=arow0; j<n2; j++, dst0+=2, pv+=2)
        {
            dst0[0] += pu[0]*pv[0];
            dst0[1] += pu[0]*pv[1];
        }

        //
        // final update
        //
        if( n%2!=0 )
            dst0[0] += pu[0]*pv[0];
    }
    return true;
}

/********************************************************************
This is real GEMM kernel
********************************************************************/
bool ialglib::_i_rmatrixgemmf(int m,
     int n,
     int k,
     double alpha,
     const ap::real_2d_array& _a,
     int ia,
     int ja,
     int optypea,
     const ap::real_2d_array& _b,
     int ib,
     int jb,
     int optypeb,
     double beta,
     ap::real_2d_array& _c,
     int ic,
     int jc)
{
    if( m>alglib_r_block || n>alglib_r_block || k>alglib_r_block )
        return false;

    int i, stride, cstride;
    double *crow;
    double __abuf[alglib_r_block+alglib_simd_alignment];
    double __b[alglib_r_block*alglib_r_block+alglib_simd_alignment];
    double * const abuf = (double * const) alglib_align(__abuf,alglib_simd_alignment);
    double * const b    = (double * const) alglib_align(__b,   alglib_simd_alignment);

    //
    // copy b
    //
    if( optypeb==0 )
        mcopyblock(k, n, &_b(ib,jb), 1, _b.getstride(), b);
    else
        mcopyblock(n, k, &_b(ib,jb), 0, _b.getstride(), b);

    //
    // multiply B by A (from the right, by rows)
    // and store result in C
    //
    crow  = &_c(ic,jc);
    stride = _a.getstride();
    cstride = _c.getstride();
    if( optypea==0 )
    {
        const double *arow = &_a(ia,ja);
        for(i=0; i<m; i++)
        {
            vcopy(k, arow, 1, abuf, 1);
            if( beta==0 )
                vzero(n, crow, 1);
            mv(n, k, b, abuf, crow, 1, alpha, beta);
            crow += cstride;
            arow += stride;
        }
    }
    else
    {
        const double *acol = &_a(ia,ja);
        for(i=0; i<m; i++)
        {
            vcopy(k, acol, stride, abuf, 1);
            if( beta==0 )
                vzero(n, crow, 1);
            mv(n, k, b, abuf, crow, 1, alpha, beta);
            crow += cstride;
            acol++;
        }
    }
    return true;
}

/********************************************************************
real SYRK kernel
********************************************************************/
bool ialglib::_i_rmatrixsyrkf(int n,
     int k,
     double alpha,
     const ap::real_2d_array& a,
     int ia,
     int ja,
     int optypea,
     double beta,
     ap::real_2d_array& c,
     int ic,
     int jc,
     bool isupper)
{
    if( n>alglib_r_block || k>alglib_r_block )
        return false;
    if( n==0 )
        return true;
    
    //
    // local buffers
    //
    double *arow, *crow;
    int i;
    double __abuf[alglib_r_block*alglib_r_block+alglib_simd_alignment];
    double __cbuf[alglib_r_block*alglib_r_block+alglib_simd_alignment];
    double * const abuf   = (double * const) alglib_align(__abuf,  alglib_simd_alignment);
    double * const cbuf   = (double * const) alglib_align(__cbuf,  alglib_simd_alignment);

    //
    // copy A and C, task is transformed to "A*A^T"-form.
    // if beta==0, then C is filled by zeros (and not referenced)
    //
    // alpha==0 or k==0 are correctly processed (A is not referenced)
    //
    if( alpha==0 )
        k = 0;
    if( k>0 )
    {
        if( optypea==0 )
            mcopyblock(n, k, &a(ia,ja), 0, a.getstride(), abuf);
        else
            mcopyblock(k, n, &a(ia,ja), 1, a.getstride(), abuf);
    }
    mcopyblock(n, n, &c(ic,jc), 0, c.getstride(), cbuf);
    if( beta==0 )
    {
        for(i=0,crow=cbuf; i<n; i++,crow+=alglib_r_block)
            if( isupper )
                vzero(n-i, crow+i, 1);
            else
                vzero(i+1, crow, 1);
    }


    //
    // update C
    //
    if( isupper )
    {
        for(i=0,arow=abuf,crow=cbuf; i<n; i++,arow+=alglib_r_block,crow+=alglib_r_block)
        {
            mv(n-i, k, arow, arow, crow+i, 1, alpha, beta);
        }
    }
    else
    {
        for(i=0,arow=abuf,crow=cbuf; i<n; i++,arow+=alglib_r_block,crow+=alglib_r_block)
        {
            mv(i+1, k, abuf, arow, crow, 1, alpha, beta);
        }
    }

    //
    // copy back
    //
    mcopyunblock(n, n, cbuf, 0, &c(ic,jc), c.getstride());

    return true;
}

/********************************************************************
real TRSM kernel
********************************************************************/
bool ialglib::_i_rmatrixlefttrsmf(int m,
     int n,
     const ap::real_2d_array& a,
     int i1,
     int j1,
     bool isupper,
     bool isunit,
     int optype,
     ap::real_2d_array& x,
     int i2,
     int j2)
{
    if( m>alglib_r_block || n>alglib_r_block )
        return false;
    
    //
    // local buffers
    //
    double *pdiag, *arow;
    int i;
    double __abuf[alglib_r_block*alglib_r_block+alglib_simd_alignment];
    double __xbuf[alglib_r_block*alglib_r_block+alglib_simd_alignment];
    double __tmpbuf[alglib_r_block+alglib_simd_alignment];
    double * const abuf   = (double * const) alglib_align(__abuf,  alglib_simd_alignment);
    double * const xbuf   = (double * const) alglib_align(__xbuf,  alglib_simd_alignment);
    double * const tmpbuf = (double * const) alglib_align(__tmpbuf,alglib_simd_alignment);

    //
    // Prepare
    // Transpose X (so we may use mv, which calculates A*x, but not x*A)
    //
    bool uppera;
    mcopyblock(m, m, &a(i1,j1), optype, a.getstride(), abuf);
    mcopyblock(m, n, &x(i2,j2), 1, x.getstride(), xbuf);
    if( isunit )
        for(i=0,pdiag=abuf; i<m; i++,pdiag+=alglib_r_block+1)
            *pdiag = 1.0;
    if( optype==0 )
        uppera = isupper;
    else
        uppera = !isupper;

    //
    // Solve A^-1*Y^T=X^T where A is upper or lower triangular
    //
    if( uppera )
    {
        for(i=m-1,pdiag=abuf+(m-1)*alglib_r_block+(m-1); i>=0; i--,pdiag-=alglib_r_block+1)
        {
            double beta = 1.0/(*pdiag);
            double alpha = -beta;
            vcopy(m-1-i, pdiag+1, 1, tmpbuf, 1);
            mv(n, m-1-i, xbuf+i+1, tmpbuf, xbuf+i, alglib_r_block, alpha, beta);
        }
        mcopyunblock(m, n, xbuf, 1, &x(i2,j2), x.getstride());
    }
    else
    {   for(i=0,pdiag=abuf,arow=abuf; i<m; i++,pdiag+=alglib_r_block+1,arow+=alglib_r_block)
        {
            double beta = 1.0/(*pdiag);
            double alpha = -beta;
            vcopy(i, arow, 1, tmpbuf, 1);
            mv(n, i, xbuf, tmpbuf, xbuf+i, alglib_r_block, alpha, beta);
        }
        mcopyunblock(m, n, xbuf, 1, &x(i2,j2), x.getstride());
    }
    return true;
}

//.........这里部分代码省略.........
            return;
        }
        
        //
        // Special case: degenerate ST matrix, multicollinearity found.
        // Since we know ST eigenvalues/vectors we can translate task to
        // non-degenerate form.
        //
        // Let WG is orthogonal basis of the non zero variance subspace
        // of the ST and let WZ is orthogonal basis of the zero variance
        // subspace.
        //
        // Projection on WG allows us to use LDA on reduced M-dimensional
        // subspace, N-M vectors of WZ allows us to update reduced LDA
        // factors to full N-dimensional subspace.
        //
        m = 0;
        for(k = 0; k <= nvars-1; k++)
        {
            if( ap::fp_less_eq(d(k),1000*ap::machineepsilon*d(nvars-1)) )
            {
                m = k+1;
            }
        }
        ap::ap_error::make_assertion(m!=0, "FisherLDAN: internal error #1");
        xyproj.setbounds(0, npoints-1, 0, nvars-m);
        matrixmatrixmultiply(xy, 0, npoints-1, 0, nvars-1, false, z, 0, nvars-1, m, nvars-1, false, 1.0, xyproj, 0, npoints-1, 0, nvars-m-1, 0.0, work);
        for(i = 0; i <= npoints-1; i++)
        {
            xyproj(i,nvars-m) = xy(i,nvars);
        }
        fisherldan(xyproj, npoints, nvars-m, nclasses, info, wproj);
        if( info<0 )
        {
            return;
        }
        matrixmatrixmultiply(z, 0, nvars-1, m, nvars-1, false, wproj, 0, nvars-m-1, 0, nvars-m-1, false, 1.0, w, 0, nvars-1, 0, nvars-m-1, 0.0, work);
        for(k = nvars-m; k <= nvars-1; k++)
        {
            ap::vmove(&w(0, k), w.getstride(), &z(0, k-(nvars-m)), z.getstride(), ap::vlen(0,nvars-1));
        }
        info = 2;
    }
    else
    {
        
        //
        // General case: no multicollinearity
        //
        tm.setbounds(0, nvars-1, 0, nvars-1);
        a.setbounds(0, nvars-1, 0, nvars-1);
        matrixmatrixmultiply(sw, 0, nvars-1, 0, nvars-1, false, z, 0, nvars-1, 0, nvars-1, false, 1.0, tm, 0, nvars-1, 0, nvars-1, 0.0, work);
        matrixmatrixmultiply(z, 0, nvars-1, 0, nvars-1, true, tm, 0, nvars-1, 0, nvars-1, false, 1.0, a, 0, nvars-1, 0, nvars-1, 0.0, work);
        for(i = 0; i <= nvars-1; i++)
        {
            for(j = 0; j <= nvars-1; j++)
            {
                a(i,j) = a(i,j)/sqrt(d(i)*d(j));
            }
        }
        if( !smatrixevd(a, nvars, 1, true, d2, z2) )
        {
            info = -4;
            return;
        }
        for(k = 0; k <= nvars-1; k++)
        {
            for(i = 0; i <= nvars-1; i++)
            {
                tf(i) = z2(i,k)/sqrt(d(i));
            }
            for(i = 0; i <= nvars-1; i++)
            {
                v = ap::vdotproduct(&z(i, 0), 1, &tf(0), 1, ap::vlen(0,nvars-1));
                w(i,k) = v;
            }
        }
    }
    
    //
    // Post-processing:
    // * normalization
    // * converting to non-negative form, if possible
    //
    for(k = 0; k <= nvars-1; k++)
    {
        v = ap::vdotproduct(&w(0, k), w.getstride(), &w(0, k), w.getstride(), ap::vlen(0,nvars-1));
        v = 1/sqrt(v);
        ap::vmul(&w(0, k), w.getstride(), ap::vlen(0,nvars-1), v);
        v = 0;
        for(i = 0; i <= nvars-1; i++)
        {
            v = v+w(i,k);
        }
        if( ap::fp_less(v,0) )
        {
            ap::vmul(&w(0, k), w.getstride(), ap::vlen(0,nvars-1), -1);
        }
    }
}

void internalldltrcond(const ap::real_2d_array& l,
     const ap::integer_1d_array& pivots,
     int n,
     bool isupper,
     bool isnormprovided,
     double anorm,
     double& rcond)
{
    int i;
    int kase;
    int k;
    int km1;
    int km2;
    int kp1;
    int kp2;
    double ainvnm;
    ap::real_1d_array work0;
    ap::real_1d_array work1;
    ap::real_1d_array work2;
    ap::integer_1d_array iwork;
    double v;

    ap::ap_error::make_assertion(n>=0, "");
    
    //
    // Check that the diagonal matrix D is nonsingular.
    //
    rcond = 0;
    if( isupper )
    {
        for(i = n; i >= 1; i--)
        {
            if( pivots(i)>0&&ap::fp_eq(l(i,i),0) )
            {
                return;
            }
        }
    }
    else
    {
        for(i = 1; i <= n; i++)
        {
            if( pivots(i)>0&&ap::fp_eq(l(i,i),0) )
            {
                return;
            }
        }
    }
    
    //
    // Estimate the norm of A.
    //
    if( !isnormprovided )
    {
        kase = 0;
        anorm = 0;
        while(true)
        {
            iterativeestimate1norm(n, work1, work0, iwork, anorm, kase);
            if( kase==0 )
            {
                break;
            }
            if( isupper )
            {
                
                //
                // Multiply by U'
                //
                k = n;
                while(k>=1)
                {
                    if( pivots(k)>0 )
                    {
                        
                        //
                        // P(k)
                        //
                        v = work0(k);
                        work0(k) = work0(pivots(k));
                        work0(pivots(k)) = v;
                        
                        //
                        // U(k)
                        //
                        km1 = k-1;
                        v = ap::vdotproduct(&work0(1), 1, &l(1, k), l.getstride(), ap::vlen(1,km1));
                        work0(k) = work0(k)+v;
                        
                        //
                        // Next k
                        //
                        k = k-1;
                    }
                    else
                    {
                        
                        //
                        // P(k)
                        //
//.........这里部分代码省略.........

/*************************************************************************
Application of a sequence of  elementary rotations to a matrix

The algorithm post-multiplies the matrix by a sequence of rotation
transformations which is given by arrays C and S. Depending on the value
of the IsForward parameter either 1 and 2, 3 and 4 and so on (if IsForward=true)
rows are rotated, or the rows N and N-1, N-2 and N-3 and so on are rotated.

Not the whole matrix but only a part of it is transformed (rows from M1
to M2, columns from N1 to N2). Only the elements of this submatrix are changed.

Input parameters:
    IsForward   -   the sequence of the rotation application.
    M1,M2       -   the range of rows to be transformed.
    N1, N2      -   the range of columns to be transformed.
    C,S         -   transformation coefficients.
                    Array whose index ranges within [1..N2-N1].
    A           -   processed matrix.
    WORK        -   working array whose index ranges within [M1..M2].

Output parameters:
    A           -   transformed matrix.

Utility subroutine.
*************************************************************************/
void applyrotationsfromtheright(bool isforward,
     int m1,
     int m2,
     int n1,
     int n2,
     const ap::real_1d_array& c,
     const ap::real_1d_array& s,
     ap::real_2d_array& a,
     ap::real_1d_array& work)
{
    int j;
    int jp1;
    double ctemp;
    double stemp;
    double temp;

    
    //
    // Form A * P'
    //
    if( isforward )
    {
        if( m1!=m2 )
        {
            
            //
            // Common case: M1<>M2
            //
            for(j = n1; j <= n2-1; j++)
            {
                ctemp = c(j-n1+1);
                stemp = s(j-n1+1);
                if( ap::fp_neq(ctemp,1)||ap::fp_neq(stemp,0) )
                {
                    jp1 = j+1;
                    ap::vmove(&work(m1), 1, &a(m1, jp1), a.getstride(), ap::vlen(m1,m2), ctemp);
                    ap::vsub(&work(m1), 1, &a(m1, j), a.getstride(), ap::vlen(m1,m2), stemp);
                    ap::vmul(&a(m1, j), a.getstride(), ap::vlen(m1,m2), ctemp);
                    ap::vadd(&a(m1, j), a.getstride(), &a(m1, jp1), a.getstride(), ap::vlen(m1,m2), stemp);
                    ap::vmove(&a(m1, jp1), a.getstride(), &work(m1), 1, ap::vlen(m1,m2));
                }
            }
        }
        else
        {
            
            //
            // Special case: M1=M2
            //
            for(j = n1; j <= n2-1; j++)
            {
                ctemp = c(j-n1+1);
                stemp = s(j-n1+1);
                if( ap::fp_neq(ctemp,1)||ap::fp_neq(stemp,0) )
                {
                    temp = a(m1,j+1);
                    a(m1,j+1) = ctemp*temp-stemp*a(m1,j);
                    a(m1,j) = stemp*temp+ctemp*a(m1,j);
                }
            }
        }
    }
    else
    {
        if( m1!=m2 )
        {
            
            //
            // Common case: M1<>M2
            //
            for(j = n2-1; j >= n1; j--)
            {
                ctemp = c(j-n1+1);
                stemp = s(j-n1+1);
                if( ap::fp_neq(ctemp,1)||ap::fp_neq(stemp,0) )
//.........这里部分代码省略.........