任何矩阵库顺序无关

Any matrix library order-neutral?

本文关键字：顺序 任何矩      更新时间：2023-10-16

抱歉，如果这太长了，但我觉得这个问题需要澄清:

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

  • 对于表格函数，本机Excel-C接口将按行为主顺序(类型为FP12)的矩阵传递给CExcel SDK);

  • 对于vba，矩阵是一个LPSAFEARRAY，它将数据按列主序组织。

对于一维数据(向量)，有一个可以追溯到BLAS(30年前)的解决方案，可以翻译为C在

这样的结构中
struct VECTOR { 
    int n;
    int stride;
    double * data;
    double & operator [] (int i) { return data[(i - 1)*stride]; }   
}

换句话说，我们使用一个中间结构进行计算，它不拥有数据，但可以映射两个相邻的结构数据或数据线性间隔与一个恒定的间隙(步长)。Struct的数据可以按顺序处理，但是它们可以也可以转换为数组部分(如果使用click):

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开始，但即使这也可以使用一个(可修改的)的0-1参数抽象硬编码1的位置。all()函数/FULL()宏只对一个完整的、连续的矩阵是正确的，而不是子矩阵。但这也可以调整，在这种情况下切换到行上的循环。或多或少类似于上面的copyFrom方法()，它可以在行主要表示和列主要表示之间转换(复制)矩阵。

还请注意TranspInPlace方法:通过交换行/cols和rowstride/colstride，我们对相同的未动数据进行了逻辑转换，这意味着不再需要将逻辑转换传递给通用例程(例如，为矩阵乘法，即使(公平地说)这并不包括共轭转置的情况。

鉴于上述情况，我正在寻找一个实现这种方法的库，以便我可以使用(hook)它，但我的搜索到现在还不满意:

GSLGsl使用行主顺序。停止。

LAPACK本机代码是列为主的。c接口提供了处理行为主数据的可能性，但只能添加定制的数据(在某些例程中)在对矩阵进行操作之前对其进行物理转置。

特征作为一个模板化的库，它可以处理这两种情况。不幸的是，矩阵的顺序是模板的一个参数意味着每个矩阵方法都是重复的。它可以工作，但不是理想的。

请让我知道如果图书馆更接近我所追求的。Thx .

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

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

然后可以对mat进行任何操作，如访问mat.row(i)行，mat.col(j)列，进行乘积，求解线性方程组等。

这种方法的主要缺点是失去了SIMD向量化。

相关文章：