(57)【要約】      （修正有）【課題】大規模な固有値を計算機で求める、処理速度や記憶容量で有利な多階層のモード合成法等の近似解法において、これまでに知られていた実空間とモード空間を交互に解くＰ行列処理の非現実的な処理時間を改善し、数百万自由度以上の問題に適用できるようにする。【解決手段】モード結合による多階層のモード合成法において、解析領域全体を木の形式で多階層に領域分割し、その領域において実空間とモード空間を完全に分離する。木に領域分割した場合の根元を上の階層、枝の方を下の階層と呼ぶ。まず、実空間を下の階層から上の階層に順番に固有値問題を解く。その階層毎の実空間の固有ベクトルを基底として、上の階層から下の階層にモード空間での固有値問題を解き、最終的に全体領域の固有値・固有ベクトルを求めるアルゴリズムが構成でき、そのアルゴリズムを実現するシステムを構築する。【選択図】図１
2009151735.tif 000002
【特許請求の範囲】
【請求項１】
モード結合による多階層のモード合成法に基づき固有値問題を解くために、実空間を下の階層から上の階層に順番に固有値問題を解き、その階層毎の実空間の固有ベクトルを基底として、上の階層から下の階層にモード空間での固有値問題を解き、最終的に全体領域の固有値・固有ベクトルを求める高速に処理する手段として機能させるための固有値計算プログラム。
【請求項２】
モード結合による多階層のモード合成法に基づき固有値問題を解くために、実空間を下の階層から上の階層に順番に固有値問題を解き、その階層毎の実空間の固有ベクトルを基底として、上の階層から下の階層にモード空間での固有値問題を解き、最終的に全体領域の固有値・固有ベクトルを求める高速に処理する手段として機能させるための固有値計算プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項３】
モード結合による多階層のモード合成法に基づき固有値問題を解くために、実空間を下の階層から上の階層に順番に固有値問題を解き、その階層毎の実空間の固有ベクトルを基底として、上の階層から下の階層にモード空間での固有値問題を解き、最終的に全体領域の固有値・固有ベクトルを求める高速に処理する手段として機能させるための固有値計算プログラムを実行するための装置。
【発明の詳細な説明】
【技術分野】
【０００１】  本発明は、実ｎ次元空間における一般化固有値問題Ｋｘ＝λＭｘを対象とする。すなわち、ｎ×ｎ行列ＫとＭが与えられたときに、実固有値λおよびｎ次元固有ベクトルｘを求める問題である。行列ＫおよびＭは、スパース行列（行列の多くの要素が０である行列）である場合に有効なアルゴリズムである。大規模な工学的な振動問題では、スパース率は９９．９％以上である。ここでは、ｎが数千万に対して、数千個の固有値を求める問題を対象とするものである。ここで述べる方法は、上記の一般化固有値問題すべてに適用可能な方法であるが、先に述べた規模の大規模な問題に対して特に有効な高速演算処理方法に関する。
【背景技術】
【０００２】  構造設計において、製品性能と環境対策やユーザ満足度の要求を満たすための設計技術として、最適な音響／構造連成を求めるシミュレーション技術が要求されている。現行の音響／構造連成システムでは、動的応答（変位、応力、最大音圧、音圧積分、モード共振ピークの分散など）を最適化する構造変更を求めることが困難な機能上の問題と、使い勝手の点で利用上の制約がある。その解決案としては、構造全体を最適化変更の対象として、固有値・モード振幅を最適に調整することが必須である。このために、汎用的な音響／構造連成最適化定式の導入、大規模な設計変数規模の最適設計計算アルゴリズム、大規模な固有値計算が必要となる。【０００３】このような背景のもと、大規模構造モデルの音響構造連成応答（変形、応力、音圧、最大音圧や音圧積分などの音響特性、モード共振ピーク応答の分散低減など）を最適化解析するシミュレーションのための固有値問題を解く計算機プログラムが必要となってきている。【０００４】
【非特許文献１】長松  昭男，大熊  政明；″部分構造合成法，″培風館（１９９１）
【発明の開示】
【発明が解決しようとする課題】
【０００５】  通常の多階層モード空間モード合成法のアルゴリズムはつぎの通りである。（１）最下位の階層における内部領域群の固有値問題の固有値・固有ベクトルを求める。（２）その親の静縮合された固有値問題の固有値・固有ベクトルを求める。（３）この２つの固有値問題の固有ベクトルを基底として、その空間（モード空間）で、固有値・固有ベクトルを求める。（４）モード空間で求まった固有ベクトルを実空間における固有ベクトルに変換する。（５）ここの固有値問題の解を改めて、内部領域の解と考え、さらにひとつ上の階層の固有値問題を解く。すなわち、その親の静縮合された固有値問題を解くステップ（２）に戻る。【０００６】  ここで、φｍＢ，ｊ（ｉ）＝−Ｋｍｉ−１Ｋｍｉ，ｊφｍＢ，ｊを求める必要があり、ＫｍｉはφｍＢ，ｊ以下の領域に関するＫ行列全体であり、Ｋｍｉ−１を計算するにはつぎのような手順で求められる。
Ｐｍｉ＝Ｋｍｉ−１−Ｋｍｉ−１
このＰ行列に多大な処理時間がかかるので高速演算処理法を提供することが本発明の課題である。
【課題を解決するための手段】
【０００７】  上記の課題を解決するために、本発明は、領域分割；ベクトルｘのｎ成分をノードとし、行列ＫおよびＭの非零成分をエッジとするグラフＤ１を考える。このグラフをｐ個に分離されたグラフＤ１ｉ（ｉ＝１，ｐ１）その境界のノードで構成されるグラフＤ１Ｂ，ｊ（ｊ＝１）を作る。すなわち、Ｄ１ｉとＤ１Ｂ，ｊのノード間にはエッジＤ１ｉ，ｊは存在するが、Ｄ１ｉとＤ１ｊ（ｉ≠ｊ）のノード間にエッジは存在しない。また、Ｄ１Ｂ，ｊ（ｊ＝１）とＤ１ｉ（ｉ＝１，ｐ１）およびＤ１Ｂ，ｊで、もとのグラフＤ１全体を構成する。分割されたＤ１ｉ（ｉ＝１，ｐ１）を子、その子とエッジを共有するＤ１Ｂ，ｊ（ｊ＝１）を親と呼ぶ。この操作を、Ｄ１ｉ（ｉ＝１，ｐ１）に対して再帰的に行うことにより、Ｄ２ｉ（ｉ＝１，ｐ２）とＤ２Ｂ，ｊ（ｊ＝１，ｐ１）を構成する。注意すべきことは、Ｄ２Ｂ，ｊの個数がＤ１ｉの個数となることである。さらに、Ｄｍｉ（ｉ＝１，ｐｍ）に対して、同様の操作を行えば、Ｄｍ＋１ｉ（ｉ＝１，ｐｍ＋１）とＤｍ＋１Ｂ，ｊ（ｊ＝１，ｐｍ）を構成できる。子と親の関係も同様に定義できる。また、ひとつの親ＤｍＢ，ｊの子たちをＤｍｉ（ｉ＝１，ｐｍ（ｊ））と書くことにする。【０００８】  階層；領域分割の手順により、つぎの階層のグラフ（以下、領域とも呼ぶ）が定義できる。
Ｄ１→Ｄ１Ｂ，ｊ（ｊ＝１）→Ｄ２Ｂ，ｊ（ｊ＝１，ｐ１）→・・・→ＤｍＢ，ｊ（ｊ＝１，ｐｍ−１）→Ｄｍｊ（ｊ＝１，ｐｍ）
ここで、Ｄ１Ｂ，ｊ（ｊ＝１）を階層１、Ｄ２Ｂ，ｊ（ｊ＝１，ｐ１）を階層２、ＤｍＢ，ｊ（ｊ＝１，ｐｍ−１）を階層ｍの（境界）領域と呼び、最後のＤｍｊ（ｊ＝１，ｐｍ）をｍ＋１階層の領域または内部領域と呼ぶ。【０００９】  縮合；ある内部領域または境界領域およびその子すべてを含む領域Ｄにおいて、Ｋｘ＝λＭｘを考える。ここでは、他の領域の成分は、０と考える。このような問題を、Ｋｘ＝λＭｘ；Ｄと書くことにする。この問題は、Ｄのノードの個数の次元ｄｉｍ（Ｄ）を持つ。ここで、
2009151735.tif 000003
した固有値問題と呼ぶ。これも内部領域と同様に、Ｋｘ＝λＭｘ；ＤｍＢ，ｊと書くものとする。これは、一般に知られているＧｙａｎの静縮合を含む概念である。【００１０】縮合の演算；ＫとＭに関する縮合から成る。まず、１階層を対象としたＫの縮合は、
2009151735.tif 000004
から
2009151735.tif 000005
を構成する手続きである。全体のＫの縮合は、この手続きを再帰的に実行する手順である。つぎに、１階層を対象としたＭの縮合は、
2009151735.tif 000006
から
2009151735.tif 000007
を構成する手続きである。全体のＭの縮合は、この手続きを再帰的に実行する手順である。すなわち、静縮合とは、
2009151735.tif 000008
を
2009151735.tif 000009
2009151735.tif 000010
の固有値問題に順次変換する再帰的な手順である。【００１１】
2009151735.tif 000011
縮合された固有値問題Ｋｘ＝λＭｘ；ＤｍＢ，ｊの固有ベクトル、および、（２）内部領域Ｋｘ＝λＭｘ；Ｄｍｉの固有ベクトルを基底とする関数空間（これらの固有ベクトルの線形結合で作られる空間）を考え、その空間の範囲で、
2009151735.tif 000012
を解く。これをモード空間の固有値問題と呼ぶ。すなわち基底空間の係数ηを利用して、任意の座標はｘ＝Ｔηと表現でき、Ｋｘ＝λＭｘを（ＴＴＫＴ）η＝λ（ＴＴＭＴ）ηと変換して、モード空間での固有値問題に帰着できる。Ｋｘ＝λＭｘ；Ｄｍｉの固有ベクトルを、φｍｉ（ｉ＝１，ｐｍ（ｉ））とおき、Ｋｘ＝λＭｘ；ＤｍＢ，ｊの固有ベクトルを、φｍＢ，ｊとおくと、Ｔ＝ＰφとしてＴ行列を求めることができる。Ｐ行列の定式化は、モード空間の基底をどのように選択するかに依存しており、その定式化アルゴリズムによって異なる。【００１２】  モード空間の演算；モード空間の基底については、つぎのような方法で求める。
Ｋｘ＝λＭｘ；Ｄｍｉの固有ベクトルを、φｍｉ（ｉ＝１，ｐｍ（ｉ））とおき、Ｋｘ＝λＭｘ；ＤｍＢ，ｊの固有ベクトルを、φｍＢ，ｊとおく。モード空間の基底は、
2009151735.tif 000013
である。ここで、最後の基底は、Ｇｙａｎの静縮合の条件
2009151735.tif 000014
を満たす。すなわち、
φｍＢ，ｊ（ｉ）＝−Ｋｍｉ−１Ｋｍｉ，ｊφｍＢ，ｊ
として決定できる。これらの基底を利用して
2009151735.tif 000015
とすると、Ｋｘ＝λＭｘは、ｘ＝Ｔηを利用して、
（ＴＴＫＴ）η＝λ（ＴＴＭＴ）η
となって、モード空間での固有値問題に帰着できる。【００１３】  高速な多階層モード空間モード合成法のアルゴリズムはつぎの通りである。（１）まず、最下位の階層ｍにおける内部領域群の固有値問題Ｋｘ＝λＭｘ；Ｄｍｉで固有値・固有ベクトルを求める。（２）その親の静縮合された固有値問題Ｋｘ＝λＭｘ；ＤｍＢ，ｊで固有値・固有ベクトルを求める。ｍをひとつずつ小さくして、境界領域ＤｌＢ，ｊ（ｌ＝１，ｍ）および最下位の内部領域Ｄｍｊ（ｌ＝１，ｐｍ−１）における固有問題を解き、固有値・固有ベクトルを求める。（３）まず、ｍ＝１として、この２つの固有値問題の固有ベクトルを基底として、
2009151735.tif 000016
ける固有ベクトルに変換する。Ｔ＝Ｐφとして（ＴＴＫＴ）η＝λ（ＴＴＭＴ）ηに変換する場合のＰ行列の定式化は下記の通りである。（５）ここの固有値問題の解を改めて、境界領域の解と考え、さらにｍ＋１階層の固有値問題を解く。すなわち、その子を含めたモード空間での固有値問題を解くステップ（３）に戻る。【００１４】ここで、φｍＢ，ｊ（ｉ）＝−Ｋｍｉ−１Ｋｍｉ，ｊφｍＢ，ｊを求める必要がある。ＫｍｉはφｍＢ，ｊ以下の領域に関するＫ行列全体であり、Ｋｍｉ−１を計算するにはつぎのような手順で求められる。
Ｐｍｉ＝Ｋｍｉ−１−Ｋｍｉ−１
である。具体的には、
2009151735.tif 000017
を解くことになる。ここで、Ｐ行列の演算は、逆行列がすでに求められているため、高速に処理することができる。【００１５】  一般化固有値問題から標準固有値問題への変換；  モード空間では、一般化固有値問題Ｋｘ＝λＭｘを解く必要があるが、
2009151735.tif 000018
の行列になっている。したがって、上記の一般化固有値問題を等価な標準固有値問題へ変換すると、
2009151735.tif 000019
ができる。ここで、Ｋ−１／２ＭＫ−１／２は、Ｍ＝（ｍｉｊ）とすると、
2009151735.tif 000020
したがって、これを変換して、もとの一般化固有値問題の固有値λと固有ベクトルＫ−１／２（Ｋ１／２ｘ）を求める。ただし、ここで、
2009151735.tif 000021
とおいた。【００１６】
2009151735.tif 000022
っている。【００１７】  シフト可能な一般化固有値問題への再変換；モード空間での一般化固有値問題Ｋｘ＝λＭｘ
2009151735.tif 000023
に対しては、大きい固有値から計算できるため、この方法は都合がよかった。ただし、小
2009151735.tif 000024
有値の数が多い場合にはシフトして計算が必要であるが、（１）ＢＬＺＰＡＣＫの標準固有値問題では、固有値範囲を指定できない。（２）Ａｘ＝λｘのままシフトして解くとマイナスの固有値も計算される。（３）上記の理由からＳＬＩＣＥオプションが有効でない。の３つの理由から、多くの固有値を求めると、固有値が求まらないか、または、多くの処理時間を要することになる。したがって、標準固有値問題の解きやすいという性質を保存しつつ、つぎの一般化固有値として、ｓｌｉｃｅオプションを利用できる形式で利用する。
2009151735.tif 000025
Ａ＝Ｋ−１／２ＭＫ−１／２
ｘ′＝Ｋ１／２ｘ
として、Ａ−１ｘ′＝λｘ′として、シフトした系
（Ａ−１−λ０Ｉ）ｘ′＝（λ−λ０）ｘ′
を計算する。Ｌａｎｃｚｏｓ法の過程では、上式の左辺の行列の逆を、与えられたＵに対して
（Ａ−１−λ０Ｉ）Ｖ＝Ｕ
を満たすＶを計算する必要がある。Ａのスパース性を利用するため、
Ｖ＝（Ｉ−λ０Ａ）−１ＡＵ
で計算する。ここで、（Ｉ−λ０Ａ）−１はＡのＭ行列と同じスパース性を保ったまま、ｓｐｏｏｌｅｓで計算できる。
【発明の効果】
【００１８】    式φｍＢ，ｊ（ｉ）＝−Ｋｍｉ−１Ｋｍｉ，ｊφｍＢ，ｊにおいて、−Ｋｍｉ−１は、対角ブロック行列であり、縮合過程ですでに計算されているため、本アルゴリズムでは、高速に処理できる。【００１９】  一般化固有値問題を標準固有値問題に変換したことで短時間に固有値を求めることができ、さらに、その標準固有値問題をその固有値を高速に、かつ、安定に求めることができる。
【発明を実施するための最良の形態】
【００２０】  図１は本発明のデータ処理の実施手順を示したものである。【００２１】  まず、一般化固有値問題Ｋｘ＝λＭｘを解くための、ｎ×ｎ行列ＫとＭのデータを読み込み、領域分割を実行し、Ｄ１→Ｄ１Ｂ，ｊ（ｊ＝１）→Ｄ２Ｂ，ｊ（ｊ＝１，ｐ１）→・・・→ＤｍＢ，ｊ（ｊ＝１，ｐｍ−１）→Ｄｍｊ（ｊ＝１，ｐｍ）のデータを作成する。【００２２】  内部領域Ｄｍｊ（ｊ＝１，ｐｍ）の固有値を計算する。【００２３】
2009151735.tif 000026
らＤｍＢ，ｊに静縮合した固有値問題を解く。ｍを１だけ減じながら、ｍが１になるまで処理を繰り返す。【００２４】  Ｄ１Ｂ，ｊ（ｊ＝１）におけるモード空間での固有値を求める。つぎにＤ２Ｂ，ｊ（ｊ＝１，ｐ１）を求める。以下、ｍを１だけ増加させながら、Ｄｍｊ（ｊ＝１，ｐｍ）のモード空間での固有値を求める。
【図面の簡単な説明】
【００２５】
【図１】本発明のモード空間による多階層のモード合成法のフローチャートである。
【図２】非特許文献１の公知のモード空間による多階層のモード合成法のフローチャートである。
【図１】
2009151735.tif 000027
【図２】
2009151735.tif 000028