rugestuebencoarsen.m

五点差分型多重网格方法：各种插值算子的比较）

function RugeStuebenCoarsen(level)

% Ryan McKenzie
% Department of Computational Sciences
% University of Kentucky

%Adapted from a code by Dr. Gundolf Haase

amg_globals;

%------------------------- INITIALIZATION PHASE -------------------------%

%get the dimensions of the current level matrix
n = size( A(level).matrix , 1 );

V_set = 1:n; %set of all nodes in the current (fine) level
C_set = []; %set of coarse nodes, initialize to empty set
%NOTE: If you wish to pre-select coarse points, then the C_set must
%      be initialized accordingly

%determine and store the sets of strongly connected nodes S^(i)
for i = 1:n %for each node in the fine domain
    S(i).S = strong(i,V_set,level); %determine and store its strong set
end

%determine and store the transpose of the strongly connected sets ST^(i)
for i = 1:n %for each node in the fine domain
    S(i).ST = strongT(i,V_set,S,level); %determine and store its transpose strong set
end

F_set = []; %set of fine nodes, initialize to empty set
%initialize the fine point set according to pre-selected coarse points
for i = C_set %for each node in the coarse set
    %add any point that is not strongly connected to a pre-selected coarse
    %point into the set of fine points
    F_set = union( F_set , setdiff( S(i).ST , C_set ) );
    %if there are no pre-selected coarse points, F_set remains empty
end

%----------------------- INITIAL COARSENING PHASE -----------------------%

%generate the set of points to process as all points on the fine level
%which have not yet been pre-selected as fine or coarse (working set)
W_set = setdiff( V_set , union(F_set, C_set) );

%select a point with a maximal heuristic value and add it to the coarse set
while W_set %while there are still points to be processed
    M_Max=0; %reset maximum cardinality (number of strong connections)
    I_Max=-1; %reset the point ID for the maximal strong set
    
    for i=W_set %for each node in the set of nodes to process
        %here, the heuristic adds to a node's cardinality in two categories
        %   1)number of strong connections of node i to any other node
        %   2)number of weak connections of node i to fine points
        M_i = size(S(i).S , 2); %heuristic part 1
        M_i = M_i + size(intersect(S(i).ST, F_set),2); %heuristic part 2
        if M_i > M_Max; %if the heuristic score for this node is better than the max
            M_Max=M_i; I_Max = i; %make this node the new max
        end
    end
    if M_Max == 0 %if there were no appropriate points left to process
        F_set = setdiff( V_set, C_set ); %all remaining points are considered fine
    else %with an appropriate point discovered
        C_set = union(C_set, I_Max); %add the new point to the coarse set
        %add to the fine set any points weakly connected to the new node
        %which are not already in the coarse set
        F_set = union(F_set, setdiff(S(I_Max).ST, C_set) );
    end
    %now redefine the working set to reflect changes
    W_set = setdiff(V_set, union(F_set, C_set));
end

%------------------------ FINAL COARSENING PHASE ------------------------%

%skipped in this version...

%-------------------- CALCULATE INTERPOLATION WEIGHTS --------------------%

%Done by using a threshold strategy with the global variable SW_BOUND
%The SW_BOUND defines how strong in relation to the maximum off-diagonal
%element a connection must be to be considered in building the weight matrix

%Interpolation matrix W(level).weight will be quadratic in original numbering

%initialize temporary wieght matrix to the n by n identity matrix
W_Temp = eye(n,n); 

for i=F_set %for every fine point
    A_Max = max( abs( A(level).matrix(i,C_set) ) ); %get maximum off-diagonal entry in A
    C_i = []; %initialize the coarse coefficient set to be empty
    for j=C_set %for every coarse point
        %if the connection between fine point I and coarse point J is "strong enough"
        if abs( A(level).matrix(i,j) ) >= (SW_BOUND * A_Max);
            C_i = union(C_i, j); %add coarse point J to the coefficient set
        end
    end
    
    %get the sum of all coefficients operating on each coarse point
    %this will be used to scale the coefficients to the coarse level
    SS = sum( abs( A(level).matrix(i,C_i) ) ); 
    
    for j=C_i %for all coarse point coefficients in A
        %place the coefficient in the weight matrix after scaling
        W_Temp(i,j) = abs( A(level).matrix(i,j) / SS );
    end
end

%perform the appropriate permutations on the weight matrix
NC = size( C_set, 2); %get the number of coarse grid points
perm = [C_set , F_set]; %generate permutation set according to fine/coarse distrobution
W(level).weight = W_Temp(:,perm); %append permutation set to temporary weight matrix
%now stored in actual weight matrix

%finally, eliminate coefficients dealing with fine grid points
W(level).weight(:,NC+1:n) = [];

%we are left with an (n by NC) matrix in coarse-fine numbering

%-------------------------- COARSE LEVEL SETUP --------------------------%

%generate coarse grid matrix via the Galerkin approach
A(level+1).matrix = W(level).weight' * A(level).matrix * W(level).weight;