JP2005115497A

JP2005115497A - 連立一次方程式反復求解計算機、及び、連立一次方程式反復求解計算方法

Info

Publication number: JP2005115497A
Application number: JP2003346402A
Authority: JP
Inventors: Takumi Washio; 巧鷲尾
Original assignee: NEC Corp
Current assignee: NEC Corp
Priority date: 2003-10-03
Filing date: 2003-10-03
Publication date: 2005-04-28
Anticipated expiration: 2023-10-03
Also published as: JP3865247B2

Abstract

【課題】全体系方程式の求解の高速化と記憶容量の低減。
【解決手段】連立一次方程式を作成する作成計算器５と、連立一次方程式を前処理つき反復解法で解く反復求解器５と、反復過程の中の前処理演算において、与えられた右辺ベクトルＦとＧに対して、主未知数ベクトルＸｖを近似的に消去する消去計算器５と、主未知数ベクトルＸｖが近似的に消去され拘束的未知数ベクトルＸｐが未知数として残る主未知数消去後拘束方程式を内部の反復求解により解いて拘束的未知数ベクトルの近似解Ｖを求める第１近似計算器５と、近似解Ｖを主方程式に代入して主未知数ベクトルの近似解Ｖを求める第２近似計算器５とから構成されている。全体系の行列に対して近似精度の良い前処理行列を構成する際に現れる密な係数行列をもつ主未知数消去後拘束方程式を疎行列演算のみで効率的に解くことにより、収束解を得るまでの演算量を削減して高速化を実現する。
【選択図】図１

Description

本発明は、連立一次方程式反復求解計算機に関し、特に、流体、機械構造、電磁場のような物理状態を数学的に解析する混合型有限要素離散化法を適用する際に現れる大規模連立一次方程式を安定的に、且つ、高速に解く連立一次方程式反復求解計算機、及び、連立一次方程式反復求解計算方法に関する。

拘束条件付き線形・非線形方程式の離散化により現れる大規模連立一次方程式、又は、流体、構造、電磁場、それらの相互作用を扱う連成解析で適用される混合型有限要素離散化法で現れる大規模連立一次方程式を安定的に、且つ、高速に解くことは、高度な科学技術の発展のために求められている。

拘束条件付き方程式が離散化される場合に生じる係数行列の拘束条件式に対応する行の優対角性が保証されないので、離散化による連立方程式は標準的な反復解法では解くことができないことが知られている。質量保存、電荷量保存のような保存則に基づく連続の式が拘束条件である流体問題、構造問題、又は、電磁場問題を解く際に、主方程式と拘束方程式が合わさる混合型有限要素法により離散化が行われる場合に、その連立方程式に現れる係数行列の拘束方程式に対応する行の対角成分が零又は零に非常に近い値になるので、その離散化により生ずる膨大な数の連立方程式は標準的な反復解法では解くことができないことが知られている。

このような困難性を排除するために、離散化の際に、inf-sup条件と呼ばれる安定条件を満たす離散化手法が適用され、Uzawaアルゴリズム、又は、これを修飾するInexact Uzawaアルゴリズムが用いられる手法（特許文献１参照）、又は、拘束方程式が最後尾に並べられて係数行列が作成され不完全ＬＵ分解の際に生じる非零の対角成分フィルインを取り込み、その不完全ＬＵ分解行列を前処理付き反復法の前処理演算における前処理行列として用いる手法（非特許文献２）が採り入れられている。

非圧縮性流体解析では、ＰＳＰＧ（Pressure Stabilization Petrov Galcrkin）法のようなinf-sup条件が満たされない圧力変数に安定性をもたらしつつ非零成分の密度を適度に抑えることができる離散化法、又は、ＭＡＣ（Marker and Cell）法のように流体の運動方程式に対応する主方程式と連続の式に対応する拘束方程式の組合せを主方程式と圧力に関するポアッソン方程式の組合せに変換した後に離散化する手法が知られている。これらの２手法は、通常、２組の方程式を交互に解く準陰解法が適用されている。

inf-sup条件を満たす有限要素離散化法が適用される場合に、単純な有限要素離散化法との比較で、要素内の節点の数が増加するため離散化時の演算量が大幅に増加し、且つ、係数行列の中の非零成分の密度が増加する。このような増加は、反復法の演算量の増加を招いて高速性の点で好ましくない。ＰＳＰＧ法、又は、ＭＡＣ法により問題を離散化する離散化手法では、２組の方程式を同時に高速に解くことができる反復解法は提案されていない。２組の方程式を交互に解く準陰解法が適用される場合には、流体と構造が強く相互作用する連成問題では収束解を得ることが困難である。主方程式の中の主部未知数を近似的に消去し、拘束的未知数のみからなる変換された拘束方程式を近似的に解くことにより、前処理演算を実装する手法が提案的に考えられているが、そのような手法では、変換後の拘束方程式の係数行列が密行列になって、その方程式を近似的に解くことは困難であることが知られている。

主方程式と拘束方程式から形成される全体系方程式の収束解を適度な演算量で得ることができることが求められ、特に、inf-sup条件のような安定化条件を満たさない簡易な混合型有限要素離散化を拘束条件付き物理問題に対して実行する場合、又は、ＰＳＰＧ法のような簡易な安定化手法を拘束条件付き物理問題に対して適用する場合に、主方程式と拘束方程式から形成される全体系方程式の収束解を適度な演算量で得ることができることが求められる。

James H. Bramble, Joseph E. Pasciak and Apostol T. Vassilev, Analysis of inexact UZAWA algorithm for saddle point problems, SLAM J. Numer, Anal. 34, (1997), pp. 1072-1092 G. Segal and K. Vuik. A Simple iterative linear solver for the 3D incompressible Navier-Stokes equations discretized by the finite element method, Faculty of Technical Mathematics and Informatics, Delft University of Technology, 1995, TUD Report 95-37

本発明の課題は、主方程式と拘束方程式から形成される全体系方程式の収束解を適度な演算量で得る連立一次方程式反復求解計算機を提供することにある。
本発明の他の課題は、特に、inf-sup条件のような安定化条件を満たさない簡易な混合型有限要素離散化を拘束条件付き物理問題に対して実行する場合、又は、ＰＳＰＧ法のような簡易な安定化手法を拘束条件付き物理問題に対して適用する場合に、主方程式と拘束方程式から形成される全体系方程式の近似解を適度な演算量で計算し、それを前処理付き反復法の前処理演算として適用することにより収束解を適度な演算量で得る連立一次方程式反復求解計算機を提供することにある。

本発明による連立一次方程式反復求解計算機は、現象を記述する主方程式と前記現象を拘束する拘束方程式とを分離し、且つ、前記主方程式と前記拘束方程式をそれぞれに多変数次元メッシュ上で離散化して下記の２式：
物理方程式：ＡｖｖＸｖ＋ＡｖｐＸｐ＝Ｂｖ
拘束方程式：ＡｐｖＸｖ＋ＡｐｐＸｐ＝Ｂｐ
Ａｖｖ，Ａｖｐ，Ａｐｖ，Ａｐｐ：係数行列
Ｂｖ，Ｂｐ：既知ベクトル
Ｘｖ，Ｘｐ：求解対象未知ベクトル
として連立一次方程式を作成する作成計算器（５）と、前記方程式を前処理つき反復解法で解く反復求解計算器で構成され、反復過程の中の前処理演算において前記１式および前記２式の右辺ベクトルＢｖおよびＢｐの代わりに与えられる右辺ベクトルＦおよびにＧに対して、主未知数ベクトルＸｖを近似的に消去し、主未知数消去後拘束方程式の右辺ベクトルＱを作成する消去計算器（５）と、２式から主未知数ベクトルＸｖが近似的に消去され拘束的未知数ベクトルＸｐが未知数として残る主未知数消去後拘束方程式を内側において前処理つき反復解法を実行することにより拘束的未知数ベクトルの近似解Ｐを求める第１近似計算器（５）と、近似解Ｐを主方程式に代入して主未知数ベクトルＸｖの近似解Ｖを求める第２近似計算器（５）とから構成されている。

全体系の方程式が主方程式と拘束方程式とに分離され、主方程式に対応する主未知数ベクトルが近似的に消去され、拘束方程式が拘束的未知数ベクトルのみが未知数として残存する主未知数消去後拘束方程式に変換される。このような変換後の主未知数消去後拘束方程式の係数行列は一般に密行列となるが、適当な疎行列をもとに作成したその近似行列を前処理行列とする前処理つき反復解法を適用することにより、この密行列の成分を陽に保持することなく、主未知数消去後拘束方程式を近似的に解くことができる。このことが、精度の高い全体系方程式に対する近似解法を適切な演算量で実現できる理由である。混合型有限要素離散化法が適用される際に現れる大規模疎行列連立一次方程式を解く際の前処理つき反復法の前処理演算として、前記近似解法を適用することにより、全体系の方程式は高速に、且つ、安定に解かれ得る。主方程式としては、物理現象を記述する物理方程式に限られず、経済現象を記述する経済学方程式が好適に例示される。

消去計算器（５）は、主方程式を下記の２式：
ＭｖｖＶ＋ＡｖｐＰ＝Ｆ
ＡｐｖＶ＋ＡｐｐＰ＝Ｇ
Ｆ，Ｇ：全体系に対する前処理つき反復解法の反復過程の中で前処理演算の右辺ベクトルとして作成される既知ベクトル
のうちの第１式の初期近似解を下記式：
ＭｖｖＶ’＝Ｆ
の解としてＶ’を求め、主未知数消去後拘束方程式の右辺ベクトルＱを下記式：
Ｑ＝Ｇ−ＡｐｖＶ’
より作成する。

第１近似計算器は、下記式：
Ｃｐｐ＝Ａｐｐ−ＡｐｖＭｖｖ∧（−１）Ａｖｐ
で定義される係数行列により定義される主未知数消去後拘束方程式：
ＣｐｐＰ＝Ｑ
を近似的に解くことにより前記拘束的未知数ベクトルの近似解Ｐを求める。

主未知数消去後拘束方程式の係数行列Ｃｐｐに対して、Ｍｖｖの代わりに対角行列Ｄを代入し、Ｃｐｐを近似する疎行列Ｑｐｐ：
Ｑｐｐ＝Ａｐｐ−ＡｐｖＤ∧（−１）Ａｖｐ
の不完全ＬＵ分解行列を前処理行列として前処理つき反復解法を適用することにより、主未知数消去後拘束方程式の近似解を計算する。これにより、一般に密行列となってしまう係数行列Ｃｐｐの成分を陽に保持することなく、適度な演算量でその近似解を得ることができる。

主未知数消去後拘束方程式の係数行列Ｃｐｐに対して、Ｃｐｐを近似する疎行列Ｑｐｐ：
Ｑｐｐ＝Ａｐｐ−ＴｐｖＤ∧（−１）Ｔｖｐ
をＡｖｖの不完全ＬＵ分解行列：
Ｍｖｖ＝（Ｌｖ＋Ｄｖ）Ｄｖ∧（−１）（Ｄｖ＋Ｕｖ）
Ｕｖ：狭義上三角行列
Ｌｖ：狭義下三角行列
Ｄｖ：対角行列
に対して次式：
Ｔｐｖ（Ｄｖ＋Ｕｖ）＝Ａｐｖ
（Ｌｖ＋Ｄｖ）Ｔｖｐ＝Ａｖｐ
を近似的に満たす疎行列ＴｐｖおよびＴｖｐにより作成し、Ｑｐｐの不完全ＬＵ分解行列を前処理行列として前処理つき反復解法することにより、主未知数消去後拘束方程式の近似解を計算する。これにより、一般に密行列となってしまう係数行列Ｃｐｐの成分を陽に保持することなく、適度な演算量でその解を得ることができる。
このような作成は、異方性がある現象に対して有効である。

本発明による連立一次方程式反復求解計算方法は、既述の計算器を用いることにより下記のステップ：
主方程式を離散化する主離散化式と拘束方程式を離散化する拘束離散化式を計算器の中に記述して全体系を記述するステップ、主離散化式に対応する主未知数のうちの主未知数を近似的に計算器の中で消去することにより拘束離散化式を拘束離散化式に対応する拘束的未知数で表される主未知数消去後拘束離散化式で前記計算器の中で記述するステップ、主未知数消去後拘束離散化式の係数行列を疎行列化することにより疎行列化主未知数消去後拘束離散化式を計算器の中で記述するステップ、主未知数消去後拘束離散化式を前処理つき反復解法により解くときの前処理演算として疎行列化主未知数消去後拘束離散化式を計算器の中で記述するステップにより、流体、超弾性体に現れる運動方程式を表現する主方程式と連続の式を表現する拘束方程式を同時に満たす近似解を疎行列演算のみ高速な演算で計算でき、これを前処理つき反復解法の前処理演算として適用することにより、適度な演算量および少ない反復回数で収束が達成され、高度科学技術計算に寄与することができる。詳細に具体的に後述されるように、前処理つき反復解法の反復過程の中の前処理演算において主未知数を近似的に消去した後の主未知数消去後拘束離散化式を解くにあたって、内側においても前処理つき反復解法を用いる。これにより、主未知数消去後拘束離散化式の密な係数行列を陽に保持する必要がなくなり適度な演算量で近似精度の高い全体系の前処理演算が実現でき、外側の前処理つき反復解法に対して高速な収束性が実現される。

本発明による連立一次方程式反復求解計算機は、前処理つき反復解法の反復過程の中で用いる近似精度の良い前処理演算を実現するにあたって現れる密な係数行列を持つ方程式を再び内側の疎行列演算のみの前処理つき反復解法で解くものである。これにより全体系の係数行列に対して近似精度の良い前処理演算が実現でき、外側反復が高速に収束する。特に、inf-sup条件を満たさない混合型有限要素離散化において主方程式と拘束方程式が一体化した全体系方程式に対して有効な前処理演算を実現することができる。

本発明による連立一次方程式反復求解計算機の実現態は、図に対応して、更に詳細に記述される。その実現態は、図１に示されるように、ベクトルデータ保存器１が反復求解計算器２とともに設けられている。ベクトルデータ保存器１は、反復演算に必要であるベクトルデータ３を保存する。反復求解計算器２は、行列ベクトル積演算器４と、近似行列求解器５と、ベクトル積和演算器６と、内積演算器７とから構成されている。行列ベクトル積演算器４には、係数行列データ８を記憶する係数行列データ保持器９を備えている。近似行列求解器５は、近似行列に対して逆行列演算を実行するためのデータ１１を保存する近似行列保存器１２を備えている。

行列ベクトル積演算器４は、ベクトルデータ保存器１から与えられるベクトルデータ３と係数行列データ８との積として行列ベクトル積１３を計算により求める。係数行列データ８は、係数行列データ保持器９に蓄積されている複数の行列データのうちから指定されている行列データである。行列ベクトル積演算器４により計算された行列ベクトル積１３は、ベクトルデータ保存器１に戻されてベクトルデータ３で記憶される。近似行列求解器５は、ベクトルデータ保存器１から与えられるベクトルデータを右辺ベクトルとし、近似行列保存器１２に蓄積されている近似行列に関わるデータにより指定される近似行列を係数とする連立一次方程式を解いてその連立一次方程式解１４を求める。連立一次方程式解１４は、近似行列求解器５からベクトルデータ保存器１に戻される。ベクトル積和演算器６は、ベクトルデータ保存器１から与えられる単数又は複数のベクトルデータと単数又は複数のスカラデータに対して、指定されるベクトル積和演算（ベクトルのスカラ倍とベクトルの足し算との組合せ）を実行して、ベクトル積和１５を求める。ベクトル積和１５は、ベクトル積和演算器６からベクトルデータ保存器１に戻される。内積演算器７は、ベクトルデータ保存器１から与えられる２つのベクトルデータの内積１６を計算する。内積１６は、内積演算器７からベクトルデータ保存器１に戻される。

連立一次方程式を作成する作成計算器５で記述されている連立一次方程式の主方程式は、下記のように表される。
ＡｖｖＸｖ＋ＡｖｐＸｐ＝Ｂｖ・・・（１）
作成計算器５で記述されている連立一次方程式の拘束方程式は、下記のように表される。
ＡｐｖＸｖ＋ＡｐｐＸｐ＝Ｂｐ・・・（２）
式（１）と式（２）とから、主未知数ベクトルＸｖと拘束的未知数ベクトルＸｐの組み合わせが前処理つき反復解法により求められ得る。係数Ａｖｖと係数Ａｖｐは、主方程式の係数行列を表し、係数Ａｐｖと係数Ａｐｐは、拘束方程式の係数行列を表わしている。Ｂｖは主方程式の右辺ベクトルを表し、Ｂｐは拘束方程式の既述の右辺ベクトルを表わしている。

図２は、式（１）と式（２）を同時に満たす主未知数ベクトルＸｖと拘束的未知数ベクトルＸｐの組み合わせを求める前処理つき反復解法の計算フローを示している。初期ステップＳ０では、式（１）と式（２）の右辺ベクトルＢｖ，Ｂｐと未知数ベクトルの適当な初期解Ｘｖ，Ｘｐとが反復求解計算器２に入力される。外部から適当な初期解が与えられない場合には、未知数ベクトルＸｖ，Ｘｐの適正な初期解としてゼロベクトルが代入される。

前処理つき反復解法の一反復は、ベクトルの積和と、内積演算と、係数行列Ａｖｖ，Ａｖｐ，Ａｐｖ，Ａｐｐに対する行列ベクトル積演算と、係数行列の近似行列に対して連立一次方程式を解く前処理演算とで構成される。その前処理演算は、それ以前の積和と内積演算より算出されるＦ，Ｇを右辺とし、式（１）の係数行列Ａｖｖをその逆行列演算を少ない演算量で実行することができる近似行列Ｍｖｖに置き換えた下記の連立方程式を近似的に解く演算処理である。
ＭｖｖＶ＋ＡｖｐＰ＝Ｆ・・・（１’）
ＡｐｖＶ＋ＡｐｐＰ＝Ｇ・・・（２’）

式（１’）にＭｖｖの逆行列Ｍｖｖ∧（−１）が掛けられて未知数ベクトルＶが未知数Ｐで表される。以下、”∧（−１）”は、Ａ∧（−１）がＡの逆行列を表すことを示す記号として用いられる。
Ｖ＝Ｍｖｖ∧（−１）Ｆ−Ｍｖｖ∧（−１）ＡｖｐＰ
これが式（２’）に代入されて、主未知数ベクトルＶが消去され、拘束的未知数ベクトルＰのみからなる主未知数消去後の拘束方程式が得られる。
（Ａｐｐ−ＡｐｖＭｖｖ∧（−１）Ａｖｐ）Ｐ
＝Ｇ−ＡｐｖＭｖｖ∧（−１）Ｆ・・・（３）
この式の左辺の第１因数がＣｐｐで表される。
Ｃｐｐ＝Ａｐｐ−ＡｐｖＭｖｖ∧（−１）Ａｖｐ・・・（４）

式（３）を解くにあたって、その右辺ベクトルを計算するため、まず式（１’）の近似初期解Ｖ’を
Ｖ’＝Ｍｖｖ∧（−１）Ｆ・・・（５）
により計算する。
式（３）と式（４）と式（５）とから、
ＣｐｐＰ＝Ｇ−ＡｐｖＭｖｖ∧（−１）Ｆ
＝Ｇ−ＡｐｖＶ’・・・（６）
式（６）の解Ｐを式（１’）に代入すると、式（５）から、
Ｖ＝Ｍｖｖ∧（−１）Ｆ−Ｍｖｖ∧（−１）ＡｖｐＰ
＝Ｖ’−Ｍｖｖ∧（−１）ＡｖｐＰ・・・（７）

前処理つき反復解法においては、前処理演算の後、その前処理演算の結果得られるベクトルＶ、Ｐにもとの全体系に対する方程式（１）（２）の係数行列を行列ベクトル積演算器４により掛ける。
Ｚｖ＝ＡｖｖＶ＋ＡｖｐＰ・・・（１”）
Ｚｐ＝ＡｐｖＶ＋ＡｐｐＰ・・・（２”）
このように求められるＺｖ，Ｚｐと反復過程の中で生成されるその他のベクトルに対して内積と積和演算がベクトル積和演算器６と内積演算器７で実行され、ステップＳ４で、求解対象未知ベクトルＸｖ，Ｘｐが更新される。更新された求解対象未知ベクトルＸｖ，Ｘｐに対してステップＳ５で収束判定が行われ、収束していなければ再びステップＳ１にもどる。このような反復計算ステップＳ１〜Ｓ４が繰り返され、解Ｘｖ，Ｘｐが収束すれば（ステップＳ５）、ステップＳ６で、その収束解Ｘｖ，Ｘｐが式（１）と式（２）の未知数解Ｘｖ，Ｘｐとして収束的に定められる。

図３は、ステップＳ２において、式（１’）および式（２’）の近似解Ｖ，Ｐを求める計算過程を示している。図３のステップＳ７では、式（５）が計算され、ステップＳ８では、式（６）の右辺が計算され、ステップＳ９では、式（６）の近似解Ｐが計算され、ステップＳ１０では、式（７）の右辺第２項の因数（Ｗ＝−ＡｖｐＰ）が計算され、ステップＳ１１では、式（７）の右辺第２項のＵ＝−Ｍｖｖ∧（−１）ＡｖｐＰが計算される。このＵが式（７）に代入されて、ステップＳ１２でＶが更新される。ここで計算されたＶとＰを図２のステップＳ２での近似解Ｖ，Ｐとする。ステップＳ７、Ｓ１１でのＭｖｖに対する逆行列演算は図１の近似行列求解器５を用いて計算する。ステップＳ８、Ｓ１０での係数行列Ａｐｖ，Ａｖｐに対する行列ベクトル積演算は図１の行列ベクトル積演算器を用いて計算する。その他のベクトルの和は、図１のベクトル積和演算器６を用いて計算する。

式（４）で定義される行列Ｃｐｐの右辺第２項の因数Ｍｖｖ∧（−１）は一般に密行列であるので、行列Ｃｐｐは大概の場合で密行列である。このため、Ｃｐｐの成分を陽に求めることは記憶容量と演算量の両面で得策ではない。記憶容量と演算量の増大を抑制するために、ステップＳ９でＣｐｐに関する連立一次方程式を解く際に、Ｃｐｐの成分が陽に現れない前処理付き反復解法が内部で用いられる。

図４は、そのような内部前処理付き反復解法の計算フローを示している。ここでは、内部反復で用いられる内部反復用前処理行列Ｍｐｐとして、式（４）で定義される行列Ｃｐｐに代えられて、下記式：
Ｑｐｐ＝Ａｐｐ−ＡｐｖＤ∧（−１）Ａｖｐ・・・（８）
で定義される疎行列Ｑｐｐを計算し、それの不完全ＬＵ分解行列が用いられる。式（８）では、式（４）のＭｖｖ∧（−１）がＤ∧（−１）に変更されている。行列Ｄ∧（−１）は、行列Ｍｖｖ∧（−１）が適当に近似されていて、特に、対角行列である。前処理演算のためにＣｐｐを近似する行列Ｍｐｐは、このようなＱｐｐの不完全ＬＵ分解行列として与えられる。

式（６）の右辺ベクトルとして与えられたＱを入力したのち、ステップＳ１３においてＰの初期解としてゼロベクトルを設定し、内部反復回数ＩＴをゼロに初期化する。ステップ１４において、内部反復回数ＩＴをインクリメントした後、ステップＳ１５において内部反復で生成されるベクトルに対して内積および積和演算を行いステップＳ１６の前処理演算に用いられるベクトルＹを計算する。

ステップＳ１６では、Ｘ＝Ｍｐｐ∧（−１）Ｙを計算し、得られたベクトルＸに対してステップＳ１７では、式（６）の係数行列Ｃｐｐを掛ける。ここで、Ｃｐｐを掛けるにあたってその行列成分を用いるのではなく、Ｃｐｐの定義式（４）より得られる次式：
Ｗ＝ＣｐｐＸ＝（Ａｐｐ−ＡｐｖＭｖｖ∧（−１）Ａｖｐ）Ｘ・・・（９）
にしたがって、Ａｐｐ、Ａｐｖ、Ａｖｐに対する行列ベクトル積演算とＭｖｖに対する逆行列演算により、Ｗを計算する。より具体的には、ステップＳ１７に示すように、
Ｚ＝ＡｖｐＸ
を計算し、次に
Ｔ＝Ｍｖｖ∧（−１）Ｚ
を求め、最後に
Ｗ＝ＡｖｐＸ−ＡｐｖＴ
を計算することにより、式（９）のＷを求める。このようにして求めたＷと内部反復過程で生成されるその他のベクトルに対して内積と積和演算がベクトル積和演算器６と内積演算器７で実行され、ステップＳ１８で、内部反復における求解対象ベクトルＰが更新される。ステップＳ１６でのＭｐｐに対する逆行列演算、ステップＳ１７でのＭｖｖに対する逆行列演算は、図１の近似行列求解器５により実行する。また、ステップＳ１７でのＡｖｐ、Ａｐｐ、Ａｐｖに対する行列ベクトル積演算は、行列ベクトル積演算器４にて実行する。

一般にＭｖｖ∧（−１）が密行列となるために、Ｃｐｐも密行列となり、Ｃｐｐの行列成分を陽に用いてそれを係数行列とする連立方程式の求解を行うことは、大規模問題に対して記憶容量および計算量の点から現実的ではない。そこでステップＳ１７では、Ｃｐｐに対する行列ベクトル積がその行列成分を直接用いることなしに実行される。

ここでの内側反復は、図２の外側反復においてステップＳ２が実行されるたびに実行する必要があるので、演算量を適度な量に抑える必要がる。ステップＳ１４で設定される反復回数ＩＴが設定数ＩＴｍａｘに達すれば、内部反復計算のステップＳ１５〜Ｓ１８は強制終了され、その演算量の過度の増大が抑制される。

図５は、本発明による連立一次方程式反復求解計算機の具体的な実施例を示している。本実施例では、非圧縮性流体方程式が混合型有限要素法で離散化されて線形化される。そのように線形化される非圧縮性流体方程式は、既述の式（１）と式（２）で表される。本実施例では、式（１）で表される主方程式は流体の運動方程式の離散化式に一致し、その式の主未知数ベクトルＸｖは流速ベクトルＶ（＝（Ｖｘ，Ｖｙ））に一致し（対応し）ている。式（２）で表される拘束方程式は、流体の非圧縮性に起因する体積保存を表す連続の式の離散化式に一致していて、完全非圧縮の場合には、Ａｐｐ＝０になる。拘束的未知数ベクトルＶｐは圧力ベクトルＰ（＝（Ｐｘ，Ｐｙ））に一致する（対応する）。一般的には、流速と圧力は３次元化される。前処理つき反復求解計算器２により実行されるステップＳ２の前処理演算では、そのＭｖｖとして、Ａｖｖの不完全ＬＵ分解行列が用いられる。ステップＳ２の内部で実行するステップＳ９において式（６）をステップＳ１３〜Ｓ１９の内部の前処理つき反復解法で解く際に、式（８）の対角行列Ｄは、行列Ａｖｖの対角成分が取り出され、又は、行列Ａｖｖの各行の絶対値の和を対応する行の対角成分に代置することにより定められ、その不完全ＬＵ分解行列により、ステップＳ１６の計算で用いられるＣｐｐの近似行列Ｍｐｐが定められる。

図６は、本発明による連立一次方程式反復求解計算機の具体的な他の実施例を示している。本実施例では、メッシュ上で、ゴム、生体の筋肉のような柔らかい材質を有する超弾性体に関する超弾性体方程式が離散化されて線形化される。そのように線形化される超弾性体方程式は、既述の式（１）と式（２）で表される。本実施例では、式（１）で表される主方程式は超弾性体の運動方程式の離散化式に対応し、式（２）の拘束方程式は非圧縮性条件に対応し、その主方程式の主未知数ベクトルＸｖは変位ベクトルＶ（＝（Ｕｘ，Ｕｙ））に対応し、拘束的未知数ベクトルＸｐは圧力ベクトルＰ（＝（Ｐｘ，Ｐｙ））に対応する。ステップＳ２の前処理では、そのＭｖｖとして、Ａｖｖの不完全ＬＵ分解行列が用いられる。ステップＳ１６で用いられるＣｐｐに対する前処理行列Ｍｐｐの作成では、筋肉のように繊維に沿う異方性がある場合には、式（８）に示されるように行列Ａｖｖから定められる対角行列を用いるだけではＣｐｐにも引き継がれる異方性を近似する前処理行列Ｍｐｐを作成することができない。

このような場合には、式（８）のＱｐｐの不完全ＬＵ分解によりＭｐｐを作成するＭｐｐ作成方法は、下記のように改善される。Ａｖｖの不完全ＬＵ分解Ｍｖｖは、
Ｍｖｖ＝（Ｌｖ＋Ｄｖ）Ｄｖ∧（−１）（Ｄｖ＋Ｕｖ）・・・（１０）
により設定的に与えられる。ここで、Ｌｖは狭義下側三角行列を示し、Ｕｖは狭義上側三角行列を示し、Ｄｖは共通の対角行列を示す。この場合には、式（４）は、下記式で表さ
れる。
Ｃｐｐ
＝Ａｐｐ−Ａｐｖ（Ｄｖ＋Ｕｖ）∧（−１）Ｄｖ（Ｌｖ＋Ｄｖ）∧（−１）Ａｖｐ・・・（１１）
ここで、式（１１）の中の因数Ａｐｖ（Ｄｖ＋Ｕｖ）∧（−１），（Ｌｖ＋Ｄｖ）∧（−１）Ａｖｐを、
Ｓｐｖ＝Ａｐｖ（Ｄｖ＋Ｕｖ）∧（−１）
Ｓｖｐ＝（Ｌｖ＋Ｄｖ）∧（−１）Ａｖｐ
とおけば、因数Ｓｐｖ，Ｓｖｐは、
Ｓｐｖ（Ｄｖ＋Ｕｖ）＝Ａｐｖ・・・（１２）
（Ｌｖ＋Ｄｖ）Ｓｖｐ＝Ａｖｐ・・・（１３）
で表される。このような未知行列Ｓｐｖ，Ｓｖｐの近似解行列Ｔｐｖ，Ｔｖｐは、不完全ＬＵ分解により求められる。式（１０）は、ＴｐｖとＴｖｐにより下記式：
Ｑｐｐ＝Ａｐｐ−ＴｐｖＤｖＴｖｐ・・・（１４）
で表され、このＱｐｐはＡｖｖの異方性が考慮されてＣｐｐの近似になっている。ＴｐｖとＴｖｐを求める際に、ＡｖｐとＡｐｖの非零成分以外へのフィルインが無視され、式（１４）のＱｐｐは、疎行列である点で式（８）に同じである。

本発明は、詳細に記述されているように、与えられた拘束条件付き方程式を離散化した際に、又は、流体解析、構造解析、電磁場解析、それらの相互作用を扱う連成解析に混合型有限要素離散化法を適用する際に生じる大規模疎行列連立一次方程式を前処理付き反復法で解く際に、その前処理演算において、全体系の方程式を主方程式と拘束方程式に分離し、主方程式に対応する主未知数を近似的に消去し、残りの拘束方程式をそれに対応する拘束的未知数のみからなる主未知数消去後拘束方程式に変換し、拘束的未知数のみを未知数とする変換後の主未知数消去後拘束方程式をその係数行列を陽に保持することなしに内部の前処理付き反復解法により近似的に解くことにより、その近似解を求め、その近似解を主方程式に代入して、主未知数の近似解を計算するごとに、全体系に対する外側の反復過程の中の前処理演算における主未知数と拘束的未知数の近似解を求めている。

既述の前処理つき反復解法の実施例：
解くべき連立一次方程式は、係数行列をＡ、右辺ベクトルをＢ、未知数ベクトルをＸとして、以下のように表されているものとする。
ＡＸ＝Ｂ・・・（イ）
更に、Ａを近似し、与えられた右辺ベクトルＦに対して、Ｖを未知数ベクトルとする連立一次方程式：
ＭＶ＝Ｆ・・・（ロ）
を容易に解くことができる行列Ｍが与えられている。以下では、Richardson反復法に前処理を適用した実施例を示す。Richardson反復では、単位行列Ｉを用いて、係数行列ＡをＡ＝Ｉ＋（Ａ−Ｉ）と分離し、以下の漸化式に基づき反復解Ｘ｛ｋ｝を更新する。
Ｘ｛ｋ＋１｝＋（Ａ−Ｉ）Ｘ｛ｋ｝＝Ｂ・・・（ハ）
式（ハ）の近似更新は、次式で示されるアルゴリズムにより実行される。
Ｆ｛ｋ｝＝Ｂ−ＡＸ｛ｋ｝
Ｘ｛ｋ＋１｝＝Ｘ｛ｋ｝＋Ｆ｛ｋ｝
Richardson反復法に近似行列Ｍにより前処理を適用する際は、もとの方程式（イ）を
ＡＭ∧（−１）（ＭＸ）＝Ｂ
と書き換え、その係数行列をＡＭ∧（−１）、未知数ベクトルをＭＸとみなして、Richardson反復法を適用する。
ＭＸ｛ｋ＋１｝＋（ＡＭ∧｛−１｝−Ｉ）ＭＸ｛ｋ｝＝Ｂ
これは、下記の反復アルゴリズムで計算される。
Ｆ｛ｋ｝＝Ｂ−ＡＸ｛ｋ｝
Ｖ｛ｋ｝＝Ｍ∧｛−１｝Ｆ｛ｋ｝
Ｘ｛ｋ＋１｝＝Ｘ｛ｋ｝＋Ｖ｛ｋ｝
このような反復計算で解が収束すれば（Ｘ｛ｋ＋１｝−Ｘ｛ｋ｝が十分に小さくなれば）、式（ハ）は、
Ｘ｛ｋ＋１｝−Ｘ｛ｋ｝＋ＡＸ｛ｋ｝＝Ｂ
であるので、もとの式（イ）にＸ｛ｋ｝を代入した式：
ＡＸ｛ｋ｝＝Ｂ
は、十分な精度で満足される。収束までの反復回数ｋはもとの係数行列Ａに対するＭの近似度に依存し、ＭがＡをよく近似しているほど小さくなる。

図１は、本発明による連立一次方程式反復求解計算機の実現態を示す回路ブロック図である。図２は、本発明による連立一次方程式反復求解計算機の実現態を示す計算フロー図である。図３は、本発明による連立一次方程式反復求解計算機の他の実現態を示す計算フロー図である。図４は、本発明による連立一次方程式反復求解計算機の更に他の実現態を示す計算フロー図である。本発明による連立一次方程式反復求解計算機の実現態の実施例を示す有限要素メッシュ図である。本発明による連立一次方程式反復求解計算機の実現態の他の実施例を示す有限要素メッシュ図である。

符号の説明

５…作成計算器、反復求解計算器（消去計算器、第１近似計算器、第２近似計算器）

Claims

現象を記述する主方程式と前記現象を拘束する拘束方程式とを分離し、且つ、前記主方程式と前記拘束方程式をそれぞれに多変数次元メッシュ上で離散化して下記の２式：
主方程式：ＡｖｖＸｖ＋ＡｖｐＸｐ＝Ｂｖ
拘束方程式：ＡｐｖＸｖ＋ＡｐｐＸｐ＝Ｂｐ
Ａｖｖ，Ａｖｐ，Ａｐｖ，Ａｐｐ：係数行列
Ｂｖ，Ｂｐ：既知ベクトル
Ｘｖ，Ｘｐ：求解対象未知ベクトル
として連立一次方程式を作成する作成計算器と、
前記連立一次方程式を前処理つき反復解法で解く反復求解計算器とを構成し、
前記反復求解計算器は、
前記前処理つき反復解法の反復過程の中の前処理演算で前記２式の右辺ベクトルＢｖおよびＢｐの代わりに与えられる右辺ベクトルＦ，Ｇに対して主未知数ベクトルＸｖを近似的に消去し主未知数消去後の拘束方程式の右辺ベクトルＱを作成する消去計算器と、
前記２式から主未知数ベクトルＸｖが近似的に消去され拘束的未知数ベクトルＸｐが未知数として残る前記主未知数消去後の前記拘束方程式を内側で前記前処理つき反復解法を実行することにより前記拘束的未知数ベクトルの近似解Ｐを求める第１近似計算器と、
前記近似解Ｐを前記主方程式に代入して前記主未知数ベクトルＸｖの近似解Ｖを求める第２近似計算器とを形成する
連立一次方程式反復求解計算機。
前記消去計算器は、前記主方程式の係数行列Ａｖｖを逆行列演算を容易に行うことができる近似行列Ｍｖｖで置き換えた下記の２式：
ＭｖｖＶ＋ＡｖｐＰ＝Ｆ
ＡｐｖＶ＋ＡｐｐＰ＝Ｇ
Ｆ，Ｇ：全体系に対する前処理つき反復解法の反復過程の中で前処理演算の右辺ベクトルとして作成される既知ベクトル
のうちの第１式の初期近似解を下記式：
ＭｖｖＶ’＝Ｆ
の解としてＶ’を求め、前記主未知数消去後の前記拘束方程式の右辺ベクトルＱを下記式：
Ｑ＝Ｇ−ＡｐｖＶ’
より作成する消去計算器部分を形成する
請求項１の連立一次方程式反復求解計算機。
前記第１近似計算器は、下記式：
Ｃｐｐ＝Ａｐｐ−ＡｐｖＭｖｖ∧（−１）Ａｖｐ
で定義される係数行列により定義される前記主未知数消去後の拘束方程式：
ＣｐｐＰ＝Ｑ
を近似的に解くことにより前記拘束的未知数ベクトルの近似解Ｐを求める
請求項２の連立一次方程式反復求解計算機。
前記第２近似計算器は、前記主未知数消去後の前記拘束方程式の近似解Ｐに対して、次式：
ＭｖｖＵ＝−ＡｐｖＰ
を解き、前記第１式の初期近似解Ｖ’を式：
Ｖ＝Ｖ’＋Ｕ
により更新して主未知数ベクトルの近似解Ｖを求める
請求項２または３の連立一次方程式反復求解計算機。
前記主未知数消去後の前記拘束方程式の係数行列Ｃｐｐに対して前記Ｍｖｖの代わりに対角行列Ｄを代入し、前記Ｃｐｐを近似する疎行列Ｑｐｐ：
Ｑｐｐ＝Ａｐｐ−ＡｐｖＤ∧（−１）Ａｖｐ
の不完全ＬＵ分解行列を前処理行列として前記前処理つき反復解法により前記主未知数消去後の前記拘束方程式が解かれる
請求項５の連立一次方程式反復求解計算機。
前記主未知数消去後の前記拘束方程式の係数行列Ｃｐｐに対して、
前記Ｃｐｐを近似する疎行列Ｑｐｐ：
Ｑｐｐ＝Ａｐｐ−ＴｐｖＤｖ’∧（−１）Ｔｖｐ
をＡｖｖの不完全ＬＵ分解行列：
Ｍｖｖ＝（Ｌｖ＋Ｄｖ）Ｄｖ∧（−１）（Ｄｖ＋Ｕｖ）
Ｕｖ：狭義上三角行列
Ｌｖ：狭義下三角行列
Ｄｖ：対角行列
に対して次式：
Ｔｐｖ（Ｄｖ＋Ｕｖ）＝Ａｐｖ
（Ｌｖ＋Ｄｖ）Ｔｖｐ＝Ａｖｐ
を近似的に満たす疎行列ＴｐｖとＴｖｐを用いることにより作成し、前記Ｑｐｐの不完全ＬＵ分解行列を前処理行列として前記前処理つき反復解法により前記主未知数消去後の前記拘束方程式が解かれる
請求項３の連立一次方程式反復求解計算機。
前記現象は物理現象である
請求項１〜６から選択される１請求項の連立一次方程式反復求解計算機。
前記現象は経済学現象である
請求項１〜６から選択される１請求項の連立一次方程式反復求解計算機。
主方程式を離散化する主離散化式と拘束方程式を離散化する拘束離散化式を計算器の中に記述して全体系を記述するステップと、
前記主離散化式に対応する主未知数のうちの主未知数を近似的に前記計算器の中で消去することにより前記拘束離散化式を前記拘束離散化式に対応する拘束的未知数で表される主未知数消去後拘束離散化式で前記計算器の中で記述するステップと、
前記主未知数消去後拘束離散化式の係数行列を疎行列化することにより疎行列化主未知数消去後拘束離散化方程式を前記計算器の中で記述するステップと、
前記主未知数消去後拘束離散化式を前処理つき反復解法により解くときの前処理演算として前記疎行列化主未知数消去後拘束離散化式を前記計算器の中で記述するステップとを構成する
連立一次方程式反復求解計算方法。