mkl blas and std::complex

37 wyświetleń
Przejdź do pierwszej nieodczytanej wiadomości

Dmitry Ganyushin

nieprzeczytany,
31 gru 2015, 13:39:4431.12.2015
do matrixprogramming
Hello,

Does anybody know how to feed correctly std::complex matrices to zgemm function in mkl BLAS? I can compile the code below,  but the result is wrong:
====================================================================================
#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[i][j] = Complex(i+1,j+1);
           A1[i][j] = Complex(i+1,j+1);     
         }
       }
      
       for (int i = 0; i < Brows; i++){
         for(int j = 0; j < Bcols; j++){
           B[i][j] = -Complex(j+1,i+1);
           B1[i][j] = -Complex(j+1,i+1);
         }
       }
      
       for (int i = 0; i < Crows; i++){
         for(int j = 0; j < Ccols; j++){
           C[i][j] = Complex(0.0, 0.0);
           C1[i][j] = Complex(0.0, 0.0);
         }
       }
       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[i][j];
         }
         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[i][j];
         }
         std::cout << std::endl;
       }
      
       int iArows = Arows, iAcols = Acols, iBrows = Brows, iBcols = Bcols, iCrows = Crows, iCcols = Ccols;
      
       zgemm( notrans, notrans, &iArows, &iBcols, &iAcols, &alpha, pma, &iArows,  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[i][j] += A1[i][k] * B1[k][j];
          };
        };
      };
     

       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[i][j];
         }
         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[i][j];
         }
         std::cout << std::endl;
       }
      
      
       return 0;
}
=============================================================================================================

Result:

dmitry@linux-4doj:~/tests> ./zgemm-clean
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)(-4,-8)
(-5,-10)(4,-12)
From direct
(5,-22)(2,-26)
(1,-9)(-3,-12)


Any help is very appreciated.

Dmitry 
Odpowiedz wszystkim
Odpowiedz autorowi
Przekaż
Nowe wiadomości: 0