Базовый алгоритм линейной алгебры (MATLAB и C)

Фон: я конвертирую определенные функции MATLAB в C для конкретного приложения, чтобы уменьшить общие накладные расходы на память и повысить производительность для определенного алгоритма с помощью MEX. Некоторыми из различных включенных функций являются, например, алгоритмы декомпозиции/решателя линейных систем и их реализация в иерархии, аналогичной иерархии оператора деления слева в MATLAB.

Я понимаю, что MATLAB использует компилятор JIT, и многие люди сообщают, что трудно превзойти производительность MATLAB и т. д. и т. д., однако цель здесь состоит в том, чтобы на самом деле проверить эту теорию и действительно выяснить, можно ли добиться какого-либо увеличения производительности или эффективности.

В настоящее время я использую библиотеку GNU GSL.

К проблеме...

Я не получаю правильных результатов для базовой задачи обращения диагональной матрицы.

В MATLAB код:

x = A\B

где A — разреженная диагональная матрица, в которой каждый элемент на диагонали отличен от нуля, а B — просто обычная заполненная матрица, отвечающая требованиям размера для умножения матриц на A.

Чтобы решить эту проблему, не должен ли это быть просто базовый цикл for, который заменяет каждый элемент на главной диагонали обратным этому элементу или 1/элемент?

А потом умножить на матрицу B, конечно. Это кажется простой задачей, но она не дает таких же результатов, как MATLAB, хотя входные матрицы одинаковы.

Код С:

    gsl_matrix* matrix1;  
    gsl_matrix* matrix1_copy;
    int index; 

    ...  //  A gets initialized and filled as a sparse diagonal matrix

    gsl_matrix_memcpy(matrix1_copy, matrix1);
    gsl_vector diag = gsl_matrix_diagonal(matrix1_copy).vector;
    for (index = 0; index < diag.size; index++) {
        gsl_matrix_set(matrix1_copy, index, index, 1/gsl_vector_get(&diag, index));
    }

Это подтверждается заменой всех диагональных элементов их обратными. Эта матрица возвращается и затем используется как первая матрица в этой функции:

gsl_matrix* multiply(gsl_matrix* matrix1, gsl_matrix* matrix2) {
gsl_matrix* final_matrix = gsl_matrix_calloc(matrix1->size1, matrix2->size2);
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, matrix1, matrix2, 0.0, final_matrix);
return final_matrix;
}

Возвращаемая матрица записывается в файл, поэтому я могу легко ее прочитать. Весь этот процесс не выдает ни ошибок, ни предупреждений, и кажется, что каждый шаг выполнен, но окончательная матрица не соответствует своей копии в MATLAB; сначала трудно сказать, так как это очень большая разреженная матрица, но начало обеих сравниваемых матриц правильное, однако, вложив нижний правый угол матрицы в MATLAB, я наблюдаю числа, которые просто 0 в моем C код. Моя точность настроена на правильный вывод.

Любые идеи?


c matlab matrix gsl
person Community    schedule 05.08.2014    source источник
comment
Пара вещей, чтобы проверить. 1) Использовать тот же алгоритм в Matlab, т.е. сделать x=AA*B, где AA=1./A. Я почти уверен, что x=A/B не выполняет inv(A)*B, а решает систему A*x=B без вычисления обратной матрицы. Поскольку выполняемые операции не идентичны из-за ошибок округления, вы можете получить несколько разные значения. 2) Чтобы исключить неправильное использование библиотеки GSL, кодируйте умножение самостоятельно, используя циклы (т.е. x(i,k)=sum(1/A(i,j)*B(j,k), где сумма превышает j) 3) Возможно значения становятся нулями при записи в файл. Проверьте, как вы записываете/читаете из файла.   -  person PetrH    schedule 05.08.2014

Ответы (1)


arrow_upward
0
arrow_downward

Вы должны использовать обратную подстановку и рассматривать диагональную матрицу A как особый случай. Таким образом, вы охватите более широкий спектр проблем с помощью одного и того же алгоритма (и тех же модульных тестов!)

Первый дополнительный комментарий: поскольку вы разрабатываете алгоритм с нуля, вы должны использовать простые тестовые примеры. Тестирование со сложными случаями будет затруднено. Я имею в виду, что ваши тестовые случаи должны быть настолько малы, чтобы вы могли проверить результат вручную.

Второй дополнительный комментарий: Matlab использует библиотеки BLAS/LAPACK, оптимизированные Intel, также известные как Библиотека ядра Intel Math. Если ваша цель — превзойти Matlab, вам следует использовать BLAS или BLAS-подобную библиотеку: перейти к BLAS, OpenBLAS или FLAME.

Последний дополнительный комментарий: вы не первый человек, которого я знаю, который хочет победить Matlab. Учитесь на ошибках других. Выберите целевую платформу и оптимизируйте как черт. Знай свой процессор до мозга костей. Используйте многопоточность. Используйте МПИ. Возьмите это как свою библию: методы оптимизации Agner

person gire    schedule 05.08.2014
comment
Спасибо за ответ, особенно за ссылку на руководство по оптимизации. - person ; 09.08.2014