導入

有限要素法で全体係数行列が疎行列となるときは疎行列向けのソルバを使用します．前回の記事で紹介した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の対応する要素の行列中での列インデックスと行インデックスを格納します．

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形式と言ったら上記のようにソート済みで重複もないことを前提とします．

ここで，非零成分の各配列がまず行についてソートされていることに注目すると，row_indの代わりに次のような配列row_ptrを採用しても行列$[K]$を復元できることに気づきます²．

row_ptr

    row_ptr = (/ 1, 4, 6, 7, 8, 11 /)

ここで，row_ptrのi番目の要素の値は行列$[K]$のi行目を第一列から数えていった時の最初の非零成分がvalの何要素目に格納されているかを指し示す配列です．つまり

row_ptr

    i行目の第一列から数えて最初の非零成分の値=val(row_ptr(i))

ということです．row_ptrを行ポインタを格納する配列と呼ぶことにします．
row_indの値の代わりにrow_ptrを採用しても行列$[K]$を復元できるというのは，以下の手順で行列のi行目を復元できるということです．

行列$[K]$の次元をnとします（例に挙げた行列の場合は6です）
もしi<nなら，行列のi行目の非零成分は左端から順に，val(row_ptr(i))～val(row_ptr(i+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は先に挙げた行列の場合

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を追加で格納する方の行ポインタの配列が採用されます³．

全体係数行列の(i,j)成分を更新したいとき

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

先ほどの説明を踏まえると，これは以下の手順で達成できることがわかります．

$[K]_{ij}$の値がvalのどの要素に格納されているか知りたいとします
$[K]_{ij}$の値は，val(row_ptr(i))～val(row_ptr(i+1)-1)の中のどこかに格納されているはずです
kをrow_ptr(i)からrow_ptr(i+1)-1の間で動かして，col_ind(k)==jとなった時のval(k)にアクセスしたかった$[K]_{ij}$の値が格納されています

雑に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)成分にはどのようにアクセスしたらよいのか，についてまとめました．

Hishimuma_t, "疎行列とベクトルを掛けたい貴方に Chapter06 疎行列の格納形式：COO", (https://zenn.dev/hishinuma_t/books/sparse-matrix-and-vector-product/viewer/coo) ↩
COO形式とCSR形式の比較やrow_indの代わりにrow_ptrを使うことによる利点などは他の記事ですでにいろいろと説明されているのでここでは触れません． ↩
Fortranの配列は1始まりなのでこのようになります．配列が0始まりの場合にはここでの説明を色々調整して，row_ptrの末尾にはnnzが格納されることになります． ↩