SLISC 密矩阵的切割

                     

贡献者： addis

　　 如果我们想把矩阵的一部分（例如一行）输入到某个函数中怎么办呢？如果我们不使用任何容器而是用指针，函数的输入输出也使用指针，那么可以通过 step 参数实现：

// 矢量求和
Doub sum(Doub_I *x, Long_I N, Long_I step)
{
    Doub s = 0;
    for (Long i = 0; i < N; ++i)
        s += x[step*i];
    return s;
}

Long N0 = 10, N1 = 10; row = 3;
Doub *a; // 列主序矩阵
new Doub a[N0*N1];
for (Long i = 0; i < N0*N1; ++i)
    a[i] = i;
// 对第 3 行求和
cout << "sum a(3,all) = " << sum(a+3, N1, N0) << endl;
// 对第 3 列求和
cout << "sum a(all,3) = " << sum(a+3*N0, N0, 1) << endl;
// 对第 3 列的前 3 个元求和
cout << "sum(0-2,3) = " << sum(a+3*N0, 3, 1) << endl;
delete a[];

再看一个例子，如何对子矩阵求和呢？我们需要一个参数叫 leading dimension，该参数表示列标（行主序时：行标）增加 1 时单索引增加的数量。

// N0, N1 时子矩阵的尺寸， lda 是原矩阵的行数
Doub sum(Doub_I *a, Long_I N0, Long_I N1, Long_I lda)
{
    Doub s = 0;
    for (Long j = 0; i < N1; ++i) {
        for (Long i = 0; i < N0; ++i)
            s += a[i];
        a += lda;
    }
    return s;
}

使用方法：

Long N0 = 10, N1 = 10;
Doub *a;
new Doub a[N0*N1]; // 列主序
for (Long i = 0; i < N0*N1; ++i)
    a[i] = i;
// 对 3-5 行， 4-6 列的子矩阵求和
cout << "sum(3-5,4-6) = " << sum(a+4*N0+3, 3, 3, N0) << endl;

可见直接使用裸指针给函数传递数组是非常灵活的，在 C 语言中这也是很常见的做法，见 BLAS 和 LAPACK。然而这么做的缺点是容易出错，例如索引超出分配的内存，或者忘记用 delete 释放内存。这样的错误调试起来会非常困难。

1. view 类型

　　 在 SLISC 库中，函数的输入或输出除了各种容器外，也可以是 view 类型的一种。例如上面的 sum 函数也可以定义为

Doub sum(DvecDoub_I x)
{
    Doub s = x[0];
    for (Long i = 1; i < x.size(); ++i)
        s += x[x.step()*i];
    return s;
}

其中 DvecDoub 是 view 类型的一种，它包含一个指针 Doub *m_p，一个长度 Long *m_N 以及一个步长 Long *m_step。它生成和毁灭时不会分配或释放内存，它只用于表示另一个容器的一部分数据，所以叫做 view，也叫 slicing。

　　 SLISC 提供一系列开头为 cut 的函数，可以从容器中截取一部分然后返回 view 类型的一种。例如以上 sum 函数可以这样使用

CmatDoub a(10, 10);
// 对第 3 行求和
cout << "sum row 3 = " << sum(cut1(a, 3)) << endl;

其中 cut1(a, 3) 从矩阵 a 中切下第 3 行（1 表示第 1 个维度，即 “行”），然后返回一个临时的 const DvecDoub & 类型传给 sum 函数。

　　 类似地，cut0(a, 3) 从矩阵 a 中切下第 3 列（0 表示第 0 个维度，即 “列”），然后返回一个临时的 const SvecDoub & 类型。SvecDoub 是一个密矢量的 view，没有 m_step 成员，所有元素在内存中都是相邻的。在 view 类型命名中，S 代表 smooth，D 代表 dash。

　　 同理，如果要切割一个子矩阵，用 cut(a, i, N0, j, N1) 可以切出一个 DcmatDoub 类型，它的成员包括一个起始指针 Doub *m_p，子矩阵尺寸 Long *m_N0, *m_N1，总长度 Long m_N，以及 leading dimension m_lda。

　　 view 类型和指针一样，也分为 low level const 和 high level const，const SvecDoub 代表 high level const，也就是它一旦初始化后数据成员就都是 const，其中 m_p 是 high level，即 const Doub *。view 的 low level const 则有专门的类型如 SvecDoub_c，它的成员可以随意改变，但所代表的数据不能被修改。

　　 在函数参数中，SvecDoub_I 的定义是 const SvecDoub_c & 用于输入；SvecDoub_O, SvecDoub_IO 的定义是 const SvecDoub &，用于输出。