(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2023142272
(43)【公開日】2023-10-05
(54)【発明の名称】計算プログラム及び計算方法
(51)【国際特許分類】
G06F 17/12 20060101AFI20230928BHJP
【ＦＩ】
G06F17/12
【審査請求】未請求
【請求項の数】5
【出願形態】ＯＬ
(21)【出願番号】P 2022049092
(22)【出願日】2022-03-24
(71)【出願人】
【識別番号】000005223
【氏名又は名称】富士通株式会社
(74)【代理人】
【識別番号】100121083
【弁理士】
【氏名又は名称】青木 宏義
(74)【代理人】
【識別番号】100138391
【弁理士】
【氏名又は名称】天田 昌行
(74)【代理人】
【識別番号】100074099
【弁理士】
【氏名又は名称】大菅 義之
(72)【発明者】
【氏名】大山 洋介
【テーマコード（参考）】
5B056
【Ｆターム（参考）】
5B056BB02
(57)【要約】
【課題】並列に動作する複数の処理単位を用いた反復法の計算時間を短縮する。
【解決手段】コンピュータは、１つ又は複数のループ処理各々において、並列に動作する複数の処理単位を用いて解の更新を反復する反復法の計算を実行する。コンピュータは、１つ又は複数のループ処理よりも後の所定のループ処理において、複数の処理単位を用いて反復法の計算を実行する。コンピュータは、１つ又は複数のループ処理各々の反復法の計算において解が更新された回数に基づき、所定のループ処理の反復法の計算において更新終了を判定する判定処理のタイミングを決定する。判定処理は、複数の処理単位の間における通信の結果に基づいて、更新終了を判定する処理を含む。
【選択図】図１７
【解決手段】コンピュータは、１つ又は複数のループ処理各々において、並列に動作する複数の処理単位を用いて解の更新を反復する反復法の計算を実行する。コンピュータは、１つ又は複数のループ処理よりも後の所定のループ処理において、複数の処理単位を用いて反復法の計算を実行する。コンピュータは、１つ又は複数のループ処理各々の反復法の計算において解が更新された回数に基づき、所定のループ処理の反復法の計算において更新終了を判定する判定処理のタイミングを決定する。判定処理は、複数の処理単位の間における通信の結果に基づいて、更新終了を判定する処理を含む。
【選択図】図１７
【特許請求の範囲】
【請求項1】
コンピュータのための計算プログラムであって、
前記計算プログラムは、
１つ又は複数のループ処理各々において、並列に動作する複数の処理単位を用いて解の更新を反復する反復法の計算を実行し、
前記１つ又は複数のループ処理よりも後の所定のループ処理において、前記複数の処理単位を用いて前記反復法の計算を実行し、
前記１つ又は複数のループ処理各々の前記反復法の計算において前記解が更新された回数に基づき、前記所定のループ処理の前記反復法の計算において更新終了を判定する判定処理のタイミングを決定する、
処理を前記コンピュータに実行させ、
前記判定処理は、前記複数の処理単位の間における通信の結果に基づいて、前記更新終了を判定する処理を含むことを特徴とする計算プログラム。
前記計算プログラムは、
１つ又は複数のループ処理各々において、並列に動作する複数の処理単位を用いて解の更新を反復する反復法の計算を実行し、
前記１つ又は複数のループ処理よりも後の所定のループ処理において、前記複数の処理単位を用いて前記反復法の計算を実行し、
前記１つ又は複数のループ処理各々の前記反復法の計算において前記解が更新された回数に基づき、前記所定のループ処理の前記反復法の計算において更新終了を判定する判定処理のタイミングを決定する、
処理を前記コンピュータに実行させ、
前記判定処理は、前記複数の処理単位の間における通信の結果に基づいて、前記更新終了を判定する処理を含むことを特徴とする計算プログラム。
【請求項2】
前記判定処理のタイミングを決定する処理は、前記１つ又は複数のループ処理各々の前記反復法の計算において前記解が更新された回数の統計値に基づき、前記所定のループ処理の前記反復法の計算において前記判定処理を開始する開始タイミングを決定する処理を含むことを特徴とする請求項１記載の計算プログラム。
【請求項3】
前記判定処理のタイミングを決定する処理は、前記所定のループ処理の前記反復法の計算において前記解の更新が複数回行われる間に、前記判定処理が１回行われるように、前記判定処理のタイミングを決定する処理を含むことを特徴とする請求項１又は２記載の計算プログラム。
【請求項4】
前記反復法の計算は、連立一次方程式の解を求める計算であり、
前記複数の処理単位各々には、前記連立一次方程式の解の一部を計算する処理が割り当てられ、
前記通信の結果に基づいて前記更新終了を判定する処理は、
前記複数の処理単位の間で、前記連立一次方程式の解の一部に対する残差に関する情報を送受信する処理と、
前記連立一次方程式の解の一部に対する残差に関する情報を用いて、前記連立一次方程式の解に対する残差を求める処理と、
前記連立一次方程式の解に対する残差に基づいて、前記解の更新を終了するか否かを判定する処理と、
を含むことを特徴とする請求項１又は２記載の計算プログラム。
前記複数の処理単位各々には、前記連立一次方程式の解の一部を計算する処理が割り当てられ、
前記通信の結果に基づいて前記更新終了を判定する処理は、
前記複数の処理単位の間で、前記連立一次方程式の解の一部に対する残差に関する情報を送受信する処理と、
前記連立一次方程式の解の一部に対する残差に関する情報を用いて、前記連立一次方程式の解に対する残差を求める処理と、
前記連立一次方程式の解に対する残差に基づいて、前記解の更新を終了するか否かを判定する処理と、
を含むことを特徴とする請求項１又は２記載の計算プログラム。
【請求項5】
１つ又は複数のループ処理各々において、並列に動作する複数の処理単位を用いて解の更新を反復する反復法の計算を実行し、
前記１つ又は複数のループ処理よりも後の所定のループ処理において、前記複数の処理単位を用いて前記反復法の計算を実行し、
前記１つ又は複数のループ処理各々の前記反復法の計算において前記解が更新された回数に基づき、前記所定のループ処理の前記反復法の計算において更新終了を判定する判定処理のタイミングを決定する、
処理をコンピュータが実行し、
前記判定処理は、前記複数の処理単位の間における通信の結果に基づいて、前記更新終了を判定する処理を含むことを特徴とする計算方法。
前記１つ又は複数のループ処理よりも後の所定のループ処理において、前記複数の処理単位を用いて前記反復法の計算を実行し、
前記１つ又は複数のループ処理各々の前記反復法の計算において前記解が更新された回数に基づき、前記所定のループ処理の前記反復法の計算において更新終了を判定する判定処理のタイミングを決定する、
処理をコンピュータが実行し、
前記判定処理は、前記複数の処理単位の間における通信の結果に基づいて、前記更新終了を判定する処理を含むことを特徴とする計算方法。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、計算技術に関する。
【背景技術】
【0002】
対称連立一次方程式の解を求める反復法として、共役勾配法（Conjugate Gradient Method，ＣＧ法）及び前処理付き共役勾配法（Preconditioned CG Method，ＰＣＧ法）が知られている（例えば、非特許文献１を参照）。非対称連立一次方程式の解を求める反復法として、双共役勾配法（Biconjugate Gradient Method，ＢｉＣＧ法）及び前処理付き双共役勾配法（Preconditioned BiCG Method，ＰＢｉＣＧ法）も知られている（例えば、非特許文献２を参照）。
【0003】
数値解析演算の収束判定、収束までの計算回数等の情報を計算中に取得することができる解析装置も知られている（例えば、特許文献１を参照）。反復計算において、計算完了までの計算量をより精度よく推定する反復計算量推定装置も知られている（例えば、特許文献２を参照）。関数を近似する方法も知られている（例えば、特許文献３を参照）。
【先行技術文献】
【特許文献】
【0004】
【特許文献1】特開２００７－２８７０５５号公報
【特許文献2】特開２０１７－９７３９２号公報
【特許文献3】米国特許出願公開第２０１０／０２９３２１３号明細書
【非特許文献】
【0005】
【非特許文献1】E. F. Kaasschieter, “Preconditioned conjugate gradients for solving singular systems”, Journal of Computational and Applied mathematics 24, 1988, pages 265－275.
【非特許文献2】“Biconjugate Gradient Method”、Wolfram MathWorld、September 19, 2021、［online］、［令和３年１０月１日検索］、インターネット＜ＵＲＬ：https://mathworld.wolfram.com/BiconjugateGradientMethod.html＞
【発明の概要】
【発明が解決しようとする課題】
【0006】
連立一次方程式の解を求める反復法の計算を、複数のプロセスを用いて並列化する場合、プロセス間通信にかかる時間がボトルネックになることがある。
【0007】
なお、かかる問題は、反復法の計算を複数のプロセスを用いて並列化する場合に限らず、反復法の計算を様々な処理単位を用いて並列化する場合において生ずるものである。
【0008】
１つの側面において、本発明は、並列に動作する複数の処理単位を用いた反復法の計算時間を短縮することを目的とする。
【課題を解決するための手段】
【0009】
１つの案では、計算プログラムは、以下の処理をコンピュータに実行させる。
【0010】
コンピュータは、１つ又は複数のループ処理各々において、並列に動作する複数の処理単位を用いて解の更新を反復する反復法の計算を実行する。コンピュータは、１つ又は複数のループ処理よりも後の所定のループ処理において、複数の処理単位を用いて反復法の計算を実行する。
【0011】
コンピュータは、１つ又は複数のループ処理各々の反復法の計算において解が更新された回数に基づき、所定のループ処理の反復法の計算において更新終了を判定する判定処理のタイミングを決定する。判定処理は、複数の処理単位の間における通信の結果に基づいて、更新終了を判定する処理を含む。
【発明の効果】
【0012】
１つの側面によれば、並列に動作する複数の処理単位を用いた反復法の計算時間を短縮することができる。
【図面の簡単な説明】
【0013】
【図1】ＣＧ法の疑似コードを示す図である。
【図2】ＰＣＧ法の疑似コードを示す図である。
【図3】ＰＢｉＣＧ法の疑似コードを示す図である。
【図4】並列化されたＣＧ法の第１の疑似コードを示す図である。
【図5】ｇＳｕｍＰｒｏｄの疑似コードを示す図である。
【図6】ｇＳｐＭＶの疑似コードを示す図である。
【図7】ｇＳｕｍＭａｇの疑似コードを示す図である。
【図8】並列化されたＣＧ法の第２の疑似コードを示す図である。
【図9】実施形態の情報処理装置の機能的構成図である。
【図10】第１の計算処理のフローチャートである。
【図11】情報処理装置の具体例を示す機能的構成図である。
【図12】判定処理の回数を調整可能なＣＧ法の疑似コードを示す図である。
【図13】ｍｉｎ（ｓｊ）の変化を示す図である。
【図14】ループ処理における第１の動作を示す図である。
【図15】ループ処理における第２の動作を示す図である。
【図16】ループ処理における第３の動作を示す図である。
【図17】残差ノルムの計算回数を示す図である。
【図18】第２の計算処理のフローチャートである。
【図19】解更新処理のフローチャートである。
【図20】履歴更新処理のフローチャートである。
【図21】解更新処理を示す図である。
【図22】複数のノードを含む情報処理装置のハードウェア構成図である。
【図23】流体解析シミュレーションにおける計算時間を示す図である。
【図24】流体解析シミュレーションにおけるδｉの計算回数及び解更新回数を示す図である。
【図25】情報処理装置のハードウェア構成図である。
【発明を実施するための形態】
【0014】
以下、図面を参照しながら、実施形態を詳細に説明する。
【0015】
まず、次式のような連立一次方程式の解法について説明する。
【0016】
Ａｘ＝ｂ （１）
【0017】
行列Ａは、係数行列を表し、ベクトルｂは、定数項を表す。ベクトルｘは、連立一次方程式の解を表す。Ａとしては、例えば、疎行列が用いられる。
【0018】
Ａが対称行列である場合、式（１）は対称連立一次方程式となり、Ａが非対称行列である場合、式（１）は非対称連立一次方程式となる。対称連立一次方程式の反復法としては、例えば、ＣＧ法又はＰＣＧ法が用いられ、非対称連立一次方程式の反復法としては、例えば、ＢｉＣＧ法又はＰＢｉＣＧ法が用いられる。
【0019】
図１は、ＣＧ法の疑似コードの例を示している。ベクトルｘｉは、反復法において更新される解を表し、ベクトルｘ０は、初期解を表す。ベクトルｒｉは、反復法において更新される残差を表し、ベクトルｒ０は、初期残差を表す。
【0020】
ｉｆ文１０１のｒｉ＝０は、ＣＧ法のｆｏｒループの終了条件を表す。実用上は、ｒｉのノルムが適当な閾値εよりも小さくなった場合に、反復法の計算が打ち切られ、最後に計算されたｘｉが解として出力される。
【0021】
図２は、非特許文献１に記載されたＰＣＧ法の疑似コードの例を示している。式２０１のＭは、前処理行列を表し、Ｍ－１は、Ｍの逆行列を表す。式２０１は、ＰＣＧ法において、Ｍ－１とｒｉの疎行列ベクトル積を計算する前処理を表す。ＰＣＧ法によれば、ｒｉに前処理を適用することで解の収束性が改善される。
【0022】
図３は、ＰＢｉＣＧ法の疑似コードの例を示している。式３０１は、図２の式２０１に対応する。式３０２のＭ－Ｔは、Ｍ－１の転置行列を表す。Ｍ－Ｔは、Ｍの転置行列ＭＴの逆行列に一致する。式３０１及び式３０２は、ＰＢｉＣＧ法における前処理を表す。
【0023】
行列Ａが対称行列である場合、ＰＢｉＣＧ法のアルゴリズムは、ＰＣＧ法のアルゴリズムに一致する。前処理を適用しない場合、Ｍ－１は単位行列となり、ＰＢｉＣＧ法のアルゴリズムは、ＢｉＣＧ法のアルゴリズムに一致する。
【0024】
複数のプロセスを用いて反復法の計算を並列化する場合、変数が複数の部分に分割され、各部分を用いた計算が各プロセスに割り当てられる。
【0025】
図４は、並列化されたＣＧ法の第１の疑似コードの例を示している。ベクトルｘ（ｐ）は、分割されたベクトルｘのうち、プロセスｐによって更新されるベクトルを表す。例えば、ｘの要素数がプロセス数で等分され、各部分の要素が順番に各プロセスｐに割り当てられる。ベクトルｒ（ｐ）は、分割されたベクトルｒのうち、プロセスｐによって更新されるベクトルを表す。ｘ（ｐ）は、連立一次方程式の解の一部の一例であり、ｒ（ｐ）は、解の一部に対する残差の一例である。
【0026】
式４０１は、プロセスｐに割り当てられたｒ０
（ｐ）を求める計算式である。式４０１のｇＳｐＭＶ（Ａ，ｘ０
（ｐ））は、Ａとｘ０の行列ベクトル積を表すベクトルの要素のうち、プロセスｐが使用する部分を表す。式４０３のｇＳｐＭＶ（Ａ，ｐｉ
（ｐ））は、Ａとｐｉの行列ベクトル積を表すベクトルの要素のうち、プロセスｐが使用する部分を表す。
【0027】
式４０２のｇＳｕｍＰｒｏｄ（ｒｉ
（ｐ），ｒｉ
（ｐ））は、ベクトルｒｉとベクトルｒｉの内積を表す。式４０４のｇＳｕｍＰｒｏｄ（ｐｉ
（ｐ），ｑｉ
（ｐ））は、ベクトルｐｉとベクトルｑｉの内積を表す。
【0028】
図５は、図４のｇＳｕｍＰｒｏｄの疑似コードの例を示している。図５の疑似コードは、分割されたベクトルｘ（ｐ）と分割されたベクトルｙ（ｐ）を用いて計算された、ベクトルｘとベクトルｙの内積（ｘ，ｙ）を出力する。
【0029】
式５０１のａｌｌｒｅｄｕｃｅ（ａ）は、ＭＰＩ（Message Passing Interface）のａｌｌｒｅｄｕｃｅに対応し、各プロセスｐによって計算されたａの値の総和を表す。ａｌｌｒｅｄｕｃｅ通信は、集団通信の一例である。集団通信は、プロセスのような複数の通信主体の間で行われる、一対多、多対一、又は多対多の通信である。
【0030】
【0031】
行集合Ｒ（ｐ）は、入出力を表す列ベクトルの要素のうち、プロセスｐが担当する要素の行番号の集合を表す。行列Ａ（ｐ）は、ＡのＲ（ｐ）×Ｒ（ｐ）成分を要素とする行列を表す。Ｒ（ｐ）×Ｒ（ｐ）成分は、Ａの要素のうち、Ｒ（ｐ）の各番号が示す行とＲ（ｐ）の各番号が示す列に対応するＡの要素を表す。
【0032】
Ａ（ｐ，ｑ）は、ＡのＲ（ｐ）×Ｒ（ｑ）成分のうち、非ゼロ要素からなるベクトルを表す。ｎｎｚ（ｐ，ｑ）は、Ｒ（ｐ）×Ｒ（ｑ）成分に含まれる非ゼロ要素の個数を表す。ｐ（ｐ，ｑ）は、プロセスｑが行う（Ａｘ）（ｑ）の計算において、Ａ（ｑ，ｐ）に含まれる非ゼロ要素の列位置に対応するように、ｘ（ｐ）を変換する写像を表す。式６０１のＳｐＭＶ（Ａ（ｐ），ｘ（ｐ））は、Ａ（ｐ）とｘ（ｐ）の行列ベクトル積を表す。
【0033】
一例として、Ａが次のような４行４列の疎行列であり、プロセス０及びプロセス１の２個のプロセスを用いてｙ（ｐ）を計算する場合について説明する。
【0034】
【数1】
【0035】
ａ、ｂ、ｃ、ｄ、ｅ、及びｆは、非ゼロ要素である。ｘ＝（ｘ０，ｘ１，ｘ２，ｘ３）Ｔ、Ｒ（０）＝｛０，１｝、Ｒ（１）＝｛２，３｝である場合、ｘ（０）＝（ｘ０，ｘ１）Ｔ、ｘ（１）＝（ｘ２，ｘ３）Ｔとなる。また、Ｒ（０）×Ｒ（０）＝｛（０，０），（０，１），（１，０），（１，１）｝となり、Ｒ（１）×Ｒ（１）＝｛（２，２），（２，３），（３，２），（３，３）｝となる。したがって、Ａ（０）及びＡ（１）は、次のような２行２列の行列となる。
【0036】
【数2】
【0037】
この場合、ｎｎｚ（０，１）＝２、ｎｎｚ（１，０）＝１、Ａ（０，１）＝（ｂ，ｄ）Ｔ、Ａ（１，０）＝（ｅ）Ｔである。また、ｐ（０，１）（ｘ（０））＝ｐ（０，１）（（ｘ０，ｘ１）Ｔ）＝（ｘ０）Ｔであり、ｐ（１，０）（ｘ（１））＝ｐ（１，０）（（ｘ２，ｘ３）Ｔ）＝（ｘ２，ｘ２）Ｔである。
【0038】
プロセス０は、ｐ（０，１）（ｘ（０））＝（ｘ０）Ｔをプロセス１へ送信し、次式により、ｙ（０）を計算する。
【0039】
ｙ（０）＝ＳｐＭＶ（Ａ（０），ｘ（０））
＝（ａｘ０，ｃｘ１）Ｔ （５）
＝（ａｘ０，ｃｘ１）Ｔ （５）
【0040】
次に、プロセス０は、ｚ＝ｐ（１，０）（ｘ（１））＝（ｘ２，ｘ２）Ｔをプロセス１から受信し、次式により、ｙ（０）を更新する。
【0041】
ｙ（０）＝ｙ（０）＋（Ａ（０，１））Ｔｚ
＝（ａｘ０＋ｂｘ２，ｃｘ１＋ｄｘ２）Ｔ （６）
＝（ａｘ０＋ｂｘ２，ｃｘ１＋ｄｘ２）Ｔ （６）
【0042】
一方、プロセス１は、ｐ（１，０）（ｘ（１））＝（ｘ２，ｘ２）Ｔをプロセス０へ送信し、次式により、ｙ（１）を計算する。
【0043】
ｙ（１）＝ＳｐＭＶ（Ａ（１），ｘ（１））
＝（０，ｆｘ３）Ｔ （７）
＝（０，ｆｘ３）Ｔ （７）
【0044】
次に、プロセス１は、ｚ＝ｐ（０，１）（ｘ（０））＝（ｘ０）Ｔをプロセス０から受信し、次式により、ｙ（１）を更新する。
【0045】
ｙ（１）＝ｙ（１）＋（Ａ（１，０））Ｔｚ
＝（ｅｘ０，ｆｘ３）Ｔ （８）
＝（ｅｘ０，ｆｘ３）Ｔ （８）
【0046】
この場合、Ａｘは、式（６）のｙ（０）及び式（８）のｙ（１）を用いて、次式により表される。
【0047】
【数3】
【0048】
式（９）のＡｘは、Ａ及びｘを分割することなく、Ａとｘの行列ベクトル積を計算した結果と一致している。
【0049】
図４のｉｆ文では、ｆｏｒループの終了条件としてｒｉ
（ｐ）＝０が用いられている。しかし、数値計算ではｒｉ
（ｐ）＝０となることがないため、例えば、以下のような終了条件が設定される。
【0050】
（ａ）解更新回数に基づく終了条件
ｉ＝Ｉｍａｘ
ｉ＝Ｉｍａｘ
【0051】
（ｂ）残差ノルムに基づく終了条件
ｎ（ｒｉ）＜ε
ｎ（ｒｉ）＜ε
【0052】
（ｃ）相対残差に基づく終了条件
ｎ（ｒｉ）／ｎ（ｒ０）＜τ
ｎ（ｒｉ）／ｎ（ｒ０）＜τ
【0053】
ｎ（ｒｉ）は、ｒｉのＬ１ノルムを表し、εは、ｎ（ｒｉ）に対する閾値を表す。ｎ（ｒｉ）／ｎ（ｒ０）は、ｎ（ｒ０）に対するｎ（ｒｉ）の比率を表し、τは、ｎ（ｒｉ）／ｎ（ｒ０）に対する閾値を表す。ｎ（ｒｉ）は、ｒｉ
（ｐ）を用いて、次式により表される。
【0054】
ｎ（ｒｉ）＝ｇＳｕｍＭａｇ（ｒｉ
（ｐ）） （１０）
【0055】
図７は、ｇＳｕｍＭａｇの疑似コードの例を示している。図７の疑似コードは、ｎ（ｘ（ｐ））の総和をｎ（ｘ）として出力する。ｘ（ｐ）をｒｉ
（ｐ）に置き換えた場合、ｎ（ｒｉ
（ｐ））の総和がｎ（ｒｉ）として出力される。この場合、ａｌｌｒｅｄｕｃｅ（ａ）によってプロセス間で送受信されるａは、解の一部に対する残差に関する情報の一例である。
【0056】
ｇＳｐＭＶ、ｇＳｕｍＰｒｏｄ、及びｇＳｕｍＭａｇは、ＯｐｅｎＦＯＡＭ（Open source Field Operation And Manipulation）において使用されている名称である。
【0057】
図８は、並列化されたＣＧ法の第２の疑似コードの例を示している。図８の疑似コードは、図４のｉｆ文において残差ノルムに基づく終了条件を用いた疑似コードに対応する。式８０２のδｉは、式（１０）のｎ（ｒｉ）に対応し、ｉｆ文８０３のδｉ＜εは、残差ノルムに基づく終了条件に対応する。式８０２及びｉｆ文８０３は、ｘｉ
（ｐ）の更新の終了条件を判定する処理を表す。
【0058】
式８０２のδｉが計算される時点では、ｘ１
（ｐ）～ｘｉ
（ｐ）が計算されているため、解がｉ回更新されている。ｉｆ文８０３の終了条件が満たされて、解の更新が打ち切られた場合、その時点におけるイテレータｉは、ＣＧ法の反復計算における解更新回数を表す。ｉｆ文８０３の判定は解の更新よりも先に行われるため、解更新回数が０回になることもあり得る。
【0059】
式８０１及び式８０５のｇＳｐＭＶを実行する際、プロセス間通信にかかる時間は、ＳｐＭＶの計算時間によって隠蔽することができる。一方、式８０２のｇＳｕｍＭａｇと式８０４及び式８０６のｇＳｕｍＰｒｏｄを実行する際、ａｌｌｒｅｄｕｃｅ通信によって計算がブロックされる。このため、ａｌｌｒｅｄｕｃｅ通信にかかる時間は隠蔽されない。
【0060】
行列ベクトル積、ベクトル演算等の計算時間は、プロセス数Ｐに応じてスケールするのに対して、ａｌｌｒｅｄｕｃｅ通信にかかる時間は、ｌｏｇＰのオーダーである。したがって、多数のプロセスが使用される場合、ＣＧ法におけるａｌｌｒｅｄｕｃｅ通信はボトルネックとなる。
【0061】
例えば、流体解析シミュレーションのように、時間発展系を扱うアプリケーションでは、同じ係数行列に対してＣＧ法を用いたループ処理が繰り返し実行され、各ループ処理における解更新回数が毎回類似する傾向が見られる。また、残差ノルムは急峻には変化せず、徐々に小さくなることが多いため、各ループ処理において、残差ノルムを必ずしも毎回計算する必要はない。
【0062】
この場合、残差ノルムが閾値よりも小さくなって、終了条件が満たされると予測されるループ処理の付近のみで残差ノルムを計算することで、残差ノルムの計算回数を削減することができる。さらに、そのループ処理の付近における残差ノルムの計算頻度自体を少なくすることで、残差ノルムの計算回数がさらに削減される。残差ノルムの計算回数を削減することで、式８０２におけるａｌｌｒｅｄｕｃｅ通信の回数が減少するため、各ループ処理の計算が高速化される。
【0063】
【0064】
図１０は、図９の情報処理装置９０１が行う第１の計算処理の例を示すフローチャートである。まず、演算処理部９１１は、１つ又は複数のループ処理各々において、並列に動作する複数の処理単位を用いて解の更新を反復する反復法の計算を実行する（ステップ１００１）。次に、演算処理部９１１は、１つ又は複数のループ処理よりも後の所定のループ処理において、複数の処理単位を用いて反復法の計算を実行する（ステップ１００２）。
【0065】
次に、演算処理部９１１は、１つ又は複数のループ処理各々の反復法の計算において解が更新された回数に基づき、所定のループ処理の反復法の計算において更新終了を判定する判定処理のタイミングを決定する（ステップ１００３）。判定処理は、複数の処理単位の間における通信の結果に基づいて、更新終了を判定する処理を含む。
【0066】
図９の情報処理装置９０１によれば、並列に動作する複数の処理単位を用いた反復法の計算時間を短縮することができる。
【0067】
図１１は、図９の情報処理装置９０１の具体例を示している。図１１の情報処理装置１１０１は、ＣＰＵ（Central Processing Unit）１１１１、記憶部１１１２、及び出力部１１１３を含む。ＣＰＵ１１１１は、図９の演算処理部９１１に対応する。
【0068】
情報処理装置１１０１は、様々なアプリケーションの数値計算において、式（１）の連立一次方程式の解を反復法により求めることができる。アプリケーションは、流体解析シミュレーション、気候シミュレーション、材料科学、生化学等の分野における分子動力学シミュレーション等である。反復法としては、例えば、ＣＧ法、ＰＣＧ法、ＢｉＣＧ法、又はＰＢｉＣＧ法が用いられる。
【0069】
例えば、流体解析シミュレーションは、ＰＩＳＯ（Pressure-Implicit with Splitting of Operators）法を用いて流体の定常解析を行うシミュレーションであってもよい。この場合、式（１）のベクトルｘは、物理量を表す。物理量は、圧力又は速度であってもよい。ｘが圧力を表す場合、ＣＧ法又はＰＣＧ法を用いてもよく、ｘが速度を表す場合、ＢｉＣＧ法又はＰＢｉＣＧ法を用いてもよい。
【0070】
ＣＰＵ１１１１は、複数のコアを含み、コア数以下の個数のプロセスを並列に動作させることができる。各コアは、演算回路であり、プロセスは、処理単位の一例である。ＣＰＵ１１１１は、所定時間間隔で、１つ又は複数のプロセスを用いてアプリケーションのループ処理を繰り返す。ＣＰＵ１１１１は、各時間ステップのループ処理において、反復法の計算を行うことで解ｘを求める。出力部１１１３は、各ループ処理において求められた解ｘと、そのループ処理における解更新回数とを出力する。
【0071】
記憶部１１１２は、ＣＰＵ１１１１が用いるデータとして、プロセス毎にアプリケーションオブジェクト１１２１及び反復法オブジェクト１１２２を記憶する。
【0072】
アプリケーションオブジェクト１１２１は、ｘ（ｐ）、更新履歴、及びアプリケーションデータを含む。ｘ（ｐ）は、ベクトルｘを分割することでプロセスｐに割り当てられたベクトルを表す。単一のプロセス０のみが動作する場合は、ベクトルｘは分割されず、ｘ（ｐ）＝ｘ（０）＝ｘとなる。更新履歴は、実行済みの各ループ処理における解更新回数を表す。解更新回数は、解ｘｉ
（ｐ）が更新された回数を表す。更新履歴は、プロセス間で一致している。アプリケーションデータは、アプリケーション特有のパラメータ等を表す。
【0073】
反復法オブジェクト１１２２は、ｘｉ
（ｐ）、ｒｉ
（ｐ）、ε、ＮＩ、ＮＨ、η、及びその他の反復法データを含む。イテレータｉは、反復法の計算におけるｆｏｒループのループ変数を表す。ｘｉ
（ｐ）は、ｉ番目のｆｏｒループにおける更新対象の解を表し、ｒｉ
（ｐ）は、ｉ番目のｆｏｒループにおける更新対象の残差を表す。εは、残差ノルムに対する閾値を表す。
【0074】
ＮＩは、反復法の計算において更新終了を判定する判定処理の周期を表す。判定処理では、ｘｉ
（ｐ）の更新を終了するか否かが判定される。ＮＨは、更新履歴に含まれる解更新回数の個数を表す。ηは、判定処理を開始する開始タイミングを決定するための係数を表す。ＮＩ及びＮＨは、１以上の整数であり、ηは、０以上１未満の実数である。ＮＩ、ＮＨ、及びηは、例えば、ユーザにより指定される。その他の反復法データは、係数行列Ａ、定数項を表すベクトルｂ、イテレータｉ等を含む。
【0075】
反復法がＢｉＣＧ法又はＰＢｉＣＧ法である場合、反復法オブジェクトは、ｉ番目のｆｏｒループにおける更新対象である転置版の残差をさらに含む。
【0076】
【0077】
ｉｆ文１２０１のｃｈｅｃｋＲｅｓｉｄｕａｌ（ｉ，｛ｓｊ｝）は、ｉ番目のｆｏｒループにおいて判定処理を行うか否かを決定する関数を表す。ｓｊは、ｊ回前（ｊ＝１，２，．．．，ＮＨ）のループ処理における解更新回数を表し、｛ｓｊ｝は、ｓ１～ｓＮＨを含む更新履歴を表す。ｃｈｅｃｋＲｅｓｉｄｕａｌ（ｉ，｛ｓｊ｝）は、次式により記述される。
【0078】
ｃｈｅｃｋＲｅｓｉｄｕａｌ（ｉ，｛ｓｊ｝）＝Ｔ
（ｉ≧ηｍｉｎ（ｓｊ）かつｉ％ＮＩ＝０） （１１）
ｃｈｅｃｋＲｅｓｉｄｕａｌ（ｉ，｛ｓｊ｝）＝Ｆ （otherwise） （１２）
（ｉ≧ηｍｉｎ（ｓｊ）かつｉ％ＮＩ＝０） （１１）
ｃｈｅｃｋＲｅｓｉｄｕａｌ（ｉ，｛ｓｊ｝）＝Ｆ （otherwise） （１２）
【0079】
ｍｉｎ（ｓｊ）は、｛ｓｊ｝に含まれるｓ１～ｓＮＨの最小値を表す。最小値は、統計値の一例である。１回前までに終了しているループ処理の実行回数がｊ回未満の場合、そのｓｊはｍｉｎ（ｓｊ）の計算対象から除外される。最初のループ処理（ｉ＝０）では、ｍｉｎ（ｓｊ）として０が用いられる。ｉ％ＮＩは、ｉをＮＩで除算したときの剰余を表す。
【0080】
Ｔは、真理値“真”を表し、Ｆは、真理値“偽”を表す。ｃｈｅｃｋＲｅｓｉｄｕａｌ（ｉ，｛ｓｊ｝）＝Ｔの場合、ｉ番目のｆｏｒループにおいて判定処理が行われ、ｃｈｅｃｋＲｅｓｉｄｕａｌ（ｉ，｛ｓｊ｝）＝Ｆの場合、ｉ番目のｆｏｒループにおいて判定処理がスキップされる。
【0081】
η＝０の場合、ηｍｉｎ（ｓｊ）＝０となる。したがって、ｆｏｒループがＮＩ回実行される間に１回の頻度で、ｃｈｅｃｋＲｅｓｉｄｕａｌ（ｉ，｛ｓｊ｝）＝Ｔとなる。ηが１に近いほど、反復法の開始直後の多数のｆｏｒループにおいて、判定処理を省略することができる。しかし、予測よりも速くδｉがε未満になった場合、解の収束の検出が遅れて、不要なｆｏｒループの反復が行われる可能性がある。η＝０かつＮＩ＝１の場合、ｆｏｒループにおいて毎回判定処理が行われる。
【0082】
図１３は、ＮＨ＝３の場合におけるｍｉｎ（ｓｊ）の変化の例を示している。ｓｊ（ｊ＝１～３）が未定義である場合、そのループ処理においてｓｊは存在しない。ｓは、各ループ処理における解更新回数を表し、次のループ処理においてｓ１として用いられる。ｋ回目のループ処理において、ｓ＝ξｋ－１である。各ループ処理は、所定のループ処理に対応する。
【0083】
１回目のループ処理では、ｓ１～ｓ３が未定義であるため、ｍｉｎ（ｓｊ）＝０、ｓ＝ξ０となる。２回目のループ処理では、ｓ２及びｓ３が未定義であるため、ｍｉｎ（ｓｊ）＝ｍｉｎ｛ξ０｝＝ξ０、ｓ＝ξ１となる。３回目のループ処理では、ｓ３が未定義であるため、ｍｉｎ（ｓｊ）＝ｍｉｎ｛ξ０，ξ１｝、ｓ＝ξ２となる。
【0084】
４回目のループ処理では、ｍｉｎ（ｓｊ）＝ｍｉｎ｛ξ０，ξ１，ξ２｝、ｓ＝ξ３となる。５回目のループ処理では、ｍｉｎ（ｓｊ）＝ｍｉｎ｛ξ１，ξ２，ξ３｝、ｓ＝ξ４となる。
【0085】
複数のプロセスが動作する場合、基本的に各プロセスは独立して動作するが、ｇＳｐＭＶ、ｇＳｕｍＭａｇ、ｇＳｕｍＰｒｏｄ等の実行時にはプロセス間で通信が行われ、データが交換される。
【0086】
ＣＰＵ１１１１は、ｉｆ文１２０１を実行することで、更新履歴に基づいて、各ループ処理において判定処理を開始する開始タイミングを決定する。開始タイミングは、各ループ処理において、最初にｉ≧ηｍｉｎ（ｓｊ）が満たされるｆｏｒループに対応する。更新履歴に基づいて開始タイミングを決定することで、反復法の開始から判定処理の開始タイミングまでの間のｆｏｒループにおいて、判定処理を省略することができる。
【0087】
ＣＰＵ１１１１は、判定処理の開始タイミング以降において、ｆｏｒループがＮＩ回実行される間に判定処理が１回行われるように、判定処理のタイミングを決定する。これにより、判定処理の開始タイミング以降のｆｏｒループにおける判定処理の回数を削減することができる。
【0088】
判定処理において、ＣＰＵ１１１１は、ａｌｌｒｅｄｕｃｅ通信を含むｇＳｕｍＭａｇ（ｒｉ
（ｐ））を実行することでδｉを計算し、計算されたδｉをεと比較することで、ｘｉ
（ｐ）の更新を終了するか否かを決定する。
【0089】
図１４は、図１２のループ処理における第１の動作の例を示している。図１４のループ処理では、η＝０かつＮＩ＝１である。横軸は、イテレータｉを表し、縦軸は、δｉを表す。曲線は、δｉが毎回計算されたと仮定した場合のδｉの変化を表す。縦の直線は、δｉが実際に計算されたタイミングを表し、横軸のｍは、０番目～ｉ番目のｆｏｒループにおいて、δｉが実際に計算された回数を表す。
【0090】
この場合、０番目～５番目のｆｏｒループにおいてδｉが毎回計算され、５番目のｆｏｒループにおいてδｉ＜εとなるため、解更新回数は５回となる。更新終了時のδｉの計算回数は６回である。
【0091】
【0092】
この場合、０番目、２番目、４番目、及び６番目のｆｏｒループにおいてδｉが計算され、６番目のｆｏｒループにおいてδｉ＜εとなるため、解更新回数は６回となる。更新終了時のδｉの計算回数は４回である。したがって、δｉの計算回数は、図１４のループ処理よりも２回減少し、解更新回数は、図１４のループ処理よりも１回増加する。
【0093】
【0094】
この場合、ηｍｉｎ（ｓｊ）＝１．５であり、２番目、４番目、及び６番目のｆｏｒループにおいてδｉが計算される。そして、６番目のｆｏｒループにおいてδｉ＜εとなるため、解更新回数は６回となる。更新終了時のδｉの計算回数は３回である。したがって、δｉの計算回数は、図１４のループ処理よりも３回減少し、解更新回数は、図１４のループ処理よりも１回増加する。
【0095】
このように、δｉの計算回数の削減と解更新回数の増加は、トレードオフの関係にある。η及びＮＩを適切な値に調整することで、解の不要な更新を抑えつつ、δｉの計算回数を効果的に削減することができる。
【0096】
図１７は、複数の時間ステップのループ処理における残差ノルムの計算回数の例を示している。図１７（ａ）は、１回目のループ処理における動作の例を示している。この場合、０番目～５番目のｆｏｒループにおいてδｉが毎回計算され、５番目のｆｏｒループにおいてδｉ＜εとなるため、解更新回数は５回となる。更新終了時のδｉの計算回数は６回である。
【0097】
図１７（ｂ）は、２回目のループ処理における動作の例を示している。この場合、０番目～４番目のｆｏｒループにおいてδｉが毎回計算され、４番目のｆｏｒループにおいてδｉ＜εとなるため、解更新回数は４回となる。更新終了時のδｉの計算回数は５回である。
【0098】
図１７（ｃ）は、３回目のループ処理における動作の例を示している。この場合、３番目及び５番目のｆｏｒループにおいてδｉが計算され、５番目のｆｏｒループにおいてδｉ＜εとなるため、解更新回数は５回となる。更新終了時のδｉの計算回数は２回である。したがって、δｉの計算回数は、解更新回数よりも３回少ない。
【0099】
図１２の疑似コードによれば、判定処理のみにおいて残差ノルムの計算が行われるため、残差ノルムの計算は、解及び残差を更新する計算には影響しない。したがって、図１５及び図１６に示したように、判定処理をスキップすることで、解更新回数が増加することはあるが、解の収束性が悪化することはない。判定処理をスキップすることで増加した解の更新時間よりも、削減されたプロセス間通信の通信時間の方が長い場合、連立一次方程式の解を求める処理全体の計算時間が短縮される。
【0100】
図１８は、図１１の情報処理装置１１０１が行う第２の計算処理の例を示すフローチャートである。ＣＰＵ１１１１は、アプリケーションオブジェクト１１２１及び反復法オブジェクト１１２２を用いて、アプリケーションプログラムを実行することで、図１８の計算処理を行う。計算処理は、流体解析シミュレーション、気候シミュレーション、分子動力学シミュレーション等に対応する。
【0101】
ＣＰＵ１１１１は、ステップ１８０１～ステップ１８０３の時間ステップループのループ処理を繰り返す。時間ステップループのループ処理において、ＣＰＵ１１１１は、アプリケーションオブジェクト１１２１を用いて、アプリケーションの計算を行う（ステップ１８０１）。アプリケーションの計算は、例えば、物理量を求める計算のうち、ＣＧ法を使用しない計算に対応する。
【0102】
次に、ＣＰＵ１１１１は、ステップ１８０２のＣＧ法ループの処理を繰り返す。ＣＧ法ループの処理は、図１２の疑似コードにおけるｆｏｒループの処理に対応し、ＣＰＵ１１１１は、各プロセスｐにｆｏｒループの処理を実行させる。ＣＧ法ループの処理において、ＣＰＵ１１１１は、反復法オブジェクト１１２２を用いて、プロセスｐ毎にｘ（ｐ）を更新することで、解ｘを更新する（ステップ１８０２）。ＣＧ法ループの処理は、残差ノルムが終了条件を満たすまで繰り返される。
【0103】
ＣＧ法ループの処理の繰り返しが終了した後、ＣＰＵ１１１１は、アプリケーションオブジェクト１１２１に含まれる更新履歴を更新する（ステップ１８０３）。そして、出力部１１１３は、求められた解と解更新回数とを出力する。時間ステップループの処理は、計算処理における最後の時間ステップに到達するまで繰り返される。
【0104】
【0105】
まず、ＣＰＵ１１１１は、ｃｈｅｃｋＲｅｓｉｄｕａｌ（ｉ，｛ｓｊ｝）がＴであるか否かをチェックする（ステップ１９０１）。ｃｈｅｃｋＲｅｓｉｄｕａｌ（ｉ，｛ｓｊ｝）がＦである場合（ステップ１９０１，ＮＯ）、ＣＰＵ１１１１は、ｘｉ
（ｐ）及びｒｉ
（ｐ）を更新して、ｘｉ＋１
（ｐ）及びｒｉ＋１
（ｐ）を計算する（ステップ１９０５）。そして、ＣＰＵ１１１１は、ｉを１だけインクリメントして、次のｆｏｒループの処理を行う。
【0106】
一方、ｃｈｅｃｋＲｅｓｉｄｕａｌ（ｉ，｛ｓｊ｝）がＴである場合（ステップ１９０１，ＹＥＳ）、ＣＰＵ１１１１は、δｉを計算し（ステップ１９０２）、計算されたδｉをεと比較する（ステップ１９０３）。
【0107】
δｉがε以上である場合（ステップ１９０３，ＮＯ）、ＣＰＵ１１１１は、ステップ１９０５の処理を行う。そして、ＣＰＵ１１１１は、ｉを１だけインクリメントして、次のｆｏｒループの処理を行う。一方、δｉがε未満である場合（ステップ１９０３，ＹＥＳ）、ＣＰＵ１１１１は、その時点におけるイテレータｉを解更新回数ｓに設定して（ステップ１９０４）、ｆｏｒループの処理を終了する。
【0108】
図２０は、図１８のステップ１８０３における履歴更新処理の例を示すフローチャートである。まず、ＣＰＵ１１１１は、ｊ＝ＮＩ～２について、ｊの降順に、ステップ２００１のｓｊ更新ループの処理を繰り返す。ｓｊ更新ループの処理において、ＣＰＵ１１１１は、更新履歴｛ｓｊ｝に含まれるｓｊ－１の値をｓｊに上書きする（ステップ２００１）。
【0109】
ｓｊ更新ループの処理の繰り返しが終了した後、ＣＰＵ１１１１は、ステップ１９０４で設定されたｓの値をｓ１に上書きする（ステップ２００２）。
【0110】
一例として、Ａが次のような４行４列の疎行列である場合に、ある時間ステップのループ処理で行われる解更新処理の例について説明する。
【0111】
【数4】
【0112】
ｂ＝（１，２，３，４）Ｔ、ｘ０＝（１，３，５，－３）Ｔである場合、解析解はｘ＝（０．８，２，３，３．７）Ｔである。プロセス０及びプロセス１の２個のプロセスを用いて連立方程式の解ｘを計算する場合、ｘｉ
（０）及びｘｉ
（１）として２次元ベクトルが用いられる。
【0113】
ε＝１０-５、ＮＩ＝２、ＮＨ＝２、η＝０．５、ｓ１＝３、ｓ２＝４である場合、ηｍｉｎ（ｓｊ）は、次式のように計算される。
【0114】
ηｍｉｎ（ｓｊ）＝０．５ｍｉｎ｛３，４｝＝１．５ （１４）
【0115】
図２１は、この場合の解更新処理の例を示している。ｃｈｅｃｋＲｅｓｉｄｕａｌの全体は、ｃｈｅｃｋＲｅｓｉｄｕａｌ（ｉ，｛ｓｊ｝）の値を表し、δｉ＜εは、終了条件の判定結果を表す。δｉ及びδｉ＜εの「計算しない」は、判定処理がスキップされることを表す。
【0116】
δｉは、倍精度浮動小数点数型で計算され、指数表示により小数点以下第３位まで表示されている。ｘｉ＋１
（０）及びｘｉ＋１
（１）は、倍精度浮動小数点数型で計算され、小数点以下第５位まで表示されている。
【0117】
ｉ＝０のｆｏｒループでは、式（１１）のｉ≧ηｍｉｎ（ｓｊ）の真理値がＦとなり、ｉ％ＮＩ＝０の真理値がＴとなる。したがって、全体の真理値はＦとなるため、判定処理がスキップされる。そして、プロセス０によりｘ１
（０）が計算され、プロセス１によりｘ１
（１）が計算される。
【0118】
ｉ＝１のｆｏｒループでは、ｉ≧ηｍｉｎ（ｓｊ）及びｉ％ＮＩ＝０の真理値がＦとなる。したがって、全体の真理値はＦとなるため、判定処理がスキップされる。そして、プロセス０によりｘ２
（０）が計算され、プロセス１によりｘ２
（１）が計算される。
【0119】
ｉ＝２のｆｏｒループでは、ｉ≧ηｍｉｎ（ｓｊ）及びｉ％ＮＩ＝０の真理値がＴとなる。したがって、全体の真理値はＴとなるため、判定処理が行われる。この場合、δｉ＜εの判定結果はＦであるため、解更新処理が継続される。そして、プロセス０によりｘ３
（０）が計算され、プロセス１によりｘ３
（１）が計算される。
【0120】
ｉ＝３のｆｏｒループでは、ｉ≧ηｍｉｎ（ｓｊ）の真理値がＴとなり、ｉ％ＮＩ＝０の真理値がＦとなる。したがって、全体の真理値はＦとなるため、判定処理がスキップされる。そして、プロセス０によりｘ４
（０）が計算され、プロセス１によりｘ４
（１）が計算される。
【0121】
ｉ＝４のｆｏｒループでは、ｉ≧ηｍｉｎ（ｓｊ）及びｉ％ＮＩ＝０の真理値がＴとなる。したがって、全体の真理値はＴとなるため、判定処理が行われる。この場合、δｉ＜εの判定結果はＦであるため、解更新処理が継続される。そして、プロセス０によりｘ５
（０）が計算され、プロセス１によりｘ５
（１）が計算される。
【0122】
ｉ＝５のｆｏｒループでは、ｉ≧ηｍｉｎ（ｓｊ）の真理値がＴとなり、ｉ％ＮＩ＝０の真理値がＦとなる。したがって、全体の真理値はＦとなるため、判定処理がスキップされる。そして、プロセス０によりｘ６
（０）が計算され、プロセス１によりｘ６
（１）が計算される。
【0123】
ｉ＝６のｆｏｒループでは、ｉ≧ηｍｉｎ（ｓｊ）及びｉ％ＮＩ＝０の真理値がＴとなる。したがって、全体の真理値はＴとなるため、判定処理が行われる。この場合、δｉ＜εの判定結果はＴであるため、解更新処理が打ち切られる。したがって、ｘ７
（０）及びｘ７
（１）は計算されない。そして、出力部１１１３は、ｘ６
（０）及びｘ６
（１）を解として出力し、ｓ＝６を解更新回数として出力する。
【0124】
図１２の疑似コードでは、残差ノルムに基づく終了条件が用いられているが、終了条件は、解更新回数又は相対残差に基づく終了条件であってもよく、複数の終了条件の組み合わせであってもよい。複数の終了条件の組み合わせは、複数の終了条件の論理和であってもよい。
【0125】
例えば、相対残差に基づく終了条件を用いる場合、ＣＰＵ１１１１は、ＣＧ法ループの開始時にｎ（ｒ０）＝ｇＳｕｍＭａｇ（ｒ０
（ｐ））を計算しておき、判定処理において、δｉ／ｎ（ｒ０）＜τが満たされるか否かをチェックする。
【0126】
ＣＧ法の代わりにＰＣＧ法、ＢｉＣＧ法、又はＰＢｉＣＧ法を用いた場合も、同様にして判定処理の回数を削減することで、連立一次方程式の解を求める処理の計算時間を短縮することができる。
【0127】
図２２は、複数のノードを含む情報処理装置のハードウェア構成例を示している。図２２の情報処理装置は、ノード２２０１－１～ノード２２０１－Ｕ（Ｕは２以上の整数）を含む。各ノード２２０１－ｕ（ｕ＝１～Ｕ）は、ＣＰＵ２２１１及びメモリ２２１２を含む。ＣＰＵ２２１１及びメモリ２２１２は、ハードウェアである。ＣＰＵ２２１１は、図９の演算処理部９１１に対応する。
【0128】
【0129】
ノード２２０１－１～ノード２２０１－Ｕは、通信ネットワーク２２０２を介して互いに通信することができる。通信ネットワーク２２０２は、例えば、Ethernet（登録商標）、InfiniBand、Tofuインターコネクト等の規格によるインターコネクトである。通信ネットワーク２２０２の接続方式は、メッシュ、トーラス、又はfat treeであってもよい。２つのノード２２０１－ｕそれぞれにおいて動作しているプロセス同士の通信は、通信ネットワーク２２０２を介して行われる。
【0130】
【0131】
図２３は、１００×１００×１００構造格子の物理量を計算する流体解析シミュレーションにおける計算時間の例を示している。この例では、OpenFOAM HPC Benchmark suiteのcavityと呼ばれる問題が用いられ、３２個のノードそれぞれにおいて４０個のプロセスを動作させている。したがって、プロセスの総数は１２８０個である。εは１０-４である。
【0132】
Ｃ１は、η＝０かつＮＩ＝１の場合のシミュレーション結果を表す。Ｃ２は、η＝０かつＮＩ＝１０の場合のシミュレーション結果を表す。Ｃ３は、η＝０．５かつＮＩ＝１０かつＮＨ＝１０の場合のシミュレーション結果を表す。Ｃ２及びＣ３の計算時間は、Ｃ１の計算時間よりも短縮されており、Ｃ３の計算処理は、Ｃ１の計算処理と比較して約１．２４倍高速化されている。
【0133】
【0134】
図２４（ａ）は、各時間ステップのループ処理毎のδｉの計算回数の例を示している。Ｃ１～Ｃ３各々の箱ひげ図において、箱内の×印は、計算回数の平均値を示し、箱内の横線は、計算回数の中央値を示す。箱の下端は、計算回数の第一四分位数を示し、箱の上端は、計算回数の第三四分位数を示し、ひげの下端は、計算回数の最小値を示し、ひげの上端は、計算回数の最大値を示す。Ｃ２及びＣ３の計算回数は、Ｃ１の計算回数よりも大幅に減少している。
【0135】
図２４（ｂ）は、各時間ステップのループ処理毎の解更新回数の例を示している。Ｃ１～Ｃ３各々の箱ひげ図において、箱内の×印は、解更新回数の平均値を示し、箱内の横線は、解更新回数の中央値を示す。箱の下端は、解更新回数の第一四分位数を示し、箱の上端は、解更新回数の第三四分位数を示し、ひげの下端は、解更新回数の最小値を示し、ひげの上端は、解更新回数の最大値を示す。Ｃ２及びＣ３の解更新回数は、Ｃ１の解更新回数と同程度である。
【0136】
図９の情報処理装置９０１、図１１の情報処理装置１１０１、及び図２２の情報処理装置の構成は一例に過ぎず、情報処理装置の用途又は条件に応じて一部の構成要素を省略又は変更してもよい。例えば、ＣＰＵの代わりに、ＧＰＵ（Graphics Processing Unit）等のアクセラレータを用いて、並列処理を実行することも可能である。この場合、プロセスの代わりに、スレッド等の別の処理単位を用いてもよい。
【0137】
【0138】
図１～図８及び図１２に示した疑似コードは一例に過ぎず、反復法のアルゴリズムは、別の形式で記述することもできる。図１３に示したｍｉｎ（ｓｊ）は一例に過ぎず、ｍｉｎ（ｓｊ）は、ＮＨに応じて変化する。図１４～図１７に示したループ処理の動作は一例に過ぎず、ループ処理の動作は、η、ＮＨ、及びＮＩに応じて変化する。
【0139】
図２１に示した解更新処理は一例に過ぎず、解更新処理は、η、ＮＨ、及びＮＩに応じて変化する。図２４に示したシミュレーション結果は一例に過ぎず、シミュレーション結果は、シミュレーションで用いられる連立方程式に応じて変化する。
【0140】
式（１）～式（１４）は一例に過ぎず、情報処理装置は、別の計算式を用いて計算処理を行うこともできる。例えば、式（１１）のｍｉｎ（ｓｊ）の代わりに、ｓ１～ｓＮＨの平均値、中央値、最頻値、最大値等の別の統計値を用いてもよい。
【0141】
図２５は、図９の情報処理装置９０１及び図１１の情報処理装置１１０１のハードウェア構成例を示している。図２５の情報処理装置は、ＣＰＵ２５０１、メモリ２５０２、入力装置２５０３、出力装置２５０４、補助記憶装置２５０５、媒体駆動装置２５０６、及びネットワーク接続装置２５０７を含む。これらの構成要素はハードウェアであり、バス２５０８により互いに接続されている。
【0142】
メモリ２５０２は、例えば、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）等の半導体メモリであり、処理に用いられるプログラム及びデータを記憶する。メモリ２５０２は、図１１の記憶部１１１２として動作してもよい。
【0143】
ＣＰＵ２５０１は、例えば、メモリ２５０２を利用してプログラムを実行することにより、図９の演算処理部９１１として動作する。ＣＰＵ２５０１は、図１１のＣＰＵ１１１１としても動作する。ＣＰＵ２５０１は、プロセッサと呼ばれることもある。
【0144】
入力装置２５０３は、例えば、キーボード、ポインティングデバイス等であり、ユーザ又はオペレータからの指示又は情報の入力に用いられる。出力装置２５０４は、例えば、表示装置、プリンタ等であり、ユーザ又はオペレータへの問い合わせ又は指示、及び処理結果の出力に用いられる。処理結果は、時間ステップ毎の解を含む計算結果であってもよい。出力装置２５０４は、図１１の出力部１１１３として動作してもよい。
【0145】
補助記憶装置２５０５は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク装置、テープ装置等である。補助記憶装置２５０５は、ハードディスクドライブ又はＳＳＤ（Solid State Drive）であってもよい。情報処理装置は、補助記憶装置２５０５にプログラム及びデータを格納しておき、それらをメモリ２５０２にロードして使用することができる。補助記憶装置２５０５は、図１１の記憶部１１１２として動作してもよい。
【0146】
媒体駆動装置２５０６は、可搬型記録媒体２５０９を駆動し、その記録内容にアクセスする。可搬型記録媒体２５０９は、メモリデバイス、フレキシブルディスク、光ディスク、光磁気ディスク等である。可搬型記録媒体２５０９は、ＣＤ－ＲＯＭ（Compact Disk Read Only Memory）、ＤＶＤ（Digital Versatile Disk）、ＵＳＢ（Universal Serial Bus）メモリ等であってもよい。ユーザ又はオペレータは、可搬型記録媒体２５０９にプログラム及びデータを格納しておき、それらをメモリ２５０２にロードして使用することができる。
【0147】
このように、処理に用いられるプログラム及びデータを格納するコンピュータ読み取り可能な記録媒体は、メモリ２５０２、補助記憶装置２５０５、又は可搬型記録媒体２５０９のような、物理的な（非一時的な）記録媒体である。
【0148】
ネットワーク接続装置２５０７は、ＬＡＮ（Local Area Network）、ＷＡＮ（Wide Area Network）等の通信ネットワークに接続され、通信に伴うデータ変換を行う通信インタフェース回路である。情報処理装置は、プログラム及びデータを外部の装置からネットワーク接続装置２５０７を介して受信し、それらをメモリ２５０２にロードして使用することができる。ネットワーク接続装置２５０７は、図１１の出力部１１１３として動作してもよい。
【0149】
なお、情報処理装置が図２５のすべての構成要素を含む必要はなく、情報処理装置の用途又は条件に応じて一部の構成要素を省略することも可能である。例えば、ユーザ又はオペレータとのインタフェースが不要である場合は、入力装置２５０３及び出力装置２５０４を省略してもよい。可搬型記録媒体２５０９又は通信ネットワークを使用しない場合は、媒体駆動装置２５０６又はネットワーク接続装置２５０７を省略してもよい。
【0150】
開示の実施形態とその利点について詳しく説明したが、当業者は、特許請求の範囲に明確に記載した本発明の範囲から逸脱することなく、様々な変更、追加、省略をすることができるであろう。
【0151】
図１乃至図２５を参照しながら説明した実施形態に関し、さらに以下の付記を開示する。
（付記１）
コンピュータのための計算プログラムであって、
前記計算プログラムは、
１つ又は複数のループ処理各々において、並列に動作する複数の処理単位を用いて解の更新を反復する反復法の計算を実行し、
前記１つ又は複数のループ処理よりも後の所定のループ処理において、前記複数の処理単位を用いて前記反復法の計算を実行し、
前記１つ又は複数のループ処理各々の前記反復法の計算において前記解が更新された回数に基づき、前記所定のループ処理の前記反復法の計算において更新終了を判定する判定処理のタイミングを決定する、
処理を前記コンピュータに実行させ、
前記判定処理は、前記複数の処理単位の間における通信の結果に基づいて、前記更新終了を判定する処理を含むことを特徴とする計算プログラム。
（付記２）
前記判定処理のタイミングを決定する処理は、前記１つ又は複数のループ処理各々の前記反復法の計算において前記解が更新された回数の統計値に基づき、前記所定のループ処理の前記反復法の計算において前記判定処理を開始する開始タイミングを決定する処理を含むことを特徴とする付記１記載の計算プログラム。
（付記３）
前記判定処理のタイミングを決定する処理は、前記所定のループ処理の前記反復法の計算において前記解の更新が複数回行われる間に、前記判定処理が１回行われるように、前記判定処理のタイミングを決定する処理を含むことを特徴とする付記１又は２記載の計算プログラム。
（付記４）
前記反復法の計算は、連立一次方程式の解を求める計算であり、
前記複数の処理単位各々には、前記連立一次方程式の解の一部を計算する処理が割り当てられ、
前記通信の結果に基づいて前記更新終了を判定する処理は、
前記複数の処理単位の間で、前記連立一次方程式の解の一部に対する残差に関する情報を送受信する処理と、
前記連立一次方程式の解の一部に対する残差に関する情報を用いて、前記連立一次方程式の解に対する残差を求める処理と、
前記連立一次方程式の解に対する残差に基づいて、前記解の更新を終了するか否かを判定する処理と、
を含むことを特徴とする付記１又は２記載の計算プログラム。
（付記５）
１つ又は複数のループ処理各々において、並列に動作する複数の処理単位を用いて解の更新を反復する反復法の計算を実行し、
前記１つ又は複数のループ処理よりも後の所定のループ処理において、前記複数の処理単位を用いて前記反復法の計算を実行し、
前記１つ又は複数のループ処理各々の前記反復法の計算において前記解が更新された回数に基づき、前記所定のループ処理の前記反復法の計算において更新終了を判定する判定処理のタイミングを決定する、
処理をコンピュータが実行し、
前記判定処理は、前記複数の処理単位の間における通信の結果に基づいて、前記更新終了を判定する処理を含むことを特徴とする計算方法。
（付記６）
前記判定処理のタイミングを決定する処理は、前記１つ又は複数のループ処理各々の前記反復法の計算において前記解が更新された回数の統計値に基づき、前記所定のループ処理の前記反復法の計算において前記判定処理を開始する開始タイミングを決定する処理を含むことを特徴とする付記５記載の計算方法。
（付記７）
前記判定処理のタイミングを決定する処理は、前記所定のループ処理の前記反復法の計算において前記解の更新が複数回行われる間に、前記判定処理が１回行われるように、前記判定処理のタイミングを決定する処理を含むことを特徴とする付記５又は６記載の計算方法。
（付記８）
前記反復法の計算は、連立一次方程式の解を求める計算であり、
前記複数の処理単位各々には、前記連立一次方程式の解の一部を計算する処理が割り当てられ、
前記通信の結果に基づいて前記更新終了を判定する処理は、
前記複数の処理単位の間で、前記連立一次方程式の解の一部に対する残差に関する情報を送受信する処理と、
前記連立一次方程式の解の一部に対する残差に関する情報を用いて、前記連立一次方程式の解に対する残差を求める処理と、
前記連立一次方程式の解に対する残差に基づいて、前記解の更新を終了するか否かを判定する処理と、
を含むことを特徴とする付記５又は６記載の計算方法。
（付記１）
コンピュータのための計算プログラムであって、
前記計算プログラムは、
１つ又は複数のループ処理各々において、並列に動作する複数の処理単位を用いて解の更新を反復する反復法の計算を実行し、
前記１つ又は複数のループ処理よりも後の所定のループ処理において、前記複数の処理単位を用いて前記反復法の計算を実行し、
前記１つ又は複数のループ処理各々の前記反復法の計算において前記解が更新された回数に基づき、前記所定のループ処理の前記反復法の計算において更新終了を判定する判定処理のタイミングを決定する、
処理を前記コンピュータに実行させ、
前記判定処理は、前記複数の処理単位の間における通信の結果に基づいて、前記更新終了を判定する処理を含むことを特徴とする計算プログラム。
（付記２）
前記判定処理のタイミングを決定する処理は、前記１つ又は複数のループ処理各々の前記反復法の計算において前記解が更新された回数の統計値に基づき、前記所定のループ処理の前記反復法の計算において前記判定処理を開始する開始タイミングを決定する処理を含むことを特徴とする付記１記載の計算プログラム。
（付記３）
前記判定処理のタイミングを決定する処理は、前記所定のループ処理の前記反復法の計算において前記解の更新が複数回行われる間に、前記判定処理が１回行われるように、前記判定処理のタイミングを決定する処理を含むことを特徴とする付記１又は２記載の計算プログラム。
（付記４）
前記反復法の計算は、連立一次方程式の解を求める計算であり、
前記複数の処理単位各々には、前記連立一次方程式の解の一部を計算する処理が割り当てられ、
前記通信の結果に基づいて前記更新終了を判定する処理は、
前記複数の処理単位の間で、前記連立一次方程式の解の一部に対する残差に関する情報を送受信する処理と、
前記連立一次方程式の解の一部に対する残差に関する情報を用いて、前記連立一次方程式の解に対する残差を求める処理と、
前記連立一次方程式の解に対する残差に基づいて、前記解の更新を終了するか否かを判定する処理と、
を含むことを特徴とする付記１又は２記載の計算プログラム。
（付記５）
１つ又は複数のループ処理各々において、並列に動作する複数の処理単位を用いて解の更新を反復する反復法の計算を実行し、
前記１つ又は複数のループ処理よりも後の所定のループ処理において、前記複数の処理単位を用いて前記反復法の計算を実行し、
前記１つ又は複数のループ処理各々の前記反復法の計算において前記解が更新された回数に基づき、前記所定のループ処理の前記反復法の計算において更新終了を判定する判定処理のタイミングを決定する、
処理をコンピュータが実行し、
前記判定処理は、前記複数の処理単位の間における通信の結果に基づいて、前記更新終了を判定する処理を含むことを特徴とする計算方法。
（付記６）
前記判定処理のタイミングを決定する処理は、前記１つ又は複数のループ処理各々の前記反復法の計算において前記解が更新された回数の統計値に基づき、前記所定のループ処理の前記反復法の計算において前記判定処理を開始する開始タイミングを決定する処理を含むことを特徴とする付記５記載の計算方法。
（付記７）
前記判定処理のタイミングを決定する処理は、前記所定のループ処理の前記反復法の計算において前記解の更新が複数回行われる間に、前記判定処理が１回行われるように、前記判定処理のタイミングを決定する処理を含むことを特徴とする付記５又は６記載の計算方法。
（付記８）
前記反復法の計算は、連立一次方程式の解を求める計算であり、
前記複数の処理単位各々には、前記連立一次方程式の解の一部を計算する処理が割り当てられ、
前記通信の結果に基づいて前記更新終了を判定する処理は、
前記複数の処理単位の間で、前記連立一次方程式の解の一部に対する残差に関する情報を送受信する処理と、
前記連立一次方程式の解の一部に対する残差に関する情報を用いて、前記連立一次方程式の解に対する残差を求める処理と、
前記連立一次方程式の解に対する残差に基づいて、前記解の更新を終了するか否かを判定する処理と、
を含むことを特徴とする付記５又は６記載の計算方法。
【符号の説明】
【0152】
１０１、８０３、１２０１ ｉｆ文
２０１、３０１、３０２、４０１～４０４、５０１、６０１、８０１、８０２、８０４～８０６ 式
９０１、１１０１ 情報処理装置
９１１ 演算処理部
１１０１ 情報処理装置
１１１１、２２１１、２５０１ ＣＰＵ
１１１２ 記憶部
１１１３ 出力部
１１２１ アプリケーションオブジェクト
１１２２ 反復法オブジェクト
２２０１－１～２２０１－Ｕ ノード
２２０２ 通信ネットワーク
２２１２、２５０２ メモリ
２５０３ 入力装置
２５０４ 出力装置
２５０５ 補助記憶装置
２５０６ 媒体駆動装置
２５０７ ネットワーク接続装置
２５０８ バス
２５０９ 可搬型記録媒体
２０１、３０１、３０２、４０１～４０４、５０１、６０１、８０１、８０２、８０４～８０６ 式
９０１、１１０１ 情報処理装置
９１１ 演算処理部
１１０１ 情報処理装置
１１１１、２２１１、２５０１ ＣＰＵ
１１１２ 記憶部
１１１３ 出力部
１１２１ アプリケーションオブジェクト
１１２２ 反復法オブジェクト
２２０１－１～２２０１－Ｕ ノード
２２０２ 通信ネットワーク
２２１２、２５０２ メモリ
２５０３ 入力装置
２５０４ 出力装置
２５０５ 補助記憶装置
２５０６ 媒体駆動装置
２５０７ ネットワーク接続装置
２５０８ バス
２５０９ 可搬型記録媒体
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】
【図22】
【図23】
【図24】
【図25】