DPF

/*  Function Rosenbrock
F(x,y)=(1-x)^2+100(y-x^2)^2
partial derivative x : f = -2+2*x-400*y*x+400*x^3
partial derivative y : f = 200*(y-x^2)
*/

#include <conio.h>
#include <math.h>
#include <iostream>

using namespace std;

double function_Rosenbroke(double * vector) {
    return pow((1 - vector[0]), 2) + 100 * pow((vector[1] - pow(vector[0], 2)), 2);
}
//zeroing matrix

double ** zeroing(double **A,int N) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++)
            A[i][j] = 0;
    }
    return A;
}
//multiply matrixes

double** product(double **A, double **U,int N) {

    double **c = new double*[N];
    for (int i = 0; i < N; i++)
        c[i] = new double[N];

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            c[i][j] = 0;
            for (int t = 0; t < N; t++)
                c[i][j] += A[i][t] * U[t][j];
        }
    }
    return c;
}
// multiply matrix on vector

double* product_vector(double **A, double *b, int N) {

    double *a = new double[N];

    for (int i = 0; i < N; i++) {
        a[i] = 0;
        for (int j = 0; j < N; j++)
            a[i] += A[i][j] * b[j];
    }
    return a;
}
// find gradients
double grad1(double* vector, int N) {
    return -2 + 2 * vector[0] - 400 * vector[1] * vector[0] + 400 * pow(vector[0], 3);  //partial derivative x
}
double grad2(double* vector, int N) {
    return 200 * (vector[1] - pow(vector[0], 2)); //partial derivative y
}

//find direction

double * func_direction(double **A, double *b, int N) {

    double *direction = new double[N];
    for (int i = 0; i < N; i++)
        direction[i] = 0;

    direction = product_vector(A, b,  N);

    for (int i = 0; i < N; i++)
        direction[i] = -direction[i];
    return direction;
}

double **make_matrix(double * a, double * b, int N) {
    double **matrix = new double*[N];
    for (int i = 0; i < N; i++)
        matrix[i] = new double[N];
    zeroing(matrix,N);

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++)
            matrix[i][j] = a[i] * b[j];
    }
    return matrix;
}

double make_number(double * a, double *b, int N) {
    double number = 0;
    for (int i = 0; i < N; i++)
        number += a[i] * b[i];
    return number;
}
//division matrix by number
double ** division_matrix_by_number(double ** matrix, double number, int N) {

    double **matrix1 = new double*[N];
    for (int i = 0; i < N; i++)
        matrix1[i] = new double[N];
    zeroing(matrix1,N);

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++)
            matrix1[i][j] = matrix[i][j] / number;
    }
    return matrix1;
}
 //find reverse Hesse's  matrix
double ** formula(double ** matrix, double *alpha, double *beta, int N) {
    double number1 = 0, number2 = 0;

    double **matrix1 = new double*[N];
    for (int i = 0; i <N; i++)
        matrix1[i] = new double[N];
    zeroing(matrix1,N);

    double **matrix2 = new double*[N];
    for (int i = 0; i < N; i++)
        matrix2[i] = new double[N];
    zeroing(matrix2,N);

    double **B = new double*[N];
    for (int i = 0; i < N; i++)
        B[i] = new double[N];
    zeroing(B,N);

    double **C = new double*[N];
    for (int i = 0; i < N; i++)
        C[i] = new double[N];
    zeroing(C,N);

    double *a = new double[N];
    for (int i = 0; i < N; i++)
        a[i] = 0;

    double *b = new double[N];
    for (int i = 0; i < N; i++)
        b[i] = 0;

    matrix1 = make_matrix(alpha, alpha,N);
    number1 = make_number(alpha, beta,N);
    B = division_matrix_by_number(matrix1, number1,N);
    a = product_vector(matrix, beta,N);
    matrix2 = make_matrix(a, beta,N);
    matrix2 = product(matrix2, matrix,N);
    b = product_vector(matrix, beta,N);
    number2 = make_number(b, beta,N);
    C = division_matrix_by_number(matrix2, number2,N);

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++)
            matrix[i][j] = matrix[i][j] + B[i][j] - C[i][j];
    }

    return matrix;
}

double length_of_grad(double *gradient) {
    return sqrt(pow(gradient[0], 2) + pow(gradient[1], 2));
}

int main() {
    int N = 2;
    double k = 0, k_max, eps = 0.00016, t, number1, number2;
    cout << "enter max approximation: ";
    cin >> k_max;
    cout << "enter point of min 0<t<1: ";
    cin >> t;

    double **A = new double*[N]; //identity matrix
    for (int i = 0; i < N; i++)
        A[i] = new double[N];

    double *direction = new double[N];
    for (int i = 0; i < N; i++)
        direction[i] = 0;

    for (int i = 0; i < N; i++) { // enter identity matrix
        for (int j = 0; j < N; j++) {
            if (i == j)
                A[i][j] = 1;
            else
                A[i][j] = 0;
        }
    }

    double *vector = new double[N];
    for (int i = 0; i < N; i++)
        vector[i] = 0;
    cout << "enter points x and y: ";
    cin >> vector[0] >> vector[1];

    double *gradient = new double[N];
    for (int i = 0; i < N; i++)
        gradient[i] = 0;

    double *alpha = new double[N];
    for (int i = 0; i < N; i++)
        alpha[i] = 0;

    double *beta = new double[N];
    for (int i = 0; i < N; i++)
        beta[i] = 0;

    gradient[0] = grad1(vector,N);    // function's gradient evaluation at a point X
    gradient[1] = grad2(vector,N);    // function's gradient evaluation at a point Y

    do {
        number1 = 0, number2 = 0;

        direction = func_direction(A, gradient,N);            //to find direction

        for (int i = 0; i < N; i++)                        // old points
            alpha[i] = vector[i];

        for (int i = 0; i < N; i++)                        // new points
            vector[i] = vector[i] + t * direction[i];

        for (int i = 0; i < N; i++)                        // new points - old points
            alpha[i] = vector[i] - alpha[i];

        for (int i = 0; i < N; i++)                       //old values of grad
            beta[i] = gradient[i];

        gradient[0] = grad1(vector,N);                      // new values of grads
        gradient[1] = grad2(vector,N);

        for (int i = 0; i < N; i++)                         //new grad - old grad
            beta[i] = gradient[i] - beta[i];

        A = formula(A, alpha, beta,N);                          //Ak=Ak-1 + B - C;

        k++;

        cout << endl << "grad " << length_of_grad(gradient) << " " << k << endl;

    } while ((length_of_grad(gradient) > eps) && (k < k_max));

    cout << vector[0] << " " << vector[1] << endl;
    cout << function_Rosenbroke(vector);

    _getch();
    return 0;
}