function [x,U] = IDRstabg3(A,b,x0,M1,M2,tol,maxit,s,Utol2)

    % The IDRstab algorithm, robust GMRES-style version
    % [x] = IDRstabg(A,b) performs 400 iterations of IDR(4)stab(L) with L
    %  adaptively either 1 or 2 to solve A*x=b for x up to a relative
    %  residual of 1e-8 for an initial guess of x=0*b. 
    %  No preconditioning is used.
    %
    % [x] = IDRstabg(A,b,x0,M1,M2,tol,maxit,s) uses IDR(s)stab(L), with L
    %  as above, an initial guess x0, left preconditioner M1, right 
    %  preconditioner M2, relative tolerance tol, terminates after at 
    %  most maxit iterations.
    %
    % Author:   Martin Neuenhofen
    % Version:  26.10.2016
    % 
    % Copyrigtht (c) 2016 by Martin Neuenhofen.
    % This program is distributed without any warranty.
    
    % init user inputs
    if( nargin < 1 || isempty(A) || ~ismatrix(A) || size(A,1)~=size(A,2) )
        fprintf('First input must be a quadratic system matrix A.\n');
        x           = NaN;
        return;
    end
    if( nargin < 2 || isempty(b) )
        fprintf('Second input must be a n-by-1 right-hand side b.\n');
        x           = NaN;
        return;
    else
        N           = size(b,1);
    end
    if(nargin < 3 || isempty(x0) || size(x0,1)~=N)
        if(nargin>2 && ~isempty(x0) && size(x0,1)~=N)
            fprintf('Warning: Initial guess x0 (third input) does not match the dimension N=%i of b.\n',N);
            fprintf('Notice: Initial guess x0 is set to zero vector.\n');
        end
        x0          = zeros(N,1);
    end
    if(nargin < 4 || isempty(M1) || ~ismatrix(M1) || size(M1,1)~=size(M1,2) )
        M1          = speye(N,N);
    end
    if(nargin < 5 || isempty(M2) || ~ismatrix(M2) || size(M2,1)~=size(M2,2) )
        M2          = speye(N,N);
    end
    if(nargin < 6 || isempty(tol))
        tol         = 1e-8;
    end
    if(tol<eps || tol>1)
        if(tol < 1e3*eps)
            fprintf('Warning: Relative residual tolerance tol=%1.3g (sixth input)\nis assumed to not be achievable.\n',tol);
        elseif(tol < eps)
            fprintf('Warning: Relative residual tolerance tol=%1.3g (sixth input)\nis most probably not achievable.\n',tol);
        elseif(tol > 1)
            fprintf('Hint: Relative residual tolerance tol=%1.3g (sixth input) should \nlive in the open interval (eps,1).\n',tol);
        end
    end
    if(nargin < 7 || isempty(maxit))
        maxit       = 400;
    end
    if(maxit < 100)
        if(maxit < 0)
            fprintf('Warning: Maximum number of iterations maxit=%i (seventh input)\nmust be larger then 0.\n',tol);
        elseif(s<maxit)
            fprintf('Hint: Below 100 iterations one should not use an iterative solver.\nTry GMRES instead by choosing s=%i.\n',maxit);
        end
    end
    if(nargin < 8 || isempty(s))
        s           = 4;
    elseif(s<=0)
        fprintf('Warning: Parameter s (eighth input) must be greater or equal to 1.\n');
        fprintf('Notice: s is set to 4.\n');
        s           = 4;
    end
    if(nargin < 9 || isempty(Utol2))
        recycle_U       = false;
        flag_written_U  = true;   % beta < tol2 can never be achieved --> no write-back
        tol2            = 0;
    elseif(size(Utol2,1)==N && size(Utol2,2)==s)
        recycle_U       = true;
        flag_written_U  = true;
        tol2            = 0;
    elseif(size(Utol2,1)==1 && size(Utol2,2)==1)
        recycle_U       = false;
        tol2    = Utol2;
        if(tol2 > 1e0)
            tol2            = 1e-2;
        elseif(tol2 < 1e-4)
            tol2            = 1e-4;
        end
        flag_written_U  = false;
        fprintf('Notice: Write out data for M(%i)stab(L) at tol2=%1.3g\n',s,tol2);
    end
    U               = NaN;
    
    % init r and x
    tic;
    x           = M2*x0;
    r           = M1\(b - A*x0);
    b           = M1\b;
    norm_b      = norm(b);
    beta        = norm(r);
    
    abstol      = tol  * norm_b;
    abstol2     = tol2 * norm_b;
    if(beta < abstol)
        fprintf('Initial guess satisfies residual requirements:\n\tnorm(b-A*x0)/norm(b)=%1.3g < %1.3g\n',beta/norm_b,tol);
        return;
    else
        fprintf('#matvec\texec. time\t|r|\n');
        fprintf('\t%i\t%3.3f\t%1.3g\n',0,0,norm_b);
        fprintf('\t%i\t%3.3f\t%1.3g\n',1,toc,beta);
    end
    
    iter        = 1;
    ITER        = [0,1];
    RESNORM     = [norm_b,beta];
    
    % Initialize the direct sub lin solver
    opts_triu   = []; opts_triu.UT = true;
    opts_tril   = []; opts_tril.LT = true;
    
    % Initialize all large arrays
%     V0          = zeros(N,s);
%     V1          = zeros(N,s);
    W           = zeros(N,2*s+3);
    H           = zeros(2*s+3,2*s+2);
    tH          = zeros(2*s+3,2*s+2);
    Y           = zeros(s,2*s+2);
    ce        	= zeros(2*s+2,1);
    se         	= zeros(2*s+2,1);
    gamma       = zeros(2*s+3,1);
    
    if(recycle_U)
        rng('default');
        P           = orth(randn(N,s));     % initialize P in the same way as before
        V0          = Utol2;
        V1          = M1\(A*(M2\V0));
        [V1,R]      = qr(V1,0);
        V0          = V0/R;
        Z           = P'*V1;
        [Qz,Z]      = qr(Z');
        Z           = Z';
        V0          = V0 * Qz;
        V1          = V1 * Qz;
        % bio res
        zeta        = linsolve(Z,P'*r,opts_tril);
        x           = x + V0 * zeta;
        r           = r - V1 * zeta;
        beta        = norm(r);
        
        iter        = iter+s;
        ITER        = [ITER,iter];
        RESNORM     = [RESNORM,beta];
        
    else
        
        % GMRES for the first s iters
        W(:,1)      = 1/beta * r;
        gamma(1)    = beta;
        for j=1:s
            w           = M1\(A * (M2\W(:,j)));
            for i=1:j
                H(i,j)      = W(:,i)'*w;
                w           = w - H(i,j) * W(:,i);
            end
            H(j+1,j)    = norm(w);
            tH(1:j+1,j) = H(1:j+1,j);
            for i=1:j-1
                H(i:i+1,j)  = [ce(i+1) , se(i+1) ; -se(i+1) , ce(i+1)] * H(i:i+1,j);
            end
            geta        = sqrt( H(j,j)*H(j,j) + H(j+1,j)*H(j+1,j) );
            se(j+1)   	= H(j+1,j)/geta;
            ce(j+1)   	= H(j  ,j)/geta;
            H(j,j)      = geta;
            gamma(j+1)  = -se(j+1)*gamma(j);
            gamma(j)    =  ce(j+1)*gamma(j);

            iter        = iter+1;
            ITER        = [ITER,iter];
            RESNORM     = [RESNORM,abs(gamma(j+1))];
            fprintf('\t%i\t%3.3f\t%1.3g\n',iter,toc,abs(gamma(j+1)));
            if( abs(gamma(j+1))<=abstol || iter>=maxit)
                y           = linsolve(triu(H(1:j,1:j)),gamma(1:j),opts_triu);
                x           = x + W(:,1:j)*y;
                x           = M2 * x;
                interior_plotter_function(ITER,RESNORM);
                return;
            end
            W(:,j+1)    = 1/H(j+1,j) * w;
        end

        % Initialize IDR projector
        %   Understanding this is not important for the following.
        %   Description:
        %       From above GMRES we have
        %           A*W(:,1:s) = W(:,1:s+1) * Qh * triu(H)
        %       where Qh is (s+1)-by-(s) and only given from the above Givens
        %       rotations.
        %
        %       In addition we have Z=P'*W*Qh at line {eqn:Z_init}. Z is 
        %       overwritten such that Z*Q' := Z and Z becomes lower triangular.
        %
        %       We compute as follows:
        %               V0  = W(:,1:s  ) * inv( triu(H) ) * Q   % \label{eqn:V0_init}
        %               V1  = W(:,1:s+1) *      Qh        * Q   % \label{eqn:V1_init}
        %
        %   So afterwards it holds:
        %       P'* V1  = Z     is lower triangular
        %       A * V0  = V1    is orthogonal
        % + Compute P
        rng('default');
        P           = orth(randn(N,s));
        % + Compute Z
        Y(:,1:s+1)  = P'*W(:,1:s+1);
        % ++ Compute linear combination for oblique projected initial residual
        e           = beta * Y(:,1);
        y           = (Y(:,1:s+1)*tH(1:s+1,1:s))\e;
        x           = x + W(:,1:s)*y;
        e           = [beta;zeros(s,1)] - tH(1:s+1,1:s)*y;
        beta        = norm(e);
        r           = W(:,1:s+1) * e;
        ITER        = [ITER,iter];      % residual for identical iteration as before 
        RESNORM     = [RESNORM,beta];   % but residual grows (saw profile)
        fprintf('\t%i\t%3.3f\t%1.3g\n',iter,toc,beta);
        % ++ Apply Qh from the right to Z
        for i=1:s
            Y(:,i:i+1)  = Y(:,i:i+1) * [ce(i+1) , se(i+1) ; -se(i+1) , ce(i+1)]';
        end
        Z           = Y(:,1:s);     % \label{eqn:Z_init}
        % ++ Make the LQ' decomposition for Z
        [Q,Z]       = qr(Z');
        Z           = Z';           % Z is lower left triangular
        H(1:s,1:s)  = linsolve(triu(H(1:s,1:s)),Q,opts_triu);
        Q           = [Q;zeros(1,s)];
        for i=s:-1:1
            Q(i:i+1,:) 	= [ce(i+1) , se(i+1) ; -se(i+1) , ce(i+1)]' * Q(i:i+1,:);
        end
        % + Compute the right linear factor matrices from the equations 
        %   {eqn:V0_init}-{eqn:V1_init}
        V0          = W(:,1:s  ) * H(1:s,1:s);
        V1          = W(:,1:s+1) * Q;

        % Validation of the initialization: 
        % The following should be < C * eps_machine * (cond(A)+cond(Z))
    %     cons_Z  = norm(P'*V1-Z)/norm(Z)
    %     cons_V1 = norm(V1-A*V0,'inf')/norm(V1,'inf')  
    %     bior_r  = norm(P'*r)/norm(r)
    end
    
    % initialize parameters for alternative cycles of L
    x_global        = zeros(N,1);           % it holds: b - A*(x_global+x) = r
    b_local         = b;
    norm_b_local    = norm(b_local);
    repl_up_gap     = 1e-2;
    repl_down_gap   = 1e-2;
    interm_res_max  = beta;
    
    cycle2_count    = 0;
    cycle2_max      = 3;            % maximum number of subsequent L=2-cycles before L=1-cycle will be performed
    
    % Start the iteration loop
    while(true)
        
        % plan the flying restart & L=1
        if(beta < repl_down_gap * norm_b_local)
            replace_down    = true;
            replace_up      = false;
        else
            replace_down    = false;
            if(beta < repl_up_gap * interm_res_max)
                replace_up      = true;
            else
                replace_up      = false;
                interm_res_max  = max( interm_res_max,beta );
            end
        end
        
        % induce the cycle type
        if( replace_up || replace_down || cycle2_count>=cycle2_max)
            CYCLE_TYPE      = 'IDR(s)stab(1)';
            cycle2_count    = 0;
        else
            CYCLE_TYPE      = 'IDR(s)stab(2)';
            cycle2_count    = cycle2_count+1;
        end
        
%         CYCLE_TYPE      = 'IDR(s)stab(1)';
        % Perform the cycle
        switch CYCLE_TYPE
            
            case 'IDR(s)stab(1)'
                r1          = M1\(A * (M2\r));
                % choose omega
                kappa       = 0.7;
                r1Tr        = r1'*r;
                norm_r1     = norm(r1);
                omg         = r1Tr / ( norm_r1*norm_r1 );
                rho         = r1Tr / ( norm_r1 * beta );
                if(abs(rho)<kappa)
                    omg         = omg * kappa / abs(rho);
                end
                x           = x + omg * r;
                r           = r - omg * r1;
                beta        = sqrt( beta*beta - 2 * omg * real(r1Tr) + (omg*norm_r1)^2 );
                iter        = iter+1;
                ITER        = [ITER,iter];
                RESNORM     = [RESNORM,abs(gamma(j+1))];
                fprintf('\t%i\t%3.3f\t%1.3g\n',iter,toc,abs(gamma(j+1)));
                if( abs(gamma(j+1))<=abstol || iter>=maxit)
                    x           = x + x_global;
                    x           = M2 \ x;
                    interior_plotter_function(ITER,RESNORM);
                    return;
                end
                        
                % Projector modification
                % + here is the expensive step
                %   we do [Q,R] = qr([V1 V0],0) by appending the V0 on the
                %   QR of V1 (notice that V1 is already orthogonal).
                R           = zeros(2*s,s);
%                 V1t     = V1;     % temporary for later validation
                for i=1:s
                    for j=1:s
                        R(j,i)      = V1(:,j)'*V0(:,i);
                        V0(:,i)     = V0(:,i) - R(j,i)   * V1(:,j);
                    end
                    for j=1:i-1
                        R(s+j,i)    = V0(:,j)'*V0(:,i);
                        V0(:,i)     = V0(:,i) - R(s+j,i) * V0(:,j);
                    end
                    R(s+i,i)   	= norm( V0(:,i) );
                    V0(:,i)     = 1/R(s+i,i) * V0(:,i);
                end
                % + the relation [V1 V0] = Q*R (where V1 and V0 are
                %   already overwritten) can be utilized to do the update
                %       " V0 := -1/omg*V0 + V1 "
                %   in a robust way.
                F           = -1/omg * R;
                F(1:s,1:s)  = F(1:s,1:s) + eye(s,s);
                [Qf,Rf]     = qr(F,0);
                % reflip the Z to lower triangle from right
                Z           = Z/Rf;
                [Qz,Z]      = qr(Z');
                Z           = Z';
                % Do the above update for V0 in the robust way.
                V0          = [V1 V0] * (Qf*Qz);
%               errZ = norm( Z - ((P'*V1t)/Rf)*Qz ) / norm(Z)
                
                % GMRES for s steps
                %   This is:
                %       A * (W-V*Z\Y) = bar{W}*\underline{tH}
                %
                %       In the above decomposition it is Y = P'*W.
                %       The matrix Y is different for L=1 and L=2 !!
                W(:,1)      = 1/beta * r;
                gamma       = zeros(s+1,1);
                gamma(1)    = beta;
                for j=1:s
                    Y(:,j)      = P'*W(:,j);
                    zeta        = linsolve(Z,Y(:,j),opts_tril);
                    w           = W(:,j) - V0 * zeta;
                    w           = M1\(A * (M2\w));
                    for i=1:j
                        H(i,j)      = W(:,i)'*w;
                        w           = w - H(i,j) * W(:,i);
                    end
                    H(j+1,j)    = norm(w);
                    tH(1:j+1,j) = H(1:j+1,j);
                    for i=1:j-1
                        H(i:i+1,j)  = [ce(i+1) , se(i+1) ; -se(i+1) , ce(i+1)] * H(i:i+1,j);
                    end
                    geta        = sqrt( H(j,j)*H(j,j) + H(j+1,j)*H(j+1,j) );
                    se(j+1)   	= H(j+1,j)/geta;
                    ce(j+1)   	= H(j  ,j)/geta;
                    H(j,j)      = geta;
                    gamma(j+1)  = -se(j+1)*gamma(j);
                    gamma(j)    =  ce(j+1)*gamma(j);

                    iter        = iter+1;
                    ITER        = [ITER,iter];
                    RESNORM     = [RESNORM,abs(gamma(j+1))];
                    fprintf('\t%i\t%3.3f\t%1.3g\n',iter,toc,abs(gamma(j+1)));
                    if( abs(gamma(j+1))<=abstol || iter>=maxit)
                        y           = linsolve(triu(H(1:j,1:j)),gamma(1:j),opts_triu);
                        x           = x + W(:,1:j)*y - V0 * linsolve(Z,Y(:,1:j)*y,opts_tril);
                        x           = x + x_global;
                        x           = M2 \ x;
                        interior_plotter_function(ITER,RESNORM);
                        return;
                    end
                    W(:,j+1)    = 1/H(j+1,j) * w;
                end
                Y(:,s+1)    = P'*W(:,s+1);
                
                % bi-orthogonalize the residual
                y           = (Y(:,1:s+1)*tH(1:s+1,1:s)) \ (beta*Y(:,1));
                x           = x + W(:,1:s) * y - V0 * linsolve(Z,Y(:,1:s)*y,opts_tril);
                e           = [beta;zeros(s,1)] - tH(1:s+1,1:s)*y;
                beta        = norm(e);
                ITER        = [ITER,iter];      % again: same iteration
                RESNORM     = [RESNORM,beta];
                fprintf('\t%i\t%3.3f\t%1.3g\n',iter,toc,beta);
                r           = W(:,1:s+1) * e;
%                 bioRes      = norm(P'*r)/beta
                
                % build new projector
                [Q,R]       = qr(tH(1:s+1,1:s),0);
                hZ          = Y(:,1:s+1)*Q;
                [Qz,hZ]     = qr(hZ');
                V0          = W(:,1:s) * (R\Qz) - V0 * linsolve(Z,Y(:,1:s)*(R\Qz),opts_tril);
                V1          = W(:,1:s+1)*(Q*Qz);
                Z           = hZ';
%                 err_Z       = norm(Z-P'*V1)/norm(Z)
                
                % Validation:
%                 cond_V1     = cond(V1)-1
%                 cons_V1     = norm(V1-A*V0)
%                 cons_r      = norm(b-A*x-r)/beta
%                 cons_beta   = (beta-norm(r))/beta
                
            %end case
            
            case 'IDR(s)stab(2)'
                % GMRES for 2*s+2 steps
                %   This is:
                %       A * (W-V*Z\Y) = bar{W}*\underline{tH}
                %       tH is again decomposed into Qg * triu(H)
                %       where Qg = G{2*s+2} * ... * G{1} is the bunch of Givens
                %       rotations.
                %
                %       In the above decomposition it is Y = P'*A*W
%                 beta
%                 norm(r)
                W(:,1)      = 1/beta * r;
                gamma       = zeros(2*s+3,1);
                gamma(1)    = beta;
                for j=1:2*s+2
                    w           = M1\(A * (M2\W(:,j)));
                    Y(:,j)      = P'* w;
                    zeta        = linsolve(Z,Y(:,j),opts_tril);
                    w           = w - V1 * zeta;
                    for i=1:j
                        H(i,j)      = W(:,i)'*w;
                        w           = w - H(i,j) * W(:,i);
                    end
                    H(j+1,j)    = norm(w);
                    tH(1:j+1,j) = H(1:j+1,j);
                    for i=1:j-1
                        H(i:i+1,j)  = [ce(i+1) , se(i+1) ; -se(i+1) , ce(i+1)] * H(i:i+1,j);
                    end
                    geta        = sqrt( H(j,j)*H(j,j) + H(j+1,j)*H(j+1,j) );
                    se(j+1)   	= H(j+1,j)/geta;
                    ce(j+1)   	= H(j  ,j)/geta;
                    H(j,j)      = geta;
                    gamma(j+1)  = -se(j+1)*gamma(j);
                    gamma(j)    =  ce(j+1)*gamma(j);

                    iter        = iter+1;
                    ITER        = [ITER,iter];
                    RESNORM     = [RESNORM,abs(gamma(j+1))];
                    fprintf('\t%i\t%3.3f\t%1.3g\n',iter,toc,abs(gamma(j+1)));
                    if( abs(gamma(j+1))<=abstol || iter>=maxit )
                        y           = linsolve(triu(H(1:j,1:j)),gamma(1:j),opts_triu);
                        x           = x + W(:,1:j)*y - V0 * linsolve(Z,Y(:,1:j)*y,opts_tril);
                        x           = x + x_global;
                        x           = M2 \ x;
                        interior_plotter_function(ITER,RESNORM);
                        return;
                    end
                    W(:,j+1)    = 1/H(j+1,j) * w;
                end
%                 cond_W = cond(W)
%                 err_arnoldi = norm(A*W(:,1:2*s+2) - [V1 W] * [Z\Y;H])/norm(A*W(:,1:2*s+2))
                
                % Compute biorthogonal residuals
                [Qy1,~]     = qr( Y(:,1:2*s+1)',0 );
                [Qy2,~]     = qr( Y(:,1:2*s+2)',0 );
                y           = ( [Qy1';Qy2'*tH(1:2*s+2,1:2*s+1)]*tH(1:2*s+1,1:2*s) ) \ [ Qy1'*[beta;zeros(2*s,1)] ; Qy2' * tH(1:2*s+2,1:2*s+1)*[beta;zeros(2*s,1)] ];
%                 res1        = Y(:,1:2*s+1)                       * ( [beta;zeros(2*s,1)]-tH(1:2*s+1,1:2*s)*y )
%                 res2        = Y(:,1:2*s+2) * tH(1:2*s+2,1:2*s+1) * ( [beta;zeros(2*s,1)]-tH(1:2*s+1,1:2*s)*y )
                % + update x according to bi-orthogonal residual r2 (unique under the given degrees of freedom)
                x           = x + W(:,1:2*s)*y - V0 * linsolve(Z,Y(:,1:2*s)*y,opts_tril);
%                 rel_orth_r = norm(P'*(b-A*x))/norm(b-A*x)
%                 norm_r = norm(b-A*x)

                % Compute stability polynomial coefficients tau and apply the polynomial to x and r
                % + compute linear factors for r r1 r2
                e           = [beta;zeros(2*s,1)];
                c0          = e - tH(1:2*s+1,1:2*s) * y;
                c1          = tH(1:2*s+2,1:2*s+1) * c0;
                c2          = tH(1:2*s+3,1:2*s+2) * c1;
                % + compute tau
                S           = [c0'*c0, 0, 0; c1(1:2*s+1)'*c0, c1'*c1, 0; c2(1:2*s+1)'*c0, c2(1:2*s+2)'*c1, c2'*c2];
                S(1:2,2:3)  = S(2:3,1:2)';
                if(1)
                    % + +   Use Alg. 2 from Fokkema's "Enhanced BiCGstab(L)" paper
                    %       to replace the above tau by one which yields a
                    %       larger biorthogonal fraction in the subsequent res.
                    y0          = [-1; S(2,2)\S(2,1); 0];
                    y2          = [ 0; S(2,2)\S(2,3); -1];
                    k0          = sqrt(y0'*S*y0);
                    k2          = sqrt(y2'*S*y2);
                    rho         = (y2'*S*y0) / ( k0 * k2 );
                    gamma_hat   = rho/abs(rho) * max( abs(rho) , 0.7 );
                    y0          = y0 - gamma_hat * k0/k2 * y2;
                    tau         = y0(2:3,1);
                else
                    tau         = S(2:3,2:3)\S(2:3,1);
                end
                % + update x and r residual-minimal (under the given degrees of freedom
                %   with respect to that obtained r is bi-orthogonal)
                e           = [c0;0;0] - tau(1) * [c1;0] - tau(2) * c2;
                x           = x + W(:,1:2*s+2) * ( tau(1)*[c0;0]+tau(2)*c1 );
                beta        = norm(e);
                ITER        = [ITER,iter];      % again: same iteration
                RESNORM     = [RESNORM,beta];
                fprintf('\t%i\t%3.3f\t%1.3g\n',iter,toc,beta);
                if( beta<=abstol || iter>=maxit)
                    x           = x + x_global;
                    x           = M2 \ x;
                    interior_plotter_function(ITER,RESNORM);
                    return;
                end
                r           =     W(:,1:2*s+3) * e;
        %         rel_orth_r  = norm(P'*(b-A*x))/norm(b-A*x)
        %         norm_r      = norm(b-A*x)
        %         cons_r      = norm(r - (b-A*x))/norm_r

                % Compute ortho-bi-ortho subsequent projectors V0 V1
                [Q,~]       = qr(tH(1:2*s+1,1:2*s)'*Qy1);
                Q2          = Q(:,s+1:2*s);
        %         Y(:,1:2*s+1)*tH(1:2*s+1,1:2*s) * Q2
                % + compute linear factors of new V1
                C           = tH(1:2*s+3,1:2*s) - tau(1) * tH(1:2*s+3,1:2*s+1) * tH(1:2*s+1,1:2*s) - tau(2) * tH(1:2*s+3,1:2*s+2) * tH(1:2*s+2,1:2*s+1) * tH(1:2*s+1,1:2*s);
                D           = - tau(2) * linsolve(Z,Y*tH(1:2*s+2,1:2*s+1)*tH(1:2*s+1,1:2*s),opts_tril);
                % + compute matrix R from QR decomposition of current [W,V1]
                B           = V1'*W;
                    
                [R,p]       = chol( [ W'*W, B' ; B , V1'*V1 ] );
                if(p==0)    % the cholesky decomposition breaks down
                    % This can mean either of the following
                    %   -   The condition of V1 from the last step is quite different
                    %       from 1. This would not happen without round-off.
                    %   -   The condition of the new W is bad. We do not assume
                    %       this at this step.
                    %   -   The minimal angle between W1 and W is zero and
                    %       round-off affected
                    [~,R]       = qr([W,V1],0);
                end
        %         cond([W V1]/R)
                % + compute matrx F, which inherits the non-orthogonality/illness of the 
                %   sum (V0-tau1*V1-tau2*V2) from the new basis matrices if they
                %   were computed straight forward.
                F           = R * [C;D] * Q2;
                [~,F]       = qr(F,0);
        %         F
                % + compute first version of G
                G           = Q2 / F;
                % + compute subsequent Z
                hZ          = Y(:,1:2*s+2) * tH(1:2*s+2,1:2*s+1) * tH(1:2*s+1,1:2*s) * G;
                [Qz,hZ]     = qr(hZ');
                G           = G * Qz;
                % + obtain V0, V1 from linear factors
                V0          = W(:,1:2*s+2) * ( [G;zeros(2,s)] - tau(1)*tH(1:2*s+2,1:2*s)*G - tau(2)*tH(1:2*s+2,1:2*s+1)*tH(1:2*s+1,1:2*s)*G ) - V0 * linsolve(Z,Y(:,1:2*s)*G,opts_tril);
                V1          = [W V1] * ( [C;D]*G );
        %         cond_V1 = cond(V1)
        %         cons_V1 = norm(V1-A*V0)
        %         cons_r  = norm(b-A*x-r)/beta
                % + write according Z
                Z           = -tau(2) * hZ';
        %         P'*V1
        %         Z
            %end case
        end % switch
        
%         replace_up  = false;    replace_down = false;
        % perform the actual restart
        if(replace_down)
            b_local         = b_local   - M1\(A * (M2\x));
            x_global        = x_global  + x;
            x               = zeros(N,1);
            r               = b_local;
            norm_b_local    = norm(b_local);
            % re-biorthogonalize the residual
            zeta            = linsolve(Z,P'*r,opts_tril);
            x               = x         + V0*zeta;
            r               = r         - V1*zeta;
            iter            = iter + 1;
            beta            = norm(r);
        elseif(replace_up)
            r               = b_local   - M1\(A * (M2\x));
            interm_res_max  = norm(r);
            % re-biorthogonalize the residual
            zeta            = linsolve(Z,P'*r,opts_tril);
            x               = x         + V0*zeta;
            r               = r         - V1*zeta;
            iter            = iter + 1;
            beta            = norm(r);
        end
        
        flag_write_U        = beta>abstol2 && strcmp(CYCLE_TYPE,'IDR(s)stab(1)');
        
        % write U
        if(flag_write_U)
            U               = V0;
            fprintf('Written out recycling data U for M(%i)stab(L).\n',s);
        end
        
    end
    
end

function [] = interior_plotter_function(ITER,RESNORM)
    hold all;
    plot(ITER,log10(RESNORM./RESNORM(1)));
end