?GEMM_

List

sgemm_cblas_sgemm単精度実　数一般行列と一般行列の積
dgemm_cblas_dgemm倍精度実　数一般行列と一般行列の積
cgemm_cblas_cgemm単精度複素数一般行列と一般行列の積
zgemm_cblas_zgemm倍精度複素数一般行列と一般行列の積

概略

一般行列と一般行列の積を計算します。結果を別途渡した行列にスカラ倍したものを加算します（詳しくは計算式参照）

計算式

C := alpha * AB + beta * C

beta=0のときは、Cへは単なる代入を行います。

サンプルコード

通常の実数の例 
#include <cblas.h>
#include <stdio.h>

int main(void)
{
  double A[2*3]={1.0,2.0,3.0, 
                 4.0,5.0,6.0};
  double B[3*4]={1.0,2.0,3.0,4.0,
                 5.0,6.0,7.0,8.0,
                 9.0,0.0,1.0,2.0};
  double C[2*4]; //x * y = 38 14 20 26
                 //        83 38 53 68 となることを確かめる
  
  //二次元配列の場合
  //double A[2][3]={{1.0,2.0,3.0},
  //                {4.0,5.0,6.0}}; //のように二次元配列で定義しても以下は同じコードで動く
  //double B[3][4]={{1.0,2.0,3.0,4.0},
  //                {5.0,6.0,7.0,8.0},
  //                {9.0,0.0,1.0,2.0}};
  //double C[2][4];
  
  int i,j;
  cblas_dgemm(CblasRowMajor,//通常のCの行列であればCBlasRowMajorを指定
              CblasNoTrans, //Aについて転置しない場合　CblasNoTrans 転置する場合　 CblasTrans
              CblasNoTrans, //Bについて転置しない場合　CblasNoTrans 転置する場合　 CblasTrans
              2,            //Aの行数
              4,            //Bの列数
              3,            //Aの列数　＝　Bの行数（一致していないと掛け算できない）
              1,            //alpha の値
              A,            //A
              3,            //leading dimension (通常はAの列数)
              B,            //B
              4,            //leading dimension (通常はBの列数)
              0,            //beta の値
              C,            //C
              4             //leading dimension (通常はCの列数)
              );

  for(i=0;i<2;i++){
    for(j=0;j<4;j++){
                                 //二次元配列の場合でもこの書き方でOK
        printf("%lf ",C[i*4+j]); //(i,j)成分は i*ldC+jの形で書ける ldCはCのleading dimension :Cの列数
        //printf("%lf ",C[i][j]);//もちろん、二次元配列で定義すればこれでも動く
    }
    printf("\n");
  }

  return 0;
}
転置する場合の例 
#include <cblas.h>
#include <stdio.h>

int main(void)
{
  double A[3*2]={1.0,2.0,
                 3.0,4.0,
                 5.0,6.0};
  double B[4*3]={1.0,2.0,3.0,
                 4.0,5.0,6.0,
                 7.0,8.0,9.0,
                 0.0,1.0,2.0};
  double C[2*4]; //x * y = 22 49  76 13
                 //        28 64 100 16 となることを確かめる
  int i,j;
  cblas_dgemm(CblasRowMajor,//通常のCの行列であればCBlasRowMajorを指定
              CblasTrans,   //Aについて転置しない場合　CblasNoTrans 転置する場合　 CblasTrans
              CblasTrans,   //Bについて転置しない場合　CblasNoTrans 転置する場合　 CblasTrans
              2,            //Aの行数(転置後)
              4,            //Bの列数(転置後)
              3,            //Aの列数(転置後)　＝　Bの行数(転置後)（一致していないと掛け算できない）
              1,            //alpha の値
              A,            //A
              2,            //leading dimension (通常はAの列数(転置前))
              B,            //B
              3,            //leading dimension (通常はBの列数(転置前))
              0,            //beta の値
              C,            //C
              4             //leading dimension (通常はCの列数)
              );

  for(i=0;i<2;i++){
    for(j=0;j<4;j++){
        printf("%lf ",C[i*4+j]); //(i,j)成分は i*ldC+jの形で書ける ldCはCのleading dimension :Cの列数
    }
    printf("\n");
  }

  return 0;
}

複素数の例
#include <cblas.h>
#include <stdio.h>

typedef struct
{
  double r;
  double i;
} complex_;

int main(void)
{
  //複素数を構造体無しで扱う場合はaxpyを参照
  complex_ A[2][2]={{{1.0,2.0},{3.0,4.0}},  // 1+2i 3+4i
                    {{5.0,6.0},{7.0,8.0}}}; // 5+6i 7+8i
  complex_ B[2][2]={{{8.0,7.0},{6.0,5.0}},  // 8+7i 6+5i
                    {{4.0,3.0},{2.0,1.0}}}; // 4+3i 2+1i
  complex_ C[2][2]; //x * y =-6 +  48 i, -2 + 28 i
                    //        2 + 136 i,  6 + 84 i となることを確かめる
  complex_ alpha={1,0};
  complex_ beta={0,0};
  int i,j;
  cblas_zgemm(CblasRowMajor,//通常のCの行列であればCBlasRowMajorを指定
              CblasNoTrans, //A 転置しない場合 CblasNoTrans 転置 CblasTrans 共役転置 CblasConjTrans
              CblasNoTrans, //B 転置しない場合 CblasNoTrans 転置 CblasTrans 共役転置 CblasConjTrans
              2,            //Aの行数
              2,            //Bの列数
              2,            //Aの列数　＝　Bの行数（一致していないと掛け算できない）
              &alpha,       //alpha の値：構造体へのポインタになる
              A,            //A
              2,            //leading dimension (通常はAの列数)
              B,            //B
              2,            //leading dimension (通常はBの列数)
              &beta,        //beta の値：構造体へのポインタになる
              C,            //C
              2             //leading dimension (通常はCの列数)
              );

  for(i=0;i<2;i++){
    for(j=0;j<2;j++){
        printf("%lf + %lf i, ",C[i][j].r,C[i][j].i);
    }
    printf("\n");
  }
  return 0;
}

引数／戻り値

cblas_dgemm

変数名型 上書き 概要
Order enum CBLAS_ORDER C言語では通常はCBlasRowMajor(列方向に並べる)
Fortran式ならCBlasColMajor(行方向に並べる)
Transaenum CBLAS_TRANSPOSE 行列Aの転置を指定
転置しない CblasNoTrans
転置 CblasTrans
共役転置 CblasConjTrans
Transbenum CBLAS_TRANSPOSE 行列Bの転置を指定
転置しない CblasNoTrans
転置 CblasTrans
共役転置 CblasConjTrans
m int 行列Aの行数
n int 行列Bの列数
k int 行列Aの列数、行列Bの行数
alpha double スカラーalpha
A double* 行列Aの先頭ポインタ
ldA int Aのleading dimension
CBlasRowMajorを指定したら列数
CBlasColMajor指定なら行数
B double* 行列Bの先頭ポインタ
ldB int Bのleading dimension
CBlasRowMajorを指定したら列数
CBlasColMajor指定なら行数
beta double スカラーbeta
C double* 上書き 行列Cの先頭ポインタ
ldC int Cのleading dimension
CBlasRowMajorを指定したら列数
CBlasColMajor指定なら行数
戻り値void

dgemm_

変数名型(dgemm_)概要
transachar* 行列Aの転置を指定 ("N"(そのまま）,"T"（転置）,"C"（共役転置）から選択)
transbchar* 行列Bの転置を指定 ("N"(そのまま）,"T"（転置）,"C"（共役転置）から選択)
m int* 行列Aの行数
n int* 行列Bの列数
k int* 行列Aの列数、行列Bの行数
alpha double*スカラーalpha
A double*行列Aの先頭ポインタ
ldA int* Aのleading dimension (通常は行数を指定すれば良い）
B double*行列Bの先頭ポインタ
ldB int* Bのleading dimension (通常は行数を指定すれば良い）
beta double*スカラーbeta
C double*行列Cの先頭ポインタ
ldC int* Cのleading dimension (通常は行数を指定すれば良い）
戻り値void

プロトタイプ宣言

void sgemm_(char *transa, char *transb, int *m, int *n, int *k,
float *alpha, float *A, int *ldA, float *B, int *ldB,
float *beta , float *C, int *ldC);

void dgemm_(char *transa, char *transb, int *m, int *n, int *k,
double *alpha, double *A, int *ldA, double *B, int *ldB,
double *beta , double *C, int *ldC);

void cgemm_(char *transa, char *transb, int *m, int *n, int *k,
complex *alpha, complex *A, int *ldA, complex *B, int *ldB,
complex *beta , complex *C, int *ldC);

void zgemm_(char *transa, char *transb, int *m, int *n, int *k,
doublecomplex *alpha, doublecomplex *A, int *ldA, doublecomplex *B, int *ldB,
doublecomplex *beta , doublecomplex *C, int *ldC);

void cblas_sgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB,
const blasint M, const blasint N, const blasint K,
const float alpha, const float *A, const blasint lda, const float *B, const blasint ldb, const float beta, float *C, const blasint ldc);

void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB,
const blasint M, const blasint N, const blasint K,
const double alpha, const double *A, const blasint lda, const double *B, const blasint ldb, const double beta, double *C, const blasint ldc);

void cblas_cgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB,
const blasint M, const blasint N, const blasint K,
const void *alpha, const void *A, const blasint lda, const void *B, const blasint ldb, const void *beta, void *C, const blasint ldc);

void cblas_zgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB,
const blasint M, const blasint N, const blasint K,
const void *alpha, const void *A, const blasint lda, const void *B, const blasint ldb, const void *beta, void *C, const blasint ldc);


※OPENBLAS_CONSTはcblas.hにて
# define OPENBLAS_CONST const
で定義されているため、ここではconstで置き換えている

※blasintはcommon.hにて
#if defined(OS_WINDOWS) && defined(__64BIT__)
typedef long long BLASLONG;
typedef unsigned long long BLASULONG;
#else
typedef long BLASLONG;
typedef unsigned long BLASULONG;
#endif

#ifdef USE64BITINT
typedef BLASLONG blasint;
#else
typedef int blasint;
#endif
で定義される(不要なコードは除いている)ため、環境により異なる。
初学者は通常はintの別名として定義されると考えて使用しても概ね差し支えない
(Cにおいてはintの値を渡しても自動的に昇格されるので)