C++ real_1d_array::getvector方法代码示例

void lbfgslincomb(const int& n,
     const double& da,
     const ap::real_1d_array& dx,
     int sx,
     ap::real_1d_array& dy,
     int sy)
{
    int fx;
    int fy;

    fx = sx+n-1;
    fy = sy+n-1;
    ap::vadd(dy.getvector(sy, fy), dx.getvector(sx, fx), da);
}

double lbfgsdotproduct(const int& n,
     const ap::real_1d_array& dx,
     int sx,
     const ap::real_1d_array& dy,
     int sy)
{
    double result;
    double v;
    int fx;
    int fy;

    fx = sx+n-1;
    fy = sy+n-1;
    v = ap::vdotproduct(dx.getvector(sx, fx), dy.getvector(sy, fy));
    result = v;
    return result;
}

/*************************************************************************
Dense solver.

Similar to RMatrixSolveM() but solves task with one right part  (where b/x
are vectors, not matrices).

See RMatrixSolveM()  description  for  more  information  about subroutine
parameters.

  -- ALGLIB --
     Copyright 24.08.2009 by Bochkanov Sergey
*************************************************************************/
void rmatrixsolve(const ap::real_2d_array& a,
     int n,
     const ap::real_1d_array& b,
     int& info,
     densesolverreport& rep,
     ap::real_1d_array& x)
{
    ap::real_2d_array bm;
    ap::real_2d_array xm;

    if( n<=0 )
    {
        info = -1;
        return;
    }
    bm.setlength(n, 1);
    ap::vmove(bm.getcolumn(0, 0, n-1), b.getvector(0, n-1));
    rmatrixsolvem(a, n, bm, 1, info, rep, xm);
    x.setlength(n);
    ap::vmove(x.getvector(0, n-1), xm.getcolumn(0, 0, n-1));
}

/*************************************************************************
Obsolete 1-based subroutine
*************************************************************************/
void shermanmorrisonupdateuv(ap::real_2d_array& inva,
     int n,
     const ap::real_1d_array& u,
     const ap::real_1d_array& v)
{
    ap::real_1d_array t1;
    ap::real_1d_array t2;
    int i;
    int j;
    double lambda;
    double vt;

    t1.setbounds(1, n);
    t2.setbounds(1, n);
    
    //
    // T1 = InvA * U
    // Lambda = v * T1
    //
    for(i = 1; i <= n; i++)
    {
        vt = ap::vdotproduct(&inva(i, 1), &u(1), ap::vlen(1,n));
        t1(i) = vt;
    }
    lambda = ap::vdotproduct(&v(1), &t1(1), ap::vlen(1,n));
    
    //
    // T2 = v*InvA
    //
    for(j = 1; j <= n; j++)
    {
        vt = ap::vdotproduct(v.getvector(1, n), inva.getcolumn(j, 1, n));
        t2(j) = vt;
    }
    
    //
    // InvA = InvA - correction
    //
    for(i = 1; i <= n; i++)
    {
        vt = t1(i)/(1+lambda);
        ap::vsub(&inva(i, 1), &t2(1), ap::vlen(1,n), vt);
    }
}

/*************************************************************************
Inverse matrix update by the Sherman-Morrison formula

The algorithm computes the inverse of matrix A+u*v’ by using the given matrix
A^-1 and the vectors u and v.

Input parameters:
    InvA    -   inverse of matrix A.
                Array whose indexes range within [0..N-1, 0..N-1].
    N       -   size of matrix A.
    U       -   the vector modifying the matrix.
                Array whose index ranges within [0..N-1].
    V       -   the vector modifying the matrix.
                Array whose index ranges within [0..N-1].

Output parameters:
    InvA - inverse of matrix A + u*v'.

  -- ALGLIB --
     Copyright 2005 by Bochkanov Sergey
*************************************************************************/
void rmatrixinvupdateuv(ap::real_2d_array& inva,
     int n,
     const ap::real_1d_array& u,
     const ap::real_1d_array& v)
{
    ap::real_1d_array t1;
    ap::real_1d_array t2;
    int i;
    int j;
    double lambda;
    double vt;

    t1.setbounds(0, n-1);
    t2.setbounds(0, n-1);
    
    //
    // T1 = InvA * U
    // Lambda = v * T1
    //
    for(i = 0; i <= n-1; i++)
    {
        vt = ap::vdotproduct(&inva(i, 0), &u(0), ap::vlen(0,n-1));
        t1(i) = vt;
    }
    lambda = ap::vdotproduct(&v(0), &t1(0), ap::vlen(0,n-1));
    
    //
    // T2 = v*InvA
    //
    for(j = 0; j <= n-1; j++)
    {
        vt = ap::vdotproduct(v.getvector(0, n-1), inva.getcolumn(j, 0, n-1));
        t2(j) = vt;
    }
    
    //
    // InvA = InvA - correction
    //
    for(i = 0; i <= n-1; i++)
    {
        vt = t1(i)/(1+lambda);
        ap::vsub(&inva(i, 0), &t2(0), ap::vlen(0,n-1), vt);
    }
}

/*************************************************************************
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( ctemp!=1||stemp!=0 )
                {
                    jp1 = j+1;
                    ap::vmove(work.getvector(m1, m2), a.getcolumn(jp1, m1, m2), ctemp);
                    ap::vsub(work.getvector(m1, m2), a.getcolumn(j, m1, m2), stemp);
                    ap::vmul(a.getcolumn(j, m1, m2), ctemp);
                    ap::vadd(a.getcolumn(j, m1, m2), a.getcolumn(jp1, m1, m2), stemp);
                    ap::vmove(a.getcolumn(jp1, m1, m2), work.getvector(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( ctemp!=1||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( ctemp!=1||stemp!=0 )
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    g.setbounds(1, n);
    xold.setbounds(1, n);
    diag.setbounds(1, n);
    (*funcgrad)(x, f, g, params);
    iter = 0;
    info = 0;
    if( n<=0||m<=0||m>n||epsg<0||epsf<0||epsx<0||maxits<0 )
    {
        info = -1;
        return;
    }
    nfun = 1;
    point = 0;
    for(i = 1; i <= n; i++)
    {
        diag(i) = 1;
    }
    xtol = 100*ap::machineepsilon;
    gtol = 0.9;
    stpmin = pow(double(10), double(-20));
    stpmax = pow(double(10), double(20));
    ispt = n+2*m;
    iypt = ispt+n*m;
    for(i = 1; i <= n; i++)
    {
        w(ispt+i) = -g(i);
    }
    gnorm = sqrt(lbfgsdotproduct(n, g, 1, g, 1));
    stp1 = 1/gnorm;
    ftol = 0.0001;
    maxfev = 20;
    while(true)
    {
        ap::vmove(xold.getvector(1, n), x.getvector(1, n));
        fold = f;

        iter = iter+1;
        info = 0;
        bound = iter-1;
        if( iter!=1 )
        {
            if( iter>m )
            {
                bound = m;
            }
            ys = lbfgsdotproduct(n, w, iypt+npt+1, w, ispt+npt+1);
            yy = lbfgsdotproduct(n, w, iypt+npt+1, w, iypt+npt+1);
            for(i = 1; i <= n; i++)
            {
                diag(i) = ys/yy;
            }
            cp = point;
            if( point==0 )
            {
                cp = m;
            }
            w(n+cp) = 1/ys;
            for(i = 1; i <= n; i++)
            {
                w(i) = -g(i);
            }
            cp = point;
            for(i = 1; i <= bound; i++)
            {
                cp = cp-1;
                if( cp==-1 )

/*************************************************************************
Obsolete 1-based subroutine
*************************************************************************/
void totridiagonal(ap::real_2d_array& a,
     int n,
     bool isupper,
     ap::real_1d_array& tau,
     ap::real_1d_array& d,
     ap::real_1d_array& e)
{
    int i;
    int ip1;
    int im1;
    int nmi;
    int nm1;
    double alpha;
    double taui;
    double v;
    ap::real_1d_array t;
    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);
    tau.setbounds(1, ap::maxint(1, n-1));
    d.setbounds(1, n);
    e.setbounds(1, ap::maxint(1, n-1));
    if( isupper )
    {
        
        //
        // Reduce the upper triangle of A
        //
        for(i = n-1; i >= 1; i--)
        {
            
            //
            // Generate elementary reflector H(i) = I - tau * v * v'
            // to annihilate A(1:i-1,i+1)
            //
            // DLARFG( I, A( I, I+1 ), A( 1, I+1 ), 1, TAUI );
            //
            ip1 = i+1;
            im1 = i-1;
            if( i>=2 )
            {
                ap::vmove(t.getvector(2, i), a.getcolumn(ip1, 1, im1));
            }
            t(1) = a(i,ip1);
            generatereflection(t, i, taui);
            if( i>=2 )
            {
                ap::vmove(a.getcolumn(ip1, 1, im1), t.getvector(2, i));
            }
            a(i,ip1) = t(1);
            e(i) = a(i,i+1);
            if( taui!=0 )
            {
                
                //
                // Apply H(i) from both sides to A(1:i,1:i)
                //
                a(i,i+1) = 1;
                
                //
                // Compute  x := tau * A * v  storing x in TAU(1:i)
                //
                // DSYMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO, TAU, 1 );
                //
                ip1 = i+1;
                ap::vmove(t.getvector(1, i), a.getcolumn(ip1, 1, i));
                symmetricmatrixvectormultiply(a, isupper, 1, i, t, taui, tau);
                
                //
                // Compute  w := x - 1/2 * tau * (x'*v) * v
                //
                ip1 = i+1;
                v = ap::vdotproduct(tau.getvector(1, i), a.getcolumn(ip1, 1, i));
                alpha = -0.5*taui*v;
                ap::vadd(tau.getvector(1, i), a.getcolumn(ip1, 1, i), alpha);
                
                //
                // Apply the transformation as a rank-2 update:
                //    A := A - v * w' - w * v'
                //
                // DSYR2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A, LDA );
                //
                ap::vmove(t.getvector(1, i), a.getcolumn(ip1, 1, i));
                symmetricrank2update(a, isupper, 1, i, t, tau, t2, double(-1));
                a(i,i+1) = e(i);
            }
            d(i+1) = a(i+1,i+1);
            tau(i) = taui;
        }
        d(1) = a(1,1);
//.........这里部分代码省略.........

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

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

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;
}

//.........这里部分代码省略.........
                    grow = grow/xj;
                    j = j+jinc;
                }
            }
        }
    }
    if( ap::fp_greater(grow*tscal,smlnum) )
    {
        
        //
        // Use the Level 2 BLAS solve if the reciprocal of the bound on
        // elements of X is not too small.
        //
        if( upper&&notran||!upper&&!notran )
        {
            if( nounit )
            {
                vd = a(n,n);
            }
            else
            {
                vd = 1;
            }
            x(n) = x(n)/vd;
            for(i = n-1; i >= 1; i--)
            {
                ip1 = i+1;
                if( upper )
                {
                    v = ap::vdotproduct(&a(i, ip1), &x(ip1), ap::vlen(ip1,n));
                }
                else
                {
                    v = ap::vdotproduct(a.getcolumn(i, ip1, n), x.getvector(ip1, n));
                }
                if( nounit )
                {
                    vd = a(i,i);
                }
                else
                {
                    vd = 1;
                }
                x(i) = (x(i)-v)/vd;
            }
        }
        else
        {
            if( nounit )
            {
                vd = a(1,1);
            }
            else
            {
                vd = 1;
            }
            x(1) = x(1)/vd;
            for(i = 2; i <= n; i++)
            {
                im1 = i-1;
                if( upper )
                {
                    v = ap::vdotproduct(a.getcolumn(i, 1, im1), x.getvector(1, im1));
                }
                else
                {