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

Question

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

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

• для функций листа собственный интерфейс Excel-C передает C матрицу в порядке следования строк (тип FP12 в 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 ]

Когда мы переходим от векторов к матрицам, я читал, что можно абстрагироваться от порядка матриц, используя как шаг строки, так и шаг столбца. Моя наивная попытка может быть:

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: меняя строки / столбцы и строки-строки / столбцы-строки, мы получаем логическое преобразование одних и тех же нетронутых данных, что означает, что больше не нужно передавать логические переключатели в универсальную подпрограмму (например, GEMM для умножения матриц)., хотя (чтобы быть справедливым) это не покрывает случай сопряженной транспозиции).

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

GSL Gsl использует порядок строк-мажоров. Стоп.

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

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

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

Answer 1

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

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

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

Основным недостатком этого подхода является потеря SIMD-векторизации.