おことわり
いかにズラズラとサンプルとその解説を載せていきます。この実験環境はCygwin上やLinux上でGotoBLASを使ったものです。基本的に環境依存では書いてないので、どの環境でも動くと信じていますが、仮に動作しなかったとしても当方では一切責任を負いませんのであしからず。
ちなみに、当方環境で動いたもののみ掲載しています。
ヘッダファイル
このサイトで紹介しているBLASルーチンの全プロトタイプ宣言を記載したファイルを置いておきます。以下では、このファイルを使用しています。
blas.h
基本的な使い方 - dgemvを例に
ソースコード
一例としてこのように用います。
#include<stdio.h>
#include<stdlib.h>
#include "blas.h"
main()
{
int i,j;
double *A;
double *x,*y;
int m=3,n=3;int ONE=1;
double alpha=1,beta=2;
A=(double*)malloc(3 * 3 * sizeof(double));
x=(double*)malloc( 3 * sizeof(double));
y=(double*)malloc( 3 * sizeof(double));
A[0]=1; A[3]=2; A[6]=3;
A[1]=0; A[4]=2; A[7]=3;
A[2]=0; A[5]=0; A[8]=3;
x[0]=1; y[0]=1;
x[1]=1; y[1]=1;
x[2]=1; y[2]=1;
dgemv_("N",&m,&n,&alpha,A,&m,x,&ONE,&beta,y,&ONE);
for(j=0;j<n;j++){
printf("%f ",y[j]);
}
printf("\n");
}
コンパイル
これをtestdgemv.cとすると、コンパイルは(/usr/lib にlibgoto2.aがあるとすると）
gcc testdgemv.c -lgoto2
とすると、コンパイルが通ります。
もし、
undefined reference to '_dgemv_'
ld returned 1 exit status
などのエラーが出てきた場合は、コンパイル自体には成功しているが、リンクに失敗した状態です。
libgoto2.aなどのファイルへのパスが通っていない状態となります
この時は、-Lオプションでパスを指定する必要があり
gcc testdgemv.c -L/home/Azalea/GotoBLAS2 -lgoto2
のように、パスを指定します。
これでもうまくコンパイルできない場合はBLASの生成時につけられた名前が異なる可能性があります当サイトでは、-DAdd_を使った場合を仮定しています。
Makefile等に -DNoChange, -DAdd__, -DUpCase,が指定されている場合は、-DAdd_に書き換えてコンパイルし直してください
メモリの格納
メモリへの格納は、列方向に行われます。詳しくは、
使い方概略を参考にしてください。
関数の呼び出し方と引数の渡し方
今回は、
dgemv_(一般行列とベクトルの積）を計算させています。
結果は、
8.000000 7.000000 5.000000
のようになるはずです。
さて、呼び出し方ですが、通常のCの関数のように呼び出せます。プロトタイプ宣言を書いているので、引数のチェックも行なってくれます。（余談ですが、
Cではコレがなくてもコンパイル通りますし、リンクも出来ます。試しに、blas.hをインクルードしないで実験してみてください(警告は出るかもしれませんが）。しかしこの場合、
引数のチェックは行わず、関数はすべてintを返すよう解釈してコードが読み込まれることに注意してください。C++ではこれができなくなりました）
引数は全てポインタになっている点に注意してください
これはFORTRANの仕様に依るものです。
そのため、定数を渡す際には、一度変数に詰め込んでから、渡す必要があります。
int m=3,n=3;int ONE=1;
double alpha=1,beta=2;
のようにしているのは、
dgemv_("N",&m,&n,&alpha,A,&m,x,&ONE,&beta,y,&ONE);
で、格納されているアドレスを渡す為です。
ただし、
char*だけは特別で、これだけは、変数に別途詰め込まず、直接""でくくって書くことができます。
Cでは文字列リテラルは全てポインタ(char*)で、これは読み込み専用の空間に格納されるという仕様になっている為です
"N"は転置を行わないことを指示する命令ですので、このままの形でAは演算に供されます
ldAとincx
ldAはLeading Dimensionの略称です。
配列Aの列数について聞いています。これは、
行列Aの列数とは無関係のものです。
FORTRANのルーチンでは、メモリアクセスの際、i行j列（i,jは0からm-1,n-1までとする）にアクセスする際に
A[i+j*ldA]という形でアクセスを行います。この時のldAを引数で欲しがっているのです
これは、
通常、Aの行数に相当します。そのため、このサイトでは、「通常は行数を指定すれば良い」と記しています。
このような引数がある理由は、例えば、Aの行列の一部のみを行列とみなして演算させるような状況があるから…だと思います。
そのため、行列の一部のみを演算に使う、など特殊な状況ではない限りは行数を格納した変数へのポインタを渡せば問題ありません
incxはベクトルxのインクリメント幅です。
これは、
xのi番目の成分が *(x+(i*incx)) でアクセスされるという意味になります。
例えば、xの要素の奇数番目だけ二倍したい、などという場合にはdscal_の引数incxに2を入れれば良いということになります。
もちろんこのような例は一般的ではないので、このサイトでは「通常1を指定すれば良い」としています
なお、
GotoBLASでは、incx=1の場合以外は最適化されていませんので、動作速度が低下する可能性があります