Продукт BLAS dgemm неожиданно ведет себя с опцией CblasTrans

Я хотел бы задать вам вопрос новичка BLAS.
Казалось бы, простая задача о матричном умножении матрицы А на ее собственную
транспонировать: C: = A ‘* A

Мой пример (2×3): A: = [1 2 3; 4 5 6].
Следовательно, A ‘есть (3×2) и C должно быть (3×3).

В Row Major и планировании использовать опцию CblasTrans я бы ожидал
lda = ldb = 3 в обоих случаях A и A ‘.

К сожалению, нижняя демонстрационная программа все еще генерирует совершенно неправильный продукт
и простая перестановка параметров с моей стороны пока не достигла цели.
На самом деле полученные значения смехотворно высоки, и я
сбит с толку 6-элементной структурой результата.

Что мне здесь не хватает?

/**
* transposeMat.cpp, compile using: g++ -lcblas transposeMat.cpp
*/

#include <cstdlib>
#include <cblas.h>
#include <iostream>
#include <sstream>
#include <string>

using namespace std;

string matrix2string(int m, int n, double* A, CBLAS_ORDER order)
{
ostringstream oss;
for (int j=0;j<m;j++)
{
for (int k=0;k<n;k++)
{
switch (order)
{
case CblasRowMajor:
oss << A[j*n+k];
break;
case CblasColMajor:
oss << A[j+k*m];
break;
default:
return "[matrix2string(..): Unknown order.]";
}
if (k < n-1) oss << '\t';
}
if (j < m-1) oss << endl;
}
return oss.str();
}

int main(int argc, char** argv)
{
int m=2;
int n=3;
// RowMajor matrix [ 1,2,3 ; 4,5,6 ]
double A[6] = { 1,2,3,4,5,6 };
// Using A for both xgemm-Parameters brings no luck! This is not enough though.
double B[6] = { 1,2,3,4,5,6 };
// Container for the result which will be 3x3.
double C[9] = { 0,0,0,0,0,0,0,0,0 };
// C:=A'A
// Params: (Majority,TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
cblas_dgemm(CblasRowMajor,CblasTrans,CblasNoTrans,m,n,n,1,&A[0],n,&B[0],n,0,&C[0],n);
//> ADDED COMMENT AFTER aka.nice ANSWERED THE QUESTION. ----------
// 1.: "MxN" really are the dimensions of matrix C and K is the "in-between"//   dimension shared by the factors of the product.
// 2.: The op(A) on the BLAS reference card actually seems to read "after
//   the internal transpose of A".
// 3.: Taken this into the code the above matrix B also becomes unnecessary.
// Hence this programm runs expectedly if you
//   replace the upper line by:
// cblas_dgemm(CblasRowMajor,CblasTrans,CblasNoTrans,n,n,m,1,&A[0],n,&A[0],n,0,&C[0],n);
//< --------------------------------------------------------------
cout << "A:" << endl << matrix2string(m,n,&A[0],CblasRowMajor).c_str() << endl <<
"C:" << endl << matrix2string(n,n,&C[0],CblasRowMajor).c_str() << endl;
/** Output:
A:
1       2       3
4       5       6
C:
34      44      54
90      117     144
0       0       0
*/
return EXIT_SUCCESS;
}

3

Решение

Посмотрите на DGEMM из netlib: http://www.netlib.org/blas/dgemm.f

Вы увидите, что:

*  DGEMM  performs one of the matrix-matrix operations
*
*     C := alpha*op( A )*op( B ) + beta*C,

и это:

*  M      - INTEGER.
*           On entry,  M  specifies  the number  of rows  of the  matrix
*           op( A )  and of the  matrix  C.  M  must  be at least  zero.
*           Unchanged on exit.

Таким образом, если A есть (2,3), то op (A) = A ‘есть (3,2).

Если вы посмотрите на определение для других аргументов, вы увидите, что вы должны передать M = 3, N = 3, K = 2

2

Другие решения

Других решений пока нет …

По вопросам рекламы ammmcru@yandex.ru
Adblock
detector