SLISC 密矩阵的切割

                     

贡献者: addis

  • 本文处于草稿阶段。
预备知识 SLISC 的矢量和矩阵

   如果我们想把矩阵的一部分(例如一行)输入到某个函数中怎么办呢?如果我们不使用任何容器而是用指针,函数的输入输出也使用指针,那么可以通过 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 &,用于输出。

  

未完成:具体介绍一下所有类型

                     

© 小时科技 保留一切权利