test.cpp

数学计算程序

/* We illustrate the operation of the code using the problem

Final convergence status = 0
Convergence tolerance for gradient satisfied
projected gradient max norm: 7.007567e-009
function value:              -4.234085e+000

Total cg  iterations:                    6
Total cg  function evaluations:         13
Total cg  gradient evaluations:          7
Total cbb iterations:                    3
Total cbb function evaluations:          8
Total cbb gradient evaluations:          4
------------------------------------------
Total function evaluations:             21
Total gradient evaluations:             11
==========================================

0  0  4.67391  0  0  2.63256
2.04134
 
See the drivers in the distribution of cg_descent Version 3.0
for more examples */

#include <math.h> 
#include "asa_user.h" /* needed by the program which calls asa_cg */
#include "../PackingGA/pair_sc.h"
#include "../PackingGA/pair_sc.cpp"

/* prototypes for the function and gradient evaluation routines */
double myvalue(double *x, INT32 n) ;

void mygrad(double  *g, double  *x, INT32 n) ;

double myvalgrad(double  *g, double  *x, INT32 n) ;

int main (void)
{
    double *x, *lo, *hi ;
	double grad_tol;
    INT32 i, n ;

    /* if you want to change parameter value, you need the following: */
    //asacg_parm cgParm ;
    //asa_parm asaParm ;

	/* tolerance */
	grad_tol = 1.e-8;
	//grad_tol = 0.01;
    /* allocate arrays for problem solution and bounds */
    n = 6; /* problem dimension */
    x  = (double *) malloc (n*sizeof (double)) ;
    lo = (double *) malloc (n*sizeof (double)) ;
    hi = (double *) malloc (n*sizeof (double)) ;
    for (i = 0; i < n; i++) lo [i] = (double) 0 ;
    for (i = 0; i < n; i++) hi [i] = (double) 10 ;

    ///* if you want to change parameter value, initialize strucs with default */
    //asa_cg_default (&cgParm) ;
    //asa_default (&asaParm) ;

    ///* if you want to change parameters, change them here: */
    //cgParm.PrintParms = TRUE ;
    //cgParm.PrintLevel = 1 ;
    //asaParm.PrintParms = TRUE ;
    //asaParm.PrintLevel = 1 ;

    ///* starting guess */
    //for (i = 0; i < n; i++) x[i] = 1 ;
    /////* run the code */
    //asa_cg (x, lo, hi, n, NULL, &cgParm, &asaParm, grad_tol, myvalue, mygrad, myvalgrad, NULL) ;

    /* if no change in parameters, you could replace Parm arguments by NULL*/
	/* starting guess */
    //for (i = 0; i < n; i++) x [i] = 1 ;
    //asa_cg (x, lo, hi, n, NULL, NULL, NULL, grad_tol, myvalue, mygrad, myvalgrad, NULL) ;

    /* with some loss of efficiency, you could omit the valgrad routine */
	/* starting guess */
	//for (i = 0; i < n/2; i++) x [i] = 0 ;
	for (i = 0; i < n-1; i++) x [i] = 0 ; 
	x[n-1] = 5;
	//for (i = n/2; i < n; i++) x[i] = 5.0; 
	//setup_param(string("Ni"));
	//Rc = 100;
    asa_cg (x, lo, hi, n, NULL, NULL, NULL, grad_tol, myvalue, mygrad, NULL, NULL);

	for(i = 0; i < n; i++)
		cout << x[i] << "  ";
	cout << endl;
	cout << sqrt((x[0]-x[3])*(x[0]-x[3]) + (x[1]-x[4])*(x[1]-x[4]) + (x[2]-x[5])*(x[2]-x[5])) << endl;

    free (x) ;
    free (lo) ;
    free (hi) ;
}

/* evaluate the objective function */
double myvalue(double *xx, INT32 n)
{
    double f;
    double *x = new double[n/3];
	double *y = new double[n/3];
	double *z = new double[n/3];
	for(int i = 0; i < n;)
	{
		x[i/3] = xx[i];
		y[i/3] = xx[i+1];
		z[i/3] = xx[i+2];
		i += 3;
	}
	setup_param(string("Ni"));
	//Rc = 100;
	f = EnergySC(x, y, z, n/3, string("Ni"));
	delete []x;
	delete []y;
	delete []z;
    return (f) ;
}

/* evaluate the gradient of the objective function */
void mygrad(double  *g, double  *xx, INT32 n)
{
	int i;
    double *x = new double[n/3];
	double *y = new double[n/3];
	double *z = new double[n/3];
	double **f = new double*[3];
	for(i = 0; i < 3; i++) f[i] = new double[n/3];
	for(i = 0; i < n;)
	{
		x[i/3] = xx[i];
		y[i/3] = xx[i+1];
		z[i/3] = xx[i+2];
		i += 3;
	}
	setup_param(string("Ni"));
	//Rc = 100;
	ForceSC(x, y, z, n/3, f, string("Ni"));
	for(i = 0; i < n;)
	{
		g[i] = -f[0][i/3];
		g[i+1] = -f[1][i/3];
		g[i+2] = -f[2][i/3];
		i += 3;
	}
	delete []y;
	delete []z;
	delete []x;
	for(i = 0; i < 3; i++) delete []f[i];
	delete []f;
}

/* value and gradient of the objective function */
double myvalgrad (double  *g, double  *x, INT32 n)
{
    double f, xi, t ;
    INT32 i ;
    f = 0 ;
    for (i = 0; i < n; i++)
    {
        t = i + 1 ;
        t = sqrt (t) ;
        xi = x [i] ;
        f += exp (xi) - t*xi ;
        g [i] = exp (xi) -  t ;
    }
    return (f) ;
}