連載講座111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111

行列演算用言語 LAMAX-S(3 )内田智史，本郷茂，八巻直一
1.はじめにLAMAX白S 入門の最後として、数理計画法の分野から、改訂シンプレックス法と非線形方程式の求解の 2 つの例

題を紹介します。その後、 LAMAX-S の現状と今後についてさまざまな観点からふれたいと思います。

LAMAX-S を用いてこのアルゴリズムを記述するには、少し問題があります。それは LAMAX-S では、大文字と小文字を文法上区別しないので、右辺ベクトル (b) と基底行列 (B) をプログラム上の変数としてそのまま表現できないのです。そこで、アルゴリズムとプログラム中のbiE

耐

Bpc

Bπ

NbeM

ccb

18

久一山

N
CCTbbnM

句

改訂シンプレックス法は、その手順が行列で表現されJているため、 LAMAX-S で記述するには最適な例題です。

ANZAN

2. 改訂シンプレックス法

VE‑ｭ
‑bca

各変数の対応を次のようにしました。

最大化z ニ ctx

条件

Ax= b

LAMAX-S によって書かれたプログラムは、ほとんどx~O

アルゴリズムそのままですので、ポイントとなるところ

このアルゴリズムは文献 [3] から引用すると次のようだけを説明することにします。になります。

このプログラムでは、要素が整数型のベクトル Bindex 、区三百Nindex を宣言しています。これは、基底行列および非制約行列 A を、 B(実行可能基底行列)、 N( 非基底行列)に複写する。

基底行列の変数の番号を保存するためのものです。プログラムの先頭部分で、データを各変数に代入しています

利益ベクトル r を、 CB ， CN に複写するつ17‑7'1:1'1が、このデータは参考文献 [3] から引用しました。

それでは、このプログラムに関連した LAMAX-S の文

方程式 7rB = C}ヲを解いて、 Tを求めるの法事項について説明することにしましょう。 LAMAX-S

被約費用ベクトル cJv ==c
Jv‑7rN を計算する内では、行列の部分をA[ i:j ， k:l] と表現し、これは行列もし、 CN::: 0 ならステップ 4 へ

A の I 行から l 行、 k 列から I 列までの部分行列を意味しでNの最大の要素の依を ëjs とする oます。 i:j や k: 工の部分に*と記述するとそれはすべてのMテップ引範囲を意味し、ユ:キと記述するとそれは ri 以降J を意味

方程式 Bãj" ==αJ" を解いて ãjs を計算する。

します。プログラム中で、 A の部分行列をそれぞれ B 、もし、石弘三 O ならステップ 5 へ

N に代入していますが、それは次の意味です。

方程式 Bb ==b を解いて b を求める。

1 mpmp+l

br/ãrj"_ー min{b;/ãij ， lãij ， >O} なる Tを求める o

IAMテップ引B のけ;1]目と N の js?il] 目を入れ換える。

CB のパ悶目と CN の js1岡目の 11è:を入れ換える。
ステップ 1 へ=B- 1 b,XN ==0 が最適解となり終了

区五ヨIB

B料 (-1 )という表現は、 B の逆行列を意味します。しかし、 LAMAX-S では計算量と精度の観点から、逆行列の計算を可能な限り避ける [1] ように最適化しますn

無限解となり終了

うちだ,

図 2-1 行列 A の部分行列のコピー

区7'1 :1' 41
Xn

NIf I I

pi=cB' 柑林(ー 1) の計算では、実際には逆行列を求めず、さとし神奈川大学〒 221 横浜市神奈川区六角橋 3 -27ー 1

行列 B を三角分解して連立一次方程式の解として求めまほんごう

す n また、 solve 文を用いれば、明示的に連立一時方程やまきしげる専修大学

なおかず綿システム計画研究所

204 (40)

式を解くことができます。

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オベレーションズ・リ+ーチ

solve A*x

医亨7ヂτ1bで述立一次方程式 Ax= b を解いて z を求めます。

アルゴリズム '1" の|ステップ 2 |で、b,.ja'-j, = min{丸 /aij， laij， >O}
なる r を求める部分がありますが、これを LAMAX-S
では次のように allmin 関数を用いて処理します。 allminl現数は行列の最小の値を返す関数です。 allmin(A) とす

ると、行列 A の最小の要素の偵を返しますが、これにいろいろな条件を付加することができます。たとえば、allmin(A , +diag叩al) では、行列 A の対角要素の中から最小の要素を探します。それでは、次の命令は何でしょうか。ここで、 A と B は同じ寸法です。

r =allmin( A, +pos , B .
gt.0.0 )k ニ o とし、初期点を ;Co とする。

|ステップ 1

If (Xk) が十分 O に近ければステ、y プ 5 へ。

|ステ、y プ 2I

連立ー一次方程式 g(xkld =‑J(x ，.) を解いて d を求める。

ZIJ

IAT';
l+d を計算し近似解を更新する。
IAT';I::t 4コ
ハ +1 =Xkk= 止 +1 とし、ステッブ 1 へ

IAT';I::t ヨ解 ;Cj を表示して終了
プログラムは、次のようになります。ところで、 f(x)はベクトル関数であり、 g(x) は行列関数です。たとえば、+pos の場合、第 3 引数に記述された論理演算が行列のこのプログラムでは、 f(x) = (xi+x~-25 ， X1 十 X2 ̲7)t

要素単位ごとに演算され、その結果が真となる要素の {，'rとしています。

ii'i:と同じ位置の行列 A の要素が検索対象となります。つ行列(ベクトル)を返す関数では、構造の指定を func

まり、この例では、 A の最小要素を検索しますが、検索tion 文 iこ記述しますが、その寸法には必ず*を記述します内

対象となるのは、 B の要素が 0 よりも大きい位置ーだけとなります。したがって、 γ を求めるには、

r=allmin( bb Y. ajs , +pos , ajs.gt.O.O )とすればよいのです内ここで日は、要素ごとの除算を意味

この例では、ニュートン法の数学的表現の言語上での再現性がきわめて良いことが示されています。行列表現の導入と、行列関数の概念の導入がアルゴリズムの表現には欠かせない要素であることが認識される例といえま
しょう 3

します。

<=>は、行列やベクトル、要素を交換する命令です内
A<=>B で A と B が交換されます。

このプロクラムでは、変数の番号と値を並べて表示するために、

I

Bindex , xB I という行列を合成しています内これによって 1 列円にはベクトル Bindex ， 2 列目には xB からなる行列が合成されます。

4. LAMAX・ S の今後と関連研究
この連載で、 6 つのプログラム(ヤコビ法(岡有値・阿釘ベクトル)、重回帰分析、主成分分析、改訂シンプレックス法、ニュートン法(非線形方程式))の例を紹介してきました。かなり難しい内容の例題でも、J， AAIAX-S で記述すればひじように簡潔に記述できることがご理解いただけたのではな L 、か ι思います。今までにかなりの例
題を LAMAX-S で記述してきましたが、ほとんど 1 枚3. ニュートン法

OR モデルには、非線形方程式に帰着するものがしばの A4 用紙にプログラムが納まりました。また、プログしば現れます。そのとき、求める解を得るための有効なラミングに要する時間ですが、他の言語で行うのに比べ手法として、ニュートン法があります。ここでは、非線て 5 倍から 10 倍程度のスピード‘でプログラミンクーを行形方程式 J(x) = 0 の解を求めるプログラムを紹介します。ただし、 J， x は n 次元ベクトルとします。また、以下の説明において、 g(x) は f の z におけるヤコビ行列とします。このアルゴリズムは次のようになります門

うことができました。また、何よりもアルゴリズムそのままの形でプログラム上に記述できるため、アルゴリズ
ムのイメージが強烈に浮び上がり、教育的効果が高いことがわかりました。実際、 LAMAX-S を用いて学部学生にプログラムを組ませたところかなり好評でした。
このように教育的な立場から、 LAMAX-S の利用は大いに効果があると考えられますが、我々は実務家や OR

1993 年 4 月号

© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

(41)205

研究者の実際の研究活動にも LAMAX-S は十分に活用できると考えています。現在、実際に我々が提供できる

このような技術を結集し、私たちは、実用的な LAMAX­
S 処理系を目指してゆこうと考えています。

LAMAX-S 処理系は、パソコン版だけであり、多くの文さらに、 LAMAX-S の「その後」のプログラミング環法上の制限があります。近いうちに、我々はこの制限を極境についても研究が進められています。1\'IEX で書かれ力除いてさらにワークステーションに移植する予定です。

たアルゴリズムを LAMAX-S のプログラムに変換する処

このような言語の設計・開発では利用者の意見が大き理系や LAMAX-S のソースコードを生成する数学ソフなウエイトを占めますので、このワークステーション版

トウェアプログラミング支援のための CASE ツールなど処理系を用いて、実際に多くの方々に利用していただき、が、プロトタイプではありますが完成しています。次のLAMAX-S の諸機能をさらに高めていこうと考えていま

アルゴリズムは I阜，TEX によって書かれたヤコビ法によるす。現在の LAMAX-S の文法検討事項として、 (1) 集合連立一次方程式の解法ですが、特殊なトランスレータにの導入、 (2) インデックス行列の導入、 (3) 関数の微分機

よって、この弘司X のファイルを LAMAX-S プログラ能の導入、などがあります0

ムに変換し実行させることができます。

・方、スーパコンビュータに対応可能な本格的な LAMAX­

近い将来、現在あるようなアスキー文字をベースにし円処理系もその開発が進められています。この新処理系

たテキスト中心のプログラミング言語は姿を消すかも知

では、多様なライブラリを利用でき、さらに自動チューれません。私たちの将来の研究は、プログラミング言語

ニング機能および数学的特性や数学的記述をベースにしそのものではなく、数学的ソフトウェアのプログラミンて、計算量や記憶量を減らす処理(これを私達は数学的グ環境やその実現技術に関するものになって行くと考え最適化機能と呼んでいます)も含まれています。この新ています。

処理系は、現在の PC98 版よりも生成される FORTRAN
ヤコビ法(弘TEX → LAMAX回 S)

コードの質を大幅に向上させていて、スーパコンビュー
タだけでなくワークステーション・パソコンにも利用可能になります。

医霊ヨ

現在、この新処理系の開発と平行して、自動チューニ

A:行列 [n xnJ
B: ベクトル [nJ
D:対角行列 [n X n]
L:狭義下三角行列 [n X n]

ング機能を実現するために、各種コンビュータに対するチューニング特性測定システムを開発しています。これは、 LAMAX-S が各コンビュータの特性に合わせたい、わU:狭義上三角行列 [n X 時]
X， Ne切X: ベクトル [n]ゆるチューニングされた )FORTRAN コードを生成できるようにするために、さまざまなサンプルテストプログ

ラムを生成し、実際にコンピコータにかけ、チューニングの特性を調べるものです。この測定結果を元に LAMAX­S のチューニング知識ベースを構築し、自動チューニン

π: 整定数=5ε: 実定数=10- 61Step.1
1A , Bを入力する。

何色ep. 2
1 A を L ， D ， U に分割する。

1
Step. 叶初期点の設定，/

FORTRAN プログラ

1EElli

チューニングされたプリプロセッサ、・、

LAMAX‑S

ソースプログラム・・・唱 ALAMAX‑S唱A

チューニング知識ベース

噌i

チューニング特性iIIlJ 定システム

A
噌

X

=

種々のコンピュータ，，Z，I，BEl--11グを実現する予定です。

図 4・1 チューエング知識ペースの利用

206 (42)© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

オベレーションズ・リサーチ

医百方程式を解〈Cc ステップ 0

方程式 DNe却X=B+(L+U) Xを NewXについて解く。

=diJf (Ne切X-X) の絶対値が最大の行列要素の値も(.， diff ざ叫ま、I Step.三ト行く。

各種データの初期化

B = A(* , mp+l:*]
N = A(*.l:mp]
cB= c(mp+l:*]cN= c[l:mp]c ステップ 1

10000 continue
pi = cB'*Bホキ(-1)
cNb (c耳， ‑ pi牟N )'もし mc> maxcnt ならば、恒三子行く。
X=NewX
mc=mc+lユÏ(allmax(cNb ).le.O.O) g
oto99999

js= almaxr(cNb)

IStep. 子行く。

c ステップ 2

solveB 本 ajs = a[* , js]ii(allmax(ajs).le.O.O) goto 99998
solveB ホ bbbvr = allmi且( bb.
Yajs , +pos , ajs.g七 .0.0)「収束せずJ と表示し、プログラムを終了する。
NcwX を表示し、プログラムを終了する。

なお、 PC98 版 LAMAX-S 処理系入手希望の方はシスc ステザプ S

基底の入れ替え

B[*.rJ<=> N[ホ， js]cB[r]<=> cN[js]Bindex[r] <時 Nindex[js]go to 10000

テム計画研究所 (03-5489-0211 )にお問い合わせください(詳細は本連載 (1) ，Vol.38 ， No.2 , pp.88-91 をご覧下さい)。

c ステップ 4

参考文献

99999 continue
write( キ，*)理、 Vo1. 2 ， No.1 , 1992

2. 内田智史、八巻直一、本郷茂、井上美明、唐津豊:

c ステップ 599998 continue

行列の数学的特性を重視した言語 LAMAX-S の設write(川本)

計思想とその処理系について、日本応用数理学会

end

論文誌、

cVol. 2 , No.3 , pp.155-168 , 1992

改訂シンプレックス法

parameter(mp=3 , np=6.1p=3)
real:matrix [mp , np] A
併制約行列
real:vector[mp]bv <*右辺ベクト lv
real:vectorﾚlp]c
仰利益ベクトル
real:matrix~ ， mp] B
俳基底行列
real:matrix [mp , lp] N<*非基底行列real:vecto~~]cB <*real:ve orQp]
cN <本
real:rvector[mp] pi 併単体乗数ベクトル
real:vector[mp]xb 併最適解
real:vector[mp]ajsp]bb
real:vector[mreal:vector[lpJcNb <*披約費用ベクトル
integer:vector[mp] mask
integer:vector[mpJ Bindex <*基底行列のインc

デックス

cc

ニュートン法x=15 , 81'
10 continueiﾏ( absmax( ﾏ(x) ) .
lt. eps )go to 20
solve G(x) ネ d = ‑ﾏ(x)x x+ dk=k+1
goto 10c
20 continue
callmprint(x)stopendc- ー--・回・開園嗣ーーー---・・・圃・ー・圃ー・・・.._----剛園田・・圃・・・・・・圃--real:vector 凶 Ïunction Ïω<*非基底行列のイ〉〆

デックス

k =0c

“

integer:vector[lp] Nindex
integer js , r

real:vector[2] xc
まデータの設定

1x' * x ・ 26 ，x[1] +x [2] 同7

l'retur耳A = 11 , 1 , 2 , 1 , 0 , 0
2 , 0 , 2 , 0, 1, 0

end

2 , 1 , 3 , 0, 0, 1 I

bv 14 , 5 , 7 l'
13.2.4 , 0.0 , 01'cBindex= 14.5, 6 l'Nindex= 11 , 2 , 3 l'
1993 年 4 月号

，解:無限解'parameter(n=2 , eps=l.e・ 3)
real:vectorﾚlJx
rea工 :vector [n] dreal:vector[*] ﾏ
real:matrix(*β] G

3. 今野浩:線形計画法、日科技連、 pp.60-63

c，終了'solve B 本 xb
bv
callmprint(1 Bindex , xb 1)stop

L 二宮市三:科学技術計算への二つの提案、応用数

c- ・"・・・ー・・・ー・・開園・・・・"・・・・・・・・ー・"ー・・・画--...--------・開園

ﾏunctionG(玄}real:vector[2] x

real:matrix[* ，ヰ]c
G = 12*x[1] ,1,
returnend© 日本オペレーションズ・リサーチ学会. 無断複写・複製・転載を禁ず.

2ヰ x[2J

1107(43)2