c++ - 矩阵库是否有顺序中性？

Question

抱歉，如果时间太长，但是我认为这个问题需要澄清：

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


对于工作表函数，本机Excel-C接口将C以行优先顺序（矩阵的FP12类型）传递给C
Excel SDK）；
对于vba，矩阵是LPSAFEARRAY，其数据以列大顺序组织。


对于一维数据（矢量），有一个可以追溯到BLAS（30年前）的解决方案，可以将其转换为
C具有类似的结构

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

换句话说，我们在计算中使用了一个中间结构，该结构不拥有数据，但可以映射两个连续的结构
数据或以恒定间隔线性排列的数据（步幅）。 Struct的数据可以顺序处理，但是可以
还转换为数组部分（如果使用cilk）：

data [ 0 : n : stride ]

当我们从向量移到矩阵时，我已经读到可以使用
行跨步和列跨步。我的naif尝试可能是：

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开始，但是即使这样，也可以使用（可修改的）0-1参数来抽象
硬编码1的位置。all（）函数/ FULL（）宏仅对整个连续矩阵正确，而不对
一个子矩阵。但这也可以调整，在这种情况下，切换到行上的循环。与上面的copyFrom方法（）大致相同，该方法可以在行为主和列主表示之间转换（复制）矩阵。

还要注意TranspInPlace方法：通过交换行/列和行跨度/同行距，我们可以对相同的原始数据进行逻辑转换，这意味着不再需要将逻辑开关传递给通用例程（例如，GEMM用于
矩阵乘法，即使（公平地说）这并不涵盖共轭转置的情况）。

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

GSL
Gsl使用行优先顺序。停止。

拉帕克
本机代码是专栏主要的。 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向量化功能。