最適化(線形計画法)

/****************************/
/* 線形計画法               */
/*      coded by Y.Suganuma */
/****************************/
import java.io.*;
import java.util.StringTokenizer;

/*************/
/* クラス LP */
/*************/
class LP
{
	private int n;   // 制約条件の数
	private int m;   // 変数の数
	private int m_s;   // スラック変数の数
	private int m_a;   // 人為変数の数
	private int mm;   // m + m_s + m_a
	private int a[];   // 人為変数があるか否か
	private int cp[];   // 比較演算子(-1:左辺<右辺, 0:左辺=右辺, 1:左辺>右辺)
	private int row[];   // 各行の基底変数の番号
	private double z[];   // 目的関数の係数
	private double s[][];   // 単体表
	private double eps;   // 許容誤差
	private int err;   // エラーコード (0:正常終了, 1:解無し)

	/******************************/
	/* コンストラクタ             */
	/*      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][mm+1];
		for (i1 = 0; i1 <= n; i1++) {
			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)
				System.out.printf("\nphase 1\n");
			else
				System.out.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)
				System.out.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) {
				System.out.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)
				System.out.printf("x%d", row[i1]+1);
			else
				System.out.printf(" z");
			for (i2 = 0; i2 <= mm; i2++)
				System.out.printf(" %f", s[i1][i2]);
			System.out.printf("\n");
		}
	}

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

		if (err > 0)
			System.out.printf("\n解が存在しません\n");
		else {
			System.out.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)
					System.out.printf("%f", x);
				else
					System.out.printf(", %f", x);
			}
			System.out.printf(") のとき,最大値 %f\n", s[n][0]);
		}
	}
}

/****************/
/* main program */
/****************/
public class Test
{
	public static void main (String[] args) throws IOException
	{
		BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
		double eq_l[][], eq_r[], z[];
		int i1, i2, n, m, cp[], sw;
		String c;
		StringTokenizer str;
					// 入力
							// 変数の数と式の数
		str  = new StringTokenizer(in.readLine(), " ");
		m    = Integer.parseInt(str.nextToken());
		n    = Integer.parseInt(str.nextToken());
		cp   = new int [n];
		z    = new double [m];
		eq_l = new double [n][m];
		eq_r = new double [n];
							// 目的関数の係数
		str = new StringTokenizer(in.readLine(), " ");
		for (i1 = 0; i1 < m; i1++)
			z[i1] = Double.parseDouble(str.nextToken());
							// 制約条件
		for (i1 = 0; i1 < n; i1++) {
			str = new StringTokenizer(in.readLine(), " ");
			for (i2 = 0; i2 < m; i2++)
				eq_l[i1][i2] = Double.parseDouble(str.nextToken());
			c = str.nextToken();
			if (c.compareTo("<") == 0)
				cp[i1] = -1;
			else if (c.compareTo(">") == 0)
				cp[i1] = 1;
			else
				cp[i1] = 0;
			eq_r[i1] = Double.parseDouble(str.nextToken());
		}
					// 実行
		LP lp = new LP(n, m, z, eq_l, eq_r, cp);
		sw = lp.optimize(1);
					// 結果の出力
		lp.result();
	}
}