ちゃんと、_mm256_addやら_mm256_fmaddも使います。
が、その前に、２次元の複素数配列をもう一回attributeを使って行ってみます。
複素数は、
std::complex
を使うのが一般的でしょうが、これを使うのはあまりうまいやり方ではありません。レジスタ上に実数、虚数、実数、虚数と並ぶ事になるからです。
もちろん、AVX系命令を使えば、複素数×複素数のヴェクトル計算も出来ますが、今はあくまでattributeを使ってみると言う趣旨なので、実数は実数だけでベクトルレジスタ上に、虚数は虚数だけでヴェクトルレジスタ上に並ぶようにしてみましょう。
２次元配列なので、XY平面上の格子点に複素数が乗っているとイメージしてください。レジスタ上に乗る単精度数は８個なので、Y方向を４つ、X方向を２つに分割する事にします。
このとき、各配列と、SIMDのヴェクトルレジスタにのる配列を格納するAccess関数が変わっている事に注意してください。
#include < stdio.h>
#include < stdlib.h>
#include < math.h>
#include < sys/time.h>
#include < iostream>
#define L 512
#define Ncmplx 2
#define Re 0
#define Im 1
#ifdef SSE
#define NSIMD (4)
#define NBT (2)
#define NBX (2)
typedef float f4 __attribute__ ((vector_size (16))); // SSE vector, 16 bytes = 4 floats
typedef f4 SIMDfloat;
#endif
#ifdef AVX
#define NSIMD (8)
#define NBT (4)
#define NBX (2)
typedef float f8 __attribute__ ((vector_size (32))); // AVX vector, 32 bytes = 8 floats
typedef f8 SIMDfloat;
#endif
#define BT (L/NBT)
#define BX (L/NBX)
float U [L][L][Ncmplx];
float V [L][L][Ncmplx];
float W [L][L][Ncmplx];
SIMDfloat Usimd [L/NBT][L/NBX][Ncmplx];
SIMDfloat Vsimd [L/NBT][L/NBX][Ncmplx];
SIMDfloat Wsimd [L/NBT][L/NBX][Ncmplx];
#define Access(A,t,x,reim) ((float *)&A[t%BT][x%BX][reim])[ (x/BX)+NBX*(t/BT) ]
#define Access_A(t,x,reim) std::cout << t << " " << x << " : " << t%BT << " " << x%BX << " " << x/BX+NBX*(t/BT) << std::endl;
void report(char *routine, struct timeval *t1, struct timeval * t2)
{
struct timeval diff;
timersub(t2,t1,&diff);
printf("%s took %ld s %6ld usec\n",routine,(long)diff.tv_sec,(long)diff.tv_usec);
}
void scalar_multiply(void){
for(int t=0; t< L; t++) {
for(int x=0; x< L; x++) {
U[t][x][Re]= V[t][x][Re] * W[t][x][Re] - V[t][x][Im] * W[t][x][Im] ;
U[t][x][Im]= V[t][x][Re] * W[t][x][Im] + V[t][x][Im] * W[t][x][Re] ;
}
}
}
void simd_multiply(void){
for(int t=0; t< BT; t++) {
for(int x=0; x< BX; x++) {
Usimd[t][x][Re]= Vsimd[t][x][Re] * Wsimd[t][x][Re] - Vsimd[t][x][Im] * Wsimd[t][x][Im];
Usimd[t][x][Im]= Vsimd[t][x][Re] * Wsimd[t][x][Im] + Vsimd[t][x][Im] * Wsimd[t][x][Re] ;
}
}
}
int main(void){
for(int t=0; t< L; t++){
for(int x=0; x< L; x++){
Access(Vsimd,t,x,Re)=V[t][x][Re] = drand48();
Access(Vsimd,t,x,Im)=V[t][x][Im] = drand48();
Access(Wsimd,t,x,Re)=W[t][x][Re] = drand48();
Access(Wsimd,t,x,Im)=W[t][x][Im] = drand48();
Access_A(t,x,Im);
}
}
int i;
struct timeval tv1,tv2,tv3,tv4;
gettimeofday(&tv1,(struct timezone *)NULL);
for(i=0; i< 100; i++){
scalar_multiply();
}
gettimeofday(&tv2,(struct timezone *)NULL);
for(i=0; i< 100; i++){
simd_multiply();
}
gettimeofday(&tv3,(struct timezone *)NULL);
// Check results
for(int t=0; t< L; t++){
for(int x=0; x< L; x++){
if(Access(Usimd,t,x,Re)!=U[t][x][Re] )
printf(" Wrong %d %d %f %f\n",t,x,Access(Usimd,t,x,Re),U[t][x][Re]);
if(Access(Usimd,t,x,Im)!=U[t][x][Im] )
printf(" Wrong %d %d %f %f\n",t,x,Access(Usimd,t,x,Im),U[t][x][Im]);
}
}
report((char *)"scalar multiply ",&tv1,&tv2);
report((char *)"simd multiply ",&tv2,&tv3);
}
このプログラムでは、SIMDベクトル上に配置する関数Accessは、以下のように代わります。
#define Access(A,t,x,reim) ((float *)&A[t%BT][x%BX][reim])[ (x/BX)+NBX*(t/BT) ]
このように定義すると、Usimd[0][0][0]のレジスタ上には、以下の図の８点が入る事になります。
なぜ、こんな分割をするかと言うと、話が長くなるので、さっくり行かせる為に、ここでattributeから普通にAVX2や我らがXeonPhiのAVX-523系の演算を使ったプログラムに話を変える事にします。
次回、待望の_mm256_fmaddの登場です。