function [solution, kuhn_tucker, iterations, end_condition] = ...
    dykstra(data,covariance,constraint_array,constraint_constant,equality_flag)

zero = 0.0;
two = 2.0;
end_condition = 0;
iterations = 0;
n = length(data);
k = size(constraint_array,1);
eps = 1.0e-6;
eps2 = 1.0e-10;
itermax = 1.0e+04;

nk = n*k;
nktwo = nk + nk;
nktwok = nktwo + k;
nwork = nktwok + n;

work = zeros(nwork,1);
solution = zeros(n,1);
kuhn_tucker = zeros(k,1);

if(n <= 0)
    end_condition = 2;
    return
end

if(k <= 0)
    end_condition = 2;
    return
end


for i = 1:k
    if((equality_flag(i) ~= 0)  & (equality_flag(i) ~= 1))
        end_condition = 2;
        return
    end
end


for i = 1:k
    for j = 1:n
        index = nk + (j-1)*k + i;
        work(index) = zero;
         for h = 1:n
            work(index) = work(index) + ...
            constraint_array(i,h)*covariance(h,j);
         end
    end
end

for i = 1:k
    index = nktwo + i;
    work(index) = zero;
    for j = 1:n
        indtwo = nk + (j-1)*k + i;
        work(index) = work(index) + constraint_array(i,j)*work(indtwo);
    end
    if (abs(work(index)) <= eps2)
        end_condition = 3;
        return
    end
end

for j = 1:n
    index = nktwok + j;
    work(index) = data(j);
    solution(j) = data(j);
    for i = 1:k
        index = (i-1)*n + j;
        work(index) = zero;
    end
end

iterations = 0;

while (iterations <= itermax)
    
    iterations = iterations + 1;
    supdif = zero;
    
    for i = 1:k
        for j = 1:n
            index = (i-1)*n + j;
            solution(j) = solution(j) + work(index);
        end
        sum = zero;
        for j = 1:n
            sum = sum + constraint_array(i,j)*solution(j);
        end
        
        sum = sum - constraint_constant(i);
        
        if((equality_flag(i) == 0) & (sum <= zero)) 
            for j = 1:n
                index = (i-1)*n + j;
                work(index) = zero;
            end
        else
            index = nktwo + i;
            temp = sum/work(index);
            for j = 1:n
                indone = (i-1)*n + j;
                indtwo = nk + (j-1)*k + i;
                work(indone) = work(indtwo)*temp;
                solution(j) = solution(j) - work(indone);
            end
        end
        
        for j = 1:n
            index = nktwok + j;
            abdif = abs(solution(j) - work(index));
            if(supdif < abdif) 
                supdif = abdif;
            end
        end
    end
    
    if (supdif <= eps)
        for i = 1:k
            for j = 1:n
                index = nk + (j-1)*k + i;
                if (abs(work(index)) > eps2)
                    indtwo = (i-1)*n + j;
                    kuhn_tucker(i) = two*work(indtwo)/work(index);
                end
            end
        end
        
        return
        
    end
    
   for j = 1:n
       index = nktwok + j;
       work(index) = solution(j);
   end
end

end_condition = 1;

return

            
            
        
        
            
            
