simplx.cpp

Visual C++ 常用数值算法集 源代码

void simplx(double a[], int m, int n, int mp, int np, int m1,
			int m2,	int m3, int& icase, int izrov[], int iposv[])
{
	int kp,ip,m12,iq,kh;
    double q1,bmax,eps = 0.000001;
    int l1[101], l2[101], l3[101];
    if (m != m1 + m2 + m3)
	{
		cout<<" bad input constraint counts"<<endl;
        return;
    }
    int i,k,nl1 = n;
    for (k = 1; k<=n; k++)
	{
        l1[k] = k;
        izrov[k] = k;
    }
    int nl2 = m;
    for (i = 1; i<=m; i++)
	{
        if (a[i*np+1] < 0.0)
		{
            cout<<" bad input tableau."<<endl;
            return;
        }
        l2[i] = i;
        iposv[i] = n + i;
    }
    for (i = 1; i<=m2; i++)
	{
        l3[i] = 1;
    }
    int ir = 0;
    if (m2 + m3 == 0)
	{
		goto r3;
	}
	ir = 1;
	for (k = 1; k<=n + 1; k++)
	{
		double q1 = 0.0;
		for (i = m1 + 1; i<=m; i++)
		{
			q1 = q1 + a[i*np+k];
		}
		a[(m + 1)*np+k] = -q1;
	}
	do
	{
		simp1(a, mp, np, m + 1, l1, nl1, 0, kp, bmax);
		if (bmax <= eps && a[(m + 1)*np+1] < -eps)
		{
			icase = -1;
			return;
		}
		else
		{
			if (bmax <= eps && a[(m + 1)*np+1] <= eps)
			{
				m12 = m1 + m2 + 1;
				if (m12 <= m)
				{
					for (ip = m12; ip<=m; ip++)
					{
						if (iposv[ip] == ip + n)
						{
							simp1(a,mp,np,ip,l1,nl1,1,kp,bmax);
						}
						if (bmax > 0.0)
						{
							goto r1;
						}
					}
				}
				ir = 0;
				m12 = m12 - 1;
				if (m1 + 1 > m12)
				{
					break;
				}
				for (i = m1 + 1; i<=m12; i++)
				{
					if (l3[i - m1] == 1)
					{
						for (k = 1; k=n + 1; k++)
						{
							a[i*np+k] = -a[i*np+k];
						}
					}
				}
				break;
			}
		}
		simp2(a, m, n, mp, np, l2, nl2, ip, kp, q1);
		if (ip == 0)
		{
			icase = -1;
			return;
		}
r1:         simp3(a, mp, np, m + 1, n, ip, kp);
		if (iposv[ip] >= n + m1 + m2 + 1)
		{
			for (k = 1; k<=nl1; k++)
			{
				if (l1[k] == kp)
				{
					break;
				}
			}
			nl1 = nl1 - 1;
			for (iq = k; iq<=nl1; iq++)
			{
				l1[iq] = l1[iq + 1];
			}
		}
		else
		{
			if (iposv[ip] < n + m1 + 1)
			{
				goto r2;
			}
			kh = iposv[ip] - m1 - n;
			if (l3[kh] == 0)
			{
				goto r2;
			}
			l3[kh] = 0;
		}
		a[(m + 1)*np+kp + 1] = a[(m + 1)*np+kp + 1] + 1.0;
		for (i = 1; i<=m + 2; i++)
		{
			a[(i-1)*np+kp + 1] = -a[(i-1)*np+kp + 1];
		}
r2:     iq = izrov[kp];
		izrov[kp] = iposv[ip];
		iposv[ip] = iq;
	}while (ir != 0);
r3: simp1(a, mp, np, 0, l1, nl1, 0, kp, bmax);
    if (bmax <= 0.0)
	{
        icase = 0;
        return;
    }
    simp2(a, m, n, mp, np, l2, nl2, ip, kp, q1);
    if (ip == 0)
	{
        icase = 1;
        return;
    }
    simp3(a, mp, np, m, n, ip, kp);
    goto r2;
}