function [ alpha, FoundValue, exitflag ] = solve_quadprog_both( M, V, AlphaStart, sDM, TShadow, Accuracy, Options, ZeroVec, OneVec )

% Function to solve for alpha in the quadratic optimisation problem with
% ONBC
%
%
% INPUTS 
% ------
% M:            Matrix formed from horizontally concatenating column vectors of the relative impulse responses of the bounded variable
% V:            Matrix of the unconstrained absolute IRFs of the constrained variable
% AlphaStart:   Start value for alpha
% sDM:			Original signs of the diagonal of M
% TShadow:      Number of periods, the constraint is non-binding
% Accuracy:		[Optional] Permitted amount of constraint violation
% Options:		[Optional] quadprog options structure
% ZeroVec:		[Optional] Optional argument to save on unnecessary memory allocations
% OneVec:		[Optional] Optional argument to save on unnecessary memory allocations
%
% OUTPUTS
% -------
% alpha:        Vector that determines the size for each shadow shock to keep the constrained variable bounded by zero

	% Get an initial approximation for alpha to increase the stability and the speed of the computation

	if nargin < 7
		ZeroVec = zeros( TShadow, 1 );
		OneVec = ones( size( ZeroVec ) );
		if nargin < 6
			Options = optimset;
			if nargin < 5
				Accuracy = sqrt( eps );
			end
		end
	end
	M_Sub = M( 1:2*TShadow, : );
	V_Sub = V( 1:2*TShadow );
 	alpha = AlphaStart * max( 0, max( -V_Sub ) );
    V_Candidate = V_Sub + M_Sub * alpha;
	indices = [ ];
	for conv_iter = 1:TShadow
		%indices = union( indices, find( V_Candidate < -Accuracy ) );
		if isempty( indices )
			break;
		end
		oldalpha = alpha;
		alpha( indices ) = max( ZeroVec( indices ), alpha( indices ) - M_Sub( indices, indices ) \ V_Candidate( indices ) );
		indices = setdiff( indices, find( alpha < Accuracy ) );
		V_Candidate = V_Sub + M_Sub * alpha;
		if min( V_Candidate ) >= -Accuracy || norm( oldalpha - alpha ) < Accuracy
			break;
		end
	end
    
	if min( V_Candidate ) < -Accuracy
		alpha = max( ZeroVec, linprog( OneVec, -M_Sub, V_Sub, [ ], [ ], ZeroVec, [ ], [ ], optimset( 'LargeScale', 'on', 'MaxIter', 1000, 'TolFun', Accuracy ) ) ); 
		V_Candidate = V_Sub + M_Sub * alpha;
	end
	

	% Solve the quadratic optimisation problem
    warning('off','optim:quadprog:WillBeRemoved');
    [ alpha, FoundValue, exitflag ] = quadprog( M + M', V, -M, V, [ ], [ ], ZeroVec, [ ], alpha, Options );
    
	if FoundValue >= Accuracy || exitflag <= 0
		error( 'Model7:quadprog', 'Failed to solve quadratic programming problem. Try increasing TShadow.' );
	end
	if alpha( TShadow ) > (Accuracy*100000000)
		warning( 'Model7:Inaccuracy', 'The final component of alpha is equal to %e > 0. This is indicative of TshadowI being too small.', alpha( end ) );
	end
	alpha = alpha .* sDM;	
end