void dfpmin(double p[],int n,double ftol,int& iter,double& fret)
{
    int i,j,its,itmax = 200;
    double eps = 0.0000000001;
    double hessin[51][51], xi[51], g[51], dg[51], hdg[51];
    double fp = func2(p, n);
    dfunc(p, g);
    for (i = 1; i<=n; i++)
	{
        for (j = 1; j<=n; j++)
		{
            hessin[i][j] = 0.0;
        }
        hessin[i][i] = 1.0;
        xi[i] = -g[i];
    }
    for (its = 1; its<=itmax; its++)
	{
        iter = its;
        linmin(p, xi, n, fret);
        if (2.0 * fabs(fret - fp)<=ftol*(fabs(fret)+fabs(fp)+eps))
		{
            return;
        }
        fp = fret;
        for (i = 1; i<=n; i++)
		{
            dg[i] = g[i];
        }
        fret = func2(p, n);
        dfunc(p, g);
        for (i = 1; i<=n; i++)
		{
            dg[i] = g[i] - dg[i];
        }
        for (i = 1; i<=n; i++)
		{
            hdg[i] = 0.0;
            for (j = 1; j<=n; j++)
			{
                hdg[i] = hdg[i] + hessin[i][j] * dg[j];
            }
        }
        double fac = 0.0;
        double fae = 0.0;
        for (i = 1; i<=n; i++)
		{
            fac = fac + dg[i] * xi[i];
            fae = fae + dg[i] * hdg[i];
        }
        fac = 1.0 / fac;
        double fad = 1.0 / fae;
        for (i = 1; i<=n; i++)
		{
            dg[i] = fac * xi[i] - fad * hdg[i];
        }
        for (i = 1; i<=n; i++)
		{
            for (j = 1; j<=n; j++)
			{
				double aaa = fac*xi[i]*xi[j]-fad*hdg[i]*hdg[j];
				double bbb = fae*dg[i]*dg[j];
				hessin[i][j] = hessin[i][j] + aaa + bbb;
            }
        }
        for (i = 1; i<=n; i++)
		{
            xi[i] = 0.0;
            for (j = 1; j<=n; j++)
			{
                xi[i] = xi[i] - hessin[i][j] * g[j];
            }
        }
    }
    if (its > itmax)
	{
		cout<<" too many iterations in dfpmin"<<endl;
	}
}