(57)【要約】回路シミュレーションにおいて並列方程式を解くためのコンピュータ実装された方法を説明する。この方法は、回路ヤコビ行列を疎結合のパーティションに分割するステップと、電圧ベクトルと行列をパーティションにしたがってリオーダーするステップと、ヤコビ行列を二つの行列ＭとＮにスプリットするステップとを含み、Ｍは並列処理に適した行列であり、Ｎは結合行列である。ＭとＮは、Ｍ−１Ｊｘ＝（Ｉ＋Ｍ−１Ｎ）ｘ＝Ｍ−１ｒを形成するようにプレコンディションされ、ヤコビ行列Ｊは反復解法を用いて解かれる。【選択図】図４
2009520306.tif 000002
【特許請求の範囲】
【請求項１】
集積回路のシミュレーションオペレーションのためのコンピュータ実装された方法であって、
集積回路のオペレーションをモデル化する微分代数方程式（ＤＡＥ）の系を生成するステップと、
前記ＤＡＥの系を離散化するステップと、
前記ＤＡＥが非線形である場合、離散化されたＤＡＥを線形化して回路ヤコビ行列をもつ線形系を形成するステップと、
線形化された回路ヤコビ行列ソルバーを用いて前記線形系を解くステップとを含み、
前記線形系を解くステップは、
前記回路ヤコビ行列を二つの行列ＭとＮにスプリットするステップと、
対角成分が１で非対角成分がゼロの恒等行列をＩとして、前記二つの行列をプレコンディションしてＩ＋Ｍ−１Ｎの形の行列をもつプレコンディションされた方程式を形成するステップと、
直接・反復組み合わせ解法を用いて前記プレコンディションされた方程式におけるＩ＋Ｍ−１Ｎに対する解を求めるステップとを含むことを特徴とする方法。
【請求項２】
前記ＤＡＥの系を生成するステップは、集積回路に対する設計をＤＡＥの系にマッピングするステップを含み、ＤＡＥの系は集積回路の複数のノードのそれぞれにおける電圧の変化を定義し、前記電圧の変化は定義可能な期間にわたって発生するものであることを特徴とする請求項１の方法。
【請求項３】
前記ＤＡＥの系を解くステップはさらに、
前記解を前記離散化された方程式に代入して複数のノードのそれぞれに対して電圧の変化を求めるステップと、
前記ＤＡＥの有限差分離散化と非線形反復法を用いて電圧の変化を解くための前記線形系を求めるステップと、
複数のノードのそれぞれに対して電圧の変化を現在の電圧ベクトルに足すことにより、各ノードの新しい電圧値を与える新しい電圧ベクトルを求めるステップと、
次のタイムステップでＤＡＥの新しい系を解くために１タイムステップ進め、続くタイムステップに対して当該解くステップを繰り返すことにより、集積回路の過渡的な振る舞いをモデル化するステップとを含むことを特徴とする請求項２の方法。
【請求項４】
Ｍは並列処理に適した行列であり、Ｎは結合行列であることを特徴とする請求項１の方法。
【請求項５】
行列Ｉ＋Ｍ−１Ｎを解くステップは、
Ｕを上三角行列、Ｌを下三角行列として、Ｍ−１をＵ−１Ｌ−１に置き換えることにより項Ｕ−１Ｌ−１Ｎを形成するステップと、
Ｕ−１Ｌ−１Ｎを第１、第２、第３の三つの部分にスプリットするステップとを含み、
前記第１の部分は行列の相対的に大きい値であるエントリを含み、前記第２の部分と前記第３の部分は、計算過程で相対的に小さな値であるエントリを含むことを特徴とする請求項１の方法。
【請求項６】
前記第１の部分は行列Ｅであり、前記第２の部分は行列（ＬＵ）−１Ｆ１であり、前記第３の部分は行列Ｆ２であり、Ｕ−１Ｌ−１Ｎ＝Ｅ＋（ＬＵ）−１Ｆ１＋Ｆ２を満たすことを特徴とする請求項５の方法。
【請求項７】
Ｆ１はＬ−１Ｎ＝Ｅ０＋Ｌ−１Ｆ１の計算過程で形成され、
Ｆ２は、Ｕ−１Ｌ−１Ｎ＝Ｕ−１（Ｅ０＋Ｌ−１Ｆ１）＝Ｅ＋Ｕ−１Ｆ２＋Ｕ−１Ｌ−１Ｆ１の計算過程で形成されることを特徴とする請求項６の方法。
【請求項８】
Ｍ−１が部分系（Ｉ＋Ｅ＋Ｕ−１Ｆ２＋Ｕ−１Ｌ−１Ｆ１）ｘ＝Ｕ−１Ｌ−１ｒにより与えられ、ｘは、回路ノードにおける以前の近似からの電圧の変化を含むベクトルを表すΔｖとして定義され、
2009520306.tif 000040
であり、ｖはノードの電圧ベクトルであり、Ｑ（ｖ）はノードの電荷であり、ｉ（ｖ）はノードの電流であり、ｕ０は回路の電源を表すことを特徴とする請求項５の方法。
【請求項９】
部分系行列（Ｉ＋Ｅ）｜Ｓを形成するために、前記部分系はプレコンディショナー行列Ｉ＋Ｅを射影することによってプレコンディションされ、このプレコンディションは、（（Ｉ＋Ｅ）｜Ｓ）−１（Ｉ＋Ｅ＋Ｕ−１Ｆ２＋Ｕ−１Ｌ−１Ｆ１）｜Ｓｘ｜Ｓ＝（（Ｉ＋Ｅ）｜Ｓ）−１（Ｕ−１Ｌ−１ｒ）｜Ｓを与えることを特徴とする請求項８の方法。
【請求項１０】
回路ヤコビ行列をＭ＋Ｎにスプリットするステップ、Ｍ−１をＵ−１Ｌ−１に置き換えるステップ、およびＵ−１Ｌ−１Ｎを三つの部分にスプリットするステップは、再帰的に部分系行列（Ｉ＋Ｅ）｜Ｓに対して実行されることを特徴とする請求項８の方法。
【請求項１１】
最終的な部分系のサイズが所定の閾値より小さくなるか、前記部分系が行列特性解析を通じて直接解法にとって効率的であると判定されたとき、プレコンディションおよび再帰は停止し、最終的な部分系を解くために直接法が利用されることを特徴とする請求項１０の方法。
【請求項１２】
回路ヤコビ行列は頂点セパレータの集合にしたがってリオーダーされることを特徴とする請求項１の方法。
【請求項１３】
回路ヤコビ行列は、頂点セパレータに対応する列をスプリットアウトすることにより、行列ＭとＮにスプリットされることを特徴とする請求項１２の方法。
【請求項１４】
行列Ｍは上行列Ｕと下行列Ｌに分解され、プレコンディションするステップは、前記下行列Ｌを用いて行列ＭとＮをプレコンディションし、プレコンディションされた回路方程式を形成することを特徴とする請求項１３の方法。
【請求項１５】
プレコンディションされた回路方程式は、Ｌ−１ＪΔｖ＝Ｌ−１（Ｍ＋Ｎ）Δｖ＝（Ｕ＋Ｌ−１Ｎ）Δｖ＝Ｌ−１ｒにより定義され、Δｖは回路ノードにおける以前の近似からの電圧の変化を含むベクトルを表し、
2009520306.tif 000041
であり、ｖはノードの電圧ベクトルであり、Ｑ（ｖ）はノードの電荷であり、ｉ（ｖ）はノードの電流であり、ｕ０は回路の電源を表すことを特徴とする請求項１４の方法。
【請求項１６】
プレコンディションされた回路方程式（Ｕ＋Ｌ−１Ｎ）Δｖ＝Ｌ−１ｒは、Ｌ−１ＮをＬ−１Ｎ＝Ｅ＋Ｌ−１Ｆを満たすような第１の部分Ｅと第２の部分Ｌ−１Ｆにスプリットすることにより解かれ、第１の部分ＥはＬ−１Ｎの計算の後、相対的に大きいエントリを含み、第２の部分Ｌ−１Ｆは相対的に小さいエントリを含むことを特徴とする請求項１５の方法。
【請求項１７】
Ｌ−１Ｎのスプリッティングは、エントリの値を比較することによりＬ−１Ｎの計算をする過程でなされることを特徴とする請求項１６の方法。
【請求項１８】
前記第１の部分と前記第２の部分を最初の方程式に代入した結果、（Ｕ＋Ｅ＋Ｌ−１Ｆ）Δｖ＝Ｌ−１ｒによって定義される解くべき線形系が得られることを特徴とする請求項１５の方法。
【請求項１９】
前記解くべき線形系はプレコンディショナー（Ｕ＋Ｅ）を用いてプレコンディションされ、（Ｕ＋Ｅ）−１（Ｕ＋Ｅ＋Ｌ−１Ｆ）Δｖ＝（Ｉ＋（Ｕ＋Ｅ）−１Ｌ−１Ｆ）Δｖ＝（Ｕ＋Ｅ）−１Ｌ−１ｒで定義されるプレコンディションされた線形系が形成されることを特徴とする請求項１８の方法。
【請求項２０】
回路ヤコビ行列をリオーダーするステップ、回路ヤコビ行列を行列ＭとＮにスプリットしてプレコンディションされた回路方程式を形成するステップ、プレコンディションされた回路方程式を形成するために下行列Ｌを利用するステップ、および、Ｌ−１Ｎを二つの部分にスプリットしプレコンディショナー（Ｕ＋Ｅ）で系をプレコンディションすることにより、プレコンディションされた回路方程式を解くステップを含む一連のオペレーションを再帰的に実行するステップと、
縮小された部分系のサイズが所定の閾値より小さくなるか、縮小された部分系が行列特性解析を通じて直接解法にとって効率的であると判定されたとき、再帰を停止するステップとをさらに含むことを特徴とする請求項１９の方法。
【請求項２１】
Ｋｒｙｌｏｖ部分空間反復法はプレコンディションされた線形系の部分系に適用され、この部分系は頂点セパレータに対応することを特徴とする請求項１９の方法。
【請求項２２】
前記直接・反復組み合わせ解法はＫｒｙｌｏｖ部分空間反復法を含むことを特徴とする請求項１の方法。
【請求項２３】
集積回路のすべてのアクティブなパーティションのグループは一緒に解かれ、すべてのアクティブなパーティションのグループは、各タイムステップで各非線形の反復過程で動的に成長しうることを特徴とする請求項１の方法。
【請求項２４】
ＤＡＥの系は集積回路の静的な回路オペレーションをモデル化することを特徴とする請求項１の方法。
【請求項２５】
ＤＡＥの系は集積回路の周期的または準周期的な安定状態のオペレーションをモデル化することを特徴とする請求項１の方法。
【請求項２６】
疎結合のパーティションに分割された回路ヤコビ行列を用いて、集積回路のオペレーションをモデル化する微分代数方程式（ＤＡＥ）の系を解くことをコンピュータシステムに実行させるプログラムインストラクションが記録されたマシーン読み取り可能な媒体であって、当該マシーン読み取り可能な媒体は、
複数のノードのそれぞれにおける電圧の値を含む電圧ベクトルと前記回路ヤコビ行列を前記パーティションにしたがってリオーダーするためのインストラクションと、
前記回路ヤコビ行列を二つの行列ＭとＮにスプリットするためのインストラクションと、
対角成分が１で非対角成分がゼロの恒等行列をＩとして、前記二つの行列をプレコンディションしてＩ＋Ｍ−１Ｎの形の行列をもつプレコンディションされた方程式を形成するためのインストラクションと、
直接・反復組み合わせ解法を用いて前記プレコンディションされた方程式におけるＩ＋Ｍ−１Ｎに対する解を求めるためのインストラクションとを含むことを特徴とするマシーン読み取り可能媒体。
【請求項２７】
Ｍは並列処理に適した行列であり、Ｎは結合行列であることを特徴とする請求項２６のマシーン読み取り可能媒体。
【請求項２８】
行列Ｉ＋Ｍ−１Ｎを解くためのインストラクションは、
Ｕを上三角行列、Ｌを下三角行列として、Ｍ−１をＵ−１Ｌ−１に置き換えることにより項Ｕ−１Ｌ−１Ｎを形成するためのインストラクションと、
Ｕ−１Ｌ−１Ｎを第１、第２、第３の三つの部分にスプリットするためのインストラクションとを含み、
前記第１の部分は行列の相対的に大きい値であるエントリを含み、前記第２の部分と前記第３の部分は、計算過程で相対的に小さな値であるエントリを含むことを特徴とする請求項２６のマシーン読み取り可能媒体。
【請求項２９】
前記第１の部分は行列Ｅであり、前記第２の部分は行列（ＬＵ）−１Ｆ１であり、前記第３の部分は行列Ｆ２であり、Ｕ−１Ｌ−１Ｎ＝Ｅ＋（ＬＵ）−１Ｆ１＋Ｆ２を満たすことを特徴とする請求項２８のマシーン読み取り可能媒体。
【請求項３０】
Ｆ１はＬ−１Ｎ＝Ｅ０＋Ｌ−１Ｆ１の計算過程で形成され、
Ｆ２は、Ｕ−１Ｌ−１Ｎ＝Ｕ−１（Ｅ０＋Ｌ−１Ｆ１）＝Ｅ＋Ｕ−１Ｆ２＋Ｕ−１Ｌ−１Ｆ１の計算過程で形成されることを特徴とする請求項２９のマシーン読み取り可能媒体。
【請求項３１】
Ｍ−１が部分系（Ｉ＋Ｅ＋Ｕ−１Ｆ２＋Ｕ−１Ｌ−１Ｆ１）ｘ＝Ｕ−１Ｌ−１ｒにより与えられ、ｘは、回路ノードにおける以前の近似からの電圧の変化を含むベクトルを表すΔｖとして定義され、
2009520306.tif 000042
であり、ｖはノードの電圧ベクトルであり、Ｑ（ｖ）はノードの電荷であり、ｉ（ｖ）はノードの電流であり、ｕ０は回路の電源を表すことを特徴とする請求項２９のマシーン読み取り可能媒体。
【請求項３２】
部分系行列（Ｉ＋Ｅ）｜Ｓを形成するために、前記部分系はプレコンディショナー行列Ｉ＋Ｅを射影することによってプレコンディションされ、このプレコンディションは、（（Ｉ＋Ｅ）｜Ｓ）−１（Ｉ＋Ｅ＋Ｕ−１Ｆ２＋Ｕ−１Ｌ−１Ｆ１）｜Ｓｘ｜Ｓ＝（（Ｉ＋Ｅ）｜Ｓ）−１（Ｕ−１Ｌ−１ｒ）｜Ｓを与えることを特徴とする請求項３１のマシーン読み取り可能媒体。
【請求項３３】
回路ヤコビ行列をＭ＋Ｎにスプリットするためのインストラクション、Ｍ−１をＵ−１Ｌ−１に置き換えるためのインストラクション、およびＵ−１Ｌ−１Ｎを三つの部分にスプリットするためのインストラクションを、再帰的に部分系行列（Ｉ＋Ｅ）｜Ｓに対して実行するためのインストラクションをさらに含むことを特徴とする請求項３１のマシーン読み取り可能媒体。
【請求項３４】
縮小された部分系のサイズ、あるいはその縮小された部分系が行列特性解析を通じて直接解法にとって効率的であるかどうかの判定にもとづいて、再帰を停止するかどうかを決定するためのインストラクションと、
再帰の停止の後、直接法によって前記縮小された部分系を解くためのインストラクションとをさらに含むことを特徴とする請求項３３のマシーン読み取り可能媒体。
【請求項３５】
回路ヤコビ行列は頂点セパレータの集合にしたがってリオーダーされることを特徴とする請求項２６のマシーン読み取り可能媒体。
【請求項３６】
回路ヤコビ行列は、頂点セパレータに対応する列をスプリットアウトすることにより、行列ＭとＮにスプリットされることを特徴とする請求項３５のマシーン読み取り可能媒体。
【請求項３７】
行列Ｍは上行列Ｕと下行列Ｌに分解するためのインストラクションをさらに含み、プレコンディションするためのインストラクションは、前記下行列Ｌを用いて行列ＭとＮをプレコンディションし、Ｌ−１ＪΔｖ＝Ｌ−１（Ｍ＋Ｎ）Δｖ＝（Ｕ＋Ｌ−１Ｎ）Δｖ＝Ｌ−１ｒにより定義されるプレコンディションされた回路方程式を形成するためのインストラクションを含み、
Δｖは回路ノードにおける以前の近似からの電圧の変化を含むベクトルを表し、
2009520306.tif 000043
であり、ｖはノードの電圧ベクトルであり、Ｑ（ｖ）はノードの電荷であり、ｉ（ｖ）はノードの電流であり、ｕ０は回路の電源を表すことを特徴とする請求項３６のマシーン読み取り可能媒体。
【請求項３８】
プレコンディションされた回路方程式（Ｕ＋Ｌ−１Ｎ）Δｖ＝Ｌ−１ｒは、Ｌ−１ＮをＬ−１Ｎ＝Ｅ＋Ｌ−１Ｆを満たすような第１の部分Ｅと第２の部分Ｌ−１Ｆにスプリットすることにより解くためのインストラクションをさらに含み、第１の部分ＥはＬ−１Ｎの計算の後、相対的に大きいエントリを含み、第２の部分Ｌ−１Ｆは相対的に小さいエントリを含むことを特徴とする請求項３７のマシーン読み取り可能媒体。
【請求項３９】
Ｌ−１Ｎのスプリッティングは、エントリの値を比較することによりＬ−１Ｎの計算をする過程でなされることを特徴とする請求項３８のマシーン読み取り可能媒体。
【請求項４０】
前記第１の部分と前記第２の部分を最初の方程式に代入した結果、（Ｕ＋Ｅ＋Ｌ−１Ｆ）Δｖ＝Ｌ−１ｒによって定義される解くべき線形系が得られることを特徴とする請求項３８のマシーン読み取り可能媒体。
【請求項４１】
前記解くべき線形系を（Ｕ＋Ｅ）を用いてプレコンディションし、（Ｕ＋Ｅ）−１（Ｕ＋Ｅ＋Ｌ−１Ｆ）Δｖ＝（Ｉ＋（Ｕ＋Ｅ）−１Ｌ−１Ｆ）Δｖ＝（Ｕ＋Ｅ）−１Ｌ−１ｒで定義されるプレコンディションされた線形系を形成するためのインストラクションをさらに含むことを特徴とする請求項４０のマシーン読み取り可能媒体。
【請求項４２】
回路ヤコビ行列をリオーダーするためのインストラクション、回路ヤコビ行列を行列ＭとＮにスプリットしてプレコンディションされた回路方程式を形成するためのインストラクション、プレコンディションされた回路方程式を形成するために下行列Ｌを利用するためのインストラクション、および、Ｌ−１Ｎを二つの部分にスプリットしプレコンディショナー（Ｕ＋Ｅ）で系をプレコンディションすることにより、プレコンディションされた回路方程式を解くためのインストラクションを再帰的に実行するためのインストラクションステップと、
縮小された部分系のサイズ、あるいはその縮小された部分系が行列特性解析を通じて直接解法にとって効率的であるかどうかの判定にもとづいて、再帰を停止するかどうかを決定するためのインストラクションと、
再帰の停止の後、直接法によって前記縮小された部分系を解くためのインストラクションとをさらに含むことを特徴とする請求項４１のマシーン読み取り可能媒体。
【請求項４３】
Ｋｒｙｌｏｖ部分空間反復法はプレコンディションされた線形系の部分系に適用され、この部分系は頂点セパレータに対応することを特徴とする請求項４１のマシーン読み取り可能媒体。
【請求項４４】
前記直接・反復組み合わせ解法はＫｒｙｌｏｖ部分空間反復法を含むことを特徴とする請求項２６のマシーン読み取り可能媒体。
【請求項４５】
集積回路のすべてのアクティブなパーティションのグループは一緒に解かれ、すべてのアクティブなパーティションのグループは、各タイムステップで各非線形の反復過程で動的に成長しうることを特徴とする請求項２６のマシーン読み取り可能媒体。
【請求項４６】
ＤＡＥの系は集積回路の静的な回路オペレーションをモデル化することを特徴とする請求項２６のマシーン読み取り可能媒体。
【請求項４７】
ＤＡＥの系は集積回路の周期的または準周期的な安定状態のオペレーションをモデル化することを特徴とする請求項２６のマシーン読み取り可能媒体。
【発明の詳細な説明】
【背景技術】
【０００１】
電子設計自動化（ＥＤＡ）では、集積回路設計者が自分の設計を探求し、検証するために利用することができるソフトウェアが開発されている。半導体産業がナノメートルの時代に入ったため、設計者はＥＤＡツールを使って、製造される回路にナノメートルレベルの影響があるかどうかを調べることが必要になってきている。設計を検証しなかったとしたら、製造される集積回路は高い確率で機能不全に陥り、再設計を余儀なくされ、費用がかかる。
【０００２】
回路シミュレーションは回路設計が動作するかどうかを検証する最も正確な方法である。図１は、初期設計が実現されたときから回路を製造する際に行われる典型的な変換を例示するフローダイアグラム５０である。ブロック５２は初期回路設計を示す。この設計は、一連の接続されたデバイスエレメント５４に分解される。当該設計における各デバイスエレメントは、半導体製造工場や製造者によって検証されてきた正確な解析的モデルでモデル化される。解析的なモデルで表されたデバイスエレメントを用いれば、波線で示すブロック５５で示される回路シミュレータが、電圧値や電流値をある期間にわたってシミュレートすることができる。回路シミュレータ５５は、データに対して回路シミュレーションの演算を行うようにプログラムされたコンピュータシステムを含む。回路シミュレーションにおいて、当該回路のすべてのノードの電圧値と電流値は、フローダイアグラム５０のブロック５６に表された微分代数方程式（ＤＡＥ）の系（システム）を解くことにより得られる。
【０００３】
ＤＡＥは有限差分法を用いて離散化することができ、ニュートン−ラフソン法のような非線形反復法を用いて、反復処理で方程式を解く。各反復において、非線形方程式が以前に得られた解の周りで線形化され、ブロック５８で示される線形化された方程式系が生成される。そして、この線形系を解く必要がある。
【０００４】
行列解法技術が多くの科学、工学の分野で広く用いられている。ＥＤＡにおいて、行列解法技術は、回路シミュレーションのような分野で線形方程式系を解く上で欠かすことができない。
【０００５】
回路方程式は次のような形を取る。
【０００６】
【数１】
2009520306.tif 000003
【０００７】
ここでｖは当該回路でシミュレートされるすべてのノードにおける電圧のベクトルである。Ｑ（ｖ）はこれらのノードにおける電荷である。ｉ（ｖ）はこれらのノードにおける電流である。ｕ０は当該回路の電源である。
【０００８】
上記の方程式を解く際、有限微分法を用いて最初に微分演算子を近似する。ここでは説明のため、後退オイラー法を説明するが、他の有限微分法が使われているならば、他の方法を同様に適用することができる。離散方程式系は次のような形になる。
【０００９】
【数２】
2009520306.tif 000004
【００１０】
ここで時間ステップｔｍは既知であると仮定し、時刻ｔｍ＋１における解を求める。非線形離散方程式を解くためにニュートン−ラフソン反復法を用いることができる。反復処理において、初期の推定値ｖ０が与えられる。その後、解ｖ１、ｖ２、…のシーケンスが得られ、図２のグラフに示されるような非線形方程式の解に収束する。各反復において、非線形方程式は、以前の反復で得られた既知の解ｖｎの周りに線形化される。そして、次の線形方程式系を解いてΔｖを得る。
【００１１】
【数３】
2009520306.tif 000005
【００１２】
ここで
【００１３】
【数３ａ】
2009520306.tif 000006
【００１４】
【数３ｂ】
2009520306.tif 000007
【００１５】
【数３ｃ】
2009520306.tif 000008
【００１６】
Δｖを求めた後、更新された解ｖｎ＋１＝ｖｎ＋Δｖを求める。このプロセスは次の条件が満たされるまで続く。
【００１７】
【数４】
2009520306.tif 000009
【００１８】
ここでｔｏｌはある小さな誤差許容値である。式４の条件が満たされると、解は図２に示すように収束したと考えられ、ｔｍ＋１における解は次のようになる。
【００１９】
【数５】
2009520306.tif 000010
【００２０】
従来のＳｐｉｃｅ回路シミュレーションにおいて、誤差チェックはすべてのノードに対して行われ、｜Δｖ｜は、ベクトルΔｖすなわち全ノードにおける解の変化の中で、すべてのエントリの内の最大絶対値である。式４が満たされると、手続きは１タイムステップ進み、フローチャート５０のブロック５６に戻り、次の時間に対応する次の非線形方程式系を解く。このようにして、モデル化された回路において過渡的な変動をシミュレートすることができる。
【００２１】
各時刻において生成された電圧値を解析し、回路シミュレーションが当該回路の期待された演算と合致しているかどうかどうかを判定する。たとえば、電圧値は回路設計の論理的な解析と比較される。別の例では、電圧値と電流値は、回路設計の電源およびタイミングの分析を実行するために使うことができる。演算６８において、シミュレーションが期待された演算に合致しないならば、設計が変更され、欠陥が正される。そして、手続きは新しい回路設計のもとで演算５２に戻る。しかし、電圧が期待された演算と比べて遜色がないならば、回路設計はプロトタイプとして製造される。すなわち、製造回路７０が生成される。
【００２２】
線形系を解くことの主要な部分は、回路ヤコブ行列を解くことであり、これは、ＬＵ分解ステップと従来の直接的アプローチにおける前向きおよび後ろ向き代入ステップとを含む。高速で正確な行列解法は、高速で正確な回路シミュレーションにおいてたいへん重要である。
【００２３】
行列解法技術は二つのカテゴリに分けることができる。一つのカテゴリは直接行列解法と呼ばれる。他方のカテゴリは反復行列解法と呼ばれる。
【００２４】
回路シミュレーションにおいて次の線形系を考える。
【００２５】
【数６】
2009520306.tif 000011
【００２６】
ここでＪは回路のヤコビ行列であり、ｒは残差ベクトルであり、Δｖはノード電圧解の更新ベクトルである。
【００２７】
直接行列解法は、行列Ｊを最初にＬＵ分解する
【００２８】
【数７】
2009520306.tif 000012
【００２９】
ここでＬは下三角行列であり、Ｕは上三角行列である。その後、次式を解き、線形系に対する解Δｖを得る。
【００３０】
【数８】
2009520306.tif 000013
【００３１】
かつ
【００３２】
【数９】
2009520306.tif 000014
【００３３】
以下の説明において、便宜上、Δｖを表すためにｘを用いる。
【００３４】
【数１０】
2009520306.tif 000015
【００３５】
反復行列解法は、反復しながら解を得るか、解に近づこうとする。初期推定値ｘ０のもとで、反復法は次の処理で解ｘに近づく。
【００３６】
【数１１】
2009520306.tif 000016
【００３７】
ここでｎ≧０であり、Ｂは反復解法のためにプレコンディション（precondition）された行列である。反復解法が効率的であるためには、この処理は、相対的に小さい反復回数で許容できる正確さで解に近づくことができるものでなければならず、Ｂｘｎの計算は高速でなければならない。Ｋｒｙｌｏｖ部分空間反復法は十分な収束特性をもち、合理的な速さで回路解に近づくことができる。プレコンディションの目的は、反復処理を高速化するために、行列Ｂを恒等行列に近づけることである。
【００３８】
回路シミュレーションにおいて、線形系のサイズは非常に大きくなり、行列の要素のパターンが疎（スパース）になることがある。このため、標準的な回路シミュレータは、行列Ａの疎（スパース）行列表現を採用する。行列の疎（スパース）表現では、非零である行列の要素だけを保存するのが一般的である。疎（スパース）行列の非零パターンはＬＵ分解の効率化のために重要である。これは、オペレーションはＬおよびＵに非零のエントリを生成することがある一方、同じ場所のＪのエントリはゼロであるからである。これらのエントリはフィル−イン（fill-in）と呼ばれる。フィルインの数を減らすため、線形系をリオーダー（reorder）することができる。リオーダーリングの処理は線形系においてベクトルｘとｂのパーミュテーション（入れ替え）である。
【００３９】
【数１２】
2009520306.tif 000017
【００４０】
ここでＪｒはＪの行と列をパーミュテーション／リオーダーリングしたものであり、ｒ’はｒの行パーミュテーションであり、ｘ’はｘの列パーミュテーションである。さらなる例示の便宜のために、下付文字および上付文字は今後の説明では省略する。リオーダーリングの目的は、行列のＬＵ分解の過程で生成されるフィル−インの数を最小にすることであり、これによりシミュレーションが高速化する。
【００４１】
典型的な並列回路シミュレーションにおいて、このタイプのリオーダーリングは別の目的で実行される。回路行列を図３に示すいわゆるダブルボーダー（double-bordered）系にリオーダーしたいとする。
【００４２】
リオーダーされた行列ＪｒのＬＵ分解を実行するにあたって、あるブロックオペレーションは並列に実行することができる。すなわち、行列ブロックＡ１、Ａ２、…、Ａｍ−１のＬＵ分解は並列に実行できる。並列計算のボトルネックは、最後のブロックのＬＵ分解を実行することである。結合（カップリング）ブロック、すなわちＣブロックやＤブロックからの寄与があるため、最後のブロックはＬＵ分解処理の過程で密になりうることに留意する。
【００４３】
Ｙ．Ｓａａｄらによって提案されたドメインベース・マルチレベル再帰ブロック不完全ＬＵ法（ＢＩＬＵＴＭ）は我々の方法とは異なる。ＢＩＬＵＴＭは一般的な行列解法をとして提案された。その方法の一つのアプローチは、Ｋｒｙｌｏｖ部分空間法をＢＩＬＵＴＭのすべてのレベルで適用することであり、ＢＩＬＵＴＭ法の出版物で報告されている。その手続きのｌ番目のレベルで、行列のブロック因数分解が次のように近似的に計算される（上付文字はレベル番号に対応する）。
【００４４】
【数１２．１】
2009520306.tif 000018
【００４５】
ここで
2009520306.tif 000019
【００４６】
回路シミュレーションにとって、上記のアプローチを取ると、トップ回路行列レベルでＫｒｙｌｏｖ部分空間反復法を適用することになる。これは、特別な特徴をもつ回路シミュレーション行列の性質のゆえに非常に非効率的である。より具体的には、回路シミュレーション行列は、非常に大きいが、疎であるかもしれない。Ｋｒｙｌｏｖ部分空間法をトップレベルで適用することは、Ｋｒｙｌｏｖ部分空間ベクトルがトップレベル行列と同じサイズであることを意味する。これらのベクトルは大きなコンピュータメモリを消費し、計算が遅くなる。ＢＩＬＵＭの別の修正した使い方は、ＢＩＬＵＴＭの特定レベル、すなわち特定の縮小（reduced）ブロックでＫｒｙｌｏｖ部分空間法を適用することである。しかし、トップレベル行列を縮小ブロックに減らしたとしても、縮小ブロックは近似技術を適用しなければ、非常に密になることがある。
【００４７】
Ｓａａｄらによるアプローチにおいて、不完全ＬＵ分解（ＩＬＵ）は保存される一方、近似誤差は破棄される。これは、ＩＬＵ分解をプレコンディション行列として用いるには良い。よって、各反復において、元の行列とベクトルの積、およびプレコンディション行列とベクトルの積を計算することが必要である。
【発明の開示】
【発明が解決しようとする課題】
【００４８】
上述の理由のため、直接・反復ハイブリッドアプローチが、並列計算を用いたときのＬＵ分解によってもたらされるボトルネックの問題を解決するために必要である。特に回路シミュレータ５５（図１）が、複数のプロセッサコアまたはプロセッサユニットをもつコンピュータシステムであり、ロジックに応答してコンピュータシステムに回路シミュレーションの機能を実行させるものである場合、新しいハイブリッドアプローチが必要である。このロジックはソフトウェア、ハードウェア、もしくはソフトウェアとハードウェアの組み合わせによって実装することができる。
【課題を解決するための手段】
【００４９】
大まかに言えば、本発明は、並列マルチレート回路シミュレータを提供し、直接法と反復法の両方の利点を組み合わせ、シミュレーション過程で並列性を増やすことで上記のニーズに応えるものである。
【００５０】
本発明は、プロセス、装置、システム、デバイス、方法を含め、数多くのやり方で実装することができることが理解されよう。本発明のいくつかの実施の形態を以下で説明する。
【００５１】
ある実施の形態では、回路シミュレーションにおいて並列方程式を解くためのコンピュータで実装された方法が記述される。この方法は、回路ヤコビ行列を疎結合のパーティションに分割し、このパーティションにしたがって電圧ベクトルと当該行列をリオーダーリングし、当該ヤコビ行列を二つの行列ＭとＮに分ける。ここでＭは並列処理に適した行列であり、Ｎは結合行列である。ＭとＮはその後、Ｍ−１Ｊｘ＝（Ｉ＋Ｍ−１Ｎ）ｘ＝Ｍ−１ｒとなるようにプレコンディションされ、ヤコビ行列Ｊが反復解法を使って解かれる。
【００５２】
本発明の利点は、これから述べる詳細な説明において、図面を用いながら発明の原理を例示することによって明らかになる。
【発明を実施するための最良の形態】
【００５３】
本発明は、これから述べる詳細な説明において、同じような構成要素には同様の符号を付した図面を参照することにより、容易に理解される。
【００５４】
これから述べる説明において、本発明を完全に理解するために、数多くの具体的な詳細が説明される。しかし、本発明はそのような具体的な詳細がなくても実施できることは当業者にとって明らかである。別の例では、周知の処理操作や実装の詳細については詳細には述べられていない。いたずらに発明をわかりにくくすることを避けるためである。
【００５５】
プロトタイプ回路の生成に先立ち、当該回路をモデル化するための公知の技術を用いて微分代数方程式（ＤＡＥ）の系を生成してもよい。ＤＡＥを数値評価に適したものにするために離散化する。離散化の方法としていろいろな方法が知られており、ＤＡＥの系を離散化するためにいずれの方法を用いてもよい。ＤＡＥが非線形である場合、たとえば回路がトランジスタのような非線形の構成要素を用いているなどの場合、線形系になるように線形化し、ここで説明する並列行列法を用いて解くことができるようにしてもよい。線形系に対する解は、回路の各ノードにおいて、電圧、電流および／または電荷を与えるものであるか、電圧、電流および／または電荷に対する変化を与えるものであってもよい。最終的な解を解析することで、物理試験のためにプロトタイプを生成するのに先立って、モデル化された回路が期待通りに振る舞うかどうかを確かめることができる。
【００５６】
［パーティショニングおよびスプリッティング］
【００５７】
図４は、回路シミュレーションのために並列方程式を高速に解くための実施の形態に係る方法を提示するフローチャート１００である。ここで説明する本アプローチは、節点解析、修正節点解析、分散エレメントを用いた方程式などを含む異なる定式化の回路方程式に適用することができる。この手続きは開始ブロック１０２で始まり、オペレーション１０４に進み、行列がグラフで表現される。最初、回路ヤコビ行列は疎結合ブロックに分割される。行列Ｊは、一般的な並列マルチレベルグラフ分割（パーティショニング）アルゴリズムを適用することができるようにグラフに変換される。この変換において、まず対称行列が以下のように得られる。
【００５８】
【数１３】
2009520306.tif 000020
【００５９】
ここでＪ’はＪの転置行列であり、１≦ｉ，ｊ≦Ｊｓｉｚｅである。ＪｓｉｚｅはＪのサイズである。
【００６０】
Ｊｓは次のようにしてグラフに変換される。エッジが下三角または上三角部分の各エントリに対して生成される。行ｉ、列ｊのエントリはグラフにおいてはｉ番目の頂点とｊ番目の頂点を結ぶエッジとして表される。各頂点は同じインデックスの行や列を表す。その後、並列グラフ分割アルゴリズムを適用して、グラフを負荷分散したコンポーネントに分割し、そのコンポーネントを分離するためにカットするのに必要なエッジの数を最小化する。適切な並列グラフ分割アルゴリズムが知られており、すぐに利用可能である。分割の後、すべての結合エントリはスプリットされ、グラフのエッジは、最後のブロックにリオーダーされるというよりはむしろ、第２の行列にカットされることになる。
【００６１】
このプロセスは図５および図６の例で説明される。図５において、Ｊｅｘａｍｐｌｅは行列の例であり、ＪＳｅｘａｍｐｌｅはＪｅｘａｍｐｌｅを対称行列に変換したものである。ＪＳｅｘａｍｐｌｅは図６に示すようにグラフ１５０に変換することができる。図４に戻り、手続きはオペレーション１０６および１０７に続き、グラフ分割アルゴリズムがオペレーション１０４で生成されたグラフに適用され、オペレーション１０８に示すように負荷分散されたパーティションが生成される。図５の例で言えば、結果的に得られるグラフ１５０は図６に示される。グラフ１５０は二つのパーティション１５２と１５４に分割され、ただ一つのエッジ１５６だけがカットのために必要である。
【００６２】
次に、オペレーション１１０において、行列Ｊｅｘａｍｐｌｅは、対応するグラフパーティションにしたがって、図７に示されるようにリオーダーされる。この例では、第２パーティション（２，３）は第１パーティション（１，４）のノードの後、リオーダーされる。パーティショニング（分割）の後、行列はオペレーション１１２において二つのマトリックスにスプリットされる。説明の目的のため、スプリッティングを次のように標記する。
【００６３】
【数１４】
2009520306.tif 000021
【００６４】
ここでＭは並列計算のために適した行列であり、Ｎは結合行列である。たとえば、行列をスプリットする方法として、図８Ａおよび８Ｂに示すような二通りの方法がある。
【００６５】
オペレーション１１４において、線形系をプレコンディションし、反復法を適用して次のように行列Ｊを解いてもよい。
【００６６】
【数１５】
2009520306.tif 000022
【００６７】
ここでｘ≡Δｖであり、Ｉは対角成分が１で非対角成分が０である恒等行列である。
【００６８】
プレコンディショニングの目的は、プレコンディションされた行列Ｉ＋Ｍ−１Ｎをできる限り恒等行列に近づけることである。このようにして、反復法が近似解に収束するまでの反復ステップが少なくなる。従来のＦａｓｔＭＯＳ回路シミュレータは、パーティショニングのヒューリスティクスをたくさん用いて、弱い結合（カップリング）すなわちＭＯＳＦＥＴデバイスのゲートにおける容量結合を生成する。こういったヒューリスティクスは依然としてここでも有用である。しかし、ＦａｓｔＭＯＳ回路シミュレータで使われる従来の緩和タイプの方法よりもずっと先進的な反復法を適用することができる。ここで説明する先進的な反復法は強い結合を扱うこともできる。こういった反復法の収束特性は、行列Ｎに現れる強い結合に対しても依然として十分有効である。
【００６９】
並列計算の性能を改善するため、分割されリオーダーされた行列は、まず上非対角ブロックのエントリを見つけることでスプリットしてもよい。たとえば、行列がブロック分割の後、図９に示す構造をもつと仮定する。見つかった上非対角ブロックはＣ１１、Ｃ２１、Ｃ１２およびＣ２２である。これらは行列の二つの列、列８と列１１に属している。図１０Ａおよび図１０Ｂに示すように、列８と列１１のエントリは対角エントリを除いてすべて、結合行列Ｎに移動することができる。行列のスプリットの後、直接行列解法において行列操作Ｍ−１Ｎを適用してもよい。Ｍ−１Ｎには二つだけ列があることに留意されたい。行列Ｉ＋Ｍ−１Ｎは図１１に示された形を取る。
【００７０】
図４に戻り、オペレーション１１２で行列をスプリッティングし、オペレーション１１４でプレコンディショニングした後、独立した部分系がオペレーション１１６によって次のように与えられ、解かれる。
【００７１】
【数１６】
2009520306.tif 000023
【００７２】
オペレーション１１６において、この系は結合系に縮小され、解かれる。系のサイズに依存して、あるいは、行列の特性解析を通して、系が直接解法によって効率的に解くことができるかどうかに依存して、オペレーション１１７に示されるように、オペレーション１０４から１１６は、再帰的にさらに部分系を縮小するために繰り返される。この再帰は図１２および図１５を参照して以下でさらに説明する。オペレーション１１８において、縮小された結合系の解が方程式１５に後ろ向き代入され、当該手続きは完了ブロック１２０に示すように終了する。このように、後ろ向き代入を使って、すべての残りの解を得ることができる。結合インデックスの集合は次のように表記される。
【００７３】
【数１７】
2009520306.tif 000024
【００７４】
このアプローチは結合部分系のサイズが比較的小さいときに有効である。結合部分系のサイズが大きいときは効果が薄れる。結合部分系の行列は次のように密である。
【００７５】
【数１８】
2009520306.tif 000025
【００７６】
行列Ｓのサイズが大きいならば、その逆行列を求めるのは効率的ではない。
【００７７】
並列計算において、方程式１８の行列Ｓがボトルネックとなり、並列行列解法のスケーラビリティをより多くのＣＰＵコアあるいは多くの結合をもつ系に制限してしまう。理想的には、パーティションの数はＣＰＵコアの数に比例し、パーティションが多いと、Ｎ行列により多くのエントリが生じる。
【００７８】
［直接・反復ハイブリッド法］
【００７９】
図１２は、行列Ｉ＋Ｍ−１Ｎを解くための直接・反復ハイブリッドあるいは組み合わせ法の手続きの例を示すフローチャート２００である。この手続きは開始ブロック２０２で示すように始まり、オペレーション２０４に進み、Ｍ−１がＵ−１Ｌ−１に置き換えられる。次に、オペレーション２０６において、Ｕ−１Ｌ−１Ｎが次のように３つの部分にスプリットされる。
【００８０】
【数１９】
2009520306.tif 000026
【００８１】
最初の部分Ｅは、Ｕ−１Ｌ−１Ｎの計算の後、比較的大きな値になる行列のエントリを取る。第２の部分（ＬＵ）−１Ｆ１と第３の部分Ｕ−１Ｆ２は計算過程で、比較的小さな値を取る。
【００８２】
Ｆ１は次式の計算過程で、形成される。
【００８３】
【数２０】
2009520306.tif 000027
【００８４】
Ｆ２は次式の計算過程で、形成される。
【００８５】
【数２１】
2009520306.tif 000028
【００８６】
Ｆ１とＦ２を前向き消去Ｌ−１と後ろ向き代入Ｕ−１の過程でスプリットアウトする際、フィル−インエントリの相対的な大きさがチェックされる。もしそのエントリが元のベクトルにあるなら、それはＦ１またはＦ２にスプリットアウトされることはない。その代わり、異なる大きさの許容値を用いて結合部分系の行や他の行のエントリを制御してもよい。
【００８７】
上記のアプローチはＳａａｄらのアプローチとは異なる。このアプローチでは不完全なＬＵ分解の残余は保たれるが、Ｓａａｄのアプローチでは不完全なＬＵ分解の残余は破棄される。それは「背景」の節で前に説明した通りである。
【００８８】
上記のようにＵ−１Ｌ−１Ｎを３つの部分にスプリットすると、Ｍ−１は次を与える。
【００８９】
【数２２】
2009520306.tif 000029
【００９０】
この系には独立の結合部分系があることに留意すべきである。部分系の行および列のインデックスは結合インデックス集合ＩＣに属する。行列Ｉ＋Ｅの部分系への射影は次のように表記できる。
【００９１】
【数２３】
2009520306.tif 000030
【００９２】
オペレーション２０８において、この行列は、以下のように方程式２２の部分系へのプレコンディショナー（preconditioner）として使われる。
【００９３】
【数２４】
2009520306.tif 000031
【００９４】
最後に、オペレーション２１０において、Ｋｒｙｌｏｖ部分空間反復法を適用することができる。たとえば、一般最小残余（ＧＭＲＥＳ）法を用いてプレコンディションされた系を解くことができる。
【００９５】
ここで述べるＫｒｙｌｏｖ部分空間ベクトルの長さは縮小された部分系のサイズのことであるが、Ｓａａｄの方法ではＫｒｙｌｏｖ部分空間ベクトルの長さは元の行列のサイズのことであり、それは「背景」の節で前に説明した通りである。
【００９６】
縮小された部分系を解くために直接解法を用いることは都合がよいことを留意すべきである。それは、縮小された系が非常に疎であるときに意味がある。こうすることにより、その方法は純粋な直接法になる。非常に疎な行列に対して、直接解法は反復解法よりも高速になりうる。ある実施の形態では、この柔軟性を今から述べる直接・反復ハイブリッド法に有効に組み入れている。
【００９７】
部分系に対する解が得られた後、部分系にその解を後ろ向き代入することを通して全体の系に対する解が得られ、この手続きは完了ブロック２１２で示すように終了する。
【００９８】
縮小された部分系（Ｉ＋Ｅ）｜Ｓを解くために上記の直接・反復ハイブリッド法を再帰的に適用してもよいことに留意する。各再帰ステップにおいて、縮小された部分系のサイズはさらに小さくなる。ある実施の形態では、縮小された部分系のサイズが予め定めた閾値よりも小さくなったとき、あるいは、部分系が行列特性解析を通じて直接解法にとって効率的であると判定されたとき、再帰は停止する。この時点で直接法は、最終的な小さな部分系に適用され、この手続きは完了ブロック２１２で示されるように終了する。たとえば、ある実施の形態では、直説法は、部分系の行列が１００行以下になったときに適用される。別の実施の形態では、直説法は、部分系の行列のＬＵ分解が生成するフィル−インの数が小さくなったときに適用される。
【００９９】
［最小頂点セパレータを用いた行列のリオーダーリング］
【０１００】
別の例示的な実施の形態では、回路行列が頂点セパレータの集合にしたがってリオーダーされ、より効率的な解法技術が適用される。頂点セパレータの集合は、グラフを二つの連結されていないコンポーネントすなわちサブグラフに分けるために、当該グラフから取り除くことができる頂点の集合であり、Ｓに接続するエッジを伴う。このグラフ分割アルゴリズムは頂点セパレータの最小集合を見つけることを試みる。たとえば、頂点セパレータの最小集合は図６のグラフ１５０の例では一つの頂点であり、図１３に示すように頂点４である。頂点３を頂点セパレータとして選ぶこともできる。
【０１０１】
頂点は最後にある頂点セパレータの最小集合を用いてリオーダーすることができる。図６に示すグラフ１５０の例では、頂点４は既に最後の頂点であるからリオーダーリングの必要はない。もし頂点３が頂点セパレータとして選ばれたとしたら、頂点３が最後の頂点になるように頂点をリオーダーする。頂点セパレータ（頂点４）は最後に順序付けられているため、行列Ｊｅｘａｍｐｌｅは図５に示した元の形を保つことができる。
【０１０２】
行列Ｊｅｘａｍｐｌｅは、最後に順序付けられた頂点セパレータに対応する列をスプリットアウトする（分け出す）ことでスプリットすることができる。このような列をスプリットするためには、当該列の非対角エントリだけをスプリットアウトすればよい。スプリットアウトされた行列をＮと表記し、残りのメインの行列をＭと表記する。我々の例では、図１４に示すようにＪｅｘａｍｐｌｅ＝Ｍ＋Ｎとなる。
【０１０３】
ＭがＭ＝ＬＵのようなＬＵ分解をもつとする。下行列Ｌは、次のように回路シミュレーションにおいて線形系をプレコンディションするために使うことができる。
【０１０４】
【数２５】
2009520306.tif 000032
【０１０５】
図１５は、線形化されプレコンディションされた回路方程式（Ｕ＋Ｌ−１Ｎ）ｘ＝Ｌ−１ｒを解くためのフローチャート２５０を示す。手続きは開始ブロック２５２で示されるように始まり、オペレーション２５４に進み、Ｌ−１Ｎが次のように二つの部分にスプリットされる。
【０１０６】
【数２６】
2009520306.tif 000033
【０１０７】
第１の部分ＥはＬ−１Ｎの計算の後、相対的に大きな値になるエントリからなる。第２の部分Ｌ−１Ｆは相対的に小さいエントリからなる。このスプリッティングは、Ｌ−１Ｎの計算過程でエントリの値を比較することで形成される。回路シミュレーションにおいて解くべき線形系は次のようになる。
【０１０８】
【数２７】
2009520306.tif 000034
【０１０９】
次に、オペレーション２５６において、線形系は（Ｕ＋Ｅ）でプレコンディションされ、次のプレコンディションされた系を形成するようになる。
【０１１０】
【数２８】
2009520306.tif 000035
【０１１１】
ＥとＦは頂点セパレータに対応する最後の列においてだけ非ゼロの値をもつことに留意すべきである。したがって、行列Ｕ＋Ｅは、頂点セパレータに対応する最後のブロックを除いて上三角形式である。この線形系において、頂点セパレータに対応する部分系は独立である。
【０１１２】
最後に、オペレーション２５８において、ＧＭＲＥＳのようなＫｒｙｌｏｖ部分空間反復法を適用して部分系を解く。各反復において、次の行列−ベクトル積が計算される。
【０１１３】
【数２９】
2009520306.tif 000036
【０１１４】
ここで、｜ｓｕｂｓｙｓｔｅｍは、ベクトルの部分系への射影を示す。Ｆ・ｘを計算する際、ｘ｜ｓｕｂｓｙｓｔｅｍだけが寄与する。なぜなら、Ｆの非ゼロの列だけが部分系に対応するものであるからである。このことはまた、部分系が独立である理由を説明している。行列Ｌの逆行列を求める際、分離されたパーティションに対応する対角ブロックの逆行列を並列に求めることができる。行列−ベクトル積である（Ｕ＋Ｅ）−１・（Ｌ−１Ｆｘ）を計算する際、次の系が解かれる。
【０１１５】
【数３０】
2009520306.tif 000037
【０１１６】
ここで、Ｌ−１Ｆｘは既に計算されたベクトルである。Ｕの対角は１を含み、頂点セパレータに対応するＵの最後のブロックは１の対角成分だけをもつ。たとえば、行列が４つのパーティションをもつとする。この場合、Ｕ＋Ｅは図１６に示す構造をもつ。この系を解くために、頂点セパレータに対応する最後のブロックＩ＋Ｅ５が最初に解かれる。この部分系は独立であることに留意する。
【０１１７】
我々のＫｒｙｌｏｖ部分空間ベクトルの長さは、縮小された部分系のサイズのことであるが、Ｓａａｄの方法ではＫｒｙｌｏｖ部分空間ベクトルの長さは、「背景」の節で説明したように、元の行列のサイズのことである。
【０１１８】
縮小された部分系を解くために直接解法を用いることは都合がよいことを留意すべきである。それは、縮小された系が非常に疎であるときに意味がある。こうすることにより、その方法は純粋な直接法になる。非常に疎な行列に対して、直接解法は反復解法よりも高速になりうる。ここで述べる方法がこの柔軟性をもつことは有利であり、そのことは直接・反復ハイブリッド法のユニークな特徴である。
【０１１９】
部分系ｘｓｕｂｓｙｓｔｅｍに対する解が得られた後、後ろ向き代入を用いて、その系を縮小することができる。残りの系を解く際、個別のパーティションに対応するブロックは並列に解かれる。系に対する解は次のように得られる。
【０１２０】
【数３１】
2009520306.tif 000038
【０１２１】
図１５を参照して上述した、直接・反復ハイブリッド法を含む行列解法技術は、部分行列Ｉ＋Ｅ５の逆行列を求めるまたはこれを解くときに適用することができる。もし部分行列Ｉ＋Ｅ５のサイズがまだ大きすぎるなら、部分行列が疎であることを保ち、その計算に並列性をさらに求められるように、この方法を再帰的に適用することができる。各再帰ステップにおいて、縮小された部分系のサイズはさらに小さくなる。縮小された部分系のサイズが予め定めた閾値よりも小さくなったとき、あるいは、部分系が行列特性解析を通じて直接解法にとって効率的であると判定されたとき、再帰は停止する。この時点で直接法は、最終的な小さな部分系に適用され、この手続きは完了ブロック２６０で示されるように終了する。たとえば、ある実施の形態では、直説法は、部分系の行列が１００行以下になったときに適用される。別の実施の形態では、直説法は、部分系の行列のＬＵ分解が生成するフィル−インの数が小さくなったときに適用される。
【０１２２】
シミュレーション行列が巨大であり、依然として疎に群がっていることがありうるから、Ｋｒｙｌｏｖ部分空間法が部分系レベルで適用される。ここで述べる方法において、列が結合（カップリング）に対応している行列ブロックＮに、スプリッティング技術が適用され、行列ブロックＮは元の線形系レベルになる。それに加えて、上述の方法がＥ（Ｌ−１Ｎの一部）の計算において列ベースになる。結果として、この実装は、共有メモリ型並列コンピュータ上で並列化するのが簡単である。
【０１２３】
この方法の一つの利点は、Ｋｒｙｌｏｖ部分空間法をサイズが巨大であるトップレベルの回路行列Ｊには適用しないことである。その代わりにＫｒｙｌｏｖ部分空間法をサイズがずっと小さい結合部分系に適用する。しかし、可能な限り元の回路行列を疎に保つためにスプリッティングがＬ−１Ｎの計算に適用される。これは、トップレベルの回路線形系のサイズをもつＫｒｙｌｏｖ部分空間ベクトルを記憶し、直交化する問題を解く。したがって、今述べているアプローチは巨大な回路行列を解くための直接法および反復法の両方の利点を活用することができる。
【０１２４】
この数値的方法は、以下で述べる柔軟なデータ構造と相まって、行列解法が回路シミュレーションにおいて効率的な過渡シミュレーションを提供することを保証する。過渡シミュレーションでは、異なる非ゼロの値をもつ行列を何度も解く必要がある。本アプローチは行列スプリッティングを動的に扱う点で柔軟性を有する。以下述べるように、本アプローチはまた、静的、周期的、準周期的シミュレーションにとって好適である。
【０１２５】
［スプリッティングのための動的なデータ構造］
【０１２６】
行列Ｅ、Ｆ１およびＦ２は非常に動的である。具体的には、シミュレーション過程のあらゆるステップで変わるかもしれない。効率的で柔軟なデータストレージを実現するために、ハッシュテーブルを用いてもよい。行列の列ベクトルの各エントリに対して、元の列ベクトルにおける行インデックス、コンパクトベクトル（非ゼロ）におけるインデックス、およびその値を格納してもよい。ある実施の形態では、たとえば、行列の各列を格納するためにＣ＋＋のＳＴＬライブラリのＭａｐまたはＶｅｃｔｏｒを用いる。
【０１２７】
ハッシュテーブルはまた、全体系において部分系に属するインデックスを格納するために用いてもよい。Ｖｅｃｔｏｒは結合インデックス集合ＩＣを格納するために使ってもよい。ハッシュテーブルを用いることの一つの利点は、全体系のサイズがどれだけ巨大であったとしても、格納する必要のあるインデックスの数は制限されることである。さらに、インデックスがこの集合に属するかどうかは線形時間で決定することができる。
【０１２８】
別の実施の形態において、バイナリツリーを動的なデータ構造のために用いる。バイナリサーチツリーを用いることにより、インデックスは常に順序づけが保たれ、ツリーをある順序で容易に行き来することができるようになる。このデータ構造を用いることにより、ツリーに挿入することがｌｏｇ（ｎ）の時間で行える。ここでｎはツリーのサイズである。ｎがそれほど大きくないなら、このＣＰＵ時間は受け入れ可能である。同様に、ツリーにおけるエントリやインデックスの場所もｌｏｇ（ｎ）の時間で見つけることができる。
【０１２９】
こういった動的なデータ構造は、行列Ｅ、Ｆ１およびＦ２を構築するのに役立つ。我々はこれらの行列の構築の効率性を保つことができる。
【０１３０】
［並列マルチレートシミュレーション］
【０１３１】
集積回路はその機能においてマルチレートの振る舞いをする。マルチレートの振る舞いをするのは、集積回路がオペレーションの過程で高周波と低周波の両方を取り扱うという事実による。
【０１３２】
回路シミュレーションにおいて、高周波信号は解くにはより多くのタイムステップが必要であるが、低周波信号はより少ないタイムステップで解くことができる。さらに、デジタルライクな信号に対しては、エッジは解くためにより多くのタイムステップが必要であるが、フラットな領域は解くのにほとんどタイムステップを必要としない。従来のＳｐｉｃｅタイプの回路シミュレーションでは、その特性は利用することができない。なぜなら全体回路は分割されていないし、回路方程式を離散化するにあたり、グローバルタイムステップが適用されるからである。正確さを保証するためには、最小タイムステップをシミュレーションで適用しなければならない。シミュレーションがより多くのタイムステップを取れば取るほど、速度が遅くなる。
【０１３３】
従来のＳｐｉｃｅタイプの回路シミュレーションを高速化するために、ＦａｓｔＭｏｓ回路シミュレータは全体の回路を異なるパーティションに分割し、マルチレートの振る舞いを利用するために異なるパーティションに異なるタイムステップを割り当てる。あるパーティションにおける信号が高周波であるなら、シミュレーションにおいて小さなタイムステップを取る。別のパーティションにおける信号が低周波であるなら、シミュレーションにおいて大きなタイムステップを取る。このマルチレートシミュレーション法を実装する際、ＦａｓｔＭＯＳシミュレータはイベント駆動型のスキームを採用する。イベント駆動型スキームはイベント伝搬にしたがって、解くべきパーティションすなわち部分回路を決定する。もしあるパーティションでイベントが起こるなら、そのパーティションは一つの反復または一つのタイムステップの間、シミュレートされる。
【０１３４】
イベント駆動型スキームはシミュレートすべきパーティションをシミュレートするだけである。そのため、不必要な計算を省くことができ、従来のＳｐｉｃｅシミュレータに比べてかなりスピードアップを図ることできる。しかし、深刻な欠点はイベントフローがシリアルである、すなわち、あるイベントの発生が以前のイベントに依存することである。この分野における公知の試みは、波形緩和および楽観的スケジューリングである。これらの技術は一般的な回路シミュレーションに対して積極的過ぎて成功しない。
【０１３５】
あるイベントのシミュレーションを並列化することは、関わる計算量がほんのわずかであるため、解決策とはならない。あるイベントをシミュレーションするための仕事量は基本的には、あるタイムステップの間、あるパーティションをシミュレーションする際に関わる仕事量である。一般的に、これらのパーティションはＦａｓｔＭＯＳシミュレーションではたいへん小さい。たとえば、一つのパーティションで二つのトランジスタがあるだけである。このような小さな量のシミュレーションの仕事を並列化しても、期待に沿うような並列効率性を達成することはできない。
【０１３６】
ここで例示の目的で説明した並列マルチレートシミュレーションスキームにおいて、イベント駆動型スキームは使われていない。その代わり、すべてのアクティブなパーティションが一緒に解かれる。このアクティブなパーティションのグループは、各タイムステップにおいて各非線形反復過程で動的に成長することがある。
【０１３７】
周期的に訂正ステップが実行される。この訂正ステップにおいて、複数のタイムステップが正確さを向上させて一緒に解かれる。上述の並列行列ソルバー（solver）を、多くのパーティションと数多くのステップを必要とする結合系を切り離す（decouple）ために使ってもよい。もし強引（brute-force）な直接法または反復法で系を解くなら、計算は非常に遅くなるか非効率的になる。シューティング−ニュートン（shooting-Newton）ＲＦシミュレーションアルゴリズムにおいて使われるものに似た反復法は、あまりに非効率であろう。そのような強引な反復法において、対角および部分対角ブロックをプレコンディショナーとしてもつ線形系を解くためにＫｒｙｌｏｖ部分空間法を用いるであろう。この場合、結合（カップリング）はプレコンディショニング解法において考慮されない。このように強引な方法は大きな回路をシミュレートするには遅く、不適当である。
【０１３８】
［いろいろなタイプの回路シミュレーションにおけるアプリケーション］
【０１３９】
過渡的な回路シミュレーションに加えて、並列行列ソルバーを静的な回路シミュレーション（ＤＣ回路シミュレーション、周期的あるいは準周期的な定常状態ＲＦシミュレーションと呼ばれることもある）に適用することができる。ＤＣシミュレーションのある実施の形態（疑似過渡ＤＣシミュレーションと呼ばれる）において、電圧源と電流源を０値からＤＣ値まで上昇させる。それから、シミュレーションは、ＤＣ解である定常状態に到達するまで続けられる。ＤＣシミュレーションの別の実施の形態では、我々の方法を次の静的な回路方程式を解くために適用することができる。
【０１４０】
【数３１．５】
2009520306.tif 000039
【０１４１】
これは動的な部分がない方程式１である。
【０１４２】
並列行列ソルバーはＲＦ（radio frequency）回路シミュレーションに適用することもできる。有限差分法やシューティング−ニュートン法のようなＲＦ回路シミュレーション方法は、いくつかのタイムインターバルにおいて周期的または準周期的解を求めるために解く。その目的のため、周期的または準周期的境界条件がない方程式１を解く。周期的定常状態を例に取ると、境界条件はｖ（ｔ１）＝ｖ（ｔ２）である。ここで［ｔ１，ｔ２］は周期的なタイムインターバルである。有限差分法またはシューティング−ニュートン法はこの境界条件を満たす解を求める。ここで、上述の行列ソルバーを線形化された線形系を解くために使ってもよい。
【０１４３】
［コンピュータ実装］
【０１４４】
上記の実施の形態を念頭において、本発明は、コンピュータシステムに記憶されたデータに関係するコンピュータ実装された様々なオペレーションを利用することができることが理解されよう。これらのオペレーションは物理的な量の物理的操作を要求するものである。通常、必ずしもそうとは限らないが、こういった量は、保存、転送、結合、比較、さもなければ操作が可能な電気的もしくは磁気的な信号の形を取る。さらに、実行される操作は生産、特定、決定、比較のような用語で参照されることがある。
【０１４５】
ここで記述した本発明の一部を形成するオペレーションはいずれも有用なマシーンオペレーションである。本発明はまた、これらのオペレーションを実行するためのデバイスまたは装置に関する。装置は要求された目的のために特別に構成することができる。あるいは、装置は、コンピュータに記憶されたコンピュータプログラムによって選択的に活性化したり、構成したりすることができる汎用のコンピュータであってもよい。特に、ここで述べた教示内容に沿って書かれたコンピュータプログラムをもつ様々な汎用マシーンを用いることができる。あるいは、要求されたオペレーションを実行するためにより特化した装置を構成するとより便利である。
【０１４６】
本発明はコンピュータ読み取り可能な媒体上のコンピュータ読み取り可能なコードとして具体化することができる。コンピュータ読み取り可能な媒体は、データを記憶する任意のデータストレージデバイスであり、コンピュータシステムによって読み取ることができる。コンピュータ読み取り可能な媒体はまた、コンピュータコードが具体化された電磁搬送波を含む。コンピュータ読み取り可能な媒体の例として、ハードドライブ、ネットワーク接続ストレージ（ＮＡＳ）、リードオンリーメモリ、ランダムアクセスメモリ、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、ＣＤ−ＲＷ、磁気テープ、その他のオプティカル／非オプティカルデータストレージデバイスがある。コンピュータ読み取り可能な媒体をネットワーク接続されたコンピュータシステム上に分散させ、コンピュータ読み取り可能なコードを分散方式で保存し、実行するようにもできる。
【０１４７】
本発明の実施の形態は、単一のコンピュータで処理することができ、また、複数のコンピュータまたは相互接続されたコンピュータコンポーネントを用いて処理することができる。ここで使われるコンピュータとは、プロセッサ、メモリ、ストレージをもつスタンドアロンのコンピュータ、または、ネットワーク端末にコンピュータリソースを提供する分散コンピューティングシステムを含む。ある分散コンピューティングシステムでは、コンピュータシステムのユーザは、現実には多くのユーザ間で共有されているコンポーネント部分にアクセスしている。したがって、ユーザはネットワーク上で仮想コンピュータにアクセスすることができ、ユーザには、一人のユーザに対してカスタマイズされた専用の単一のコンピュータとして見える。
【０１４８】
前述の発明は、明確に理解するために詳細に記述したが、添付の請求項の範囲内で変更や修正をすることができることは明らかである。したがって、本実施の形態は例示と考えるべきであり、限定的に捉えるべきものではない。本発明はここで説明した詳細な内容に限定されるものではなく、添付の請求項の範囲およびその均等の範囲で変更することができる。
【図面の簡単な説明】
【０１４９】
【図１】回路設計から集積回路を生成する際に実行される典型的な変換とオペレーションを例示するフローダイアグラムである。
【図２】解のシーケンスがどのようにして収束するかを例示するグラフである。
【図３】ダブルボーダーシステムの一例を示す図である。
【図４】回路シミュレーションにおいて並列に方程式を高速に解く方法を提示するフローチャートである。
【図５】行列の例と行列の例から形成される対称行列を示す図である。
【図６】図５の対称行列から形成されるグラフである。
【図７】図６に表されたグラフの部分にしたがってリオーダーされた図５の行列の例を示す。
【図８】図８Ａおよび図８Ｂは図７に示す行列をスプリットする方法の例を示す図である。
【図９】行列構造の一例を示す図である。
【図１０】図１０Ａおよび図１０Ｂは、図９に示す行列構造上で実行されたスプリッティングオペレーションの結果、得られる行列ＭとＮを示す図である。
【図１１】図１０Ａおよび図１０Ｂに示した行列から生成された行列Ｉ＋Ｍ−１Ｎを示す図である。
【図１２】行列Ｉ＋Ｍ−１Ｎを解くための直接・反復ハイブリッドあるいは組み合わせ法の手続き例を示すフローチャートである。
【図１３】図６のグラフ例をグラフから取り除かれた頂点の集合（頂点４）とともに示す図である。
【図１４】図１３で特定された頂点セパレータにしたがって二つの行列にスプリットされた図５の行列の例を示す図である。
【図１５】線形化されプレコンディションされた回路方程式（Ｉ＋Ｌ−１Ｎ）ｘ＝Ｌ−１ｒを解くためのフローチャートを示す図である。
【図１６】Ｕ＋Ｅの行列構造を示す図である。
【図１】
2009520306.tif 000044
【図２】
2009520306.tif 000045
【図３】
2009520306.tif 000046
【図４】
2009520306.tif 000047
【図５】
2009520306.tif 000048
【図６】
2009520306.tif 000049
【図７】
2009520306.tif 000050
【図８Ａ】
2009520306.tif 000051
【図８Ｂ】
2009520306.tif 000052
【図９】
2009520306.tif 000053
【図１０Ａ】
2009520306.tif 000054
【図１０Ｂ】
2009520306.tif 000055
【図１１】
2009520306.tif 000056
【図１２】
2009520306.tif 000057
【図１３】
2009520306.tif 000058
【図１４】
2009520306.tif 000059
【図１５】
2009520306.tif 000060
【図１６】
2009520306.tif 000061
【国際調査報告】
2009520306.tif 000062