0

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

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

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

  • 对于 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 ]

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

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 接口提供了处理主要行数据的可能性，但只需添加定制代码或（在某些例程中）在对矩阵进行操作之前物理转置矩阵。

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

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

4

1 回答 1

1

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

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

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

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

于 2016-11-05T13:23:54.273 回答