本プログラムは，多変数最適化のルーチンを集めたものである．主要なプログ ラムを以下に示す． １．基本演算 mat.c ベクトルと行列の演算 int.c 数値積分 eqs.c 区間縮小法による求根 それぞれ，線形計算，積分，求根のアルゴリズムを提供する． ２．線形計画法 lp.c 線形計画法（シンプレックス法） feas.c 線形計画法（二段階法） シンプレックス法は，目的関数と制約式が線形の場合の最適化手法である． 二段階法は，第一段階で実行可能解の発見，第二段階で最適化を行う． ３．非線形計画法 bfgs.c 準ニュートン法 neldermead.c Nelder-Meadのシンプレックス法 multi.c 乗数法（準ニュートン法使用　勾配要） multisimp.c 乗数法（Nelder-Mead法使用　勾配不要） 準ニュートン法(bfgs.c)とNelder-Mead法(neldermead.c)は，制約なし問題の 最適化手法である．前者は，目的関数の勾配を必要とするが，後者は不要であ る．ただし，収束の速さは前者が勝る．乗数法は，制約つき問題を制約なし問 題に変換する技法である．制約として，等式制約，不等式制約の双方を考慮す ることができる．変換により得られた制約なし問題を，準ニュートン法で解く プログラム(multi.c)，Nelder-Mead法で解くプログラム(multisimp.c)が開 発されている．前者は，目的関数と制約式の勾配を必要とするが，後者は必要 としない． それぞれのプログラムには，テストルーチン(*-t*.c)が付いている．簡単な使 用法は，それを見ること． ******************************************************************************** mat.c ベクトル，行列の演算 extern void vectorcopy(n, x, y) x := y extern dbl vectorvector(n, x, y) return x \dot y (inner product) extern void scalarvector(n, y, k, x) y = kx (k:scalar) extern void vectoradd(n, x, y, z) x = y + z extern void vectorsub(n, x, y, z) x = y - z extern void matrixvector(m, n, y, a, x) y = Ax (A:matrix) extern void matrixrowvector(m,n,y,a,i) y = i-th row of A extern void matrixcolumnvector(m,n,y,a,j) y = j-th column of A extern void matrixcopy(m, n, a, b) A = B extern void matrixunit(n, a) A = I (unit matrix) extern void matrixadd(m, n, a, b, c) A = B + C extern void matrixsub(m, n, a, b, c) A = B - C extern void matrixmatrix(m, n, l, a, b, c) A = BC extern void matrixinverse(n, ainv, a, eps) Ainv = a^{-1} extern void matrixrowscalar(m,n,a,i,p) 行列Aのi行をp倍 extern void matrixrowexchange(m,n,a,i0,i1) 行列Aのi0行とi1行を交換 extern void matrixrowadd(m,n,a,i0,p,i1) 行列Aのi0行にi1行のp倍を加算 extern void matrixrowaddtwo(ma,mb,n,a,b,ia,p,ib) (ma,n)行列Aのia行に(mb,n)行列Bのib行のp倍を加算 extern void matrixsearchcolumnmaxabs(m,n,a,j,is,ie,pimax,pmax) 行列のj列で絶対値最大の要素を探す． ただしis行からie行まで(ie行は含まない) extern void vectorfprint(fd, n, x) ベクトル x の出力 extern void vectorfscan(fd, n, x) ベクトル x の入力 extern void matrixfprint(fd, m, n, a) 行列 A の出力 extern void matrixfscan(fd, m, n, a) 行列 A の入力 extern void matrixprintformat(fmt) 出力形式の変更 mat-t0.c テストルーチン ******************************************************************************** int.c 数値積分 extern dbl GaussIntegral(f, a, b) ガウス公式による数値積分 入力: dbl (*f)() 被積分関数 f(x) dbl a, b 積分区間 返値: 積分値 \int_{a}^{b} f(x) dx extern dbl GaussIntegralParams(f, a, b, npar, params) パラメータ付き関数の積分 入力: dbl (*f)() 被積分関数 f(x;p) dbl a, b 積分区間 int npar パラメータ数 dbl *params パラメータ 返値: 積分値 \int_{a}^{b} f(x;p) dx int-t0.c テストルーチン int-t1.c テストルーチン ******************************************************************************** eqs.c 区間縮小法による求根 extern status IntervalReduction(f, a, b, z, timeout, eps) 入力: dbl (*f)() 関数 dbl a, b 初期区間[a,b] int timeout 繰り返し回数 dbl eps 計算精度 出力: dbl *z f(x) = 0 の根 返値: success 収束 failure f(a)とf(b)の符号が同じ timeout 回数オーバー eqs-t0.c テストルーチン ******************************************************************************** lp.c 線形計画法（シンプレックス法） extern status LPsimplex(m, n, a, b, c, d, eps, ox, ov) 線形計画問題の標準形を解く subject to Ax = b x >= 0 maximize cx + d 標準形：行列 A の後正方部は単位行列 ベクトル b の各要素は非負 入力: int m 等式の数 int n 変数の数 double *a (m,n)係数行列 double *b m次定数ベクトル double *c n次定数ベクトル double *d スカラー double eps 計算精度 出力: double *ox 最適解（存在するとき） double *ov 最適値 返値: success : 最適解が存在する failure : 最適値が発散する extern status LPsimplextableau(m, n, a, b, c, d, eps, base) extern status TwoPhasemethodtableau(m,n,a,b,c,d,cadj,dadj,eps,base) extern void LPsimplexsolution(m, n, b, d, base, ox, ov) lp-t0.c テストルーチン lp-t1.c テストルーチン lp-t-0.dat lp-t-1.dat データファイル lp-t-2.dat ******************************************************************************** feas.c 線形計画法（二段階法） extern status LPfeasiblesolution(m, n, a, b, eps, ox) 実行可能解の発見 find x (>=0) s.t. Ax = b 行列 A，ベクトル b は任意 入力: int m 等式の数 int n 変数の数 double *a (m,n)係数行列 double *b m次定数ベクトル double eps 計算精度 出力: double *ox 実行可能解（存在するとき） 返値: success : 実行可能解が存在する failure : 実行可能解が存在しない extern status TwoPhasemethod(m, n, a, b, c, d, eps, ox, ov) 任意の線形計画問題を解く subject to Ax = b x >= 0 maximize cx + d 行列 A，ベクトル b，ベクトルc，スカラーd は任意 入力: int m 等式の数 int n 変数の数 double *a (m,n)係数行列 double *b m次定数ベクトル double *c n次定数ベクトル double *d スカラー double eps 計算精度 出力: double *ox 最適解（存在するとき） double *ov 最適値 返値: success : 最適解が存在する failure : 最適値が発散する error : 実行可能解が存在しない feas-t0.c テストルーチン（実行可能解の発見） feas-t1.c テストルーチン（実行可能解の発見） feas-t1-0.dat feas-t1-1.dat データファイル feas-t1-2.dat feas-t2.c テストルーチン（二段階法） feas-t3.c テストルーチン（二段階法） feas-t3-0.dat feas-t3-1.dat データファイル feas-t3-2.dat ******************************************************************************** bfgs.c 準ニュートン法 extern status QuasiNewtonMethod(n, f, grad, x, fopt, timeout, eps) 制約なし関数最小化 入力: int n 次元数(変数の数) dbl (*f)() 目的関数 dbl *(*grad)() 目的関数の勾配を求める手続き dbl x[] 初期値 int timeout 繰り返し回数 dbl eps 計算精度 出力: dbl x[] 目的関数が最小値をとる変数の値 dbl *fopt 最小値 返値: success 収束 failure 発散 timeout 回数オーバー bfgs-t0.c テストルーチン bfgs-t1.c テストルーチン bfgs-t2.c テストルーチン ******************************************************************************** neldermead.c Nelder-Meadのシンプレックス法 extern status NelderMeadSimplexMethod(n, f, xinit, length, fopt, timeout, eps) 制約なし関数最小化 入力: int n 次元数(変数の数) dbl (*f)() 目的関数 dbl xinit[] 初期値 dbl length 初期シンプレックスのサイズ int timeout 繰り返し回数 dbl eps 計算精度 出力: dbl xinit[] 目的関数が最小値をとる変数の値 dbl *fopt 最小値 返値: success 収束 failure 発散 timeout 回数オーバー neldermead-t0.c テストルーチン neldermead-t1.c テストルーチン ******************************************************************************** multi.c 乗数法(準ニュートン法を使用) extern status MultiplierMethod(n, f, gradf, m, g, gradg, l, h, gradh, x, timeout, eps) 等式制約と不等式制約のもとでの関数最小化 入力: int n 次元数(変数の数) dbl (*f)() 目的関数 dbl *(*gradf)() 目的関数の勾配を求める手続き int m 不等式制約の数 dbl *(*g)() 不等式制約 dbl **(*gradg)() 不等式制約の勾配 int l 等式制約の数 dbl *(*h)() 等式制約 dbl **(*gradh)() 等式制約の勾配 dbl x[] 初期値 int timeout 繰り返し回数 dbl eps 計算精度 出力: dbl x[] 目的関数が最小値をとる変数の値 返値: success 収束 failure 発散 timeout 回数オーバー multi-t0.c テストルーチン multi-t1.c テストルーチン ******************************************************************************** multisimp.c 乗数法(Nelder-Meadのシンプレックス法を使用) extern status MultiplierMethodSimplex(n, f, m, g, l, h, x, length, timeout, eps) 等式制約と不等式制約のもとでの関数最小化 入力: int n 次元数(変数の数) dbl (*f)() 目的関数 int m 不等式制約の数 dbl *(*g)() 不等式制約 int l 等式制約の数 dbl *(*h)() 等式制約 dbl x[] 初期値 dbl length 初期シンプレックスのサイズ int timeout 繰り返し回数 dbl eps 計算精度 出力: dbl x[] 目的関数が最小値をとる変数の値 返値: success 収束 failure 発散 timeout 回数オーバー multisimp-t0.c テストルーチン ********************************************************************************