任何矩阵库顺序中性？ [关闭]

Question

对不起，如果这太长了，但我觉得这个问题需要澄清：

我正在为 Excel 开发一个 xll 库，即一个 C 库，其中包含可以在单元格中直接注册和调用的函数。理想情况下，这些函数也应该从 VBA 调用（或适应调用）， 为了为更复杂的计算（求根、颂歌、优化）提供解释环境，这些计算不适合 Excel 工作表。需要明确的是：有一种方法可以从 vba 调用工作表函数（函数 Application.Run），但它为所有参数和返回值从/到变体类型的转换付出了不可接受的代价。现在有趣的情况是，在同一个 Excel 环境中，二维矩阵以不同的方式传递：

• 对于工作表函数，本机 Excel-C 接口以行优先顺序向 C 传递一个矩阵（类型为 FP12 Excel SDK）；

• 对于 vba，矩阵是 LPSAFEARRAY，其数据按列优先顺序组织。

对于一维数据（向量），有一个可以追溯到 BLAS（30 年前）的解决方案，可以翻译为 C in 具有类似的结构

struct VECTOR { 
    int n;
    int stride;
    double * data;
    double & operator [] (int i) { return data[(i - 1)*stride]; }   
}   

换句话说，我们为计算使用了一个中间结构，它不拥有数据，但可以映射两个连续的 数据或以恒定间隙（步幅）线性间隔的数据。 struct 的数据可以按顺序处理，但它们可以 也翻译成数组部分（如果使用 cilk）：

data [ 0 : n : stride ]

当我们从向量转移到矩阵时，我已经读到可以从矩阵顺序中抽象使用 行跨度和列跨度。我的天真尝试可能是：

struct MATRIX {
    int rows;
    int cols;
    int rowstride;
    int colstride;
    double * data;
    inline double & operator () (int i, int j)  { return data[(i - 1)*rowstride + (j - 1)*colstride]; }
    MATRIX(int nrows, int ncols, int incrw, int inccl, double* dt) {rows = nrows; cols = ncols, rowstride = incrw; colstride = inccl; data = dt; }
    MATRIX(FP12 & A)        { rows = A.rows;  cols = A.cols;  data = A.array; rowstride = cols; colstride = 1;  }
    MATRIX(LPSAFEARRAY & x) { rows = ROWS(x); cols = COLS(x); data = DATA(x); rowstride = 1;    colstride = rows; }
    int els() { return rows * cols; }
    bool isRowMajor() { return rowstride > 1; }
    bool isScalar()   { return (rows == 1) & (cols == 1); }
    bool isRow()      { return (rows == 1); }
    bool isCol()      { return (cols == 1); }
    VECTOR col(int i) { return {rows, rowstride, &data[(i - 1)*colstride] }; }      // Col(1..)
    VECTOR row(int i) { return {cols, colstride, &data[(i - 1)*rowstride] }; }      // Row(1..)
    VECTOR all()      { return {els(), 1, data}; }  
    void copyFrom   (MATRIX & B) { for (int i = 1; i <= rows; i++) ROW(*this, i) = ROW(B, i); }
    MATRIX & Transp (MATRIX & B) { for (int i = 1; i <= rows; i++) ROW(*this, i) = COL(B, i); return *this; }
    void BinaryOp   (BinaryFcn f, MATRIX & B);
    MATRIX TranspInPlace() { return MATRIX(cols, rows, colstride, rowstride, data); }
    MATRIX SubMatrix(int irow, int icol, int nrows, int ncols) { return MATRIX(nrows, ncols, rowstride, colstride, &(*this)(irow, icol)); }
};  

来自 FP12 或 LPSAFEARRAY 的两个构造函数初始化指向按行优先或按列组织的数据的结构。这是订单中性吗？在我看来是的：数据访问（索引）和行/列选择都是正确的，独立于顺序。考虑到两次乘法，索引速度较慢，但​​可以非常快速地映射行和列：毕竟矩阵库的目的是最小化单个数据访问。此外，它是 很容易编写宏来为行或列返回数组部分，也为整个矩阵：

#define COL(A,i) (A).data[(i-1)*(A).colstride : (A).rows : (A).rowstride]           // COL(A,1)
#define ROW(A,i) (A).data[(i-1)*(A).rowstride : (A).cols : (A).colstride]           // ROW(A,1)
#define FULL(A)  (A).data[0 : (A).rows * (A).cols]                                  // FULL MATRIX

在上面的代码中，索引从 1 开始，但即使这样也可以使用（可修改的）ofs 0-1 参数在 硬编码的位置 1. all() 函数/FULL() 宏仅对整个连续矩阵是正确的，而不是 一个子矩阵。但是这也可以调整，在这种情况下切换到行上的循环。或多或少类似于上面的 copyFrom 方法（），它可以在行优先和列优先表示之间转换（复制）矩阵。

还要注意 TranspInPlace 方法：通过交换 rows/cols 和 rowstride/colstride，我们可以对相同的未触及数据进行逻辑转置，这意味着不再需要将逻辑开关传递给通用例程（例如 GEMM为了 矩阵乘法，尽管（公平地说）这不包括共轭转置的情况）。

鉴于上述情况，我正在寻找一个实现这种方法的库，以便我可以使用（挂钩）它，但我的搜索到目前为止并不令人满意：

GSL Gsl 使用行优先顺序。停下来。

LAPACK 本机代码是列优先的。 C 接口提供了处理主要行数据的可能性，但只能添加定制的 代码或（在某些例程中）在对矩阵进行操作之前对矩阵进行物理转置。

本征 作为一个模板库，它可以同时处理两者。不幸的是，矩阵顺序是模板的一个参数，它 意味着每个矩阵方法都将被复制。它有效，但并不理想。

如果图书馆更接近我所追求的，请告诉我。谢谢。

Answer 1

在 Eigen 中，您可以使用 Map 来模拟这一点，其中包含运行时内部和外部跨步。例如，如果您坚持使用ColumnMajor，则内部跨度对应于行跨度，外部跨度对应于列跨度：

Map<MatrixXd,0,Stride<> > mat(ptr, rows, cols, Stride<>(colStride,rowStride));

然后你可以对mat进行任何操作，比如访问mat.row(i)的行、mat.col(j)的列、执行产品、求解线性系统等。

这种方法的主要缺点是您松散了 SIMD 矢量化。