/* linear.c */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define EPS 1.0e-6 /* 相対精度 */
void simplx(double **a, int m, int n, int m1, int m2, int m3, int *icase,
int izrov[], int iposv[]);
void simp1(double **a, int mm, int ll[], int nll, int iabf, int *kp, double *bmax);
void simp2(double **a, int n, int l2[], int nl2, int *ip, int kp, double *q1);
void simp3(double **a, int i1, int k1, int ip, int kp);
int *ivector(long n);
void free_ivector(int *v);
double **matrix(long nr, long nc);
void free_matrix(double **m);
/*
線形計画法のシンプレックス法。入力パラメータは、a,m,n,mp,np,m1,m2,m3, 出力
パラメータは、a,icase,izrov,iposv。
*/
void simplx(double **a, int m, int n, int m1, int m2, int m3, int *icase,
int izrov[], int iposv[])
{
int i, ip, ir, is, k, kh, kp, m12, nl1, nl2;
int *l1, *l2, *l3;
double q1, bmax;
if(m != (m1 + m2 + m3))
{
fprintf(stderr, "Error : Bad input constraint counts in simplx.\n");
return;
}
l1 = ivector(n + 1);
l2 = ivector(m);
l3 = ivector(m);
nl1 = n; /* 最初変数は全部右辺に */
for(k = 1; k <= n; k++) l1[k] = izrov[k] = k; /* 添字の表を初期化 */
nl2 = m; /* 人工変数は全部右辺に */
for(i = 1; i <= m; i++) /* 添字の表の初期化 */
{
if(a[i + 1][1] < 0.0)
{
fprintf(stderr, "Error : Bad input tableau in simplx.\n");
return;
}
l2[i] = i;
iposv[i] = n + i;
}
for(i = 1; i <= m2; i++) l3[i] = 1; /* 後で使う配列の初期化 */
ir = 0;
/*
フラグ ir の値 0 はフェーズ２（可能な初期解が得られた状態）を意味する。
座標原点が可能解ならフェーズ２に直行する。
*/
if(m2 + m3)
{
ir = 1; /* フラグ ir の値 1 はフェーズ１から始めることを意味する。 */
for(k = 1; k <= (n + 1); k++) /* 補助目的関数を計算 */
{
q1 = 0.0;
for(i = m1 + 1; i <= m; i++) q1 += a[i + 1][k];
a[m + 2][k] = -q1;
}
do
{
/* 補助目的関数の最大の係数を探す */
simp1(a, m + 1, l1, nl1, 0, &kp, &bmax);
if(bmax <= EPS && a[m + 2][1] < -EPS)
{
*icase = -1;
/*
補助目的関数は負のままであり、これ以上改良できない。つまり、可能解は存在しない。
*/
free_ivector(l3);
free_ivector(l2);
free_ivector(l1);
return;
}
else if(bmax <= EPS && a[m + 2][1] <= EPS)
{
m12 = m1 + m2 + 1;
/*
補助目的関数の最大値 0 が得られた。つまり、可能な初期解が見つかった。
"goto one" で人工変数を一掃し、フェーズ２に進む。
*/
if(m12 <= m)
{
for(ip = m12; ip <= m; ip++)
{
if(iposv[ip] == ip + n)
{
simp1(a, ip, l1, nl1, 1, &kp, &bmax);
if(bmax > 0.0) goto one;
}
}
}
ir = 0; /* フェーズ２に来たことを表すフラグをセットする。 */
--m12;
if(m1 + 1 <= m12)
{
for(i = m1 + 1; i <= m12; i++)
{
if(l3[i - m1] == 1)
for(k = 1; k <= n + 1; k++)
a[i + 1][k] = - a[i + 1][k];
}
}
break; /* フェーズ２に進む。 */
}
/* ピボット要素を見つける。 */
simp2(a, n, l2, nl2, &ip, kp, &q1);
if(ip == 0) /* 補助目的関数は上限がなく、可能解は存在しない。 */
{
*icase = -1;
free_ivector(l3);
free_ivector(l2);
free_ivector(l1);
return;
}
/* 左辺・右辺変数を交換し（フェーズ１）、表を更新する。 */
one: simp3(a, m + 1, n, ip, kp);
if(iposv[ip] >= (n + m1 + m2 + 1))
{
for(k = 1; k <= nl1; k++) if(l1[k] == kp) break;
--nl1;
for(is = k; is <= nl1; is++) l1[is] = l1[is + 1];
++a[m + 2][kp + 1];
for(i = 1; i <= m + 2; i++) a[i][kp + 1] = - a[i][kp + 1];
}
else
{
if(iposv[ip] >= n + m1 + 1)
{
kh = iposv[ip] - m1 - n;
if(l3[kh])
{
l3[kh] = 0;
++a[m + 2][kp + 1];
for(i = 1; i <= m + 2; i++) a[i][kp + 1] = - a[i][kp + 1];
}
}
}
is = izrov[kp];
izrov[kp] = iposv[ip];
iposv[ip] = is;
}while(ir); /* まだフェーズ１なら do に戻る。 */
}
/*
フェーズ１の終わり。初期可能解が得られた。フェーズ２に移り、解を最適化する。
*/
for(;;)
{
simp1(a, 0, l1, nl1, 0, &kp, &bmax); /* z の行を調べる。 */
if(bmax <= 0.0) /* 解は見つかった。 */
{
*icase = 0;
free_ivector(l3);
free_ivector(l2);
free_ivector(l1);
return;
}
simp2(a, n, l2, nl2, &ip, kp, &q1); /* ピボット要素を探す。（フェーズ２） */
if(ip == 0) /* 目的関数は上限がない。 */
{
*icase = 1;
free_ivector(l3);
free_ivector(l2);
free_ivector(l1);
return;
}
simp3(a, m, n, ip, kp); /* 変数を左辺・右辺交換する。（フェーズ２） */
is = izrov[kp];
izrov[kp] = iposv[ip];
iposv[ip] = is;
}
}
/*
表 ll[ ] で与えた添字の要素の最大値を求める。その際に絶対値をとれるかどうかは
フラグ iabf で指示する。
*/
void simp1(double **a, int mm, int ll[], int nll, int iabf, int *kp, double *bmax)
{
int k;
double test;
*kp = ll[1];
*bmax = a[mm + 1][*kp + 1];
for(k = 2; k <= nll; k++)
{
if(iabf == 0) test = a[mm + 1][ll[k] + 1] - (*bmax);
else test = fabs(a[mm + 1][ll[k] + 1]) - fabs(*bmax);
if(test > 0.0)
{
*bmax = a[mm + 1][ll[k] + 1];
*kp = ll[k];
}
}
}
/*
退化（degeneracy）を考慮に入れてピボット要素を探す。
*/
void simp2(double **a, int n, int l2[], int nl2, int *ip, int kp, double *q1)
{
int k, ii, i;
double qp, q0, q;
*ip = 0;
for(i = 1; i <= nl2; i++) /* 可能なピボットはあるか？ */
{
if(a[l2[i] + 1][kp + 1] < -EPS)
{
*q1 = -a[l2[i] + 1][1] / a[l2[i] + 1][kp + 1];
*ip = l2[i];
for(i = i + 1; i <= nl2; i++)
{
ii = l2[i];
if(a[ii + 1][kp + 1] < -EPS)
{
q = -a[ii + 1][1] / a[ii + 1][kp + 1];
if(q < *q1)
{
*ip = ii;
*q1 = q;
}
else if(q == *q1) /* 退化している。 */
{
for(k = 1; k <= n; k++)
{
qp = -a[*ip + 1][k + 1] / a[*ip + 1][kp + 1];
q0 = -a[ii + 1][k + 1] / a[ii + 1][kp + 1];
if(q0 != qp) break;
}
if(q0 < qp) *ip = ii;
}
}
}
}
}
}
/*
左辺変数・右辺変数を置き換える行列演算。
*/
void simp3(double **a, int i1, int k1, int ip, int kp)
{
int kk, ii;
double piv;
piv = 1.0 / a[ip + 1][kp + 1];
for(ii = 1; ii <= i1 + 1; ii++)
{
if(ii - 1 != ip)
{
a[ii][kp + 1] *= piv;
for(kk = 1; kk <= k1 + 1; kk++)
if(kk - 1 != kp)
a[ii][kk] -= a[ip + 1][kk] * a[ii][kp + 1];
}
}
for(kk = 1; kk <= k1 + 1; kk++)
if(kk - 1 != kp) a[ip + 1][kk] *= -piv;
a[ip + 1][kp + 1] = piv;
}
/* allocate a int vector with subscript range v[nl..nh] */
int *ivector(long n)
{
int *v;
v = (int *)malloc((size_t)((n + 1) * sizeof(int)));
if(v == NULL) fprintf(stderr, "Error : allocation failure in ivector().\n");
return v;
}
/* free a int vector allocated with ivector() */
void free_ivector(int *v)
{
free((char *)v);
}
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
double **matrix(long nr, long nc)
{
double **m;
long i;
m = (double **)malloc((size_t)((nr + 1) * sizeof(double *)));
if(m == NULL)
fprintf(stderr, "Error : allocation failure 1 in matrix().\n");
m[1] = (double *)malloc((size_t)((nr * nc + 1) * sizeof(double)));
if(m[1] == NULL)
fprintf(stderr, "Error : allocation failure 2 in matrix().\n");
for(i = 2; i <= nr; i++) m[i] = m[i - 1] + nc;
return m;
}
/* free a double matrix allocated by matrix() */
void free_matrix(double **m)
{
free((char *)m[1]);
free((char *)m);
}
int main(void)
{
double **a;
int *izrov, *iposv;
int i, j, m, n, m1, m2, m3, icase;
m = 4;
n = 4;
m1 = 2; /* <= 条件数 */
m2 = 1; /* >= 条件数 */
m3 = 1; /* = 条件数 */
a = matrix(m + 2, n + 1);
izrov = ivector(n);
iposv = ivector(m);
/*
最大化条件式 : z = x1 + x2 + 3 * x3 - x4 / 2
制約条件 :
x1 + 2 * x3 <= 740
2 * x2 - 7 * x4 <= 0
x2 - x3 + 2 * x4 >= 0.5
x1 + x2 + x3 + x4 = 9
*/
a[1][1] = 0.; a[1][2] = 1.; a[1][3] = 1.; a[1][4] = 3.; a[1][5] = -0.5;
a[2][1] = 740.; a[2][2] = -1.; a[2][3] = 0.; a[2][4] = -2.; a[2][5] = 0.;
a[3][1] = 0.; a[3][2] = 0.; a[3][3] = -2.; a[3][4] = 0.; a[3][5] = 7.;
a[4][1] = 0.5; a[4][2] = 0.; a[4][3] = -1.; a[4][4] = 1.; a[4][5] = -2.;
a[5][1] = 9.; a[5][2] = -1.; a[5][3] = -1.; a[5][4] = -1.; a[5][5] = -1.;
printf("[ input ]\n");
printf("*** a ***\n");
for(i = 1; i <= m + 1; i++)
{
for(j = 1; j <= n + 1; j++) printf("%f ", a[i][j]);
printf("\n");
}
simplx(a, m, n, m1, m2, m3, &icase, izrov, iposv);
printf("[ output ]\n");
printf("icase = %d\n", icase);
printf("*** a ***\n");
for(i = 1; i <= m + 1; i++)
{
for(j = 1; j <= n + 1; j++) printf("%f ", a[i][j]);
printf("\n");
}
printf("*** izrov ***\n");
for(j = 1; j <= n; j++) printf("%d ", izrov[j]);
printf("\n");
printf("*** iposv ***\n");
for(i = 1; i <= m; i++) printf("%d ", iposv[i]);
printf("\n\n");
if(icase == 1) printf("目的関数が有界でない。\n");
else if(icase == -1) printf("解は存在しない。\n");
else
{
printf("解 : \n");
for(i = 1; i <= m; i++)
{
printf("\t x%d\t= ", i);
if(iposv[i] > n) printf("0.0\n");
else printf("%f\n", a[iposv[i] + 1][1]);
}
}
free_ivector(iposv);
free_ivector(izrov);
free_matrix(a);
return 1;
}
|
（トップに戻る）
戻る（メインメニュー）
戻る（その他メニュー）
注意事項