cg_main.cpp

C++源文件 没有建立工程 下载后请自行建立

/***************************************************
*    Conjugate Gradient Solver         v0.2        *
*    Programmed by Xiangyu Dong                    *
*--------------------------------------------------*
*  v0.2 is a prototype, and the target function is *
*  fixed. For convenience, it is set as below:     *
*      f(x)=(6+x1+x2)^2+(2-3x1-3x2-x1x2)^2         *
*  This function is described in double func().    *
*  And, the initial point is also given:           *
*      x0=(-4,6)                                   *
*  Of course, this limitation will be eliminated   *
*  in further versions.                            *
***************************************************/

#include <stdio.h>
#include <math.h>
#include "matrix.h"

#define SEARCH_RANGE 1000
#define GOLDEN_RATIO 0.61803398874989
#define ERROR 1e-6

// func: target function, here, result=(6+x1+x2)^2+(2-3x1-3x2-x1x2)^2
double func(double *x, long n)
{
	double result;
	double x1, x2;
	x1 = x[1];
	x2 = x[2];
	result = (6+x1+x2)*(6+x1+x2)+(2-3*x1-3*x2-x1*x2)*(2-3*x1-3*x2-x1*x2);
	return result;
}

double func1(double *x, long n)
{
	double result;
	double x1, x2;
	x1 = x[1];
	x2 = x[2];
	result = (x1-2)*(x1-2)+2*(x2-1)*(x2-1);
	return result;
}

// d1: derivative calculation, result=df(x)/dxi
double d1(double (*func)(double *, long), double *x, long n, long i)
{
	double h=5e-15;
	double *new_x;
	double result_l, result_r;
	long j;

	new_x = dvector(1, n);
	for (j=1; j<=n; j++)
		new_x[j] = x[j];
	
	new_x[i] -= h;
	result_l = (*func)(new_x,n);
	new_x[i] += 2*h;
	result_r = (*func)(new_x,n);

	free_dvector(new_x, 1, n);

	return (result_r-result_l)/2/h;
}

// 1-dimension search using 0.168 method
double search_0p618(double (*func)(double *, long), double *y, double *d, long n)
{
	double a, b;				// lower bound and upper bound of serach range
	double l, h;				// Two test points
	double fl, fh;				// The function value at two test points
	double *new_y;				// temp
	long i;
	double original;

	original = (*func)(y,n);
	
	a = 0; b = SEARCH_RANGE;	// Initialize the search bound
	do{
		l = a+(1-GOLDEN_RATIO)*(b-a);
		h = a+(GOLDEN_RATIO)*(b-a);
		new_y = dvector(1, n);

		for (i=1; i<=n; i++)
			new_y[i] = y[i]+l*d[i];
		fl = (*func)(new_y,n);

		for (i=1; i<=n; i++)
			new_y[i] = y[i]+h*d[i];
		fh = (*func)(new_y,n);

		while (abs(fl-fh)>ERROR){
			if (fl>fh){
				a = l;
				l = h;
				h = a+(GOLDEN_RATIO)*(b-a);
			}else{
				b = h;
				h = l;
				l = a+(1-GOLDEN_RATIO)*(b-a);
			}
			for (i=1; i<=n; i++)
				new_y[i] = y[i]+l*d[i];
			fl = (*func)(new_y,n);
			for (i=1; i<=n; i++)
				new_y[i] = y[i]+h*d[i];
			fh = (*func)(new_y,n);
		}
		a = 0;
		b = l;
	}while (fl>original);

	free_dvector(new_y, 1, n);

	return l;
}

void main()
{
	int c;
	long n;						// The dimension of problem space
	double *x;					// This is a vector storing current solution point
	double *y;					// Intermediate variables in CG algorithm, refer to textbook page 300
	double *dy;					// The derivative of y
	double *d;					// Intermediate variables in CG algorithm, refer to textbook page 300
	double e;					// Tolerance error
	long i;						// loop variable
	double norm_current;
	double norm_next;
	double beta;
	double step;				// 1-dimension search step
	long count = 0;

	e = 1e-15;
	n = 2;						// In v0.2, this value is fixed, since there are only two varibles
	x = dvector(1, n);			// Set as a n-dimension vector
	y = dvector(1, n);
	dy = dvector (1, n);
	d = dvector(1, n);
	
	// Initialization 
	x[1] = -4;
	x[2] = 6;					// Set the initial point as (-4,6)

	// Print the initial solution
	printf("The initial solution is:\n");
	for (i=1; i<=n; i++)
	printf("x%ld=%lf\t", i,x[i]);
	double result;
	result = func(x, n);
	printf("\nObjective Value = %lf\nPress ENTER to optimize ...\n",result);
	fflush(stdin);
	c=fgetc(stdin);
	printf("Processing ...\n\n");

	// Step 1
	for (i=1; i<=n; i++){
		y[i] = x[i];
		d[i] = 0 - d1(func, x, n, i);
	}
	norm_next = norm_dvector(d, n, i);
	do{
		// Step 2
		norm_current = norm_next;
		// 1-dimension search
		step = search_0p618(func, y, d, n);
		if (abs(step)<1e-15)	// d is not a descend direction
			for (i=1; i<=n; i++)
				d[i] = 0-d1(func, y, n, i);
		else{					// d is a descend direction
			// update y: y=y+step*d
			for (i=1; i<=n; i++)
				y[i] = y[i]+step*d[i];
			for (i=1; i<=n; i++)
				dy[i] = d1(func, y, n, i);
			norm_next = norm_dvector(dy, 1, n);
			beta = norm_next*norm_next/norm_current/norm_current;
			// update d: d=-d1(y)+beta*d
			for (i=1; i<=n; i++)
				d[i] = beta*d[i]-dy[i];
		}
		
		// Print the intermediate solution
		if (++count ==100){
			for (i=1; i<=n; i++)
				printf("x%ld=%lf\t", i,y[i]);
			result = func(y, n);
			printf("\nObjective Value = %lf\n",result);
			printf("Processing ...\n\n");
			count = 0;
		}
	}while (norm_current>e);

	printf("Optimization finished!\n");
	for (i=1; i<=n; i++)
		printf("x%ld=%lf\t", i,y[i]);
	result = func(y, n);
	printf("\nObjective Value = %lf\n",result);
	fflush(stdin);
	c=fgetc(stdin);

	free_dvector(x, 1, n);		// Memory free
	free_dvector(y, 1, n);		// Memory free
	free_dvector(dy, 1, n);		// Memory free
	free_dvector(d, 1, n);		// Memory free
}