LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@sukimaengineer

SPICE系回路シミュレーター向けスパースLUソルバ「KLU」を使ってみる

Last updated at Posted at 2025-01-18

はじめに

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

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

目次

  1. KLUとは
  2. 前提
  3. KLUのインストール方法
  4. KLUでスパース行列を解いてみる (C++20例)
    1. サンプル行列
    2. CCS形式の配列 (Ap, Ai, Ax)
    3. C++20コード例
    4. コンパイルと実行
    5. 連立方程式の解
  5. 応用: 時間領域解析や反復解
  6. トラブルシューティング
  7. SUNDIALSとKLUの併用
    1. SUNDIALSとKLUの位置づけ
    2. SUNDIALSがKLUを使う流れ
    3. 回路シミュレーションでの利用例
    4. 注意点
    5. まとめ
  8. SUNDIALS+KLUのサンプルコード (C++20)
    1. サンプルの概要
    2. コード例
    3. ビルドと実行
    4. コード解説 (簡単に)
  9. よくある質問(FAQ)
  10. 参考リンク
  11. まとめ

KLUとは

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

  • 電気回路向けに最適化
    一般的なスパースソルバ(UMFPACK, SuperLUなど)より、回路行列特有の疎構造を効率的に扱う工夫があります。

  • オープンソース(LGPL)
    商用/非商用を問わず、ライセンスを守れば誰でも無償で利用可能です。

  • C言語 API を提供
    C/C++やFortran、Pythonなどから呼び出せます。Pythonでの利用例には scikit-umfpack がありますが、KLU固有のラッパーは一部のプロジェクトのみ提供されています。

前提

  • 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. ソースからビルドする場合

  1. SuiteSparse公式GitHub からソースを取得します。
  2. 解凍して 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 で繰り返し解く流れが一般的です。

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

  1. SUNDIALSのビルド時にKLUサポートを有効化

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

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

    • CVODE/IDA等の初期化後、線形ソルバとして KLU をアタッチする: 
      SUNMatrix A = SUNSparseMatrix(N, N, NNZ, CSC_MAT);
SUNLinearSolver LS = SUNLinSol_KLU(y, A);
IDASetLinearSolver(ida_mem, LS, A);

  3. 微分(代数)方程式を時間刻みに沿って求解

    • 反復計算でヤコビ行列を生成し、KLUで分解。
    • スパースかつ帯状構造を持つ行列なら効率的に解ける。

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

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

7.4. 注意点

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

11. まとめ

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

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

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

以上

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?