パソコン上で大規模メッシュを張った有限要素法(FEM)モデルの構造解析ができるようになったとはいえ、ある限界を超えるとメモリネックとなってしまって解析不能に陥ることになる。連立方程式の解法に、伝統的に使用されてきたガウスの消去法を代表例とする直接法を選択した場合の多くがそうである。ここで俄然、存在価値を高めてくるのが反復法である。もちろん、反復法でもメモリ限界があるが、直接法に比べると格段に有利である。

一口に反復法と言っても、何種類もあるが、FEMで使用されているのは共役勾配法(Conjugate Gradient Method：CGM)と思ってまず間違いない。この手法そのものは昔から存在していたが、最近のFEMソルバーで、CGMがそのまま使用されることはほとんどない。解を得るまでの反復回数を減少させて、計算時間の短縮化を狙った前処理がCGMに入る前に施されるのである。前処理が施されたCGMを特に前処理付き共役勾配法(Preconditioned Conjugate Gradient Method：PCGM)と呼んでいる。

前処理方法は過去にいろいろ提案されており、今も応用数学者の研究対象になっているテーマである。多くのバリエーションの中でも、現在の一番人気と言えば、元の係数マトリックス(すなわち剛性マトリックス)を一旦、不完全にコレスキー分解してしまうアルゴリズムである。

コレスキー分解というのは、直接法での解法の一つとして有名なコレスキー法による三角分解である(本エッセイ第51話参照)。もちろん、この三角分解を完全に実施してしまうと、連立方程式を解く行動に出ることになり本末転倒となってしまう。不完全というのがミソなのである。この不完全コレスキー分解付きCGMを俗にICCGM(Incomplete Cholesky Conjugate Gradient Method)と呼んでいる。PCGMにICCGMとやたらと略語が出てきてややこしいが、FEMユーザーは、ソルバーで使用されている反復法はPCGMだと知っておけばいいのではないかと思う。

CGMからPCGMへの開発の流れは、もともと計算時間の短縮化に目的があった。直接法では、計算時間短縮化のアルゴリズムが限界に達していたのに対し、反復法では1回でも反復回数を減らせば、それだけ計算時間が短縮でき、現在の歴史段階でも、まだ減らす余地があるからこそ、研究が続けられているのである。極論をすれば、係数マトリックスを単位マトリックスに変換してしまえば、1回の反復で終了してしまう。しかし、それをするための前処理に時間が掛かり過ぎ、意味がなくなってしまう。PCGMは、前処理と解式反復という二つの処理時間の妥協点を見つけるのが大テーマなのである。

ところで、PCGMといえども、その性質はCGMのそれがそのまま反映されているものである。それゆえ、以下は、反復法＝CGMとみなして、ここでは三つの大きな性質を紹介してみる。

 

まずは、冒頭で言ったメモリ有利の件である。本エッセイ第36話でも言ったように、通常、FEMでの剛性マトリックスは非零値が対角付近に集まるバンドマトリックスとなる。実際問題では各列の高さ(もしくは各行の長さ)が違うため、図61-1のようなきれいな一定バンド幅ではないが、話を簡単にするため、ここでは、これを仮定する。

さて、図61-1のような剛性マトリックスの場合、これを直接法で解こうとすると、零値、非零値にかかわらず、バンド内の全ての要素に対応するメモリが必要となる。なんとなれば、最初、ゼロであったマトリックス要素も、計算過程で非零になる可能性があるからである。これを、“フィルイン”と呼んでいる。これに対し、反復法では、反復過程に入った係数マトリックスの中身は変化しないので、最初にゼロであったマトリックス要素はいつまでもゼロであるため、それ用のメモリが不要となるのである。

この差は実に大きい。ちょっと、具体例を計算してみよう。いま、50万自由度の剛性マトリックスを想定してみる。構造解析での経験則で言えば、バンド幅(対称マトリックスではハーフバンド幅)は、全自由度の約1割前後程度となることが多い。しかし、ここでは、もっと少なめに仮定して1万自由度としてみる。そうすると、このバンド内のマトリックス要素用のメモリとして、50万×1万=5×109ワードすなわち、5ギガワードのメモリが必要となる。FEMでは、実数計算は倍精度計算されるので、バイト単位に換算すると、実に40ギガバイトのメモリが必要となる。現在（2009年当時）のスタンダード・パソコンである32ビットマシンでは、2ギガバイトの壁があるため、全く歯が立たないことになる。

ところで、構造解析のバンドマトリックスでは、実はバンド内のマトリックス要素もほとんどがゼロなのである。実際にモニタリングしてみれば分かることだが、マトリックスが大きくなればなるほど、マトリックス要素の99％以上がゼロなのである。すると、反復法では、この問題、400メガバイトのメモリで済むことになるというわけである。

 

2番目は反復過程での収束状況のことである。これこそ反復法の性格を特徴付けているもので、実に気まぐれなのである。いや、荒くれた性格と言う方が的を射ているかもしれない。収束終了するまでの残差ノルム値(簡単に言えば、検算時の誤差指標)の履歴を見ていると、まるで試行錯誤を繰り返しているような振る舞いなのである。単調減少の現象とは程遠いありさまを示してくれる。

図61-2を見ていただきたい。これは、10年ほど前（2009年当時）、筆者が出くわした悪夢のような体験を、記憶をたどりながら記したものである。もちろん、数字はその時のものと同じではないが、反復過程の様子に変わりはないものである。この計算は、夕方から開始したもので、残差ノルム値の状況からして、もうそろそろ収束するだろうと推測して、終了を楽しみに9時ぐらいからずっと残差ノルム値をモニタリングしていた。たしか夜中の11時半を過ぎてからだと思うが、最新の残差ノルム値が出力された。それが、何と、数値が減少するどころか、大幅に上昇してしまっているではないか。反復法に全く、裏切られたのである。あきれかえった私は、そそくさと家路についた。反復法の使用に慣れてないユーザーでは、その気まぐれな性格にあきれることが、しばしばであろう。慣れてくるにつれて、つれなく付き合えばいい、というコツを覚えるであろう－結局は収束してくれるので。

ちなみに、このときの計算は、朝方の時間に見事、収束していた。CPUタイムが11時間の計算であった。ハード面もソフト面も進歩した、現在のパソコンで－私の手元のある2.27GHzノート型パソコン－同じ問題を計算してみたら、9分ぐらいで完了してしまった。

 

3番目は、ちょっと面白い性格の話である。FEMに限らず、古くから構造解析に携わっている者は、おそらく、連立方程式の解法では直接法になじんでいる人がほとんどであろう。そこで、構造物の幾何寸法に計算時間が依存するという、反復法が示す性格が全く奇妙に映るのである。図61-3を見ていただきたい。左は辺の長さが1の立方体を片持ち支持した充実構造物の上面に分布荷重が載荷された問題である。右の構造物は、左の立方体の紙面横方向辺を3倍の長さに延長したものである。なお、ここでの話題では、幾何学的単位、また、材料定数、荷重強度の値、単位は全く本質的でないので、度外視する。

上の構造モデルに対して、60×60×60という計216,000の六面体要素をメッシュ分割してみる。このとき、節点数は226,981であり、有効自由度数が669,780となる。この二つのモデルに反復法を適用した結果、収束まで要する反復回数が次のとおりであった。

左の立方体モデル ： 106回

右の直方体モデル ： 162回

二つのモデルは、完全にメッシュの位相状態は一致しているので、もちろん計算時間は、直方体モデルの方が余計にかかることになる。幾何寸法が変わるだけで、メッシュの位相状態が一致し、総自由度数が一致している二つのモデルで、計算時間が変わるというのは、直接法の世界では考えられない反復法の特質である。

上の例でみたように、反復法はマッシブな構造に比較してスレンダー系の構造で計算上不利となる大きな特徴を持つ。したがって、塔や橋梁のような細長い構造物で反復法を利用する場合、節点数が多いFEMモデルでは相応の計算時間を覚悟しなければならない。

ここで、構造が細長くなればなるほど、反復法の計算時間のかかる理由として、反復法のファンタジーと呼べるものを紹介してみよう。図61-4のように構造物に荷重が載ったとき、弾性波が発生するとする。その波が四方八方に伝達し、やがて構造物の境界に達すると反射波となって帰ってくる。それら全ての波が干渉しあって減衰した時が、反復の終了になると理解すれば、波の発生源から境界位置までの距離が遠くなるスレンダー系構造物の不利が納得できるというものである。

ところで、ファンタジーと言ったものの、全くのファンタジーでもないことが、数値解析分野で活躍された戸川隼人先生著の本で知らされる¹ 。この本の一章で、共役勾配法における隠れた数理が波動論で説明できる旨掲載されている。

最後に、やはり波動論で説明できる現象を一つ紹介して終わりとしよう。図61-3右の直方体モデルを対象に断面を図61-5のように矩形にくり抜いてみる。すなわち、ボックス桁の問題に反復法を適用してみる。

この結果が表61-1である。これを見ると、空洞部が大きくなって、節点数、または自由度数が減っていくのに反比例して、反復回数が増大していることが分かる。これも、直接法の観点では全く非常識な話である。

ところで、この現象は、載荷により発生した弾性波が断面内を横断方向に直進する際、空洞にぶつかり迂回させられた結果、遠距離を伝達していくからだと説明できる。それゆえ、空洞が大きいほど反復回数が増えていくのである。ただ、本モデルでは、dの増加に伴い節点数も減少するので、反復回数が増えたからといって計算時間まで増加していくわけではない。剛性マトリックスが小さくなったぶん、1回の反復計算での計算時間が短縮するからである。

しかし、逆に言えば、同一節点数の有限要素法モデルであれば、空洞部が大きくなるほど、計算時間で不利になることを物語っている。ここにも、コンクリート橋梁設計者にとって、反復法が不利に働く場面を見つけるのである。

そうは言っても、大メッシュモデルの橋梁問題を直接法で解くよりも、まだ反復法で解く方が、処理時間が速いことをわれわれは経験しているし、だいいち、メモリネックで直接法では解けないことも経験しているのである。

2009年梅雨の候記

1. シリーズ17新しい応用の数学 共役勾配法 戸川隼人 教育出版 1977
   この本は、共役勾配法を分かりやすく、丁寧に解説した好著である。 [↩]