最適化(線形計画法)

/****************************/
/* 線形計画法               */
/*      coded by Y.Suganuma */
/****************************/
#include <stdio.h>

/*************/
/* クラス LP */
/*************/
class LP
{
		int n;   // 制約条件の数
		int m;   // 変数の数
		int m_s;   // スラック変数の数
		int m_a;   // 人為変数の数
		int mm;   // m + m_s + m_a
		int *a;   // 人為変数があるか否か
		int *cp;   // 比較演算子(-1:左辺<右辺, 0:左辺=右辺, 1:左辺>右辺)
		int *row;   // 各行の基底変数の番号
		double *z;   // 目的関数の係数
		double **s;   // 単体表
		double eps;   // 許容誤差
		int err;   // エラーコード (0:正常終了, 1:解無し)
	public:
		/******************************/
		/* コンストラクタ             */
		/*      n1 : 制約条件の数     */
		/*      m1 : 変数の数         */
		/*      z1 : 目的関数の係数   */
		/*      eq_l : 制約条件の左辺 */
		/*      eq_r : 制約条件の右辺 */
		/*      cp1 : 比較演算子      */
		/******************************/
		LP(int n1, int m1, double *z1, double **eq_l, double *eq_r, int *cp1)
		{
			int i1, i2, k;
					// 初期設定
			eps = 1.0e-10;
			err = 0;
			n   = n1;
			m   = m1;
			a   = new int [n];
			cp  = new int [n];
			for (i1 = 0; i1 < n; i1++) {
				a[i1]  = 0;
				cp[i1] = cp1[i1];
			}
			z = new double [m];
			for (i1 = 0; i1 < m; i1++)
				z[i1] = z1[i1];
					// スラック変数と人為変数の数を数える
			m_s = 0;
			m_a = 0;
			for (i1 = 0; i1 < n; i1++) {
				if (cp[i1] == 0) {
					m_a++;
					if (eq_r[i1] < 0.0) {
						eq_r[i1] = -eq_r[i1];
						for (i2 = 0; i2 < m; i2++)
							eq_l[i1][i2] = -eq_l[i1][i2];
					}
				}
				else {
					m_s++;
					if (eq_r[i1] < 0.0) {
						cp[i1]   = -cp[i1];
						eq_r[i1] = -eq_r[i1];
						for (i2 = 0; i2 < m; i2++)
							eq_l[i1][i2] = -eq_l[i1][i2];
					}
					if (cp[i1] > 0)
						m_a++;
				}
			}
					// 単体表の作成
							// 初期設定
			mm  = m + m_s + m_a;
			row = new int [n];
			s   = new double * [n+1];
			for (i1 = 0; i1 <= n; i1++) {
				s[i1] = new double [mm+1];
				if (i1 < n) {
					s[i1][0] = eq_r[i1];
					for (i2 = 0; i2 < m; i2++)
						s[i1][i2+1] = eq_l[i1][i2];
					for (i2 = m+1; i2 <= mm; i2++)
						s[i1][i2] = 0.0;
				}
				else {
					for (i2 = 0; i2 <= mm; i2++)
						s[i1][i2] = 0.0;
				}
			}
							// スラック変数
			k = m + 1;
			for (i1 = 0; i1 < n; i1++) {
				if (cp[i1] != 0) {
					if (cp[i1] < 0) {
						s[i1][k] = 1.0;
						row[i1]  = k - 1;
					}
					else
						s[i1][k] = -1.0;
					k++;
				}
			}
							// 人為変数
			for (i1 = 0; i1 < n; i1++) {
				if (cp[i1] >= 0) {
					s[i1][k] = 1.0;
					row[i1]  = k - 1;
					a[i1]    = 1;
					k++;
				}
			}
							// 目的関数
			if (m_a == 0) {
				for (i1 = 0; i1 < m; i1++)
					s[n][i1+1] = -z[i1];
			}
			else {
				for (i1 = 0; i1 <= m+m_s; i1++) {
					for (i2 = 0; i2 < n; i2++) {
						if (a[i2] > 0)
							s[n][i1] -= s[i2][i1];
					}
				}
			}
		}

		/*******************************/
		/* 最適化                      */
		/*      sw : =0 : 中間結果無し */
		/*           =1 : 中間結果あり */
		/*      return : =0 : 正常終了 */
		/*             : =1 : 解無し   */
		/*******************************/
		int optimize(int sw)
		{
			int i1, i2, k;
					// フェーズ1
			if (sw > 0) {
				if (m_a > 0)
					printf("\nphase 1\n");
				else
					printf("\nphase 2\n");
			}
			opt_run(sw);
					// フェーズ2
			if (err == 0 && m_a > 0) {
							// 目的関数の変更
				mm -= m_a;
				for (i1 = 0; i1 <= mm; i1++)
					s[n][i1] = 0.0;
				for (i1 = 0; i1 < n; i1++) {
					k = row[i1];
					if (k < m)
						s[n][0] += z[k] * s[i1][0];
				}
				for (i1 = 0; i1 < mm; i1++) {
					for (i2 = 0; i2 < n; i2++) {
						k = row[i2];
						if (k < m)
							s[n][i1+1] += z[k] * s[i2][i1+1];
					}
					if (i1 < m)
						s[n][i1+1] -= z[i1];
				}
							// 最適化
				if (sw > 0)
					printf("\nphase 2\n");
				opt_run(sw);
			}

			return err;
		}

		/*******************************/
		/* 最適化(単体表の変形)      */
		/*      sw : =0 : 中間結果無し */
		/*           =1 : 中間結果あり */
		/*******************************/
		void opt_run(int sw)
		{
			int i1, i2, p, q, k;
			double x, min;

			err = -1;
			while (err < 0) {
					// 中間結果
				if (sw > 0) {
					printf("\n");
					output();
				}
					// 列の選択(巡回を防ぐため必ずしも最小値を選択しない,Bland の規則)
				q = -1;
				for (i1 = 1; i1 <= mm && q < 0; i1++) {
					if (s[n][i1] < -eps)
						q = i1 - 1;
				}
					// 終了(最適解)
				if (q < 0)
					err = 0;
					// 行の選択( Bland の規則を採用)
				else {
					p   = -1;
					k   = -1;
					min = 0.0;
					for (i1 = 0; i1 < n; i1++) {
						if (s[i1][q+1] > eps) {
							x = s[i1][0] / s[i1][q+1];
							if (p < 0 || x < min || x == min && row[i1] < k) {
								min = x;
								p   = i1;
								k   = row[i1];
							}
						}
					}
							// 解無し
					if (p < 0)
						err = 1;
							// 変形
					else {
						x      = s[p][q+1];
						row[p] = q;
						for (i1 = 0; i1 <= mm; i1++)
							s[p][i1] /= x;
						for (i1 = 0; i1 <= n; i1++) {
							if (i1 != p) {
								x = s[i1][q+1];
								for (i2 = 0; i2 <= mm; i2++)
									s[i1][i2] -= x * s[p][i2];
							}
						}
					}
				}
			}
		}

		/****************/
		/* 単体表の出力 */
		/****************/
		void output()
		{
			int i1, i2;

			for (i1 = 0; i1 <= n; i1++) {
				if (i1 < n)
					printf("x%d", row[i1]+1);
				else
					printf(" z");
				for (i2 = 0; i2 <= mm; i2++)
					printf(" %f", s[i1][i2]);
				printf("\n");
			}
		}

		/**************/
		/* 結果の出力 */
		/**************/
		void result()
		{
			double x;
			int i1, i2;

			if (err > 0)
				printf("\n解が存在しません\n");
			else {
				printf("\n(");
				for (i1 = 0; i1 < m; i1++) {
					x = 0.0;
					for (i2 = 0; i2 < n; i2++) {
						if (row[i2] == i1) {
							x = s[i2][0];
							break;
						}
					}
					if (i1 == 0)
						printf("%f", x);
					else
						printf(", %f", x);
				}
				printf(") のとき,最大値 %f\n", s[n][0]);
			}
		}
};

/****************/
/* main program */
/****************/
int main()
{
	double **eq_l, *eq_r, *z;
	int i1, i2, n, m, *cp, sw;
	char c[2];
					// 入力
							// 変数の数と式の数
	scanf("%d %d", &m, &n);
	cp   = new int [n];
	z    = new double [m];
	eq_l = new double * [n];
	eq_r = new double [n];
	for (i1 = 0; i1 < n; i1++)
		eq_l[i1] = new double [m];
							// 目的関数の係数
	for (i1 = 0; i1 < m; i1++)
		scanf("%lf", &z[i1]);
							// 制約条件
	for (i1 = 0; i1 < n; i1++) {
		for (i2 = 0; i2 < m; i2++)
			scanf("%lf", &eq_l[i1][i2]);
		scanf("%s", c);
		if (c[0] == '<')
			cp[i1] = -1;
		else if (c[0] == '>')
			cp[i1] = 1;
		else
			cp[i1] = 0;
		scanf("%lf", &eq_r[i1]);
	}
					// 実行
	LP lp = LP(n, m, z, eq_l, eq_r, cp);
	sw = lp.optimize(1);
					// 結果の出力
	lp.result();

	return 0;
}