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

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

Я работаю над библиотекой xll для Excel, то есть библиотекой C, содержащей функции, которые можно регистрировать и вызывать непосредственно из ячейки. В идеале эти функции должны вызываться (или адаптироваться для вызова) также из VBA,
чтобы обеспечить интерпретированную среду для более сложных вычислений (поиск корней, оды, оптимизация), которые не вписываются в лист 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]; }
}

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

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

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

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

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

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

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