d11r10.cpp

使用VC++编写的大量数学算法的源代码

#include "iostream.h"
#include "math.h"
#include "string.h"


void simp1(double a[], int mp, int np, int mm, int ll[], int nll, int iabf,
		   int& kp, double& bmax)
{
	double test;
    kp = ll[1];
    bmax = a[mm*np+kp + 1];
    for (int k = 2; k<=nll; k++)
	{
        if (iabf == 0)
		{
            test = a[mm*np+ll[k] + 1] - bmax;
		}
        else
		{
            test = fabs(a[mm*np+ll[k] + 1]) - fabs(bmax);
        }
        if (test > 0.0)
		{
            bmax = a[mm*np+ll[k] + 1];
            kp = ll[k];
        }
    }
}

void simp2(double a[], int m, int n, int mp, int np, int l2[], int nl2, int& ip, int kp, double& q1)
{
    double q,qp,q0,eps = 0.000001;
    ip = 0;
    int i,ii,k,flag = 0;
    for (i = 1; i<=nl2; i++)
	{
        if (a[l2[i]*np+kp + 1] < -eps)
		{
			flag = 1;
		}
        if (flag == 1)
		{
			break;
		}
    }
    if (flag == 0)
	{
		return;
	}
    q1 = -a[l2[i]*np+ 1] / a[l2[i]*np + kp + 1];
    ip = l2[i];
    for (i = i + 1; i<=nl2; i++)
	{
        ii = l2[i];
        if (a[ii*np+kp + 1] < -eps)
		{
            q = -a[ii*np+1] / a[ii*np+kp+1];
            if (q < q1)
			{
                ip = ii;
                q1 = q;
			}
            else
			{
				if (q == q1)
				{
					for (k = 1; k<=n; k++)
					{
						qp = -a[ip*np+k + 1] / a[ip*np+kp + 1];
						q0 = -a[ii*np+k + 1] / a[ii*np+kp + 1];
						if (q0 != qp)
						{
							break;
						}
					}
					if (q0 < qp)
					{
						ip = ii;
					}
				}
            }
        }
    }
}

void simp3(double a[], int mp, int np, int i1, int k1, int& ip, int& kp)
{
	int kk;
    double piv = 1.0 / a[ip*np+kp + 1];
    for (int ii = 1; ii<=i1 + 1; ii++)
	{
        if (ii - 1 != ip)
		{
            a[(ii-1)*np+kp + 1] = a[(ii-1)*np+kp + 1] * piv;
            for (kk = 1; kk<=k1 + 1; kk++)
			{
                if (kk - 1 != kp)
				{
                    a[(ii-1)*np+kk] = a[(ii-1)*np+kk] - a[ip*np+kk] * a[(ii-1)*np+kp + 1];
                }
            }
        }
    }
    for (kk = 1; kk<=k1 + 1; kk++)
	{
        if (kk - 1 != kp)
		{
			a[ip*np+kk] = -a[ip*np+kk] * piv;
		}
    }
    a[ip*np + kp + 1] = piv;
}

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;
			//erase l3, l2, l1
			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;
			//erase l3, l2, l1
			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;
        //erase l3, l2, l1
        return;
    }
    simp2(a, m, n, mp, np, l2, nl2, ip, kp, q1);
    if (ip == 0)
	{
        icase = 1;
        //erase l3, l2, l1
        return;
    }
    simp3(a, mp, np, m, n, ip, kp);
    goto r2;
}

void main()
{
    //program d9r10
    //driver for routine simplx
    //incorporates  examples discussed in text
    int icase,i,j,jj,jmax,nm1m2,n = 4;
	int m = 4;
    int np = 5;
    int mp = 6;
    int m1 = 2;
    int m2 = 1;
    int m3 = 1;
    nm1m2 = n + m1 + m2;
    double a[7][6],anum[6],a1[31];
	int izrov[5], iposv[5];
	char* txt[8];
	char* alpha[6];

    txt[1] = "x1"; txt[2] = "x2"; txt[3] = "x3"; txt[4] = "x4";
    txt[5] = "y1"; txt[6] = "y2"; txt[7] = "y3";
    a[1][1]=0.0;  a[1][2]=1.0; a[1][3]=1.0; a[1][4]=3.0; a[1][5]=-0.5;
    a[2][1]=740.0;a[2][2]=-1.0;a[2][3]=0.0; a[2][4]=-2.0;a[2][5]=0.0;
    a[3][1]=0.0;  a[3][2]=0.0; a[3][3]=-2.0;a[3][4]=0.0; a[3][5]=7.0;
    a[4][1]=0.5;  a[4][2]=0.0; a[4][3]=-1;  a[4][4]=1.0; a[4][5]=-2.0;
    a[5][1]=9.0;  a[5][2]=-1.0;a[5][3]=-1.0;a[5][4]=-1.0;a[5][5]=-1.0;
    a[6][1]=0.0;  a[6][2]=0.0; a[6][3]=0.0; a[6][4]=0.0; a[6][5]=0.0;
	for (i=1; i<=mp; i++)
	{
		for (j=1; j<=np; j++)
		{
			a1[(i-1)*np+j]=a[i][j];
		}
	}
    simplx(a1, m, n, mp, np, m1, m2, m3, icase, izrov, iposv);
    if (icase == 1)
	{
        cout<<"unbounded objective function"<<endl;
	}
    else
	{
		if (icase == -1)
		{
			cout<<"no solutions satisfy constraints given"<<endl;
		}
		else
		{
			jj = 1;
			for (i = 1; i<=n; i++)
			{
				if (izrov[i] <= n + m1 + m2)
				{
					alpha[jj]=txt[izrov[i]];
					jj = jj + 1;
				}
			}
			jmax = jj - 1;
			cout.setf(ios::fixed|ios::left);
			cout.precision(3);
			cout<<endl;
			cout<<"                       ";
			for (jj = 1; jj<=jmax; jj++)
			{
				cout<<alpha[jj];
				cout<<"          ";
			}
			cout<<endl;
			for (i = 1; i<=m + 1; i++)
			{
				if (i > 1)
				{
					alpha[1]=txt[iposv[i - 1]];
				}
				else
				{
					alpha[1]="  ";
				}
				anum[1] = a1[(i-1)*np+1];
				jj = 2;
				for (j = 2; j<=n + 1; j++)
				{
					if (izrov[j - 1] <= (n + m1 + m2))
					{
						anum[jj] = a1[(i-1)*np+j];
						jj = jj + 1;
					}
				}
				jmax = jj - 1;
				cout<<alpha[1]<<"      ";
				for (jj = 1; jj<=jmax; jj++)
				{
					cout.width(12);
					cout<<anum[jj];
				}
				cout<<endl;
			}
		}
	}
}