sparse.cpp

经典算法 主要用于程序员开发 和进行开始学习用

void sparse(double b[], int n, double x[], double rsq)
{
    const int nmax = 500;
    double eps2,eps = 0.000001;
    double g[500], h[500], xi[500], xj[500];
	int j;
    eps2 = n * pow(eps, 2);
    int irst = 0;
	r1:   irst = irst + 1;
    asub(x, xi);
    double rp = 0.0;
    double bsq = 0.0;
    for (j = 1; j<=n; j++)
	{
        bsq = bsq + pow(b[j],2);
        xi[j] = xi[j] - b[j];
        rp = rp + pow(xi[j],2);
    }
    atsub(xi, g);
    for (j = 1; j<=n; j++)
	{
        g[j] = -g[j];
        h[j] = g[j];
    }
    for (int iter = 1; iter<=10 * n; iter++)
	{
        asub(h, xi);
        double anum = 0.0;
        double aden = 0.0;
        for (j = 1; j<=n; j++)
		{
            anum = anum + g[j] * h[j];
            aden = aden + pow(xi[j] ,2);
        }
        if (aden == 0.0)
		{
			cout<<"very singular matrix"<<endl;
		}
        anum = anum / aden;
        for (j = 1; j<=n; j++)
		{
            xi[j] = x[j];
            x[j] = x[j] + anum * h[j];
        }
        asub(x, xj);
        double rsq = 0.0;
        for (j = 1; j<=n; j++)
		{
            xj[j] = xj[j] - b[j];
            rsq = rsq + pow(xj[j], 2);
        }
        if (rsq == rp || rsq <= bsq * eps2)
		{
			return;
		}
        if (rsq > rp)
		{
            for (j = 1; j<=n; j++)
			{
                x[j] = xi[j];
            }
            if (irst >= 3)
			{
                return;
			}
            goto r1;
        }
        rp = rsq;
        atsub(xj, xi);
        double gg = 0.0;
        double dgg = 0.0;
        for (j = 1; j<=n; j++)
		{
            gg = gg + pow(g[j],2);
            dgg = dgg + (xi[j] + g[j]) * xi[j];
        }
        if (gg == 0.0)
		{
			return;
		}
        double gam = dgg / gg;
        for (j = 1; j<=n; j++)
		{
            g[j] = -xi[j];
            h[j] = g[j] + gam * h[j];
        }
    }
    cout<<"too many iterations"<<endl;
}