Я хотел бы задать вам вопрос новичка 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;
}
Посмотрите на 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
Других решений пока нет …