JP6601222B2 - 行列演算プログラム、行列分割方法、及び並列処理装置 - Google Patents

行列演算プログラム、行列分割方法、及び並列処理装置 Download PDF

Info

Publication number
JP6601222B2
JP6601222B2 JP2016000151A JP2016000151A JP6601222B2 JP 6601222 B2 JP6601222 B2 JP 6601222B2 JP 2016000151 A JP2016000151 A JP 2016000151A JP 2016000151 A JP2016000151 A JP 2016000151A JP 6601222 B2 JP6601222 B2 JP 6601222B2
Authority
JP
Japan
Prior art keywords
row
matrix
sparse matrix
rows
processing
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2016000151A
Other languages
English (en)
Other versions
JP2017122950A (ja
Inventor
香壱 清水
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Fujitsu Ltd
Original Assignee
Fujitsu Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Fujitsu Ltd filed Critical Fujitsu Ltd
Priority to JP2016000151A priority Critical patent/JP6601222B2/ja
Priority to US15/372,921 priority patent/US20170192818A1/en
Publication of JP2017122950A publication Critical patent/JP2017122950A/ja
Application granted granted Critical
Publication of JP6601222B2 publication Critical patent/JP6601222B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F9/00—Arrangements for program control, e.g. control units
    • G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
    • G06F9/46—Multiprogramming arrangements
    • G06F9/50—Allocation of resources, e.g. of the central processing unit [CPU]
    • G06F9/5061—Partitioning or combining of resources
    • G06F9/5066—Algorithms for mapping a plurality of inter-dependent sub-tasks onto a plurality of physical CPUs
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F9/00—Arrangements for program control, e.g. control units
    • G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
    • G06F9/46—Multiprogramming arrangements
    • G06F9/48—Program initiating; Program switching, e.g. by interrupt
    • G06F9/4806—Task transfer initiation or dispatching
    • G06F9/4843—Task transfer initiation or dispatching by program, e.g. task dispatcher, supervisor, operating system
    • G06F9/4881—Scheduling strategies for dispatcher, e.g. round robin, multi-level priority queues
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F9/00—Arrangements for program control, e.g. control units
    • G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
    • G06F9/46—Multiprogramming arrangements
    • G06F9/52—Program synchronisation; Mutual exclusion, e.g. by means of semaphores

Landscapes

  • Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Complex Calculations (AREA)

Description

本発明は、行列演算プログラム、行列分割方法、及び並列処理装置に関する。
構造物の設計や電子回路の設計などを行う際には、構造解析、流体解析、電磁場解析などの解析手法に基づく数値シミュレーションが利用される。解析対象の物理法則を規定する方程式の多くは連続的な物理量を扱う偏微分方程式で表現される。そのため、数値シミュレーションの際、偏微分方程式の求解問題は、有限要素法（ＦＥＭ：Finite Element Method）などの離散化手法を用いて行列方程式の求解問題に置き換えられる。
上述した行列方程式の係数行列は、次元が大きく、行列成分のほとんどがゼロとなる大型の疎行列となる。そのため、演算負荷及び使用メモリ量を低減するために、行列方程式の解を修正しながら繰り返し計算によって正しい解を求める反復法が利用される。反復法としては、例えば、ＣＧ（Conjugate Gradient）法、ＢｉＣＧ（Bi-Conjugate Gradient）法、ＣＲ（Conjugate Residual）法、ＣＧＳ（Conjugate Gradient Squared）法、ＩＣＣＧ（Incomplete Cholesky Conjugate Gradient）法などがある。
なお、有限要素法や反復法を利用して磁場解析や構造解析を高速に実行する方法が開示されている。また、連立一次方程式の解法に関し、直接法の１つであるコレスキー法と反復法の１つであるＩＣＣＧ法とを組み合わせて用いる方法が開示されている。この方法では、剛性マトリクスに含まれる非ゼロ成分の数がコンピュータのメモリに収まる範囲内かどうかに応じてコレスキー法とＩＣＣＧ法とを使い分けている。
特開２０１０−１２２８５０号公報 特開２００５−２０７９００号公報 特開平５−７３５２７号公報 特開２０１２−２０４８３５号公報
反復法を利用する場合、係数行列を複数の行群（行の集合）に分け、行群毎にプロセス（演算処理の実行単位）を割り当てて複数のプロセスを並列に実行することができる。係数行列が対角成分とその周辺だけに非ゼロ成分を持つバンドマトリクスの場合、プロセス間における処理負荷の偏りは小さい。一方、他の行に比べて非ゼロ成分の数が突出して多い行を含む係数行列の場合、非ゼロ成分の多い行に割り当てられたプロセスの処理が律速となり、並列処理全体の速度低下をもたらす。
例えば、導体の一部に電圧を印加して導体に流れる電流に起因して発生する磁場を有限要素法及び反復法を用いて解析する場合、各有限要素におけるベクトルポテンシャルの成分及び電流の成分が未知数となる。ベクトルポテンシャルの成分に対応する係数行列の領域は有限要素間の接続関係で規定できるため、対角成分とその周辺が非ゼロとなる。
他方、電流の成分に対応する係数は多くの有限要素において非ゼロとなるため、電流の成分に対応する行列要素の各行は多くの非ゼロ成分を含むことになる。構造解析の場合は、拘束条件を付加した場合に非ゼロ成分を多く含む行が現われる。
並列処理の場合、あるプロセスの実行結果が他のプロセスで利用される。そのため、処理が遅延しているプロセスがあると、他のプロセスに処理結果がわたらずに処理が遅延する。つまり、処理が遅延しているプロセスが律速となり全体の処理速度が低下する。並列に実行できるプロセスの数（並列数）が増えるにつれて処理能力がリニアに増加するのが理想的である。しかし、律速が生じるような条件下では、ある並列数を超えると並列数が増加しても処理性能がほとんど向上しなくなる。
１つの側面によれば、本開示の目的は、並列処理のスケーラビリティを向上させることが可能な行列演算プログラム、行列分割方法、及び並列処理装置を提供することにある。
一態様によれば、係数行列として疎行列を含む行列方程式の演算に関する処理を実行するコンピュータが、記憶部から、疎行列の各行に含まれる非ゼロ成分の多さを判定するための閾値を取得し、疎行列の中で閾値より非ゼロ成分の数が大きい第１の行を特定し、当該第１の行を複数の第２の行に分割して疎行列を拡大し、拡大した疎行列を複数の行群に分けて、処理の実行単位となるプロセスを行群毎に割り当てる行列分割方法が提供される。
並列処理のスケーラビリティを向上させることが可能になる。
第１実施形態に係る並列処理装置の一例を示した図である。 第２実施形態に係るハードウェア（単一装置）の一例を示した図である。 第２実施形態に係るハードウェア（複数装置）の一例を示した図である。 第２実施形態に係る情報処理装置が有する機能の一例を示したブロック図である。 行列方程式の構造及びＩＣＣＧ法の擬似コードを示した図である。 行列方程式の構造及び並列ＩＣＣＧ法の擬似コードを示した図である。 第２実施形態に係る行列拡張方法の一例を示した図である。 第２実施形態に係る領域分割方法の一例を示した図である。 第２実施形態に係る並列ＩＣＣＧ法の擬似コードを示した図である。 第２実施形態に係る情報処理装置による処理の流れを示した第１の図である。 第２実施形態に係る情報処理装置による処理の流れを示した第２の図である。 第２実施形態に係る行列情報のデータ構造を示した図である。 第２実施形態に係る通信用情報のデータ構造を示した図である。 第２実施形態の一応用例に係る係数行列（節点の接続関係）の作成方法を示した図である。 第２実施形態の一応用例に係る係数行列（電流未知数を追加）の作成方法を示した図である。 第２実施形態の一応用例に係る行列拡張方法及び領域分割方法の一例を示した図である。 第２実施形態の一応用例に係る行列情報の設定例を示した図である。 第２実施形態の一応用例に係る列ベクトルの非ゼロパターンを示した図である。 第２実施形態の一応用例に係るＣＰＵ間のデータコピーについて説明するための図である。 第２実施形態の一応用例に係る通信用情報の設定例を示した図である。 第２実施形態に係る関数Ｓｈａｒｅのプログラムコード例を示した図である。 第２実施形態に係る行列ベクトル積の結果代入コード及び関数Ｒｅｄｕｃｅ＿ｓｕｍのプログラムコード例を示した図である。 第２実施形態に係る技術を適用した場合の並列スケーラビリティに関する評価結果を示した図である。 第２実施形態に係る関数ＤｉｓｔＣｏｉｌＰａのプログラムコード例を示した図である。 第２実施形態に係る関数ＤｉｓｔＣｏｉｌＰａの出力に基づく割当自由度の算出例を示した図である。
以下に添付図面を参照しながら、本発明の実施形態について説明する。なお、本明細書及び図面において実質的に同一の機能を有する要素については、同一の符号を付することにより重複説明を省略する場合がある。
＜１．第１実施形態＞
図１を参照しながら、第１実施形態について説明する。図１は、第１実施形態に係る並列処理装置の一例を示した図である。なお、図１に示した並列処理装置１０は、第１実施形態に係る並列処理装置の一例である。
第１実施形態は、係数行列が疎行列になる行列方程式の求解問題を解く方法に関し、それぞれが演算の実行単位となる複数のプロセスを係数行列の部分領域に割り当てて並列に実行することで効率的に問題を処理する並列処理方法を提供する。
係数行列に含まれる非ゼロ成分のパターン（非ゼロパターン）によっては、ある一定の並列数を超えると、並列数を増加させても処理速度が向上しない状況になりうる。つまり、並列処理のスケーラビリティが低くなる非ゼロパターンが存在する。第１実施形態は、このような非ゼロパターンを有する係数行列を扱う場合であっても並列処理のスケーラビリティを向上できるようにする技術を提供する。
図１に示すように、並列処理装置１０は、記憶部１１及び演算部１２Ａ、…、１２Ｆを有する。なお、説明の都合上、演算部の数を６個（演算部１２Ａ、…、１２Ｆ）としているが、演算部の数は２以上の任意の数でよい。また、並列処理装置１０は１つのコンピュータ又は１つの筐体に含まれるコンピュータの集合であってもよいし、複数のコンピュータを通信回線で接続した分散処理システムであってもよい。
記憶部１１は、ＲＡＭ（Random Access Memory）などの揮発性記憶装置、或いは、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの不揮発性記憶装置である。演算部１２Ａ、…、１２Ｆは、ＣＰＵ（Central Processing Unit）やＤＳＰ（Digital Signal Processor）などのプロセッサである。また、演算部１２Ａ、…、１２Ｆの一部は、ＧＰＧＰＵ（General-Purpose computing on Graphics Processing Units）などのプロセッサであってもよい。演算部１２Ａ、…、１２Ｆは、記憶部１１又は他のメモリに記憶されたプログラムを実行する。
演算部１２Ａ、…、１２Ｆは、複数のプロセスを並列に実行することができる。ここで言うプロセスは演算処理の実行単位である。例えば、演算部１２Ａ、…、１２Ｆは、プロセスＰ１、…、Ｐ６を並列に実行可能である。図１の例は、係数行列Ａを含む行列方程式の演算に関する処理を並列処理装置１０で実行する様子を示している。係数行列Ａは、非ゼロ成分が少なく、多くのゼロ成分を含む疎行列である（図１（Ｂ）参照）。
図１（Ａ）に例示した有限要素モデルを考える場合、節点番号２１〜３２が付与された節点の接続関係に基づく非ゼロパターンは、図１（Ｂ）の係数行列Ａのうち最下行及び最右列以外の領域に示した非ゼロ成分の配置パターンとなる。
例えば、構造解析の場合、各節点に変位ベクトルや内力ベクトルが割り当てられ、節点毎に未知数が設定される。この設定により、節点を基準に構造物（連続体）を離散化することで、偏微分方程式で与えられる物理量の振る舞いが離散的な行列方程式で記述される。偏微分方程式を解く際に境界条件を設定するように、行列方程式を解く際にも拘束条件（一部の節点に設定した物理量を既知数に設定する条件）が追加される。
上記の拘束条件は、係数行列Ａに非ゼロ成分を多く含む行及び列を与える（図１（Ｂ）の係数行列Ａのうち最下行及び最右列を参照）。そのため、図１（Ｂ）に例示した非ゼロパターンの係数行列Ａを複数の行群に分け、各行群にプロセスを割り当てた場合、担当する非ゼロ成分の数（処理負荷に相当）が他のプロセスに比べて大幅に多くなるプロセスが生じうる。プロセスの数を増やし、係数行列Ａの分割数を増加させた場合、処理負荷の偏りが拡大する。このような状況は構造解析に限らず、電磁場解析などでも生じうる。
上記のような処理負荷の偏りを是正するため、並列処理装置１０は、係数行列Ａの中で非ゼロ成分が多い行を分割する。そのため、記憶部１１には、疎行列の各行に含まれる非ゼロ成分の多さを判定するための閾値ＴＨが格納されている。閾値ＴＨは、例えば、係数行列Ａの行数や、対角成分及びその周辺に存在する非ゼロ成分（バンド領域）がなすバンドの幅などに基づいて予め設定される。
例えば、演算部１２Ａは、疎行列（係数行列Ａ）の中で閾値ＴＨより非ゼロ成分の数が大きい第１の行Ｊを特定する。また、演算部１２Ａは、図１（Ｃ）に示すように、第１の行Ｊを複数の第２の行Ｊ１、Ｊ２に分割して疎行列を拡大する。図１（Ｃ）の例では、最下行及び最右列がそれぞれ１つずつ拡大されている。演算部１２Ａは、拡大した疎行列（係数行列Ａ）を行群Ｇ１、…、Ｇ６に分けて、行群毎にプロセスを割り当てる。
図１（Ｃ）の例では、行群Ｇ１、…、Ｇ６にそれぞれプロセスＰ１、…、Ｐ６が割り当てられている。演算部１２Ａ、…、１２Ｆは、行群Ｇ１、…、Ｇ６に割り当てられたプロセスＰ１、…、Ｐ６を並列に実行する。例えば、演算部１２Ａ、…、１２Ｆは、それぞれプロセスＰ１、…、Ｐ６（行群Ｇ１、…、Ｇ６に対応する演算処理）を実行する。
行列方程式の解法には、例えば、ＩＣＣＧ法などの反復法がある。ＩＣＣＧ法を適用する場合であれば、演算部１２Ａ、…、１２Ｆは、それぞれ行群Ｇ１、…、Ｇ６に対応する行列ベクトル積（行列とベクトルとの積）などの演算処理を実行する。
上記のように、非ゼロ成分が多い行を分割し、非ゼロ成分の数に応じて係数行列Ａを複数の行群に分けることで、各行群に割り当てられるプロセスの処理負荷をほぼ均等にすることが可能になる。その結果、一部のプロセスが実行する処理が相対的に遅いことで、他のプロセスの処理が滞るリスクが低減される。その結果、処理負荷の偏りに起因するスケーラビリティの制限を抑圧することが可能になる。
以上、第１実施形態について説明した。
＜２．第２実施形態＞
次に、第２実施形態について説明する。第２実施形態は、係数行列が疎行列になる行列方程式の求解問題を解く方法に関し、それぞれが演算の実行単位となる複数のプロセスを係数行列の部分領域に割り当てて並列に実行することで効率的に問題を処理する並列処理方法を提供する。この並列処理方法は、係数行列を適切に分割してプロセス間における処理負荷の偏りを抑制し、並列処理のスケーラビリティを向上させる技術を提供する。
［２−１．ハードウェア］
図２及び図３を参照しながら、第２実施形態に係る並列処理方法を実現することが可能な情報処理装置のハードウェアについて説明する。なお、図２に示した情報処理装置１００及び図３に示した情報処理装置１００ａ、…、１００ｆは、第２実施形態に係る並列処理方法を実現することが可能な情報処理装置の一例である。
図２は、第２実施形態に係るハードウェア（単一装置）の一例を示した図である。
情報処理装置１００は、ＣＰＵ群１０１、メモリ１０２、通信インターフェース１０３、表示インターフェース１０４、及び機器インターフェース１０５を有する。ＣＰＵ群１０１は、ＣＰＵ１０１ａ、１０１ｂ、…、１０１ｆを含む。ＣＰＵ１０１ａ、１０１ｂ、…、１０１ｆ、メモリ１０２、通信インターフェース１０３、表示インターフェース１０４、及び機器インターフェース１０５はバス１０６を解して接続されている。ＣＰＵ群１０１に含まれるＣＰＵの数は２以上の任意の数でよい。
ＣＰＵ１０１ａ、１０１ｂ、…、１０１ｆは、例えば、演算装置又は制御装置として機能し、メモリ１０２に記録された各種プログラムに基づいてハードウェア要素の動作全般又はその一部を制御する。ＣＰＵ１０１ａ、１０１ｂ、…、１０１ｆは、それぞれ複数のプロセッサコアを有していてもよい。なお、ＣＰＵ群１０１には、ＧＰＧＰＵが含まれていてもよい。
メモリ１０２は、ＣＰＵ１０１ａ、１０１ｂ、…、１０１ｆに読み込まれるプログラムや演算に用いられるデータ、そのプログラムを実行する際に変化する各種パラメータなどを一時的又は永続的に保持する記憶装置の一例である。メモリ１０２は、例えば、ＲＡＭなどの揮発性記憶装置、或いは、ＨＤＤやフラッシュメモリなどの不揮発性記憶装置である。
通信インターフェース１０３は、ネットワーク２０１に接続するための通信デバイスである。通信インターフェース１０３は、例えば、有線又は無線ＬＡＮ（Local Area Network）用の通信回路や光通信用の通信回路である。ネットワーク２０１は、有線又は無線により接続されたネットワークであり、例えば、インターネットやＬＡＮなどである。
表示インターフェース１０４は、表示装置２０２に接続するための接続デバイスである。表示装置２０２は、例えば、ＣＲＴ（Cathode Ray Tube）、ＬＣＤ（Liquid Crystal Display）、ＰＤＰ（Plasma Display Panel）、又はＥＬＤ（Electro-Luminescence Display）などである。
機器インターフェース１０５は、入力装置２０３などの外部機器を接続するための接続デバイスである。機器インターフェース１０５は、例えば、ＵＳＢ（Universal Serial Bus）ポート、ＩＥＥＥ１３９４ポート、ＳＣＳＩ（Small Computer System Interface）、ＲＳ−２３２Ｃポートなどである。機器インターフェース１０５には、着脱可能な記録媒体であるリムーバブル記録媒体（非図示）やプリンタなどの外部機器が接続されうる。リムーバブル記録媒体としては、例えば、磁気ディスク、光ディスク、光磁気ディスク、又は半導体メモリなどがある。
図２のように複数のＣＰＵが搭載された単一の情報処理装置１００を利用して第２実施形態に係る演算方法を実現することもできるが、図３に示すような分散処理システムに第２実施形態の技術を適用することもできる。図３は、第２実施形態に係るハードウェア（複数装置）の一例を示した図である。図３の例では、ＣＰＵ１０１ａ、１０１ｂ、…、１０１ｆがそれぞれ情報処理装置１００ａ、１００ｂ、…、１００ｆに搭載されている。そして、情報処理装置１００ａ、１００ｂ、…、１００ｆは通信回線で接続されている。
メモリ１０２ａ、１０２ｂ、…、１０２ｆは、上述したメモリ１０２と同じである。通信インターフェース１０３ａ、１０３ｂ、…、１０３ｆは、上述した通信インターフェース１０３と同じである。説明の都合上、情報処理装置１００ａ、１００ｂ、…、１００ｆでは表示インターフェース１０４及び機器インターフェース１０５に相当する要素の記載が省略されている。情報処理装置１００ａ、１００ｂ、…、１００ｆは、相互に演算結果を送受信しながらＣＰＵ１０１ａ、１０１ｂ、…、１０１ｆに分散した演算処理を実行することができる分散処理システムとして動作しうる。
図２に示した情報処理装置１００又は図３に示した分散処理システムのハードウェアを利用することで、第２実施形態に係る演算方法を実現することができる。なお、図２及び図３に示したハードウェアは一例である。例えば、複数のＧＰＧＰＵや複数のＣＰＵコアを並列に動作させて複数のプロセスを並列実行するシステムに第２実施形態に係る技術を適用することもできる。以下では、説明の都合上、図２に示した情報処理装置１００を利用する例について説明する。
以上、ハードウェアについて説明した。
［２−２．機能］
次に、図４を参照しながら、情報処理装置１００の機能について説明する。図４は、第２実施形態に係る情報処理装置が有する機能の一例を示したブロック図である。
図４に示すように、情報処理装置１００は、記憶部１１１、行列拡張部１１２、プロセス割当部１１３、及び並列計算部１１４を有する。
記憶部１１１の機能は、例えば、メモリ１０２を利用して実現できる。行列拡張部１１２及びプロセス割当部１１３の機能は、例えば、ＣＰＵ１０１ａの機能を利用して実現できる。並列計算部１１４の機能は、例えば、ＣＰＵ群１０１に含まれる複数のＣＰＵ（例えば、ＣＰＵ１０１ａ、…、１０１ｆの一部又は全部）を利用して実現できる。
記憶部１１１には、係数行列１１１ａ及び閾値１１１ｂが格納される。係数行列１１１ａは、演算処理の対象となる行列方程式に含まれる係数行列である。例えば、ベクトルｂ、ｘ、及び係数行列Ａを含む行列方程式「Ａｘ＝ｂ」を解く場合、係数行列Ａの情報が記憶部１１１に格納される。この場合、記憶部１１１には、ベクトルｂ（右辺ベクトル）の情報も格納される。
閾値１１１ｂは、係数行列１１１ａの各行における非ゼロ数（非ゼロ成分の数）の多さを判断するための閾値である。この閾値は、例えば、係数行列Ａのサイズや情報処理装置１００が並列に実行するプロセス（演算処理の実行単位）の数などに応じて設定される。
行列拡張部１１２は、係数行列Ａの各行に含まれる非ゼロ成分の数（非ゼロ数）を計数し、非ゼロ数が上記の閾値より大きい行を特定する。そして、行列拡張部１１２は、特定した行をそれぞれ複数の行に分割して係数行列Ａを拡張する。つまり、行列拡張部１１２は、非ゼロ数が多い行を非ゼロ数が少ない複数の行に分割する。
プロセス割当部１１３は、情報処理装置１００が並列に実行するプロセスの数に応じて、行列拡張部１１２が拡張した係数行列Ａを複数の行群に分割する。なお、行群とは、１つ又は複数の行を含む行の集合である。そして、プロセス割当部１１３は、各行群にプロセスを割り当てる。並列計算部１１４は、各行群に割り当てられたプロセスをＣＰＵ群１０１に含まれる複数のＣＰＵ（ＣＰＵ１０１ａ、…、１０１ｆの一部又は全部）に実行させる。
（適用例）
ここで、図５〜図９を参照しながら、ＩＣＣＧ法に基づく行列方程式の計算アルゴリズムに第２実施形態の技術を適用する方法について説明する。なお、ＩＣＣＧ法は一例であり、第２実施形態に係る技術は、ＩＣＣＧ法以外の反復法に基づく行列方程式の求解問題にも適用可能である。
図５は、行列方程式の構造及びＩＣＣＧ法の擬似コードを示した図である。図６は、行列方程式の構造及び並列ＩＣＣＧ法の擬似コードを示した図である。図７は、第２実施形態に係る行列拡張方法の一例を示した図である。図８は、第２実施形態に係る領域分割方法の一例を示した図である。図９は、第２実施形態に係る並列ＩＣＣＧ法の擬似コードを示した図である。
まず、ＩＣＣＧ法とは、ＣＧ法に不完全コレスキー（ＩＣ：Incomplete Cholesky）分解と呼ばれる前処理を組み合わせた演算手法である。例えば、行列Ａが対象行列の場合、下記の式（１）に示すように、対角行列Ｄ及び非対角行列Ｌを用いて行列Ａを分解することができる。このとき、対角行列Ｄ及び非対角行列Ｌは下記の式（２）及び式（３）で与えられる。
Figure 0006601222
ＩＣ分解は、「非対角行列Ｌに含まれる非ゼロ成分の配置パターン（非ゼロパターン）が行列Ａの非ゼロパターンと同じである」という制約の下で行列Ｃ（Ｃ＝ＬＤＬＴ≒Ａ）を計算する手法である。なお、上付きのＴは転置を表す。このような制限を設ける利点は、非対角行列Ｌの値を格納する配列のサイズを予め確定できる点にある。また、ＩＣＣＧ法は、係数行列Ａが図５（Ａ）に示すような疎行列の場合に演算負荷や使用メモリ量を低減できる利点がある。
ベクトルｂ、ｘ、及び係数行列Ａを含む行列方程式「Ａｘ＝ｂ」を解くＩＣＣＧ法の疑似コードは、例えば、図５（Ｂ）のようになる。なお、Ａ、Ｃは行列である。ｋなどのインデックスが付されたｂ、ｒ、ｘ、ｚ、ｐ、ｑはベクトルである。ｋなどのインデックスが付されたρ、ε、α、βはスカラーである。ｋは整数である。（・，・）はベクトルの内積を表す。Ｓｑｒｔ（・）は平方根を表す。｜・｜はベクトルの大きさを表すノルムである。なお、疑似コードの左端に表示された０１などの数字は行番号である。
図５（Ｂ）の疑似コードを実行することでベクトルｘの解（近似解）を得ることができる。この疑似コードの中で、７行目にある「ｑk＝Ａｐk」の演算は、右辺に行列Ａとベクトルｐkとの積（行列ベクトル積）を含む。このような行列ベクトル積の演算は演算負荷が高い。そのため、行列ベクトル積の演算をいかに効率良く実行できるかが、処理負荷の低減及び処理時間の短縮に影響する。
そこで、図６（Ａ）に示すように、係数行列Ａを複数の部分領域（例えば、｛Ａ（１）｝、｛Ａ（２）｝）に分割し、各部分領域の処理を異なるＣＰＵに割り当てて並列に処理を実行する方法（並列ＩＣＣＧ法）が用いられる。なお、ＩＣ分解を適用する場合、図６（Ａ）の鎖線で囲まれたブロック領域（｛Ａ（ＩＣ１）｝、｛Ａ（ＩＣ２）｝）にある非ゼロ成分が処理に利用される。
並列ＩＣＣＧ法の疑似コードは図６（Ｂ）のようになる。この疑似コードの７行目には、関数Ｓｈａｒｅ（・）の処理が含まれている。関数Ｓｈａｒｅ（・）は、ある部分領域の処理を担当するＣＰＵが、自己の演算結果の一部を他のＣＰＵに送り、他のＣＰＵの演算結果のうち、自己の演算に用いる演算結果を受け取る処理を表す。
つまり、関数Ｓｈａｒｅ（・）は、並列に処理を実行する複数のＣＰＵが、それぞれ処理に利用するデータ（演算結果）を共有する処理を表現したものである。例えば、図６（Ｂ）の疑似コードには、９行目、１３行目、１６行目に内積の演算を含む。これらの内積の演算には、他のＣＰＵが計算したデータが用いられる。そのため、ＭＰＩ（Message Passing Interface）関数を利用してデータをやり取りする処理が挿入される。なお、関数Ｓｈａｒｅ（・）のプログラムコード例は後に示す。
図６（Ａ）に例示した係数行列Ａは、対角成分及びその周辺に非ゼロ成分を有するバンドマトリクスといえる。この場合、係数行列Ａを同程度の非ゼロ成分を含む複数の行群（図６の例では２つの行群）に分けることで負荷の分散を図ることができる。しかし、図７に示すように、他の行と比べて非ゼロ数が大幅に大きい行（図７（Ａ）の例では最下行）が係数行列Ａに含まれる場合、図６（Ｂ）の疑似コードをそのまま適用すると、負荷が十分に分散されないことがある。
図７（Ａ）に示した係数行列Ａの右側に（３）、（３）、…、（１１）という数字が並んでいるが、これらの数字は、対応する行の非ゼロ数である。例えば、係数行列Ａを３個の行群に分ける場合、各行群に含まれる非ゼロ成分の数をほぼ同程度に調整することは可能である。そのため、３つのＣＰＵに同程度の処理を分担させることが可能になり、並列数（この場合は３）に見合った処理性能の向上が期待できる。
一方、全行数が１３行の係数行列Ａを１０個の行群に分ける場合、行群の集合は、非ゼロ成分が１１（又は１１以上）の行群と、非ゼロ数が５程度の行群とを含むことになる。既に述べたように、ある行群を担当するＣＰＵの演算に、他の行群を担当するＣＰＵの演算結果が利用される。そのため、非ゼロ数が大きい行群を担当するＣＰＵの演算結果を、非ゼロ数が小さい行群を担当するＣＰＵが待つことになる。つまり、並列数（この場合は１０）に見合った処理性能の向上が得られないリスクが高い。
そのため、第２実施形態では、行列拡張部１１２が、非ゼロ数が閾値Ｓ（例えば、Ｓ＝７）より大きい行（ラージ行）を特定し、図７（Ｂ）に示すように、非ゼロ数が小さい複数の行に分割する。このとき、行列拡張部１１２は、分割後の係数行列Ａが対称行列となるように列（図７（Ｂ）の例では最右列）も分割する。
また、行列拡張部１１２は、ラージ行に対応するベクトル｛ｐ｝｛ｑ｝を拡張する。例えば、ラージ行に対応するベクトル｛ｑ｝の要素がＱである場合、分割後の行に対応する要素をそれぞれＱ1、Ｑ2とする。但し、「Ｑ＝Ｑ1＋Ｑ2」である。ラージ行に対応するベクトル｛ｐ｝の要素がＰの場合、分割後の行に対応する要素の値は非ゼロ成分の配置に依存する。図７（Ｂ）の例では、係数行列Ａの右下端にある非ゼロ成分を分割後の下行右端に割り当てている。この場合、ベクトル｛ｐ｝の最下段にある要素がＰとなり、その上段が０となる。
上記のように係数行列Ａを拡張することで、図８（Ａ）に示すように、非ゼロ数のバランスを考慮して係数行列Ａを複数の行群に分割することが可能になる。図８（Ａ）の例では、係数行列Ａが行群Ｇ１、…、Ｇ６に分割されている。なお、ＩＣＣＧ法を適用する場合、プロセス割当部１１３は、演算に用いる非ゼロ成分の範囲（鎖線で囲った範囲）を考慮し、その範囲に含まれる非ゼロ数を考慮して行群のサイズを決める。そして、プロセス割当部１１３は、各行群を異なるプロセスに割り当てる。図８（Ｂ）の例では、行群Ｇ１、…、Ｇ６が実行プロセスＰa、…、Ｐfに割り当てられている。
図７（Ｂ）に示した係数行列Ａの拡張（行列拡張）に関する処理を考慮した並列ＩＣＣＧ法の疑似コードは図９のようになる。係数行列Ａを拡張した場合、上記のように、ラージ行に対応するベクトル｛ｑ｝の要素Ｑも分割されるため、分割された要素（上記の例ではＱ1、Ｑ2）の和を計算する処理が追加される。図９に示した疑似コードでは、９行目に追加された関数Ｒｅｄｕｃｅ＿ｓｕｍ（・）が、要素の和を計算する処理を担う。
但し、ｃは係数行列Ａに含まれるラージ行の数、Ｎkcは並列数（並列実行するプロセスの数）、Ｑ（ｋ，ｍ）はｋ番目のラージ行に対応するベクトル｛ｑ｝の要素を分割して得たｍ番目の分割要素を表す。例えば、図７（Ｂ）の例ではｃが１であり、Ｑ1はＱ（１，１）に対応し、Ｑ2はＱ（１，２）に対応する。図８の例ではＮkcが６となる。
以上、情報処理装置１００の機能について説明した。
［２−３．処理の流れ］
次に、図１０及び図１１を参照しながら、ベクトルｂ、ｘ及び係数行列Ａを含む行列方程式「Ａｘ＝ｂ」を解く際の情報処理装置１００による処理の流れについて説明する。図１０は、第２実施形態に係る情報処理装置による処理の流れを示した第１の図である。図１１は、第２実施形態に係る情報処理装置による処理の流れを示した第２の図である。
（Ｓ１０１）行列拡張部１１２は、並列数Ｎ、閾値Ｓ、係数行列Ａ、右辺ベクトルｂ、非ゼロ成分の総数ｎallを取得する。並列数Ｎは、並列処理に用いるプロセスの数（１つのＣＰＵに１つのプロセスを割り当てる場合にはＣＰＵの数）である。閾値Ｓは、係数行列Ａのサイズなどに応じて予め設定されている。非ゼロ成分の総数ｎallは、係数行列Ａに含まれる非ゼロ成分の総数である。
（Ｓ１０２、Ｓ１０７）行列拡張部１１２は、パラメータｋを１からｎdまで変化させながら、Ｓ１０３からＳ１０６までの処理を繰り返し実行する。ｎdは、係数行列Ａの全行数である。
（Ｓ１０３）行列拡張部１１２は、係数行列Ａのｋ行目に位置する非ゼロ成分の数ｎkをカウントする。
（Ｓ１０４）行列拡張部１１２は、Ｓ１０３でカウントした非ゼロ成分の数ｎkが閾値Ｓより大きいか否かを判定する。非ゼロ成分の数ｎkが閾値Ｓより大きい場合、処理はＳ１０５へと進む。一方、非ゼロ成分の数ｎkが閾値Ｓより大きくない場合、処理はＳ１０７へと進み、パラメータｋがｎd以下の場合にはＳ１０２以降の処理が実行される。
（Ｓ１０５）行列拡張部１１２は、ｋ行目の行をラージ行に設定する。ラージ行は、非ゼロ成分を多く含む行として分割の対象になる行である。
（Ｓ１０６）行列拡張部１１２は、パラメータｎcにｎkを加算する。パラメータｎcは、係数行列Ａのラージ行に含まれる非ゼロ成分の総数をカウントするためのパラメータである。Ｓ１０６の処理が完了すると、処理はＳ１０７へと進み、パラメータｋがｎd以下の場合にはＳ１０２以降の処理が実行される。
（Ｓ１０８）行列拡張部１１２は、ラージ行の分割数Ｎcを計算する。ラージ行の分割数Ｎcは、係数行列Ａに含まれるラージ行の領域をいくつの行群に分割するかを示すパラメータである。
例えば、図７に示すようにラージ行（図７の例では最下行）が１つの場合、ラージ行の分割数Ｎcは、下記の式（５）で与えられる。ｎaは、ｎallからｎcを引いた数である。関数ＦＬＯＯＲ（・）は、浮動小数を整数に変換する関数であり、関数に渡された値の小数点以下を切り捨てた整数を返す。
下記の式（５）は、下記の式（４）の関係式に基づいている。つまり、ラージ行の非ゼロ成分をＮc個のグループに分けた場合の非ゼロ数と、ラージ行以外の領域に含まれる非ゼロ成分を（Ｎ−Ｎc）個のグループに分けた場合の非ゼロ数とが同等になるようにＮcが決められる。このような方法でＮcを決めることで、係数行列Ａを分割して得られる複数の行群について非ゼロ数の偏りを抑制することができる。
Figure 0006601222
なお、ラージ行が複数の場合にＮcを決める方法については後述する。
（Ｓ１０９）行列拡張部１１２は、係数行列Ａに含まれるラージ行の領域をＮc個の行群に分割する（行列拡張）。例えば、ラージ行が１つの場合、行列拡張部１１２は、図７に示すような方法で係数行列Ａのラージ行をＮc個の行群に分割する（図７の例では、１つのラージ行が、それぞれ１つの行を含む２個の行群に分割されている）。また、行列拡張部１１２は、係数行列Ａが対象行列となるように係数行列Ａの列を拡張すると共に、係数行列Ａの行列ベクトル積に含まれるベクトル｛ｐ｝及び行列ベクトル積を格納するベクトル｛ｑ｝を拡張する（図７を参照）。
（Ｓ１１０）プロセス割当部１１３は、ラージ行を除く係数行列Ａの領域を（Ｎ−Ｎc）個の行群に分割する（領域分割）。例えば、図８の例では、ラージ行を除く係数行列Ａの領域が行群Ｇ１、…、Ｇ４に分割されている。
（Ｓ１１１）プロセス割当部１１３は、Ｓ１０９でラージ行を分割して得た各行にプロセスを割り当てると共に、Ｓ１１０の領域分割で得た各行群にプロセスを割り当てる。
このとき、プロセス割当部１１３は、図１２に示すようなデータ構造を有する行列情報を生成する。図１２は、第２実施形態に係る行列情報のデータ構造を示した図である。図１２に示すように、行列情報には、行のサイズ、非ゼロ成分の数、行の先頭番号、各行の配列先頭番号の配列、列番号の配列、係数の配列を示すパラメータの情報が格納される。
行列情報は、並列実行される各プロセスの担当ＣＰＵに渡される。そのため、プロセス毎に行列情報が生成される。行のサイズは、担当プロセスが処理する行のサイズ（行群をなす行の数）を表す。非ゼロ成分の数は、担当プロセスが処理する行群に含まれる非ゼロ成分の数を表す。行の先頭番号は、担当プロセスが処理する行群の先頭が係数行列Ａの何行目に位置するかを表す。
各行の配列先頭番号の配列は、係数の配列の中で担当プロセスが処理する各行の非ゼロ成分のデータが格納されている場所の先頭位置を表す配列である。列番号の配列は、係数の配列の中で担当プロセスが処理する各行に対応して、非ゼロ成分のデータが格納されている列の位置を表す配列である。係数の配列は、担当プロセスが処理する係数行列Ａの非ゼロ成分が格納される配列である。このような行列情報を各プロセスに配分することで、非ゼロ成分のデータを効率的に受け渡すことができ、メモリ使用量の節約に寄与する。
（Ｓ１１２）プロセス割当部１１３は、プロセスの割り当て内容に応じて、複数のプロセスを並列実行するために用いる通信用情報を生成する。Ｓ１１１でプロセスの割り当て内容が決定すると、どのＣＰＵがどの行群の処理を担当するのかが決まる。そして、各ＣＰＵが演算処理を進めるために、各ＣＰＵがどのＣＰＵから演算結果を取得するか、どのＣＰＵに自己の演算結果を提供するかが特定される。通信用情報は、このようなＣＰＵ間における演算結果の送受信を可能にするための情報である。
例えば、プロセス割当部１１３は、図１３に示すようなデータ構造を有する通信用情報を生成する。図１３は、第２実施形態に係る通信用情報のデータ構造を示した図である。図１３に示すように、通信用情報には、送信時に用いる情報として、送信する値の数、送信先のＣＰＵ数、送信先のＣＰＵ番号、送信する値の配列先頭番号、送信する値の配列番号、送信する値の配列が含まれる。また、通信用情報には、受信時に用いる情報として、受信する値の数、受信先のＣＰＵ数、受信先のＣＰＵ番号、受信する値の配列先頭番号、受信する値の配列番号、受信する値の配列が含まれる。
（Ｓ１１３）並列計算部１１４は、Ｓ１１１で生成された行列情報及びＳ１１２で生成された通信用情報の送受信を実行する。このとき、Ｓ１０１からＳ１１２までの処理を実行したＣＰＵが、ＣＰＵ毎に対応する行列情報及び通信用情報を該当ＣＰＵに送信し、送信された行列情報及び通信用情報を該当ＣＰＵが受信する。つまり、行列情報及び通信用情報が並列処理のプロセスを担当するＣＰＵに配信される。
（Ｓ１１４）並列計算部１１４は、複数のＣＰＵを並列に動作させ、通信用情報に基づくＣＰＵ間の連携を維持しながら、行列情報に基づく各行群の演算処理を実行させる（未知数の計算）。Ｓ１１４の処理が完了すると、図１０及び図１１に示した一連の処理は終了する。
以上、処理の流れについて説明した。
［２−４．応用例：磁場解析］
ここで、第２実施形態の技術を磁場解析に応用する例について説明する。以下では、説明の都合上、並列計算にＣＰＵ＃１、ＣＰＵ＃２、…、ＣＰＵ＃１０を利用可能な場合を例に説明する。
（係数行列）
コイルから発生する磁場を記述する方程式は、ベクトルポテンシャルＡ、電流密度Ｊ、磁気抵抗率ν、鎖交磁束Φ、抵抗Ｒ、電流Ｉ、端子電圧Ｖ、コイルの断面積Ｓ、コイルの電流密度の単位方向ベクトルｎを用いて下記の式（６）〜式（８）で与えられる。式（７）の左辺第１項は磁束の時間変化に起因する誘導起電力を表し、第２項は抵抗による電圧降下を表す。
Figure 0006601222
例えば、図１４（Ａ）に示すように、矩形のコイルに２５個の節点を設定し、斜線部分を電流が未知数の領域、他の部分を電流が既知の領域とする有限要素モデルを考えると、節点の接続関係に基づく係数行列Ａの非ゼロパターンは図１４（Ｂ）のようになる。図１４は、第２実施形態の一応用例に係る係数行列（節点の接続関係）の作成方法を示した図である。なお、番号（１）、（１１）、（２１）は、行及び列の位置を識別しやすくするために記載した行番号及び列番号であり、節点番号に対応する。
電流密度Ｊが固定であれば係数行列Ａは図１４（Ｂ）のようなバンドマトリクスとなる。一方、電流密度Ｊに未知数（電流未知数）を含む場合、係数行列Ａの非ゼロパターンは、図１５のようになる。図１５に示すように、最下行及び最右列に電流未知数に対応する非ゼロ成分が追加される。図１４（Ａ）の有限要素モデルでは節点４から２５の電流値が未知であるため、節点４から２５に対応する成分が非ゼロとなる。図１５は、第２実施形態の一応用例に係る係数行列（電流未知数を追加）の作成方法を示した図である。
（行列拡張・領域分割）
閾値Ｓが１０の場合、電流未知数に対応する行がラージ行になる。並列数Ｎが１０の場合、ラージ行の分割数Ｎcは、上記の式（４）及び式（５）に基づいて下記の式（９）及び式（１０）のように計算される。下記の式（１０）の結果から、ラージ行の分割数Ｎcは２となる。一方、ラージ行以外の行を含む領域の分割数（Ｎ−Ｎc）は８となる。なお、ＲＯＵＮＤ（・）は四捨五入する関数である。ＭＯＤ（・）は剰余を出力する関数である。
Figure 0006601222
ラージ行を２つに分割し、残りの行を８つの行群に分割すると図１６のようになる。図１６は、第２実施形態の一応用例に係る行列拡張方法及び領域分割方法の一例を示した図である。行群は、各行群に含まれる非ゼロ成分の数が同数に近くなるように設定される。なお、図１６の鎖線領域はＩＣ分解における対象となる領域を示している。このようにして設定された各行群にはプロセスが割り当てられ、割当プロセスを実行するＣＰＵが対応付けられる。図１６の例ではＣＰＵ＃１、ＣＰＵ＃２、…、ＣＰＵ＃１０が行群Ｇ１、Ｇ２、…、Ｇ１０に対応付けられている。
（通信用情報・行列情報）
図１６のように係数行列Ａの拡張及びプロセスの割り当てが完了すると、図１７に示すような行列情報を設定することが可能になる。図１７は、第２実施形態の一応用例に係る行列情報の設定例を示した図である。
図１７に示すように、行列情報は、プロセスを実行するＣＰＵ毎に設定される。図１７の例では、行群Ｇ１を担当するＣＰＵ＃１について、行のサイズ（n＿rows）が４、非ゼロ成分の数（n＿nonzero）が１４、行の先頭番号（n＿row0）が１に設定されている。また、各行の配列先頭番号の配列（ptr＿row[]）、列番号の配列（col[]）、係数の配列（A＿mat[]）が設定されている。なお、Ａi#jは係数行列Ａのｉ行ｊ列成分を表す。ＣＰＵ＃２、…、ＣＰＵ＃１０についても図１６の割り当て内容に応じてパラメータが設定される。
図１６の割り当て内容が決まると、各ＣＰＵが自己の演算により正しい値が得られる列ベクトル（ベクトル｛ｐ｝、ベクトル｛ｑ｝）の成分が特定できる。例えば、図１６に示した係数行列Ａ（拡張後）の非ゼロパターンから、図１８に示すような列ベクトルの非ゼロパターン（ハッチング部分）がＣＰＵ毎に特定できる。図１８は、第２実施形態の一応用例に係る列ベクトルの非ゼロパターンを示した図である。
また、図１６の割り当て内容から、各ＣＰＵが自己の演算で正しい値を取得できる成分（濃いハッチング部分）と、正しい値を取得できない成分（薄いハッチング部分）とが特定できる。正しい値を取得できない成分については、図１９に矢印で示したように、他のＣＰＵが演算処理によって取得した正しい値（演算結果）がコピーされる。図１９は、第２実施形態の一応用例に係るＣＰＵ間のデータコピーについて説明するための図である。
図１９のようにＣＰＵ間で送受信される演算結果のデータ、データの送信元及び送信先が特定できると、図２０に示すような通信用情報を設定することができる。図２０は、第２実施形態の一応用例に係る通信用情報の設定例を示した図である。
図２０の例では、ＣＰＵ＃１に提供される通信用情報として、送信する値の数（n＿all＿send）が２、送信先のＣＰＵ数（num＿cpu＿send）が１、受信する値の数（n＿all＿recv）が３、受信先のＣＰＵ数（num＿cpu＿recv）が２に設定されている。また、送信先のＣＰＵ番号（n＿cpu＿send[]）、送信する値の配列先頭番号（ptr＿send[]）、送信する値の配列番号（vec＿num＿send[]）も設定される。
図１９を参照すると、ＣＰＵ＃１は、上から３段目、４段目の値をＣＰＵ＃２に送信する。そのため、送信する値の数が２、送信先のＣＰＵ数が１となる。また、ＣＰＵ＃１は、自己の演算処理に利用する列ベクトルの成分（非ゼロ成分）のうち、正しい値を得られない成分が５段目、６段目、２６段目にある。そのため、ＣＰＵ＃１は、５段目、６段目の値をＣＰＵ＃２から受信し、２６段目の値をＣＰＵ＃１０から受信する。そのため、受信する値の数が３、受信先のＣＰＵ数が２となる。他のＣＰＵについても同様である。
（プログラムコード例＃１：関数Ｓｈａｒｅ）
上記の通信用情報を取得したＣＰＵは、図２１に示すようなプログラムコードを実行して演算結果をやり取りする。図２１は、第２実施形態に係る関数Ｓｈａｒｅのプログラムコード例を示した図である。図２１の例ではＣ言語による表現を用いている。３行目から１１行目までが送信時の処理であり、１４行目から２２行目までが受信時の処理である。
送信時の処理では、４行目から７行目までの処理で送信する値を送信用の配列にコピーし、１０行目のＭＰＩ関数を利用して配列のデータを他のＣＰＵに送信している。受信時の処理では、１７行目のＭＰＩ関数を利用して配列のデータを他のＣＰＵから受信し、１８行目から２０行目までの処理で受信した値を演算処理に用いるベクトルにコピーしている。関数Ｓｈａｒｅ（・）によりデータが送受信される。
（プログラムコード例＃２：Ｒｅｄｕｃｅ＿ｓｕｍ）
上述した行列拡張を行う際、係数行列Ａの拡張に伴って列ベクトルも拡張される。そのとき、ラージ行に対応する行列ベクトル積が格納される列ベクトルの成分（Ｑ）も分割されるため、その成分を元に戻す処理が演算処理の中に組み込まれる。元に戻す処理を実行する関数が関数Ｒｅｄｕｃｅ＿ｓｕｍである。
関数Ｒｅｄｕｃｅ＿ｓｕｍを実行するにあたり、図２２（Ａ）に示すようなプログラムコードを実行してラージ行を分割して得た各行に対応する行列ベクトル積の演算結果を配列（Q[]）に代入する。そして、図２２（Ｂ）に示す関数Ｒｅｄｕｃｅ＿ｓｕｍのプログラムコードを実行して、その配列の合計値を計算する。図２２は、第２実施形態に係る行列ベクトル積の結果代入コード及び関数Ｒｅｄｕｃｅ＿ｓｕｍのプログラムコード例を示した図である。なお、図２２の例ではＣ言語の表現を用いている。
（効果：並列処理のスケーラビリティ）
これまで説明してきた第２実施形態の技術を適用すると、図２３（本手法）のように高い並列スケーラビリティを実現することができる。図２３は、第２実施形態に係る技術を適用した場合の並列スケーラビリティに関する評価結果を示した図である。図２３には比較のために従来手法の評価結果も示しているが、従来手法の場合、ＣＰＵ数が１５０を超えた辺りでＣＰＵ数の増加に対する速度比の向上が停止している。一方、本手法の場合（第２実施形態の技術を適用した場合）、ＣＰＵ数が２５０を超えてもＣＰＵ数の増加に伴う速度比の向上が見られる。このように、第２実施形態の技術を適用することで、並列スケーラビリティを向上させることが可能になる。
（補足：複数コイルの場合）
上記の応用例は、端子電圧が印加される１つのコイルから発生する磁場を解析する手法であった。ここでは説明の都合上、コイルの数を１つとしたが、例えば、芯材に複数のコイルを巻き付けたインダクタモデルなどにも応用することができる。この場合、複数のコイルに対応する複数のラージ行が係数行列Ａに含まれる。そのため、ラージ行の分割数Ｎcの計算方法は、下記の式（１１）〜式（１５）のように拡張される。
コイルの数をｙcoil、ｙ番目のコイルに関する自由度（総未知数）をｎcyとすると、全てのコイルの自由度ｎc#allは下記の式（１１）で与えられる。また、下記の式（１２）からラージ行全体の分割数Ｎcは下記の式（１３）で与えられる。また、ラージ行毎の分割数Ｎcy（ｙ番目のコイルに対応するラージ行の分割数）は、下記の式（１４）で与えられる。
なお、１つのプロセスに割り当てられるコイルの平均自由度＜ｎc＞は、下記の式（１５）で与えられる。ｎcyが＜ｎc＞に満たない場合、Ｎcyが０となる。この場合、複数のコイルに対応するラージ行をグループ化し、ラージ行のグループに対して少なくとも１つのプロセスが割り当てられるようにする。
例えば、プロセス割当部１１３は、ｎcyが小さい順（ｎcy＜ｎc(y+1)）にコイルの番号（ラージ行の順序に相当）を並べ替え、Ｎcyが１に満たないラージ行をグループ化する。このとき、プロセス割当部１１３は、グループ化したラージ行の全体で分割数が１を超えるようにラージ行の組み合わせを設定する。
Figure 0006601222
グループの識別番号をＩＤ＿Ｇ、グループの分割数をＰａ＿Ｇ、グループの自由度をＤｏｆ＿Ｇ、各ラージ行の分割数をＰａとすると、これらの値は図２４に示すプログラムコード（関数ＤｉｓｔＣｏｉｌＰａ（・））で計算することができる。図２４は、第２実施形態に係る関数ＤｉｓｔＣｏｉｌＰａのプログラムコード例を示した図である。ＩＤ＿Ｇ、Ｐａ＿Ｇ、Ｄｏｆ＿Ｇ、Ｐａはそれぞれ配列の形で出力される。
プロセス割当部１１３は、各グループに対する分割数が１であればラージ行を分割しない。また、プロセス割当部１１３は、分割数が１より大きく、グループに属するラージ行の数が１の場合に分割数Ｐａで分割する。また、プロセス割当部１１３は、分割数が１より大きく、グループに属するラージ行の数が１より大きい場合に、最後のラージ行を同じグループに追加したときに分割数が１を超えたと判断する。そして、プロセス割当部１１３は、同じグループに属する最後のラージ行だけを対象に分割数に応じた分割を行う。
４つのＣＰＵ＃１、…、＃４を５つのラージ行（５つのコイルに対応）に割り当てる例を図２５に示す。図２５は、第２実施形態に係る関数ＤｉｓｔＣｏｉｌＰａの出力に基づく割当自由度の算出例を示した図である。図２４に示した関数ＤｉｓｔＣｏｉｌＰａを実行することで、Ｐａ＿Ｇ及びＰａが得られる。図２５の例では、５つのラージ行が２つの大きなグループ（Pa＿G[0], Pa＿G[1]に対応するグループ）に分けられ、それぞれ自由度に応じてＣＰＵが割り当てられている。
なお、自由度を調整するためにＣＰＵ＃１に４番目のコイルに対応する一部の自由度Ｘ4が割り当てられている。この自由度Ｘ4は、下記の式（１６）により計算される。ＭＯＤ（・）は剰余を求める関数である。図２５のように複数のコイルに対応する複数のラージ行をグループ化し、グループ単位でプロセスを割り当てることで、複数のコイルを含むモデルに対しても柔軟に第２実施形態の技術を適用することが可能になる。
Figure 0006601222
以上、第２実施形態について説明した。
１０ 並列処理装置
１１ 記憶部
１２Ａ、…、１２Ｆ 演算部
２１、…、３２ 節点
Ａ 係数行列
Ｇ１、…、Ｇ６ 行群
Ｊ 第１の行
Ｊ１、Ｊ２ 第２の行
Ｐ１、…、Ｐ６ プロセス
ＴＨ 閾値

Claims (5)

  1. 係数行列として疎行列を含む行列方程式の演算に関する処理を実行するコンピュータに、
    記憶部から、前記疎行列の各行に含まれる非ゼロ成分の多さを判定するための閾値を取得し、前記疎行列の中で前記閾値より非ゼロ成分の数が大きい第１の行を特定し、当該第１の行を複数の第２の行に分割して前記疎行列を拡大し、拡大した前記疎行列のうち、前記第２の行の集合を複数の第１の行群に分け、前記第２の行以外の第３の行の集合を複数の第２の行群に分け、前記第１の行群毎および前記第２の行群毎に、処理の実行単位となるプロセスを割り当てる
    処理を実行させる、行列演算プログラム。
  2. 係数行列として疎行列を含む行列方程式の演算に関する処理を実行するコンピュータに、
    記憶部から、前記疎行列の各行に含まれる非ゼロ成分の多さを判定するための閾値を取得し、前記疎行列の中で前記閾値より非ゼロ成分の数が大きい第１の行を特定し、当該第１の行を複数の第２の行に分割して前記疎行列を拡大し、拡大された前記疎行列の１以上の行を含む行群に含まれる非ゼロ成分の数に応じて当該行群に含める行の組み合わせを決定して、拡大した前記疎行列を複数の行群に分けて、処理の実行単位となるプロセスを行群毎に割り当てる
    処理を実行させる、行列演算プログラム。
  3. 係数行列として疎行列を含む行列方程式の演算に関する処理を実行するコンピュータに、
    記憶部から、前記疎行列の各行に含まれる非ゼロ成分の多さを判定するための閾値を取得し、前記疎行列の中で前記閾値より非ゼロ成分の数が大きい第１の行を特定し、前記第１の行に含まれる非ゼロ成分の数に応じて、前記第１の行を分割して得られる第２の行の数を決定し、決定した数の複数の前記第２の行に前記第１の行を分割して前記疎行列を拡大し、拡大した前記疎行列を複数の行群に分けて、処理の実行単位となるプロセスを行群毎に割り当てる
    処理を実行させる、行列演算プログラム。
  4. 係数行列として疎行列を含む行列方程式の演算に関する処理を実行するコンピュータが、
    記憶部から、前記疎行列の各行に含まれる非ゼロ成分の多さを判定するための閾値を取得し、前記疎行列の中で前記閾値より非ゼロ成分の数が大きい第１の行を特定し、当該第１の行を複数の第２の行に分割して前記疎行列を拡大し、拡大した前記疎行列のうち、前記第２の行の集合を複数の第１の行群に分け、前記第２の行以外の第３の行の集合を複数の第２の行群に分け、前記第１の行群毎および前記第２の行群毎に、処理の実行単位となるプロセスを割り当てる
    行列分割方法。
  5. 処理の実行単位となる複数のプロセスを並列に実行可能な複数の演算部を有し、係数行列として疎行列を含む行列方程式の演算に関する処理を実行する並列処理装置であって、
    前記疎行列の各行に含まれる非ゼロ成分の多さを判定するための閾値が格納される記憶部と、
    前記疎行列の中で前記閾値より非ゼロ成分の数が大きい第１の行を特定し、当該第１の行を複数の第２の行に分割して前記疎行列を拡大し、拡大した前記疎行列のうち、前記第２の行の集合を複数の第１の行群に分け、前記第２の行以外の第３の行の集合を複数の第２の行群に分け、前記第１の行群毎および前記第２の行群毎に、前記プロセスを割り当てる前記演算部と
    を有する、並列処理装置。
JP2016000151A 2016-01-04 2016-01-04 行列演算プログラム、行列分割方法、及び並列処理装置 Active JP6601222B2 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
JP2016000151A JP6601222B2 (ja) 2016-01-04 2016-01-04 行列演算プログラム、行列分割方法、及び並列処理装置
US15/372,921 US20170192818A1 (en) 2016-01-04 2016-12-08 Matrix division method and parallel processing apparatus

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2016000151A JP6601222B2 (ja) 2016-01-04 2016-01-04 行列演算プログラム、行列分割方法、及び並列処理装置

Publications (2)

Publication Number Publication Date
JP2017122950A JP2017122950A (ja) 2017-07-13
JP6601222B2 true JP6601222B2 (ja) 2019-11-06

Family

ID=59226397

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2016000151A Active JP6601222B2 (ja) 2016-01-04 2016-01-04 行列演算プログラム、行列分割方法、及び並列処理装置

Country Status (2)

Country Link
US (1) US20170192818A1 (ja)
JP (1) JP6601222B2 (ja)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6645509B2 (ja) * 2015-11-04 2020-02-14 富士通株式会社 構造解析方法、及び構造解析プログラム
US9934287B1 (en) * 2017-07-25 2018-04-03 Capital One Services, Llc Systems and methods for expedited large file processing
CN107704686B (zh) * 2017-10-11 2020-02-07 大连海事大学 快速分解法潮流计算修正方程系数矩阵的矩阵运算方法
CN107658880B (zh) * 2017-11-16 2019-12-03 大连海事大学 基于关联矩阵运算的快速分解法系数矩阵计算方法
CN107834562B (zh) * 2017-11-16 2019-12-03 大连海事大学 基于Matlab矩阵运算的快速分解法系数矩阵计算法
KR102192560B1 (ko) * 2017-12-07 2020-12-17 한국과학기술원 분산 컴퓨팅 행렬 희소화 방법 및 시스템
CN111240744B (zh) * 2020-01-03 2022-03-22 支付宝(杭州)信息技术有限公司 一种提高涉及稀疏矩阵并行计算效率的方法和系统
WO2021152852A1 (ja) * 2020-01-31 2021-08-05 三菱電機株式会社 制御装置、機械学習装置、および、制御方法
CN114817845B (zh) * 2022-05-20 2023-05-30 昆仑芯(北京)科技有限公司 数据处理方法、装置、电子设备及存储介质
CN115408653B (zh) * 2022-11-01 2023-03-21 泰山学院 一种IDRstab算法高可扩展并行处理方法及系统
CN115952385B (zh) * 2023-03-10 2023-05-05 山东省计算中心（国家超级计算济南中心） 用于大规模稀疏方程组求解的并行超节点排序方法及系统

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4750818A (en) * 1985-12-16 1988-06-14 Cochran Gregory M Phase conjugation method
JP3391262B2 (ja) * 1998-05-11 2003-03-31 日本電気株式会社 記号計算システム及び方法、並びに並列回路シミュレーションシステム
US7065545B2 (en) * 2002-05-07 2006-06-20 Quintero-De-La-Garza Raul Gera Computer methods of vector operation for reducing computation time
US8751556B2 (en) * 2010-06-11 2014-06-10 Massachusetts Institute Of Technology Processor for large graph algorithm computations and matrix operations
US9122932B2 (en) * 2012-04-30 2015-09-01 Xerox Corporation Method and system for automatically detecting multi-object anomalies utilizing joint sparse reconstruction model
US9704102B2 (en) * 2013-03-15 2017-07-11 William Marsh Rice University Sparse factor analysis for analysis of user content preferences
JP6083300B2 (ja) * 2013-03-29 2017-02-22 富士通株式会社 プログラム、並列演算方法および情報処理装置
US10209685B2 (en) * 2015-10-29 2019-02-19 Mitsubishi Electric Research Laboratories, Inc. Method and apparatus for preconditioned model predictive control

Also Published As

Publication number Publication date
JP2017122950A (ja) 2017-07-13
US20170192818A1 (en) 2017-07-06

Similar Documents

Publication Publication Date Title
JP6601222B2 (ja) 行列演算プログラム、行列分割方法、及び並列処理装置
CN109919311B (zh) 生成指令序列的方法、执行神经网络运算的方法和装置
JP6900487B2 (ja) ハードウェアにおける平均プーリングの実行
KR101959376B1 (ko) 멀티 코어 최적화된 순환 신경망을 위한 시스템 및 방법
JP6083300B2 (ja) プログラム、並列演算方法および情報処理装置
US11100386B2 (en) Buffer addressing for a convolutional neural network
US10572225B1 (en) Circuit arrangements and methods for performing multiply-and-accumulate operations
US9037798B2 (en) System and method of operating a computing device to perform memoization including transforming input/output parameters to reduce redundancies and efficiently cache data
JP2019148969A (ja) 行列演算装置、行列演算方法および行列演算プログラム
Gmeiner et al. Parallel multigrid on hierarchical hybrid grids: a performance study on current high performance computing clusters
US20170206089A1 (en) Information processing apparatus and computational method
WO2010039325A1 (en) Method for solving reservoir simulation matrix equation using parallel multi-level incomplete factorizations
CN109885407B (zh) 数据处理方法和装置、电子设备、存储介质
Godwin et al. High-performance sparse matrix-vector multiplication on GPUs for structured grid computations
JP2022550730A (ja) 高速なスパースニューラルネットワーク
US11106968B1 (en) Circuit arrangements and methods for traversing input feature maps
JP2022550170A (ja) コンピュータがテンソルデータの計算を実現する方法、装置、媒体及び機器
US20090171636A1 (en) High-speed operation method for coupled equations based on finite element method and boundary element method
JP2016143210A (ja) 磁界シミュレータプログラム、磁界シミュレータ装置および磁界シミュレーション方法
JP5262177B2 (ja) ベクトル積の並列処理方法
US10210136B2 (en) Parallel computer and FFT operation method
JP6065543B2 (ja) ニューラルネットワーク設計方法、フィッティング方法、及びプログラム
CN111712811A (zh) Hd地图的可扩展图形slam
US20190012831A1 (en) Systems, methods and, media for simulating deformations of nonlinear elastic bodies
JP6981223B2 (ja) 疎行列ベクトル積演算装置及び疎行列ベクトル積演算方法

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20180912

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20190415

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20190423

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20190621

RD02 Notification of acceptance of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7422

Effective date: 20190621

RD04 Notification of resignation of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7424

Effective date: 20190621

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20190910

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20190923

R150 Certificate of patent or registration of utility model

Ref document number: 6601222

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150