CSR形式で表現された行列の成分にアクセスする方法

# 導入
有限要素法で全体係数行列が疎行列となるときは疎行列向けのソルバを使用します．[前回の記事](https://qiita.com/chopsticks/items/0bcd0efb78a5ee14be60)で紹介したPARDISOも疎行列向けのソルバであり，疎行列のデータをCompressed Sparse Row(CSR)形式で要求します．CSR形式は行列を素直に二次元配列で記憶する方法と比べて非直感的な方法です．

本記事では，有限要素法を実装する際などにおけるCSR形式との付き合い方についてまとめようと思います．まずCSR形式とは何なのか，そしてCSR形式で表現された行列の$(i,j)$成分にどのようにアクセスすればよいのか，についてです．

# Compressed Sparse Row(CSR)形式とは
まずは，CSR形式についての説明から始めます．
行列$[K]$の成分のほとんどが零であるとき$[K]$は疎行列と呼ばれます．成分のほとんどが零なので，この係数行列を表現するのに二次元配列を確保するのはもったいないです．非零成分とその成分の行列中での位置だけを覚えておけば十分というのはすんなり納得できるかと思います．これはCoordinate(COO)形式と呼ばれる格納形式の考え方です．

具体的に次のような行列を考えます．
$$ 
    [K] = 
    \begin{bmatrix}
    2 & 2 & 0 & 3 & 0 & 0 \\\\
    0 & 1 & 0 & 1 & 0 & 0 \\\\
    0 & 0 & 4 & 0 & 0 & 0 \\\\
    0 & 0 & 0 & 0 & 1 & 0 \\\\
    0 & 0 & 0 & 1 & 0 & 2 \\\\
    0 & 0 & 0 & 3 & 0 & 1 \\\\
    \end{bmatrix}
$$
このような行列を，COO形式では三つの一次元配列`val`，`col_ind`および`row_ind`で表現します．
まず`val`には非零成分の値を格納します．そして`col_ind`と`row_ind`には，`val`の対応する要素の行列中での列インデックスと行インデックスを格納します．
```fortran:COO
    val     = (/ 2, 2, 3, 1, 1, 4, 1, 1, 2, 3, 1 /)
    col_ind = (/ 1, 2, 4, 2, 4, 3, 5, 4, 5, 6, 6 /)
    row_ind = (/ 1, 1, 1, 2, 2, 3, 4, 5, 5, 5, 6 /)
```
ここでは，それぞれの配列をまず行について，そのあと列について昇順にソートして表現しましたが，COO形式と呼ばれるもの一般にそのようなソートがされていることや，重複がないという保証はないことに注意が必要です[^coo]．しかし，本記事ではCOO形式と言ったら上記のようにソート済みで重複もないことを前提とします．

[^coo]:Hishimuma_t, "疎行列とベクトルを掛けたい貴方に Chapter06 疎行列の格納形式：COO", (https://zenn.dev/hishinuma_t/books/sparse-matrix-and-vector-product/viewer/coo)

ここで，非零成分の各配列がまず行についてソートされていることに注目すると，`row_ind`の代わりに次のような配列`row_ptr`を採用しても行列$[K]$を復元できることに気づきます[^hosoku]．
[^hosoku]:COO形式とCSR形式の比較や`row_ind`の代わりに`row_ptr`を使うことによる利点などは他の記事ですでにいろいろと説明されているのでここでは触れません．
```fortran:row_ptr
    row_ptr = (/ 1, 4, 6, 7, 8, 11 /)
```
ここで，`row_ptr`の`i`番目の要素の値は行列$[K]$の`i`行目を第一列から数えていった時の最初の非零成分が`val`の何要素目に格納されているかを指し示す配列です．つまり
```fortran:row_ptr
    i行目の第一列から数えて最初の非零成分の値=val(row_ptr(i))
```
ということです．`row_ptr`を行ポインタを格納する配列と呼ぶことにします．
`row_ind`の値の代わりに`row_ptr`を採用しても行列$[K]$を復元できるというのは，以下の手順で行列の`i`行目を復元できるということです．
1. 行列$[K]$の次元を`n`とします（例に挙げた行列の場合は`6`です）
1. もし`i<n`なら，行列の`i`行目の非零成分は左端から順に，`val(row_ptr(i))`～`val(row_ptr(i+1)-1)`です
1. もし`i==n`なら，行列の`i`行目の非零成分は左端から順に，`val(row_ptr(i))`～`val(size(val))`です

ここで，`i`が`i<n`のときと`i==n`のときで場合分けが必要となるのが気持ち悪いです．これは，`row_ptr`の定義を修正することで簡単に解決できます．まず，`size(val)`，すなわち行列の非零成分の総数を`nnz`とします．そして`row_ptr`の末尾に要素を一つ追加して，`nnz+1`を格納することにします．`row_ptr`は先に挙げた行列の場合
```fortran:row_ptr改
    row_ptr = (/ 1, 4, 6, 7, 8, 11, nnz+1 /)
```
となります．こうすることで，行列の`i`行目の復元方法は
- 行列の次元を`n`とします（例に挙げた行列の場合は`6`です）
- 行列の`i`行目の非零成分は左端から順に，`val(row_ptr(i))`～`val(row_ptr(i+1)-1)`です

のようにシンプルにすることができました．実際，CSR形式といったらこちらの，末尾に`nnz+1`を追加で格納する方の行ポインタの配列が採用されます[^index]．
[^index]:Fortranの配列は1始まりなのでこのようになります．配列が0始まりの場合にはここでの説明を色々調整して，`row_ptr`の末尾には`nnz`が格納されることになります．<!--世の中は広いため，配列のインデックスを0始まりにしたい人たちもいます．そのような人たちはここでの説明を色々と調整する必要があります．`col_ind`と`row_ptr`の値も1つずつずれて最後の要素には`nnz`が格納されることになります．-->

# 全体係数行列の(i,j)成分を更新したいとき
非線形FEMなどをする際には，解析の反復中に逐一全体係数行列を更新したくなります．でも全体係数行列が今はCSR形式で表現されています．このようなとき，どうやって$[K]_{ij}$の値にアクセスしたらよいでしょうか．という問題について考えます．

先ほどの説明を踏まえると，これは以下の手順で達成できることがわかります．
1. $[K]_{ij}$の値が`val`のどの要素に格納されているか知りたいとします
1.  $[K]_{ij}$の値は，`val(row_ptr(i))`～`val(row_ptr(i+1)-1)`の中のどこかに格納されているはずです
1.  `k`を`row_ptr(i)`から`row_ptr(i+1)-1`の間で動かして，`col_ind(k)==j`となった時の`val(k)`にアクセスしたかった$[K]_{ij}$の値が格納されています

雑にFortranコード風に書いてみると，
```fortran:HowToAccess
    do k=row_ptr(i),row_ptr(i+1)-1
        if(col(k)==j)then
            !val(k)が[K]_{ij}の値
        endif
    enddo
```
といった感じになると思います．

# まとめ
今回は，CSR形式とはどのようなものか，そしてCSR形式で行列が表現されているときにその行列の(i,j)成分にはどのようにアクセスしたらよいのか，についてまとめました．