function  [m,flag] = mysstatelinear(m,opt)
% MYSSTATELINEAR  [Not a public function] Steady-state solver for linear models.
%
% Backed IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

if ~exist('opt','var')
    opt = struct();
    opt.warning = true;
end

%**************************************************************************

ny = sum(m.nametype == 1);
realsmall = getrealsmall();
nalt = size(m.Assign,3);
flag = true([1,nalt]);

[flag,index] = isnan(m,'solution');
if opt.warning && flag
    % Solutions for some parameterisations don't exist.
    warning_(55,sprintf(' #%g',find(index)));
end

for ialt = find(~index)
    doonesstate();
end

if any(m.log)
    realassign = real(m.Assign(1,m.log,:));
    imagassign = imag(m.Assign(1,m.log,:));
    m.Assign(1,m.log,:) = exp(realassign) + 1i*exp(imagassign);
end

% Nested functions.

%***********************************************************************
    function doonesstate()
        T = m.solution{1}(:,:,ialt);
        K = m.solution{3}(:,:,ialt);
        Z = m.solution{4}(:,:,ialt);
        D = m.solution{6}(:,:,ialt);
        U = m.solution{7}(:,:,ialt);
        [nx,nb] = size(T);
        nunit = sum(abs(abs(m.eigval(1,:,ialt))-1) <= realsmall);
        nf = nx - nb;
        nstable = nb - nunit;
        Tf = T(1:nf,:);
        Ta = T(nf+1:end,:);
        Kf = K(1:nf,1);
        Ka = K(nf+1:end,1);
        % Alpha vector.
        if any(any(abs(Ta(1:nunit,1:nunit) - eye(nunit)) > realsmall))
            % I(2) or higher-order systems.
            % Write the steady-state system at two different times: t and t+d.
            d = 10;
            E1 = [eye(nb),zeros(nb);eye(nb),d*eye(nb)];
            E2 = [Ta,-Ta;Ta,(d-1)*Ta];
            temp = pinv(E1 - E2) * [Ka;Ka];
            a2 = temp(nunit+(1:nstable));
            da1 = temp(nb+(1:nunit));
            utils.warning('model', ...
                ['Model is not stationary or difference stationary. ', ...
                'Some of the steady-state growth rates are not pinned down ', ...
                'to fixed numbers.']);
        else
            % I(0) or I(1) systems.
            a2 = (eye(nstable) - Ta(nunit+1:end,nunit+1:end)) ...
                \ Ka(nunit+1:end,1);
            da1 = Ta(1:nunit,nunit+1:end)*a2 + Ka(1:nunit,1);
        end
        % Transition variables.
        x = [Tf*[-da1;a2]+Kf;U(:,nunit+1:end)*a2];
        dx = [Tf(:,1:nunit)*da1;U(:,1:nunit)*da1];
        x(abs(x) <= realsmall) = 0;
        dx(abs(dx) <= realsmall) = 0;
        realid = real(m.solutionid{2});
        imagid = imag(m.solutionid{2});
        index = imagid == 0;
        realid(~index) = [];
        x(~index) = [];
        dx(~index) = [];
        m.Assign(1,realid,ialt) = x(:).' + 1i*dx(:).';
        % Measurement variables.
        if ny > 0
            realid = real(m.solutionid{1});
            y = Z(:,nunit+1:end)*a2 + D;
            dy = Z(:,1:nunit)*da1;
            m.Assign(1,realid,ialt) = y(:).' + 1i*dy(:).';
        end
    end
% doonesstate().

end