AVXで倍精度浮動小数点のでっかい行列の転置

#SIMD(AVX)で巨大な行列を転置する。
4*4の転置は[こちらの記事](http://qiita.com/nbjiao/items/db1a5ff2c0afafb5717a)をありがたく利用する。
10000*10000の行列を転置するようなケース。4*4のブロック単位で転置していく。
対角上のブロックは、そのまま転置。対角線に対して対称な位置にあるブロックどうしを転置して交換する。

```c:trans.c
+---+---+---+    +---+---+---+
| A | B | C |    | At| Dt| Gt|
+---+---+---+    +---+---+---+
| D | E | F | => | Bt| Et| Ht|
+---+---+---+    +---+---+---+
| G | H | I |    | Ct| Ft| It|
+---+---+---+    +---+---+---+
/*
1つのブロックは4*4。
AEIはただ転置。BとD、CとG、FとHは転置して交換する。
*/
void  transpose_pd(double **mat, int row, int col)
{
  // rowは4の倍数、row==colを前提とする。
  int i,j;
  for (i=0; i<row; i+=4) {
    j=i;
    __m256d ymm_mat0 = _mm256_loadu_pd(&mat[i+0][j]); // |03|02|01|00|
    __m256d ymm_mat1 = _mm256_loadu_pd(&mat[i+1][j]); // |13|12|11|10|
    __m256d ymm_mat2 = _mm256_loadu_pd(&mat[i+2][j]); // |23|22|21|20|
    __m256d ymm_mat3 = _mm256_loadu_pd(&mat[i+3][j]); // |33|32|31|30|

    __m256d ymm_mat01_lo = _mm256_unpacklo_pd(ymm_mat0, ymm_mat1); // |12|02|10|00|
    __m256d ymm_mat01_hi = _mm256_unpackhi_pd(ymm_mat0, ymm_mat1); // |13|03|11|01|
    __m256d ymm_mat23_lo = _mm256_unpacklo_pd(ymm_mat2, ymm_mat3); // |32|22|30|20|
    __m256d ymm_mat23_hi = _mm256_unpackhi_pd(ymm_mat2, ymm_mat3); // |33|23|31|21|
    __m256d ymm_swap0 = _mm256_permute2f128_pd(ymm_mat01_lo, ymm_mat23_lo, 0b00100000); // |30|20|10|00|
    __m256d ymm_swap1 = _mm256_permute2f128_pd(ymm_mat01_hi, ymm_mat23_hi, 0b00100000); // |31|21|11|01|
    __m256d ymm_swap2 = _mm256_permute2f128_pd(ymm_mat01_lo, ymm_mat23_lo, 0b00110001); // |32|22|12|02|
    __m256d ymm_swap3 = _mm256_permute2f128_pd(ymm_mat01_hi, ymm_mat23_hi, 0b00110001); // |33|23|13|03|

    _mm256_storeu_pd(&mat[i+0][j], ymm_swap0);
    _mm256_storeu_pd(&mat[i+1][j], ymm_swap1);
    _mm256_storeu_pd(&mat[i+2][j], ymm_swap2);
    _mm256_storeu_pd(&mat[i+3][j], ymm_swap3);
    for (j=i+4; j<col; j+=4) {
      __m256d ymm_matA0 = _mm256_loadu_pd(&mat[i+0][j]); // |03|02|01|00|
      __m256d ymm_matA1 = _mm256_loadu_pd(&mat[i+1][j]); // |13|12|11|10|
      __m256d ymm_matA2 = _mm256_loadu_pd(&mat[i+2][j]); // |23|22|21|20|
      __m256d ymm_matA3 = _mm256_loadu_pd(&mat[i+3][j]); // |33|32|31|30|
      __m256d ymm_matA01_lo = _mm256_unpacklo_pd(ymm_matA0, ymm_matA1); // |12|02|10|00|
      __m256d ymm_matA01_hi = _mm256_unpackhi_pd(ymm_matA0, ymm_matA1); // |13|03|11|01|
      __m256d ymm_matA23_lo = _mm256_unpacklo_pd(ymm_matA2, ymm_matA3); // |32|22|30|20|
      __m256d ymm_matA23_hi = _mm256_unpackhi_pd(ymm_matA2, ymm_matA3); // |33|23|31|21|
      __m256d ymm_swapA0 = _mm256_permute2f128_pd(ymm_matA01_lo, ymm_matA23_lo, 0b00100000); // |30|20|10|00|
      __m256d ymm_swapA1 = _mm256_permute2f128_pd(ymm_matA01_hi, ymm_matA23_hi, 0b00100000); // |31|21|11|01|
      __m256d ymm_swapA2 = _mm256_permute2f128_pd(ymm_matA01_lo, ymm_matA23_lo, 0b00110001); // |32|22|12|02|
      __m256d ymm_swapA3 = _mm256_permute2f128_pd(ymm_matA01_hi, ymm_matA23_hi, 0b00110001); // |33|23|13|03|

      __m256d ymm_matB0 = _mm256_loadu_pd(&mat[j+0][i]);
      __m256d ymm_matB1 = _mm256_loadu_pd(&mat[j+1][i]);
      __m256d ymm_matB2 = _mm256_loadu_pd(&mat[j+2][i]);
      __m256d ymm_matB3 = _mm256_loadu_pd(&mat[j+3][i]);
      __m256d ymm_matB01_lo = _mm256_unpacklo_pd(ymm_matB0, ymm_matB1);
      __m256d ymm_matB01_hi = _mm256_unpackhi_pd(ymm_matB0, ymm_matB1);
      __m256d ymm_matB23_lo = _mm256_unpacklo_pd(ymm_matB2, ymm_matB3);
      __m256d ymm_matB23_hi = _mm256_unpackhi_pd(ymm_matB2, ymm_matB3);
      __m256d ymm_swapB0 = _mm256_permute2f128_pd(ymm_matB01_lo, ymm_matB23_lo, 0b00100000);
      __m256d ymm_swapB1 = _mm256_permute2f128_pd(ymm_matB01_hi, ymm_matB23_hi, 0b00100000);
      __m256d ymm_swapB2 = _mm256_permute2f128_pd(ymm_matB01_lo, ymm_matB23_lo, 0b00110001);
      __m256d ymm_swapB3 = _mm256_permute2f128_pd(ymm_matB01_hi, ymm_matB23_hi, 0b00110001);

      _mm256_storeu_pd(&mat[j+0][i], ymm_swapA0);
      _mm256_storeu_pd(&mat[j+1][i], ymm_swapA1);
      _mm256_storeu_pd(&mat[j+2][i], ymm_swapA2);
      _mm256_storeu_pd(&mat[j+3][i], ymm_swapA3);
      _mm256_storeu_pd(&mat[i+0][j], ymm_swapB0);
      _mm256_storeu_pd(&mat[i+1][j], ymm_swapB1);
      _mm256_storeu_pd(&mat[i+2][j], ymm_swapB2);
      _mm256_storeu_pd(&mat[i+3][j], ymm_swapB3);
    }
  }
  return;
}
```

##Cで書いたのと比べる
あんまり速くないな。1.5倍弱の高速化。

```c:trans.c
void transpose(double **mat, int row, int col)
{
  int i,j;
  for (i=0; i<row; i++) {
    for (j=i+1; j<col; j++) {
      double swap = mat[i][j];
      mat[i][j] = mat[j][i];
      mat[j][i] = swap;
    }
  }
  return;
}
```