JP2956800B2 - 連立一次方程式に関する計算装置 - Google Patents
連立一次方程式に関する計算装置Info
- Publication number
- JP2956800B2 JP2956800B2 JP3239260A JP23926091A JP2956800B2 JP 2956800 B2 JP2956800 B2 JP 2956800B2 JP 3239260 A JP3239260 A JP 3239260A JP 23926091 A JP23926091 A JP 23926091A JP 2956800 B2 JP2956800 B2 JP 2956800B2
- Authority
- JP
- Japan
- Prior art keywords
- matrix
- coefficient matrix
- diagonal
- preprocessing
- elements
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
Description
【０００１】
【産業上の利用分野】本発明は、連立一次方程式の数値
解を計算する装置に係わり、特に連立一次方程式の係数
行列に対する前処理行列を計算する装置であって、大規
模な数値シミュレーションを行うベクトル計算機、並列
処理方式の計算機、ワークステーションなどに関する。
解を計算する装置に係わり、特に連立一次方程式の係数
行列に対する前処理行列を計算する装置であって、大規
模な数値シミュレーションを行うベクトル計算機、並列
処理方式の計算機、ワークステーションなどに関する。
【０００２】
【従来の技術】従来の技術としては、村田健郎他著「ス
ーパーコンピュータ、科学技術計算への応用」（丸善）
に論じられている不完全三角分解を前処理に使用する連
立一次方程式の解析方式がある。この方法でもベクトル
計算機向き前処理方式の工夫が記述されているが、複数
のベクトル計算処理装置を持つコンピュータや超並列計
算機には適用が困難である。
ーパーコンピュータ、科学技術計算への応用」（丸善）
に論じられている不完全三角分解を前処理に使用する連
立一次方程式の解析方式がある。この方法でもベクトル
計算機向き前処理方式の工夫が記述されているが、複数
のベクトル計算処理装置を持つコンピュータや超並列計
算機には適用が困難である。
【０００３】行列Ａの三角分解はＡを下三角行列Ｌと上
三角行列Ｕの積ＬＵ（＝Ａ）に分解するものであり、有
限要素法で離散近似して得られる連立一次方程式では係
数行列は疎である（全行列要素に占める非ゼロ要素の割
合が低い）のに対し、ＬＵ分解された行列はそれよりも
非ゼロ要素が多くなる。通常この非ゼロ要素の発生をfi
ll−inと呼んでいるが、上記論文の不完全ＬＵ分解法で
は、このfill−in部分を無視して（ゼロに近似して）Ｌ
Ｕ分解を行う。従って、係数行列の全非ゼロ要素に対し
て不完全ＬＵ分解を行って得られるＬ，ＵはそれぞれＡ
の下三角部分、上三角部分と同じ構造を持つ。
三角行列Ｕの積ＬＵ（＝Ａ）に分解するものであり、有
限要素法で離散近似して得られる連立一次方程式では係
数行列は疎である（全行列要素に占める非ゼロ要素の割
合が低い）のに対し、ＬＵ分解された行列はそれよりも
非ゼロ要素が多くなる。通常この非ゼロ要素の発生をfi
ll−inと呼んでいるが、上記論文の不完全ＬＵ分解法で
は、このfill−in部分を無視して（ゼロに近似して）Ｌ
Ｕ分解を行う。従って、係数行列の全非ゼロ要素に対し
て不完全ＬＵ分解を行って得られるＬ，ＵはそれぞれＡ
の下三角部分、上三角部分と同じ構造を持つ。
【０００４】さらに上記の方式を発展させて、複数のベ
クトル計算処理装置をもつコンピュータへ適用する方式
がＶan der Ｖorst（1982）に発表されているが、不
完全三角分解部分は複数のベクトル計算処理装置への適
用が困難で、反復計算部分でもmx×myに分割した領域で
はベクトル長がmxに制限される。
クトル計算処理装置をもつコンピュータへ適用する方式
がＶan der Ｖorst（1982）に発表されているが、不
完全三角分解部分は複数のベクトル計算処理装置への適
用が困難で、反復計算部分でもmx×myに分割した領域で
はベクトル長がmxに制限される。
【０００５】また他の従来の技術としては、村田健郎他
著「大型数値シミュレーション」（岩波書店）に論じら
れている不完全三角分解を前処理に使用して共役勾配法
系の解法で連立一次方程式の解を反復計算する方法があ
る。この方法では不完全三角分解にＧustafsson流の補
正を加えているが、補正係数αは、風上係数ωと移流拡
散問題におけるセル・ペクレ数Ｐeを使用して、α＝0.9
5max(0,1−Ｐe／(4＋6ω)）としている。一般に連立一
次方程式計算の汎用ライブラリでは、風上係数ωやセル
・ペクレ数をプログラムに渡すのが困難であり、またこ
の方法で補正係数αを計算しただけでは、セル・ペクレ
数が10を超える問題では共役勾配法系の収束が遅くな
り、多くの場合収束しなくなる。
著「大型数値シミュレーション」（岩波書店）に論じら
れている不完全三角分解を前処理に使用して共役勾配法
系の解法で連立一次方程式の解を反復計算する方法があ
る。この方法では不完全三角分解にＧustafsson流の補
正を加えているが、補正係数αは、風上係数ωと移流拡
散問題におけるセル・ペクレ数Ｐeを使用して、α＝0.9
5max(0,1−Ｐe／(4＋6ω)）としている。一般に連立一
次方程式計算の汎用ライブラリでは、風上係数ωやセル
・ペクレ数をプログラムに渡すのが困難であり、またこ
の方法で補正係数αを計算しただけでは、セル・ペクレ
数が10を超える問題では共役勾配法系の収束が遅くな
り、多くの場合収束しなくなる。
【０００６】
【発明が解決しようとする課題】上記、不完全三角分解
法では、どのように工夫しても並列化の程度は次元数ｎ
より相当小さくなる。また、従来不完全三角分解を使用
しないで乗算行列を前処理行列とする方法も提案されて
いたが、いずれも不完全三角分解法に比較して連立一次
方程式の解の収束が遅く、係数行列の性質が悪くなる
（悪条件性を有する）と収束が不安定になるという欠点
を持っていた。
法では、どのように工夫しても並列化の程度は次元数ｎ
より相当小さくなる。また、従来不完全三角分解を使用
しないで乗算行列を前処理行列とする方法も提案されて
いたが、いずれも不完全三角分解法に比較して連立一次
方程式の解の収束が遅く、係数行列の性質が悪くなる
（悪条件性を有する）と収束が不安定になるという欠点
を持っていた。
【０００７】また共役勾配法系の計算手法を用いて連立
一次方程式の解を反復計算するとき、従来の不完全三角
分解法を前処理に使用したのでは、移流拡散方程式を差
分法で離散化した行列の性質が悪くなる（セル・ペクレ
数が大きい場合等）と収束が遅く、パラメータの設定方
法を少し誤ると収束しなくなるケースが多かった。
一次方程式の解を反復計算するとき、従来の不完全三角
分解法を前処理に使用したのでは、移流拡散方程式を差
分法で離散化した行列の性質が悪くなる（セル・ペクレ
数が大きい場合等）と収束が遅く、パラメータの設定方
法を少し誤ると収束しなくなるケースが多かった。
【０００８】本発明の目的は、複数台のベクトル演算器
を有する計算機や並列処理計算機によく適用できるよう
な連立一次方程式の前処理行列を計算する方式を提供す
ることにある。
を有する計算機や並列処理計算機によく適用できるよう
な連立一次方程式の前処理行列を計算する方式を提供す
ることにある。
【０００９】本発明の他の目的は、連立一次方程式の係
数行列の性質が悪条件であっても、その数値解の収束性
を向上させることにある。
数行列の性質が悪条件であっても、その数値解の収束性
を向上させることにある。
【００１０】
（１）有限要素で離散近似して得られる連立一次方程式
を共役勾配法系の計算手法で、高速でかつ安定的に、し
かも並列処理向き前処理を行うため下記のような技術的
手段を採用した。
を共役勾配法系の計算手法で、高速でかつ安定的に、し
かも並列処理向き前処理を行うため下記のような技術的
手段を採用した。
【００１１】まず、並列処理向きにするため、前処理は
行列Ａの不完全三角分解を使用せず、行列Ａの非対角非
ゼロ要素をｍ個の部分行列Ｅ1，Ｅ2，・・・，Ｅmに分
け、前処理行列を（Ｉ−Ｅi）,ｉ＝１，２，・・・，ｍ
の乗算で構成する。これにより、連立一次方程式の解析
手法の全部分で行列の次元数ｎの並列度を有することを
可能にする。
行列Ａの不完全三角分解を使用せず、行列Ａの非対角非
ゼロ要素をｍ個の部分行列Ｅ1，Ｅ2，・・・，Ｅmに分
け、前処理行列を（Ｉ−Ｅi）,ｉ＝１，２，・・・，ｍ
の乗算で構成する。これにより、連立一次方程式の解析
手法の全部分で行列の次元数ｎの並列度を有することを
可能にする。
【００１２】次に、共役勾配法系の計算手法で収束性を
向上させ、計算を高速化させるために、非対角の非ゼロ
要素を行ごとに列番号インデックスを付けて記憶手段に
記憶させ、絶対値の大きい順にソートして、行列Ａの部
分行列、Ｅ1，Ｅ2，Ｅ3，・・・，Ｅmを作るとき、Ｅ1，
Ｅ2，Ｅ3，Ｅ4と最初の方は行列の一行当り一要素で分
割し、要素の絶対値の値が小さくなるに従って行列の一
行当りに複数個の要素を含めて分割する方式を用いる。
これにより（Ｉ−Ｅi），ｉ＝１，２，・・・,ｍと分割
した各行列のノルムσi＝‖Ｅi‖を１より相当小さくす
ることができ、共役勾配法系の収束性を向上させ、連立
一次方程式の数値解を高速に求めることができる。
向上させ、計算を高速化させるために、非対角の非ゼロ
要素を行ごとに列番号インデックスを付けて記憶手段に
記憶させ、絶対値の大きい順にソートして、行列Ａの部
分行列、Ｅ1，Ｅ2，Ｅ3，・・・，Ｅmを作るとき、Ｅ1，
Ｅ2，Ｅ3，Ｅ4と最初の方は行列の一行当り一要素で分
割し、要素の絶対値の値が小さくなるに従って行列の一
行当りに複数個の要素を含めて分割する方式を用いる。
これにより（Ｉ−Ｅi），ｉ＝１，２，・・・,ｍと分割
した各行列のノルムσi＝‖Ｅi‖を１より相当小さくす
ることができ、共役勾配法系の収束性を向上させ、連立
一次方程式の数値解を高速に求めることができる。
【００１３】次に、離散近似した連立一次方程式の係数
行列が悪条件（この場合は係数行列Ａの対角優位性が大
きくくずれた場合）にも、共役勾配法系の解法の収束を
安定化させるため、通常行列Ａと右辺ベクトルを対角行
列でスケーリングするところを、対角要素と非対角要素
の絶対値の和とでスケーリングする方法を採用した。こ
れにより、行列Ａが悪条件でも、分割した前処理行列
（Ｉ−Ｅi），ｉ＝１，２，・・・，ｍのいずれにおい
ても、各分割行列のノルムσi＝‖Ｅi‖を１より小さく
することができ、連立一次方程式を安定に解析すること
を可能にする。
行列が悪条件（この場合は係数行列Ａの対角優位性が大
きくくずれた場合）にも、共役勾配法系の解法の収束を
安定化させるため、通常行列Ａと右辺ベクトルを対角行
列でスケーリングするところを、対角要素と非対角要素
の絶対値の和とでスケーリングする方法を採用した。こ
れにより、行列Ａが悪条件でも、分割した前処理行列
（Ｉ−Ｅi），ｉ＝１，２，・・・，ｍのいずれにおい
ても、各分割行列のノルムσi＝‖Ｅi‖を１より小さく
することができ、連立一次方程式を安定に解析すること
を可能にする。
【００１４】（２）差分法で離散近似して得られる連立
一次方程式を共役勾配法系の計算手法で、高速かつ安定
に、しかも並列処理向き前処理を行うため下記のような
技術的手段を採用した。
一次方程式を共役勾配法系の計算手法で、高速かつ安定
に、しかも並列処理向き前処理を行うため下記のような
技術的手段を採用した。
【００１５】まず、並列処理向きにするために、前処理
は行列Ａの不完全三角分解を使用せず、行列Ａの非対角
非ゼロ帯行列を一本の帯ずつに分け、係数行列右下の帯
よりＡ1，Ａ2，・・・，Ａmとするｍ個の部分行列を作
成し、前処理行列を、（Ｉ−Ａi），ｉ＝１，２，・・
・，ｍの乗算で構成する。このときｍの値は、２次元５
点差分法のとき４、３次元７点差分法のとき６となる。
さらに周期境界条件があるとさらに２又は４だけｍの値
は大きくなる。これにより、連立一次方程式の解析手法
の全部分で行列の次元数ｎの並列度を有することが可能
である。また非対角行列をｍ個の部分行列に分けること
で、共役勾配法系の収束性を向上させる。
は行列Ａの不完全三角分解を使用せず、行列Ａの非対角
非ゼロ帯行列を一本の帯ずつに分け、係数行列右下の帯
よりＡ1，Ａ2，・・・，Ａmとするｍ個の部分行列を作
成し、前処理行列を、（Ｉ−Ａi），ｉ＝１，２，・・
・，ｍの乗算で構成する。このときｍの値は、２次元５
点差分法のとき４、３次元７点差分法のとき６となる。
さらに周期境界条件があるとさらに２又は４だけｍの値
は大きくなる。これにより、連立一次方程式の解析手法
の全部分で行列の次元数ｎの並列度を有することが可能
である。また非対角行列をｍ個の部分行列に分けること
で、共役勾配法系の収束性を向上させる。
【００１６】次に、離散近似した連立一次方程式の係数
行列が悪条件（この場合は係数行列Ａの対角優位性が大
きく崩れた場合）でも、共役勾配法系の解法の収束を安
定化させるため、通常、行列Ａと右辺ベクトルを対角行
列でスケーリングするところを、対角要素と非対角要素
の絶対値の和でスケーリングする方法を採用した。これ
により、行列Ａが悪条件でも、分割した前処理行列（Ｉ
−Ａi）,ｉ＝１，２，・・・，ｍのいずれにおいても、
各分割行列のノルムσi＝‖Ａi‖を１より小さくするこ
とができ、連立一次方程式を安定に解析することを可能
にする。
行列が悪条件（この場合は係数行列Ａの対角優位性が大
きく崩れた場合）でも、共役勾配法系の解法の収束を安
定化させるため、通常、行列Ａと右辺ベクトルを対角行
列でスケーリングするところを、対角要素と非対角要素
の絶対値の和でスケーリングする方法を採用した。これ
により、行列Ａが悪条件でも、分割した前処理行列（Ｉ
−Ａi）,ｉ＝１，２，・・・，ｍのいずれにおいても、
各分割行列のノルムσi＝‖Ａi‖を１より小さくするこ
とができ、連立一次方程式を安定に解析することを可能
にする。
【００１７】（３）差分法で離散近似して得られる連立
一次方程式を、不完全三角分解付き共役勾配法系の計算
手法で、高速かつ安定に計算するため、不完全三角分解
に下記のような技術手段を採用した。
一次方程式を、不完全三角分解付き共役勾配法系の計算
手法で、高速かつ安定に計算するため、不完全三角分解
に下記のような技術手段を採用した。
【００１８】従来、不完全三角分解を行うとき、係数行
列Ａに対してそのまま不完全三角分解していたのに対
し、本発明では、対角要素と非対角要素の絶対値の和と
を合わせて作成した値を対角要素と見なして不完全三角
分解を行う方式を採用した。これにより係数行列の性質
の悪さにあまり関係なく、Ｇustafsson流の補正におけ
る補正係数αを一定値に固定することができ、収束の安
定性を保つことができた。
列Ａに対してそのまま不完全三角分解していたのに対
し、本発明では、対角要素と非対角要素の絶対値の和と
を合わせて作成した値を対角要素と見なして不完全三角
分解を行う方式を採用した。これにより係数行列の性質
の悪さにあまり関係なく、Ｇustafsson流の補正におけ
る補正係数αを一定値に固定することができ、収束の安
定性を保つことができた。
【００１９】また、上記不完全三角分解を前処理に使用
した場合、係数行列の性質が悪い（対角要素に対して非
対角要素の絶対値の和の値が大きい）場合でも収束の安
定性が非常に良いという利点と、複数台のベクトル演算
器での処理効率が悪いという欠点がある。一方、上記
（２）の（Ｉ−Ａi），ｉ＝１，２，・・・，ｍの乗算
で構成する行列を前処理に使用すると、不完全三角分解
法と逆の利点と欠点がある。このためどのような条件で
も高速でかつ安定な収束をさせる技術的手段として、使
用可能なベクトル演算器の台数と対角要素と非対角要素
の絶対値の和との比の二つをパラメータとして、二つの
前処理方式を選定する方式を採用した。
した場合、係数行列の性質が悪い（対角要素に対して非
対角要素の絶対値の和の値が大きい）場合でも収束の安
定性が非常に良いという利点と、複数台のベクトル演算
器での処理効率が悪いという欠点がある。一方、上記
（２）の（Ｉ−Ａi），ｉ＝１，２，・・・，ｍの乗算
で構成する行列を前処理に使用すると、不完全三角分解
法と逆の利点と欠点がある。このためどのような条件で
も高速でかつ安定な収束をさせる技術的手段として、使
用可能なベクトル演算器の台数と対角要素と非対角要素
の絶対値の和との比の二つをパラメータとして、二つの
前処理方式を選定する方式を採用した。
【００２０】
（１）共役勾配法系の解法の前処理として、不完全三角
解法ではベクトル計算機向きにはできても、超並列処理
向きにするのは困難である。上記（１）の発明では、前
処理行列を（Ｉ−Ｅi），ｉ＝１，２，・・・，ｍの乗
算で構成することで、複数台のベクトル演算機を有する
コンピュータにも一部の超並列計算機にも向く計算方式
とすることができる。
解法ではベクトル計算機向きにはできても、超並列処理
向きにするのは困難である。上記（１）の発明では、前
処理行列を（Ｉ−Ｅi），ｉ＝１，２，・・・，ｍの乗
算で構成することで、複数台のベクトル演算機を有する
コンピュータにも一部の超並列計算機にも向く計算方式
とすることができる。
【００２１】共役勾配法系の解法で、収束性を向上させ
る手段として、非対角の非ゼロ要素を行ごとに列番号イ
ンデックスを付けて記憶させ、絶対値の大きい順にソー
トし、行列Ａの部分行列、Ｅ1，Ｅ2，・・・，Ｅmを作
るとき最初の方は一要素ずつで分割し、要素の絶対値が
小さくなるに従って複数個の要素を含めて分割する技術
的手段を用いることにより、分割した各行列のノルムσ
i＝‖Ｅi‖，ｉ＝１，２，・・・，ｍをすべて１より相
当小さくし、共役勾配法系の収束性を向上させ、連立一
次方程式の解析を高速化することができる。
る手段として、非対角の非ゼロ要素を行ごとに列番号イ
ンデックスを付けて記憶させ、絶対値の大きい順にソー
トし、行列Ａの部分行列、Ｅ1，Ｅ2，・・・，Ｅmを作
るとき最初の方は一要素ずつで分割し、要素の絶対値が
小さくなるに従って複数個の要素を含めて分割する技術
的手段を用いることにより、分割した各行列のノルムσ
i＝‖Ｅi‖，ｉ＝１，２，・・・，ｍをすべて１より相
当小さくし、共役勾配法系の収束性を向上させ、連立一
次方程式の解析を高速化することができる。
【００２２】さらに、共役勾配法系の収束を安定化させ
る技術的手段として、対角要素と非対角要素の絶対値の
和とで行列Ａと右辺ベクトルｂをスケーリングする方式
を採用し、係数行列Ａの対角優位性が大きく崩れた場合
も、連立一次方程式の数値解が収束するようにできる。
る技術的手段として、対角要素と非対角要素の絶対値の
和とで行列Ａと右辺ベクトルｂをスケーリングする方式
を採用し、係数行列Ａの対角優位性が大きく崩れた場合
も、連立一次方程式の数値解が収束するようにできる。
【００２３】また非ゼロ要素を下三角と上三角に分けて
記憶する必要がないため、記憶容量は不完全三角分解法
より少なくて済む。
記憶する必要がないため、記憶容量は不完全三角分解法
より少なくて済む。
【００２４】（２）共役勾配法系の解法の前処理とし
て、不完全三角解法ではベクトル計算機向きにはできて
も、超並列処理向きにするのは困難である。上記（２）
の発明では、前処理行列を（Ｉ−Ａi），ｉ＝１，２，
・・・，ｍの乗算で構成することで、複数台のベクトル
演算器を有するコンピュータにも、超並列計算機にも向
く計算方式とすることができる。また前処理行列を（Ｉ
−Ａi），ｉ＝１，２，・・・，ｍの乗算で構成するこ
とで、従来行われていたＥ＝Ａ1＋Ａ2＋・・・＋Ａmと
して、Ｉ−Ｅ＋Ｅ2−Ｅ3＋Ｅ4・・・とするより、共役
勾配法系の収束性を向上させることができる。
て、不完全三角解法ではベクトル計算機向きにはできて
も、超並列処理向きにするのは困難である。上記（２）
の発明では、前処理行列を（Ｉ−Ａi），ｉ＝１，２，
・・・，ｍの乗算で構成することで、複数台のベクトル
演算器を有するコンピュータにも、超並列計算機にも向
く計算方式とすることができる。また前処理行列を（Ｉ
−Ａi），ｉ＝１，２，・・・，ｍの乗算で構成するこ
とで、従来行われていたＥ＝Ａ1＋Ａ2＋・・・＋Ａmと
して、Ｉ−Ｅ＋Ｅ2−Ｅ3＋Ｅ4・・・とするより、共役
勾配法系の収束性を向上させることができる。
【００２５】さらに共役勾配法系の収束を安定化させる
技術的手段として、対角要素と非対角要素の絶対値の和
とで行列Ａと右辺ベクトルｂをスケーリングする方式を
採用し、係数行列Ａの対角優位性が大きく崩れた場合
も、連立一次方程式の数値解が安定に収束する。
技術的手段として、対角要素と非対角要素の絶対値の和
とで行列Ａと右辺ベクトルｂをスケーリングする方式を
採用し、係数行列Ａの対角優位性が大きく崩れた場合
も、連立一次方程式の数値解が安定に収束する。
【００２６】（３）従来の係数行列Ａに対してそのまま
不完全三角分解する方法では対角優位性が大きく崩れる
と、前処理付き共役勾配法系では収束しなくなるという
欠点があった。上記（３）の発明のように、対角要素と
非対角要素の和とを合わせて作成した値を対角要素と見
なして不完全三角分解を行うと、不完全三角分解の途中
で分解後の対角要素が極端に小さくなるという状態も回
避でき、係数行列の優位性が大きく崩れても収束は安定
するという効果がある。さらにＧustafsson流の補正係
数αを0.95程度に固定しても収束の速さが大きく変化し
ないという利点を持つ。
不完全三角分解する方法では対角優位性が大きく崩れる
と、前処理付き共役勾配法系では収束しなくなるという
欠点があった。上記（３）の発明のように、対角要素と
非対角要素の和とを合わせて作成した値を対角要素と見
なして不完全三角分解を行うと、不完全三角分解の途中
で分解後の対角要素が極端に小さくなるという状態も回
避でき、係数行列の優位性が大きく崩れても収束は安定
するという効果がある。さらにＧustafsson流の補正係
数αを0.95程度に固定しても収束の速さが大きく変化し
ないという利点を持つ。
【００２７】また複数台のベクトル演算器を持つベクト
ル並列計算機において、共役勾配法系の解法で連立一次
方程式の解を反復計算するとき、上記不完全三角分解を
前処理にするのと上記（２）の発明により（Ｉ−Ａ
i），ｉ＝１，２，・・・，ｍの乗算で構成される行列
を前処理にするのとでは、それぞれ利点と欠点があり、
それはベクトル演算器の台数と対角要素と非対角要素の
絶対値の和との比の二つのパラメータに大きく左右され
る。そのためこの二つのパラメータにより、二つの前処
理方法のうちより良い方式を自動選定し、より高速な計
算を可能にする。
ル並列計算機において、共役勾配法系の解法で連立一次
方程式の解を反復計算するとき、上記不完全三角分解を
前処理にするのと上記（２）の発明により（Ｉ−Ａ
i），ｉ＝１，２，・・・，ｍの乗算で構成される行列
を前処理にするのとでは、それぞれ利点と欠点があり、
それはベクトル演算器の台数と対角要素と非対角要素の
絶対値の和との比の二つのパラメータに大きく左右され
る。そのためこの二つのパラメータにより、二つの前処
理方法のうちより良い方式を自動選定し、より高速な計
算を可能にする。
【００２８】
（１）実施例１ 図１は本実施例における共役勾配法系の計算手法に用い
る前処理行列Ｍの作成手順を示す図である。図２は前処
理行列Ｍを構成する行列の具体的分割方法を示す図であ
る。図３は前処理付き共役勾配法系の計算手順を示す図
である。図４は有限要素法による数値シミュレーション
の手順を示す図である。図５は有限要素分割を示す図で
ある。図６は図５に示す分割の結果作成される連立一次
方程式の係数行列の記憶方法を示す図である。
る前処理行列Ｍの作成手順を示す図である。図２は前処
理行列Ｍを構成する行列の具体的分割方法を示す図であ
る。図３は前処理付き共役勾配法系の計算手順を示す図
である。図４は有限要素法による数値シミュレーション
の手順を示す図である。図５は有限要素分割を示す図で
ある。図６は図５に示す分割の結果作成される連立一次
方程式の係数行列の記憶方法を示す図である。
【００２９】図１は図３の前処理行列を作成するステッ
プ９の具体的な作成手順を示した図である。図２は図１
の行列分割を行うステップ５の具体的方法を示したもの
である。図３は図４の連立一次方程式の数値解を計算す
るステップ16の具体例を示したものである。
プ９の具体的な作成手順を示した図である。図２は図１
の行列分割を行うステップ５の具体的方法を示したもの
である。図３は図４の連立一次方程式の数値解を計算す
るステップ16の具体例を示したものである。
【００３０】図１の１は有限要素法で作成される連立一
次方程式の係数行列と右辺ベクトルを記憶手段に記憶す
るステップである。係数行列Ａは対角行列Ｄと非対角行
列、波つきのＡ、とは分けて記憶される。また非対角行
列、波つきのＡ、は非ゼロ要素だけ記憶させるために、
非ゼロ要素について行ごとに列番号インデックスを格納
した列番号テーブルＬと対にして格納される。図１では
行列の次元数をＮとし、対角要素を除いた一行当りの最
大非ゼロ要素数をＮＤで示している。２は係数行列と右
辺ベクトルをスケーリングするため値を計算するステッ
プである。いま対角要素Ｄ(i)，ｉ＝１，・・・，Ｎは
すべて正としておく。いままではこのスケーリングに１
／Ｄ(i)，ｉ＝１，・・・，Ｎが使用されていたが、対
角優位性が大きく崩れた場合は解が収束しなくなる。こ
のため、ここでは各行の非対角要素の絶対値の和Ｕi
（ｉ＝１，・・・，Ｎ）と対角要素Ｄ(i)（ｉ＝１，・
・・，Ｎ）より計算する方式を採用する。パラメータα
は一般に4.0程度に固定しておけばよい。但しそれでも
収束しない場合にαの値を変更できるようにしておく。
３は係数行列と右辺ベクトルのスケーリングを行うステ
ップである。ステップ２で計算したωi（ｉ＝１，・・
・，Ｎ)を乗算することでスケーリングを行う。４は係
数行列要素をソートするステップである。ここではその
絶対値が大きい順にソートしている。なおこのスケーリ
ングとソーティングの順序はこの逆でもよい。５はソー
トした行列要素を使用してｍ個の部分行列Ｅ1，・・
・，Ｅmに分割するステップである。この部分について
は図２でさらに詳しく述べる。６は部分行列Ｅ1，・・
・，Ｅmを使用して前処理行列Ｍを作成するステップで
ある。上記各ステップは、記憶手段に格納されたデータ
を参照しながら上記の計算及び処理を実行する手段によ
って実現される。
次方程式の係数行列と右辺ベクトルを記憶手段に記憶す
るステップである。係数行列Ａは対角行列Ｄと非対角行
列、波つきのＡ、とは分けて記憶される。また非対角行
列、波つきのＡ、は非ゼロ要素だけ記憶させるために、
非ゼロ要素について行ごとに列番号インデックスを格納
した列番号テーブルＬと対にして格納される。図１では
行列の次元数をＮとし、対角要素を除いた一行当りの最
大非ゼロ要素数をＮＤで示している。２は係数行列と右
辺ベクトルをスケーリングするため値を計算するステッ
プである。いま対角要素Ｄ(i)，ｉ＝１，・・・，Ｎは
すべて正としておく。いままではこのスケーリングに１
／Ｄ(i)，ｉ＝１，・・・，Ｎが使用されていたが、対
角優位性が大きく崩れた場合は解が収束しなくなる。こ
のため、ここでは各行の非対角要素の絶対値の和Ｕi
（ｉ＝１，・・・，Ｎ）と対角要素Ｄ(i)（ｉ＝１，・
・・，Ｎ）より計算する方式を採用する。パラメータα
は一般に4.0程度に固定しておけばよい。但しそれでも
収束しない場合にαの値を変更できるようにしておく。
３は係数行列と右辺ベクトルのスケーリングを行うステ
ップである。ステップ２で計算したωi（ｉ＝１，・・
・，Ｎ)を乗算することでスケーリングを行う。４は係
数行列要素をソートするステップである。ここではその
絶対値が大きい順にソートしている。なおこのスケーリ
ングとソーティングの順序はこの逆でもよい。５はソー
トした行列要素を使用してｍ個の部分行列Ｅ1，・・
・，Ｅmに分割するステップである。この部分について
は図２でさらに詳しく述べる。６は部分行列Ｅ1，・・
・，Ｅmを使用して前処理行列Ｍを作成するステップで
ある。上記各ステップは、記憶手段に格納されたデータ
を参照しながら上記の計算及び処理を実行する手段によ
って実現される。
【００３１】図２は非対角非ゼロ要素行列、波つきの
Ａ、をｍ個の行列に分割する具体的な手法を示す図であ
る。７は一行当りの最大非ゼロ要素数ＮＤが大きい場合
の分割方法である。｜波つきのＡ（ｉ，ｊ）｜≦｜波つ
きのＡ（ｉ，ｊ＋１）｜（ｉ＝１，・・・，ｎ，ｊ＝
１，・・・，ＮＤ−１）のようにｊでソートされている
ため、ｊの値が小さいところは１行１要素で分割し、ｊ
が大きくなると少しずつ多くの要素を集めて部分行列と
している。８は一行当りの最大非ゼロ要素数ＮＤが小さ
い場合の分割方法を示したものである。
Ａ、をｍ個の行列に分割する具体的な手法を示す図であ
る。７は一行当りの最大非ゼロ要素数ＮＤが大きい場合
の分割方法である。｜波つきのＡ（ｉ，ｊ）｜≦｜波つ
きのＡ（ｉ，ｊ＋１）｜（ｉ＝１，・・・，ｎ，ｊ＝
１，・・・，ＮＤ−１）のようにｊでソートされている
ため、ｊの値が小さいところは１行１要素で分割し、ｊ
が大きくなると少しずつ多くの要素を集めて部分行列と
している。８は一行当りの最大非ゼロ要素数ＮＤが小さ
い場合の分割方法を示したものである。
【００３２】なお図２は分割方法の一例を示すものであ
るが、一般的には｜波つきのＡ（ｉ，ｊ）｜単独かある
いはいくつかの｜波つきのＡ（ｉ，ｊ）｜を加えたもの
がほぼ同じ程度の値になるように部分行列のグループ分
けをすればよく、これは計算機によって処理可能であ
る。
るが、一般的には｜波つきのＡ（ｉ，ｊ）｜単独かある
いはいくつかの｜波つきのＡ（ｉ，ｊ）｜を加えたもの
がほぼ同じ程度の値になるように部分行列のグループ分
けをすればよく、これは計算機によって処理可能であ
る。
【００３３】図３は前処理付きＣＧＰＭ法（共役勾配最
小多項式法、Ｖander ＶorstのＢi−ＣＧＳＴＡＢ法
（1990）と同様な計算手法）で連立一次方程式Ａ・ベク
トルｘ＝ベクトルｂの解であるベクトルｘを反復計算し
て求める計算手順である。９は入力された行列Ａから前
処理行列Ｍを作成するステップである。この具体的手順
は図１で示した通りである。このとき、行列Ａも右辺ベ
クトルｂもスケーリング処理される。10は反復計算のた
めの前準備をするステップである。ここでベクトルｘは
解の初期ベクトルで、分からない場合はオールゼロを入
れておけばよい。ベクトルｒは残差ベクトルである。ま
た（ベクトルｒ，ベクトルｒ）は内積計算を示し、以下
も同じ意味で使用する。11は反復計算を制御するステッ
プである。12は前処理付きＣＧＰＭ法の反復計算を実行
するステップである。計算解はベクトルｘに求まる。Ａ
とＭは行列、ベクトルｐ，ｑ，ｒ，ｒ0，ｅおよびｖは
次元数Ｎのベクトルを示し、それ以外の記号はスカラを
示す。ステップ12は公知技術であり、詳細説明を省略す
る。上記各ステップは、記憶手段に格納されたデータを
参照しながら計算および処理を実行する手段によって実
現される。
小多項式法、Ｖander ＶorstのＢi−ＣＧＳＴＡＢ法
（1990）と同様な計算手法）で連立一次方程式Ａ・ベク
トルｘ＝ベクトルｂの解であるベクトルｘを反復計算し
て求める計算手順である。９は入力された行列Ａから前
処理行列Ｍを作成するステップである。この具体的手順
は図１で示した通りである。このとき、行列Ａも右辺ベ
クトルｂもスケーリング処理される。10は反復計算のた
めの前準備をするステップである。ここでベクトルｘは
解の初期ベクトルで、分からない場合はオールゼロを入
れておけばよい。ベクトルｒは残差ベクトルである。ま
た（ベクトルｒ，ベクトルｒ）は内積計算を示し、以下
も同じ意味で使用する。11は反復計算を制御するステッ
プである。12は前処理付きＣＧＰＭ法の反復計算を実行
するステップである。計算解はベクトルｘに求まる。Ａ
とＭは行列、ベクトルｐ，ｑ，ｒ，ｒ0，ｅおよびｖは
次元数Ｎのベクトルを示し、それ以外の記号はスカラを
示す。ステップ12は公知技術であり、詳細説明を省略す
る。上記各ステップは、記憶手段に格納されたデータを
参照しながら計算および処理を実行する手段によって実
現される。
【００３４】図４は有限要素法による数値シミュレーシ
ョンの手順を示す図であり、単純化するため、線形定常
解析の場合の一例を示す。13は解析の前処理を行うステ
ップである。14は有限要素近似するための領域のメッシ
ュ分割を行うステップである。15は有限要素法によって
連立一次方程式を作成するステップを示す。16はその連
立一次方程式の数値解を計算するステップである。ステ
ップ16は図３に示す前処理付き共役勾配法系の解法で計
算される。
ョンの手順を示す図であり、単純化するため、線形定常
解析の場合の一例を示す。13は解析の前処理を行うステ
ップである。14は有限要素近似するための領域のメッシ
ュ分割を行うステップである。15は有限要素法によって
連立一次方程式を作成するステップを示す。16はその連
立一次方程式の数値解を計算するステップである。ステ
ップ16は図３に示す前処理付き共役勾配法系の解法で計
算される。
【００３５】図５は有限要素分割の一例を示す。18は４
節点４辺形要素及び一部３節点３辺形要素での分割を示
す。19は分割した要素の節点番号を示す。
節点４辺形要素及び一部３節点３辺形要素での分割を示
す。19は分割した要素の節点番号を示す。
【００３６】図６は、図５の有限要素分割で作成される
連立一次方程式の係数行列の非対角非ゼロ要素の記憶方
法の一例を示す。20は非ゼロ要素の列番号テーブルＬ
（Ｎ，ＮＤ）を示す。21は非ゼロ要素行列、波つきのＡ
（Ｎ，ＮＤ）、を示す。ｊ＝Ｌ（ｉ，ｋ）とするとき、
Ａ（ｉ，ｊ）は密行列表示のａi,jの要素を示す。22は
列番号インデックスで、この部分が空白の場合は対応す
る非ゼロ要素が存在しないことを示す。23は行列の非ゼ
ロ要素を示し、24は記憶領域は持つが値がゼロの要素を
示す。この形の非ゼロ要素行列、波つきのＡ（Ｎ，Ｎ
Ｄ）、を使用して図１に示す前処理行列Ｍを作成する。
連立一次方程式の係数行列の非対角非ゼロ要素の記憶方
法の一例を示す。20は非ゼロ要素の列番号テーブルＬ
（Ｎ，ＮＤ）を示す。21は非ゼロ要素行列、波つきのＡ
（Ｎ，ＮＤ）、を示す。ｊ＝Ｌ（ｉ，ｋ）とするとき、
Ａ（ｉ，ｊ）は密行列表示のａi,jの要素を示す。22は
列番号インデックスで、この部分が空白の場合は対応す
る非ゼロ要素が存在しないことを示す。23は行列の非ゼ
ロ要素を示し、24は記憶領域は持つが値がゼロの要素を
示す。この形の非ゼロ要素行列、波つきのＡ（Ｎ，Ｎ
Ｄ）、を使用して図１に示す前処理行列Ｍを作成する。
【００３７】本実施例において、前処理行列を不完全三
角分解でなく、（Ｉ−Ｅi），ｉ＝１，・・・，ｍの乗
算行列にしたのは、並列ベクトル計算機及び一部の超並
列計算機に良く適用できるようにするためである。不完
全三角分解を使用すると並列度は行列の次元数Ｎより大
幅に小さくなるが、本方式では常に次元数Ｎの並列度を
持つ。
角分解でなく、（Ｉ−Ｅi），ｉ＝１，・・・，ｍの乗
算行列にしたのは、並列ベクトル計算機及び一部の超並
列計算機に良く適用できるようにするためである。不完
全三角分解を使用すると並列度は行列の次元数Ｎより大
幅に小さくなるが、本方式では常に次元数Ｎの並列度を
持つ。
【００３８】非対角非ゼロ行列を絶対値の大きい順にソ
ートし、ｍ個の部分行列Ｅ1,・・・，Ｅmに分割し、前
処理行列を（Ｉ−Ｅi），ｉ＝１，２，・・・，ｍの乗
算で構成することにより、不完全三角分解の前処理と同
等か、それ以上の解の収束性を共役勾配法系の解法で得
ることができる。
ートし、ｍ個の部分行列Ｅ1,・・・，Ｅmに分割し、前
処理行列を（Ｉ−Ｅi），ｉ＝１，２，・・・，ｍの乗
算で構成することにより、不完全三角分解の前処理と同
等か、それ以上の解の収束性を共役勾配法系の解法で得
ることができる。
【００３９】また係数行列と右辺ベクトルのスケーリン
グに対角行列だけでなく非対角要素の行ごとの絶対値の
和を合わせて使用することにより、共役勾配法系の収束
を安定化させることができる。特に、対角優位性が大幅
に崩れた行列を係数とする場合に、従来収束しなかった
ものが収束するようになる。
グに対角行列だけでなく非対角要素の行ごとの絶対値の
和を合わせて使用することにより、共役勾配法系の収束
を安定化させることができる。特に、対角優位性が大幅
に崩れた行列を係数とする場合に、従来収束しなかった
ものが収束するようになる。
【００４０】また、前処理行列と係数行列がスケーリン
グした同じ値の要素を使用することで、主記憶容量を減
少させ、ワークステーションでも大次元行列が扱える。 （２）実施例２ 図７は本実施例における共役勾配法系の計算に用いる前
処理行列Ｍの作成手順を示す図である。図８は差分法に
よる数値シミュレーションの手順を示す図である。図９
は３次元７点差分法における非ゼロ帯行列の構成を示す
図である。
グした同じ値の要素を使用することで、主記憶容量を減
少させ、ワークステーションでも大次元行列が扱える。 （２）実施例２ 図７は本実施例における共役勾配法系の計算に用いる前
処理行列Ｍの作成手順を示す図である。図８は差分法に
よる数値シミュレーションの手順を示す図である。図９
は３次元７点差分法における非ゼロ帯行列の構成を示す
図である。
【００４１】図７は図３の前処理行列を作成するステッ
プ９の具体的な作成手順を示した図である。図３は図８
の連立一次方程式の数値解を計算するステップ43の具体
例を示した図となる。
プ９の具体的な作成手順を示した図である。図３は図８
の連立一次方程式の数値解を計算するステップ43の具体
例を示した図となる。
【００４２】図７の31は差分法で作成される連立一次方
程式の係数行列と右辺ベクトルを記憶手段に記憶するス
テップであり、係数行列Ａは対角行列Ｄと非対角行列、
波つきのＡ、とに分けて記憶される。ここでは３次元７
点差分法で離散化した場合の例が示されている。２次元
５点差分法では、非対角非ゼロ要素帯行列は波つきのＡ
（Ｎ，４）となる。また周期境界条件がある場合は非対
角非ゼロ要素帯行列の本数が２又は４増加する。ここで
次元数はＮとする。32は係数行列と右辺ベクトルをスケ
ーリングするための値を計算するステップである。いま
対角要素Ｄ（ｉ），ｉ＝１，・・・，Ｎはすべて正とし
ておく。従来はこのスケーリングに１／Ｄ（ｉ），ｉ＝
１，・・・，Ｎが使用されていたが、対角優位性が大き
く崩れた場合は解が収束しなくなる。このため、ここで
は各行の非対角要素の絶対値の和Ｕi,（ｉ＝１，・・
・，Ｎ）と対角要素Ｄ(i)(ｉ＝１，・・・，Ｎ）より計
算する方式を採用する。パラメータαは一般に4.0程度
に固定しておけばよい。但しそれでも収束しない場合に
αの値を変更できるようにしておく。33は係数行列と右
辺ベクトルのスケーリングを行うステップである。ステ
ップ32で計算したωiを乗算することでスケーリングを
行う。34は非対角行列要素をｍ個の部分行列Ａ1，Ａ2，
・・・，Ａmに分割するステップである。ここでは３次
元７点差分のため６個の部分行列Ａ1，Ａ2，・・・，Ａ
6に分割している。35は部分行列Ａ1，・・・，Ａ6を使
用して前処理行列Ｍを作成するステップである。上記各
ステップは、記憶手段に格納されたデータを参照しなが
ら上記の計算及び処理を実行する手段によって実現され
る。
程式の係数行列と右辺ベクトルを記憶手段に記憶するス
テップであり、係数行列Ａは対角行列Ｄと非対角行列、
波つきのＡ、とに分けて記憶される。ここでは３次元７
点差分法で離散化した場合の例が示されている。２次元
５点差分法では、非対角非ゼロ要素帯行列は波つきのＡ
（Ｎ，４）となる。また周期境界条件がある場合は非対
角非ゼロ要素帯行列の本数が２又は４増加する。ここで
次元数はＮとする。32は係数行列と右辺ベクトルをスケ
ーリングするための値を計算するステップである。いま
対角要素Ｄ（ｉ），ｉ＝１，・・・，Ｎはすべて正とし
ておく。従来はこのスケーリングに１／Ｄ（ｉ），ｉ＝
１，・・・，Ｎが使用されていたが、対角優位性が大き
く崩れた場合は解が収束しなくなる。このため、ここで
は各行の非対角要素の絶対値の和Ｕi,（ｉ＝１，・・
・，Ｎ）と対角要素Ｄ(i)(ｉ＝１，・・・，Ｎ）より計
算する方式を採用する。パラメータαは一般に4.0程度
に固定しておけばよい。但しそれでも収束しない場合に
αの値を変更できるようにしておく。33は係数行列と右
辺ベクトルのスケーリングを行うステップである。ステ
ップ32で計算したωiを乗算することでスケーリングを
行う。34は非対角行列要素をｍ個の部分行列Ａ1，Ａ2，
・・・，Ａmに分割するステップである。ここでは３次
元７点差分のため６個の部分行列Ａ1，Ａ2，・・・，Ａ
6に分割している。35は部分行列Ａ1，・・・，Ａ6を使
用して前処理行列Ｍを作成するステップである。上記各
ステップは、記憶手段に格納されたデータを参照しなが
ら上記の計算及び処理を実行する手段によって実現され
る。
【００４３】前処理付き共役勾配法の計算手順について
は、図３について上述した通りである。ステップ９は入
力された行列Ａから前処理行列Ｍを作成するステップで
あり、この具体的手順が図７について上述したものであ
る。
は、図３について上述した通りである。ステップ９は入
力された行列Ａから前処理行列Ｍを作成するステップで
あり、この具体的手順が図７について上述したものであ
る。
【００４４】図８は差分法における数値シミュレーショ
ンの手順を示す図であり、単純にするため線形定常解析
の場合の一例を示す。40は解析の前処理を行うステップ
である。41は差分近似するための領域の分割を行うステ
ップである。42は差分法によって連立一次方程式を作成
するステップを示す。43はその連立一次方程式の数値解
を計算するステップである。ステップ43は図３に示す前
処理付き共役勾配法系の解法で計算される。
ンの手順を示す図であり、単純にするため線形定常解析
の場合の一例を示す。40は解析の前処理を行うステップ
である。41は差分近似するための領域の分割を行うステ
ップである。42は差分法によって連立一次方程式を作成
するステップを示す。43はその連立一次方程式の数値解
を計算するステップである。ステップ43は図３に示す前
処理付き共役勾配法系の解法で計算される。
【００４５】図９は３次元７点差分法における非ゼロ帯
行列の構成を示す図である。45は非ゼロ帯行列の構成に
おける位置関係を示す。46は対角行列Ｄの位置を示す。
47は６本の非対角帯行列の位置を示す。48はこの行列が
次元数Ｎの行列であることを示す。
行列の構成を示す図である。45は非ゼロ帯行列の構成に
おける位置関係を示す。46は対角行列Ｄの位置を示す。
47は６本の非対角帯行列の位置を示す。48はこの行列が
次元数Ｎの行列であることを示す。
【００４６】本実施例において、前処理行列を不完全三
角分解でなく（Ｉ−Ａi），ｉ＝１，２，・・・，ｍの
乗算行列にしたのは、並列ベクトル計算機及び超並列計
算機に良く適用できるようにするためである。不完全三
角分解を使用すると、並列度は行列の次元数Ｎより大幅
に小さくなるが、本方式では常に次元数Ｎの並列度を有
す。また前処理行列を（Ｉ−Ａｉ），ｉ＝１，２，・・
・，ｍの乗算で構成したことにより、従来、行われてい
たＥ＝Ａ1＋Ａ2＋・・・＋ＡmとしＩ−Ｅ＋Ｅ2−Ｅ3＋
Ｅ4・・・とするより共役勾配法系の収束性を向上さ
せ、不完全三角分解とほぼ同等の収束速度が得られる。
角分解でなく（Ｉ−Ａi），ｉ＝１，２，・・・，ｍの
乗算行列にしたのは、並列ベクトル計算機及び超並列計
算機に良く適用できるようにするためである。不完全三
角分解を使用すると、並列度は行列の次元数Ｎより大幅
に小さくなるが、本方式では常に次元数Ｎの並列度を有
す。また前処理行列を（Ｉ−Ａｉ），ｉ＝１，２，・・
・，ｍの乗算で構成したことにより、従来、行われてい
たＥ＝Ａ1＋Ａ2＋・・・＋ＡmとしＩ−Ｅ＋Ｅ2−Ｅ3＋
Ｅ4・・・とするより共役勾配法系の収束性を向上さ
せ、不完全三角分解とほぼ同等の収束速度が得られる。
【００４７】また係数行列と右辺ベクトルのスケーリン
グに対角行列だけでなく非対角要素の行ごとの絶対値の
和を合わせて使用することにより、共役勾配法系の収束
を安定化させることができる。特に、対角優位性が大幅
に崩れた行列を係数とする場合に、従来収束しなかった
ものが収束するようになる。
グに対角行列だけでなく非対角要素の行ごとの絶対値の
和を合わせて使用することにより、共役勾配法系の収束
を安定化させることができる。特に、対角優位性が大幅
に崩れた行列を係数とする場合に、従来収束しなかった
ものが収束するようになる。
【００４８】（３）実施例３ 図10は本実施例における共役勾配法系の計算に用いる前
処理用不完全三角分解行列の作成手順を示す図である。
図11は３次元差分法における係数行列と不完全三角分解
行列の構成を示す図である。図12は前処理付き共役勾配
法系の計算手順を示す図である。図13は前処理方式の選
定処理の流れを示すフローチャートである。図14は別の
前処理付き共役勾配法系の計算手順を示す図である。
処理用不完全三角分解行列の作成手順を示す図である。
図11は３次元差分法における係数行列と不完全三角分解
行列の構成を示す図である。図12は前処理付き共役勾配
法系の計算手順を示す図である。図13は前処理方式の選
定処理の流れを示すフローチャートである。図14は別の
前処理付き共役勾配法系の計算手順を示す図である。
【００４９】図10は図12の前処理行列を作成するステッ
プ57の具体的な作成手順の一例を示した図である。図13
は前処理方式の選定を行い結果をキーの値で示したもの
で、図14ではこのキーの値に従った前処理方式を利用し
た共役勾配法系の計算手順の一例を示したものである。
図８の43に示す連立一次方程式の数値解を計算するステ
ップはスカラ計算機のとき図12に示す計算方式で、並列
・ベクトル計算機のとき図14に示す計算方式を適用す
る。
プ57の具体的な作成手順の一例を示した図である。図13
は前処理方式の選定を行い結果をキーの値で示したもの
で、図14ではこのキーの値に従った前処理方式を利用し
た共役勾配法系の計算手順の一例を示したものである。
図８の43に示す連立一次方程式の数値解を計算するステ
ップはスカラ計算機のとき図12に示す計算方式で、並列
・ベクトル計算機のとき図14に示す計算方式を適用す
る。
【００５０】図10の51は係数行列Ａを不完全三角分解す
るための準備処理の一例を示すステップで、Ｕiで行列
各行の非対角要素の絶対値の和を計算し、ωiは不完全
三角分解するときに対角要素と見なす値を計算してい
る。ここでＮは行列の次元数とし、Ｕi，ωiを除く記号
は図11で示した記号に対応する。ｄiは対角要素を、
ａi，ｂi，ｃi，ｅi，ｆi，ｇiは非対角要素を示す。ま
たｄi＞０としておく。52はＧustafssonの提案した補正
方法に従って補正を行った不完全三角分解処理を行うス
テップである。従来は波つきのｄi＝１｛ｄi−・・・｝
と計算していたものをここでは波つきのｄi＝１／｛ωi
−・・・｝と計算する方式を採用する。これにより係数
行列の性質が悪い(対角要素に対して非対角要素の絶対
値の和の値が大きい)場合でも収束が安定し、しかも、
補正係数αの値は行列の性質の悪さにかかわらず、α＝
0.95程度に固定できる。さらに収束性を速めるには、次
元数に比例してαの値を大きくすればよい。但しαは１
より小さい値とする。またａi，ｇiは対角よりｍだけ離
れ、ｂi，ｆiは対角よりｌだけ離れ、ｃi，ｅiは対角よ
り１だけ離れた位置にあるものとする。ステップ52の計
算をベクトル計算機で行うときにはｉ＝１，２，３，・
・・，Ｎの順序で計算せず、ハイパープレーン法と呼ば
れる方向に順序付けて計算すれば良く、その方式は概に
公知技術となっている。
るための準備処理の一例を示すステップで、Ｕiで行列
各行の非対角要素の絶対値の和を計算し、ωiは不完全
三角分解するときに対角要素と見なす値を計算してい
る。ここでＮは行列の次元数とし、Ｕi，ωiを除く記号
は図11で示した記号に対応する。ｄiは対角要素を、
ａi，ｂi，ｃi，ｅi，ｆi，ｇiは非対角要素を示す。ま
たｄi＞０としておく。52はＧustafssonの提案した補正
方法に従って補正を行った不完全三角分解処理を行うス
テップである。従来は波つきのｄi＝１｛ｄi−・・・｝
と計算していたものをここでは波つきのｄi＝１／｛ωi
−・・・｝と計算する方式を採用する。これにより係数
行列の性質が悪い(対角要素に対して非対角要素の絶対
値の和の値が大きい)場合でも収束が安定し、しかも、
補正係数αの値は行列の性質の悪さにかかわらず、α＝
0.95程度に固定できる。さらに収束性を速めるには、次
元数に比例してαの値を大きくすればよい。但しαは１
より小さい値とする。またａi，ｇiは対角よりｍだけ離
れ、ｂi，ｆiは対角よりｌだけ離れ、ｃi，ｅiは対角よ
り１だけ離れた位置にあるものとする。ステップ52の計
算をベクトル計算機で行うときにはｉ＝１，２，３，・
・・，Ｎの順序で計算せず、ハイパープレーン法と呼ば
れる方向に順序付けて計算すれば良く、その方式は概に
公知技術となっている。
【００５１】図11の53は元の係数行列Ａの非ゼロ要素位
置の構成を示したものであり、ｄは対角行列を、ａ，
ｂ，ｃ，ｅ，ｆ，ｇは非対角行列を示す。ａi，ｂi，ｃ
i，ｄi，ｅi，ｆi，ｇiは各行列のｉ行の要素を示すも
のである。54は不完全三角分解された下三角行列Ｌの非
ゼロ要素の位置構成を、55は対角行列Ｄを、56は上三角
行列Ｕの非ゼロ要素の位置構成を示す。ここで行列Ａと
同じ記号ａ，ｂ等は同じ値の行列であることを示し、波
つきのｄだけが新たに計算される。
置の構成を示したものであり、ｄは対角行列を、ａ，
ｂ，ｃ，ｅ，ｆ，ｇは非対角行列を示す。ａi，ｂi，ｃ
i，ｄi，ｅi，ｆi，ｇiは各行列のｉ行の要素を示すも
のである。54は不完全三角分解された下三角行列Ｌの非
ゼロ要素の位置構成を、55は対角行列Ｄを、56は上三角
行列Ｕの非ゼロ要素の位置構成を示す。ここで行列Ａと
同じ記号ａ，ｂ等は同じ値の行列であることを示し、波
つきのｄだけが新たに計算される。
【００５２】図12は前処理付きＣＧＰＭ法（共役勾配最
小多項式法、Ｖon der ＶorstのＢi−ＣＧＳＴＡＢ法
(1990)と同様な計算手法）で連立一次方程式の解ｘを反
復計算して求める計算手順を示す図である。57は入力さ
れた行列Ａから不完全三角分解行列ＬＤＵを作成するス
テップである。ステップ57の具体的手順は図10で示した
通りである。58は反復計算の前準備をするステップであ
る。ここでベクトルｘは解の初期ベクトルで、分からな
い場合はオール・ゼロを入れておけばよい。ベクトルｒ
は残差ベクトルである。また（ベクトルｒ，ベクトル
ｒ）はベクトルの内積計算を示し、結果はスカラ値とな
る。59は反復計算の制御を行うステップである。60は不
完全三角分解付きＣＧＰＭ法の反復計算を実行するステ
ップであり、計算解はベクトルｘに求まる。またｑ＝
〔ＬＤＵ〕~1・波つきのベクトルＰ、（波つきのベクト
ルＰ＝Ａ・波つきのベクトルＰ）は波つきのベクトルＰ
にＬＤＵ行列の逆行列を乗算するのではなく、一般に知
られている前進、後退代入計算を使用して波つきのベク
トルＰからベクトルｑを計算する。ＡとＬ，Ｄ，Ｕは次
元数Ｎの行列、ベクトルｂ，ｘ，ｐ，ｑ，ｒ，ｒ0，ｅ
およびｖは次元数Ｎのベクトルを示し、それ以外の記号
はスカラを示す。
小多項式法、Ｖon der ＶorstのＢi−ＣＧＳＴＡＢ法
(1990)と同様な計算手法）で連立一次方程式の解ｘを反
復計算して求める計算手順を示す図である。57は入力さ
れた行列Ａから不完全三角分解行列ＬＤＵを作成するス
テップである。ステップ57の具体的手順は図10で示した
通りである。58は反復計算の前準備をするステップであ
る。ここでベクトルｘは解の初期ベクトルで、分からな
い場合はオール・ゼロを入れておけばよい。ベクトルｒ
は残差ベクトルである。また（ベクトルｒ，ベクトル
ｒ）はベクトルの内積計算を示し、結果はスカラ値とな
る。59は反復計算の制御を行うステップである。60は不
完全三角分解付きＣＧＰＭ法の反復計算を実行するステ
ップであり、計算解はベクトルｘに求まる。またｑ＝
〔ＬＤＵ〕~1・波つきのベクトルＰ、（波つきのベクト
ルＰ＝Ａ・波つきのベクトルＰ）は波つきのベクトルＰ
にＬＤＵ行列の逆行列を乗算するのではなく、一般に知
られている前進、後退代入計算を使用して波つきのベク
トルＰからベクトルｑを計算する。ＡとＬ，Ｄ，Ｕは次
元数Ｎの行列、ベクトルｂ，ｘ，ｐ，ｑ，ｒ，ｒ0，ｅ
およびｖは次元数Ｎのベクトルを示し、それ以外の記号
はスカラを示す。
【００５３】図13は前処理方式の選定処理の流れを示す
フローチャートで、ＫＥＹ＝０のとき前処理乗算行列Ｍ
の作成を示し、ＫＥＹ＝１のとき不完全三角分解ＬＤＵ
を前処理として使用することを示す。61は対角要素ｄi
と非対角要素の絶対値の和の比の最大値をＳにセットす
るステップを示す。Ｓの値が小さいと行列の性質は良
く、大きいと悪いことを示す。62は非定常計算等で係数
行列Ａが完全な対角優位行列となる場合、すなわちＳ≦
１−ε1のときＫＥＹ＝０をセットするよう判定するス
テップである。一般にε1＝0.1程度の値としておく。63
は係数行列Ａの対角優位性が大きく崩れた場合、すなわ
ちＳ≧ε2のときＫＥＹ＝１をセットするよう判定する
ステップである。一般にε2は２程度の値としておく。6
4は係数行列の性質が特別に良くも悪くもないとき、複
数のベクトル演算器が使用できるかどうかを判定するス
テップで、使用可能のときＫＥＹ＝０を、一台しか使用
できない場合ＫＥＹ＝１をセットする。ＫＥＹ＝０をセ
ットした場合は複数のベクトル演算器を効率良く使用す
る前処理を選定する。
フローチャートで、ＫＥＹ＝０のとき前処理乗算行列Ｍ
の作成を示し、ＫＥＹ＝１のとき不完全三角分解ＬＤＵ
を前処理として使用することを示す。61は対角要素ｄi
と非対角要素の絶対値の和の比の最大値をＳにセットす
るステップを示す。Ｓの値が小さいと行列の性質は良
く、大きいと悪いことを示す。62は非定常計算等で係数
行列Ａが完全な対角優位行列となる場合、すなわちＳ≦
１−ε1のときＫＥＹ＝０をセットするよう判定するス
テップである。一般にε1＝0.1程度の値としておく。63
は係数行列Ａの対角優位性が大きく崩れた場合、すなわ
ちＳ≧ε2のときＫＥＹ＝１をセットするよう判定する
ステップである。一般にε2は２程度の値としておく。6
4は係数行列の性質が特別に良くも悪くもないとき、複
数のベクトル演算器が使用できるかどうかを判定するス
テップで、使用可能のときＫＥＹ＝０を、一台しか使用
できない場合ＫＥＹ＝１をセットする。ＫＥＹ＝０をセ
ットした場合は複数のベクトル演算器を効率良く使用す
る前処理を選定する。
【００５４】図14は二つの前処理手法のどちらかを選定
して、前処理付きＣＧＰＭ法で連立一次方程式の解、ベ
クトルｘ、を反復計算して求める計算手順を示す図であ
る。67は図13に従ってＫＥＹの値を設定するステップで
ある。68はＫＥＹ＝０か１かにより異なった前処理を実
行するステップである。ＫＥＹ＝０のとき図７に示す手
順によって前処理行列Ｍを作成する。69は残差ベクトル
ｒを設定するステップで、ここでもＫＥＹ＝０か１で処
理方式を変える。70は反復計算の前準備をするステップ
である。71は反復計算の制御を行うステップである。7
2，73，74及び75は前処理付きＣＧＰＭ法の反復計算を
実行するステップであり、計算解はベクトルｘに求ま
る。ステップ72と74はＫＥＹ＝０か１により異なった計
算を行う。ステップ23と25は共通な計算を行う。ＡとＭ
及びＬ，Ｄ，Ｕは次元数Ｎの行列、ベクトルｂ，ｘ，
ｐ，ｑ，ｒ，ｒ0，ｅおよびｖは次元数Ｎのベクルを示
し、それ以外の記号はスカラとする。
して、前処理付きＣＧＰＭ法で連立一次方程式の解、ベ
クトルｘ、を反復計算して求める計算手順を示す図であ
る。67は図13に従ってＫＥＹの値を設定するステップで
ある。68はＫＥＹ＝０か１かにより異なった前処理を実
行するステップである。ＫＥＹ＝０のとき図７に示す手
順によって前処理行列Ｍを作成する。69は残差ベクトル
ｒを設定するステップで、ここでもＫＥＹ＝０か１で処
理方式を変える。70は反復計算の前準備をするステップ
である。71は反復計算の制御を行うステップである。7
2，73，74及び75は前処理付きＣＧＰＭ法の反復計算を
実行するステップであり、計算解はベクトルｘに求ま
る。ステップ72と74はＫＥＹ＝０か１により異なった計
算を行う。ステップ23と25は共通な計算を行う。ＡとＭ
及びＬ，Ｄ，Ｕは次元数Ｎの行列、ベクトルｂ，ｘ，
ｐ，ｑ，ｒ，ｒ0，ｅおよびｖは次元数Ｎのベクルを示
し、それ以外の記号はスカラとする。
【００５５】図８は、すでに説明した通り、差分法にお
ける数値シミュレーションの計算手順を示す図である。
ステップ43はその連立一次方程式の数値解を計算するス
テップであり、この部分はスカラ計算機のとき図12に示
す計算手順で、並列ベクトル計算機のとき図14で示す計
算手順で計算する。
ける数値シミュレーションの計算手順を示す図である。
ステップ43はその連立一次方程式の数値解を計算するス
テップであり、この部分はスカラ計算機のとき図12に示
す計算手順で、並列ベクトル計算機のとき図14で示す計
算手順で計算する。
【００５６】本実施例によれば、不完全三角分解付き共
役勾配法系の反復解法で、連立一次方程式の数値解を安
定に計算する効果を得る。特に対角優位性が大きく崩れ
た行列では従来の不完全三角分解では解が収束しないケ
ースが多発していたのをなくすことができる。またＧus
tafsson流の補正における補正係数αを、行列の対角優
位性の良し悪しにかかわらずα＝0.95程度に固定するこ
とで、100×100分割の２次元問題で３〜５倍の収束性の
向上ができる。さらにαを領域の分割数に比例して大き
くする（但しα＜1.0）ことでさらに収束を速くするこ
とができる。
役勾配法系の反復解法で、連立一次方程式の数値解を安
定に計算する効果を得る。特に対角優位性が大きく崩れ
た行列では従来の不完全三角分解では解が収束しないケ
ースが多発していたのをなくすことができる。またＧus
tafsson流の補正における補正係数αを、行列の対角優
位性の良し悪しにかかわらずα＝0.95程度に固定するこ
とで、100×100分割の２次元問題で３〜５倍の収束性の
向上ができる。さらにαを領域の分割数に比例して大き
くする（但しα＜1.0）ことでさらに収束を速くするこ
とができる。
【００５７】また複数台のベクトル演算器を持つスーパ
ーコンピュータにおいて、使用可能なベクトル演算器の
台数と係数行列の性質の良し悪しで二種類の前処理方式
を自動選定させ、連立一次方程式の数値解をより高速に
求めることができる。
ーコンピュータにおいて、使用可能なベクトル演算器の
台数と係数行列の性質の良し悪しで二種類の前処理方式
を自動選定させ、連立一次方程式の数値解をより高速に
求めることができる。
【００５８】
【発明の効果】本発明によれば、複数台のベクトル演算
器を有する計算機や並列処理計算機によく適用できるよ
うな連立一次方程式の前処理行列を作成できるため、こ
れらの計算機に関して連立一次方程式の数値解を高速に
求めることができる。
器を有する計算機や並列処理計算機によく適用できるよ
うな連立一次方程式の前処理行列を作成できるため、こ
れらの計算機に関して連立一次方程式の数値解を高速に
求めることができる。
【００５９】また、連立一次方程式の係数行列の性質が
悪条件であっても、その数値解の収束性を向上させるこ
とができる。
悪条件であっても、その数値解の収束性を向上させるこ
とができる。
【図１】第１の実施例において共役勾配法系の計算手法
に用いる前処理行列Ｍの作成手順を示す図である。
に用いる前処理行列Ｍの作成手順を示す図である。
【図２】前処理行列Ｍを構成する行列の具体的分割方法
を示す図である。
を示す図である。
【図３】前処理付き共役勾配法系の計算手順を示す図で
ある。
ある。
【図４】有限要素法による数値シミュレーションの手順
を示す図である。
を示す図である。
【図５】有限要素分割を示す図である。
【図６】図５に示す分割の結果作成される連立一次方程
式の係数行列の記憶方法を示す図である。
式の係数行列の記憶方法を示す図である。
【図７】第２の実施例において共役勾配法系の計算手法
に用いる前処理行列Ｍの作成手順を示す図である。
に用いる前処理行列Ｍの作成手順を示す図である。
【図８】差分法による数値シミュレーションの手順を示
す図である。
す図である。
【図９】３次元７点差分法における非ゼロ帯行列の構成
を示す図である。
を示す図である。
【図１０】第３の実施例において共役勾配法系の計算手
法に用いる前処理行列Ｍの作成手順を示す図である。
法に用いる前処理行列Ｍの作成手順を示す図である。
【図１１】３次元差分法における係数行列と不完全三角
分解行列の構成を示す図である。
分解行列の構成を示す図である。
【図１２】前処理付き共役勾配法系の計算手順を示す図
である。
である。
【図１３】前処理方式の選定処理の流れを示すフローチ
ャートである。
ャートである。
【図１４】別の前処理付き共役勾配法系の計算手順を示
す図である。
す図である。
フロントページの続き (56)参考文献 特開 平３−28961（ＪＰ，Ａ) 特開 平１−219951（ＪＰ，Ａ) 特開 平１−125667（ＪＰ，Ａ) 特開 昭63−127366（ＪＰ，Ａ) 特開 昭63−95568（ＪＰ，Ａ) 特開 昭62−172461（ＪＰ，Ａ) 情報処理学会研究報告 Ｖｏｌ．93, Ｎｏ．89（1993−10−14）93−ＨＰＣ− 49 ｐｐ．17−24 情報処理学会論文誌 Ｖｏｌ．27，Ｎ ｏ．１（1986−１−10）ｐｐ．11−19 (58)調査した分野(Int.Cl.6，ＤＢ名) G06F 17/12 ＪＩＣＳＴファイル（ＪＯＩＳ)
Claims (4)
- 【請求項１】疎行列を係数行列とする連立１次方程式の
解を計算する装置であり、 前記係数行列のゼロでない要素を格納する記憶手段と、 前記記憶手段に格納された係数行列の各行について、非
対角要素値の絶対値の和または対角要素の値のいづれか
大きいほうの値と対角要素の値との平均で前記係数行列
の各要素をスケーリングするスケーリング手段と、 前記スケーリング手段によりスケーリングされた前記係
数行列の対角要素を除く要素に対して、各行ごとに要素
の値の絶対値が降順に並ぶようにソートするソート手段
と、前記ソート手段によりソートされた前記係数行列を
任意の列で分割し、ｍ個の部分行列Ｅ１、Ｅ２…Ｅｍに
分割する分割手段と、 前記分割手段により分割された前記部分行列の各々を単
位行列Ｉから減算した結果を（Ｉ‐Ｅ１）×（Ｉ‐Ｅ
２）×……×（Ｉ‐Ｅｍ）の形式で乗算する演算手段
と、前記乗算する演算手段によって得られた行列Mを複
数行単位で分割し並列計算機の各プロセッサに割り当て
る割り当て手段とを備えることを特徴とする連立１次方
程式の解を計算する装置。 - 【請求項２】疎行列を係数行列とする連立１次方程式の
解を計算する装置であり、 前記係数行列のゼロでない要素を格納する記憶手段と、 前記記憶手段に格納された係数行列の各行について、非
対角要素値の絶対値の和または対角要素の値のいづれか
大きいほうの値と対角要素の値との平均で前記係数行列
の各要素をスケーリングするスケーリング手段と、 前記スケーリング手段によりスケーリングされた前記係
数行列を1列単位に分割し、ｍ個の部分行列Ａ１、Ａ２
…Ａｍに分割する分割手段と、 前記分割手段により分割された前記部分行列の各々を単
位行列Ｉから減算した結果を（Ｉ‐Ａ１）×（Ｉ‐Ａ
２）×……×（Ｉ‐Ａｍ）の形式で乗算する演算手段
と、前記乗算する演算手段によって得られた行列Mを複
数行単位で分割し並列計算機の各プロセッサに割り当て
る割り当て手段とを備えることを特徴とする連立１次方
程式の解を計算する装置。 - 【請求項３】前記疎行列を係数行列とする連立１次方程
式の解を計算する装置において、 対角要素の値及び非対角要素の絶対値の和から、係数行
列の対角優位性を判定する判定手段と、 前記係数行列が対角優位であった場合に請求項1及び請
求項2に示す連立1次方程式解法装置を選択し、前記係数
行列が対角優位でなかった場合に不完全LU分解を前処理
行列として使用する連立1次方程式解法装置とを選択す
る選択手段を備えることを特徴とする連立1次方程式解
法装置。 - 【請求項４】前記疎行列を係数行列とする連立１次方程
式の解を計算する装置において、 複数個の演算器を利用できるか否かを判定する判定手段
と、 前記複数個のベクトル計算機を使用できる場合に請求項
1及び請求項2に示す連立1次方程式解法装置を使用し、
単一のベクトル演算器しか使用できない場合に不完全LU
分解を前処理行列として使用する連立1次方程式解法装
置を選択する選択手段を備えることを特徴とする連立1
次方程式解法装置。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP3239260A JP2956800B2 (ja) | 1991-09-19 | 1991-09-19 | 連立一次方程式に関する計算装置 |
| US07/947,801 US5604911A (en) | 1991-09-19 | 1992-09-21 | Method of and apparatus for preconditioning of a coefficient matrix of simultaneous linear equations |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP3239260A JP2956800B2 (ja) | 1991-09-19 | 1991-09-19 | 連立一次方程式に関する計算装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPH0581310A JPH0581310A (ja) | 1993-04-02 |
| JP2956800B2 true JP2956800B2 (ja) | 1999-10-04 |
Family
ID=17042120
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP3239260A Expired - Fee Related JP2956800B2 (ja) | 1991-09-19 | 1991-09-19 | 連立一次方程式に関する計算装置 |
Country Status (2)
| Country | Link |
|---|---|
| US (1) | US5604911A (ja) |
| JP (1) | JP2956800B2 (ja) |
Families Citing this family (35)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3639323B2 (ja) * | 1994-03-31 | 2005-04-20 | 富士通株式会社 | メモリ分散型並列計算機による連立１次方程式計算処理方法および計算機 |
| US6078938A (en) * | 1996-05-29 | 2000-06-20 | Motorola, Inc. | Method and system for solving linear systems |
| US6182270B1 (en) * | 1996-12-04 | 2001-01-30 | Lucent Technologies Inc. | Low-displacement rank preconditioners for simplified non-linear analysis of circuits and other devices |
| US5844821A (en) * | 1997-04-29 | 1998-12-01 | Lucent Technologies Inc. | Systems and methods for determining characteristics of a singular circuit |
| TW388921B (en) * | 1997-11-28 | 2000-05-01 | Nippon Electric Co | Semiconductor process device simulation method and storage medium storing simulation program |
| US6154716A (en) * | 1998-07-29 | 2000-11-28 | Lucent Technologies - Inc. | System and method for simulating electronic circuits |
| JP3154406B2 (ja) * | 1998-11-02 | 2001-04-09 | 日本電気株式会社 | シミュレーション方法および記録媒体 |
| US6810370B1 (en) * | 1999-03-31 | 2004-10-26 | Exxonmobil Upstream Research Company | Method for simulation characteristic of a physical system |
| JP3768375B2 (ja) * | 2000-04-04 | 2006-04-19 | Ｎｅｃエレクトロニクス株式会社 | 計算装置および電子回路シミュレーション装置 |
| US7006951B2 (en) * | 2000-06-29 | 2006-02-28 | Object Reservoir, Inc. | Method for solving finite element models using time slabbing |
| US7369973B2 (en) | 2000-06-29 | 2008-05-06 | Object Reservoir, Inc. | Method and system for representing reservoir systems |
| JP3827941B2 (ja) * | 2000-11-16 | 2006-09-27 | 株式会社日立製作所 | 連立一次方程式求解方法及びその実施装置並びにその処理プログラムを記録した記録媒体 |
| JP3809062B2 (ja) * | 2000-11-21 | 2006-08-16 | 富士通株式会社 | マルチレベル不完全ブロック分解による前処理を行う処理装置 |
| JP3639206B2 (ja) * | 2000-11-24 | 2005-04-20 | 富士通株式会社 | 共有メモリ型スカラ並列計算機における並列行列処理方法、及び記録媒体 |
| US20030149962A1 (en) * | 2001-11-21 | 2003-08-07 | Willis John Christopher | Simulation of designs using programmable processors and electronically re-configurable logic arrays |
| US7328195B2 (en) * | 2001-11-21 | 2008-02-05 | Ftl Systems, Inc. | Semi-automatic generation of behavior models continuous value using iterative probing of a device or existing component model |
| US20030182518A1 (en) * | 2002-03-22 | 2003-09-25 | Fujitsu Limited | Parallel processing method for inverse matrix for shared memory type scalar parallel computer |
| GB0208329D0 (en) * | 2002-04-11 | 2002-05-22 | Univ York | Data processing particularly in communication systems |
| US20040236806A1 (en) * | 2002-06-24 | 2004-11-25 | Turner James D. | Method, apparatus and articles of manufacture for computing the sensitivity partial derivatives of linked mechanical systems |
| US7418370B2 (en) * | 2004-03-31 | 2008-08-26 | International Business Machines Corporation | Method, apparatus and computer program providing broadband preconditioning based on reduced coupling for numerical solvers |
| CN100489558C (zh) * | 2004-06-07 | 2009-05-20 | 埃克森美孚上游研究公司 | 用于求解隐式储层仿真矩阵的方法 |
| US20070255778A1 (en) * | 2006-04-27 | 2007-11-01 | Jean-Paul Theis | Software method for solving systems of linear equations having integer variables |
| WO2008026261A1 (fr) | 2006-08-30 | 2008-03-06 | Fujitsu Limited | Procédé de traitement de calcul haute vitesse d'une équation de combinaison à base d'un procédé à éléments finis et d'un procédé à éléments limites |
| US20080120357A1 (en) * | 2006-11-22 | 2008-05-22 | Jean-Paul Theis | Software method for solving systems of linear equations having integer variables |
| WO2009008086A1 (ja) * | 2007-07-12 | 2009-01-15 | Fujitsu Limited | 計算装置、計算方法および計算プログラム |
| BRPI0820870A2 (pt) | 2007-12-13 | 2015-06-16 | Exxonmobil Upstream Res Co | Método para simular um modelo de reservatório. |
| US8175853B2 (en) * | 2008-03-28 | 2012-05-08 | International Business Machines Corporation | Systems and methods for a combined matrix-vector and matrix transpose vector multiply for a block-sparse matrix |
| US8306798B2 (en) * | 2008-12-24 | 2012-11-06 | Intel Corporation | Fluid dynamics simulator |
| US9176928B2 (en) * | 2009-07-07 | 2015-11-03 | L3 Communication Integrated Systems, L.P. | System for convergence evaluation for stationary method iterative linear solvers |
| US8577949B2 (en) * | 2009-07-07 | 2013-11-05 | L-3 Communications Integrated Systems, L.P. | System for conjugate gradient linear iterative solvers |
| US8712738B2 (en) * | 2010-04-30 | 2014-04-29 | International Business Machines Corporation | Determining ill conditioning in square linear system of equations |
| WO2014104093A1 (ja) * | 2012-12-28 | 2014-07-03 | 株式会社村田製作所 | 連立１次方程式の解の算出装置及び算出方法、並びにそれに適用されるプログラム |
| US10242136B2 (en) * | 2015-05-20 | 2019-03-26 | Saudi Arabian Oil Company | Parallel solution for fully-coupled fully-implicit wellbore modeling in reservoir simulation |
| CN106997333B (zh) * | 2016-01-22 | 2020-07-28 | 阿里巴巴集团控股有限公司 | 逻辑回归梯度的计算方法和装置 |
| CN113609720B (zh) * | 2021-07-07 | 2022-10-25 | 广州中望龙腾软件股份有限公司 | 一种有限元分析的主从自由度处理方法、设备及存储介质 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4697247A (en) * | 1983-06-10 | 1987-09-29 | Hughes Aircraft Company | Method of performing matrix by matrix multiplication |
| US4787057A (en) * | 1986-06-04 | 1988-11-22 | General Electric Company | Finite element analysis method using multiprocessor for matrix manipulations with special handling of diagonal elements |
| US5301342A (en) * | 1990-12-20 | 1994-04-05 | Intel Corporation | Parallel processing computer for solving dense systems of linear equations by factoring rows, columns, and diagonal, inverting the diagonal, forward eliminating, and back substituting |
-
1991
- 1991-09-19 JP JP3239260A patent/JP2956800B2/ja not_active Expired - Fee Related
-
1992
- 1992-09-21 US US07/947,801 patent/US5604911A/en not_active Expired - Lifetime
Non-Patent Citations (2)
| Title |
|---|
| 情報処理学会研究報告 Ｖｏｌ．９３，Ｎｏ．８９（１９９３−１０−１４）９３−ＨＰＣ−４９ ｐｐ．１７−２４ |
| 情報処理学会論文誌 Ｖｏｌ．２７，Ｎｏ．１（１９８６−１−１０）ｐｐ．１１−１９ |
Also Published As
| Publication number | Publication date |
|---|---|
| JPH0581310A (ja) | 1993-04-02 |
| US5604911A (en) | 1997-02-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP2956800B2 (ja) | 連立一次方程式に関する計算装置 | |
| Krishnan et al. | Efficient preconditioning of laplacian matrices for computer graphics | |
| US20170061279A1 (en) | Updating an artificial neural network using flexible fixed point representation | |
| US7933752B2 (en) | Method, apparatus and computer program providing broadband preconditioning based on a reduced coupling for numerical solvers | |
| KR20190062303A (ko) | 폴드된 특징 데이터에 대한 컨볼루션 연산을 수행하기 위한 방법 및 장치 | |
| JP2017130036A (ja) | 情報処理装置、演算方法、および演算プログラム | |
| WO2010109359A2 (en) | Processing of linear systems of equations | |
| KR101990735B1 (ko) | 사전 그래프 분할 기반 행렬 벡터 곱을 이용한 대규모 그래프 마이닝 방법 및 장치 | |
| Xu et al. | Towards a scalable hierarchical high-order CFD solver | |
| JP2024028901A (ja) | ハードウェアにおけるスパース行列乗算 | |
| Oliker et al. | Parallel implementation of an adaptive scheme for 3D unstructured grids on the SP2 | |
| Merkel et al. | Iterative solution of large-scale 3D-BEM industrial problems | |
| JP3809062B2 (ja) | マルチレベル不完全ブロック分解による前処理を行う処理装置 | |
| Rauber et al. | Algorithms for Systems of Linear Equations | |
| JP2007133710A (ja) | 連立一次方程式反復解法における前処理方法および行列リオーダリング方法 | |
| Penzl | A cyclic low rank Smith method for large, sparse Lyapunov equations with applications in model reduction and optimal control | |
| Bar-On et al. | A fast parallel Cholesky decomposition algorithm for tridiagonal symmetric matrices | |
| JPH06119368A (ja) | 緩和法による連立方程式解析計算機 | |
| JP3300388B2 (ja) | 並列プロセス生成方法 | |
| Bevilacqua et al. | Orthogonal iterations on Structured Pencils | |
| Nepomniaschaya et al. | Associative parallel algorithm of checking spanning trees for optimality | |
| Torun et al. | Enhancing Block Cimmino for Sparse Linear Systems with Dense Columns via Schur Complement | |
| Lau | Numerical solution of skew-symmetric linear systems | |
| CN116702857A (zh) | 计算机可读记录介质、机器学习方法和信息处理设备 | |
| JPH07325806A (ja) | マルチプロセッサシステムのプロセス割当方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20070723 Year of fee payment: 8 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20080723 Year of fee payment: 9 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20080723 Year of fee payment: 9 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090723 Year of fee payment: 10 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090723 Year of fee payment: 10 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100723 Year of fee payment: 11 |
|
| LAPS | Cancellation because of no payment of annual fees |