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

//.........这里部分代码省略.........
        ll = ll+1;
        
        //
        // E(LL) through E(M-1) are nonzero, E(LL-1) is zero
        //
        if( ll==m-1 )
        {
            
            //
            // 2 by 2 block, handle separately
            //
            svdv2x2(d(m-1), e(m-1), d(m), sigmn, sigmx, sinr, cosr, sinl, cosl);
            d(m-1) = sigmx;
            e(m-1) = 0;
            d(m) = sigmn;
            
            //
            // Compute singular vectors, if desired
            //
            if( ncvt>0 )
            {
                mm0 = m+(vstart-1);
                mm1 = m-1+(vstart-1);
                ap::vmove(&vttemp(vstart), &vt(mm1, vstart), ap::vlen(vstart,vend), cosr);
                ap::vadd(&vttemp(vstart), &vt(mm0, vstart), ap::vlen(vstart,vend), sinr);
                ap::vmul(&vt(mm0, vstart), ap::vlen(vstart,vend), cosr);
                ap::vsub(&vt(mm0, vstart), &vt(mm1, vstart), ap::vlen(vstart,vend), sinr);
                ap::vmove(&vt(mm1, vstart), &vttemp(vstart), ap::vlen(vstart,vend));
            }
            if( nru>0 )
            {
                mm0 = m+ustart-1;
                mm1 = m-1+ustart-1;
                ap::vmove(utemp.getvector(ustart, uend), u.getcolumn(mm1, ustart, uend), cosl);
                ap::vadd(utemp.getvector(ustart, uend), u.getcolumn(mm0, ustart, uend), sinl);
                ap::vmul(u.getcolumn(mm0, ustart, uend), cosl);
                ap::vsub(u.getcolumn(mm0, ustart, uend), u.getcolumn(mm1, ustart, uend), sinl);
                ap::vmove(u.getcolumn(mm1, ustart, uend), utemp.getvector(ustart, uend));
            }
            if( ncc>0 )
            {
                mm0 = m+cstart-1;
                mm1 = m-1+cstart-1;
                ap::vmove(&ctemp(cstart), &c(mm1, cstart), ap::vlen(cstart,cend), cosl);
                ap::vadd(&ctemp(cstart), &c(mm0, cstart), ap::vlen(cstart,cend), sinl);
                ap::vmul(&c(mm0, cstart), ap::vlen(cstart,cend), cosl);
                ap::vsub(&c(mm0, cstart), &c(mm1, cstart), ap::vlen(cstart,cend), sinl);
                ap::vmove(&c(mm1, cstart), &ctemp(cstart), ap::vlen(cstart,cend));
            }
            m = m-2;
            continue;
        }
        
        //
        // If working on new submatrix, choose shift direction
        // (from larger end diagonal element towards smaller)
        //
        // Previously was
        //     "if (LL>OLDM) or (M<OLDLL) then"
        // fixed thanks to Michael Rolle < [email protected] >
        // Very strange that LAPACK still contains it.
        //
        bchangedir = false;
        if( idir==1&&fabs(d(ll))<1.0E-3*fabs(d(m)) )
        {
            bchangedir = true;

/*************************************************************************
Reduction of a rectangular matrix to  bidiagonal form

The algorithm reduces the rectangular matrix A to  bidiagonal form by
orthogonal transformations P and Q: A = Q*B*P.

Input parameters:
    A       -   source matrix. array[0..M-1, 0..N-1]
    M       -   number of rows in matrix A.
    N       -   number of columns in matrix A.

Output parameters:
    A       -   matrices Q, B, P in compact form (see below).
    TauQ    -   scalar factors which are used to form matrix Q.
    TauP    -   scalar factors which are used to form matrix P.

The main diagonal and one of the  secondary  diagonals  of  matrix  A  are
replaced with bidiagonal  matrix  B.  Other  elements  contain  elementary
reflections which form MxM matrix Q and NxN matrix P, respectively.

If M>=N, B is the upper  bidiagonal  MxN  matrix  and  is  stored  in  the
corresponding  elements  of  matrix  A.  Matrix  Q  is  represented  as  a
product   of   elementary   reflections   Q = H(0)*H(1)*...*H(n-1),  where
H(i) = 1-tau*v*v'. Here tau is a scalar which is stored  in  TauQ[i],  and
vector v has the following  structure:  v(0:i-1)=0, v(i)=1, v(i+1:m-1)  is
stored   in   elements   A(i+1:m-1,i).   Matrix   P  is  as  follows:  P =
G(0)*G(1)*...*G(n-2), where G(i) = 1 - tau*u*u'. Tau is stored in TauP[i],
u(0:i)=0, u(i+1)=1, u(i+2:n-1) is stored in elements A(i,i+2:n-1).

If M<N, B is the  lower  bidiagonal  MxN  matrix  and  is  stored  in  the
corresponding   elements  of  matrix  A.  Q = H(0)*H(1)*...*H(m-2),  where
H(i) = 1 - tau*v*v', tau is stored in TauQ, v(0:i)=0, v(i+1)=1, v(i+2:m-1)
is    stored    in   elements   A(i+2:m-1,i).    P = G(0)*G(1)*...*G(m-1),
G(i) = 1-tau*u*u', tau is stored in  TauP,  u(0:i-1)=0, u(i)=1, u(i+1:n-1)
is stored in A(i,i+1:n-1).

EXAMPLE:

m=6, n=5 (m > n):               m=5, n=6 (m < n):

(  d   e   u1  u1  u1 )         (  d   u1  u1  u1  u1  u1 )
(  v1  d   e   u2  u2 )         (  e   d   u2  u2  u2  u2 )
(  v1  v2  d   e   u3 )         (  v1  e   d   u3  u3  u3 )
(  v1  v2  v3  d   e  )         (  v1  v2  e   d   u4  u4 )
(  v1  v2  v3  v4  d  )         (  v1  v2  v3  e   d   u5 )
(  v1  v2  v3  v4  v5 )

Here vi and ui are vectors which form H(i) and G(i), and d and e -
are the diagonal and off-diagonal elements of matrix B.
*************************************************************************/
void rmatrixbd(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 maxmn;
    int i;
    double ltau;

    
    //
    // Prepare
    //
    if( n<=0||m<=0 )
    {
        return;
    }
    maxmn = ap::maxint(m, n);
    work.setbounds(0, maxmn);
    t.setbounds(0, maxmn);
    if( m>=n )
    {
        tauq.setbounds(0, n-1);
        taup.setbounds(0, n-1);
    }
    else
    {
        tauq.setbounds(0, m-1);
        taup.setbounds(0, m-1);
    }
    if( m>=n )
    {
        
        //
        // Reduce to upper bidiagonal form
        //
        for(i = 0; i <= n-1; i++)
        {
            
            //
            // Generate elementary reflector H(i) to annihilate A(i+1:m-1,i)
            //
            ap::vmove(t.getvector(1, m-i), a.getcolumn(i, i, m-1));
            generatereflection(t, m-i, ltau);
            tauq(i) = ltau;
            ap::vmove(a.getcolumn(i, i, m-1), t.getvector(1, m-i));
            t(1) = 1;
//.........这里部分代码省略.........

/*************************************************************************
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;
//.........这里部分代码省略.........

/*************************************************************************
Inversion of a matrix given by its LU decomposition.

Input parameters:
    A       -   LU decomposition of the matrix (output of RMatrixLU subroutine).
    Pivots  -   table of permutations which were made during the LU decomposition
                (the output of RMatrixLU subroutine).
    N       -   size of matrix A.

Output parameters:
    A       -   inverse of matrix A.
                Array whose indexes range within [0..N-1, 0..N-1].

Result:
    True, if the matrix is not singular.
    False, if the matrix is singular.

  -- LAPACK routine (version 3.0) --
     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
     Courant Institute, Argonne National Lab, and Rice University
     February 29, 1992
*************************************************************************/
bool rmatrixluinverse(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( !rmatrixtrinverse(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), &work(j+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.getvector(0, n-1), a.getcolumn(j, 0, n-1));
            ap::vmove(a.getcolumn(j, 0, n-1), a.getcolumn(jp, 0, n-1));
            ap::vmove(a.getcolumn(jp, 0, n-1), work.getvector(0, n-1));
        }
    }
    return result;
}

//.........这里部分代码省略.........
    {
        for(l = ai1; l <= ai2; l++)
        {
            for(r = bi1; r <= bi2; r++)
            {
                v = alpha*a(l,aj1+r-bi1);
                k = ci1+l-ai1;
                ap::vadd(&c(k, cj1), &b(r, bj1), ap::vlen(cj1,cj2), v);
            }
        }
        return;
    }
    
    //
    // A*B'
    //
    if( !transa&&transb )
    {
        if( arows*acols<brows*bcols )
        {
            for(r = bi1; r <= bi2; r++)
            {
                for(l = ai1; l <= ai2; l++)
                {
                    v = ap::vdotproduct(&a(l, aj1), &b(r, bj1), ap::vlen(aj1,aj2));
                    c(ci1+l-ai1,cj1+r-bi1) = c(ci1+l-ai1,cj1+r-bi1)+alpha*v;
                }
            }
            return;
        }
        else
        {
            for(l = ai1; l <= ai2; l++)
            {
                for(r = bi1; r <= bi2; r++)
                {
                    v = ap::vdotproduct(&a(l, aj1), &b(r, bj1), ap::vlen(aj1,aj2));
                    c(ci1+l-ai1,cj1+r-bi1) = c(ci1+l-ai1,cj1+r-bi1)+alpha*v;
                }
            }
            return;
        }
    }
    
    //
    // A'*B
    //
    if( transa&&!transb )
    {
        for(l = aj1; l <= aj2; l++)
        {
            for(r = bi1; r <= bi2; r++)
            {
                v = alpha*a(ai1+r-bi1,l);
                k = ci1+l-aj1;
                ap::vadd(&c(k, cj1), &b(r, bj1), ap::vlen(cj1,cj2), v);
            }
        }
        return;
    }
    
    //
    // A'*B'
    //
    if( transa&&transb )
    {
        if( arows*acols<brows*bcols )
        {
            for(r = bi1; r <= bi2; r++)
            {
                for(i = 1; i <= crows; i++)
                {
                    work(i) = 0.0;
                }
                for(l = ai1; l <= ai2; l++)
                {
                    v = alpha*b(r,bj1+l-ai1);
                    k = cj1+r-bi1;
                    ap::vadd(&work(1), &a(l, aj1), ap::vlen(1,crows), v);
                }
                ap::vadd(c.getcolumn(k, ci1, ci2), work.getvector(1, crows));
            }
            return;
        }
        else
        {
            for(l = aj1; l <= aj2; l++)
            {
                k = ai2-ai1+1;
                ap::vmove(work.getvector(1, k), a.getcolumn(l, ai1, ai2));
                for(r = bi1; r <= bi2; r++)
                {
                    v = ap::vdotproduct(&work(1), &b(r, bj1), ap::vlen(1,k));
                    c(ci1+l-aj1,cj1+r-bi1) = c(ci1+l-aj1,cj1+r-bi1)+alpha*v;
                }
            }
            return;
        }
    }
}

//.........这里部分代码省略.........
        ap::vmove(&da(i, 0), &a(i, 0), ap::vlen(0,n-1));
    }
    rmatrixlu(da, n, n, p);
    rep.r1 = rmatrixlurcond1(da, n);
    rep.rinf = rmatrixlurcondinf(da, n);
    if( ap::fp_less(rep.r1,10*ap::machineepsilon)||ap::fp_less(rep.rinf,10*ap::machineepsilon) )
    {
        for(i = 0; i <= n-1; i++)
        {
            for(j = 0; j <= m-1; j++)
            {
                x(i,j) = 0;
            }
        }
        rep.r1 = 0;
        rep.rinf = 0;
        info = -3;
        return;
    }
    info = 1;
    
    //
    // solve
    //
    for(k = 0; k <= m-1; k++)
    {
        
        //
        // First, non-iterative part of solution process:
        // * pivots
        // * L*y = b
        // * U*x = y
        //
        ap::vmove(bc.getvector(0, n-1), b.getcolumn(k, 0, n-1));
        for(i = 0; i <= n-1; i++)
        {
            if( p(i)!=i )
            {
                v = bc(i);
                bc(i) = bc(p(i));
                bc(p(i)) = v;
            }
        }
        y(0) = bc(0);
        for(i = 1; i <= n-1; i++)
        {
            v = ap::vdotproduct(&da(i, 0), &y(0), ap::vlen(0,i-1));
            y(i) = bc(i)-v;
        }
        xc(n-1) = y(n-1)/da(n-1,n-1);
        for(i = n-2; i >= 0; i--)
        {
            v = ap::vdotproduct(&da(i, i+1), &xc(i+1), ap::vlen(i+1,n-1));
            xc(i) = (y(i)-v)/da(i,i);
        }
        
        //
        // Iterative improvement of xc:
        // * calculate r = bc-A*xc using extra-precise dot product
        // * solve A*y = r
        // * update x:=x+r
        //
        // This cycle is executed until one of two things happens:
        // 1. maximum number of iterations reached
        // 2. last iteration decreased error to the lower limit
        //

/*************************************************************************
Obsolete 1-based subroutine.
*************************************************************************/
bool choleskydecomposition(ap::real_2d_array& a, int n, bool isupper)
{
    bool result;
    int i;
    int j;
    double ajj;
    double v;
    int jm1;
    int jp1;

    
    //
    //     Test the input parameters.
    //
    ap::ap_error::make_assertion(n>=0, "Error in CholeskyDecomposition: incorrect function arguments");
    
    //
    //     Quick return if possible
    //
    result = true;
    if( n==0 )
    {
        return result;
    }
    if( isupper )
    {
        
        //
        // Compute the Cholesky factorization A = U'*U.
        //
        for(j = 1; j <= n; j++)
        {
            
            //
            // Compute U(J,J) and test for non-positive-definiteness.
            //
            jm1 = j-1;
            v = ap::vdotproduct(a.getcolumn(j, 1, jm1), a.getcolumn(j, 1, jm1));
            ajj = a(j,j)-v;
            if( ajj<=0 )
            {
                result = false;
                return result;
            }
            ajj = sqrt(ajj);
            a(j,j) = ajj;
            
            //
            // Compute elements J+1:N of row J.
            //
            if( j<n )
            {
                for(i = j+1; i <= n; i++)
                {
                    jm1 = j-1;
                    v = ap::vdotproduct(a.getcolumn(i, 1, jm1), a.getcolumn(j, 1, jm1));
                    a(j,i) = a(j,i)-v;
                }
                v = 1/ajj;
                jp1 = j+1;
                ap::vmul(&a(j, jp1), ap::vlen(jp1,n), v);
            }
        }
    }
    else
    {
        
        //
        // Compute the Cholesky factorization A = L*L'.
        //
        for(j = 1; j <= n; j++)
        {
            
            //
            // Compute L(J,J) and test for non-positive-definiteness.
            //
            jm1 = j-1;
            v = ap::vdotproduct(&a(j, 1), &a(j, 1), ap::vlen(1,jm1));
            ajj = a(j,j)-v;
            if( ajj<=0 )
            {
                result = false;
                return result;
            }
            ajj = sqrt(ajj);
            a(j,j) = ajj;
            
            //
            // Compute elements J+1:N of column J.
            //
            if( j<n )
            {
                for(i = j+1; i <= n; i++)
                {
                    jm1 = j-1;
                    v = ap::vdotproduct(&a(i, 1), &a(j, 1), ap::vlen(1,jm1));
                    a(i,j) = a(i,j)-v;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    ap::real_1d_array t2;
    ap::real_1d_array t3;

    if( n<=0 )
    {
        return;
    }
    t.setbounds(1, n);
    t2.setbounds(1, n);
    t3.setbounds(1, n);
    if( n>1 )
    {
        tau.setbounds(0, n-2);
    }
    d.setbounds(0, n-1);
    if( n>1 )
    {
        e.setbounds(0, n-2);
    }
    if( isupper )
    {
        
        //
        // Reduce the upper triangle of A
        //
        for(i = n-2; i >= 0; i--)
        {
            
            //
            // Generate elementary reflector H() = E - tau * v * v'
            //
            if( i>=1 )
            {
                ap::vmove(t.getvector(2, i+1), a.getcolumn(i+1, 0, i-1));
            }
            t(1) = a(i,i+1);
            generatereflection(t, i+1, taui);
            if( i>=1 )
            {
                ap::vmove(a.getcolumn(i+1, 0, i-1), t.getvector(2, i+1));
            }
            a(i,i+1) = t(1);
            e(i) = a(i,i+1);
            if( taui!=0 )
            {
                
                //
                // Apply H from both sides to A
                //
                a(i,i+1) = 1;
                
                //
                // Compute  x := tau * A * v  storing x in TAU
                //
                ap::vmove(t.getvector(1, i+1), a.getcolumn(i+1, 0, i));
                symmetricmatrixvectormultiply(a, isupper, 0, i, t, taui, t3);
                ap::vmove(&tau(0), &t3(1), ap::vlen(0,i));
                
                //
                // Compute  w := x - 1/2 * tau * (x'*v) * v
                //
                v = ap::vdotproduct(tau.getvector(0, i), a.getcolumn(i+1, 0, i));
                alpha = -0.5*taui*v;
                ap::vadd(tau.getvector(0, i), a.getcolumn(i+1, 0, i), alpha);
                
                //

/*************************************************************************
Unpacking matrix Q which reduces symmetric matrix to a tridiagonal
form.

Input parameters:
    A       -   the result of a SMatrixTD subroutine
    N       -   size of matrix A.
    IsUpper -   storage format (a parameter of SMatrixTD subroutine)
    Tau     -   the result of a SMatrixTD subroutine

Output parameters:
    Q       -   transformation matrix.
                array with elements [0..N-1, 0..N-1].

  -- ALGLIB --
     Copyright 2005-2008 by Bochkanov Sergey
*************************************************************************/
void smatrixtdunpackq(const ap::real_2d_array& a,
     const int& n,
     const bool& isupper,
     const ap::real_1d_array& tau,
     ap::real_2d_array& q)
{
    int i;
    int j;
    ap::real_1d_array v;
    ap::real_1d_array work;

    if( n==0 )
    {
        return;
    }
    
    //
    // init
    //
    q.setbounds(0, n-1, 0, n-1);
    v.setbounds(1, n);
    work.setbounds(0, n-1);
    for(i = 0; i <= n-1; i++)
    {
        for(j = 0; j <= n-1; j++)
        {
            if( i==j )
            {
                q(i,j) = 1;
            }
            else
            {
                q(i,j) = 0;
            }
        }
    }
    
    //
    // unpack Q
    //
    if( isupper )
    {
        for(i = 0; i <= n-2; i++)
        {
            
            //
            // Apply H(i)
            //
            ap::vmove(v.getvector(1, i+1), a.getcolumn(i+1, 0, i));
            v(i+1) = 1;
            applyreflectionfromtheleft(q, tau(i), v, 0, i, 0, n-1, work);
        }
    }
    else
    {
        for(i = n-2; i >= 0; i--)
        {
            
            //
            // Apply H(i)
            //
            ap::vmove(v.getvector(1, n-i-1), a.getcolumn(i, i+1, n-1));
            v(1) = 1;
            applyreflectionfromtheleft(q, tau(i), v, i+1, n-1, 0, n-1, work);
        }
    }
}

/*************************************************************************
Solving a system of linear equations with a system  matrix  given  by  its
Cholesky decomposition.

The algorithm solves systems with a square matrix only.

Input parameters:
    A       -   Cholesky decomposition of a system matrix (the result of
                the SMatrixCholesky subroutine).
    B       -   right side of a system.
                Array whose index ranges within [0..N-1].
    N       -   size of matrix A.
    IsUpper -   points to the triangle of matrix A in which the Cholesky
                decomposition is stored. If IsUpper=True,  the  Cholesky
                decomposition has the form of U'*U, and the upper triangle
                of matrix A stores matrix U (in  that  case,  the  lower
                triangle isn’t used and isn’t changed by the subroutine)
                Similarly, if IsUpper = False, the Cholesky decomposition
                has the form of L*L',  and  the  lower  triangle  stores
                matrix L.

Output parameters:
    X       -   solution of a system.
                Array whose index ranges within [0..N-1].

Result:
    True, if the system is not singular. X contains the solution.
    False, if the system is singular (there is a zero element on the main
diagonal). In this case, X doesn't contain a solution.

  -- ALGLIB --
     Copyright 2005-2008 by Bochkanov Sergey
*************************************************************************/
bool spdmatrixcholeskysolve(const ap::real_2d_array& a,
     ap::real_1d_array b,
     int n,
     bool isupper,
     ap::real_1d_array& x)
{
    bool result;
    int i;
    double v;

    ap::ap_error::make_assertion(n>0, "Error: N<=0 in SolveSystemCholesky");
    
    //
    // det(A)=0?
    //
    result = true;
    for(i = 0; i <= n-1; i++)
    {
        if( ap::fp_eq(a(i,i),0) )
        {
            result = false;
            return result;
        }
    }
    
    //
    // det(A)<>0
    //
    x.setbounds(0, n-1);
    if( isupper )
    {
        
        //
        // A = U'*U, solve U'*y = b first
        //
        b(0) = b(0)/a(0,0);
        for(i = 1; i <= n-1; i++)
        {
            v = ap::vdotproduct(a.getcolumn(i, 0, i-1), b.getvector(0, i-1));
            b(i) = (b(i)-v)/a(i,i);
        }
        
        //
        // Solve U*x = y
        //
        b(n-1) = b(n-1)/a(n-1,n-1);
        for(i = n-2; i >= 0; i--)
        {
            v = ap::vdotproduct(&a(i, i+1), &b(i+1), ap::vlen(i+1,n-1));
            b(i) = (b(i)-v)/a(i,i);
        }
        ap::vmove(&x(0), &b(0), ap::vlen(0,n-1));
    }
    else
    {
        
        //
        // A = L*L', solve L'*y = b first
        //
        b(0) = b(0)/a(0,0);
        for(i = 1; i <= n-1; i++)
        {
            v = ap::vdotproduct(&a(i, 0), &b(0), ap::vlen(0,i-1));
            b(i) = (b(i)-v)/a(i,i);
        }
        
        //
//.........这里部分代码省略.........

/*************************************************************************
Обращение треугольной матрицы

Подпрограмма обращает следующие типы матриц:
    * верхнетреугольные
    * верхнетреугольные с единичной диагональю
    * нижнетреугольные
    * нижнетреугольные с единичной диагональю

В случае, если матрица верхне(нижне)треугольная, то  матрица,  обратная  к
ней, тоже верхне(нижне)треугольная, и после  завершения  работы  алгоритма
обратная матрица замещает переданную. При этом элементы расположенные ниже
(выше) диагонали не меняются в ходе работы алгоритма.

Если матрица с единичной  диагональю, то обратная к  ней  матрица  тоже  с
единичной  диагональю.  В  алгоритм  передаются  только    внедиагональные
элементы. При этом в результате работы алгоритма диагональные элементы  не
меняются.

Входные параметры:
    A           -   матрица. Массив с нумерацией элементов [0..N-1,0..N-1]
    N           -   размер матрицы
    IsUpper     -   True, если матрица верхнетреугольная
    IsunitTriangular-   True, если матрица с единичной диагональю.

Выходные параметры:
    A           -   матрица, обратная к входной, если задача не вырождена.

Результат:
    True, если матрица не вырождена
    False, если матрица вырождена

  -- LAPACK routine (version 3.0) --
     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
     Courant Institute, Argonne National Lab, and Rice University
     February 29, 1992
*************************************************************************/
bool rmatrixtrinverse(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( 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.getvector(0, j-1), a.getcolumn(j, 0, j-1));
                for(i = 0; i <= j-1; i++)
                {
                    if( i<j-1 )
                    {
                        v = ap::vdotproduct(&a(i, i+1), &t(i+1), ap::vlen(i+1,j-1));
                    }
                    else
                    {
                        v = 0;
                    }
                    if( nounit )
                    {
                        a(i,j) = v+a(i,i)*t(i);
                    }
//.........这里部分代码省略.........

bool solvesystemcholesky(const ap::real_2d_array& a,
     ap::real_1d_array b,
     int n,
     bool isupper,
     ap::real_1d_array& x)
{
    bool result;
    int i;
    int im1;
    int ip1;
    double v;

    ap::ap_error::make_assertion(n>0, "Error: N<=0 in SolveSystemCholesky");
    
    //
    // det(A)=0?
    //
    result = true;
    for(i = 1; i <= n; i++)
    {
        if( ap::fp_eq(a(i,i),0) )
        {
            result = false;
            return result;
        }
    }
    
    //
    // det(A)<>0
    //
    x.setbounds(1, n);
    if( isupper )
    {
        
        //
        // A = U'*U, solve U'*y = b first
        //
        b(1) = b(1)/a(1,1);
        for(i = 2; i <= n; i++)
        {
            im1 = i-1;
            v = ap::vdotproduct(a.getcolumn(i, 1, im1), b.getvector(1, im1));
            b(i) = (b(i)-v)/a(i,i);
        }
        
        //
        // Solve U*x = y
        //
        b(n) = b(n)/a(n,n);
        for(i = n-1; i >= 1; i--)
        {
            ip1 = i+1;
            v = ap::vdotproduct(&a(i, ip1), &b(ip1), ap::vlen(ip1,n));
            b(i) = (b(i)-v)/a(i,i);
        }
        ap::vmove(&x(1), &b(1), ap::vlen(1,n));
    }
    else
    {
        
        //
        // A = L*L', solve L'*y = b first
        //
        b(1) = b(1)/a(1,1);
        for(i = 2; i <= n; i++)
        {
            im1 = i-1;
            v = ap::vdotproduct(&a(i, 1), &b(1), ap::vlen(1,im1));
            b(i) = (b(i)-v)/a(i,i);
        }
        
        //
        // Solve L'*x = y
        //
        b(n) = b(n)/a(n,n);
        for(i = n-1; i >= 1; i--)
        {
            ip1 = i+1;
            v = ap::vdotproduct(a.getcolumn(i, ip1, n), b.getvector(ip1, n));
            b(i) = (b(i)-v)/a(i,i);
        }
        ap::vmove(&x(1), &b(1), ap::vlen(1,n));
    }
    return result;
}

//.........这里部分代码省略.........
        ap::vmove(&wr1s(0), &wr1(0), ap::vlen(0,n-1));
        ap::vmove(&wi1s(0), &wi1(0), ap::vlen(0,n-1));
        for(i = 0; i <= n-1; i++)
        {
            for(j = 0; j <= n-2-i; j++)
            {
                if( wr1s(j)>wr1s(j+1) )
                {
                    tmp = wr1s(j);
                    wr1s(j) = wr1s(j+1);
                    wr1s(j+1) = tmp;
                    tmp = wi1s(j);
                    wi1s(j) = wi1s(j+1);
                    wi1s(j+1) = tmp;
                }
            }
        }
        for(i = 0; i <= n-1; i++)
        {
            valonlydiff = ap::maxreal(valonlydiff, fabs(wr0s(i)-wr1s(i)));
            valonlydiff = ap::maxreal(valonlydiff, fabs(wi0s(i)-wi1s(i)));
        }
        
        //
        // Test right vectors
        //
        if( needr )
        {
            k = 0;
            while(k<=n-1)
            {
                if( wi1(k)==0 )
                {
                    ap::vmove(vec1r.getvector(0, n-1), vr.getcolumn(k, 0, n-1));
                    for(i = 0; i <= n-1; i++)
                    {
                        vec1i(i) = 0;
                    }
                    curwr = wr1(k);
                    curwi = 0;
                }
                if( wi1(k)>0 )
                {
                    ap::vmove(vec1r.getvector(0, n-1), vr.getcolumn(k, 0, n-1));
                    ap::vmove(vec1i.getvector(0, n-1), vr.getcolumn(k+1, 0, n-1));
                    curwr = wr1(k);
                    curwi = wi1(k);
                }
                if( wi1(k)<0 )
                {
                    ap::vmove(vec1r.getvector(0, n-1), vr.getcolumn(k-1, 0, n-1));
                    ap::vmoveneg(vec1i.getvector(0, n-1), vr.getcolumn(k, 0, n-1));
                    curwr = wr1(k);
                    curwi = wi1(k);
                }
                for(i = 0; i <= n-1; i++)
                {
                    vt = ap::vdotproduct(&a(i, 0), &vec1r(0), ap::vlen(0,n-1));
                    vec2r(i) = vt;
                    vt = ap::vdotproduct(&a(i, 0), &vec1i(0), ap::vlen(0,n-1));
                    vec2i(i) = vt;
                }
                ap::vmove(&vec3r(0), &vec1r(0), ap::vlen(0,n-1), curwr);
                ap::vsub(&vec3r(0), &vec1i(0), ap::vlen(0,n-1), curwi);
                ap::vmove(&vec3i(0), &vec1r(0), ap::vlen(0,n-1), curwi);
                ap::vadd(&vec3i(0), &vec1i(0), ap::vlen(0,n-1), curwr);

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

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

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.
    TAUQ        -   scalar factors which are used to form Q.
                    Output of ToBidiagonal subroutine.
    Z           -   multiplied matrix.
                    array[0..ZRows-1,0..ZColumns-1]
    ZRows       -   number of rows in matrix Z. If FromTheRight=False,
                    ZRows=M, otherwise ZRows can be arbitrary.
    ZColumns    -   number of columns in matrix Z. If FromTheRight=True,
                    ZColumns=M, otherwise ZColumns can be arbitrary.
    FromTheRight -  pre- or post-multiply.
    DoTranspose -   multiply by Q or Q'.

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

  -- ALGLIB --
     Copyright 2005 by Bochkanov Sergey
*************************************************************************/
void rmatrixbdmultiplybyq(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 i1;
    int i2;
    int istep;
    ap::real_1d_array v;
    ap::real_1d_array work;
    int mx;

    if( m<=0||n<=0||zrows<=0||zcolumns<=0 )
    {
        return;
    }
    ap::ap_error::make_assertion(fromtheright&&zcolumns==m||!fromtheright&&zrows==m, "RMatrixBDMultiplyByQ: 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);
    if( m>=n )
    {
        
        //
        // setup
        //
        if( fromtheright )
        {
            i1 = 0;
            i2 = n-1;
            istep = +1;
        }
        else
        {
            i1 = n-1;
            i2 = 0;
            istep = -1;
        }
        if( dotranspose )
        {
            i = i1;
            i1 = i2;
            i2 = i;
            istep = -istep;
        }
        
        //
        // Process
        //
        i = i1;
        do
        {
            ap::vmove(v.getvector(1, m-i), qp.getcolumn(i, i, m-1));
            v(1) = 1;
            if( fromtheright )
            {
                applyreflectionfromtheright(z, tauq(i), v, 0, zrows-1, i, m-1, work);
            }
            else
//.........这里部分代码省略.........

/*************************************************************************
Linear regression

Variant of LRBuild which uses vector of standatd deviations (errors in
function values).

INPUT PARAMETERS:
    XY          -   training set, array [0..NPoints-1,0..NVars]:
                    * NVars columns - independent variables
                    * last column - dependent variable
    S           -   standard deviations (errors in function values)
                    array[0..NPoints-1], S[i]>0.
    NPoints     -   training set size, NPoints>NVars+1
    NVars       -   number of independent variables

OUTPUT PARAMETERS:
    Info        -   return code:
                    * -255, in case of unknown internal error
                    * -4, if internal SVD subroutine haven't converged
                    * -1, if incorrect parameters was passed (NPoints<NVars+2, NVars<1).
                    * -2, if S[I]<=0
                    *  1, if subroutine successfully finished
    LM          -   linear model in the ALGLIB format. Use subroutines of
                    this unit to work with the model.
    AR          -   additional results


  -- ALGLIB --
     Copyright 02.08.2008 by Bochkanov Sergey
*************************************************************************/
void lrbuilds(const ap::real_2d_array& xy,
     const ap::real_1d_array& s,
     int npoints,
     int nvars,
     int& info,
     linearmodel& lm,
     lrreport& ar)
{
    ap::real_2d_array xyi;
    ap::real_1d_array x;
    ap::real_1d_array means;
    ap::real_1d_array sigmas;
    int i;
    int j;
    double v;
    int offs;
    double mean;
    double variance;
    double skewness;
    double kurtosis;

    
    //
    // Test parameters
    //
    if( npoints<=nvars+1||nvars<1 )
    {
        info = -1;
        return;
    }
    
    //
    // Copy data, add one more column (constant term)
    //
    xyi.setbounds(0, npoints-1, 0, nvars+1);
    for(i = 0; i <= npoints-1; i++)
    {
        ap::vmove(&xyi(i, 0), &xy(i, 0), ap::vlen(0,nvars-1));
        xyi(i,nvars) = 1;
        xyi(i,nvars+1) = xy(i,nvars);
    }
    
    //
    // Standartization
    //
    x.setbounds(0, npoints-1);
    means.setbounds(0, nvars-1);
    sigmas.setbounds(0, nvars-1);
    for(j = 0; j <= nvars-1; j++)
    {
        ap::vmove(x.getvector(0, npoints-1), xy.getcolumn(j, 0, npoints-1));
        calculatemoments(x, npoints, mean, variance, skewness, kurtosis);
        means(j) = mean;
        sigmas(j) = sqrt(variance);
        if( ap::fp_eq(sigmas(j),0) )
        {
            sigmas(j) = 1;
        }
        for(i = 0; i <= npoints-1; i++)
        {
            xyi(i,j) = (xyi(i,j)-means(j))/sigmas(j);
        }
    }
    
    //
    // Internal processing
    //
    lrinternal(xyi, s, npoints, nvars+1, info, lm, ar);
    if( info<0 )
    {
//.........这里部分代码省略.........