getProlong_rs.m

%> Finds prolongation based on a CF-splitting, generated by the Ruge-Stueben
%> coarsening process for example.
%>
%> @param[in]   A           The fine level matrix
%> @param[in]   splitting   Vector of length nnu, whereas splitting(i)=1 indicates
%>                          a coarse grid variable and =-1 indicates a fine grid
%>                          variable.
%> @param[in]   num_cg_vars Number of variables on coarse grid.
%>
%> @param[out] P            Prolongation matrix.
%> @param[out] Inj          Injection matrix (for x in the multigrid cycle).
%> @param[out] coarse2fine  Mapping which maps coarse level numbering to fine level numbering
function [P, Inj, coarse2fine] = getProlong_rs(A, splitting, num_cg_vars)
    
    n = length(splitting);
    pointer = 1;
    P = sparse(n,num_cg_vars);
    Prow = zeros(n*num_cg_vars,1);
    Pcol = zeros(n*num_cg_vars,1);
    Pval = zeros(n*num_cg_vars,1);
    fine2coarse = -ones(1,n);
    coarse2fine = zeros(1,num_cg_vars);

    for i = 1:n
        if ( splitting(i) == 1 )
            coarse2fine(pointer) = i;
            fine2coarse(i) = pointer;
            pointer = pointer+1;
        end
    end
    
    Pptr = 1;
    
    for i=1:n
        if ( A(i,i) == 0 ) 
            disp('Prolongation encountered A(i,i)==0');
        end
        [~,col,val] = find(A(i,:));
        if ( splitting(i) == -1 )
            sum_neg_n = 0;
            sum_neg_p = 0;
            sum_pos_n = 0;
            sum_pos_p = 0;
            for j=1:numel(col)
                if ( val(j) < 0 )
                    sum_neg_n = sum_neg_n + val(j);
                    if ( splitting(col(j)) == 1 )
                        sum_neg_p = sum_neg_p + val(j);
                    end
                elseif ( val(j) > 0 )
                    sum_pos_n = sum_pos_n + val(j);
                    if ( splitting(col(j)) == 1 )
                        sum_pos_p = sum_pos_p + val(j);
                    end
                end
            end
            %if ( sum_neg_p == 0 && sum_neg_n ~= 0 )
            %    disp('Prolongation encountered sum_neg_p==0');
            %    A(i,:)
            %    splitting(col(:))
            %    sum_neg_n
            %end
            if ( sum_neg_n == 0 ) 
                alpha = 0;
            else
                alpha = sum_neg_n / sum_neg_p;
            end 
            if ( sum_pos_n == 0 ) 
                beta = 0;
            else
                beta = sum_pos_n / sum_pos_p;
            end
            tmp_alpha = -alpha / A(i,i);
            tmp_beta = -beta / A(i,i);
            for k=1:numel(col)
                if ( splitting(col(k)) == 1 )
                    if ( val(k) < 0 ) 
%                         P(i,fine2coarse(col(k))) =  val(k)*tmp_alpha;
                        Prow(Pptr) = i;
                        Pcol(Pptr) = fine2coarse(col(k));
                        Pval(Pptr) = val(k)*tmp_alpha;
                        Pptr = Pptr +1;
                        if ( val(k)*tmp_alpha == 0 ) 
                            disp('Prolongation (-) assigned weight 0');
                            A(i,:)
                            i
                            col(k)
                        end
                    else
%                         P(i,fine2coarse(col(k))) =  val(k)*tmp_beta;
                        Prow(Pptr) = i;
                        Pcol(Pptr) = fine2coarse(col(k));
                        Pval(Pptr) = val(k)*tmp_beta;
                        Pptr = Pptr +1;
                        if ( val(k)*tmp_beta == 0 ) 
                            disp('Prolongation (+) assigned weight 0');
                            i
                            col(k)
                        end
                    end
                end
            end
        end
    end
    
    for i=1:n
        if ( splitting(i) == 1 )
%             P(i,fine2coarse(i)) = 1;
            Prow(Pptr) = i;
            Pcol(Pptr) = fine2coarse(i);
            Pval(Pptr) = 1;
            Pptr = Pptr +1;
        end
    end
    
    P = sparse(Prow(1:Pptr-1),Pcol(1:Pptr-1),Pval(1:Pptr-1));

    Inj = P';

    % TEST
    % for i=1:size(P,1)
    %     sum=0;
    %     for j=1:size(P,2)
    %         sum = sum + P(i,j);
    %     end
    %     P(i,:) = P(i,:) .* 1/sum;
    % end
            
end