solvde.cpp

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

void solvde(int itmax, double conv, double slowc, double scalv[],
			int indexv[], int ne, int nb, int m, double y[],
			int nyj, int nyk, double c[], int nci, int ncj, 
			int nck, double s[], int nsi, int nsj)
{
    double ermax[11];
	int kmax[11];
    int k1 = 1;
    int k2 = m;
    int nvars = ne * m;
    int j1 = 1;
    int j2 = nb;
    int j3 = nb + 1;
    int j4 = ne;
    int j5 = j4 + j1;
    int j6 = j4 + j2;
    int j7 = j4 + j3;
    int j8 = j4 + j4;
    int j9 = j8 + j1;
    int ic1 = 1;
    int ic2 = ne - nb;
    int ic3 = ic2 + 1;
    int ic4 = ne;
    int jc1 = 1;
    int k,j,kp,it,jcf = ic3;
    for (it = 1; it<=itmax; it++)
	{
        k = k1;
        difeq(k, k1, k2, j9, ic3, ic4, indexv, ne, s, nsi, nsj, y,
			  nyj, nyk);
        pinvs(ic3, ic4, j5, j9, jc1, k1, c, nci, ncj, nck, s, nsi,
			  nsj);
        for (k = k1 + 1; k<=k2; k++)
		{
            kp = k - 1;
            difeq(k, k1, k2, j9, ic1, ic4, indexv, ne, s, nsi, nsj,
				  y, nyj, nyk);
            red(ic1, ic4, j1, j2, j3, j4, j9, ic3, jc1, jcf, kp, c,
				nci, ncj, nck, s, nsi, nsj);
			pinvs(ic1, ic4, j3, j9, jc1, k, c, nci, ncj, nck, s,
				  nsi, nsj);
        }
        k = k2 + 1;
        difeq(k, k1, k2, j9, ic1, ic2, indexv, ne, s, nsi, nsj, y,
			  nyj, nyk);
        red(ic1, ic2, j5, j6, j7, j8, j9, ic3, jc1, jcf, k2, c,
			nci, ncj, nck, s, nsi, nsj);
        pinvs(ic1, ic2, j7, j9, jcf, k2 + 1, c, nci, ncj, nck, s,
			  nsi, nsj);
        bksub(ne, nb, jcf, k1, k2, c, nci, ncj, nck);
        double erq = 0.0;
		int jv,km;
		double errj,vmax;
        for (j = 1; j<=ne; j++)
		{
            jv = indexv[j];
            ermax[j] = 0.0;
            errj = 0.0;
            kmax[j] = 0;
            vmax = 0.0;
            for (k = k1; k<=k2; k++)
			{
				double temp;
                double vz = fabs(c[((j-1)*ncj)*nck+k]);
				temp=vz-vmax;
                if (vz > vmax)
				{
                    vmax = vz;
                    km = k;
                }
                errj = errj + vz;
            }
            erq = erq + errj / scalv[jv];
            ermax[j] = c[(j-1)*ncj*nck+km] / scalv[jv];
            kmax[j] = km;
        }
        erq = erq / nvars;
		double dum;
        if (erq > slowc)
		{
            dum = erq;
		}
        else
		{
            dum = slowc;
        }
		double fac;
        fac = slowc / dum;
        for (jv = 1; jv<=ne; jv++)
		{
            j = indexv[jv];
            for (k = k1; k<=k2; k++)
			{
                y[(j-1)*nyk+k] = y[(j-1)*nyk+k] - fac *
					               c[((jv-1)*ncj)*nck+k];
            }
        }
		cout<<"  ";
		cout.width(1);
        cout<<it;
		cout.width(12);
		cout<<erq;
		cout.width(14);
        cout<<fac;
        cout<<endl;
        for (j = 1; j<=ne; j++)
		{
            cout<<"   ";
			cout.width(12);
			cout<<kmax[j];
			cout.width(14);
			cout<<ermax[j]<<endl;
        }
        cout<<endl;
        if (erq < conv)
		{
			break;
		}
    }
}