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

/*************************************************************************
LU-разложение матрицы общего вида размера M x N

Использует  LUDecomposition.   По  функциональности  отличается  тем,  что
выводит  матрицы  L  и  U не в компактной форме, а в виде отдельных матриц
общего вида, заполненных в соответствующих местах нулевыми элементами.

Подпрограмма приведена исключительно для демонстрации того, как
"распаковывается" результат работы подпрограммы LUDecomposition

  -- ALGLIB --
     Copyright 2005 by Bochkanov Sergey
*************************************************************************/
void ludecompositionunpacked(ap::real_2d_array a,
                             int m,
                             int n,
                             ap::real_2d_array& l,
                             ap::real_2d_array& u,
                             ap::integer_1d_array& pivots)
{
    int i;
    int j;
    int minmn;

    if( m==0||n==0 )
    {
        return;
    }
    minmn = ap::minint(m, n);
    l.setbounds(1, m, 1, minmn);
    u.setbounds(1, minmn, 1, n);
    ludecomposition(a, m, n, pivots);
    for(i = 1; i <= m; i++)
    {
        for(j = 1; j <= minmn; j++)
        {
            if( j>i )
            {
                l(i,j) = 0;
            }
            if( j==i )
            {
                l(i,j) = 1;
            }
            if( j<i )
            {
                l(i,j) = a(i,j);
            }
        }
    }
    for(i = 1; i <= minmn; i++)
    {
        for(j = 1; j <= n; j++)
        {
            if( j<i )
            {
                u(i,j) = 0;
            }
            if( j>=i )
            {
                u(i,j) = a(i,j);
            }
        }
    }
}

void qrdecompositionunpacked(ap::real_2d_array a,
     int m,
     int n,
     ap::real_2d_array& q,
     ap::real_2d_array& r)
{
    int i;
    int k;
    ap::real_1d_array tau;
    ap::real_1d_array work;
    ap::real_1d_array v;

    k = ap::minint(m, n);
    if( n<=0 )
    {
        return;
    }
    work.setbounds(1, m);
    v.setbounds(1, m);
    q.setbounds(1, m, 1, m);
    r.setbounds(1, m, 1, n);
    
    //
    // QRDecomposition
    //
    qrdecomposition(a, m, n, tau);
    
    //
    // R
    //
    for(i = 1; i <= n; i++)
    {
        r(1,i) = 0;
    }
    for(i = 2; i <= m; i++)
    {
        ap::vmove(&r(i, 1), &r(1, 1), ap::vlen(1,n));
    }
    for(i = 1; i <= k; i++)
    {
        ap::vmove(&r(i, i), &a(i, i), ap::vlen(i,n));
    }
    
    //
    // Q
    //
    unpackqfromqr(a, m, n, tau, m, q);
}

/*************************************************************************
Unpacking of matrix R from the QR decomposition of a matrix A

Input parameters:
    A       -   matrices Q and R in compact form.
                Output of RMatrixQR subroutine.
    M       -   number of rows in given matrix A. M>=0.
    N       -   number of columns in given matrix A. N>=0.

Output parameters:
    R       -   matrix R, array[0..M-1, 0..N-1].

  -- ALGLIB --
     Copyright 2005 by Bochkanov Sergey
*************************************************************************/
void rmatrixqrunpackr(const ap::real_2d_array& a,
     int m,
     int n,
     ap::real_2d_array& r)
{
    int i;
    int k;

    if( m<=0||n<=0 )
    {
        return;
    }
    k = ap::minint(m, n);
    r.setbounds(0, m-1, 0, n-1);
    for(i = 0; i <= n-1; i++)
    {
        r(0,i) = 0;
    }
    for(i = 1; i <= m-1; i++)
    {
        ap::vmove(&r(i, 0), &r(0, 0), ap::vlen(0,n-1));
    }
    for(i = 0; i <= k-1; i++)
    {
        ap::vmove(&r(i, i), &a(i, i), ap::vlen(i,n-1));
    }
}

bool in_out_variable(const ap::boolean_1d_array& in, const ap::real_2d_array& X, ap::real_2d_array& x, bool io)
{
	//////////////////////////////////////////////////////////////////
	// Section: Define variables
	int rows = in.gethighbound(0) + 1;
	bool flag;
	vector<int> stdVector;
	//////////////////////////////////////////////////////////////////
	// Section: Identify how many variables are in or out

	for (int i=0; i<rows; i++)
	{
		if (in(i)==io) 
			stdVector.push_back(i);
	}
	if (stdVector.size()>0)
	{
		// Routine to extract the in/out variables
		x.setbounds(0,X.gethighbound(1),0,static_cast<int>(stdVector.size())-1);
		for (size_t i=0; i<stdVector.size(); i++)
			ap::vmove(x.getcolumn(static_cast<int>(i),0,X.gethighbound(1)), X.getcolumn(stdVector[i],0,X.gethighbound(1)));
	
        flag=TRUE;
	}
	else
		flag=FALSE;

	return flag;
}

void unpackqfromqr(const ap::real_2d_array& a,
     int m,
     int n,
     const ap::real_1d_array& tau,
     int qcolumns,
     ap::real_2d_array& q)
{
    int i;
    int j;
    int k;
    int minmn;
    ap::real_1d_array v;
    ap::real_1d_array work;
    int vm;

    ap::ap_error::make_assertion(qcolumns<=m, "UnpackQFromQR: QColumns>M!");
    if( m==0||n==0||qcolumns==0 )
    {
        return;
    }
    
    //
    // init
    //
    minmn = ap::minint(m, n);
    k = ap::minint(minmn, qcolumns);
    q.setbounds(1, m, 1, qcolumns);
    v.setbounds(1, m);
    work.setbounds(1, qcolumns);
    for(i = 1; i <= m; i++)
    {
        for(j = 1; j <= qcolumns; j++)
        {
            if( i==j )
            {
                q(i,j) = 1;
            }
            else
            {
                q(i,j) = 0;
            }
        }
    }
    
    //
    // unpack Q
    //
    for(i = k; i >= 1; i--)
    {
        
        //
        // Apply H(i)
        //
        vm = m-i+1;
        ap::vmove(v.getvector(1, vm), a.getcolumn(i, i, m));
        v(1) = 1;
        applyreflectionfromtheleft(q, tau(i), v, i, m, 1, qcolumns, work);
    }
}

/*************************************************************************
Copy
*************************************************************************/
static void rmatrixmakeacopy(const ap::real_2d_array& a,
     int m,
     int n,
     ap::real_2d_array& b)
{
    int i;
    int j;

    b.setbounds(0, m-1, 0, n-1);
    for(i = 0; i <= m-1; i++)
    {
        for(j = 0; j <= n-1; j++)
        {
            b(i,j) = a(i,j);
        }
    }
}

/*************************************************************************
Generation of random NxN symmetric positive definite matrix with given
condition number and norm2(A)=1

INPUT PARAMETERS:
    N   -   matrix size
    C   -   condition number (in 2-norm)

OUTPUT PARAMETERS:
    A   -   random SPD matrix with norm2(A)=1 and cond(A)=C

  -- ALGLIB routine --
     04.12.2009
     Bochkanov Sergey
*************************************************************************/
void spdmatrixrndcond(int n, double c, ap::real_2d_array& a)
{
    int i;
    int j;
    double l1;
    double l2;

    
    //
    // Special cases
    //
    if( n<=0||ap::fp_less(c,1) )
    {
        return;
    }
    a.setbounds(0, n-1, 0, n-1);
    if( n==1 )
    {
        a(0,0) = 1;
        return;
    }
    
    //
    // Prepare matrix
    //
    l1 = 0;
    l2 = log(1/c);
    for(i = 0; i <= n-1; i++)
    {
        for(j = 0; j <= n-1; j++)
        {
            a(i,j) = 0;
        }
    }
    a(0,0) = exp(l1);
    for(i = 1; i <= n-2; i++)
    {
        a(i,i) = exp(ap::randomreal()*(l2-l1)+l1);
    }
    a(n-1,n-1) = exp(l2);
    
    //
    // Multiply
    //
    smatrixrndmultiply(a, n);
}

/*************************************************************************
Generation of random NxN symmetric matrix with given condition number  and
norm2(A)=1

INPUT PARAMETERS:
    N   -   matrix size
    C   -   condition number (in 2-norm)

OUTPUT PARAMETERS:
    A   -   random matrix with norm2(A)=1 and cond(A)=C

  -- ALGLIB routine --
     04.12.2009
     Bochkanov Sergey
*************************************************************************/
void smatrixrndcond(int n, double c, ap::real_2d_array& a)
{
    int i;
    int j;
    double l1;
    double l2;

    ap::ap_error::make_assertion(n>=1&&ap::fp_greater_eq(c,1), "SMatrixRndCond: N<1 or C<1!");
    a.setbounds(0, n-1, 0, n-1);
    if( n==1 )
    {
        
        //
        // special case
        //
        a(0,0) = 2*ap::randominteger(2)-1;
        return;
    }
    
    //
    // Prepare matrix
    //
    l1 = 0;
    l2 = log(1/c);
    for(i = 0; i <= n-1; i++)
    {
        for(j = 0; j <= n-1; j++)
        {
            a(i,j) = 0;
        }
    }
    a(0,0) = exp(l1);
    for(i = 1; i <= n-2; i++)
    {
        a(i,i) = (2*ap::randominteger(2)-1)*exp(ap::randomreal()*(l2-l1)+l1);
    }
    a(n-1,n-1) = exp(l2);
    
    //
    // Multiply
    //
    smatrixrndmultiply(a, n);
}

/*************************************************************************
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);
}

/*************************************************************************
Unpacking matrix Q which reduces a matrix to bidiagonal form.

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.
    QColumns    -   required number of columns in matrix Q.
                    M>=QColumns>=0.

Output parameters:
    Q           -   first QColumns columns of matrix Q.
                    Array[0..M-1, 0..QColumns-1]
                    If QColumns=0, the array is not modified.

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

    ap::ap_error::make_assertion(qcolumns<=m, "RMatrixBDUnpackQ: QColumns>M!");
    ap::ap_error::make_assertion(qcolumns>=0, "RMatrixBDUnpackQ: QColumns<0!");
    if( m==0||n==0||qcolumns==0 )
    {
        return;
    }
    
    //
    // prepare Q
    //
    q.setbounds(0, m-1, 0, qcolumns-1);
    for(i = 0; i <= m-1; i++)
    {
        for(j = 0; j <= qcolumns-1; j++)
        {
            if( i==j )
            {
                q(i,j) = 1;
            }
            else
            {
                q(i,j) = 0;
            }
        }
    }
    
    //
    // Calculate
    //
    rmatrixbdmultiplybyq(qp, m, n, tauq, q, m, qcolumns, false, false);
}

static void fillidentity(ap::real_2d_array& a, int n)
{
    int i;
    int j;

    a.setbounds(0, n-1, 0, n-1);
    for(i = 0; i <= n-1; i++)
    {
        for(j = 0; j <= n-1; j++)
        {
            if( i==j )
            {
                a(i,j) = 1;
            }
            else
            {
                a(i,j) = 0;
            }
        }
    }
}

/*************************************************************************
Generation of a random uniformly distributed (Haar) orthogonal matrix

INPUT PARAMETERS:
    N   -   matrix size, N>=1
    
OUTPUT PARAMETERS:
    A   -   orthogonal NxN matrix, array[0..N-1,0..N-1]

  -- ALGLIB routine --
     04.12.2009
     Bochkanov Sergey
*************************************************************************/
void rmatrixrndorthogonal(int n, ap::real_2d_array& a)
{
    int i;
    int j;

    ap::ap_error::make_assertion(n>=1, "RMatrixRndOrthogonal: N<1!");
    a.setbounds(0, n-1, 0, n-1);
    for(i = 0; i <= n-1; i++)
    {
        for(j = 0; j <= n-1; j++)
        {
            if( i==j )
            {
                a(i,j) = 1;
            }
            else
            {
                a(i,j) = 0;
            }
        }
    }
    rmatrixrndorthogonalfromtheright(a, n, n);
}

/*************************************************************************
Generation of random NxN matrix with given condition number and norm2(A)=1

INPUT PARAMETERS:
    N   -   matrix size
    C   -   condition number (in 2-norm)

OUTPUT PARAMETERS:
    A   -   random matrix with norm2(A)=1 and cond(A)=C

  -- ALGLIB routine --
     04.12.2009
     Bochkanov Sergey
*************************************************************************/
void rmatrixrndcond(int n, double c, ap::real_2d_array& a)
{
    int i;
    int j;
    double l1;
    double l2;

    ap::ap_error::make_assertion(n>=1&&ap::fp_greater_eq(c,1), "RMatrixRndCond: N<1 or C<1!");
    a.setbounds(0, n-1, 0, n-1);
    if( n==1 )
    {
        
        //
        // special case
        //
        a(0,0) = 2*ap::randominteger(2)-1;
        return;
    }
    l1 = 0;
    l2 = log(1/c);
    for(i = 0; i <= n-1; i++)
    {
        for(j = 0; j <= n-1; j++)
        {
            a(i,j) = 0;
        }
    }
    a(0,0) = exp(l1);
    for(i = 1; i <= n-2; i++)
    {
        a(i,i) = exp(ap::randomreal()*(l2-l1)+l1);
    }
    a(n-1,n-1) = exp(l2);
    rmatrixrndorthogonalfromtheleft(a, n, n);
    rmatrixrndorthogonalfromtheright(a, n, n);
}

/*************************************************************************
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);
        }
    }
}

bool obsoletesvddecomposition(ap::real_2d_array& a,
     int m,
     int n,
     ap::real_1d_array& w,
     ap::real_2d_array& v)
{
    bool result;
    int nm;
    int minmn;
    int l;
    int k;
    int j;
    int jj;
    int its;
    int i;
    double z;
    double y;
    double x;
    double vscale;
    double s;
    double h;
    double g;
    double f;
    double c;
    double anorm;
    ap::real_1d_array rv1;
    bool flag;

    rv1.setbounds(1, n);
    w.setbounds(1, n);
    v.setbounds(1, n, 1, n);
    result = true;
    if( m<n )
    {
        minmn = m;
    }
    else
    {
        minmn = n;
    }
    g = 0.0;
    vscale = 0.0;
    anorm = 0.0;
    for(i = 1; i <= n; i++)
    {
        l = i+1;
        rv1(i) = vscale*g;
        g = 0;
        s = 0;
        vscale = 0;
        if( i<=m )
        {
            for(k = i; k <= m; k++)
            {
                vscale = vscale+fabs(a(k,i));
            }
            if( ap::fp_neq(vscale,0.0) )
            {
                for(k = i; k <= m; k++)
                {
                    a(k,i) = a(k,i)/vscale;
                    s = s+a(k,i)*a(k,i);
                }
                f = a(i,i);
                g = -extsign(sqrt(s), f);
                h = f*g-s;
                a(i,i) = f-g;
                if( i!=n )
                {
                    for(j = l; j <= n; j++)
                    {
                        s = 0.0;
                        for(k = i; k <= m; k++)
                        {
                            s = s+a(k,i)*a(k,j);
                        }
                        f = s/h;
                        for(k = i; k <= m; k++)
                        {
                            a(k,j) = a(k,j)+f*a(k,i);
                        }
                    }
                }
                for(k = i; k <= m; k++)
                {
                    a(k,i) = vscale*a(k,i);
                }
            }
        }
        w(i) = vscale*g;
        g = 0.0;
        s = 0.0;
        vscale = 0.0;
        if( i<=m&&i!=n )
        {
            for(k = l; k <= n; k++)
            {
                vscale = vscale+fabs(a(i,k));
            }
            if( ap::fp_neq(vscale,0.0) )
//.........这里部分代码省略.........