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

bool generalizedsymmetricdefiniteevdreduce(ap::real_2d_array& a,
        int n,
        bool isuppera,
        const ap::real_2d_array& b,
        bool isupperb,
        int problemtype,
        ap::real_2d_array& r,
        bool& isupperr)
{
    bool result;
    ap::real_2d_array t;
    ap::real_1d_array w1;
    ap::real_1d_array w2;
    ap::real_1d_array w3;
    int i;
    int j;
    double v;

    ap::ap_error::make_assertion(n>0, "GeneralizedSymmetricDefiniteEVDReduce: N<=0!");
    ap::ap_error::make_assertion(problemtype==1||problemtype==2||problemtype==3, "GeneralizedSymmetricDefiniteEVDReduce: incorrect ProblemType!");
    result = true;

    //
    // Problem 1:  A*x = lambda*B*x
    //
    // Reducing to:
    //     C*y = lambda*y
    //     C = L^(-1) * A * L^(-T)
    //     x = L^(-T) * y
    //
    if( problemtype==1 )
    {

        //
        // Factorize B in T: B = LL'
        //
        t.setbounds(1, n, 1, n);
        if( isupperb )
        {
            for(i = 1; i <= n; i++)
            {
                ap::vmove(t.getcolumn(i, i, n), b.getrow(i, i, n));
            }
        }
        else
        {
            for(i = 1; i <= n; i++)
            {
                ap::vmove(&t(i, 1), &b(i, 1), ap::vlen(1,i));
            }
        }
        if( !choleskydecomposition(t, n, false) )
        {
            result = false;
            return result;
        }

        //
        // Invert L in T
        //
        if( !invtriangular(t, n, false, false) )
        {
            result = false;
            return result;
        }

        //
        // Build L^(-1) * A * L^(-T) in R
        //
        w1.setbounds(1, n);
        w2.setbounds(1, n);
        r.setbounds(1, n, 1, n);
        for(j = 1; j <= n; j++)
        {

            //
            // Form w2 = A * l'(j) (here l'(j) is j-th column of L^(-T))
            //
            ap::vmove(&w1(1), &t(j, 1), ap::vlen(1,j));
            symmetricmatrixvectormultiply(a, isuppera, 1, j, w1, 1.0, w2);
            if( isuppera )
            {
                matrixvectormultiply(a, 1, j, j+1, n, true, w1, 1, j, 1.0, w2, j+1, n, 0.0);
            }
            else
            {
                matrixvectormultiply(a, j+1, n, 1, j, false, w1, 1, j, 1.0, w2, j+1, n, 0.0);
            }

            //
            // Form l(i)*w2 (here l(i) is i-th row of L^(-1))
            //
            for(i = 1; i <= n; i++)
            {
                v = ap::vdotproduct(&t(i, 1), &w2(1), ap::vlen(1,i));
                r(i,j) = v;
            }
        }

        //
//.........这里部分代码省略.........

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

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

Входные параметры:
    A       -   матрица A. Нумерация элементов: [1..M, 1..N]
    M       -   число строк в матрице A
    N       -   число столбцов в матрице A

Выходные параметры:
    A       -   матрицы L и U в компактной форме (см. ниже).
                Нумерация элементов: [1..M, 1..N]
    Pivots  -   матрица перестановок в компактной форме (см. ниже).
                Нумерация элементов: [1..Min(M,N)]

Матрица A представляется, как A = P * L * U, где P - матрица перестановок,
матрица L - нижнетреугольная (или нижнетрапецоидальная, если M>N) матрица,
U - верхнетреугольная (или верхнетрапецоидальная, если M<N) матрица.

Рассмотрим разложение более подробно на примере при M=4, N=3:

                   (  1          )    ( U11 U12 U13  )
A = P1 * P2 * P3 * ( L21  1      )  * (     U22 U23  )
                   ( L31 L32  1  )    (         U33  )
                   ( L41 L42 L43 )

Здесь матрица L  имеет  размер  M  x  Min(M,N),  матрица  U  имеет  размер
Min(M,N) x N, матрица  P(i)  получается  путем  перестановки  в  единичной
матрице размером M x M строк с номерами I и Pivots[I]

Результатом работы алгоритма являются массив Pivots  и  следующая матрица,
замещающая  матрицу  A,  и  сохраняющая  в компактной форме матрицы L и U
(пример приведен для M=4, N=3):

 ( U11 U12 U13 )
 ( L21 U22 U23 )
 ( L31 L32 U33 )
 ( L41 L42 L43 )

Как видно, единичная диагональ матрицы L  не  сохраняется.
Если N>M, то соответственно меняются размеры матриц и расположение
элементов.

  -- LAPACK routine (version 3.0) --
     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
     Courant Institute, Argonne National Lab, and Rice University
     June 30, 1992
*************************************************************************/
void ludecomposition(ap::real_2d_array& a,
                     int m,
                     int n,
                     ap::integer_1d_array& pivots)
{
    int i;
    int j;
    int jp;
    ap::real_1d_array t1;
    double s;

    pivots.setbounds(1, ap::minint(m, n));
    t1.setbounds(1, ap::maxint(m, n));
    ap::ap_error::make_assertion(m>=0&&n>=0);

    //
    // Quick return if possible
    //
    if( m==0||n==0 )
    {
        return;
    }
    for(j = 1; j <= ap::minint(m, n); j++)
    {

        //
        // Find pivot and test for singularity.
        //
        jp = j;
        for(i = j+1; i <= m; i++)
        {
            if( fabs(a(i,j))>fabs(a(jp,j)) )
            {
                jp = i;
            }
        }
        pivots(j) = jp;
        if( a(jp,j)!=0 )
        {

            //
            //Apply the interchange to rows
            //
            if( jp!=j )
            {
                ap::vmove(t1.getvector(1, n), a.getrow(j, 1, n));
                ap::vmove(a.getrow(j, 1, n), a.getrow(jp, 1, n));
                ap::vmove(a.getrow(jp, 1, n), t1.getvector(1, n));
            }

//.........这里部分代码省略.........

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

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

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

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

Входные параметры:
    A           -   матрица. Массив с нумерацией элементов [1..N,1..N]
    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 invtriangular(ap::real_2d_array& a,
                   int n,
                   bool isupper,
                   bool isunittriangular)
{
    bool result;
    bool nounit;
    int i;
    int j;
    int nmj;
    int jm1;
    int jp1;
    double v;
    double ajj;
    ap::real_1d_array t;

    result = true;
    t.setbounds(1, n);

    //
    // Test the input parameters.
    //
    nounit = !isunittriangular;
    if( isupper )
    {

        //
        // Compute inverse of upper triangular matrix.
        //
        for(j = 1; j <= n; 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>1 )
            {
                jm1 = j-1;
                ap::vmove(t.getvector(1, jm1), a.getcolumn(j, 1, jm1));
                for(i = 1; i <= j-1; i++)
                {
                    if( i<j-1 )
                    {
                        v = ap::vdotproduct(a.getrow(i, i+1, jm1), t.getvector(i+1, jm1));
                    }
                    else
                    {
                        v = 0;
                    }
//.........这里部分代码省略.........

//************************************************************************
//Обращение матрицы, заданной LU-разложением
//
//Входные параметры:
//    A       -   LU-разложение  матрицы   (результат   работы  подпрограммы
//                LUDecomposition).
//    Pivots  -   таблица перестановок,  произведенных в ходе LU-разложения.
//                (результат работы подпрограммы LUDecomposition).
//    N       -   размерность матрицы
//
//Выходные параметры:
//    A       -   матрица, обратная к исходной. Массив с нумерацией
//                элементов [1..N, 1..N]
//
//Результат:
//    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 inverselu(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;
    int jp1;
    double v;

    result = true;

    //
    // Quick return if possible
    //
    if( n==0 )
    {
        return result;
    }
    work.setbounds(1, n);

    //
    // Form inv(U)
    //
    if( !invtriangular(a, n, true, false) )
    {
        result = false;
        return result;
    }

    //
    // Solve the equation inv(A)*L = inv(U) for inv(A).
    //
    for(j = n; j >= 1; j--)
    {
        //
        // Copy current column of L to WORK and replace with zeros.
        //
        for(i = j+1; i <= n; i++)
        {
            work(i) = a(i,j);
            a(i,j) = 0;
        }

        //
        // Compute current column of inv(A).
        //
        if( j<n )
        {
            jp1 = j+1;
            for(i = 1; i <= n; i++)
            {
                v = ap::vdotproduct(a.getrow(i, jp1, n), work.getvector(jp1, n));
                a(i,j) = a(i,j)-v;
            }
        }
    }

    //
    // Apply column interchanges.
    //
    for(j = n-1; j >= 1; j--)
    {
        jp = pivots(j);
        if( jp!=j )
        {
            ap::vmove(work.getvector(1, n), a.getcolumn(j, 1, n));
            ap::vmove(a.getcolumn(j, 1, n), a.getcolumn(jp, 1, n));
            ap::vmove(a.getcolumn(jp, 1, n), work.getvector(1, n));
        }
    }
    return result;
//.........这里部分代码省略.........

static void naivematrixmatrixmultiply(const ap::real_2d_array& a,
     int ai1,
     int ai2,
     int aj1,
     int aj2,
     bool transa,
     const ap::real_2d_array& b,
     int bi1,
     int bi2,
     int bj1,
     int bj2,
     bool transb,
     double alpha,
     ap::real_2d_array& c,
     int ci1,
     int ci2,
     int cj1,
     int cj2,
     double beta)
{
    int arows;
    int acols;
    int brows;
    int bcols;
    int i;
    int j;
    int k;
    int l;
    int r;
    double v;
    ap::real_1d_array x1;
    ap::real_1d_array x2;

    
    //
    // Setup
    //
    if( !transa )
    {
        arows = ai2-ai1+1;
        acols = aj2-aj1+1;
    }
    else
    {
        arows = aj2-aj1+1;
        acols = ai2-ai1+1;
    }
    if( !transb )
    {
        brows = bi2-bi1+1;
        bcols = bj2-bj1+1;
    }
    else
    {
        brows = bj2-bj1+1;
        bcols = bi2-bi1+1;
    }
    ap::ap_error::make_assertion(acols==brows, "NaiveMatrixMatrixMultiply: incorrect matrix sizes!");
    if( arows<=0||acols<=0||brows<=0||bcols<=0 )
    {
        return;
    }
    l = arows;
    r = bcols;
    k = acols;
    x1.setbounds(1, k);
    x2.setbounds(1, k);
    for(i = 1; i <= l; i++)
    {
        for(j = 1; j <= r; j++)
        {
            if( !transa )
            {
                if( !transb )
                {
                    v = ap::vdotproduct(b.getcolumn(bj1+j-1, bi1, bi2), a.getrow(ai1+i-1, aj1, aj2));
                }
                else
                {
                    v = ap::vdotproduct(&b(bi1+j-1, bj1), &a(ai1+i-1, aj1), ap::vlen(bj1,bj2));
                }
            }
            else
            {
                if( !transb )
                {
                    v = ap::vdotproduct(b.getcolumn(bj1+j-1, bi1, bi2), a.getcolumn(aj1+i-1, ai1, ai2));
                }
                else
                {
                    v = ap::vdotproduct(b.getrow(bi1+j-1, bj1, bj2), a.getcolumn(aj1+i-1, ai1, ai2));
                }
            }
            if( ap::fp_eq(beta,0) )
            {
                c(ci1+i-1,cj1+j-1) = alpha*v;
            }
            else
            {
                c(ci1+i-1,cj1+j-1) = beta*c(ci1+i-1,cj1+j-1)+alpha*v;
//.........这里部分代码省略.........