Любая матричная библиотека порядка-нейтральная?

Question

Извините, если это слишком долго, но я чувствую, что вопрос должен быть уточнен:

Я работаю над библиотекой xll для Excel, то есть библиотекой C, содержащей функции, которые можно зарегистрировать и вызывать непосредственно из клетка. В идеале эти функции должны быть вызваны (или адаптированы для вызова) также от VBA, , чтобы обеспечить интерпретируемую среду для более сложных вычислений (поиск корней, ode, оптимизация), которые не очень хорошо вписываются в лист Excel. Чтобы быть ясным: есть способ вызвать функцию листа из vba (функция Application.Run), но она платит недопустимую цену конверсий из/в тип варианта для всех параметров и возвращаемого значения. Сейчас забавная ситуация в том, что, в той же среде Excel, двумерный матрица передается по-разному:

• для листовых функций, нативный интерфейс Excel-C переходит в C матрицу в строке-мажорный порядке (тип FP12 Excel SDK);

• для vba, матрица представляет собой LPSAFEARRAY, которая имеет данные, упорядоченные по порядку столбцов.

Для одномерных данных (векторы) существует решение восходит к BLAS (30 лет назад), который может быть переведен на С в имеющий структуру, как

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

Другими словами, мы используем для расчетов промежуточную структуру, которая не владеет данными, но может отображать как смежные данные или данные, линейно расположенные с постоянным зазором (шаг). Данные структуры могут быть обработаны последовательно, но они могут также переведены на секцию массива (при использовании 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, но даже это может быть отведенным используя (модифицируемый) параметр 0-1 в место жестко закодированного 1. Макрос функции all()/FULL() верны только для целой, смежной матрицы, а не подматрицы. Но также это можно было бы скорректировать, переключив в этом случае на цикл над строками. Более или менее подобный выше метод copyFrom(), который может преобразовывать (копировать) матрицу между строкой-строкой и представлением-столбцом.

Обратите также внимание на метод TranspInPlace: путем замены строк/столбцов и rowstride/colstride у нас есть логическая перестановка одних и тех же нетронутых данных, что означает, что больше не нужно передавать логические коммутаторы в универсальную подпрограмму (например, GEMM для матричное умножение, хотя (если быть справедливым) это не распространяется на случай сопряженной транспозиции).

Учитывая вышеизложенное, я ищу для библиотеки, реализующей этот подход, так что я могу использовать (крюк), но мой поиск до сих пор не является удовлетворительным:

GSL Gsl использует строки-мажорных заказ. Стоп.

LAPACK Нативный код является основным. C-интерфейс дает возможность обрабатывать данные основных строк, но только добавляет индивидуальный код или (в некоторой рутинной) физическую перенос матрицы перед его работой.

EIGEN Будучи темплатной библиотекой, он может обрабатывать и то, и другое. К сожалению, порядок матрицы является параметром шаблона, который означает, что каждый метод матрицы будет дублироваться. Он работает, но не идеален.

Пожалуйста, дайте мне знать, если библиотекари ближе к тому, что я хочу. Спасибо.

Answer 1

В Eigen вы можете имитировать это с помощью Map с внутренними и внешними шагами выполнения. Например, если придерживаться ColumnMajor то внутрипартийной шаг соответствует пропашной походке, и внешний шаг соответствует колонной походке:

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

После этого вы можете делать какую-либо операцию по mat, как доступу к строки mat.row(i), столбцы mat.col(j), выполняет продукты, решает линейные системы и т. д.

Основным недостатком этого подхода является то, что вы теряете векторную схему SIMD.