C++ ap::real_2d_array类代码示例

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

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

void qrdecomposition(ap::real_2d_array& a,
     int m,
     int n,
     ap::real_1d_array& tau)
{
    ap::real_1d_array work;
    ap::real_1d_array t;
    int i;
    int k;
    int mmip1;
    int minmn;
    double tmp;

    minmn = ap::minint(m, n);
    work.setbounds(1, n);
    t.setbounds(1, m);
    tau.setbounds(1, minmn);
    
    //
    // Test the input arguments
    //
    k = ap::minint(m, n);
    for(i = 1; i <= k; i++)
    {
        
        //
        // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
        //
        mmip1 = m-i+1;
        ap::vmove(t.getvector(1, mmip1), a.getcolumn(i, i, m));
        generatereflection(t, mmip1, tmp);
        tau(i) = tmp;
        ap::vmove(a.getcolumn(i, i, m), t.getvector(1, mmip1));
        t(1) = 1;
        if( i<n )
        {
            
            //
            // Apply H(i) to A(i:m,i+1:n) from the left
            //
            applyreflectionfromtheleft(a, tau(i), t, i, m, i+1, n, work);
        }
    }
}

/*************************************************************************
QR decomposition of a rectangular matrix of size MxN

Input parameters:
    A   -   matrix A whose indexes range within [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 and R in compact form (see below).
    Tau -   array of scalar factors which are used to form
            matrix Q. Array whose index ranges within [0.. Min(M-1,N-1)].

Matrix A is represented as A = QR, where Q is an orthogonal matrix of size
MxM, R - upper triangular (or upper trapezoid) matrix of size M x N.

The elements of matrix R are located on and above the main diagonal of
matrix A. The elements which are located in Tau array and below the main
diagonal of matrix A are used to form matrix Q as follows:

Matrix Q is represented as a product of elementary reflections

Q = H(0)*H(2)*...*H(k-1),

where k = min(m,n), and each H(i) is in the form

H(i) = 1 - tau * v * (v^T)

where tau is a scalar stored in Tau[I]; v - real vector,
so that v(0:i-1) = 0, v(i) = 1, v(i+1:m-1) stored in A(i+1:m-1,i).

  -- 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.
     Translation from FORTRAN to pseudocode (AlgoPascal)
     by Sergey Bochkanov, ALGLIB project, 2005-2007.
*************************************************************************/
void rmatrixqr(ap::real_2d_array& a, int m, int n, ap::real_1d_array& tau)
{
    ap::real_1d_array work;
    ap::real_1d_array t;
    int i;
    int k;
    int minmn;
    double tmp;

    if( m<=0||n<=0 )
    {
        return;
    }
    minmn = ap::minint(m, n);
    work.setbounds(0, n-1);
    t.setbounds(1, m);
    tau.setbounds(0, minmn-1);
    
    //
    // Test the input arguments
    //
    k = minmn;
    for(i = 0; i <= k-1; i++)
    {
        
        //
        // Generate elementary reflector H(i) to annihilate A(i+1:m,i)
        //
        ap::vmove(t.getvector(1, m-i), a.getcolumn(i, i, m-1));
        generatereflection(t, m-i, tmp);
        tau(i) = tmp;
        ap::vmove(a.getcolumn(i, i, m-1), t.getvector(1, m-i));
        t(1) = 1;
        if( i<n )
        {
            
            //
            // Apply H(i) to A(i:m-1,i+1:n-1) from the left
            //
            applyreflectionfromtheleft(a, tau(i), t, i, m-1, i+1, n-1, work);
        }
    }
}

/*************************************************************************
Obsolete 1-based subroutine
*************************************************************************/
void shermanmorrisonupdaterow(ap::real_2d_array& inva,
     int n,
     int updrow,
     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
    //
    ap::vmove(t1.getvector(1, n), inva.getcolumn(updrow, 1, n));
    
    //
    // T2 = v*InvA
    // Lambda = v * InvA * U
    //
    for(j = 1; j <= n; j++)
    {
        vt = ap::vdotproduct(v.getvector(1, n), inva.getcolumn(j, 1, n));
        t2(j) = vt;
    }
    lambda = t2(updrow);
    
    //
    // 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 updates matrix A^-1 when adding a vector to a row
of matrix A.

Input parameters:
    InvA    -   inverse of matrix A.
                Array whose indexes range within [0..N-1, 0..N-1].
    N       -   size of matrix A.
    UpdRow  -   the row of A whose vector V was added.
                0 <= Row <= N-1
    V       -   the vector to be added to a row.
                Array whose index ranges within [0..N-1].

Output parameters:
    InvA    -   inverse of modified matrix A.

  -- ALGLIB --
     Copyright 2005 by Bochkanov Sergey
*************************************************************************/
void rmatrixinvupdaterow(ap::real_2d_array& inva,
     int n,
     int updrow,
     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
    //
    ap::vmove(t1.getvector(0, n-1), inva.getcolumn(updrow, 0, n-1));
    
    //
    // T2 = v*InvA
    // Lambda = v * InvA * U
    //
    for(j = 0; j <= n-1; j++)
    {
        vt = ap::vdotproduct(v.getvector(0, n-1), inva.getcolumn(j, 0, n-1));
        t2(j) = vt;
    }
    lambda = t2(updrow);
    
    //
    // 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);
    }
}

/*************************************************************************
Tests Z*Z' against diag(1...1)
Returns absolute error.
*************************************************************************/
static double testort(const ap::real_2d_array& z, int n)
{
    double result;
    int i;
    int j;
    double v;

    result = 0;
    for(i = 0; i <= n-1; i++)
    {
        for(j = 0; j <= n-1; j++)
        {
            v = ap::vdotproduct(z.getcolumn(i, 0, n-1), z.getcolumn(j, 0, n-1));
            if( i==j )
            {
                v = v-1;
            }
            result = ap::maxreal(result, fabs(v));
        }
    }
    return result;
}

/*************************************************************************
Generate matrix with given condition number C (2-norm)
*************************************************************************/
static void rmatrixgenzero(ap::real_2d_array& a0, int n)
{
    int i;
    int j;

    a0.setlength(n, n);
    for(i = 0; i <= n-1; i++)
    {
        for(j = 0; j <= n-1; j++)
        {
            a0(i,j) = 0;
        }
    }
}

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

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

The algorithm updates matrix A^-1 when adding a number to an element
of matrix A.

Input parameters:
    InvA    -   inverse of matrix A.
                Array whose indexes range within [0..N-1, 0..N-1].
    N       -   size of matrix A.
    UpdRow  -   row where the element to be updated is stored.
    UpdColumn - column where the element to be updated is stored.
    UpdVal  -   a number to be added to the element.


Output parameters:
    InvA    -   inverse of modified matrix A.

  -- ALGLIB --
     Copyright 2005 by Bochkanov Sergey
*************************************************************************/
void rmatrixinvupdatesimple(ap::real_2d_array& inva,
     int n,
     int updrow,
     int updcolumn,
     double updval)
{
    ap::real_1d_array t1;
    ap::real_1d_array t2;
    int i;
    int j;
    double lambda;
    double vt;

    ap::ap_error::make_assertion(updrow>=0&&updrow<n, "RMatrixInvUpdateSimple: incorrect UpdRow!");
    ap::ap_error::make_assertion(updcolumn>=0&&updcolumn<n, "RMatrixInvUpdateSimple: incorrect UpdColumn!");
    t1.setbounds(0, n-1);
    t2.setbounds(0, n-1);
    
    //
    // T1 = InvA * U
    //
    ap::vmove(t1.getvector(0, n-1), inva.getcolumn(updrow, 0, n-1));
    
    //
    // T2 = v*InvA
    //
    ap::vmove(&t2(0), &inva(updcolumn, 0), ap::vlen(0,n-1));
    
    //
    // Lambda = v * InvA * U
    //
    lambda = updval*inva(updcolumn,updrow);
    
    //
    // InvA = InvA - correction
    //
    for(i = 0; i <= n-1; i++)
    {
        vt = updval*t1(i);
        vt = vt/(1+lambda);
        ap::vsub(&inva(i, 0), &t2(0), ap::vlen(0,n-1), vt);
    }
}

static void mheapresize(ap::real_2d_array& heap,
     int& heapsize,
     int newheapsize,
     int heapwidth)
{
    ap::real_2d_array tmp;
    int i;

    tmp.setlength(heapsize, heapwidth);
    for(i = 0; i <= heapsize-1; i++)
    {
        ap::vmove(&tmp(i, 0), &heap(i, 0), ap::vlen(0,heapwidth-1));
    }
    heap.setlength(newheapsize, heapwidth);
    for(i = 0; i <= heapsize-1; i++)
    {
        ap::vmove(&heap(i, 0), &tmp(i, 0), ap::vlen(0,heapwidth-1));
    }
    heapsize = newheapsize;
}