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

void complexludecompositionunpacked(ap::complex_2d_array a,
     int m,
     int n,
     ap::complex_2d_array& l,
     ap::complex_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);
    complexludecomposition(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);
            }
        }
    }
}

/*************************************************************************
Generation of random NxN Hermitian 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 HPD matrix with norm2(A)=1 and cond(A)=C

  -- ALGLIB routine --
     04.12.2009
     Bochkanov Sergey
*************************************************************************/
void hpdmatrixrndcond(int n, double c, ap::complex_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
    //
    hmatrixrndmultiply(a, n);
}

/*************************************************************************
Copy
*************************************************************************/
static void makeacopy(const ap::complex_2d_array& a,
     int m,
     int n,
     ap::complex_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 Hermitian 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 hmatrixrndcond(int n, double c, ap::complex_2d_array& a)
{
    int i;
    int j;
    double l1;
    double l2;

    ap::ap_error::make_assertion(n>=1&&ap::fp_greater_eq(c,1), "HMatrixRndCond: 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
    //
    hmatrixrndmultiply(a, n);
}

/*************************************************************************
Generation of random NxN complex matrix with given condition number C 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 cmatrixrndcond(int n, double c, ap::complex_2d_array& a)
{
    int i;
    int j;
    double l1;
    double l2;
    hqrndstate state;
    ap::complex v;

    ap::ap_error::make_assertion(n>=1&&ap::fp_greater_eq(c,1), "CMatrixRndCond: N<1 or C<1!");
    a.setbounds(0, n-1, 0, n-1);
    if( n==1 )
    {
        
        //
        // special case
        //
        hqrndrandomize(state);
        hqrndunit2(state, v.x, v.y);
        a(0,0) = v;
        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);
    cmatrixrndorthogonalfromtheleft(a, n, n);
    cmatrixrndorthogonalfromtheright(a, n, n);
}

/*************************************************************************
Generation of a random Haar distributed orthogonal complex 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 cmatrixrndorthogonal(int n, ap::complex_2d_array& a)
{
    int i;
    int j;

    ap::ap_error::make_assertion(n>=1, "CMatrixRndOrthogonal: 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;
            }
        }
    }
    cmatrixrndorthogonalfromtheright(a, n, n);
}

/*************************************************************************
Unsets 2D array.
*************************************************************************/
static void unset2dc(ap::complex_2d_array& a)
{

    a.setbounds(0, 0, 0, 0);
    a(0,0) = 2*ap::randomreal()-1;
}

/*************************************************************************

  -- ALGLIB --
     Copyright 2005, 2007 by Bochkanov Sergey
*************************************************************************/
void unpackqfromhermitiantridiagonal(const ap::complex_2d_array& a,
     const int& n,
     const bool& isupper,
     const ap::complex_1d_array& tau,
     ap::complex_2d_array& q)
{
    int i;
    int j;
    ap::complex_1d_array v;
    ap::complex_1d_array work;
    int i_;
    int i1_;

    if( n==0 )
    {
        return;
    }
    
    //
    // init
    //
    q.setbounds(1, n, 1, n);
    v.setbounds(1, n);
    work.setbounds(1, n);
    for(i = 1; i <= n; i++)
    {
        for(j = 1; j <= n; j++)
        {
            if( i==j )
            {
                q(i,j) = 1;
            }
            else
            {
                q(i,j) = 0;
            }
        }
    }
    
    //
    // unpack Q
    //
    if( isupper )
    {
        for(i = 1; i <= n-1; i++)
        {
            
            //
            // Apply H(i)
            //
            for(i_=1; i_<=i;i_++)
            {
                v(i_) = a(i_,i+1);
            }
            v(i) = 1;
            complexapplyreflectionfromtheleft(q, tau(i), v, 1, i, 1, n, work);
        }
    }
    else
    {
        for(i = n-1; i >= 1; i--)
        {
            
            //
            // Apply H(i)
            //
            i1_ = (i+1) - (1);
            for(i_=1; i_<=n-i;i_++)
            {
                v(i_) = a(i_+i1_,i);
            }
            v(1) = 1;
            complexapplyreflectionfromtheleft(q, tau(i), v, i+1, n, 1, n, work);
        }
    }
}

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

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

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

  -- ALGLIB --
     Copyright 2005, 2007, 2008 by Bochkanov Sergey
*************************************************************************/
void hmatrixtdunpackq(const ap::complex_2d_array& a,
     const int& n,
     const bool& isupper,
     const ap::complex_1d_array& tau,
     ap::complex_2d_array& q)
{
    int i;
    int j;
    ap::complex_1d_array v;
    ap::complex_1d_array work;
    int i_;
    int i1_;

    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)
            //
            i1_ = (0) - (1);
            for(i_=1; i_<=i+1;i_++)
            {
                v(i_) = a(i_+i1_,i+1);
            }
            v(i+1) = 1;
            complexapplyreflectionfromtheleft(q, tau(i), v, 0, i, 0, n-1, work);
        }
    }
    else
    {
        for(i = n-2; i >= 0; i--)
        {
            
            //
            // Apply H(i)
            //
            i1_ = (i+1) - (1);
            for(i_=1; i_<=n-i-1;i_++)
            {
                v(i_) = a(i_+i1_,i);
            }
            v(1) = 1;
            complexapplyreflectionfromtheleft(q, tau(i), v, i+1, n-1, 0, n-1, work);
        }
    }
}

/*************************************************************************
Finding the eigenvalues and eigenvectors of a Hermitian matrix

The algorithm finds eigen pairs of a Hermitian matrix by  reducing  it  to
real tridiagonal form and using the QL/QR algorithm.

Input parameters:
    A       -   Hermitian matrix which is given  by  its  upper  or  lower
                triangular part.
                Array whose indexes range within [0..N-1, 0..N-1].
    N       -   size of matrix A.
    IsUpper -   storage format.
    ZNeeded -   flag controlling whether the eigenvectors  are  needed  or
                not. If ZNeeded is equal to:
                 * 0, the eigenvectors are not returned;
                 * 1, the eigenvectors are returned.

Output parameters:
    D       -   eigenvalues in ascending order.
                Array whose index ranges within [0..N-1].
    Z       -   if ZNeeded is equal to:
                 * 0, Z hasn’t changed;
                 * 1, Z contains the eigenvectors.
                Array whose indexes range within [0..N-1, 0..N-1].
                The eigenvectors are stored in the matrix columns.

Result:
    True, if the algorithm has converged.
    False, if the algorithm hasn't converged (rare case).

Note:
    eigen vectors of Hermitian matrix are defined up to multiplication  by
    a complex number L, such as |L|=1.

  -- ALGLIB --
     Copyright 2005, 23 March 2007 by Bochkanov Sergey
*************************************************************************/
bool hmatrixevd(ap::complex_2d_array a,
     int n,
     int zneeded,
     bool isupper,
     ap::real_1d_array& d,
     ap::complex_2d_array& z)
{
    bool result;
    ap::complex_1d_array tau;
    ap::real_1d_array e;
    ap::real_1d_array work;
    ap::real_2d_array t;
    ap::complex_2d_array q;
    int i;
    int k;
    double v;

    ap::ap_error::make_assertion(zneeded==0||zneeded==1, "HermitianEVD: incorrect ZNeeded");
    
    //
    // Reduce to tridiagonal form
    //
    hmatrixtd(a, n, isupper, tau, d, e);
    if( zneeded==1 )
    {
        hmatrixtdunpackq(a, n, isupper, tau, q);
        zneeded = 2;
    }
    
    //
    // TDEVD
    //
    result = smatrixtdevd(d, e, n, zneeded, t);
    
    //
    // Eigenvectors are needed
    // Calculate Z = Q*T = Re(Q)*T + i*Im(Q)*T
    //
    if( result&&zneeded!=0 )
    {
        work.setbounds(0, n-1);
        z.setbounds(0, n-1, 0, n-1);
        for(i = 0; i <= n-1; i++)
        {
            
            //
            // Calculate real part
            //
            for(k = 0; k <= n-1; k++)
            {
                work(k) = 0;
            }
            for(k = 0; k <= n-1; k++)
            {
                v = q(i,k).x;
                ap::vadd(&work(0), &t(k, 0), ap::vlen(0,n-1), v);
            }
            for(k = 0; k <= n-1; k++)
            {
                z(i,k).x = work(k);
            }
            
            //
//.........这里部分代码省略.........