sse3 メモ

＊複素数計算を SSE3 で行う。
  Fortran と C を混ぜて使う。
  SSE3 の計算は intel C の intrinsic 命令を使う。


=================================================================================
(倍精度)
=================================================================================
＊複素数の対応
fortran 側、複素数
  complex(8) :: a

は

C 側 構造体
  typedef struct {
    double re;
    double im;
  } cmplx;

  cmplx *a;

と対応させる。

＊ C 側での複素数計算
 (0)変数の宣言

  cmplx *a,*b;          // cmplx 構造体へのポインタ , ２つの double(64bit) を持つ
  __m128d  t1,t2,t3;    //  128bit SSE2 レジスタ変数, ２つの double(64bit) を持つ
  
 (1)変数のロード

   t1 = _mm_load_pd((double *)a);   
または
   t1 = *(__m128d *)a;
   // アドレス a から 2 つの 64bit double を SSE2 レジスタ t1 にロード

 (2)変数のストア

   _mm_store_pd((double *)a,t1);   
または、
   *(__m128d *)a = t1;
   // SSE2 レジスタ t1 の内容を アドレス a へ 2 つの 64bit double をストア

 (3)２つの複素数の和(2-comp. vector add)

   t3 = _mm_add_pd(t2,t1);  
   // t3 = t2 + t1;    t1,t2,t3 はそれぞれ２成分ベクトルと考える。

 (4)２つの複素数の差(2-comp. vector sub)

   t3 = _mm_sub_pd(t2,t1);  
   // t3 = t2 - t1;    t1,t2,t3 はそれぞれ２成分ベクトルと考える。

 (5)２つの複素数の計算。実部は差を取り、虚部は和を取る(SSE3のみ)

   t3 = _mm_addsub_pd(t2,t1);  
   // t3.re = t2.re - t1.re;    
   // t3.im = t2.im + t1.im;    t1,t2,t3 はそれぞれ複素数と考える。

 (6)２つの sse2 レジスタの積

   t3 = _mm_mul_pd(t2,t1);  
   // t3.re = t2.re*t1.re;    
   // t3.im = t2.im*t1.im;      t1,t2,t3 はそれぞれ複素数と考える。

 (7) 実部を取り出す

   t2 = _mm_unpacklo_pd(t1,t1);  
   // t2.re = t1.re;    
   // t2.im = t1.re;

 (8) 虚部を取り出す

   t2 = _mm_unpackhi_pd(t1,t1);  
   // t2.re = t1.im;    
   // t2.im = t1.im;

 (9) 実部と虚部を入れ換える

   t2 = _mm_shuffle_pd(t1,t1,1);  
   // t2.re = t1.im;    
   // t2.im = t1.re;

(10) 実部の符号を反転する。

   t2 = _mm_xor_pd(t1,_mm_set_pd(0.0,-0.0));  
または
   t2 = _mm_xor_pd(t1,_mm_set_sd(-0.0));  
   // t2.re =-t1.re;
   // t2.im = t1.im;

(11) 虚部の符号を反転する。

   t2 = _mm_xor_pd(t1,_mm_set_pd(-0.0,0.0));  
   // t2.re = t1.re;
   // t2.im =-t1.im;


=================================================================================
(単精度)
=================================================================================
＊複素数の対応
fortran 側、複素数
  complex(4) :: a(2)

は

C 側 構造体
  typedef struct {
    float re;
    float im;
  } fcmplx;

  fcmplx a[2];

と対応させる。

＊ C 側での複素数計算
 (0)変数の宣言

  fcmplx  a[2],b[2];    //  fcmplx[2], ４つの  float(32bit) を持つ
  __m128   t1,t2,t3;    //  128bit SSE2 レジスタ変数, ４つの float(32bit) を持つ
  
 (1)変数のロード

   t1 = _mm_load_ps((float *)a);   
または、
   t1 = *(__m128 *)a;
   // アドレス a から 4 つの 32bit float を SSE2 レジスタ t1 にロード

 (2)変数のストア

   _mm_store_ps((float *)a,t1);   
または、
   *(__m128 *)a = t1;
   // SSE2 レジスタ t1 の内容を アドレス a へ ４ つの 32bit float をストア

 (3)２つの複素数x2 の和(4-comp. vector add)

   t3 = _mm_add_ps(t2,t1);  
   // t3 = t2 + t1;    t1,t2,t3 はそれぞれ４成分ベクトルと考える。

 (4)２つの複素数x2の差(4-comp. vector sub)

   t3 = _mm_sub_ps(t2,t1);  
   // t3 = t2 - t1;    t1,t2,t3 はそれぞれ４成分ベクトルと考える。

 (5)２つの複素数x2の計算。実部は差を取り、虚部は和を取る(SSE3のみ)

   t3 = _mm_addsub_ps(t2,t1);  
   // t3_1.re = t2_1.re - t1_1.re;    
   // t3_1.im = t2_1.im + t1_1.im;
   // t3_2.re = t2_2.re - t1_2.re;    
   // t3_2.im = t2_2.im + t1_2.im;    t1,t2,t3 はそれぞれ複素数x2と考える。

 (6)２つの sse2 レジスタの積

   t3 = _mm_mul_ps(t2,t1);  
   // t3_1.re = t2_1.re*t1_1.re;    
   // t3_1.im = t2_1.im*t1_1.im;
   // t3_2.re = t2_2.re*t1_2.re;    
   // t3_2.im = t2_2.im*t1_2.im;      t1,t2,t3 はそれぞれ複素数x2と考える。

 (7) 実部を取り出す

   t2 = _mm_moveldup_ps(t1);    //  SSE3 のみ  
または
   t2 = _mm_shuffle_ps(t1,t1,_MM_SHUFFLE(2,2,0,0));    //   SSE2
   // t2_1.re = t1_1.re;    
   // t2_1.im = t1_1.re;
   // t2_2.re = t1_2.re;    
   // t2_2.im = t1_2.re;

 (8) 虚部を取り出す

   t2 = _mm_movehdup_ps(t1);    //   SSE3 のみ
または
   t2 = _mm_shuffle_ps(t1,t1,_MM_SHUFFLE(3,3,1,1));    //   SSE2
   // t2_1.re = t1_1.im;
   // t2_1.im = t1_1.im;
   // t2_2.re = t1_2.im;
   // t2_2.im = t1_2.im;

 (9) 実部と虚部を入れ換える

   t2 = _mm_shuffle_ps(t1,t1,_MM_SHUFFLE(2,3,0,1));  
   // t2_1.re = t1_1.im;    
   // t2_1.im = t1_1.re;
   // t2_2.re = t1_2.im;    
   // t2_2.im = t1_2.re;

(10) 実部の符号を反転する。

   t2 = _mm_xor_ps(t1,_mm_set_ps(0.0,-0.0,0.0,-0.0));  

   // t2_1.re =-t1_1.re;
   // t2_1.im = t1_1.im;
   // t2_2.re =-t1_2.re;
   // t2_2.im = t1_2.im;

(11) 虚部の符号を反転する。

   t2 = _mm_xor_ps(t1,_mm_set_ps(-0.0,0.0,-0,0,0,0));  
   // t2_1.re = t1_1.re;
   // t2_1.im =-t1_1.im;
   // t2_2.re = t1_2.re;
   // t2_2.im =-t1_2.im;