Decimation.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - Decimation
%    Copyright (C) 2009-2015 Peter Rakyta, Ph.D.
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see http://www.gnu.org/licenses/.
%
%
%> @addtogroup basic Basic Functionalities
%> @{
%> @file Decimation.m
%> @brief A class providing function handle to reduce the number of sites in the Hamiltonian via decimation procedure.
%> @} 
%> @brief A class providing function handle to reduce the number of sites in the Hamiltonian via decimation procedure.
%> The Decimation parameter in structure #Opt picks the following decimation procedures \n
%> Opt.Decimation=1 or Opt.Decimation=2 picks procedure #decimation_1_and_2 \n
%> Opt.Decimation=3 picks procedure #decimation_3 \n
%> Opt.Decimation=4 picks procedure #decimation_4 \n
%%
classdef Decimation < handle & Messages

    
    properties ( Access = public )  
%> The function handle of the decimation method
        DecimationFunc 
    end

methods ( Access = public )
%% Contructor of the class
%> @brief Constructor of the class.
%> @param Opt An instance of the structure #Opt.
%> @param varargin Cell array of optional parameters:
%> @param 'ForcedDecimation' Overrides the parameter value Opt.Decimation in structure #Opt.
%> @return An instance of the class
    function obj = Decimation( Opt, varargin)
        obj = obj@Messages( Opt );
        
        p = inputParser;
        p.addParameter('ForcedDecimation', []);
        p.parse(varargin{:});
        
        ForcedDecimation = p.Results.ForcedDecimation;

        if ~isempty(ForcedDecimation)
            obj.Opt.Decimation = ForcedDecimation;
        end   
        
        if isempty( Opt.Decimation )
           error(['EQuUs:', class(obj), ':Decimation:Wrong_Decimation_Parameter'], 'Empty Decimation parameter given in Opt.Decimation.') 
        end
            

        if Opt.Decimation == 3    
            obj.DecimationFunc = @(E, Way2Hamiltonian, NameOfH, NameOfRegular_sites)(obj.decimation_3(E, Way2Hamiltonian, NameOfH, NameOfRegular_sites));
        elseif Opt.Decimation == 4    
            obj.DecimationFunc = @( E, Way2Hamiltonian, NameOfH, NameOfRegular_sites )(obj.decimation_4(E, Way2Hamiltonian, NameOfH, NameOfRegular_sites ));
        elseif  Opt.Decimation == 1 || Opt.Decimation == 2
            obj.DecimationFunc = @(E, Way2Hamiltonian, NameOfH, NameOfRegular_sites)(obj.decimation_1_and_2(E, Way2Hamiltonian, NameOfH, NameOfRegular_sites));
        else
            warning(['EQuUs:', class(obj), ':Decimation:Wrong_Decimation_Parameter'], 'Wrong Decimation parameter given! Using default 4.')
            obj.DecimationFunc = @(E, Way2Hamiltonian, NameOfH, NameOfRegular_sites)(obj.decimation_4(E, Way2Hamiltonian, NameOfH, NameOfRegular_sites));
        end

    end  
    
%% decimation_4
%> @brief Algorithm to reduce the number of the sites in the Hamiltonian via decimation. The number of the decimated sites in one turn varies dynamically in order to avoid badly conditioned matrices. Equally fast as #decimation_3, but more stable and more accurate.
%> @param E The energy value
%> @param Way2Hamiltonian An instance of class CreateHamiltonians (or a derived class).
%> @param NameOfH Name of the attribute storing the Hamiltonian to be decimate  (typically use 'Hscatter')
%> @param NameOfRegular_sites Name of the attribute storing the idexes of sites to be kept after the decimation.  (typically use 'kulso_szabfokok')
    function decimation_4(obj, E, Way2Hamiltonian, NameOfH, NameOfRegular_sites )
    
     if Way2Hamiltonian.Read('HamiltoniansDecimated')
          obj.display(['EQuUs:', class(obj), ':decimation_4: The Hamiltonian is already decimated'],1);
          return
        elseif ~Way2Hamiltonian.Read('HamiltoniansCreated')
          obj.display(['EQuUs:', class(obj), ':decimation_4: The Hamiltonian needs to be created before the decimation'],1);
          return
        end
        
        %blokk_size_max = 50;
        blokk_size_max = obj.Opt.Decimation_block;
                
        Hamiltoni       = Way2Hamiltonian.Read( NameOfH );
        regular_sites = Way2Hamiltonian.Read( NameOfRegular_sites );
        
        if length(regular_sites) >= size(Hamiltoni,1)
            obj.display(['EQuUs:', class(obj), ':decimation_4: Nothing to be decimated']);
            return
        end
        
        Way2Hamiltonian.Clear( NameOfH );
        
        sorting_sites();
        blocks = blocks_to_be_decimated();        
        
        for idx = length(blocks):-1:2
            decimation_of_block2( blocks(idx)-blocks(idx-1) );
        end
                       
                
        Way2Hamiltonian.Write( NameOfH, Hamiltoni );
        Way2Hamiltonian.Write('HamiltoniansDecimated', true)
        
        

    %-----------------------------------------------------------------
    % nested functions        
        
        %------------------------------------------------------------------
        % determines the blocks of the Hamiltonian to be decimated out
        function blocks = blocks_to_be_decimated()
            blocks = length(regular_sites):blokk_size_max:size(Hamiltoni,1);
            if blocks(end) <  size(Hamiltoni,1)  
                blocks(end+1) = size(Hamiltoni,1);
            end
            
            if length(blocks) < 2
                error(['EQuUs:', class(obj), ':decimation_4'], 'Error in determining the blocks to be decimated');
            end
        end   
        
               
        %------------------------------------------------------------------
        % sites to be decimated out are sorted to the end
        function sorting_sites()
            indexes = false( size(Hamiltoni,1), 1);
            indexes(regular_sites) = true;
            
            Hamiltoni = [Hamiltoni(indexes,indexes),  Hamiltoni(indexes,~indexes);
                                Hamiltoni(~indexes,indexes), Hamiltoni(~indexes,~indexes)];           
            
            
            regular_sites = 1:length(regular_sites);
            
            
        end
        
        
        %------------------------------------------------------------------
        % decimates out the given block in the Hamiltonian
        function decimation_of_block2( num )
            
            while num > 0
            
                nevezo  = E*eye(num) - Hamiltoni(end-num+1:end, end-num+1:end);
                if issparse(nevezo)
                    cond_number = rcond(full(nevezo));
                else
                    cond_number = rcond(nevezo);
                end
                tolerance   = 1e-4;
            
                %if abs(cond_number-1) > tolerance
                if cond_number < tolerance && num  > 1
                    num2 = increase_condition_number();
                    decimation_of_block2( num2 );
                    num = num - num2;
                else
                    % here the given block with good condition number is
                    % decimated out
                    [~, col_idx] = find(  Hamiltoni(end-num+1:end, 1:end-num)   );
                    col_blocks = unique(floor(col_idx/blokk_size_max));
                    for jdx = 1:length(col_blocks)
                        col_min = max( [col_blocks(jdx)*blokk_size_max, 1] );
                        col_max = min( [size(Hamiltoni,2)-num, (col_blocks(jdx)+1)*blokk_size_max-1] );
                        cols = col_min:col_max;                        
                        H_21 = Hamiltoni( end-num+1:end, cols);
                        nevezo2   = nevezo\H_21;
                        
                        [row_idx2, ~] = find(  Hamiltoni(1:end-num, end-num+1:end)   );
                        row_blocks = unique(floor(row_idx2/blokk_size_max));
                        for kdx = 1:length(row_blocks)
                            row_min = max( [row_blocks(kdx)*blokk_size_max, 1] );
                            row_max = min( [size(Hamiltoni,1)-num, (row_blocks(kdx)+1)*blokk_size_max-1] );
                            rows = row_min:row_max;
                            H_12 = Hamiltoni( rows, end-num+1:end );
                            H_11 = H_12*nevezo2;
                            
                            if issparse(Hamiltoni)
                                size_H = size(Hamiltoni,1);
                                
                                [rows_H11, cols_H11, elements_H11] = find(H_11);
                                if size(H_11,1) == 1
                                    rows_H11 = transpose(rows_H11);
                                    cols_H11 = transpose(cols_H11);
                                    elements_H11 = transpose(elements_H11);
                                end
                                rows_H11 = rows_H11 + row_min - 1 ;
                                cols_H11 = cols_H11 + col_min - 1;                          
                                
                                [rows_H, cols_H, elements] = find( Hamiltoni );
                                Hamiltoni = [];
                                Hamiltoni = sparse( [rows_H; rows_H11], [cols_H; cols_H11], [elements; elements_H11], size_H, size_H );
                            else
                                Hamiltoni(rows, cols) = Hamiltoni( rows, cols ) + H_11;
                            end
                            
                        end
                        
                    end 
                    
                    % delete the decimated sites from the Hamiltonian
                    if issparse(Hamiltoni)
                                size_H = size(Hamiltoni,1);   
                                [rows, cols, elements] = find( Hamiltoni(1:end-num, 1:end-num) );
                                Hamiltoni = sparse( rows, cols, elements, size_H-num, size_H-num );
                    else
                                Hamiltoni(:, end-num+1:end ) = [];
                                Hamiltoni(end-num+1:end, :) = [];
                    end
                    
                    
                    % convert to full matrix if sparse is not memory
                    % efficient anymore
                    whos_H = whos('Hamiltoni');
                    if whos_H.size(1)*whos_H.size(2)*16 <= whos_H.bytes && whos_H.sparse == 1
                        Hamiltoni = full(Hamiltoni);
                        obj.display(['EQuUs:', class(obj), 'Sparse Hamiltonian was converted to full.']);
                    end 
                            
                    
                    num = 0;
                    
                end
            end
            
            
            
            %-----------------------------------
            % determine block size for which the matrix can be inverted
            function num_tmp = increase_condition_number()
                num_tmp = num;
                while cond_number < tolerance
                    num_tmp = floor(num_tmp/2);
                    if num_tmp == 1
                        break
                    end
                    nevezo = E*eye(num_tmp) - Hamiltoni(end-num_tmp+1:end, end-num_tmp+1:end);
                    if issparse(nevezo)
                        cond_number = rcond(full(nevezo));
                    else
                        cond_number = rcond(nevezo);
                    end
                end 
                
                nevezo = [];
            end
            
        end        
        
    %-----------------------------------------------------------------
    % end nested functions       
                
        
    end

%% decimation_3
%> @brief Algorithm to reduce the number of the sites in the Hamiltonian via decimation. The number of the decimated sites in one turn varies dynamically in order to avoid badly conditioned matrices.
%> @param E The energy value
%> @param Way2Hamiltonian An instance of class CreateHamiltonians (or a derived class).
%> @param NameOfH Name of the attribute storing the Hamiltonian to be decimate  (typically use 'Hscatter')
%> @param NameOfRegular_sites Name of the attribute storing the idexes of sites to be kept after the decimation.  (typically use 'kulso_szabfokok')
    function decimation_3(obj, E, Way2Hamiltonian, NameOfH, NameOfRegular_sites )
    
     if Way2Hamiltonian.Read('HamiltoniansDecimated')
          obj.display(['EQuUs:', class(obj), ':decimation_3: The Hamiltonian is already decimated'],1);
          return
        elseif ~Way2Hamiltonian.Read('HamiltoniansCreated')
          obj.display(['EQuUs:', class(obj), ':decimation_3: The Hamiltonian needs to be created before the decimation'],1);
          return
          end
        
        %blokk_size_max = 50;
        blokk_size_max = obj.Opt.Decimation_block;
        
        Hamiltoni       = Way2Hamiltonian.Read( NameOfH );
        regular_sites = Way2Hamiltonian.Read( NameOfRegular_sites );
        Way2Hamiltonian.Clear( NameOfRegular_sites );
        
        obj.setCoordinates( Way2Hamiltonian, regular_sites );
        
        megmaradt_szabfokok = sites_after_decimation(regular_sites);
        Way2Hamiltonian.Write( NameOfRegular_sites, megmaradt_szabfokok );
        
        Way2Hamiltonian.Clear( NameOfH );
        
        H_size          = size(Hamiltoni ,1);
        blokkok         = indexing_of_regular_sites();
        belso_szabfokok = find(blokkok{2} == 'b');
        
        whos_H          = whos('Hamiltoni');
        if whos_H.size(1)*whos_H.size(2)*16 <= whos_H.bytes && whos_H.sparse == 1
            Hamiltoni = full(Hamiltoni);
            obj.display('A Hamltoni mar nem ritka');
        end 

        Idx            = length(belso_szabfokok);
        iidx           = 1;
        
        while iidx <= Idx
            idx = belso_szabfokok(iidx);
            decimation_of_block();
            
            whos_H = whos('Hamiltoni');
            if whos_H.size(1)*whos_H.size(2)*16 <= whos_H.bytes && whos_H.sparse == 1
                Hamiltoni = full(Hamiltoni);
                obj.display('A Hamltoni mar nem ritka');
            end 
            
            iidx = iidx + 1;
        end
        
        
                
        Way2Hamiltonian.Write( NameOfH, Hamiltoni );
        Way2Hamiltonian.Write('HamiltoniansDecimated', 1)
        
    %-----------------------------------------------------------------
    % nested functions
        
        %------------------------------------------------------------------
        function decimation_of_block()
            
            oszlop  = [sum(blokkok{1}(1:idx-1))+1,sum(blokkok{1}(1:idx))];
            nevezo  = E*eye(blokkok{1}(idx)) - Hamiltoni(oszlop(1):oszlop(2),oszlop(1):oszlop(2));
            cond_number = rcond(nevezo); 
            tolerance   = 1e-4;
            
            %if abs(cond_number-1) > tolerance
            if cond_number < tolerance
                clear('nevezo'); 
                [oszlop, nevezo] = increase_condition_number();
            end
            
            if cond_number < tolerance
                obj.display(num2str(cond(nevezo)) )
            end
     
            [sor_idx,oszlop_idx] = find(Hamiltoni(oszlop(1):oszlop(2),:));
            oszlop_szamok        = oszlop_szamolo(oszlop_idx', oszlop, blokk_size_max);
            Kdx                  = size(oszlop_szamok,1);
            kdx                  = 1;
            while kdx <= Kdx
                oszlop_k  = oszlop_szamok(kdx,:); 
                jdx       = 1;
                H_darab21 = Hamiltoni(oszlop(1):oszlop(2),oszlop_k(1):oszlop_k(2));
                nevezo2   = nevezo\H_darab21;
                while jdx <= Kdx
                    oszlop_j  = oszlop_szamok(jdx,:);
                    H_darab11 = Hamiltoni(oszlop_j(1):oszlop_j(2),oszlop_k(1):oszlop_k(2));
                    H_darab12 = Hamiltoni(oszlop_j(1):oszlop_j(2),oszlop(1):oszlop(2));
                    H_darab11 = H_darab11 + H_darab12*nevezo2;
                    
                    Hamiltoni(oszlop_j(1):oszlop_j(2),oszlop_k(1):oszlop_k(2)) = H_darab11;
                    jdx       = jdx + 1;
                end
                kdx = kdx + 1;
            end 
        
            if issparse(Hamiltoni)
                [sor_idx,oszlop_idx,elemek] = find(Hamiltoni); 
                Hamiltoni           = [];
                indexek             = logical(((oszlop(1)<=sor_idx) & (sor_idx <=oszlop(2))) | ((oszlop(1)<=oszlop_idx) & (oszlop_idx <=oszlop(2))));
                sor_idx(indexek)    = [];
                oszlop_idx(indexek) = [];
                elemek(indexek)     = [];
                indexek             = logical(sor_idx>oszlop(2));
                sor_idx(indexek)    = sor_idx(indexek) - blokkok{1}(idx);
                indexek                  = logical(oszlop_idx>oszlop(2));
                oszlop_idx(indexek) = oszlop_idx(indexek) - blokkok{1}(idx);
                clear indexek;
                blokkok{1}(idx) = [];
                blokkok{2}(idx) = [];
                Hamiltoni       = sparse(sor_idx,oszlop_idx,elemek,sum(blokkok{1}),sum(blokkok{1}));
                clear sor_idx oszlop_idx elemek;
            else   
                Hamiltoni(oszlop(1):oszlop(2),:) = [];
                Hamiltoni(:,oszlop(1):oszlop(2)) = [];
                blokkok{1}(idx) = [];
                blokkok{2}(idx) = [];
                
            end                
                belso_szabfokok = belso_szabfokok -1;
                
                
            %------------------------------------------------------------------               
            % kiszamolja hogy a kidecimalando szabfokokat hogy kell beoszlopozni
            function oszlop_szamok = oszlop_szamolo(oszlop_idx, oszlop_decimalo, blokk_size_max)
    
            blokk_size_max = 150*blokk_size_max;
            ldx = 1;
            Ldx = length(oszlop_idx);

            if oszlop_idx(1) >= oszlop_decimalo(1)
                while oszlop_idx(ldx) <= oszlop_decimalo(2)
                    ldx = ldx + 1;
                end
                oszlop_szamok = [oszlop_idx(ldx),oszlop_idx(ldx)];    
            else
                oszlop_szamok = [oszlop_idx(1),oszlop_idx(1)];
            end

            ldx = ldx + 1;

            while (ldx <= Ldx) && (oszlop_decimalo(1) > oszlop_idx(ldx))
                if (oszlop_idx(ldx)-oszlop_szamok(end,1) > blokk_size_max)
                    oszlop_szamok(end,2) = oszlop_idx(ldx-1);
                    oszlop_szamok        = [oszlop_szamok;[oszlop_idx(ldx),oszlop_idx(ldx)]];       
                end    
                ldx = ldx + 1;
            end

            oszlop_szamok(end,2) = oszlop_idx(ldx-1);

            while (ldx < Ldx) && (oszlop_decimalo(2) > oszlop_idx(ldx))
                ldx = ldx + 1;
            end

            oszlop_szamok        = [oszlop_szamok;[oszlop_idx(ldx),oszlop_idx(ldx)]];

            while ldx <= Ldx
                if (oszlop_idx(ldx)-oszlop_szamok(end,1) > blokk_size_max)
                    oszlop_szamok(end,2) = oszlop_idx(ldx-1);
                    oszlop_szamok        = [oszlop_szamok;[oszlop_idx(ldx),oszlop_idx(ldx)]];     
                end
                ldx = ldx + 1;
            end

            oszlop_szamok(end,2) = oszlop_idx(ldx-1);
            %oszlop_szamok = [1,oszlop_decimalo(1)-1;oszlop_decimalo(2)+1,H_size];

            end
            %------------------------------------------------------------------
        
        
            %------------------------------------------------------------------
            function [oszlop,nevezo] = increase_condition_number()
        
            %display(['condition number = ',num2str(cond_number)]); 
            blokkok{1} = [blokkok{1}(1:idx);blokkok{1}(idx:end)];
            blokkok{2} = [blokkok{2}(1:idx);blokkok{2}(idx:end)];
            egyik_fele = blokkok{2}(idx);
            blokkok{1}(idx+1) = 0;
            
            %while abs(cond_number-1) >= tolerance && egyik_fele >= 2
            while cond_number <= tolerance && egyik_fele >= 2
                egyik_fele        = round(blokkok{1}(idx)*2.99/4);
                masik_fele        = blokkok{1}(idx) - egyik_fele;
                blokkok{1}(idx)   = egyik_fele;
                blokkok{1}(idx+1) = masik_fele + blokkok{1}(idx+1);
                oszlop            = [sum(blokkok{1}(1:idx-1))+1,sum(blokkok{1}(1:idx))];
                nevezo            = E*eye(blokkok{1}(idx)) - Hamiltoni(oszlop(1):oszlop(2),oszlop(1):oszlop(2));
                cond_number       = rcond(nevezo);
            end
            %display(num2str(blokkok{1}(idx)))
            %display(['atjavitott condition number = ',num2str(cond_number),' valamint kidecimalando szabfokok = ',num2str(blokkok{1}(idx))]);
            belso_szabfokok = [belso_szabfokok(1:iidx);belso_szabfokok(iidx:end)];
            belso_szabfokok(iidx+1:end) = belso_szabfokok(iidx+1:end) + 1;
            Idx = length(belso_szabfokok);
        
        
        
            end
            %------------------------------------------------------------------
                
                
        end
        %------------------------------------------------------------------
        
        
        %------------------------------------------------------------------
        % Determines the blocks of regular sites.
        function blokkok = indexing_of_regular_sites()
        
        
            regular_sites = sort(regular_sites);
            iii_max         = length(regular_sites);
            iiidx           = 1;
            blokkok         = cell(2,1);
        
            while iiidx <= iii_max
            
                szabfok_tipus = [];
                blokk_size    = 1;
                if (sum(blokkok{1}) + blokk_size) == regular_sites(iiidx)
                    while (iiidx == iii_max || regular_sites(iiidx+1) == regular_sites(iiidx)+1) && blokk_size <= blokk_size_max
                        blokk_size = blokk_size + 1;
                        iiidx      = iiidx + 1;
                        if iiidx >= iii_max; blokk_size = blokk_size-1; break; end
                    end
                    if blokk_size == 1 && iiidx < iii_max ; iiidx = iiidx + 1; end
                    szabfok_tipus = 'k';  % kulso szabfok 
                elseif (sum(blokkok{1}) + blokk_size) < regular_sites(iiidx)
                    while (sum(blokkok{1}) + blokk_size + 1 < regular_sites(iiidx)) && blokk_size <= blokk_size_max
                        blokk_size = blokk_size + 1;
                    end                
                    szabfok_tipus = 'b'; % kidecimalando belso szabfokok
                elseif isempty(szabfok_tipus) && (sum(blokkok{1}) + blokk_size) > regular_sites(iiidx)
                    iiidx = iiidx + 1;
                    szabfok_tipus = [];
                    blokk_size    = [];
                end
            
                blokkok{1} = [blokkok{1};blokk_size];
                blokkok{2} = [blokkok{2};szabfok_tipus];
     
            
            end
        
            szabfokok = H_size-sum(blokkok{1});
        
            if szabfokok > 0
                blokk_size_max     = blokk_size_max + 1;
                blokkok_kieg       = ones((szabfokok-mod(szabfokok,blokk_size_max))/blokk_size_max,1)*blokk_size_max;
                szabfok_tipus_kieg = char(ones((szabfokok-mod(szabfokok,blokk_size_max))/blokk_size_max,1)*'b');
                blokkok{1}         = [blokkok{1};blokkok_kieg;mod(szabfokok,blokk_size_max)];
                blokkok{2}         = [blokkok{2};szabfok_tipus_kieg;'b'];
            end
        
        end
        %------------------------------------------------------------------
        
        
        %------------------------------------------------------------------
        % Determines the sites after decimation
        function remaining_sites = sites_after_decimation(regular_sites)
        
            L          = length(regular_sites);
            idx        = 1;
            ertek_temp = [];
            Maximum    = max(regular_sites);
            remaining_sites = zeros(L,1);
            while idx <= L
                szabfok_max   = max(remaining_sites);
                [ertek,index] = min(regular_sites);
                if ertek_temp == ertek
                    remaining_sites(index) = szabfok_max;
                else
                    ertek_temp    = ertek;
                    remaining_sites(index) = szabfok_max+1;
                end
                regular_sites(index) = Maximum+1;
                idx = idx + 1;
            end
    
    
   
        end
        
    % end nested functions        
    %------------------------------------------------------------------       
    
        
        
    end
    
%% decimation_1_and_2
%> @brief Algorithm to reduce the number of sites in the Hamiltonian via decimation. The number of the decimated sites is fixed during the calculations.
%> @param E The energy value
%> @param Way2Hamiltonian An instance of class CreateHamiltonians (or a derived class).
%> @param NameOfH Name of the attribute storing the Hamiltonian to be decimate (typically use 'Hscatter')
%> @param NameOfRegular_sites Name of the attribute storing the idexes of sites to be kept after the decimation  (typically use 'kulso_szabfokok').
function decimation_1_and_2(obj, E, Way2Hamiltonian, NameOfH, NameOfRegular_sites )

    if Way2Hamiltonian.Read('HamiltoniansDecimated')
     obj.display(['EQuUs:', class(obj), ':decimation_1_and_2: The Hamiltonian is already decimated'],1);
     return
    elseif ~Way2Hamiltonian.Read('HamiltoniansCreated')
          obj.display(['EQuUs:', class(obj), ':decimation_1_and_2: The Hamiltonian needs to be created before the decimation'],1);
          return
    end
    
    Hamiltoni       = Way2Hamiltonian.Read( NameOfH );
    regular_sites = Way2Hamiltonian.Read( NameOfRegular_sites );
    Way2Hamiltonian.Clear( NameOfH );
    Way2Hamiltonian.Clear( NameOfRegular_sites );
    
    obj.setCoordinates( Way2Hamiltonian, regular_sites );
    
    Idx = size(Hamiltoni, 1);
    idx = 1;
    [minidx,minkoord]   = min(regular_sites);
    minidxtemp          = 0;
    szabfokok_szama     = length(regular_sites); 
    megmaradt_szabfokok = zeros(szabfokok_szama,1);



    while idx <= Idx
    
        if idx < minidx
            if obj.Opt.Decimation == 1
                sitewise_decimation;
                regular_sites   = regular_sites - 1;
                minidx            = minidx - 1;
                Idx               = size(Hamiltoni, 1);
            elseif obj.Opt.Decimation == 2
                blokk_size_max    = obj.Opt.Decimation_block;
                blokk_size        = minidx - idx;
                if blokk_size > blokk_size_max
                    blokk_size = blokk_size_max;
                end
                decimalas_blokkonkent(blokk_size);
                s = whos('Hamiltoni');
                if s.size(1)*s.size(2)*16 <= s.bytes && s.sparse == 1
                    Hamiltoni = full(Hamiltoni);
                else
                    Hamiltoni = sparse(Hamiltoni);
                end
                regular_sites   = regular_sites - blokk_size;
                minidx            = minidx - blokk_size;
                Idx               = size(Hamiltoni, 1);
            else
                error(['EQuUs:', class(obj), ':decimation_1_and_2'], 'Wrong decimation parameter: Opt_Decimation');
            end       
        
        else
            if minidx ~= minidxtemp
                idx = idx + 1;
            end   
  
            if minidx <= Idx
                megmaradt_szabfokok(minkoord) = idx-1;
                regular_sites(minkoord)     = Idx+1;
            end
                minidxtemp        = minidx;
                [minidx,minkoord] = min(regular_sites);
        end    

      
    end

    if megmaradt_szabfokok(minkoord) == 0
        megmaradt_szabfokok(minkoord) = minidx;
    end
    
    Way2Hamiltonian.Write( NameOfH, Hamiltoni );
    Way2Hamiltonian.Write( NameOfRegular_sites, megmaradt_szabfokok );
    Way2Hamiltonian.Write( 'HamiltoniansDecimated', 1 );
    
    
    %-----------------------------------------------------------------
    % nested functions
    
    %------------------------------------------------------------------
    function decimalas_blokkonkent(blokk_size)


        Hsize = Idx;

        E_matrix = eye(blokk_size)*E;

        if idx == 1
            Hamiltoni = Hamiltoni(blokk_size+1:Hsize,blokk_size+1:Hsize) + ...
                Hamiltoni(blokk_size+1:Hsize,1:blokk_size) /(E_matrix-Hamiltoni(1:blokk_size,1:blokk_size) * Hamiltoni(1:blokk_size,blokk_size+1:Hsize));
    
        elseif idx+blokk_size-1 == Hsize
    
            Hamiltoni = Hamiltoni(1:Hsize-blokk_size,1:Hsize-blokk_size) + ...
                Hamiltoni(1:Hsize-blokk_size,Hsize-blokk_size+1:Hsize) /(E_matrix-Hamiltoni(Hsize-blokk_size+1:Hsize,Hsize-blokk_size+1:Hsize)) * Hamiltoni(Hsize-blokk_size+1:Hsize,1:Hsize-blokk_size);
       
        else
            blokk_size = blokk_size-1;
            nevezo  = (E_matrix-Hamiltoni(idx:idx+blokk_size,idx:idx+blokk_size))\eye(size(E_matrix));
            Hdeci11 = Hamiltoni(1:idx-1,1:idx-1) + Hamiltoni(1:idx-1,idx:idx+blokk_size)*nevezo*Hamiltoni(idx:idx+blokk_size,1:idx-1);
            Hdeci12 = Hamiltoni(1:idx-1,idx+blokk_size+1:Hsize) + Hamiltoni(1:idx-1,idx:idx+blokk_size)*nevezo*Hamiltoni(idx:idx+blokk_size,idx+blokk_size+1:Hsize);
            Hdeci21 = Hamiltoni(idx+blokk_size+1:Hsize,1:idx-1) + Hamiltoni(idx+blokk_size+1:Hsize,idx:idx+blokk_size)*nevezo*Hamiltoni(idx:idx+blokk_size,1:idx-1);
            Hdeci22 = Hamiltoni(idx+blokk_size+1:Hsize,idx+blokk_size+1:Hsize) + Hamiltoni(idx+blokk_size+1:Hsize,idx:idx+blokk_size)*nevezo*Hamiltoni(idx:idx+blokk_size,idx+blokk_size+1:Hsize);
            Hamiltoni = [];
            Hamiltoni = [Hdeci11,Hdeci12;Hdeci21,Hdeci22];
        end
        
    end
    %------------------------------------------------------------------
    
    
    %------------------------------------------------------------------
    % decimalas pontonkent
    function sitewise_decimation()
        
        
        
        
        nem_nulla_idx = find(Hamiltoni(:,idx));
        rekurzioszam  = max(size(nem_nulla_idx));

        sor    = 1;
        while sor <= rekurzioszam
       
            if ~(nem_nulla_idx(sor) == idx)
                oszlop = 1;
                while oszlop <= rekurzioszam
                    if nem_nulla_idx(oszlop) == idx
                        oszlop = oszlop + 1;
                    else
                        Hamiltoni(nem_nulla_idx(sor),nem_nulla_idx(oszlop)) = Hamiltoni(nem_nulla_idx(sor),nem_nulla_idx(oszlop)) + ...
                            Hamiltoni(nem_nulla_idx(sor),idx)*Hamiltoni(idx,nem_nulla_idx(oszlop)) / (E-Hamiltoni(idx,idx));
                        oszlop = oszlop + 1;
                    end
                end
            end
            sor = sor + 1; 
        end
        Hamiltoni(idx,:) = [];
        Hamiltoni(:,idx) = [];

    end
    %------------------------------------------------------------------        
    
    % end nested functions
            
end  
    
end % methods public

methods ( Access = protected )
%% setCoordinates
%> @brief Sets the coordinates of the remaining sites after the decimation is finshed.
%> @param Way2Hamiltonian An instance of class CreateHamiltonians (or a derived class).
%> @param regular_sites Array of idexes of sites kept after the decimation
    function setCoordinates( obj, Way2Hamiltonian, regular_sites )
        
        coordinates     = Way2Hamiltonian.Read( 'coordinates' );
        
        if isempty(coordinates)
            obj.display(['EQuUs:', class(obj), ':setCoordinates: structure coordinates is empty!'])
            return
        end
        
        indexes = logical( regular_sites <= length(coordinates.x) );
        
        coordinates.x     = coordinates.x(regular_sites(indexes));
        coordinates.y     = coordinates.y(regular_sites(indexes));   
        
        if isprop( coordinates, 'spinup' ) && ~isempty( coordinates.spinup )
            coordinates.spinup = coordinates.spinup(regular_sites(indexes));
        end
        
        if isprop( coordinates, 'BdG_u' ) && ~isempty( coordinates.BdG_u )
            coordinates.BdG_u = coordinates.BdG_u(regular_sites(indexes));
        end
        
        Way2Hamiltonian.Write('coordinates', coordinates );        
        
    end        

end %methods protected

end % Global end