BLAS/CBLASの簡単な使い方

?AXPY_

List

saxpy_	cblas_saxpy	単精度実数ベクトル同士の加算
daxpy_	cblas_daxpy	倍精度実数ベクトル同士の加算
caxpy_	cblas_caxpy	単精度複素数ベクトル同士の加算
zaxpy_	cblas_zaxpy	倍精度複素数ベクトル同士の加算

概略

ベクトル同士の加算を行います。行列も大きさが非常に長いベクトルだと思えば使えます。与えたベクトルYの内容は破壊され、計算結果が書きこまれます。

計算式

Y := alpha * X + Y

サンプルコード

実数の例

#include <cblas.h>
#include <stdio.h>

int main(void)
{
  int i=0;
  double x[2] = {3.0,4.0};
  double y[2] = {7.0,8.0};
  double a = 2.0;
  cblas_daxpy(2, a, x, 1, y, 1); // 2*(3,4)+(7,8)
  for(i=0;i<2;i++){
    printf("%lf ",y[i]);         // (13,16)
  }
  printf("\n");
  
  return 0;
}

FORTRAN版へリンクする場合

#include <cblas.h>
#include <stdio.h>

int main(void)
{
  int i=0;
  double x[2] = {3.0,4.0};
  double y[2] = {7.0,8.0};
  double a = 2.0;
  int ONE = 1;
  int len = 2;
  
  daxpy_(&len,&a,x,&ONE,y,&ONE);
  
  for(i=0;i<2;i++){
    printf("%lf ",y[i]);
  }
  printf("\n");
  
  return 0;
}

複素数の例

#include <cblas.h>
#include <stdio.h>

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

int main(void)
{
  int i=0;
  
  //複素数を構造体無しで扱う場合
  double x[4] = {3.0,4.0,5.0,6.0};// x = (3+4i , 5+6i)
  double y[4] = {7.0,8.0,9.0,10.0}// y = (7+8i , 9+10i);
  double a[2] = {1.0,2.0};        // a =  1+2i
  cblas_zaxpy(2, a, x, 1, y, 1);  //n=2で指定するので注意
  for(i=0;i<4;i++)
    printf("%lf ",y[i]);          //(2+18i , 2+26i)
  printf("\n");
  
  //複素数を構造体で扱う場合
  complex_ xx[2]={{3.0,4.0},{5.0,6.0}};  // x = (3+4i , 5+6i)
  complex_ yy[2]={{7.0,8.0},{9.0,10.0}}; // y = (7+8i , 9+10i)
  complex_ aa={1.0,2.0};                 // a =  1+2i
  cblas_zaxpy(2,&aa,xx,1,yy,1);  // ( (1+2i)*(3+4i) + 7+8i ,(1+2i)*(5+6i)+(9+10i) )

  for(i=0; i<2; i++){
    printf("%lf + %lfi ", yy[i].r,yy[i].i); //(2+18i , 2+26i)
  }
  printf("\n");
  
  return 0;
}

引数／戻り値

cblas_daxpy

変数名	型	上書き	概要
n	int		ベクトルX,Yの大きさ（長さ）
alpha	double		Xベクトルのスカラ倍の係数
X	double*		ベクトルXの先頭ポインタ
incX	int		Xのインクリメント幅（通常1を指定すれば良い）
Y	double*	上書き	ベクトルYの先頭ポインタ
incY	int		Yのインクリメント幅（通常1を指定すれば良い）
戻り値	void

daxpy_

変数名	型	概要
n	int*	ベクトルX,Yの大きさ（長さ）
alpha	double*	Xベクトルのスカラ倍の係数
X	double*	ベクトルXの先頭ポインタ
incX	int*	Xのインクリメント幅（通常1を指定すれば良い）
Y	double*	ベクトルYの先頭ポインタ
incY	int*	Yのインクリメント幅（通常1を指定すれば良い）
戻り値	void

プロトタイプ宣言

void saxpy_(int *n, float *alpha, float *x, int *incx, float *y, int *incy);
void daxpy_(int *n, double *alpha, double *x, int *incx, double *y, int *incy);
void caxpy_(int *n, complex *alpha, complex *x, int *incx, complex *y, int *incy);
void zaxpy_(int *n, doublecomplex *alpha, doublecomplex *x, int *incx, doublecomplex *y, int *incy);

void cblas_saxpy(const blasint n, const float alpha, const float *x, const blasint incx, float *y, const blasint incy);
void cblas_daxpy(const blasint n, const double alpha, const double *x, const blasint incx, double *y, const blasint incy);
void cblas_caxpy(const blasint n, const void *alpha, const void *x, const blasint incx, void *y, const blasint incy);
void cblas_zaxpy(const blasint n, const void *alpha, const void *x, const blasint incx, void *y, const blasint incy);

※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の値を渡しても自動的に昇格されるので)

BLAS/CBLASの簡単な使い方 AXPY

menu