SIMD演算を用いた高精度疎行列計算ソフトウェアの高速化

 物理シミュレーションでは大規模な疎行列に対する連立一次方程式の求解が求められる．このような連立一次方程式を解くにはKrylov部分空間法とよばれる解法を用いることが一般的である．Krylov部分空間法は連立一次方程式の近似解を求める反復法アルゴリズムの1つのカテゴリで，スカラ，ベクトル，疎行列に対する演算（基本演算とよぶ）の組み合わせによって構成される．Krylov部分空間法は丸め誤差の蓄積により反復回数が増大したり，近似解が求まらなくなったりする．

 このような問題を解決するため，プログラム全体や誤差の蓄積しやすい一部の演算を高精度化する手法がある．高精度演算は計算時間やメモリデータ量が多く必要になることや標準的なプログラミング環境から高精度演算を行うためのソフトウェアが少ないことが問題となる．現実的な計算時間で高精度演算を実現するために，まずは倍精度と比べて2倍の精度である4倍精度に着目した．倍精度数演算を組み合わせて4倍精度演算を実現できるDouble-Double精度（DD）演算は特殊なハードウェアを必要とせずに4倍精度演算を実現できる．

 本研究ではKrylov部分空間法を実装するために必要なDDの基本演算を高速化することを目的とした．高速化のアプローチとして汎用プロセッサでサポートされている並列化機構を効率的に活用することを考え，現在主流なアーキテクチャに実装されているFMA命令，SIMD命令，マルチスレッドを用いた際に効率的に計算を行えるメモリレイアウトや実装方法を検討した．プロセッサのさまざまな高速化機構を効率的に利用するためにはメモリレイアウトや実装を工夫することが重要である．一方で高速化した演算をユーザに利用してもらうためには独自のメモリレイアウトや実装を隠し，簡潔なインタフェースを提供することが必要になる．そのため高速化した基本演算を利用するためのインタフェースを設計した．

 疎行列はほとんどの要素が零の行列で，メモリ使用量を節約するために零要素を記憶せず，行番号や列番号のインデックスを用いた格納形式を用いる．一般的によく使われる格納形式であるCRS形式はSIMD命令を用いた並列化の効率が悪いことが実験によって明らかになったため，SIMD命令を効率的に利用するためにBCRS形式とよばれる形式を用いた手法を提案し，CRS形式と比べて平均で約2.61倍高速に疎行列ベクトル積を計算できるという結果を得た．

 実装したこれらの演算をユーザから簡単に利用できるようにするため，倍精度型とDD型を組み合わせて利用できるインタフェースを設計した．開発したソフトウェアを利用することで変数の宣言部を変更するだけでプログラム全体の計算精度を変更でき，簡単に高速な高精度Krylov部分空間法を実現できることを示した．