最速下降 最优化

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
#include <string.h>
#define m 2
#define epsilon   pow(10,-5)

double norm_1(double *x);
double fun_obj(double *x);//原函数
void fun_grad(double *x,double *grad);//求梯度：
void armijograd(double *x0,double *grad);
int main()
{
    int i;
    double *x;
    x = (double*)malloc(sizeof(double)*m);
    x[0] = 3.0;
    x[1] = 1.0;
    double *grad;
    grad = (double*)malloc(sizeof(double)*m);
    grad[0] = 0.0;
    grad[1] = 0.0;

    armijograd(x,grad);
    for(i = 0;i<m;i++)
        printf("%f ",x[i]);
    free(x);free(grad);
    return 0;
}

double fun_obj(double *x)//原函数
{
    double f = 0.0;
    f = pow(x[0]-1,2)+pow(x[1]-1,2);
    return f;
}

void fun_grad(double *x,double *grad)//求梯度：
{
    grad[0] = 2*(x[0]-1);
    grad[1] = 2*(x[1]-1);

}

double norm_1(double *x)
{
    double sum = 0.0;
    int i;
    for(i=0;i<m;i++)
        sum+=x[i]*x[i];
    sum = sqrt(sum);
    return sum;
}

void armijograd(double *x0,double *grad)
{
    int max_iter = 5000;   // max number of iterations
    double EPS = 1e-6;        //threshold of gradient norm

    double rho = 0.45; double sigma = 0.2;  // Armijo parameters

    double *xk;xk = (double *)malloc(sizeof(double)*m);
    int i;
    for(i=0;i<m;i++)
        xk[i] = x0[i]; // initialization
    int k = 0;
    double *d;d= (double *)malloc(sizeof(double)*m);
    int mk = 0; int max_mk = 100;
    while(k < max_iter)
        {
            k = k + 1;
            fun_grad( xk ,grad);   // gradient vector
            for(i=0;i<m;i++)
                d[i] = -1 * grad[i];              // search direction
            if( norm_1(d) < EPS )  //precision
                break;
            while (mk <= max_mk)//armijo line search
            {
                fun_grad(xk,grad);//grad存放xk的梯度向量
                for(i=0;i<m;i++)
                    d[i] = -1 * grad[i];
                for(i = 0;i<m;i++)
                    x0[i] = xk[i] + pow(rho,mk) * d[i];
                double temp = 0.0;
                for(i = 0;i<m;i++)
                    temp+=grad[i]*d[i];
                if( fun_obj( x0 ) <=  fun_obj( xk ) + sigma * pow(rho,mk) *temp)
                    break;
                mk = mk + 1;
            }
            for(i = 0;i<m;i++)
                xk[i] = xk[i] + pow(rho,mk)*d[i];  //update
        }

    for(i = 0;i<m;i++)
        x0[i] = xk[i];
}