powell.cpp

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

void erase(double pbar[], double prr[], double pr[])
{
	for (int i=0; i<=20; i++)
	{
		pbar[i]=0;
		prr[i]=0;
		pr[i]=0;
	}
}

void powell(double p[], double xi[], int n, int np, double ftol,
			int& iter, double& fret)
{
    int i,j,ibig,itmax = 200;
    double pt[21], ptt[21], xit[21];
	double t,dum,fp,del,fptt;
    fret = func2(p, n);
    for (j = 1; j<=n; j++)
	{
        pt[j] = p[j];
    }
    iter = 0;
    do
	{
		do
		{
			do
			{
				iter = iter + 1;
				fp = fret;
				ibig = 0;
				del = 0.0;
				for (i = 1; i<=n; i++)
				{
					for (j = 1; j<=n; j++)
					{
						xit[j] = xi[(j-1)*np+i];
					}
					fptt = fret;
					linmin(p, xit, n, fret);
					if (fabs(fptt - fret) > del)
					{
						del = fabs(fptt - fret);
						ibig = i;
					}
				}
				if (2.0*fabs(fp-fret)<=ftol*(fabs(fp)+fabs(fret)))
				{
					erase(xit, ptt, pt);
					return;
				}
				if (iter == itmax)
				{
					cout<<"powell exceeding maximum iterations"<<endl;
					return;
				}
				for (j = 1; j<=n; j++)
				{
					ptt[j] = 2.0 * p[j] - pt[j];
					xit[j] = p[j] - pt[j];
					pt[j] = p[j];
				}
				fptt = func2(ptt, n);
			}while (fptt >= fp);
			dum = fp - 2 * fret + fptt;
			t = 2.0 * dum*pow((fp-fret-del),2)-del*pow((fp-fptt),2);
		}while (t >= 0.0);
		linmin(p, xit, n, fret);
		for (j = 1; j<=n; j++)
		{
			xi[(j-1)*np+ibig] = xit[j];
		}
	}while(1);
}