nested_dissection.m

This a matlab code for nested dissection.

function perm=nested_dissection(A)

% A is the adjacency matrix which must be renumbered
% A is sparse and symmetric
% The output perm is a permutation array which defines the
% renumbering


% The nested dissection renumbering scheme uses a graph
% partitioning algorithm to find a seperator set, that splits the
% graph into two components.  Once we have the
% partition, the seperator consists of nodes connected to
% cut-edges.

% The renumbering rule then is to number each of the components
% first and the seperator set last

% So, firstly calculate the partitioning ...

partOfA=graphPartitioning(A,psize,ga);

% Next extract the two components and seperator set
% Let's not worry how to do this initially -- just think about what
% the output of this stage should be and define a function.

% The two components and seperator set can each be represented by
% their own adjacency matrix -- so the output should be 3 adjacency
% matrices which we call C1, C2 and Sep

[C1, C2, Sep]=extractComponents(A, partOfA);

% Now renumber C1 ....
% We would like to renumber C1 using the *same* algorithm --
% i.e. through a recursive call to this function.

% If the size of C1 is M1 nodes, we should number them using the
% first M1 labels.  That's no problem the *first* time this function
% is run, clearly in that case the labels are the numbers between 1
% and M1.  

% However, we have to think about the general case.  We also want
% to renumber the nodes in C2, which, if the number of nodes in C2
% is M2,  requires a renumbering of the labels between M1+1 and M1+M2
% In general, the requirement will be to renumber some contiguous
% block of numbers between some minimum label and some maximum
% label.

% At this stage we realise that more inputs are required to the
% function - it should instead be defined as 

%% function perm = nested_dissection(A, minLabel, maxLabel);

% So if A an adjacency matrix which must be renumbered between
% labels minLabel and maxLabel, C1 must be renumbered between
% minLabel and minLabel+M1

% Now we have the following code ....

M1 = size(C1, 1);
M2 = size(C2, 1);
M3 = size(Sep,1);

perm = nested_dissection(C1, minLabel, minLabel+M1);
perm = nested_dissection(C2, minLabel+M1+1, minLabel+M1+M2-1);
perm = nested_dissection(Sep, minLabel+M1+M2, minLabel+M1+M2+M3-1);

% Notice how the permutation array perm is assigned to each call to
% the nested_dissection function.  This is required so that we get
% the result of the renumbering back from the function that created
% it.  But if each function call is going to update the perm array,
% we also need to ensure that it is given the most up-do-date perm
% as input.  Hence, we need to redefine the function as follows

% function perm = nested_dissection(A, minLabel, maxLabel, perm)

% and the code becomes

perm = nested_dissection(C1, minLabel, minLabel+M1, perm);
perm = nested_dissection(C2, minLabel+M1+1, minLabel+M1+M2-1, perm);
perm = nested_dissection(Sep, minLabel+M1+M2, minLabel+M1+M2+M3-1, ...
			 perm);

% Note the first time that nested_dissection is called, the perm
% array could be set as the identity array, perm(i)=i


% Okay how are we doing ... 

% Two problems should be becoming apparent...

% (1) If A is an NxN matrix and C1 is an M1xM1 matrix, how do we
% keep a track of which nodes of A are contained in C1? Note that
% the end result should be an Nx1 permutation matrix perm and
% clearly the minLabel, maxLabel etc that we are referring to above
% are labels in the set [1..N].  So to update the perm array, we
% need to know the relationship between nodes in A and nodes in C1
% (similarly C2 and Sep)  (We talked about this same problem last
% week when discussing how to do a k-way partitioning)

% (2) If the function remains as above it will keep on calling
% itself forever.  We need some criterion to stop the recursion.

% Let's go back through the lines of code already defined and try
% to sort out (1)

% The place to gather the relationship between nodes in A and nodes
% in C1, C2 and Sep is in the function that extracts the components
% We will make extractComponents return 3 extra arrays, nodeMap1,
% nodeMap2 and nodeMap3, such that

% (a) The length of nodeMap1 is M1
% (b) nodeMap1(i) = j means that node i in C1 corresponds to the
% global label j  (i.e. j belongs to the set [1..N]). 

% Similarly, for nodeMap2 and nodeMap3.

% BUT, remember that if nested_dissection is to be called
% recursively, taking in C1 on its second call etc., then there
% should be a node map array associated with the matrix that is
% input to the function.  Hence, we should define the function as
% follows

% function perm = nested_dissection(A, minLabel, maxLabel, perm, nodeMapA)

% and on the very first call to the function, we would have
% nodeMapA defined as the identity map i.e. nodeMapA(i)=i

% Hence, we re-write the code as follows:


[C1, C2, Sep, nodeMap1, nodeMap2, nodeMap3]= ...
    extractComponents(A, partOfA, nodeMapA);

M1 = size(C1, 1);
M2 = size(C2, 1);
M3 = size(Sep,1);

perm = nested_dissection(C1, minLabel, minLabel+M1, perm, nodeMap1);
perm = nested_dissection(C2, minLabel+M1+1, minLabel+M1+M2-1, perm, nodeMap2);
perm = nested_dissection(Sep, minLabel+M1+M2, minLabel+M1+M2+M3-1, ...
			 perm, nodeMap3);


% Finally, we get to problem (2) - the fact that, as written, this
% function is just going to keep on calling itself over and over
% again.

% We should stop the recursive call when the adjacency matrix is
% 'sufficiently' small.  How small you use is a parameter of your
% program - it could be, say, 100.  So if the matrix size is smaller
% than 100, we should instead apply some 'simple' numbering scheme
% e.g. reverse cuthill mckee, or the simplest example is a random
% numbering scheme.  With this modification, the code will look
% like


[C1, C2, Sep, nodeMap1, nodeMap2, nodeMap3]= ...
    extractComponents(A, artOfA, nodeMapA);

M1 = size(C1, 1);
M2 = size(C2, 1);
M3 = size(Sep,1);

if (M1>=100)
  perm = nested_dissection(C1, minLabel, minLabel+M1, perm, nodeMap1);
else
  perm = simpleRenumbering(C1, minLabel, minLabel+M1, perm, nodeMap1);
end

if (M2>=100)
  perm = nested_dissection(C2, minLabel+M1+1, minLabel+M1+M2-1, perm, ...
			   nodeMap2);
else
  perm = simpleRenumbering(C2, minLabel+M1+1, minLabel+M1+M2-1, perm, ...
			   nodeMap2);
end
  
if (M3>=100)
  perm = nested_dissection(Sep, minLabel+M1+M2, minLabel+M1+M2+M3-1, ...
			   perm, nodeMap3);
else
  perm = simpleRenumbering(Sep, minLabel+M1+M2, minLabel+M1+M2+M3-1, ...
			   perm, nodeMap3);
end

% and it remains to write the simpleRenumbering function