function [THIS,NPATH,SING1] = mysolve(THIS,ALT,OPT)
% mysolve  [Not a public function] First-order accurate quasi-triangular state-space form.
%
% Backend IRIS function.
% No help provided.

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

% If called with just one input arguments, `solve_` operates in a fast
% mode: works only on the first parameterisations and skips the
% expansion matrices.
fast = nargin == 1;

if ~exist('ALT','var')
    ALT = 1;
end

if ~exist('OPT','var')
    % Create default options when `mysolve` is called from `myupdatemodel`.
    % This is to avoid calling `passvalopt` which is slow, and hence not
    % a good idea for iterative procedures.
    OPT = struct( ...
        'progress',false, ...
        'select',true, ...
        'symbolic',true, ...
        'tolerance',getrealsmall(), ...
        'warning',false, ...
        'expand',0);
end

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

tol = OPT.tolerance;
ny = length(THIS.systemid{1});
nx = length(THIS.systemid{2});
nb = sum(imag(THIS.systemid{2}) < 0);
nf = nx - nb;
ne = length(THIS.systemid{3});
nn = sum(THIS.nonlin);
fkeep = ~THIS.metadelete;
nfkeep = sum(fkeep);
nalt = size(THIS.Assign,3);

if islogical(ALT)
    ALT = find(ALT);
elseif isequal(ALT,Inf)
    ALT = 1 : nalt;
end

% Pre-allocate solution matrices.
doallocsolution();

% Set `NPATH` to 1 initially to handle correctly the cases when only a
% subset of parameterisations is solved for.
NPATH = ones(1,nalt);
ALT = ALT(:).';

if OPT.progress
    progress = progressbar('IRIS model.solve progress');
end

SING1 = false(sum(THIS.eqtntype == 1),nalt);

for ialt = ALT
    % Select only the equations in which at least one parameter has
    % changed since last differentiation.
    eqselect = affectedeqtn(THIS,ialt,OPT.select);
    eqselect(THIS.eqtntype >= 3) = false;
    [THIS,deriv,nanderiv] = ...
        myderiv(THIS,eqselect,ialt,OPT.symbolic);
    if any(nanderiv)
        donanderiv();
    end
    [THIS,system] = mysystem(THIS,deriv,eqselect,ialt);
    % Check system matrices for complex numbers.
    if ~isreal(system.K{1}) ...
            || ~isreal(system.K{2}) ...
            || ~isreal(system.A{1}) ...
            || ~isreal(system.A{2}) ...
            || ~isreal(system.B{1}) ...
            || ~isreal(system.B{2}) ...
            || ~isreal(system.E{1}) ...
            || ~isreal(system.E{2})
        NPATH(ialt) = 1i;
        continue;
    end
    % Check system matrices for NaNs.
    if any(isnan(system.K{1})) ...
            || any(isnan(system.K{2})) ...
            || any(any(isnan(system.A{1}))) ...
            || any(any(isnan(system.A{2}))) ...
            || any(any(isnan(system.B{1}))) ...
            || any(any(isnan(system.B{2}))) ...
            || any(any(isnan(system.E{1}))) ...
            || any(any(isnan(system.E{2})))
        NPATH(ialt) = NaN;
        continue;
    end
    [SS,TT,QQ,ZZ,THIS.eigval(1,:,ialt),eqorder] = doschur();
    if NPATH(ialt) == 1
        if ~dosspace()
            NPATH(ialt) = -1;
        end
    end
    if OPT.progress
        update(progress,ialt/length(ALT));
    end
end

if OPT.expand > 0
    THIS = expand(THIS,OPT.expand);
end

% Nested functions.

%**********************************************************************
    function [SS,TT,QQ,ZZ,eigval,eqorder] = doschur()
        % Ordered real QZ decomposition.
        fA = full(system.A{2});
        fB = full(system.B{2});
        eqorder = 1 : size(fA,1);
        % If the QZ re-ordering fails, change the order of equations --
        % place the first equation last, and repeat.
        warning('off','MATLAB:ordqz:reorderingFailed');
        while true
            AA = fA(eqorder,:);
            BB = fB(eqorder,:);
            [SS,TT,QQ,ZZ] = qz(AA,BB,'real');
            % Ordered inverse eigvals.
            eigval = -ordeig(SS,TT);
            eigval = eigval(:).';
            issevn2 = dosevn2patch();
            stable = abs(eigval) >= 1 + tol;
            unit = abs(abs(eigval)-1) < tol;
            % Clusters of unit, stable, and unstable eigenvalues.
            clusters = zeros(size(eigval));
            % Unit roots first.
            clusters(unit) = 2;
            % Stable roots second.
            clusters(stable) = 1;
            % Unstable roots last.
            % Re-order by the clusters.
            lastwarn('');
            [SS,TT,QQ,ZZ] = ordqz(SS,TT,QQ,ZZ,clusters);
            isemptywarn = isempty(lastwarn());
            % If the first equations is ordered second, it indicates the
            % next cycle would bring the equations to their original order.
            % We stop and throw an error.
            if isemptywarn || eqorder(2) == 1
                break
            else
                eqorder = eqorder([2:end,1]);
            end
        end
        warning('on','MATLAB:ordqz:reorderingFailed');
        if ~isemptywarn
            utils.error('model', ...
                ['QZ re-ordering failed because ', ...
                'some eigenvalues are too close to swap, and ', ...
                'equation re-ordering does not help.']);
        end
        if OPT.warning && eqorder(1) ~= 1
            utils.warning('model', ...
                ['Numerical instability in QZ decomposition. ', ...
                'Equations re-ordered %g time(s).'], ...
                eqorder(1)-1);
        end
        
        % Re-order the inverse eigvals.
        eigval = -ordeig(SS,TT);
        eigval = eigval(:).';
        issevn2 = dosevn2patch() | issevn2;
        if OPT.warning && issevn2
            utils.warning('model', ...
                ['Numerical instability in QZ decomposition. ', ...
                'SEVN2 patch applied.'])
        end
        
        % Undo the eigval inversion.
        infeigval = eigval == 0;
        eigval(~infeigval) = 1./eigval(~infeigval);
        eigval(infeigval) = Inf;
        nunit = sum(unit);
        nstable = sum(stable);
        
        % Check BK saddle-path condition.
        if any(isnan(eigval))
            NPATH(ialt) = -2;
        elseif nb == nstable + nunit
            NPATH(ialt) = 1;
        elseif nb > nstable + nunit
            NPATH(ialt) = 0;
        else
            NPATH(ialt) = Inf;
        end

        %******************************************************************
        function flag = dosevn2patch()
            % Sum of two eigvals near to 2 may indicate inaccuracy.
            % Largest eigval less than 1.
            flag = false;
            eigval0 = eigval;
            eigval0(abs(eigval) >= 1-tol) = 0;
            eigval0(imag(eigval) ~= 0) = 0;
            if any(eigval0 ~= 0)
                [ans,below] = max(abs(eigval0)); %#ok<*NOANS,*ASGLU>
            else
                below = [];
            end
            % Smallest eig greater than 1.
            eigval0 = eigval;
            eigval0(abs(eigval) <= 1+tol) = Inf;
            eigval0(imag(eigval) ~= 0) = Inf;
            if any(~isinf(eigval0))
                [ans,above] = min(abs(eigval0));
            else
                above = [];
            end
            if ~isempty(below) && ~isempty(above) ...
                    && abs(eigval(below) + eigval(above) - 2) <= tol ...
                    && abs(eigval(below) - 1) <= 1e-6
                eigval(below) = sign(eigval(below));
                eigval(above) = sign(eigval(above));
                TT(below,below) = sign(TT(below,below))*abs(SS(below,below));
                TT(above,above) = sign(TT(above,above))*abs(SS(above,above));
                flag = true;
            end
        end
        % dosevn2patch().
        
    end
% doschur().

%**************************************************************************
    function flag = dosspace()
        
        flag = true;
        isnonlin = any(THIS.nonlin);
        S11 = SS(1:nb,1:nb);
        S12 = SS(1:nb,nb+1:end);
        S22 = SS(nb+1:end,nb+1:end);
        T11 = TT(1:nb,1:nb);
        T12 = TT(1:nb,nb+1:end);
        T22 = TT(nb+1:end,nb+1:end);
        Z11 = ZZ(fkeep,1:nb);
        Z12 = ZZ(fkeep,nb+1:end);
        Z21 = ZZ(nf+1:end,1:nb);
        Z22 = ZZ(nf+1:end,nb+1:end);
        
        % Transform the other system matrices by QQ.
        if eqorder(1) == 1
            % No equation re-ordering.
            % Constant.
            C = QQ*system.K{2};
            % Effect of transition shocks.
            D = QQ*full(system.E{2});
            if isnonlin
                % Effect of add-factors in transition equations earmarked
                % for non-linear simulations.
                N = QQ*system.N{2};
            end
        else
            % Equations have been re-ordered while computing QZ.
            % Constant.
            C = QQ*system.K{2}(eqorder,:);
            % Effect of transition shocks.
            D = QQ*full(system.E{2}(eqorder,:));
            if isnonlin
                % Effect of add-factors in transition equations earmarked
                % for non-linear simulations.
                N = QQ*system.N{2}(eqorder,:);
            end
        end

        C1 = C(1:nb,1);
        C2 = C(nb+1:end,1);
        D1 = D(1:nb,:);
        D2 = D(nb+1:end,:);
        if isnonlin
            N1 = N(1:nb,:);
            N2 = N(nb+1:end,:);
        end
        
        % Quasi-triangular state-space form.
        
        U = Z21;
        
        % Steady state for non-linear models. They are needed in non-linear
        % models to back out the constant vectors.
        if ~THIS.linear
            ysstate = ...
                mytrendarray(THIS,THIS.solutionid{1},0,false,ialt);
            xfsstate = ...
                mytrendarray(THIS,THIS.solutionid{2}(1:nfkeep),[-1,0],false,ialt);
            xbsstate = ...
                mytrendarray(THIS,THIS.solutionid{2}(nfkeep+1:end),[-1,0],false,ialt);
            asstate = U \ xbsstate;
            if any(isnan(asstate(:)))
                flag = false;
                return
            end
        end
        
        % Unstable block.
        
        G = -Z21\Z22;
        if any(isnan(G(:)))
            flag = false;
            return
        end
        
        Ru = -T22\D2;
        if any(isnan(Ru(:)))
            flag = false;
            return
        end
        
        if isnonlin
            Yu = -T22\N2;
            if any(isnan(Yu(:)))
                flag = false;
                return
            end
        end
        
        if THIS.linear
            Ku = -(S22+T22)\C2;
        else
            Ku = zeros(nfkeep,1);
        end
        if any(isnan(Ku(:)))
            flag = false;
            return
        end

        % Transform stable block == transform backward-looking variables:
        % a(t) = s(t) + G u(t+1).
        
        Ta = -S11\T11;
        if any(isnan(Ta(:)))
            flag = false;
            return
        end
        Xa0 = S11\(T11*G + T12);
        if any(isnan(Xa0(:)))
            flag = false;
            return
        end
        
        Ra = -Xa0*Ru - S11\D1;
        if any(isnan(Ra(:)))
            flag = false;
            return
        end
        
        if isnonlin
            Ya = -Xa0*Yu - S11\N1;
            if any(isnan(Ya(:)))
                flag = false;
                return
            end
        end
        
        Xa1 = G + S11\S12;
        if any(isnan(Xa1(:)))
            flag = false;
            return
        end
        if THIS.linear
            Ka = -(Xa0 + Xa1)*Ku - S11\C1;
        else
            Ka = asstate(:,2) - Ta*asstate(:,1);
        end
        if any(isnan(Ka(:)))
            flag = false;
            return
        end
        
        % Forward-looking variables.
        
        % Duplicit rows (metadelete) already deleted from Z11 and Z12.
        Tf = Z11;
        Xf = Z11*G + Z12;
        Rf = Xf*Ru;
        if isnonlin
            Yf = Xf*Yu;
        end
        if THIS.linear
            Kf = Xf*Ku;
        else
            Kf = xfsstate(:,2) - Tf*asstate(:,1);
        end
        if any(isnan(Kf(:)))
            flag = false;
            return
        end
        
        % State-space form:
        % [xf(t);a(t)] = T a(t-1) + K + R(L) e(t) + Y(L) addfactor(t),
        % U a(t) = xb(t).
        T = [Tf;Ta];
        K = [Kf;Ka];
        R = [Rf;Ra];
        if isnonlin
            Y = [Yf;Ya];
        end
        
        % Remove negligible entries from U.
        %maxsing = svds(U,1);
        %tolU = size(U,1)*eps(maxsing)^(2/3);
        %U(abs(U) <= tolU) = 0;

        % y(t) = Z a(t) + D + H e(t)
        if ny > 0
            ZZ = -full(system.A{1}\system.B{1});
            if any(isnan(ZZ(:)))
                flag = false;
                s = size(system.A{1},1);
                r = rank(full(system.A{1}));
                if r < s
                    d = s - r;
                    [U,S] = svd(full(system.A{1})); %#ok<NASGU>
                    SING1(:,ialt) = any(abs(U(:,end-d+1:end)) > tol,2);
                end
                return
            end
            H = -full(system.A{1}\system.E{1});
            if any(isnan(H(:)))
                flag = false;
                return
            end
            if THIS.linear
                D = full(-system.A{1}\system.K{1});
            else
                D = ysstate - ZZ*xbsstate(:,2);
            end
            if any(isnan(D(:)))
                flag = false;
                return
            end
            Z = ZZ*U;
            % Remove negligible entries from Z.
            %maxsing = svds(Z,1);
            %tolZ = size(Z,2)*eps(maxsing)^(2/3);
            %Z(abs(Z) <= tolZ) = 0;
        else
            Z = zeros(0,nb);
            H = zeros(0,ne);
            D = zeros(0,1);
        end
        
        % Necessary initial conditions in xb vector.
        if ~fast
            THIS.icondix(1,:,ialt) = any(abs(T/U) > tol,1);
        end

        if ~fast && ~isnan(OPT.expand)
            % Forward expansion.
            % a(t) <- -Xa J^(k-1) Ru e(t+k)
            % xf(t) <- Xf J^k Ru e(t+k)
            J = -T22\S22;
            Xa = Xa1 + Xa0*J;
            % Highest computed power of J: e(t+k) requires J^k.
            Jk = eye(size(J));

            THIS.Expand{1}(:,:,ialt) = Xa;
            THIS.Expand{2}(:,:,ialt) = Xf;
            THIS.Expand{3}(:,:,ialt) = Ru;
            THIS.Expand{4}(:,:,ialt) = J;
            THIS.Expand{5}(:,:,ialt) = Jk;
            if isnonlin
                THIS.Expand{6}(:,:,ialt) = Yu;
            end
        end
        
        THIS.solution{1}(:,:,ialt) = T;
        THIS.solution{2}(:,:,ialt) = R;
        THIS.solution{3}(:,:,ialt) = K;
        THIS.solution{4}(:,:,ialt) = Z;
        THIS.solution{5}(:,:,ialt) = H;
        THIS.solution{6}(:,:,ialt) = D;
        THIS.solution{7}(:,:,ialt) = U;
        if isnonlin
            THIS.solution{8}(:,:,ialt) = Y;
        end

    end
% dosspace().

%**************************************************************************
    function doallocsolution()
        if isempty(THIS.eigval)
            THIS.eigval = nan([1,nx,nalt]);
        else
            THIS.eigval(:,:,ALT) = NaN;
        end
        
        if isempty(THIS.icondix)
            THIS.icondix = false(1,nb,nalt);
        else
            THIS.icondix(1,:,ALT) = false;
        end
        
        if isnan(OPT.expand)
            THIS.Expand = {};
        else
            if isempty(THIS.Expand) || isempty(THIS.Expand{1})
                THIS.Expand{1} = nan(nb,nf,nalt);
                THIS.Expand{2} = nan(nfkeep,nf,nalt);
                THIS.Expand{3} = nan(nf,ne,nalt);
                THIS.Expand{4} = nan(nf,nf,nalt);
                THIS.Expand{5} = nan(nf,nf,nalt);
                THIS.Expand{6} = nan(nf,nn,nalt);
            else
                THIS.Expand{1}(:,:,ALT) = NaN;
                THIS.Expand{2}(:,:,ALT) = NaN;
                THIS.Expand{3}(:,:,ALT) = NaN;
                THIS.Expand{4}(:,:,ALT) = NaN;
                THIS.Expand{5}(:,:,ALT) = NaN;
                THIS.Expand{6}(:,:,ALT) = NaN;
            end
        end
        
        if isempty(THIS.solution) || isempty(THIS.solution{1})
            THIS.solution{1} = nan(nfkeep+nb,nb,nalt); % T
            THIS.solution{2} = nan(nfkeep+nb,ne,nalt); % R
            THIS.solution{3} = nan(nfkeep+nb,1,nalt); % K
            THIS.solution{4} = nan(ny,nb,nalt); % Z
            THIS.solution{5} = nan(ny,ne,nalt); % H
            THIS.solution{6} = nan(ny,1,nalt); % D
            THIS.solution{7} = nan(nb,nb,nalt); % U
            THIS.solution{8} = nan(nfkeep+nb,nn,nalt); % Y
        else
            THIS.solution{1}(:,:,ALT) = nan(nfkeep+nb,nb,length(ALT));
            if size(THIS.solution{2},2) > ne
                THIS.solution{2} = nan(nfkeep+nb,ne,nalt);
            else
                THIS.solution{2}(:,:,ALT) = nan(nfkeep+nb,ne,length(ALT));
            end
            THIS.solution{3}(:,:,ALT) = nan(nfkeep+nb,1,length(ALT));
            THIS.solution{4}(:,:,ALT) = nan(ny,nb,length(ALT));
            THIS.solution{5}(:,:,ALT) = nan(ny,ne,length(ALT));
            THIS.solution{6}(:,:,ALT) = nan(ny,1,length(ALT));
            THIS.solution{7}(:,:,ALT) = nan(nb,nb,length(ALT));
            if size(THIS.solution{8},2) > nn
                THIS.solution{8} = nan(nfkeep+nb,nn,nalt);
            else
                THIS.solution{8}(:,:,ALT) = nan(nfkeep+nb,nn,length(ALT));
            end
        end
    end
% doallocsolution().

%**************************************************************************
    function donanderiv()
        % donanderiv  Report equations with NaN derivatives.
        utils.error('model', ...
            'NaN derivatives in this equation: ''%s''.', ...
            THIS.eqtn{nanderiv});
    end

end