Common Lisp で行列計算するブログを追試してみた (最後)

VOPでチューニングした方をアンローリングで最適化した。

# アンローラーの書き直し

`dotimes-unroll` をもう一度書きなおした。
元の実装から、ベタ書きに変換する `dotimes-inline` をまず分割して、
残りがオフセット計算からアンロールを行う `dotimes-unroll3` になった。

```cl
(defmacro dotimes-inline ((var count &optional result-form) &body body &environment env)
  (check-type var symbol)
  (let ((count (macroexpand count env)))
    (assert (and (constantp count) (numberp count)))
    (iter (for c to count)
          (when (first-iteration-p)
            (collect 'progn))
          (collect
              (if (< c count)
                  `(symbol-macrolet ((,var ,c))
                     ,@body)
                  result-form)))))
```

`dotimes-inline` は、 `count` が定数でないとエラーを吐く。シンボルマクロも認識して、その展開結果が定数ならOK. SBCLのIR1をもう少しよく理解すれば、変数が定数畳み込みされた時にも展開できるんだろうけど、それはめんどうなのでやめた。

```cl
(defmacro dotimes-unroll3 (((var &optional (delta (gensym "D")))
                            (count unroll)
                            &optional result)
                           (&body update)
                           &body body
                           &environment env)
  (with-gensyms (quat mod nth-loop)
    (let ((unroll (macroexpand unroll env)))
      `(multiple-value-bind (,quat ,mod) (floor ,count ,unroll)
         (declare (type fixnum ,quat)
                  (type (mod ,unroll) ,mod))
         (symbol-macrolet ((,delta ,unroll))
           (dotimes (,nth-loop ,quat)
             (declare (ignorable ,nth-loop))
             (dotimes-inline (,var ,unroll)
               ,@body)
             ,@update))
         (symbol-macrolet ((,delta 1)
                           (,var 0))
           (dotimes (,nth-loop ,mod)
             (declare (type (mod ,unroll) ,var))
             ,@body
             ,@update))
         ,result))))
```

`dotimes-unroll3` は、今までの`dotimes-unroll` とちがって `update` フォームを指定できる。 これはアンロールの際のオフセットを別に計算できるようにするため。次の例を見ればわかると思う。

# dotimes-unroll3 で2段アンローリング

もとの手でやったアンローリングがこれだ。colsが偶数なのを仮定してるのでコレでただしい。でも汎用にするには余りを考えないといけない。

```cl
(dotimes (col (/ cols 2))
  (sb-kernel:data-vector-set-with-offset
   (sb-kernel:%array-data-vector mc) mc-index 0
   (+ (sb-kernel:data-vector-ref-with-offset
       (sb-kernel:%array-data-vector mc) mc-index 0)
      (* cell (sb-kernel:data-vector-ref-with-offset
               (sb-kernel:%array-data-vector mb) mb-index 0))))
  (sb-kernel:data-vector-set-with-offset
   (sb-kernel:%array-data-vector mc) mc-index 1
   (+ (sb-kernel:data-vector-ref-with-offset
       (sb-kernel:%array-data-vector mc) mc-index 1)
      (* cell (sb-kernel:data-vector-ref-with-offset
               (sb-kernel:%array-data-vector mb) mb-index 1))))
  (incf mb-index 2)
  (incf mc-index 2))
```

アンロールの余りを計算するために抽象化するとこう。

```cl
  (multiple-value-bind (quat mod) (floor cols uj)
    (dotimes (col quat)
      (dotimes-inline (offset uj)
        (sb-kernel:data-vector-set-with-offset
         (sb-kernel:%array-data-vector mc) mc-index offset
         (+ (sb-kernel:data-vector-ref-with-offset
             (sb-kernel:%array-data-vector mc) mc-index offset)
            (* cell (sb-kernel:data-vector-ref-with-offset
                     (sb-kernel:%array-data-vector mb) mb-index offset)))))
      (incf mb-index uj)
      (incf mc-index uj))
    (dotimes (col mod)
      (sb-kernel:data-vector-set-with-offset
       (sb-kernel:%array-data-vector mc) mc-index 0
       (+ (sb-kernel:data-vector-ref-with-offset
           (sb-kernel:%array-data-vector mc) mc-index 0)
          (* cell (sb-kernel:data-vector-ref-with-offset
                   (sb-kernel:%array-data-vector mb) mb-index 0))))
      (incf mb-index)
      (incf mc-index)))
```

オフセットの更新部分を `update` に移して、マクロでやるとこう。

```cl
(dotimes-unroll3 ((col delta) (cols 2))
    ((incf mb-index delta)
     (incf mc-index delta))
  (sb-kernel:data-vector-set-with-offset
   (sb-kernel:%array-data-vector mc) mc-index col
   (+ (sb-kernel:data-vector-ref-with-offset
       (sb-kernel:%array-data-vector mc) mc-index col)
      (* cell (sb-kernel:data-vector-ref-with-offset
               (sb-kernel:%array-data-vector mb) mb-index col))))))))))
```

# 結果まとめ

最終回なので、今までの性能のまとめを書いておこうとおもう。
テストの繰り返し回数が変わっているので前回と数字が違うが傾向は一緒。
行列はもういいや。

|                                | time [sec] |
|--------------------------------+------------:|
| **IJK-loop**                     |            |
| SIMPLE-GEMM                    |       22.1 |
| CACHE-GEMM                     |     14.008 |
| ROW-MAJOR-GEMM                 |      9.412 |
| RM-GEMM+STATIC-SIZE            |      8.428 |
| RM-GEMM+STATIC-SIZE+UNROLL2    |       7.98 |
| RM-GEMM+STATIC-SIZE+UNROLL8    |      7.612 |
| autotune (1 2 16)              |      6.736 |
|                                |            |
| **IKJ-loop**                     |            |
| SIMPLE-GEMM-K                  |     17.924 |
| CACHE-GEMM-K                   |     14.588 |
| RM-GEMM-K                      |       8.04 |
| RM-GEMM+STATIC-SIZE-K          |      8.108 |
| RM-GEMM+STATIC-SIZE+UNROLL2-K1 |      8.332 |
|                                |            |
| **IKJ-loop + VOP**               |            |
| RM-GEMM+STATIC-SIZE+UNROLL2-K2 |      6.052 |
| RM-GEMM+STATIC-SIZE+UNROLL8-K  |      6.232 |
| autotune (1 2 16)              |      5.632 |

コードはここ。

https://github.com/guicho271828/play-on-matrix