CDS 矩阵-向量乘积

如果 n \times n 矩阵 A 以 CDS 格式存储，仍然可以通过按行或按列的方式执行矩阵-向量乘积 y=Ax，但这并没有充分利用 CDS 格式的优势。其思路是在双重嵌套循环中进行坐标变换。通过替换 j\rightarrow i+j，我们得到

y_i\leftarrow y_i+a_{i,j}x_j\quad\Rightarrow\quad y_i\leftarrow y_i+a_{i,i+j}x_{i+j} ~.

在内层循环中使用索引 i，我们可以看到表达式 a_{i,i+j} 访问了矩阵的第 j 条对角线（其中主对角线的编号为 0）。

该算法现在将包含一个双重嵌套循环，外层循环枚举对角线 diag=-p,q，其中 p 和 q 分别是主对角线左侧和右侧的（非负）对角线数量。内层循环的边界由以下要求决定：

1\leq {\tt i},{\tt i+j}\leq n.

算法变为

   for i = 1, n
       y(i) = 0
   end;
   for diag = -diag_left, diag_right
       for loc = max(1,1-diag), min(n,n-diag)
           y(loc) = y(loc) + val(loc,diag) * x(loc+diag)
       end;
   end;

转置矩阵-向量乘积 y=A^Tx 是上述算法的一个微小变体。使用更新公式

\begin{aligned} y_i&\leftarrow y_i+a_{j,i}x_j \\ &=y_i+a_{j,i+j-j}x_{i+j} \end{aligned}

我们得到

   for i = 1, n
       y(i) = 0
   end;
   for diag = -diag_right, diag_left
       for loc = max(1,1-diag), min(n,n-diag)
           y(loc) = y(loc) + val(loc+diag, -diag) * x(loc+diag)
       end;
   end;

基于 CDS 的矩阵-向量乘积 y=Ax（或 y=A^Tx）的内存访问在内层迭代中涉及三个向量。另一方面，不存在间接寻址，并且该算法可通过向量长度实质上为矩阵阶数 n 的向量化实现。由于数据访问的规律性，大多数机器可以通过保持三个基址寄存器并使用简单的偏移寻址来高效执行此算法。