4. Matrix和Array的分块运算(block operation)

本章介绍block operation的基本操作。根据线性代数的知识，我们知道矩阵运算是可以通过划分为若干块分别实现的。 这里所说的块(block)是指Matrix或者Array中的一个长方形区域。

4.1 初识分块运算

一般Eigen通过调用.block()方法进行分块运算。生成一个从第i行j列(i,j)开始，具有p行q列(p,q)的块可以有如下的两种形式:

• dynamic-size:

          matrix.block(i,j,p,q); 
  

• fixed-size:

          matrix.block<p,q>(i,j); 
  

这两种形式在语法上是一致的，都可以用于fixed-size和dynamic-size的Matrix和Array。 不同之处只在于当分块规模较小时fixed-size版本产生的代码可以更快一些。

下面是一个用于Matrix的分块运算示例，在该例程中分别用dynamic-size和fixed-size两种方法求取并打印块。

        #include <iostream>
        #include <Eigen/Dense>

        using namespace std;
        int main()
        {
            Eigen::MatrixXf m(4,4);
            m <<  1, 2, 3, 4,
                5, 6, 7, 8,
                9,10,11,12,
               13,14,15,16;
            cout << "Block in the middle" << endl;
            cout << m.block<2,2>(1,1) << endl << endl;
            for (int i = 1; i <= 3; ++i) {
                cout << "Block of size " << i << "x" << i << endl;
                cout << m.block(0,0,i,i) << endl << endl;
            }
        } 

    Block in the middle
     6  7
    10 11

    Block of size 1x1
    1

    Block of size 2x2
    1 2
    5 6

    Block of size 3x3
     1  2  3
     5  6  7
     9 10 11

以上的例程中，block是用作右值的。block还可以用作左值进行复制运算，我们将在下一个例程中展示这种功能。 此外，下面的例程中将用于一个Array的分块运算，其工作方式与Matrix的是一样的。

        #include <iostream>
        #include <Eigen/Dense>

        using namespace std;
        using namespace Eigen;
        int main()
        {
            Array22f m;
            m << 1,2,
               3,4;
            Array44f a = Array44f::Constant(0.6);
            cout << "Here is the array a:" << endl << a << endl << endl;
            a.block<2,2>(1,1) = m;
            cout << "Here is now a with m copied into its central 2x2 block:" << endl << a << endl << endl;
            a.block(0,0,2,3) = a.block(2,1,2,3);
            cout << "Here is now a with bottom-right 2x3 block copied into top-left 2x2 block:" << endl << a << endl << endl;
        } 

其输出结果如下：

        Here is the array a:
        0.6 0.6 0.6 0.6
        0.6 0.6 0.6 0.6
        0.6 0.6 0.6 0.6
        0.6 0.6 0.6 0.6

        Here is now a with m copied into its central 2x2 block:
        0.6 0.6 0.6 0.6
        0.6   1   2 0.6
        0.6   3   4 0.6
        0.6 0.6 0.6 0.6

        Here is now a with bottom-right 2x3 block copied into top-left 2x2 block:
          3   4 0.6 0.6
        0.6 0.6 0.6 0.6
        0.6   3   4 0.6
        0.6 0.6 0.6 0.6

虽然.block()可以用于所有与分块运算相关的操作，但是Eigen还提供了一些更为高效的API。 例如，如果我们只需要访问矩阵的第j列，可以通过.block(0,j,p,1)来实现，其中p为矩阵的总行数。 Eigen推荐我们使用.col()，它将给编译器提供更多的优化空间，使得编译出的代码可以更高效。

独立的访问行或者列是两种特殊的分块操作。为此，Eigen专门提供了.row()和.col()。 下面的例程将介绍如何使用这两个函数，而且也显示了分块运算也可用于算术运算。

        #include <iostream>
        #include <Eigen/Dense>

        using namespace std;
        int main()
        {
            Eigen::MatrixXf m(3,3);
            m << 1,2,3,
               4,5,6,
               7,8,9;
            cout << "Here is the matrix m:" << endl << m << endl;
            cout << "2nd Row: " << m.row(1) << endl;
            m.col(2) += 3 * m.col(0);
            cout << "After adding 3 times the first column into the third column, the matrix m is:\n";
            cout << m << endl;
        } 

输出结果如下：

        Here is the matrix m:
        1 2 3
        4 5 6
        7 8 9
        2nd Row: 4 5 6
        After adding 3 times the first column into the third column, the matrix m is:
         1  2  6
         4  5 18
         7  8 30

Eigen也提供一些特殊的函数用于访问Matrix或者Array的Corner-related块。如下：

matrix.topLeftCorner(p,q); 

matrix.topLeftCorner<p,q>(); 

matrix.bottomLeftCorner(p,q); 

matrix.bottomLeftCorner<p,q>(); 

matrix.topRightCorner(p,q); 

matrix.topRightCorner<p,q>(); 

matrix.bottomRightCorner(p,q); 

matrix.bottomRightCorner<p,q>(); 

matrix.topRows(p,q); 

matrix.topRows<p,q>(); 

matrix.bottomRows(p,q); 

matrix.bottomRows<p,q>(); 

matrix.leftCols(p,q); 

matrix.leftCols<p,q>(); 

matrix.rightCols(p,q); 

matrix.rightCols<p,q>(); 

参考如下例程：

        #include <iostream>
        #include <Eigen/Dense>

        using namespace std;
        int main()
        {
            Eigen::Matrix4f m;
            m << 1, 2, 3, 4,
               5, 6, 7, 8,
               9, 10,11,12,
               13,14,15,16;
            cout << "m.leftCols(2) =" << endl << m.leftCols(2) << endl << endl;
            cout << "m.bottomRows<2>() =" << endl << m.bottomRows<2>() << endl << endl;
            m.topLeftCorner(1,3) = m.bottomRightCorner(3,1).transpose();
            cout << "After assignment, m = " << endl << m << endl;
        } 

输出结果如下：

        m.leftCols(2) =
         1  2
         5  6
         9 10
        13 14

        m.bottomRows<2>() =
         9 10 11 12
        13 14 15 16

        After assignment, m = 
         8 12 16  4
         5  6  7  8
         9 10 11 12
        13 14 15 16

4.2 向量的分块运算

vector.head(p,q); 

vector.head<p,q>(); 

vector.tail(p,q); 

vector.tail<p,q>(); 

vector.segment(p,q); 

vector.segment<p,q>(); 

参考如下例程：

        #include <iostream>
        #include <Eigen/Dense>

        using namespace std;
        int main()
        {
            Eigen::ArrayXf v(6);
            v << 1, 2, 3, 4, 5, 6;
            cout << "v.head(3) =" << endl << v.head(3) << endl << endl;
            cout << "v.tail<3>() = " << endl << v.tail<3>() << endl << endl;
            v.segment(1,4) *= 2;
            cout << "after 'v.segment(1,4) *= 2', v =" << endl << v << endl;
        } 

        v.head(3) =
        1
        2
        3

        v.tail<3>() = 
        4
        5
        6

        after 'v.segment(1,4) *= 2', v =
         1
         4
         6
         8
        10
         6