Этот вопрос повторяется и на него следует дать более четкие ответы, чем "Matlab использует высоко оптимизированные библиотеки" или "Matlab использует MKL" один раз в Stackoverflow.
История
Матричное умножение (вместе с матрично-векторным, векторно-векторным умножением и многими матричными разложениями) является (являются) наиболее важной проблемой в линейной алгебре. Инженеры решали эти проблемы с компьютерами с первых дней.
Я не специалист по истории, но, видимо, тогда все просто переписали его версию на Фортране с помощью простых циклов. Затем пришла некоторая стандартизация с идентификацией «ядер» (базовых процедур), которые необходимы для решения большинства задач линейной алгебры. Эти основные операции были затем стандартизированы в спецификации под названием «Основные подпрограммы линейной алгебры (BLAS)». Инженеры могли бы тогда назвать эти стандартные, хорошо проверенные подпрограммы BLAS в своем коде, делая их работу намного проще.
BLAS:
BLAS эволюционировал от уровня 1 (первая версия, которая определяла скалярно-векторные и векторно-векторные операции) до уровня 2 (операции вектор-матрица) до уровня 3 (операции матрица-матрица) и предоставлял все больше и больше «ядер» так стандартизировано все больше и больше фундаментальных операций линейной алгебры. Оригинальные реализации Fortran 77 все еще доступны на сайте Netlib .
На пути к повышению производительности:
Таким образом, за эти годы (особенно между выпусками BLAS уровня 1 и уровня 2: в начале 80-х) аппаратное обеспечение изменилось с появлением векторных операций и иерархий кэша. Эти изменения позволили существенно повысить производительность подпрограмм BLAS. Затем разные поставщики внедрили свои подпрограммы BLAS, которые становились все более и более эффективными.
Я не знаю всех исторических реализаций (тогда я не был ни ребенком, ни ребенком), но две самые заметные из них появились в начале 2000-х годов: Intel MKL и GotoBLAS. Ваш Matlab использует Intel MKL, который является очень хорошим, оптимизированным BLAS, и это объясняет великолепную производительность, которую вы видите.
Технические подробности по умножению матриц:
Так почему же Matlab (MKL) так быстро работает на dgemm (общее умножение матрицы на матрицу с двойной точностью)? Проще говоря: потому что он использует векторизацию и хорошее кэширование данных. В более сложных терминах: см. Статью , предоставленную Джонатаном Муром.
По сути, когда вы выполняете умножение в предоставленном вами коде C ++, вы совсем не дружественны к кешу. Поскольку я подозреваю, что вы создали массив указателей на массивы строк, ваш доступ во внутреннем цикле к k-му столбцу "matice2": matice2[m][k] очень медленный. Действительно, когда вы обращаетесь к matice2[0][k], вы должны получить k-й элемент массива 0 вашей матрицы. Затем на следующей итерации вы должны получить доступ к matice2[1][k], который является k-м элементом другого массива (массив 1). Затем на следующей итерации вы получаете доступ к еще одному массиву и т. Д. Так как вся матрица matice2 не может поместиться в старшие кеши (размер 8*1024*1024 байт), программа должна извлечь нужный элемент из main память, теряя много времени.
Если вы просто транспонировали матрицу так, чтобы доступ осуществлялся по смежным адресам памяти, ваш код уже работал бы намного быстрее, потому что теперь компилятор может одновременно загружать целые строки в кэш. Просто попробуйте эту модифицированную версию:
timer.start();
float temp = 0;
//transpose matice2
for (int p = 0; p < rozmer; p++)
{
for (int q = 0; q < rozmer; q++)
{
tempmat[p][q] = matice2[q][p];
}
}
for(int j = 0; j < rozmer; j++)
{
for (int k = 0; k < rozmer; k++)
{
temp = 0;
for (int m = 0; m < rozmer; m++)
{
temp = temp + matice1[j][m] * tempmat[k][m];
}
matice3[j][k] = temp;
}
}
timer.stop();
Таким образом, вы можете видеть, как простота размещения кэша значительно увеличила производительность вашего кода.В настоящее время в реальных реализациях dgemm это используется на очень обширном уровне: они выполняют умножение на блоки матрицы, определяемые размером TLB (буфер трансляции с кратким изложением, если коротко: что может эффективно кэшироваться), чтобы они передавали потокдля процессора точно количество данных, которые он может обработать.Другим аспектом является векторизация, они используют векторизованные инструкции процессора для оптимальной пропускной способности команд, чего вы не можете сделать из своего кроссплатформенного кода C ++.
Наконец, люди утверждают, что это из-за Штрассена или Копперсмита -Алгоритм Винограда неверен, оба эти алгоритма не реализуются на практике из-за аппаратных соображений, упомянутых выше.