はじめに

この記事では、電気回路シミュレーション(SPICE系)などで大規模なスパース行列を解く際に有用な KLU ソルバについて解説します。具体的には、KLUの概要からインストール方法、簡単な線形方程式のサンプルプログラム、そしてSUNDIALSとの併用例まで紹介します。さらに、よくある質問(FAQ)を追加して理解を深められるよう工夫しています。

本記事では、サンプルコードをC++20に対応した形で作成しています。出力部分に std::format (C++20) を使ったフォーマット表示を取り入れています。（利用にはC++20以降をサポートし、かつ <format> を実装している標準ライブラリが必要です。）

KLUとは

KLU (Sparse LU Solver for Circuit Simulators) は、フロリダ大学の Timothy A. Davis 氏が開発した SuiteSparse ライブラリ群の一部です。とくに電気回路シミュレーション(SPICE系)で発生する行列のスパース性や帯状構造を活かし、高速かつメモリ効率の良いLU分解を実現します。

前提

OS/ディストリビューション: 本記事では主に Linux (Ubuntu系) を例に説明
コンパイラ: C++20に対応した gcc / clang (C++20の <format> が使用可能であること)
SuiteSparse のバージョン: 7.x系 (執筆時点)

※WindowsやmacOSの場合でも、CMakeや自前のビルドスクリプトでコンパイルすれば、ほぼ同様の手順で利用できます。インストール先やビルドオプションを適宜調整してください。
※std::format を利用するには、GCC 12以上やClang 13以上など、C++20の実装が揃ったコンパイラが必要な場合が多いです。GCC 10や11では <format> ヘッダが未実装または不完全の場合があります。

1. KLUのインストール方法

1-1. Linux(Ubuntu)でAPTを使う場合

Ubuntuなどのディストリビューションによっては、SuiteSparse一式をまとめてインストールできます。

sudo apt-get update
sudo apt-get install libsuitesparse-dev

このコマンドで、klu.h を含む SuiteSparse の各種ヘッダとライブラリ (libklu.so, libumfpack.so など) がインストールされます。

1-2. ソースからビルドする場合

SuiteSparse公式GitHub からソースを取得します。
解凍して SuiteSparse/ ディレクトリに移動したら、以下のように Makefile/CMake でビルドします。

make       # あるいは
cmake .
make
sudo make install

インストール先を指定したい場合は、環境変数 INSTALL などを指定する、またはCMakeの -DCMAKE_INSTALL_PREFIX=... オプションを使って調整してください。

2. KLUでスパース行列を解いてみる (C++20例)

2-1. サンプル行列

ここでは、3×3 の簡単な線形方程式を解いてみます。
例として、次の連立方程式を考えます。

\begin{cases}
2x + 1y + 0z = 5 \\
1x + 3y + 1z = 10 \\
0x + 1y + 2z = 8
\end{cases}

行列形式で書くと、

\begin{pmatrix}
2 & 1 & 0 \\
1 & 3 & 1 \\
0 & 1 & 2
\end{pmatrix}
\cdot
\begin{pmatrix}
x \\
y \\
z
\end{pmatrix}
=
\begin{pmatrix}
5 \\
10 \\
8
\end{pmatrix}

これを 圧縮列形式(Compressed Column Storage, CCS) で表し、KLUで LU 分解 & 解の計算を行います。

どうしてCCS形式が必要？

KLUなどのスパースソルバは、多くのゼロ要素を無視して計算負荷を下げるため、圧縮表現が必須です。CCS形式は“列優先”で非ゼロ要素の情報を持つため、行列要素へのアクセス・更新が効率化されます。

2-2. CCS形式の配列 (Ap, Ai, Ax)

CCS形式(別名CSC: Compressed Sparse Column)では、列ごとに「非ゼロ要素の行インデックス」と「値」をまとめます。具体的には、

Ap: 列ポインタ配列(サイズは 列数+1)
Ai: 行インデックス配列(非ゼロ要素数分)
Ax: 実際の値(非ゼロ要素数分)

今回の行列

\begin{pmatrix}
2 & 1 & 0 \\
1 & 3 & 1 \\
0 & 1 & 2
\end{pmatrix}

を列ごとに見ると:

列0: (行0, 2), (行1, 1)
列1: (行0, 1), (行1, 3), (行2, 1)
列2: (行1, 1), (行2, 2)

非ゼロ要素数は合計 7 つ。

Ap (列ポインタ)

std::array<int, 4> Ap = {0, 2, 5, 7};

Ai, Ax (行インデックスと値)

std::array<int, 7> Ai = {0, 1, 0, 1, 2, 1, 2};
std::array<double, 7> Ax = {2, 1, 1, 3, 1, 1, 2};

b ベクトル

std::array<double, 3> b = {5, 10, 8};

2-3. C++20コード例

// filename: example_klu.cpp
#include <format>    // C++20以降
#include <iostream>
#include <array>
#include "klu.h"     // インストール先に応じてパスを調整

int main()
{
    int n = 3;  // 行列サイズ (3x3)

    // CCS形式のデータ
    std::array<int, 4> Ap = {0, 2, 5, 7};
    std::array<int, 7> Ai = {0, 1, 0, 1, 2, 1, 2};
    std::array<double, 7> Ax = {2, 1, 1, 3, 1, 1, 2};

    // 右辺ベクトル
    std::array<double, 3> b = {5, 10, 8};

    klu_common Common;
    klu_defaults(&Common);

    // シンボリック解析
    klu_symbolic* Symbolic = klu_analyze(n, Ap.data(), Ai.data(), &Common);
    if (!Symbolic) {
        std::cerr << "klu_analyze failed.\n";
        return -1;
    }

    // 数値分解
    klu_numeric* Numeric = klu_factor(Ap.data(), Ai.data(), Ax.data(), Symbolic, &Common);
    if (!Numeric) {
        std::cerr << "klu_factor failed.\n";
        klu_free_symbolic(&Symbolic, &Common);
        return -1;
    }

    // 解を求める (x を b に上書き)
    int ok = klu_solve(Symbolic, Numeric, n, 1, b.data(), &Common);
    if (!ok) {
        std::cerr << "klu_solve failed.\n";
    } else {
        std::cout << std::format("Solution: x={}, y={}, z={}\n", b[0], b[1], b[2]);
    }

    // 後処理
    klu_free_numeric(&Numeric, &Common);
    klu_free_symbolic(&Symbolic, &Common);

    return 0;
}

2-4. コンパイルと実行

# SuiteSparse(KLU) ライブラリがインストールされている前提
# C++20モードでコンパイルし、リンク時にsuitesparse関連ライブラリを指定します

g++ -std=c++20 example_klu.cpp -o example_klu \
    -I/usr/include/suitesparse \
    -L/usr/lib/x86_64-linux-gnu \
    -lklu -lamd -lcolamd -lbtf -lsuitesparseconfig

# 実行
./example_klu
# => "Solution: x=..., y=..., z=..."

注意: std::format はC++20の仕様ですが、GCC 10/11では未実装または不完全です。GCC 12以降やClang 13以降など、実装が揃った環境でご利用ください。

2-5. 連立方程式の解

連立方程式

\begin{pmatrix}
2x + y = 5 \\
x + 3y + z = 10 \\
y + 2z = 8
\end{pmatrix}

を手計算で解くと、

x = 1.625, \ y = 1.75, \ z = 3.125.

プログラム実行時に得られる解も、(数値誤差を除いて) ほぼ同じ値になります。

3. 応用: 時間領域解析や反復解

SPICE系シミュレーターを作成する場合、修正ノード解析(MNA) により大規模かつスパースな行列を生成し、KLU で繰り返し解く流れが一般的です。

ネットリスト解析 -> MNA行列 (A) と右辺ベクトル (b) を形成
klu_analyze / klu_factor -> 初期LU分解(シンボリック + 数値分解)
klu_solve -> 解を得る
非線形素子(MOS, ダイオードなど)の更新 -> 行列再設定
必要に応じて klu_refactor (数値再分解のみ) などを呼び出し、ループを回す
時間ステップごとに行列更新 -> 再度 KLU で解く

KLUはフィルイン(新たな非ゼロ要素の発生)を最小化するための順序付け（AMD など）を自動で行ってくれるため、回路規模が大きくなっても比較的スケールしやすいのが特徴です。

4. トラブルシューティング

Symbolic や Numeric が NULL になる
- 行列データ((Ap, Ai, Ax)) が正しく設定されていないか、行列サイズ((n)) が合っていない可能性があります。
リンカエラー (undefined reference to 'klu_...')
- コンパイル時に -lklu のほか、依存する -lamd, -lcolamd, -lbtf, -lsuitesparseconfig 等もリンクしているか確認してください。
数値が発散する／NaNになる
- スパース行列が特異行列（ランク落ち）になっていないか、行列要素や回路モデルを見直しましょう。
- 回路の場合、ノードが完全に浮いている(接地がない等)と、行列が特異になる可能性があります。
Windowsでビルドできない
- SuiteSparse公式リポジトリには、CMakeファイルが用意されています。Visual Studio等でCMakeプロジェクトを生成し、同様にビルド可能です。

7. SUNDIALSとKLUの併用

SUNDIALS（CVODEやIDAなど）とKLUを「併用」することは十分に可能です。
SUNDIALS側では内部で生成されるヤコビアンや線形系を解く際に、「外部のスパース直接ソルバ」と連携する仕組みが用意されており、その一つとしてSuiteSparse(KLU)がサポートされています。

7.1. SUNDIALSとKLUの位置づけ

SUNDIALS
Lawrence Livermore National Laboratoryが開発している、常微分方程式(ODE) や 微分代数方程式(DAE) の数値解法ライブラリ群。
例: CVODE, CVODES, IDA, IDAS, KINSOL, ARKODE など
KLU
SuiteSparseに含まれるスパース行列専用のLU分解ソルバ。電気回路シミュレーション向けに最適化され、修正ノード解析(MNA) 系の行列を高速に解けます。

SUNDIALS自身は、内部で線形方程式を解く必要が生じたときにユーザーが選択した線形ソルバを呼び出し、ヤコビ行列(または近似ヤコビ行列)を解きます。ここにKLUを組み込むことで、大規模かつスパースな系を高速に解くことができます。

7.2. SUNDIALSがKLUを使う流れ

SUNDIALSのビルド時にKLUサポートを有効化
- -DKLU_ENABLE=ON などを指定し、SuiteSparse(KLU)へのパスを渡してビルド。
- これにより SUNDIALS が KLU用ラッパ(SUNLinSol_KLU など)をコンパイルし、使用可能となる。

ソルバの初期化時にKLUを指定

CVODE/IDA等の初期化後、線形ソルバとして KLU をアタッチする:

SUNMatrix A = SUNSparseMatrix(N, N, NNZ, CSC_MAT);
SUNLinearSolver LS = SUNLinSol_KLU(y, A);
IDASetLinearSolver(ida_mem, LS, A);

微分(代数)方程式を時間刻みに沿って求解
- 反復計算でヤコビ行列を生成し、KLUで分解。
- スパースかつ帯状構造を持つ行列なら効率的に解ける。

7.3. 回路シミュレーションでの利用例

SPICE系シミュレーションは、連立常微分方程式(または微分代数方程式) を時間領域で解く過程がある。
MNAで得られる行列はスパースで、規模が大きいほど効率的なソルバが必須。
SUNDIALS(IDA等)とKLUを組み合わせることで、大規模・非線形回路の時間領域解析を効率化可能。

7.4. 注意点

KLUだけでは時間積分はできない
- ODE/DAEの解法はSUNDIALS等が担当。KLUは線形方程式のソルバ。
SUNDIALSには他のソルバ選択肢もある
- 反復ソルバ(SPFGMRなど)や密行列ソルバ等。ただし大規模問題でスパース行列が多い場合はKLUが有利なケースが多い。
行列更新の頻度
- 非線形反復ごとに行列を再構築・分解するため、問題の規模や構造次第でKLUの効果が変わる。

7.5. まとめ

SUNDIALS (CVODE/IDAなど) と KLU は、回路シミュレーションを含む大規模スパース問題でよく併用される組み合わせです。
特に回路行列のような帯状スパース構造を活かして、時間領域解析を効率よく行うことが可能になります。

8. SUNDIALS+KLUのサンプルコード (C++20)

8.1. サンプルの概要

SUNDIALS の IDA(微分代数方程式ソルバ) と KLU を連携させる簡単な例
2次元の微分代数方程式(DAE)を解く
IDA内部でKLUが呼ばれ、ヤコビ行列をLU分解する

8.2. コード例

// filename: example_sundials_klu.cpp
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <format>
#include <iostream>

// =========== SUNDIALS関連のヘッダ ==========
#include <idas/idas.h>                   // IDACreate, IDAInitなど
#include <nvector/nvector_serial.h>      // N_VNew_Serialなど
#include <sunmatrix/sunmatrix_sparse.h>  // SUNSparseMatrix
#include <sunlinsol/sunlinsol_klu.h>     // SUNLinSol_KLU
#include <idas/idas_ls.h>               // IDASetLinearSolver

// 定数 (実数型)
#ifndef REALTYPE
#define REALTYPE double
#endif

// ---------- 2次元のDAE(残差)関数 -----------
// r_1 = y1'(t) + 2y1(t) - cos(t)
// r_2 = y2'(t) + 1y2(t) - sin(t)
static int resfn(REALTYPE t,
                 N_Vector yy,   
                 N_Vector yp,   
                 N_Vector rr,   
                 void* user_data)
{
    const realtype* y  = NV_DATA_S(yy);
    const realtype* yp_ = NV_DATA_S(yp);
    realtype* r = NV_DATA_S(rr);

    r[0] = yp_[0] + 2.0*y[0] - std::cos(t);
    r[1] = yp_[1] + 1.0*y[1] - std::sin(t);
    return 0;
}

int main()
{
    // ---- 変数数 N=2 ----
    const int N = 2;

    // ---- IDAメモリ作成 ----
    void* ida_mem = IDACreate();
    if (!ida_mem) {
        std::cerr << "Error: IDACreate failed.\n";
        return -1;
    }

    // ---- 初期状態 (y, y') ----
    N_Vector y  = N_VNew_Serial(N);
    N_Vector yp = N_VNew_Serial(N);
    if (!y || !yp) {
        std::cerr << "Error: N_VNew_Serial failed.\n";
        return -1;
    }

    NV_Ith_S(y, 0)  = 0.0;
    NV_Ith_S(y, 1)  = 0.0;
    NV_Ith_S(yp, 0) = 1.0;  // cos(0) - 2*0
    NV_Ith_S(yp, 1) = 0.0;  // sin(0) - 1*0

    // ---- IDA初期化 ----
    int flag = IDAInit(ida_mem, resfn, 0.0, y, yp);
    if (flag != IDA_SUCCESS) {
        std::cerr << "Error: IDAInit failed.\n";
        return -1;
    }

    // ---- 誤差許容値 ----
    flag = IDASStolerances(ida_mem, 1e-6, 1e-6);
    if (flag != IDA_SUCCESS) {
        std::cerr << "Error: IDASStolerances failed.\n";
        return -1;
    }

    // ---- KLU用の行列 & ソルバ ----
    SUNMatrix A = SUNSparseMatrix(N, N, 4, CSC_MAT); // 2x2で最大4要素
    if (!A) {
        std::cerr << "Error: SUNSparseMatrix creation failed.\n";
        return -1;
    }

    SUNLinearSolver LS = SUNLinSol_KLU(y, A);
    if (!LS) {
        std::cerr << "Error: SUNLinSol_KLU creation failed.\n";
        return -1;
    }

    flag = IDASetLinearSolver(ida_mem, LS, A);
    if (flag != IDA_SUCCESS) {
        std::cerr << "Error: IDASetLinearSolver failed.\n";
        return -1;
    }

    // ---- 時間ループ (0 ~ 2.0) ----
    realtype t = 0.0;
    realtype tout = 0.1;
    while (tout <= 2.0) {
        flag = IDASolve(ida_mem, tout, &t, y, yp, IDA_NORMAL);
        if (flag < 0) {
            std::cerr << "Error: IDASolve failed at t=" << t << "\n";
            break;
        }
        double y1 = NV_Ith_S(y, 0);
        double y2 = NV_Ith_S(y, 1);
        std::cout << std::format("t={:.2f}, y1={:.6f}, y2={:.6f}\n", 
                                 (double)t, y1, y2);
        tout += 0.1;
    }

    // ---- 後処理 ----
    N_VDestroy(y);
    N_VDestroy(yp);
    SUNLinSolFree(LS);
    SUNMatDestroy(A);
    IDAFree(&ida_mem);

    return 0;
}

8.3. ビルドと実行

# 例: g++でC++20, SUNDIALS+KLUライブラリをリンク (一例)
g++ -std=c++20 example_sundials_klu.cpp -o example_sundials_klu \
    -I/usr/include \
    -L/usr/lib/x86_64-linux-gnu \
    -lsundials_ida -lsundials_nvecserial -lsundials_sunlinsollklu -lsundials_sunmatrixsparse \
    -lsuitesparseconfig -lklu -lamd -lcolamd -lbtf \
    -lm

./example_sundials_klu
# => 時間ステップごとに y1, y2 が出力される

8.4. コード解説 (簡単に)

resfn: 残差関数 $r(t, y, y') = 0$ を定義
IDAInit: IDAソルバに初期条件と残差関数を登録
SUNSparseMatrix, SUNLinSol_KLU: KLU用の行列＆ソルバを用意し、IDAへ設定
IDASolve: 時間(t)を進めながら、内部でヤコビ行列をKLUがLU分解して反復を行い、$y$ と $y'$ を更新

9. よくある質問(FAQ)

Q1. KLUはどんな場面で特に効果を発揮するのですか？

A. KLUは特に電気回路シミュレーションで発生する行列構造を想定しています。回路行列はスパースかつ帯状構造を持つ場合が多く、汎用のLU分解ソルバより高速・省メモリに計算できるケースが多いです。

Q2. KLUを他のスパースソルバ(UMFPACK, SuperLU)ではなく選ぶ理由は？

A. KLUは回路特有の行列構造に合わせた最適化が行われているため、同じ問題規模・同じ行列でもフィルインが少なく済み、分解に要するメモリ量や計算量を抑えられます。汎用ソルバより高速になる可能性が高いです。

Q3. KLUはオープンソースプロジェクトだけで使われているのですか？

A. SuiteSparse自体はLGPLライセンスで公開されていますが、商用プロジェクトでもライセンスを遵守すれば使用可能です。商用Spice系シミュレーターでは独自最適化ソルバを使っている場合が多いですが、研究・教育目的や個人・社内ツールなどで KLU が活用される例は数多くあります。

Q4. CCS形式(別名CSC形式)以外の形式は使えないのですか？

A. KLUはCSC形式(列圧縮)を基本として実装されているため、それを満たす入力データが必要です。CSR形式(行圧縮)からCSC形式へ変換することも可能なので、変換実装を挟めばほかのスパース形式からも問題なく利用できます。

Q5. 数値計算の安定性はどうですか？Pivoting(ピボット選択)は行われますか？

A. 部分ピボット法や列ピボット選択を採用しており、数値安定性にも配慮されています。ただし、KLUは回路シミュレーションを想定しているため、ある程度「回路らしい」行列を仮定しています。極端に不安定な行列の場合は、他の手法(完全ピボット選択など)が必要になるケースもありえます。

Q6. メモリ使用量や速度はどのくらいですか？

A. 問題の規模や行列の構造によって異なりますが、KLUは電気回路由来のスパース行列(帯状・局所的接続)に強い最適化が施されています。結果として、同規模の一般スパースソルバ(UMFPACK, SuperLU など)よりフィルインが少ない分、有利になる場合が多いです。ただし、行列の構造によっては効果が薄いこともあります。

Q7. 時間領域解析でKLUの分解を再利用できますか？

A. KLUには再数値分解 (klu_refactor) の仕組みがあります。一度シンボリック解析を行ったあと、行列の要素だけが変わる場合は再度の数値分解だけで済むため、高速化が期待できます。ただし、大幅に回路構造(行列パターン)が変わる場合は再度シンボリック解析が必要です。

Q8. SUNDIALS+KLUをGPU上で動かせますか？

A. 現在、KLUはCPU向けの実装が中心であり、SUNDIALSのKLU連携部分もCPU計算が前提となっています。GPU対応のスパースLUソルバやSUNDIALSオプションもありますが、KLU自身はまだGPUを直接サポートしていません。将来的にGPU向け最適化が進む可能性はありますが、現段階ではCPU上での使用が主流です。

10. 参考リンク

SuiteSparse 公式サイト
SuiteSparse GitHubリポジトリ
Ngspice + KLU (NgspiceでKLUをオプション使用する例)
Xyce (大型回路シミュレーターでKLUを採用)
SUNDIALS 公式ページ
SUNDIALS GitHubリポジトリ

11. まとめ

本記事では、KLU を使ってスパース行列を解くための手順とサンプルプログラム(C++20)を紹介しました。

KLU は回路シミュレーター向けに特化した設計で、高速・省メモリなLU分解を実現
CCS形式 $(Ap, Ai, Ax)$ を正しく用意すれば、C++20 で簡単にLU分解&解の計算が可能
SPICEのような 非線形問題 や 時間領域解析 でも、KLUを繰り返し活用可能
SUNDIALSと組み合わせれば、時間積分を伴う常/微分代数方程式(ODE/DAE)の解法を効率化できる

もしスパース行列の扱いに慣れていない方でも、一度CCS形式の概念を理解すれば、大規模回路でもKLUで安定して解くことができます。ぜひ試してみてください。

以上