function u_g=icediagnostic_core_linear(m,analysis_type,varargin);
%ICEDIAGNOSTIC_CORE_LINEAR - core solution of linear problems
%
%   Core of the vertical diagnostic solution. No need for a non-linear iteration.
%   Otherwise similar to ICEDIAGNOSTIC_CORE_NONLINEAR
%
%   Usage:
%      u_g=icediagnostic_core_linear(m,analysis_type,varargin);
%
%   See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC2D, ICEDIAGNOSTIC3D, ICEDIAGNOSTIC_CORE_NONLINEAR

%global variables
global cluster gridset

if nargin==3,
	inputs=varargin{1};
else
	inputs=struct();
end

%recover fem model fields
elements=m.elements;
grids=m.grids;
materials=m.materials;
loads=m.loads;
ys=m.ys;
gridset=m.gridset;
G_mn=m.Gmn;
params=m.params;


%recover parameters
sparsity=params.sparsity;
solver_type=params.solver_type;
debug=params.debug;

disp(sprintf('%s','   solving'));

% Generate system matrices (stiffness and load)
kflag=1; pflag=1;

%generate stiffness and loads
[K_gg,p_g]=SystemMatrices(elements,grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type);
[K_gg,p_g]=PenaltySystemMatrices(grids,loads,materials,kflag, pflag, sparsity,inputs,analysis_type,K_gg,p_g);

if cluster, 
	K_gg=distributed(gplus(K_gg),'convert');
	p_g=gplus(p_g);
end

% Reduce stiffness matrix from g set to f set (effectively rooting out the single point constraints), 
% and compute modifications to loads from single point constraints.
[K_ff K_fs]=IceReducematrixfromgtof(K_gg,G_mn); 

% Reduce load from g set to f set
p_f=IceReducerightside(p_g,G_mn,K_fs,ys);

% Solve
uf=Solver(K_ff,p_f,solver_type);
if debug,
	disp(sprintf('%s%g','      condition number of stiffness matrix: ',condest(K_ff)));
end
   

% Add single point constraints back, ie increase f-set by s-set into the global g-set.
u_g= full(Mergesolution_g(uf,G_mn,ys));
