- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello,
I modified an example from the intel website to use std::complex<double> type with the zgemm function from MKL.
Unfortunately I am getting wrong results for complex matrix multiplication. I am wondering what I am doing wrong.
Thank in advance
Dmitry
=========================================================================
#include <complex>
#include <iostream>
#include <iomanip>
#define MKL_Complex16 std::complex<double>
#include "mkl.h"
typedef std::complex<double> Complex;
const int Arows = 2, Acols = 3, Brows = 3, Bcols = 2, Crows = 2, Ccols = 2;
int main()
{
const char *notrans = "n";
const char *trans = "t";
const char *conj = "c";
Complex A[Arows][Acols], A1[Arows][Acols], B[Brows][Bcols], B1[Brows][Bcols], C[Crows][Ccols], C1[Crows][Ccols];
Complex alpha(1.0, 0.0);
Complex beta(0.0, 0.0);
for (int i = 0; i < Arows; i++){
for(int j = 0; j <Acols; j++){
A
A1
}
}
for (int i = 0; i < Brows; i++){
for(int j = 0; j < Bcols; j++){
B
B1
}
}
for (int i = 0; i < Crows; i++){
for(int j = 0; j < Ccols; j++){
C
C1
}
}
Complex* pma = &(A[0][0]);
Complex* pmb = &(B[0][0]);
Complex* pmc = &(C[0][0]);
std::cout << "Matrix A" << std::endl;
for (int i = 0; i < Arows; i++){
for(int j = 0; j < Acols; j++){
std::cout << std::setw(4) << A
}
std::cout << std::endl;
}
std::cout << "Matrix B" << std::endl;
for (int i = 0; i < Brows; i++){
for(int j = 0; j < Bcols; j++){
std::cout << std::setw(4) << B
}
std::cout << std::endl;
}
int iArows = Arows, iAcols = Acols, iBrows = Brows, iBcols = Bcols, iCrows = Crows, iCcols = Ccols;
cblas_zgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans, iArows, iBcols, iAcols, &alpha, pma, iAcols, pmb, iBrows, &beta, pmc, iArows);
for (int i = 0; i < Arows; i++){
for (int j = 0; j < Bcols; j++){
for (int k = 0; k < Acols; k++){
C1
};
};
};
std::cout << "From zgemm" << std::endl;
for (int i = 0; i < Crows; i++){
for(int j = 0; j < Ccols; j++){
std::cout <<std::setw(4) << C
}
std::cout << std::endl;
}
std::cout << "From direct" << std::endl;
for (int i = 0; i < Crows; i++){
for(int j = 0; j < Ccols; j++){
std::cout << std::setw(4) << C1
}
std::cout << std::endl;
}
return 0;
}
========================================================================
Output:
(1,1)(1,2)(1,3)
(2,1)(2,2)(2,3)
Matrix B
(-1,-1)(-2,-1)
(-1,-2)(-2,-2)
(-1,-3)(-2,-3)
From zgemm
(4,-12)(5,-15)
(0,-16)(0,-20)
From direct
(11,-12)(8,-18)
(8,-18)(2,-24)
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
You passed the wrong values of ldA, ldB and ldC in the argument list to cblas_zgemm. The correct invocation is:
cblas_zgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans, iArows, iBcols, iAcols, &alpha, pma, iAcols, pmb, iBcols, &beta, pmc, iCcols);
Since all three matrices are stored as C-native (i.e., row-major) arrays, the ldA, ldB and ldC arguments should be the column-counts of A, B and C.
You may find clarification on this issue by looking at the following lines from the file cblas_zgemmx.c in the MKL examples.
if( layout == CblasRowMajor ) { lda=cmaxa; ldb=cmaxb; ldc=cmaxc; } else { lda=rmaxa; ldb=rmaxb; ldc=rmaxc; }
I suspect that there are a number of errors in the C version of the documentation (https://software.intel.com/en-us/node/520775) in this regard. The MKL-C documentation was recently separated from the MKL-Fortran documentation and, given that Fortran stores 2-D arrays by columns and C does so by rows, such confusion is quite common. The Fortran routine zgemm (and its siblings) provide only for native (column-major) arrays, but the C routines provide for both native (row-major) and alien (column-major) arrays, so the C-documentation cannot avoid more complexity.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Add comments,
You may understand the lda, ldb, ldc as the distance of row when row-major , the value is column counts. or the distance of column when column-major, the value is row counts.
Best Regards,
Ying
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thank you very much!
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Sorry, but I still have trouble with zgemm. Now I want to calculate AT * B. My code is now
// A*B
// cblas_zgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans, iArows, iBcols, iAcols, &alpha, pma, iAcols, pmb, iBcols, &beta, pmc, iCcols);
//AT*B
cblas_zgemm( CblasRowMajor, CblasTrans, CblasNoTrans, iArows, iBcols, iAcols, &alpha, pma, iArows, pmb, iBcols, &beta, pmc, iCcols);
It looks like I have only to change the lda parameter compared to the previous case (A * B ).
But I have still different results:
Matrix A
(1,1)(1,2)
(2,1)(2,2)
(3,1)(3,2)
Matrix B
(-1,-1)(-2,-1)
(-1,-2)(-2,-2)
(-1,-3)(-2,-3)
From zgemm
(2,-8)(-1,-11)
(0,-10)(-4,-13)
From direct
(0,-28)(-11,-34)
(6,-20)(0,-26)
Could you please correct or explain this?
Best wishes
Dmitry
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi Dmitry,
I change the code as below and get the result
#include <complex>
#include <iostream>
#include <iomanip>
#define MKL_Complex16 std::complex<double>
#include "mkl.h"
typedef std::complex<double> Complex;
const int Arows = 3, Acols = 2, Brows = 3, Bcols = 2, Crows = 2, Ccols = 2;
int main()
{
const char *notrans = "n";
const char *trans = "t";
const char *conj = "c";
Complex A[Arows][Acols], A1[Arows][Acols], B[Brows][Bcols], B1[Brows][Bcols], C[Crows][Ccols], C1[Crows][Ccols];
Complex alpha(1.0, 0.0);
Complex beta(0.0, 0.0);
for (int i = 0; i < Arows; i++){
for(int j = 0; j <Acols; j++){
A
A1
}
}
for (int i = 0; i < Brows; i++){
for(int j = 0; j < Bcols; j++){
B
B1
}
}
for (int i = 0; i < Crows; i++){
for(int j = 0; j < Ccols; j++){
C
C1
}
}
Complex* pma = &(A[0][0]);
Complex* pmb = &(B[0][0]);
Complex* pmc = &(C[0][0]);
std::cout << "Matrix A" << std::endl;
for (int i = 0; i < Arows; i++){
for(int j = 0; j < Acols; j++){
std::cout << std::setw(4) << A
}
std::cout << std::endl;
}
std::cout << "Matrix B" << std::endl;
for (int i = 0; i < Brows; i++){
for(int j = 0; j < Bcols; j++){
std::cout << std::setw(4) << B
}
std::cout << std::endl;
}
int iArows = Arows, iAcols = Acols, iBrows = Brows, iBcols = Bcols, iCrows = Crows, iCcols = Ccols;
cblas_zgemm( CblasRowMajor, CblasTrans, CblasNoTrans, iAcols, iBcols, iArows, &alpha, pma, iAcols, pmb, iBcols, &beta, pmc, iCcols);
for (int i = 0; i < Acols; i++){
for (int j = 0; j < Bcols; j++){
for (int k = 0; k < Arows; k++){
C1
};
};
};
std::cout << "From zgemm" << std::endl;
for (int i = 0; i < Crows; i++){
for(int j = 0; j < Ccols; j++){
std::cout <<std::setw(4) << C
}
std::cout << std::endl;
}
std::cout << "From direct" << std::endl;
for (int i = 0; i < Crows; i++){
for(int j = 0; j < Ccols; j++){
std::cout << std::setw(4) << C1
}
std::cout << std::endl;
}
return 0;
}.
MKL C manual have some explanation and table about m, n, k, lda etc.
The key is as the m is Op(A)'s row, not A's row. So mxk * kxn = mxn in your case, (2x3)* (3x2 ) = 2 * 2 , m=2, n=2, k=3.
const int Arows = 3, Acols = 2, Brows = 3, Bcols = 2, Crows = 2, Ccols = 2;
cblas_zgemm( CblasRowMajor, CblasTrans, CblasNoTrans, iAcols, iBcols, iArows, &alpha, pma, iAcols, pmb, iBcols, &beta, pmc, iCcols);
m INTEGER. Specifies the number of rows of the matrix op(A) and of the
matrix C. The value of m must be at least zero.
n INTEGER. Specifies the number of columns of the matrix op(B) and the
number of columns of the matrix C.
The value of n must be at least zero.
k INTEGER. Specifies the number of columns of the matrix op(A) and the
number of rows of the matrix op(B).
The value of k must be at least zero.
lda Specifies the leading dimension of a as declared in the calling
(sub)program.
transa=CblasNoTrans transa=CblasTrans or
transa=CblasConjTrans
Layout =
CblasColMajor
lda must be at least
max(1, m).
lda must be at least
max(1, k)
Layout =
CblasRowMajor
lda must be at least
max(1, k)
lda must be at least
max(1, m).
Best Regards,
Ying
Matrix A
(1,1)(1,2)
(2,1)(2,2)
(3,1)(3,2)
Matrix B
(-1,-1)(-2,-1)
(-1,-2)(-2,-2)
(-1,-3)(-2,-3)
From zgemm
(0,-17)(-6,-20)
(6,-20)(0,-26)
From direct
(0,-17)(-6,-20)
(6,-20)(0,-26)
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thank you very much!
![](/skins/images/45C6C2D737ED71F4C51F0145C8CB1E9C/responsive_peak/images/icon_anonymous_message.png)
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page