(57)【要約】      （修正有）【課題】数値計算プログラムの生産性を向上することができる計算コード生成装置、方法及びプログラムを提供する。【解決手段】計算コード生成装置１は、計算対象の数学的問題に対する数値計算が実行される計算コードＣｏを自動生成する。計算コード生成装置は、記憶部と、制御部とを備える。記憶部は、計算対象の数学的問題に対する数値計算法を数式で宣言的に記述する形式記述情報Ｄｘを格納する。制御部は、記憶部に格納された形式記述情報に基づく演算処理を行う。制御部は、形式記述情報に基づいて、形式記述情報における数式及び数式が含む変数の関係を解析する。制御部は、解析結果に基づいて、所定の手続き型プログラミング言語により計算対象の関数に対する数値計算を実行する計算コードを生成する。形式記述情報は、反復型数値計算法を漸化式で規定する情報を含む。【選択図】図１
2017111749.tif 000002
【特許請求の範囲】
【請求項１】
計算対象の数学的問題に対する数値計算が実行される計算コードを自動生成する計算コード生成装置であって、
数値計算法を数式で宣言的に記述する形式記述情報を格納する記憶部と、
前記記憶部に格納された形式記述情報に基づく演算処理を行う制御部とを備え、
前記制御部は、
前記形式記述情報に基づいて、前記形式記述情報における数式及び前記数式が含む変数の関係を解析して、
解析結果に基づいて、所定のプログラミング言語により前記数値計算を実行する計算コードを生成し、
前記形式記述情報は、反復型数値計算法を漸化式で規定する情報を含む
計算コード生成装置。
【請求項２】
前記制御部は、前記形式記述情報における漸化式を、前記漸化式に含まれる変数を介して関連する一連の数式に展開することにより、前記解析を行う
請求項１に記載の計算コード生成装置。
【請求項３】
前記制御部は、前記記憶部に格納された形式記述情報における数式毎に各数式に含まれる変数を分類して、前記漸化式で規定する情報を解釈する
請求項１又は２に記載の計算コード生成装置。
【請求項４】
前記制御部は、
前記計算コードが実行される計算環境を示す計算環境情報を取得し、
前記計算環境情報に応じて最適化して、前記計算コードを出力する
請求項１〜３のいずれか１項に記載の計算コード生成装置。
【請求項５】
前記形式記述情報における漸化式は、少なくとも３回にわたる反復計算回数にそれぞれ対応する少なくとも３つの変数を含む
請求項１〜４のいずれか１項に記載の計算コード生成装置。
【請求項６】
前記形式記述情報は、前記数値計算法をモジュール化して記述される
請求項１〜５のいずれか１項に記載の計算コード生成装置。
【請求項７】
前記生成された計算コードの外部から、前記数学的問題の具体形が指定される
請求項１〜６のいずれか１項に記載の計算コード生成装置。
【請求項８】
前記生成された計算コードの内部に、前記数学的問題の具体形が組み込まれている
請求項１〜６のいずれか１項に記載の計算コード生成装置。
【請求項９】
コンピュータが、計算対象の数学的問題に対する数値計算が実行される計算コードを自動生成する方法であって、
前記コンピュータの記憶部には、前記計算対象の数学的問題に対する数値計算法を数式で宣言的に記述する形式記述情報が格納されており、
前記コンピュータの制御部が、前記形式記述情報に基づいて、前記形式記述情報における数式及び前記数式が含む変数の関係を解析するステップと、
前記制御部が、解析結果に基づいて、所定のプログラミング言語により前記計算対象の関数に対する数値計算を実行する計算コードを生成するステップとを含み、
前記形式記述情報は、反復型数値計算法を漸化式で規定する情報を含む
方法。
【請求項１０】
請求項９に記載の方法をコンピュータに実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、数値計算法が実行されるプログラムコードを自動生成する計算コード生成装置、方法及びプログラムに関する。
【背景技術】
【０００２】
従来、数値計算アルゴリズムにおいては、計算性能が重要視されるため、特定の手続型言語による参照実装が一般的に公開されている。しかし、この手法には移植性や保守性の問題があることから、機械的に数値計算法のプログラムコードを生成するための技術が提案されている（例えば特許文献１，非特許文献１）。
【０００３】
特許文献１は、数値計算アルゴリズムに関する問題定義情報とプログラム様式情報の入力により、数値計算プログラムを生成する数値計算プログラム生成システムを開示している。特許文献１の数値計算プログラム生成システムによると、生成される数値計算プログラムの数値計算アルゴリズムを疑似コードの順序列で表現したものからプログラム生成を行う。これにより、既に定義されている疑似コードを利用して新たな数値計算アルゴリズムに対応する疑似コードの順序列を作り、数値計算プログラム生成システムが生成可能な数値計算アルゴリズムを拡張できる。
【０００４】
非特許文献１は、細胞生理学モデルを記述したファイルとＯＤＥ（常微分方程式）解法スキームを記述したファイルから時間発展計算を行うシミュレーションプログラムの自動生成手法を開示している。非特許文献１の自動生成手法によると、計算機が、上記の各ファイルの入力によって時間発展を計算する数式セットを自動生成し、生成された数式セットから実行コードを生成する。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２００４−８６７６０号公報
【非特許文献】
【０００６】
【非特許文献１】山下義陽、外７名、「形式的に記述されたODE解法スキームに基づくCellMLシミュレーションコード生成システム」、生体医工学、50(1)、p.68−77、2012年
【非特許文献２】Ｂｕｎｕｓ  Ｐ．，Ｆｒｉｔｚｓｏｎ  Ｐ．，“Ｍｅｔｈｏｄｓ  ｆｏｒ  Ｓｔｒｕｃｔｕｒａｌ  Ａｎａｌｙｓｉｓ  ａｎｄ  Ｄｅｂｕｇｇｉｎｇ  ｏｆ  Ｍｏｄｅｌｉｃａ  Ｍｏｄｅｌｓ”，２ｎｄ  Ｉｎｔｅｒｎａｔｉｏｎａｌ  Ｍｏｄｅｌｉｃａ  Ｃｏｎｆｅｒｅｎｃｅ，Ｐｒｏｃｅｅｄｉｎｇｓ，ｐｐ．１５７−１６５．
【非特許文献３】Ｊ．Ｕｎｇｅｒ，ｅｔ  ａｌ．，“ＳＴＲＵＣＴＵＲＡＬ  ＡＮＡＬＹＳＩＳ  ＯＦ  ＤＩＦＦＥＲＥＮＴＩＡＬ−ＡＬＧＥＢＲＡＩＣ  ＥＱＵＡＴＩＯＮ  ＳＹＳＴＥＭＳ−ＴＨＥＯＲＹ  ＡＮＤ  ＡＰＰＬＩＣＡＴＩＯＮＳ”，Ｃｏｍｐｕｔｅｒｓ  ｃｈｅｍ．，Ｅｎｇｎｇ  Ｖｏｌ．１９，Ｎｏ．８，ｐｐ．８６７−８８２，１９９５．
【発明の概要】
【発明が解決しようとする課題】
【０００７】
特許文献１の数値計算プログラム生成システムは、システム内で既に定義された疑似コードを利用して、入力される問題定義情報に応じた数値計算プログラムを生成している。このため、所望の数値計算法のアルゴリズム自体を入力として指定することで、数値計算プログラムを生成させることはできない。
【０００８】
非特許文献１の方法では、実行コードを自動生成する対象の数値計算法が、常微分方程式の初期値問題に対する修正Euler法や陽的Runge-Kutta法などの一段階法の一部に限定されている。非特許文献１の方法は、常微分方程式の初期値問題の解法以外の数値計算法や、常微分方程式の初期値問題の解法の中でも例えば陰的Runge-Kutta法や線形多段階法には適用できず、種々の数値計算法の実行コードを自動生成することができなかった。
【０００９】
本発明の目的は、数値計算プログラムの生産性を向上することができる計算コード生成装置、方法及びプログラムを提供することである。
【課題を解決するための手段】
【００１０】
本発明に係る計算コード生成装置は、計算対象とする数学的問題に対する数値計算が実行される計算コードを自動生成する計算コード生成装置である。計算コード生成装置は、記憶部と、制御部とを備える。記憶部は、計算対象の問題に対する数値計算法を数式で宣言的に記述する形式記述情報を格納する。制御部は、記憶部に格納された形式記述情報に基づく演算処理を行う。制御部は、形式記述情報に基づいて、形式記述情報における数式及び数式が含む変数の関係を解析する。制御部は、解析結果に基づいて、所定のプログラミング言語により計算対象とする問題に対する数値計算を実行する計算コードを生成する。形式記述情報は、反復型数値計算法を漸化式で規定する情報を含む。
【発明の効果】
【００１１】
本発明に係る計算コード生成装置によると、形式記述情報に記述された漸化式に基づき、種々の反復型数値計算法を実行する計算コードが生成可能であり、数値計算プログラムの生産性を向上することができる。
【図面の簡単な説明】
【００１２】
【図１】本発明の実施形態１に係る計算コード生成装置による計算コードの自動生成方法の概念を説明するための図
【図２】実施形態１に係る計算コード生成装置の構成を示すブロック図
【図３】実施形態１に係る計算コード生成処理を示すフローチャート
【図４】実施形態１における形式記述情報の一例を示す図
【図５】図４に続く形式記述情報の例を示す図
【図６】図５に続く形式記述情報の例を示す図
【図７】計算コード生成処理における形式記述解釈処理を示すフローチャート
【図８】形式記述解釈処理で生成される各種情報を説明するための図
【図９】計算コード生成処理における数式構造解析処理を示すフローチャート
【図１０】数式構造解析処理における数式の展開を説明するための図
【図１１】数式構造解析処理における二部グラフを説明するための図
【図１２】数式構造解析処理で解析される諸情報を説明するための図
【図１３】計算コード生成処理におけるコード生成処理を示すフローチャート
【図１４】実施形態１における形式記述情報の別例を示す図
【図１５】図１４に続く形式記述情報の例を示す図
【図１６】計算コード生成処理における最適化前後の計算コードを示す図
【図１７】形式記述情報の変形例を示す図
【図１８】変形例の形式記述情報における第１の参照情報を示す図
【図１９】変形例の形式記述情報における第２の参照情報を示す図
【発明を実施するための形態】
【００１３】
以下、添付の図面を参照して本発明に係る計算コード生成装置、方法及びプログラムの実施の形態を説明する。なお、以下の各実施形態において、同様の構成要素については同一の符号を付している。
【００１４】
（実施形態１）
１．構成
１−１．概要
本発明の実施形態１に係る計算コード生成装置による計算コードの自動生成方法の概念について、図１を参照して説明する。図１は、本実施形態に係る計算コード生成装置１による計算コードの自動生成方法の概念を説明するための図である。
【００１５】
本実施形態に係る計算コード生成装置１は、例えばユーザが指定する数値計算法を記述した形式記述情報Ｄｘを入力し、指定された数値計算法を実行する計算コードＣｏを自動生成する。
【００１６】
形式記述情報Ｄｘは、任意の数値計算法の完全な定義が、形式仕様記述において宣言的に記述された情報である（図４〜６参照）。形式記述情報Ｄｘにおいて、数値計算法のアルゴリズムは数式を用いて記述され、数式は例えばＭａｔｈＭＬ形式やＯｐｅｎＭａｔｈ形式などの形式記述書式により記述される。また、数値計算法が計算対象とする数学的問題の一般形を表す方程式（方程式系）における当該問題を定める関数（以下、「問題定義関数」という。）は、形式記述情報Ｄｘ中では変数として記述される。例えば、常微分方程式の初期値問題の一般形「ｄｙ／ｄｔ＝ｆ（ｙ，ｔ）」において、問題定義関数はｆである。
【００１７】
計算コードＣｏは、所定のプログラミング言語において、形式記述情報Ｄｘで指定された数値計算法による数値計算が行われるように記述されたプログラムコードである。所定のプログラミング言語には、例えばＣ言語、Ｃ＋＋、Ｊａｖａ（登録商標），Ｆｏｒｔｒａｎ，ＣＵＤＡ  Ｃ，Ｈａｓｋｅｌｌ，Ｓｃａｌａなどが含まれる。本実施形態では、問題定義関数（数学的問題）の具体形は、別途、計算コードＣｏに対する入力として指定される。
【００１８】
本実施形態に係る計算コード生成処理の対象は、常微分方程式の初期値問題の解法や境界値問題の解法、偏微分方程式の解法、（非線形）連立方程式の求根アルゴリズム、最適化問題の解探索アルゴリズム、数値連続法、数値微分法、数値積分法など、数値計算法全般にわたる。特に、数値計算法の内部で反復計算処理を行う反復型数値計算法は、形式記述情報Ｄｘにおいて、反復計算回数を添字集合とする列を有する漸化式として明示的に表現する。計算コード生成処理において漸化式の諸構造を解析することにより、例えば線形多段階法などの複雑な解法を用いる場合でも、計算コードＣｏを自動生成することができる。
【００１９】
また、本実施形態に係る計算コード生成装置１は、計算環境情報Ｄｅに応じて、計算コードＣｏを最適化して出力する。計算環境情報Ｄｅは、計算コードＣｏが実行される計算環境における各種スペック（例えば搭載メモリ量やＣＰＵ数）に関する情報である。また、計算コードＣｏを記述するプログラミング言語も、計算環境情報Ｄｅにおいて指定される。本実施形態では、形式記述情報Ｄｘに記述された数式を解析することで、機械的に判定可能な範囲で最適化された計算コードＣｏが得られる。
【００２０】
以下、本実施形態に係る計算コード生成装置１のハードウェア構成について説明する。
【００２１】
１−２．ハードウェア構成
図２は、計算コード生成装置１の構成を示すブロック図である。計算コード生成装置１は、例えばＰＣ（パーソナルコンピュータ）などの情報処理装置で構成される。計算コード生成装置１は、図２に示すように、制御部１０と、記憶部１１と、操作部１２と、表示部１３と、機器インタフェース１４と、ネットワークインタフェース１５とを備える。
【００２２】
制御部１０は、例えばソフトウェアと協働して所定の機能を実現するＣＰＵやＧＰＧＰＵ、ＭＰＵで構成され、計算コード生成装置１の全体動作を制御する。制御部１０は、記憶部１１に格納されたデータやプログラムを読み出して種々の演算処理を行い、各種の機能を実現する。例えば、御部１０は、本実施形態に係る計算コードの自動生成方法が実現される計算コード生成処理を実行する。計算コード生成処理を実行するためのプログラムは、パッケージソフトウェアであってもよい。また、制御部１０は、所定の機能を実現するように設計された専用の電子回路や再構成可能な電子回路などのハードウェア回路であってもよい。制御部１０は、ＣＰＵ，ＭＰＵ，マイコン、ＤＳＰ、ＦＰＧＡ、ＡＳＩＣ等の種々の半導体集積回路で構成されてもよい。
【００２３】
記憶部１１は、計算コード生成装置１の機能を実現するために必要なプログラム及びデータを記憶する記憶媒体であり、例えばハードディスク（ＨＤＤ）や半導体記憶装置（ＳＳＤ）を備える。また、記憶部１１は、例えば、ＤＲＡＭやＳＲＡＭ等の半導体デバイスを備えてもよく、データを一時的に記憶するとともに制御部１０の作業エリアとしても機能する。記憶部１１において、形式記述情報Ｄｘや計算環境情報Ｄｅ、計算コード入力リストＬｉｎ、数列リストＬ１，数式変数リストＬ２（後述）などが格納される。
【００２４】
操作部１２は、ユーザが操作を行うユーザインタフェースである。操作部１２は、例えば、キーボード、タッチパッド、タッチパネル、ボタン、スイッチ、及びこれらの組み合わせで構成される。操作部１２は、ユーザによって入力される諸情報を取得する取得部の一例である。
【００２５】
表示部１３は、例えば、液晶ディスプレイや有機ＥＬディスプレイで構成される。表示部１３は、例えば操作部１２から入力された情報など、種々の情報を表示する。
【００２６】
機器インタフェース１４は、計算コード生成装置１に他の機器を接続するための回路（モジュール）である。機器インタフェース１４は、所定の通信規格にしたがい通信を行う。所定の規格には、ＵＳＢ、ＨＤＭＩ（登録商標）、ＩＥＥＥ１３９５、ＷｉＦｉ、Ｂｌｕｅｔｏｏｔｈ（登録商標）等が含まれる。
【００２７】
ネットワークインタフェース１５は、無線または有線の通信回線を介して計算コード生成装置１をネットワークに接続するための回路（モジュール）である。ネットワークインタフェース１５は所定の通信規格に準拠した通信を行う。所定の通信規格には、ＩＥＥＥ８０２．３，ＩＥＥＥ８０２．１１ａ／１１ｂ／１１ｇ／１１ａｃ等の通信規格が含まれる。
【００２８】
２．動作
以下、本実施形態に係る計算コード生成装置１の動作について説明する。
【００２９】
２−１．計算コード生成処理
本実施形態に係る計算コード生成装置１によって実行される処理について、図３〜６を参照して説明する。図３は、本実施形態に係る計算コード生成処理を示すフローチャートである。図４，５，６は、形式記述情報Ｄｘの一例を示す図である。
【００３０】
計算コード生成処理は、所望の数値計算法が形式的に記述された形式記述情報Ｄｘに基づき、数値計算法が計算環境情報Ｄｅにより指定される計算環境下で実行される計算コードＣｏを生成する処理である。本処理は、計算コード生成装置１の制御部１０によって実行される。
【００３１】
まず、計算コード生成装置１の制御部１０は、形式記述情報Ｄｘに基づいて、形式記述解釈処理を行う（Ｓ１）。形式記述解釈処理では、計算コード生成装置１に入力される形式記述情報Ｄｘを取得し、形式記述情報Ｄｘにおける数値計算法の入力／出力の指定や反復条件、数値計算法のアルゴリズムを規定する数式の記述の解釈（構文解析）を行う。
【００３２】
ここで、形式記述情報Ｄｘについて説明する。形式記述情報Ｄｘには、所望の数値計算法に対する入力や出力の指定、及び数値計算法のアルゴリズムを規定する数式などの各種情報が、宣言的に記述されている。図４〜６に、形式記述情報Ｄｘの一例を示す。
【００３３】
形式記述情報Ｄｘは、例えばＸＭＬ形式で記述される。図４〜６に示すように、形式記述情報Ｄｘにおいて、記述ＸＡは計算コードＣｏに対する入力を指定する部分であり、記述ＸＢは計算コードＣｏによる出力を指定する部分であり、記述ＸＣは計算コードＣｏで実行する反復計算の反復条件を指定する部分である。また、記述ＸＤは、形式記述情報Ｄｘにおいて反復型数値計算法のアルゴリズムを規定する数式を記述したアルゴリズム部分である。記述ＸＥは、数値計算時に使用するパラメータを指定する部分である。記述ＸＡ〜ＸＣ，ＸＥの指定内容については、後述する。
【００３４】
図４〜６においては、反復型数値計算法が常微分方程式（ｄｙ／ｄｔ＝ｆ（ｙ，ｔ））の初期値問題に対する線形多段階法の一種である２次Adams-Bashforth法（以下、「２次ＡＢ法」と略記する。）の場合の一例を示している。本例では、２次ＡＢ法にEuler法を組み合わせたアルゴリズムを用いている。本例の反復型数値計算法のアルゴリズムを規定する漸化式は、以下のように表される。
ｙｉ＋１＝ｙｉ＋ｈ×（３ｇｉ−ｇｉ−１）／２      （１）
ｇｉ＝ｆ（ｙｉ，ｔｉ）                        （２）
ｔｉ＋１＝ｔｉ＋ｈ                              （３）
ｙ１＝ｙ０＋ｈ×ｇ０                            （４）
【００３５】
上式（１）〜（３）は、２次ＡＢ法のアルゴリズムを規定している。また、上式（４）において、ｙ１の計算にEuler法を用いている。図４〜６に示すように、形式記述情報Ｄｘのアルゴリズム部分ＸＤにおいて、記述ＸＤ１に上式（１）が記述され、記述ＸＤ２に上式（２）が記述され、記述ＸＤ３に上式（３）が記述され、記述ＸＤ４に上式（４）が記述されている。
【００３６】
上記の漸化式の組（１）〜（４）において、ｉは、反復計算における反復回数に対応する添字（インデックス）である。例えば、変数ｙｉは、関数ｆ（ｙ，ｔ）の引数変数ｙに対する反復計算ｉ回目の値に対応する変数であり、変数ｙｉの組｛ｙ０，ｙ１，…，ｙｉ，…｝はｙの反復計算に関する数列を構成する。各変数は、数値（スカラー）であってもよいし、ベクトル、行列、又は集合などであってもよい。また、これに応じて、数列は、各項が数値で構成される列だけでなく、ベクトル列など、あらゆる列を含むものとする。以下、変数ｙｉのように、添字が付された変数であって添字に基づく数列の一項を構成する変数を、「添字付き変数」という。
【００３７】
上式（１）中の添字付き変数ｙｉ＋１とｙｉのように、同じ数列に属する添字付き変数は、計算コードＣｏ上では別々の変数として扱われることとなるが、計算コードＣｏの生成時には、反復計算の繰り返し前後の計算ステップ等で関連付けをする必要もある。このため、本実施形態の形式記述解釈処理では、形式記述情報Ｄｘにおける漸化式に含まれる種々の変数や、各変数が属する数列に関する情報の分類を行う。形式記述解釈処理の詳細については後述する。
【００３８】
形式記述解釈処理（Ｓ１）の後に、制御部１０は、形式記述解釈処理（Ｓ１）の解析結果に基づいて、数式構造解析処理を行う（Ｓ２）。数式構造解析処理では、数値計算法のアルゴリズムを規定する数式の組における数式構造（変数を介した数式間の関係）を解析し、数値計算法を種々のプログラミング言語で記述可能にするための情報（中間データ）を生成する。
【００３９】
例えば、２次ＡＢ法の場合、上記の漸化式の組（１）〜（４）が数式構造解析処理の対象となる。漸化式は、上式（１）〜（４）に示すように、あらゆるｉ（ｉ＝１，２，…）について共通の式である。ここで、漸化式の組における数式間で、添字の相対的な関係を考慮する必要がある場合がある。例えば、上式（１），（２）において、式（１）の右辺のｇｉは式（２）により定義されることが明示的だが、式（１）の右辺のｇｉ−１も式（２）により定義されることは字句上では明示的でない。そこで、本実施形態の数式構造解析処理では、反復型数値計算法のアルゴリズムを規定する漸化式に対して、数値計算に使用される別々の添字付き変数が明示されるように、数式の展開を行う。数式構造解析処理の詳細については後述する。
【００４０】
次に、制御部１０は、数式構造解析処理（Ｓ２）の解析結果及び計算環境情報Ｄｅに基づいて、コード生成処理を行う（Ｓ３）。コード生成処理では、数式構造解析処理で得られた中間データを、計算環境情報Ｄｅで指定されるプログラミング言語で変換し、計算コードＣｏを生成する。本実施形態のコード生成処理では、計算環境情報Ｄｅが指定する計算環境に応じて自動化可能な範囲で最適化された計算コードが生成される。
【００４１】
本実施形態では、形式記述情報Ｄｘにおいて数値計算法のアルゴリズムが数式で指定されるため、数式中で並列計算可能な部分や計算結果が共通する部分を容易に判定可能である。特に、反復型数値計算法のアルゴリズムを指定する漸化式では、上記の漸化式（１）において反復計算（ｉ−１）回目で計算されるｇｉ−１のように、反復計算の過去の計算ステップの計算結果が再利用可能な部分を特定可能である。ここで、ｇｉ−１について計算結果を記憶することと再計算することのいずれが適切であるかは、計算環境に依存する。本実施形態のコード生成処理では、計算環境に応じて上記の部分等の最適化を行う。コード生成処理の詳細については後述する。
【００４２】
以上の各処理の詳細について、以下、説明する。
【００４３】
２−１−１．形式記述解釈処理
図３のフローチャートにおける形式記述解釈処理（ステップＳ１）について、図４〜８を参照して説明する。図７は、形式記述解釈処理を示すフローチャートである。図８は、形式記述解釈処理で生成される各種情報を示す図である。
【００４４】
以下では、形式記述解釈処理について、入力の数値計算法が上述の２次ＡＢ法の一例の場合について説明する。また、以下では、形式記述情報Ｄｘ（図４〜６参照）において、漸化式の組（１）〜（４）を示す宣言的記述（ＸＤ）とともに、以下の項目（Ｉ）〜（ＩＶ）が記述されている例で説明する（ＸＡ〜ＸＣ，ＸＥ）。
（Ｉ）計算コードＣｏの入力は、変数ｙ０，ｔ０，ｆである。
（ＩＩ）計算コードＣｏの出力は、全反復における変数の組（ｙｎ，ｔｎ）である。
（ＩＩＩ）反復計算の反復条件は、継続条件ｔ≦Ｔである。
（ＩＶ）パラメータ定数は、ｈ，Ｔである。
【００４５】
なお、反復条件は、反復型数値計算アルゴリズムによる反復計算を繰り返すか否かを判定するための条件である。反復条件は、反復計算の継続を判定する継続条件で指定されてもよいし（例えばｔ≦Ｔ）、反復計算の終了を判定する終了条件で指定されてもよい（例えばｔ＞Ｔ）。
【００４６】
図７のフローチャートにおいて、まず、制御部１０は、記憶部１１から形式記述情報Ｄｘを取得する（Ｓ１１）。形式記述情報Ｄｘは、例えば本処理の開始前にユーザによって操作部１２（図２）から入力されて、記憶部１１に格納されている。なお、形式記述情報Ｄｘにおいて他の情報を参照する記述がある場合、制御部１０は参照されている情報も読み出し、形式記述情報Ｄｘの記述に従ってこれらの情報を合成する（図１７〜１９参照）。
【００４７】
次に、制御部１０は、読み出した形式記述情報Ｄｘにおいて計算コードＣｏに対する入力を指定する記述（項目（Ｉ））を解釈し、計算コードＣｏに対して入力される情報を抽出する（Ｓ１２）。例えば、上述の２次ＡＢ法の一例の形式記述情報Ｄｘには、図４に示すように、計算コードＣｏの入力の記述ＸＡがある。制御部１０は、この記述ＸＡに基づいて、図８（ａ）に示すように、計算コード入力リストＬｉｎを生成する。計算コード入力リストＬｉｎは、形式記述情報Ｄｘに基づき生成される計算コードＣｏの使用時に入力されることとなる情報を管理するリストである。
【００４８】
図８（ａ）に示すように、計算コード入力リストＬｉｎは、「入力変数」と、「変数属性」とを関連付けている。入力変数は、計算コードの使用時に入力として指定される変数である。変数属性は、変数の型を示す属性情報である。
【００４９】
図８（ａ）に示すように、計算コード入力リストＬｉｎには、入力変数ｙ０，ｔ０が「定数添字付き変数」という型の変数であり、入力変数ｆが「関数」を表す変数であることが記録されている。以下、定数だけからなる添字を「定数添字」といい、添字が定数添字である添字付き変数を「定数添字付き変数」という。また、変数を含む添字を「変数添字」といい、添字が変数添字である添字付き変数を「変数添字付き変数」という。
【００５０】
次に、制御部１０は、形式記述情報Ｄｘにおいて計算コードＣｏの出力を指定する記述（項目（ＩＩ））を解釈し、計算コードＣｏによる数値計算結果として出力される情報を抽出する（Ｓ１３）。例えば、２次ＡＢ法の形式記述情報Ｄｘには、図４に示すように、計算コードＣｏの出力の記述ＸＢがある。この記述ＸＢから、制御部１０は、２次ＡＢ法の数値計算結果として、全反復における変数の組（ｙｎ，ｔｎ）の評価値が出力されることを判断する。
【００５１】
次に、制御部１０は、形式記述情報Ｄｘにおいて計算コードＣｏの数値計算アルゴリズムで用いられる反復条件やパラメータ定数を指定する記述（項目（ＩＩＩ），（ＩＶ））を解釈する（Ｓ１４）。制御部１０は、例えば計算コード入力リストＬｉｎと同様に、形式記述情報Ｄｘに記述されたパラメータ定数をリスト化する。
【００５２】
例えば、２次ＡＢ法の形式記述情報Ｄｘでは、図４に示すように、反復条件の記述ＸＣがある。この記述ＸＣから、制御部１０は、反復計算が、ｔ≦Ｔである限り継続されることを判断する。また、図６に示すように、計算コードＣｏのパラメータ定数ｈ，Ｔの記述ＸＥがある。制御部１０はこの記述ＸＥを抽出して、２次ＡＢ法の数値計算時にパラメータ定数ｈ，Ｔが用いられることを判断する。
【００５３】
次に、制御部１０は、ステップＳ１５〜Ｓ１８において、形式記述情報Ｄｘにおいて数値計算法のアルゴリズムが記述されたアルゴリズム部分の解釈を行う。
【００５４】
まず、制御部１０は、形式記述情報Ｄｘにおけるアルゴリズム部分の記述ＸＤに基づいて、抽象構文木やＤＯＭツリーなどで数値計算法のアルゴリズムを規定する数式を表現する数式情報を生成する（Ｓ１５）。
【００５５】
次に、制御部１０は、生成された数式情報から、反復型数値計算アルゴリズムにおける数列を示す情報を抽出し、数列リストＬ１を生成する（Ｓ１６）。数列リストＬ１は、反復型数値計算アルゴリズムを規定する１つ又は複数の漸化式中の添字付き変数が構成する数列を分類するリストである。図８（ｂ）に、ステップＳ１６で生成される数列リストＬ１の一例を示す。
【００５６】
図８（ｂ）に示すように、数列リストＬ１は、「数列」と、「添字付き変数」と、「変数添字」と、「初期項フラグ」とを関連付けている。初期項は、計算コードＣｏに対する入力として初期値が設定される数列の最初の１つ又は複数の項であり、初期項フラグは、各数列について初期項の入力があるか否かを示すフラグである。
【００５７】
図８（ｂ）は、上述の２次ＡＢ法の例の数列リストＬ１を例示しており、この数列リストＬ１は漸化式の組（１）〜（４）の記述に基づき生成される。数列リストＬ１では、例えば変数添字ｉを有する添字付き変数ｙｉで構成される数列｛ｙｉ｝に関して、計算コード入力リストＬｉｎ（図８（ａ））の初期項ｙ０に応じて、初期項フラグが「ＯＮ」に設定されている。添字付き変数ｔｉ，ｇｉの数列｛ｔｉ｝，｛ｇｉ｝についても同様に管理される（例えば、計算コード入力リストＬｉｎにｇ０等がないため、数列｛ｇｉ｝の初期項フラグは「ＯＦＦ」に設定されている。）。
【００５８】
図７に戻り、制御部１０は、漸化式の組に含まれる添字付き変数の変数添字の正規化を行う（Ｓ１７）。変数添字の正規化では、漸化式中の添字付き変数の変数添字の内で最小の変数添字が「０」に設定されるとともに、同式中の変数添字間の相対的な関係（添字相対値）を維持するように他の添字付き変数の変数添字が設定される。例えば、漸化式（１）は、変数添字の最小値が「ｉ−１」であるため、以下のように正規化される。
ｙ［２］＝ｙ［１］＋ｈ（３ｇ［１］−ｇ［０］）／２      （１’）
【００５９】
次に、制御部１０は、正規化された漸化式の組における数式毎に、添字付き変数を抽出してリスト化し、数式変数リストＬ２を生成する（Ｓ１８）。数式変数リストＬ２は、数値計算法のアルゴリズム部分を記述する漸化式の組の各数式に含まれる別々の変数を分類するリストである。図８（ｃ）に、ステップＳ１８で生成される数式変数リストＬ２の一例を示す。
【００６０】
図８（ｃ）に示すように、数式変数リストＬ２は、「数式」と、複数の変数（「変数１」、「変数２」、…）とを関連付けている。数式変数リストＬ２では、正規化された漸化式の組において、各の数式に含まれる全ての変数が、添字の違いも考慮して、１つずつ記録される。例えば、漸化式（１）については、ステップＳ１７で正規化されていることから（上式（１’）参照）、図８（ｃ）に示すように、４つの変数ｙ［２］、ｙ［１］、ｇ［１］、ｇ［０］が記録されている。
【００６１】
以上のように、制御部１０は形式記述情報Ｄｘ中のアルゴリズム部分の解釈を行う（Ｓ１５〜Ｓ１８）。ステップＳ１５で生成された数式情報は、ステップＳ１６〜Ｓ１８で得られた各種情報に基づき、制御部１０によって適宜、更新される。
【００６２】
これにより、制御部１０は、図３のフローチャートにおける形式記述解釈処理（ステップＳ１）を終了し、数式構造解析処理（ステップＳ２）に進む。
【００６３】
以上の処理によると、反復計算回数に対応する添字付き変数が構成する数列や、数式毎に含まれる添字付き変数など、制御部１０が形式記述情報Ｄｘにおいて反復型数値計算法を表す漸化式を解析するために必要な情報が整理される。
【００６４】
以上の各処理において、制御部１０は、数式情報や計算コード入力リストＬｉｎ，数列リストＬ１，数式変数リストＬ２などの各種情報を生成すると、随時、記憶部１１に保持し、保持した情報を適宜、以降の処理に利用する。
【００６５】
なお、以上の説明では、数式構造解析処理（図３のステップＳ２）よりも前の形式記述解釈処理（ステップＳ１）において各種情報を抽出したが、上記情報を事前に抽出せずに、数式構造解析処理（ステップＳ２）以降の処理において必要の都度、抽出してもよい。
【００６６】
また、以上の説明では、変数添字の種類が「ｉ」の１種類のみの例について説明したが、これに限らず、複数種類の変数添字があってもよい。例えば、偏微分方程式の解法における有限差分法のように空間的方向を示す添字が、反復計算回数（時間方向）に対応する添字とは別に用いられてもよい。また、変数添字について、反復計算処理が２重ループや多重ループになっていてもよい。
【００６７】
また、以上の説明では、数列｛ｙｉ｝，｛ｔｉ｝がそれぞれ１つの初期項ｙ０，ｔ０を有する例について説明したが、初期項は、数列の最初の１項に限らず、複数項であってもよい。例えば、数列｛ｙｉ｝に関して初期項としてｙ０とｙ１の値が入力されてもよい。
【００６８】
また、以上の説明では、２次ＡＢ法にEuler法を組み合わせた一例の式（１）〜（４）が記述された形式記述情報Ｄｘに対する形式記述解釈処理について説明した。形式記述情報Ｄｘには、任意の反復型数値計算を記述することができ、例えば複数の反復型数値計算法を組み合せて記述されてもよい。また、適宜、式（２）のｇｉのように、中間変数を導入して、形式記述情報Ｄｘが記述されてもよい。
【００６９】
２−１−２．数式構造解析処理
図３のフローチャートにおける数式構造解釈処理（ステップＳ２）について、図９〜１１を参照して説明する。図９は、数式構造解析処理を示すフローチャートである。図１０は、数式構造解析処理における数式の展開を説明するための図である。図１１は、数式構造解析処理における二部グラフを説明するための図である。
【００７０】
以下では、解析対象の漸化式が陽的に記述されていても陰的に記述されていても統一的に数式構造を解析可能に実現される数式構造解析処理の一例について説明する。
【００７１】
まず、制御部１０は、数列リストＬ１及び数式変数リストＬ２に基づいて、初期項が入力として設定されない数列に属する添字付き変数の添字の正規化後の値を、数式間で揃えるように、漸化式の組を展開する（Ｓ２１）。本処理は、数列の添字の値が数式間で特定の相対的関係を有する漸化式の組を、数値計算時に使用され得る数列内の別変数が明示された一連の数式に展開するために行われる。
【００７２】
ここで、図１０（ａ）〜（ｃ）を用いて、ステップＳ２１の処理の一例を説明する。本例では、図１０（ａ）の数列リストＬ１Ａに示すように、２つの数列｛ｘｋ｝，｛ｚｋ｝があり、数列｛ｘｋ｝の初期項フラグは「ＯＮ」である一方、数列｛ｚｋ｝の初期項フラグは「ＯＦＦ」である。また、図１０（ｂ）の数式変数リストＬ２Ａに示すように、２つの数式Ｅｑ１，Ｅｑ２がある。数式変数リストＬ２Ａによると、（正規化後の）第１の数式Ｅｑ１は添字付き変数ｘ［１］，ｚ［１］，ｚ［０］間の関係を表し、（正規化後の）第２の数式Ｅｑ２は添字付き変数ｘ［０］，ｚ［１］間の関係を表していることがわかる。
【００７３】
ここで、数式変数リストＬ２Ａからは、第１の数式Ｅｑ１中に存在する数列｛ｚｋ｝に属する添字付き変数ｚ［１］，ｚ［０］のうちで、ｚ［１］が第２の数式Ｅｑ２に関連していることは即座に読み取れるが、ｚ［０］が第２の数式Ｅｑ２に関連していることは字句上からだけでは読み取れない。しかし、第２の数式Ｅｑ２の実際（正規化前）の変数はｚｋ＋１，ｘｋであり（図１０（ａ））、実際の第２の数式Ｅｑ２はｚ［０］にも関連している。そこで、ステップＳ２１の処理において、制御部１０は、第１の数式Ｅｑ１が含む初期項フラグが「ＯＦＦ」である添字付き変数ｚ［０］の添字の値「０」に揃えるように、第２の数式Ｅｑ２を展開する。
【００７４】
具体的には、数式変数リストＬ２Ａ（図１０（ｂ））において第２の数式Ｅｑ２を示すデータを、各添字の値を「−１」だけシフトして複製し、図１０（ｃ）に示すように、二つの数式Ｅｑ２−１，Ｅｑ２−２として数式変数リストＬ２Ａ’に記録する。これに合わせて、制御部１０は、形式記述処理で得られた数式情報等も適宜、更新する。数式（漸化式）Ｅｑ１，Ｅｑ２を展開した一連の数式Ｅｑ１，Ｅｑ２−１，Ｅｑ２−２の内の数式Ｅｑ２−２には、ｚ［０］とともにｘ［−１］も含まれる（図１０（ｃ））。このように、数式変数リストＬ２Ａ’で示される展開後の一連の数式には、各数列中で数値計算時に使用され得る別々の変数が、全て出揃うこととなる。
【００７５】
図９に戻り、制御部１０は、数式変数リストＬ２Ａ’で示される展開後の一連の数式に基づき、例えば各数式と添字付き変数との間の二部グラフを用いて、数式間の変数を介した依存関係を解析する（Ｓ２２）。ステップＳ２１，Ｓ２２の処理により、各数式における添字の値の範囲が解析される。
【００７６】
図１１（ａ），（ｂ）は、ステップＳ２２の処理で用いられる二部グラフの例を示す。図１１（ａ）に示すように、二部グラフは、第１の点群Ｐ１と、第２の点群Ｐ２と、第１及び第２の点群Ｐ１，Ｐ２の間に引かれる線とで構成される。図１１（ａ）に示す二部グラフにおいて、第１の点群Ｐ１は、図１０（ｃ）に示す展開後の数式変数リストＬ２Ａ’における各数式Ｅｑ１，Ｅｑ２−１，Ｅｑ２−２対応する３つの点からなる。また、第２の点群Ｐ２は、数式変数リストＬ２Ａ’において互いに異なる全ての添字付き変数ｘ［１］，ｚ［１］，ｚ［０］，ｘ［０］，ｘ［−１］に対応する５つの点からなる。
【００７７】
制御部１０は、ステップＳ２２の処理において、数式変数リストＬ２Ａ’に基づいて、図１１（ａ）の二部グラフを配列等で生成する。図１１（ａ）の二部グラフにおいて、第１及び第２の点群Ｐ１，Ｐ２の間の線は、数式Ｅｑ１，Ｅｑ２−１，Ｅｑ２−２と、各数式を構成する変数ｘ［１］，ｚ［１］，ｚ［０］，ｘ［０］，ｘ［−１］とを結んだものである。例えば、数式Ｅｑ１には、添字付き変数ｘ［１］，ｚ［１］，ｚ［０］が現れるため（図１０（ｃ）参照）、数式Ｅｑ１に対応する点と、上記変数に対応する３点との間に線が引かれている。図１１（ａ）の二部グラフによると、例えば、第１の数式Ｅｑ１におけるｘ［１］は、数式Ｅｑ２−１，Ｅｑ２−２を介して、ｘ［０］及びｘ［−１］（計算ステップが１ステップ前及び２ステップ前のデータ）を用いて計算されることが判明する。
【００７８】
次に、制御部１０は、形式記述処理で得た計算コード入力リストＬｉｎ等の情報に基づいて、図１１（ｂ）に示すように、第１の点群Ｐ１に、初期項を含む定数添字付き変数に対応する点（ｘ０，ｘ１）を追加する。これとともに、制御部１０は、第１の点群Ｐ１に追加した点と、追加した点に対して第２の点群Ｐ２内で所定の関係を有する点（例えば、追加した点に対応する定数添字付き変数と同じ数列に属する変数添字付き変数に対応する点）との間に線を引くように、二部グラフを更新する。図１１（ｂ）の例では、第１の点群Ｐ１に追加した２点の初期項ｘ０，ｘ１と同じ数列｛ｘｋ｝の変数ｘ［１］，ｘ［０］，ｘ［−１］の内の、添字の値が小さい２つ（ｘ［−１］，ｘ［０］）の点と、追加した２点（ｘ０，ｘ１）との間に線が引かれている。
【００７９】
次に、制御部１０は、更新した二部グラフの最大マッチングを求め、最大マッチングに基づき、展開後の数式の中で連立方程式として解く必要がある数式の集合を示す連立方程式情報と、計算を行う順序に関する数式間の半順序関係を示す計算順序情報とを導出する（Ｓ２３）。計算順序情報が示す数式間の半順序関係は、数値計算時に順番どおりに計算する必要がある数式間の順序の関係だけでなく、独立に計算可能な数式の組み合わせの関係も含む。上記の二部グラフの最大マッチングを求めることにより、連立方程式情報及び計算順序情報は既存の手法で導出することができる（例えば、非特許文献２，３参照）。
【００８０】
次に、制御部１０は、漸化式中の添字の範囲の解析結果に基づき、変数添字付き変数を算出可能な数式を抽出することにより、ループ部分情報を取得する（Ｓ２４）。ループ部分情報は、反復計算において計算の繰り返しになるループ部分を示す情報である。例えば、制御部１０は、数式変数リストＬ２Ａ，Ｌ２Ａ’、更新された二部グラフ（図１０（ｂ），１１（ｂ）参照）及び数式間の半順序関係に基づき、一連の数式において数式Ｅｑ１、Ｅｑ２をループ部分情報として抽出する。なお、ループ部分情報は、ステップＳ２３で計算順序情報と同時に抽出されてもよい。ループ部分情報は、既存の手法（例えば、非特許文献２，３参照）を利用し、半順序関係から添字をずらすことで反復計算できる部分として抽出できる。
【００８１】
次に、制御部１０は、連立方程式情報、計算順序情報及びループ部分情報並びに形式記述解釈処理で得た各種情報に基づき、種々のプログラミング言語に変換可能な中間形式で数値計算法を表す中間データを生成する（Ｓ２５）。中間データは、連立方程式情報、計算順序情報及びループ部分情報を示す情報を含む。
【００８２】
具体的に、制御部１０は、展開後の数式中で、連立方程式情報が示す数式の集合を特定し、中間形式において、特定した数式の集合を連立方程式として解くことを設定する。さらに、制御部１０は、計算順序情報に基づく計算順序において、一連の数式の半順序関係を中間データとして保存する。さらに、制御部１０は、ループ部分情報及び数列リストＬ１に基づき、反復計算のためのループ部分において、繰り返し前後の計算ステップ間の各変数を対応付け、中間データを生成する。
【００８３】
これにより、制御部１０は、図３のフローチャートにおける数式構造解析処理（ステップＳ２）を終了し、コード生成処理（ステップＳ３）に進む。
【００８４】
以上の処理によると、数列中の変数間で共通の漸化式の組を展開することで、その数式構造を解析でき、反復型数値計算法を実現する中間データを得ることができる。また、計算環境に依存しない中間形式の中間データを生成することで、種々の計算環境に応じた計算コードＣｏが生成可能になる。
【００８５】
例えば、ステップＳ２３で連立方程式として解く必要がある数式の集合が存在した場合、制御部１０は、コード生成処理（図３のステップＳ３）において中間データに基づいて、非線形方程式の求根アルゴリズムを用いてこの漸化式を解くことを判断できる。この場合、例えば、制御部１０は、予め記憶部１１に記録された非線形方程式の求根アルゴリズムを記述する情報を参照して、上記の漸化式を当該求根アルゴリズムの計算対象とするように、適宜、計算コードに組み込む。
【００８６】
また、ステップＳ２１の処理により漸化式が陰的であっても陽的であっても同様に展開され、陰的な漸化式を含む漸化式の組に対する数値計算が可能な中間データを生成することができる。
【００８７】
以上の説明では、漸化式が陽的に記述されていても陰的に記述されていても統一的に数式構造を解析可能な方法について説明した。数式構造解析処理では、添字が最大の添字付き変数に対して陽的に記述可能である陽的な漸化式と、陽的に記述不可能である陰的な漸化式とについて、異なる方法で数式構造が解析されてもよい。以下、この方法の一例を説明する。
【００８８】
まず、数式構造解析処理の簡略化のため、形式記述情報Ｄｘにおいて、陽的な漸化式は下記の式（５）のように記述され、陰的な漸化式は下記の式（６）のように記述されることとする（Ｆは任意の関数）。
ｘｉ＋１＝Ｆ（ｘｉ，ｘｉ−１，…）        （５）
Ｆ（ｘｉ＋１，ｘｉ，…）＝０            （６）
【００８９】
この場合、制御部１０は、例えば、形式記述情報Ｄｘにおける漸化式の組の中に、右辺が「０」の数式があるか否かによって、陰的な漸化式があるか否かを判断することができる。
【００９０】
ここで、式（５）の陽的な漸化式の記述によると、右辺を計算することにより左辺の変数の値が算出できることが明らかである。このため、漸化式の展開（ステップＳ２１）において、制御部１０は、陽的な漸化式に対しては右辺の中で、添字付き変数の添字の正規化後の値を合わせる必要がある変数があるか否かを判断すればよい。また、制御部１０は、特定の変数に揃えるように数式を展開する必要がある場合、陽的な漸化式については左辺の変数が特定の変数と添字の値だけずれている場合にのみ、添字の値をシフトするように展開すればよい。
【００９１】
例えば、図７の形式記述解釈処理（図３のＳ１）の説明で用いた２次ＡＢ法の例では、数式構造解析処理（Ｓ２）において、陽的な漸化式の組（１）〜（４）は以下のように展開される。
ｙｉ＋１＝ｙｉ＋ｈ（３ｇｉ−ｇｉ−１）／２      （１）
ｇｉ＝ｆ（ｙｉ，ｔｉ）                      （２−１）
ｇｉ−１＝ｆ（ｙｉ−１，ｔｉ−１）              （２−２）
ｔｉ＋１＝ｔｉ＋ｈ                            （３）
ｙ１＝ｙ０＋ｈｇ０                            （４）
【００９２】
上式では、式（１）の右辺のｇｉ，ｇｉ−１に合わせて、左辺がｇｉであった式（２）が式（２−１），（２−２）に展開されている。以上のように、陽的な漸化式に対しては右辺と左辺との役割を設定して展開を行うことで、数値計算に用いることがない無駄な展開を行わずに、効率的に数式構造を解析することができる。
【００９３】
図１２は、式（１）〜（４）の例において、数式構造解析処理で解析される諸情報（計算順序情報Ｄｏｒｄｅｒ及びループ部分情報Ｄｌｏｏｐ）を説明するための図である。同図が示す情報は、展開後の式（１），（２−１），（２−２），（３），（４）に対してステップＳ２３，Ｓ２４と同様に最大マッチングを求めることなどによって得られる。
【００９４】
図１２では、式（１）〜（４）における正規化後の各変数の依存関係を矢印で示している。例えば、ｙ［２］の計算（すなわち式（１）の計算）は、ｙ［１］，ｇ［１］，ｇ［０］に依存していることを表している。これらの矢印で示された（枝分かれを含む）順番が、数式間の半順序関係すなわち計算順序情報Ｄｏｒｄｅｒに対応している。例えば、ｙ［２］の計算（式（１））の前に、ｇ［１］の計算（式（２−１））及びｇ［０］の計算（式（２−２））をする必要があることなどが計算順序情報Ｄｏｒｄｅｒに含まれている。
【００９５】
また、図１２において、破線で囲んだ範囲がループ部分情報Ｄｌｏｏｐに対応している。つまり、破線で囲まれた範囲に含まれる計算が反復計算において繰り返されることとなる。ループ部分情報Ｄｌｏｏｐは、ｙ［１］，ｔ［０］，ｙ［０］が定数添字変数ｙ１，ｇ０，ｔ０、或いは以前の計算ステップにおいて計算された変数から添字をずらしたものである、ということから抽出することができる。
【００９６】
以上の解析結果を含む中間データが、制御部１０により、ステップＳ２５で生成される。
【００９７】
また、以上のように生成された中間データや各処理における解析結果は、随時、記憶部１１に保持され、以降の処理に適宜、利用される。
【００９８】
２−１−３．コード生成処理
図３のフローチャートにおけるコード生成処理（ステップＳ３）について、図１３〜１６を参照して説明する。コード生成処理（Ｓ３）では、数式構造解析処理（Ｓ２）により生成された中間データを計算環境に応じて最適化して、最終的な計算コードＣｏを生成する。以下では、手続型プログラミング言語で計算コードが生成され、最適化される例について説明する。
【００９９】
図１３は、コード生成処理を示すフローチャートである。図１４，１５は、本実施形態における形式記述情報の別例を示す図である。図１６は、コード生成処理における最適化前後の計算コードを示す図である。
【０１００】
以下では、コード生成処理について、入力の数値計算法が、非線形連立方程式の求根アルゴリズムの一種の割線法である例で説明する。図１４，１５に、数値計算法が割線法である例の形式記述情報Ｄｘａを示す。図１４，１５に示す形式記述情報Ｄｘａには、割線法を規定する以下の漸化式を示す情報が記述されている。
ｘｉ＋１＝ｘｉ−ｆ（ｘｉ）（ｘｉ−ｘｉ−１）／（ｆ（ｘｉ）−ｆ（ｘｉ−１））      （７）
【０１０１】
以下では、図１４，１５の形式記述情報Ｄｘａに基づき図３のステップＳ１，Ｓ２の形式記述解釈処理及び数式構造解析処理が行われた後に、ステップＳ３のコード生成処理が開始したこととする。
【０１０２】
図１３のフローチャートにおいて、まず、制御部１０は、計算環境情報Ｄｅ（図１）を取得する（Ｓ３０）。
【０１０３】
計算環境情報Ｄｅは、計算コードＣｏが実行される計算環境に応じた種々の指定を含む情報である。種々の指定は、例えば、計算コードＣｏを記述する手続き型プログラミング言語の指定や、計算コードＣｏの実行時に（最大）幾つに並列化して処理を実行するかの指定（並列化の指定）、メモリの使用と計算処理のいずれを優先するのかの指定（メモリ使用／計算処理の優先指定）を含む。並列化の指定は、例えば計算環境におけるＣＰＵ数や分散ネットワーク環境を考慮して設定される。メモリ使用／計算処理の優先指定は、例えばメモリの搭載量及びメモリのアクセス速度、並びに搭載ＣＰＵの性能に対するメモリ性能の相対評価等を考慮して設定される。
【０１０４】
なお、計算環境情報Ｄｅは、計算コードＣｏが実行される計算環境の各種スペックを示す情報であってもよく、制御部１０が計算環境情報Ｄｅにおける各種スペックに基づき上記の指定を判断するようにしてもよい。また、計算環境情報Ｄｅは、図３の計算コード生成処理の開始前からあらかじめ記憶部１１に格納されていてもよいし、図３のステップＳ２の数式構造解析処理後に、ユーザによって操作部１２（図２）から入力されてもよい。
【０１０５】
次に、制御部１０は、数式構造解析処理（Ｓ２）で生成された中間データ内の記述を、計算環境情報Ｄｅにおいて指定された手続き型プログラミング言語の記述に置き換えて、（最適化前の）計算コードＣｏ’を生成する（Ｓ３１）。図１６（ａ）に、ステップＳ３０の処理で生成される計算コードＣｏ’の一例を示す。
【０１０６】
次に、制御部１０は、計算環境情報Ｄｅにおいて、メモリ使用／計算処理の優先指定のうちで、メモリ使用の優先指定がされているか否かを判断する（Ｓ３２）。
【０１０７】
制御部１０は、メモリ使用の優先指定がされている場合（Ｓ３２でＹｅｓ）、数式構造解析処理の解析結果に基づいて、計算コード中で、計算結果が共通する共通部分に対する一時保存コードを生成する（Ｓ３３）。一時保存コードは、中間変数等を用いて共通部分の計算結果を一時的に保存するコードである。
【０１０８】
具体的には、制御部１０は、中間データを構成する（展開後の）全数式間または各数式の内部で、共通の部分式を検出する。共通の部分式は、添字付き変数の添字相対値まで完全に一致する部分式だけでなく、添字付き変数の添字相対値が全体的にシフトしている部分式も含む。制御部１０は、検出した部分式の評価値が格納される中間変数を有する一時保存コードを生成する。一時保存コードは、添字相対値がシフトしている共通部分に対しては、シフト幅の分、反復計算の計算ステップを跨いで一時保存するように記述される。
【０１０９】
例えば、式（７）においては２箇所のｆ（ｘｉ）が共通しており、図１６（ａ）に示すように、最適化前の計算コードＣｏ’上では、１１行目において２回、同じ評価値となる関数呼び出しを行っている。また、式（７）には、ｆ（ｘｉ）から添字の値がシフトしたｆ（ｘｉ−１）があり、この評価値は、反復計算の計算ステップが１つ前のｆ（ｘｉ）の評価値と同一である。このため、図１６（ｂ）に示すように、最適化後の計算コードＣｏ上でｆ（ｘｉ）の関数呼び出しの戻り値を一時保存する一時保存コード（１９行目）が用いられることで、使用メモリ量を増やす代わりに計算量を削減することができる。
【０１１０】
一方、実行環境や関数の内容によっては、関数呼び出しの方が有利な場合もある。このため、制御部１０は、計算環境情報Ｄｅにおいて、メモリ使用の優先指定がされていない場合（Ｓ３２でＮｏ）、ステップＳ３３の処理を行うことなく、ステップＳ３４の処理に進む。
【０１１１】
次に、制御部１０は、計算環境情報Ｄｅにおいて、（２並列以上の）並列化の指定がされているか否かを判断する（Ｓ３４）。
【０１１２】
制御部１０は、並列化の指定がされている場合（Ｓ３４でＹｅｓ）、数式構造解析処理の解析結果に基づいて、数式上で並列に計算可能な部分を判別し、計算コードを並列化する（Ｓ３５）。例えば、計算順序情報において独立に計算できる数式同士や、数式情報中の数式内で「＋」や「−」を介して結合している部分式同士は並列に計算可能であり、制御部１０は、計算順序情報や数式内の算術演算子に基づき並列計算可能な部分式を判別する。また、制御部１０は、数式構造解析処理で抽出された連立方程式情報に基づき、中間データを構成する全数式間で、連立方程式として解く必要のない、独立に計算可能な数式を判別する。制御部１０は、判別した部分式や数式について、計算環境情報Ｄｅにおいて指定された並列化の限度に応じて、計算コードの並列化を行う。
【０１１３】
一方、制御部１０は、計算環境情報において並列化の指定がない場合（Ｓ３４でＮｏ）、ステップＳ３５の処理を行うことなく、ステップＳ３６の処理に進む。
【０１１４】
次に、制御部１０は、計算環境情報Ｄｅに基づいて、計算コードＣｏ中で省略可能な部分（不要部分）の整理を行う（Ｓ３６）。例えば、制御部１０は、計算コードＣｏ中で使い回すことのない不要な中間変数の定義と、その中間変数を介した代入文を１文にまとめる。例えば、図１６（ａ）では、初期項ｘ０，ｘ１を示す入力変数ｉ０，ｉ１の値は数値計算開始後に参照されないため、この値を保存する必要はなく上書き可能である。制御部１０は、数式構造解析処理（Ｓ２）で生成された中間データに基づき変数ｉ０，ｉ１が不要（省略可能）と判断する。制御部１０は、代入文ｘ０＝ｉ０，ｘ＝ｉ１（図１６（ａ）の７，８行目）を削除して、直接、入力変数ｘ０，ｘ１を反復計算に使用するように、計算コードＣｏを整理する（図１６（ｂ）３，１４−１７行目参照）。なお、計算環境におけるプログラミング言語やコンパイラによっては、引数変数を使い回さない方が有利な場合もあるため、本処理における不要部分の整理方法は、計算環境情報Ｄｅに応じて適宜、設定される。
【０１１５】
次に、制御部１０は、計算順序情報と計算環境情報Ｄｅに基づいて、計算順序の最適化を行う（Ｓ３７）。具体的には、計算順序に任意性のある部分についてメモリ局所参照性やパイプライン処理、多重処理、並列処理などを考慮して、計算コードＣｏ上の計算順序の最適化を行う。本処理により、計算順序の最適化において実現可能な全ての計算順序について実行性能を比較評価するような手間を掛けることなく、既存の探索手法を用いて効率的に最適または準最適な計算順序を得ることができる。
【０１１６】
次に、制御部１０は、計算環境情報Ｄｅに基づき、計算コードＣｏ中で、形式記述情報Ｄｘａの反復条件に対応する反復計算の終了判定を行う箇所を最適な箇所に変更する（Ｓ３８）。例えば、図１６（ｂ）に示すように、反復計算のループ（１１〜２０行）中で終了条件の判定を可能な限り前方（１２行目）で行うことで、無駄な計算を削減することができる。なお、計算環境における数値計算の実行条件によっては、最適化前の位置（図１６（ａ）の１０〜１７行中の１３行目）での評価が有利な場合もあるため、本処理の最適化方法は、計算環境情報Ｄｅで指定される条件に応じて適宜、設定される。
【０１１７】
その後、制御部１０は、以上の処理により最適化された計算コードＣｏを記憶部１１等に出力し（Ｓ３９）、本処理を終了する。
【０１１８】
以上の処理により、計算環境情報Ｄｅで指定される種々の計算環境に応じて、数式構造解析処理の解析結果に基づき自動的に判定可能な範囲で最適化がなされた計算コードＣｏを、効率的に得ることができる。また、計算環境情報Ｄｅにおいて並列化の指定がある場合には、その指定に応じて並列化された計算コードＣｏが生成される。
【０１１９】
３．まとめ
以上のように、本実施形態に係る計算コード生成装置１は、計算対象の関数に対する数値計算が実行される計算コードＣｏを自動生成する。計算コード生成装置１は、記憶部１１と、制御部１０とを備える。記憶部１１は、計算対象の関数に対する数値計算法を数式で宣言的に記述する形式記述情報Ｄｘを格納する。制御部１０は、記憶部１１に格納された形式記述情報Ｄｘに基づく演算処理を行う。制御部１０は、形式記述情報Ｄｘに基づいて、形式記述情報Ｄｘにおける数式及び数式が含む変数の関係を解析する（Ｓ２）。制御部１０は、解析結果に基づいて、所定のプログラミング言語により計算対象の関数に対する数値計算を実行する計算コードＣｏを生成する（Ｓ３）。形式記述情報Ｄｘは、反復型数値計算法を漸化式で規定する情報を含む（図４〜６等参照）。
【０１２０】
本実施形態に係る計算コード生成装置１によると、形式記述情報Ｄｘに記述された漸化式に基づき、あらゆる反復型数値計算法を実行する計算コードＣｏが生成可能である。また、入力の形式記述情報Ｄｘでは、数値計算法全般を同一形式で統一的に記述することができる。これにより、数値計算プログラムの生産性を向上することができる。
【０１２１】
また、本実施形態において、制御部１０は、形式記述情報Ｄｘにおける漸化式を、漸化式に含まれる（数列中の）変数を介して関連する一連の数式に展開する（Ｓ２３）。数列は、反復型数値計算法における反復計算回数に対応する添字を有する変数の列で構成される。これにより、計算コード生成装置１は、数列中の変数間で共通の漸化式の組における数式と変数との関係を解析でき、反復型数値計算法を実現するための情報を得ることができる。
【０１２２】
また、本実施形態において、制御部１０は、形式記述情報Ｄｘにおける数式毎に、数式に含まれる（数列中の）変数を分類して、漸化式を示す情報を解釈する（Ｓ１）。これにより、形式記述情報Ｄｘにおいて反復型数値計算法を表す漸化式を解析するために必要な情報が整理される。
【０１２３】
また、本実施形態において、制御部１０は、計算コードＣｏが実行される計算環境を示す計算環境情報Ｄｅを取得する。制御部１０は、計算環境情報Ｄｅに応じて最適化して、計算コードＣｏを出力する。これにより、計算環境情報Ｄｅで指定される種々の計算環境に応じて、数式から自動的に判定可能な範囲で最適化がなされた計算コードＣｏを、効率的に得ることができる。
【０１２４】
また、本実施形態において、形式記述情報Ｄｘにおける漸化式（１）〜（４）は、少なくとも３回にわたる反復計算回数にそれぞれ対応する少なくとも３つの変数ｙｉ＋１，ｙｉ，ｇｉ，ｇｉー１を含む。これにより、ｉ回目の計算ステップでｙｉ＋１の評価値を計算するために、ｇｉ−１において過去の計算結果を用いるのか、再計算するのか、適宜選択して計算コードＣｏを生成することができる。
【０１２５】
また、本実施形態において、生成部Ｓ３によって生成された計算コードＣｏの外部から、数学的問題（問題定義関数関数ｆ）の具体形が指定される。これにより、生成された計算コードＣｏは、種々の問題定義関数を計算するための数値計算ライブラリ等として使用することができる。
【０１２６】
（他の実施形態）
上記の実施形態１では、計算コード生成装置１がＰＣなどの情報処理装置で構成される例について説明したが、これに限らず、例えば、計算コード生成装置１はＡＳＰサーバなどのサーバ装置であってもよい。例えば、計算コード生成装置１は、ネットワークを介して入力された形式記述情報Ｄｘや計算環境情報Ｄｅをネットワークインタフェース（取得部の一例）により取得して、計算コード生成処理を実行してもよい。また、計算コード生成装置１は、計算コード生成処理において生成した計算コードＣｏを、ネットワークを介して送信してもよい。
【０１２７】
また、計算コードの自動生成方法は、クラスタＰＣやクラウドコンピューティングなどにおいて実行されてもよい。
【０１２８】
また、計算コードの自動生成方法が実現されるプログラムは、例えば、ＣＡＥ，ＣＡＤソフトウェアに組み込まれてもよいし、数式処理ソフトウェア、数値解析ソフトウェアにおいて利用されてもよい。
【０１２９】
また、上記の実施形態１では、２次ＡＢ法や割線法などの一種の反復型数値計算法が記述された形式記述情報Ｄｘ，Ｄｘａの例について説明した。反復型数値計算法においては、その内部で別の反復型数値計算法を用いることがあり、形式記述情報Ｄｘは、数値計算法及びその一部を部分的に交換可能にするように、モジュール化して記述されてもよい。例えば、数値計算法及びその内部で用いられる数値計算法の種別毎に共通インタフェースを定義した形式言語を用いて形式記述情報Ｄｘが記述されてもよい。
【０１３０】
例えば、形式記述情報Ｄｘにおいて、非線形方程式の求根アルゴリズムの記述中で、数値微分を用いることを抽象的に記述し、具体的に利用する数値微分法は形式記述情報Ｄｘの外部から指定するようにしてもよい。また、数値計算法の内部で別の数値計算法を利用することについて、形式記述情報Ｄｘ上で明示的には記述せず、計算コード生成処理において計算コード生成装置１が判断してもよい。図１７〜１９を用いてこれを説明する。
【０１３１】
図１７は、数値計算法がニュートン法（内部で数値微分法を用いる非線形方程式の求根アルゴリズムの一種）である例の形式記述情報Ｄｘｂを示す。形式記述情報Ｄｘｂには、図１７に示すように、アルゴリズム部分を形式的に記述した第１の参照情報Ｄｘｂ１（図１８）と、反復計算の終了条件を形式的に記述した第２の参照情報Ｄｘｂ２（図１９）とを参照する記述がある。計算コード生成装置１の制御部１０は、図７の形式記述解釈処理のステップＳ１１において、上記の記述に基づき各参照情報Ｄｂ１，Ｄｂ２を読み出して、形式記述情報Ｄｘｂに合成する。
【０１３２】
また、図１８に示すように、第１の参照情報Ｄｘｂ１には、ニュートン法を規定する以下の漸化式を示す情報が記述されている。
ｘｉ＋１＝ｘｉ−ｆ（ｘｉ）／ｆ’（ｘｉ）      （８）
【０１３３】
上式（８）には、関数ｆ（ｘｉ）の導関数ｆ’（ｘｉ）（＝ｄｆ（ｘｉ）／ｄｘ）が含まれている。制御部１０は、例えば図７のステップＳ１１において、第１の参照情報Ｄｘｂ１中の導関数の記述（図１８）に基づき数値微分法を用いることを判断する。この場合、制御部１０は、記憶部１１に予め記憶された所定の数値微分法を示す形式記述情報を読み出して、形式記述情報Ｄｘｂに合成する。
【０１３４】
以上のように数値計算法の部分処理や部分解法（例えば収束判定法や数値微分）の記述をモジュール化することにより、入力する数値計算法の部分的変更を容易化することができる。また、多数の部分解法を組み合わせる複雑な数値計算法も容易に記述可能になる。
【０１３５】
また、上記の実施形態１では、計算対象の数学的問題を定める問題定義関数の具体形は、外部から計算コードＣｏに対する入力として指定された。しかし、これに限らず、問題定義関数（数学的問題）の具体形は、計算コードＣｏの内部に組み込まれてもよい。この場合、図３の計算コード生成処理において、問題定義関数の具体形は、形式記述情報Ｄｘとともに、ステップＳ１において制御部１０に入力される。また、ステップＳ２の数式構造解析処理では、図９のステップＳ２１において、問題定義関数の方程式系についても数式が展開される。
【０１３６】
例えば、常微分方程式の初期値問題における問題定義関数の具体形として、生体機能モデルの一種のFitzHugh-Nagumoモデルは、次式のように表される。
ｄｖ／ｄｔ＝ｖ−ｖ３／３−ｗ＋Ｉ（＝ｆ１（ｙ））          （９）
ｄｗ／ｄｔ＝（ｖ＋ａ−ｂ×ｗ）／τ（＝ｆ２（ｙ））      （１０）
【０１３７】
上式（９），（１０）において、Ｉ，ａ，ｂ，τは定数変数である。また、関数変数ｆ及び引数変数ｙは、２成分ベクトルであり、例えばｙ＝（ｖ，ｗ）Ｔである（Ｔは転置記号）。上式（９），（１０）を示す情報が形式記述情報Ｄｘとともに入力されることで、図９のステップＳ２１において、例えば式（２）は、次式のように展開される。
ｇｉ＝（ｖｉ−ｖｉ３／３−ｗｉ＋Ｉ，（ｖｉ＋ａ−ｂ×ｗｉ）／τ）Ｔ      （２’）
【０１３８】
上式（２’）のように、問題定義関数の方程式系についても展開された一連の数式に対して図９のステップＳ２２以降の解析が行われることで、図３のステップＳ３の計算コード生成処理では、問題定義関数の具体形に応じた最適化を行うことができる。
【０１３９】
また、上記の実施形態１では、図９の数式構造解析処理のステップＳ２２において、二部グラフを用いて展開後の数式間の変数を介した依存関係を解析し、連立方程式情報、計算順序情報及びループ部分情報を取得する例について説明したが、二部グラフを用いなくてもよい。展開後の数式構造について、既存手法の組み合わせにより様々な解析方法で解析することで、上記の各種情報を取得してもよい。
【産業上の利用可能性】
【０１４０】
本発明は、数値計算を利用するあらゆる産業分野に適用可能であり、例えば、ＣＡＥ，ＣＡＤソフトウェアへの組み込みや、数式処理ソフトウェア、数値解析ソフトウェアにおける利用が可能である。
【符号の説明】
【０１４１】
１    計算コード生成装置
１０    制御部
１１    記憶部
Ｄｘ，Ｄｘａ，Ｄｘｂ    形式記述情報
Ｄｅ    計算環境情報
Ｃｏ    計算コード
【図１】
2017111749.tif 000003
【図２】
2017111749.tif 000004
【図３】
2017111749.tif 000005
【図４】
2017111749.tif 000006
【図５】
2017111749.tif 000007
【図６】
2017111749.tif 000008
【図７】
2017111749.tif 000009
【図８】
2017111749.tif 000010
【図９】
2017111749.tif 000011
【図１０】
2017111749.tif 000012
【図１１】
2017111749.tif 000013
【図１２】
2017111749.tif 000014
【図１３】
2017111749.tif 000015
【図１４】
2017111749.tif 000016
【図１５】
2017111749.tif 000017
【図１６】
2017111749.tif 000018
【図１７】
2017111749.tif 000019
【図１８】
2017111749.tif 000020
【図１９】
2017111749.tif 000021