任何矩陣庫order-neutral？

對不起，如果這太長，但我覺得這個問題需要澄清：任何矩陣庫order-neutral？

我正在爲Excel的xll庫，即包含函數，它可以註冊，並直接調用形式細胞。理想情況下，這些函數也應該從VBA，中調用（或調整爲調用），以便爲更復雜的計算（根查找，優化，優化）提供解釋環境，這些計算不適合在Excel工作表中。要清楚：有一種方法可以從vba（函數Application.Run）調用sheet函數，但是它支付所有參數和返回值的變換類型轉換的不可接受的價格。現在有趣的情況是，同一個Excel環境中，一個二維矩陣，以不同的方式傳遞：

板材的功能，本機的Excel-C接口傳遞到C行優先的順序矩陣（Excel SDK的 FP12型）;
對於vba，矩陣是一個LPSAFEARRAY，它具有按列主序排列的數據。

對於一維的數據（矢量）有追溯到BLAS的溶液（30年前），其可以在具有像

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

換句話說的結構被翻譯成 C，我們使用計算的中間結構不擁有數據，但可以映射兩個連續的數據或以恆定間隔（跨度）線性間隔的數據。結構的數據可以按順序進行處理，但它們可以也翻譯成陣列部分（如果使用的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開始，但即使這可能是抽象使用硬編碼1的位置的（可修改的）0-1參數。all（）函數/ FULL（）宏僅對整個連續矩陣正確，而不是對子矩陣爲正確。但是這也可以調整，在這種情況下切換到行上的循環。或多或少像上面的copyFrom方法（），它可以在行 - 主和列 - 主表示之間轉換（複製）矩陣。請注意TranspInPlace方法：通過交換rows/cols和rowstride/colstride，我們對相同的，未觸及的數據進行邏輯轉置，這意味着不再需要將邏輯交換機傳遞到通用例程（例如， GEMM爲矩陣乘法，即使（公平地說）這不包括共軛轉置的情況）。

鑑於上述情況，我在尋找實施這一辦法，這樣我可以使用（掛機）是一個庫，但我的搜索至今不盡如人意：

GSL GSL使用行優先的順序。停止。

LAPACK 本地代碼是列主要的。 C接口提供了處理行主數據的可能性，但只能添加特製的代碼或（在某些例行程序中）在對其進行操作之前對矩陣進行物理移位。

EIGEN 作爲一個模板庫，它可以同時處理這兩個問題。不幸的是，矩陣順序是模板的參數，其中意味着每個矩陣方法都將被重複。它有效，但並不理想。

請讓我知道，如果圖書館員更接近我後。謝謝。

來源

2016-11-05 maurizio nardò

行和列主要通過交換A次B和B次A並解釋輸出不同，不是？ – Yakk

在這兩種情況下，相同的矩陣以不同的方式映射到**線性**存儲器中，例如， [1 2 3; 4 5 6; 7 8 9]在內存中是：row-major：[1 2 3 4 5 6 7 8 9]，col-major：[1 4 7 2 5 8 3 6 9] –

在Eigen中，您可以使用Map模擬運行時內外跨度。舉例來說，如果你堅持ColumnMajor那麼內步幅對應的行間距和外步幅對應於列步長：

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

然後，你可以做任何操作上mat，比如訪問行mat.row(i)，列mat.col(j)，執行產品，解決線性系統等。

該方法的主要缺點是您鬆散SIMD矢量化。

來源

2016-11-05 13:23:54 ggael

我會看看。對於一個主要的（如FP12），我猜這個步伐是相反的，對吧？在編碼具有上述結構的稠密三角1-3時，我沒有大問題：一個混合的，可能可行的解決方案可能是使用上述結構，並且具有諸如線性系統分解之類的東西的蹦牀圖。您如何看待這個理念？ –

即使矩陣產品更好地使用Eigen，它將充分利用您的CPU（寄存器/緩存重用，矢量化等）。 – ggael

任何矩陣庫order-neutral？

回答

相關問題