Вызов dgemm содержит два набора целых чисел. Во-первых, определить размеры матриц, как описано в математике. Глядя на документацию для BLAS по адресу http://www.netlib.org/lapack/explore-html/d7/d2b/dgemm_8f.html, мы можем увидеть аргументы подпрограммы, определенной как (и отформатировать ее так, чтобы я помнил, как она работает):
Subroutine dgemm (TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, &
B, LDB, &
BETA, C, LDC)
Набор целых чисел для определения математических измерений: M, N, K В этом наборе
- M - первое математическое измерение матрицы результатов,
C - N - второе математическое измерение матрицы результатов
C - K - длина точечного произведения, необходимая для создания отдельного элемента матрицы результатов
C
Итак M, N, K чисто по математике задачи, больше ничего. Но мы находимся на компьютере, поэтому возникает второй вопрос - как каждая из двухмерных матриц размещается в памяти компьютера, которая по своей сути является одномерным объектом? Теперь Фортран хранит свои объекты в порядке Колонна Майор. Это означает, что в памяти, если у вас есть элемент A (i, j), следующая ячейка памяти будет содержать A (i + 1, j), если мы не являемся концом столбца, или A (1, j + 1), если мы. Таким образом, матрица построена так, что «первый индекс движется быстрее всего». Чтобы определить это расположение двумерного объекта, все, что нам нужно сказать программе, это то, сколько элементов A вам нужно перепрыгнуть на go с A (i, j) на A (i, j + 1) - и это это число, которое является «ведущим измерением», LDA, LDB, LDC в документации. И, таким образом, это просто число строк в матрице, как объявлено или выделено.
Теперь давайте посмотрим на цитируемую вами документацию
LDB is INTEGER
On entry, LDB specifies the first dimension of B as declared
in the calling (sub) program. When TRANSB = 'N' or 'n' then
LDB must be at least max( 1, k ), otherwise LDB must be at
least max( 1, n ).
LDB - это количество строк в матрице B это объявлено. Таким образом,
- Если B НЕ транспонирован, каждый столбец B будет включен в точечное произведение со строкой или столбцом A, в зависимости от опции транспонирования для A. Математические измерения говорят, что точечное произведение к долго. Таким образом, число строк в B должно быть не менее k, чтобы математика была правильной, и поэтому LDB должен быть не менее k
- Если B IS транспонирован, то математическое первое измерение B также будет математическим вторым размерность матрицы результатов, и это задается N. Таким образом, LDB в этом случае должен быть не менее N
, так что это отвечает на большинство из них, и если вы всегда умножаете всю одну матрицу на объявленный всей второй матрицей, ведущее измерение будет просто самым первым измерением каждой матрицы, как объявлено. Но «по крайней мере» означает, что вы можете умножать биты массивов вместе. Подумайте о следующем, и, аккуратно разделив те целые числа, которые определяют математику, и те, которые определяют макет памяти, вы можете понять, почему он делает то, что делает? Обратите внимание, что (2, 2) в списке аргументов означает, что мы «запускаем» матрицу в этом элементе - теперь тщательно продумайте, что говорит вам ведущее измерение, и вы сможете разобраться, как это работает.
ian@eris:~/work/stack$ cat matmul2.f90
Program sirt
Integer, Parameter :: wp = Kind( 1.0d0 )
Real( wp ), Dimension( 1:5, 1:5 ) :: A
Real( wp ), Dimension( 1:3, 1:3 ) :: B
Real( wp ), Dimension( 1:4, 1:4 ) :: C
Integer :: i
A = 0.0_wp
Do i = 1, 5
A( i, i ) = Real( i, wp )
End Do
Write( *, * ) 'A = '
Do i = 1, Size( A, Dim = 1 )
Write( *, '( 100( f4.1, 1x ) )' ) A( i, : )
End Do
B = 0.0_wp
B( 1, 1 ) = 1.0_wp
B( 3, 2 ) = 1.0_wp
B( 2, 3 ) = 1.0_wp
Write( *, * ) 'B = '
Do i = 1, Size( B, Dim = 1 )
Write( *, '( 100( f4.1, 1x ) )' ) B( i, : )
End Do
! Lazy - should really just initialise only the bits of C that are NOT touched
! by the dgemm
C = 0.0_wp
Call dgemm('N', 'N', 3, 3, 3, 1.0_wp, A( 2, 2 ), Size( A, Dim = 1 ), &
B , Size( B, Dim = 1 ), &
0.0_wp, C , Size( C, Dim = 1 ) )
Write( *, * ) 'C after dgemm'
Write( *, * ) 'B = '
Do i = 1, Size( C, Dim = 1 )
Write( *, '( 100( f4.1, 1x ) )' ) C( i, : )
End Do
End Program sirt
ian@eris:~/work/stack$ gfortran-8 -std=f2008 -fcheck=all -Wall -Wextra -pedantic -O matmul2.f90 -lblas
/usr/bin/ld: warning: libgfortran.so.4, needed by //usr/lib/x86_64-linux-gnu/libopenblas.so.0, may conflict with libgfortran.so.5
ian@eris:~/work/stack$ ./a.out
A =
1.0 0.0 0.0 0.0 0.0
0.0 2.0 0.0 0.0 0.0
0.0 0.0 3.0 0.0 0.0
0.0 0.0 0.0 4.0 0.0
0.0 0.0 0.0 0.0 5.0
B =
1.0 0.0 0.0
0.0 0.0 1.0
0.0 1.0 0.0
C after dgemm
B =
2.0 0.0 0.0 0.0
0.0 0.0 3.0 0.0
0.0 4.0 0.0 0.0
0.0 0.0 0.0 0.0
Это использование для умножения части двух больших матриц на удивление распространено, особенно в алгоритмах блокированной линейной алгебры, и разделение математических и компоновочных измерений позволяет вам сделать это