function [p_f]= reducerightside( p_g, G_mn, K_fs, y_s, flag_y_s, uset )
% REDUCERIGHTSIDE [p_f]= Reducerightside( p_g, G_mn, K_fs, y_s, flag_y_s )
% Reduces right hand side vectors from g-size to f-size (p_g to p_f).
%
% Partition p_g into p_n and p_m, then modify p_n
%
% p_n = p_n + G_mn^T * p_m
%
% where G_mn is the reduction matrix from multi-point constraints
%
% Partition p_n into p_f and p_s, then modify p_f if flag_y_s > 0
%
% p_f = p_f - K_fs * y_s
%
% where y_s are enforced displacements/temperatures, 
%
% The g-set is partitioned as follows
%   g = n + m
%  		n = f + s
%   g 	all system degrees of freedom (dof)
%   m 	multi-point constraint dofs to be condensed out
%   n   remaining dof after m-set is condensed from the g-set
%   s   boundary condition dofs (single point constraints)
%   f   remaining free dofs after s-set is eliminated from the n-set
%
% Input:  p_g     	right hand side vectors (g size) 
%		  G_mn      Reduction matrix from constraints (m x n)
%		  K_fs		Partitioned stiffness matrix (f x s)
%         y_s		enforced displacements / temperatures s size
%		  flag_y_s  flag for enforced displacement / temperature
%					> 0 	K_fs * y_s is subtracted
%                   <= 0	no action
%
% Output: p_f       reduced right hand side vectors (f-size)
%
% Called by  sol101 (ms), Normalmodes (ms), Residualheat (mf)
% Calls none 
% where ms= m-script, mf= m-function, mex= executable function (c-code)
%
% All input and output variables are assumed to be in the Matlab workspace.
% The load vectors p_g and y_s must have the same number of columns. 
%
% Input from global workspace:
%
%	uset      	uset.pv_*  partitioning vector
%				uset.*size size 
%                          *=n,m etc., see Builduset
%
% uset.pv_m and uset.pv_n index vectors are wrt the g- set
% uset.pv_s and uset.pv_f index vectors are wrt the n- set
% etc.

% Reduce p_g to p_f if not empty

if ~isempty(p_g)
	
	% Reduce p_g to p_n

	if ( uset.msize > 0 )
		p_n = p_g( uset.pv_n, :);
		p_m = p_g( uset.pv_m, :);
		p_n = p_n + G_mn' * p_m;
	else
		p_n  = p_g;
	end

	% Reduce p_n to p_f

	if ( uset.ssize > 0 )
		p_f = p_n( uset.pv_f, :);
	else
		p_f = p_n;
	end

% Create a p_f if p_g is empty

else
	
	if isempty(y_s)
		disp('No right hand side found...');
		p_f=[];
		return;
	else
		ncols= size( y_s, 2);
		p_f= sparse([],[],[],uset.fsize, ncols);
	end
	
end
	
% for nonzero boundary conditions, subtract coupling forces,
% this operation is only executed if flag_y_s > 0,
% the flag is set in the calling script for linear analysis
% and the first iteration step in nonlinear analysis

if ( flag_y_s > 0 ) & ( ~isempty(y_s) )
	% check if p_f and y_s have the same number of columns

	ncolp= size(p_f,2);
	ncoly= size(y_s,2);

	if ncolp ~= ncoly
		disp('Right hand side can not be calculated...');
		return;
	end

	p_f = p_f - K_fs * y_s;
		
end


end
