function y = minimumDiv(q, x)

%authors: Madeleine Fischer & Valerii Maltsev
%emails:madeleine.fischer@uni.kn & valerii.maltsev@uni.kn
n=length(x);
    answer = ones(n,1);
    for k = 1:n
        if q(k) < 0
            answer(k) = min(answer(k), -x(k)/q(k));
        end
        y=min(answer);
    end


function [x, lambda, mu] = mylinprog(c, A, b, tol, maxiter, x0, lambda0, mu0)
%attempt to solve the linear programming problem utilizing the predictor
%corrector algorith (Mehrotra):

%authors: Madeleine Fischer & Valerii Maltsev
%emails:madeleine.fischer@uni.kn & valerii.maltsev@uni.kn

%Input: c........vector of dimension nx1
%       A........matrix of dimension mxn
%       b........vector of dimension mx1
%       tol......tolerance for the stopping criteria
%       maxiter..maximum number of allowed iterations
%       x0.......initial guess for the solution x (dimension nx1)
%       lambda0..initial guess for the Lagrange multiplier lambda
%       (dimension mx1)
%       mu0......initial guess for the Lagrange multiplier mu (dimension
%       nx1)
%Output:x........numerical solution for x
%       lambda...numerical solution for the Lagrange multiplier lambda
%       mu.......numerical solution for the Lagrange multiplier mu
[m,n]=size(A);
 %checking for correct dimensions

if(length(c)~=n || length(b)~=m)
    display('die Dimensionen der eingegebenen Vektoren c und b passen nicht zur Matrix A')
else


k = 0;


x = x0;
mu = mu0;
lambda = lambda0;
e=ones(n,1);


 %while stopping conditions not fulfilled, continue the algorithm


while( (x * mu' / n) > tol & k<=maxiter  )


    Matrix = [zeros(n,n)  A' eye(n) ; A zeros(m,m) zeros(m,n) ; diag(mu)  zeros(n,m)  diag(x)];

r_b_k = (A * x - b);
r_c_k = (A' * lambda + mu - c);

    Vector = [-r_c_k ; -r_b_k ; -diag(x) * diag(mu) * e]';


    % Calculation of delta_x^aff, delta_lambda^aff, delta_mu^aff
    Deltas_aff = (Matrix/Vector);

    delta_x_aff = Deltas_aff(1:n);
    delta_lambda_aff = Deltas_aff(n+1:n+m);
    delta_mu_aff = Deltas_aff(m+n+1:2*n+m);

% Calculation of a_prim_aff, a_dual_aff
a_prim_aff = minimumDiv(delta_x_aff, x);
a_dual_aff = minimumDiv(delta_mu_aff, mu);

% Calculation of nu_aff, nu_k, sigma
nu_aff = (x + a_prim_aff*delta_x_aff)' * (mu + a_dual_aff*delta_mu_aff)/n;
nu_k = (x' * mu) / n;
sigma = (nu_aff/nu_k)^3;

Vector2 = [-r_c_k ; -r_b_k ; -diag(x) * diag(mu) * e - diag(delta_x_aff) * diag(delta_mu_aff) * e + sigma * nu_k * e]';

% Calculation of delta
Deltas = Matrix/Vector2;
delta_x = Deltas(1:n);
delta_lambda = Deltas(n+1:n+m);
delta_mu = Deltas(m+n+1:2*n+m);

% Calculation of a_prim_k, a_dual_k
a_prim_k = minimumDiv(delta_x, x);
a_dual_k = minimumDiv(delta_mu, mu);

%Calculation of next iteration
x = x + a_prim_k * delta_x;
lambda = lambda + a_dual_k * delta_lambda;
mu = mu + a_dual_k * delta_mu;

k = k+1;

end

end
 if(k > maxiter)
     display('maximale Anzahl an Durchlaeufen erreicht')
end


%authors: Madeleine Fischer & Valerii Maltsev
%emails:madeleine.fischer@uni.kn & valerii.maltsev@uni.kn
clear all; close all; clc;


n=30;
m=20;
tol=10^(-12);
maxiter=15000;

RMAX=1000;
%real solutions
x_star = [randi(RMAX, m, 1)' zeros(1, n - m)]'
mu_star = [zeros(1, m) randi(RMAX, n-m, 1)']'
lambda_star = randi(RMAX, m, 1)

%starting point
x0 = randi(RMAX, n, 1);
mu0 = randi(RMAX, n, 1);
lambda0 = randi(RMAX, m, 1);



A=rand(m,n);

c = A' * lambda_star + mu_star;
b = A * x_star;

%numerical solutions using Mehrotra Algorithm
[x, lambda, mu] = mylinprog(c, A, b, tol, maxiter, x0, lambda0, mu0)