推荐 | A suite of drivers to compute the PageRank vector for a directed graph.
The Pagerank package provides a single driver call capable of running different pagerank algorithm.
1) PageRank with power iterations
2) PageRank with gauss-seidel iterations
3) PageRank as a linear system (bicgstab and gmres solvers)
4) PageRank with the arnoldi factorization
5) Approximate personalized PageRank
The included matlab file describes all the options and parameters.
For non-linux/win32 platforms, you must compile the included .c files
mex ?setup
mex pagerank_mult.c
mex pagerank_gs_mult.c
mex spmatvec_mult.c
mex spmatvec_transmult.c
The mex files provide optimized matrix-vector multiplication routines for the PageRank code and improve the performance by 10-15%.
There is a special driver included to run a set of algorithms on the same web-graph and report comparative results.
pagerank(A,struct('alg','eval'));
Cite As
David Gleich (2023). pagerank (https://www.mathworks.com/matlabcentral/fileexchange/11613-pagerank), MATLAB Central File Exchange. Retrieved.
function [x flag hist dt] = pagerank(A,optionsu)
% PAGERANK Compute the PageRank for a directed graph.
%
% [p flag hist dt] = pagerank(A)
%
% Compute the pagerank vector p for the directed graph A, with
% teleportation probability (1-c).
%
% flag is 1 if the method converged; hist returns the convergence history
% and dt is the total time spent solving the system
%
% The matrix A should have the outlinks represented in the rows.
%
% This driver can compute PageRank using 4 different algorithms,
% the default algorithm is the Arnoldi iteration for PageRank due to
% Grief and Golub. Other algorithms include gauss-seidel iterations,
% power iterations, a linear system formulation, or an approximate
% PageRank formulation.
%
% The output p satisfies p = c A'*D^{+} p + c d'*p v + (1-c) v and
% norm(p,1) = 1.
%
% The power method solves the eigensystem x = P''^T x.
% The linear system solves the system (I-cP^T)x = (1-c)v.
% The dense method uses "\" on I-cP^T which the LU factorization.
%
% To specify a different solver for the linear system, use an anonymous
% function wrapper around one of Matlab's solver calls. To use GMRES,
% call pagerank(..., struct('linsolver', ...
% @(f,v,tol,its) gmres(f,v,[],tol, its)))
%
% Note 1: the 'approx' algorithm is the PageRank approximate personalized
% PageRank algorithm due to Gleich and Polito. It creates a set of
% active pages and runs until either norm(p(boundary),1) < options.bp or
% norm(p(boundary),inf) < options.bp, where the boundary is defined as
% the set of pages that have a non-zero personalized PageRank but are not
% in the set of active pages. As options.bp -> 0, both of these
% approximations compute the actual personalized PageRank vector.
%
% Note 2: the 'eval' algorithm evaluates five algorithms to compute the
% PageRank vector and summarizes the results in a report. The return
% from the algorithm are a set of cell arrays where
% p = cell(5,1), flag = cell(5,1), hist = cell(5,1), dt = cell(5,1)
% and each cell contains the result from one algorithm.
% p{1} is the vector computed from the 'power' algorithm
% p{2} is the vector computed from the 'gs' algorithm
% p{3} is the vector computed from the 'arnoldi' algorithm
% p{4} is the vector computed from the 'linsys' algorithm with bicgstab
% p{5} is the vector comptued from the 'linsys' algorithm with gmres
% the other outputs all match these indices.
%
% pagerank(A,options) specifies optional parameters
% options.c: the teleportation coefficient [double | {0.85}]
% options.tol: the stopping tolerance [double | {1e-7}]
% options.v: the personalization vector [vector | {uniform: 1/n}]
% options.maxiter maximum number of iterations [integer | {500}]
% options.verbose: extra output information [{0} | 1]
% options.x0: the initial vector [vector | {options.v}]
% options.alg: force the algorithm type
% ['gs' | 'power' | 'linsys' | 'dense' | {'arnoldi'} | ...
% 'approx' | 'eval']
%
% options.linsys_solver: a function handle for the linear solver used
% with the linsys option [fh | {@(f,v,tol,its) bicgstab(f,v,tol,its)}]
% options.arnoldi_k: use a k dimensional arnoldi basis [intger | {8}]
% options.approx_bp: boundary probability to expand [float | 1e-3]
% options.approx_boundary: when to expand on the boundary [1 | {inf}]
% options.approx_subiter: number of subiterations of power iterations
% [integer | {5}]
%
% Example:
% load cs-stanford;
% p = pagerank(A);
% p = pagerank(A,struct('alg','linsys',...
% 'linsys_solver',@(f,v,tol,its) gmres(f,v,[],tol, its)));
% pagerank(A,struct('alg','eval'));
%
% pagerank.m
% David Gleich
%
%
% 21 February 2006
% -- added approximate PageRank
%
% Revision 1.10
% 28 January 2006
% -- added different computational modes and timing information
%
% Revision 1.00
% 19 Octoboer 2005
%
%
% The driver does mainly parameter checking, then sends things off to one
% of the computational routines.
%
[m n] = size(A);
if (m ~= n)
error('pagerank:invalidParameter', 'the matrix A must be square');
end;
options = struct('tol', 1e-7, 'maxiter', 500, 'v', ones(n,1)./n, ...
'c', 0.85, 'verbose', 0, 'alg', 'arnoldi', ...
'linsys_solver', @(f,v,tol,its) bicgstab(f,v,tol,its), ...
'arnoldi_k', 8, 'approx_bp', 1e-3, 'approx_boundary', inf,...
'approx_subiter', 5);
if (nargin > 1)
options = merge_structs(optionsu, options);
end;
if (size(options.v) ~= size(A,1))
error('pagerank:invalidParameter', ...
'the vector v must have the same size as A');
end;
if (~issparse(A))
A = sparse(A);
end;
% normalize the matrix
P = normout(A);
switch (options.alg)
case 'dense'
[x flag hist dt] = pagerank_dense(P, options);
case 'linsys'
[x flag hist dt] = pagerank_linsys(P, options);
case 'gs'
[x flag hist dt] = pagerank_gs(P, options);
case 'power'
[x flag hist dt] = pagerank_power(P, options);
case 'arnoldi'
[x flag hist dt] = pagerank_arnoldi(P, options);
case 'approx'
[x flag hist dt] = pagerank_approx(P, options);
case 'eval'
[x flag hist dt] = pagerank_eval(P, options);
otherwise
error('pagerank:invalidParameter', ...
'invalid computation mode specified.');
end;
% ===================================
% pagerank_linsys
% ===================================
function [x flag hist dt] = pagerank_linsys(P, options)
if (options.verbose > 0)
fprintf('linear system computation...\n');
end;
tol = options.tol;
v = options.v;
maxiter = options.maxiter;
c = options.c;
solver = options.linsys_solver;
% transpose P (see pagerank_linsys_mult docs)
P = P';
f = @(x,varargin) pagerank_linsys_mult(x,P,c,length(varargin));
tic;
[x flag ignore1 ignore2 hist] = solver(f,v,tol,maxiter);
dt = toc;
% renormalize the vector to have norm 1
x = x./norm(x,1);
function y = pagerank_linsys_mult(x,P,c,tflag)
% compute the matrix vector product for the linear system. This function
% includes the transpose flag (tflag > 0) to indicate a transpose multiply.
% Because many of the algorithms just use A*x (and not A'*x) the matrix P
% should have already been transposed.
if (tflag > 0)
%y = x - c*P'*x;
y = x - c*spmatvec_transmult(P,x);
else
%y = x - c*P*x;
y = x - c*spmatvec_mult(P,x);
end;
% ===================================
% pagerank_dense
% ===================================
function [x flag hist dt] = pagerank_dense(P, options)
% solve as a dense linear system
if (options.verbose > 0)
fprintf('dense computation...\n');
end;
v = options.v;
c = options.c;
n = size(P,1);
P = eye(n) - c*full(P)';
tic;
x = P \ v;
dt = toc;
hist = norm(P*x - v,1);
flag = 0;
% renormalize the vector to have norm 1
x = x./norm(x,1);
% ===================================
% pagerank_gs
% ===================================
function [x flag hist dt] = pagerank_gs(P, options)
% use gauss-seidel computation
if (options.verbose > 0)
fprintf('gauss-seidel computation...\n');
end;
tol = options.tol;
v = options.v;
maxiter = options.maxiter;
c = options.c;
x = v;
if (isfield(options, 'x0'))
x = options.x0;
else
% this is dumb, but we need to make sure
% we actually get x it's own memory...
% right now, Matlab just has a ``shadow copy''
x(1) = x(1)-1.0;
x(1) = x(1)+1.0;
end;
delta = 1;
iter = 0;
P = -c*P;
hist = zeros(maxiter,1);
dt = 0;
while (delta > tol && iter < maxiter)
tic;
xold = pagerank_gs_mult(P,x,(1-c)*v);
dt = dt + toc;
delta = norm(x - xold,1);
hist(iter+1) = delta;
if (options.verbose > 0)
fprintf('iter=%d; delta=%f\n', iter, delta);
end;
iter = iter + 1;
end;
% resize hist
hist = hist(1:iter);
% renormalize the vector to have norm 1
x = x./norm(x,1);
% default is convergence
flag = 0;
if (delta > tol && iter == maxiter)
warning('pagerank:didNotConverge', ...
'The PageRank algorithm did not converge after %i iterations', ...
maxiter);
flag = 1;
end;
% ===================================
% pagerank_power
% ===================================
function [x flag hist dt] = pagerank_power(P, options)
% use the power iteration algorithm
if (options.verbose > 0)
fprintf('power iteration computation...\n');
end;
tol = options.tol;
v = options.v;
maxiter = options.maxiter;
c = options.c;
x = v;
if (isfield(options, 'x0'))
x = options.x0;
end;
hist = zeros(maxiter,1);
delta = 1;
iter = 0;
dt = 0;
while (delta > tol && iter < maxiter)
tic;
y =c* spmatvec_transmult(P,x);
w = 1 - norm(y,1);
y = y + w*v;
dt = dt + toc;
delta = norm(x - y,1);
hist(iter+1) = delta;
tic;
x = y;
dt = dt + toc;
if (options.verbose > 0)
fprintf('iter=%d; delta=%f\n', iter, delta);
end;
iter = iter + 1;
end;
% resize hist
hist = hist(1:iter);
flag = 0;
if (delta > tol && iter == maxiter)
warning('pagerank:didNotConverge', ...
'The PageRank algorithm did not converge after %i iterations', ...
maxiter);
flag = 1;
end;
% ===================================
% pagerank_arnoldi
% ===================================
function [x flag hist dt] = pagerank_arnoldi(P, options)
% use the power iteration algorithm
if (options.verbose > 0)
fprintf('arnoldi method computation...\n');
end;
tol = options.tol;
v = options.v;
maxiter = options.maxiter;
c = options.c;
k = options.arnoldi_k;
x = v;
if (isfield(options, 'x0'))
x = options.x0;
end;
hist = zeros(maxiter,1);
d = dangling(P);
d = double(d);
P = P';
%f = @(x) pagerank_arnoldi_mult(x,P,c,d,v);
f = @(x) pagerank_mult(x,P,c,d,v);
iter = 0;
dt = 0;
delta = 1;
while (delta > tol && iter < maxiter)
tic;
[Q H] = pagerank_arnoldi_fact(f,x,k);
[u,s,v]=svd(H-[speye(k);zeros(1,k)]);
x=Q(:,1:k)*v(:,k);
dt = dt + toc;
% for statistics purposes only
delta=norm(f(x)-x,1)/norm(x,1);
hist(iter+1) = delta;
if (options.verbose > 0)
fprintf('iter=%d; delta=%f\n', iter, delta);
end;
iter = iter + 1;
end;
% ensure correct normalization
x = sign(sum(x))*x;
x = x/norm(x,1);
% resize hist
hist = hist(1:iter);
flag = 0;
if (delta > tol && iter == maxiter)
warning('pagerank:didNotConverge', ...
'The PageRank algorithm did not converge after %i iterations', ...
maxiter);
flag = 1;
end;
function [V,H] = pagerank_arnoldi_fact(A,V,k)
% [Q,H] = ARNOLDI7(A,Q0,K,c,d,e,v)
%
% ARNOLDI: Reduce an n x n matrix A to upper Hessenberg form.
% [Q,H] = ARNOLDI(A,Q0,K) computes (k+1) x k upper
% Hessenberg matrix H and n x k matrix Q with orthonormal
% columns and Q(:,1) = Q0/NORM(Q0), such that
% Q(:,1:k+1)'*A*Q(:,1:k) = H.
%
% A can also be a function_handle to return A*x
%
%
% Written by Chen Grief
% modified by David Gleich
%
V(:,1) = V(:,1)/norm(V(:,1));
if (~isa(A,'function_handle'))
f = @(x) A*x;
A = f;
end;
w = A(V(:,1));
alpha=V(:,1)'*w;
H(1,1)=alpha;
f(:,1)=w-V(:,1)*alpha;
for j=1:k-1
beta=norm(f(:,j));
V(:,j+1)=f(:,j)/beta;
ejt=[zeros(1,j-1) beta];
Hhat=[H; ejt];
w=A(V(:,j+1));
h=V(:,1:j+1)'*w;
f(:,j+1)=w-V(:,1:j+1)*h;
H=[Hhat h];
end
% Extend Arnoldi factorization
beta=norm(f(:,k));
V(:,k+1) = f(:,k)/beta;
ejt=[zeros(1,k-1) beta];
H=[H ;ejt];
% ===================================
% pagerank_approx
% ===================================
function [x flag hist dt] = pagerank_approx(A, options)
% use the power iteration algorithm
if (options.verbose > 0)
fprintf('approximate computation...\n');
end;
tol = options.tol;
v = options.v;
maxiter = options.maxiter;
c = options.c;
bp = options.approx_bp;
subiter = options.approx_subiter;
boundary = options.approx_boundary;
n = size(A,1);
%x = v;
%if (isfield(options, 'x0'))
% x = options.x0;
%end;
if (length(find(v)) ~= n)
global_pr = 0;
else
global_pr = 1;
error('pagerank:invalidParameter',...
'approximation computations are not implemented for global pagerank yet');
end;
hist = zeros(maxiter,1);
delta = 1;
iter = 0;
dt = 0;
% set the initial set of seed pages
if (global_pr)
if (isfield(options, 'x0'))
% the seed pages come from the x0 vector if provided
p = find(options.x0);
x = x0(p);
else
% the seed pages come from the x0 vector (otherwise, choose random)
p = unique(ceil(rand(250,1)*size(P,1)));
x = ones(length(p),1)./length(p);
end;
else
% the seed pages come from the x0 vector
p = find(v);
x = ones(length(p),1)./length(p);
v = v(p);
end;
local = [];
active = p;
frontier = p;
tic;
while (iter <= maxiter && delta > tol)
% expand all pages
if (boundary == 1)
% if we are running the boundary algorithm...
[ignore sp] = sort(-x);
cs = cumsum(x(sp));
spactive = active(sp);
allexpand_ind = cs < (1-bp);
% actually, we need to add the first 0 after the last 1 in
% allexpand_ind because we need cumsum to be larger than 1-bp
allexpand_ind(min(find(allexpand_ind == 0))) = ~0;
allexpand = spactive(allexpand_ind);
toexpand = setdiff(allexpand,local);
else
%
% otherwise, just expand all pages with a sufficient tolerance
%
allexpand = active(x > bp);
toexpand = setdiff(allexpand,local);
end;
if (length(toexpand) > 0)
xp = zeros(n,1);
xp([local frontier]) = x;
local = [local toexpand];
frontier = setdiff(find(sum(A(local,:),1)), local);
active = [local frontier];
x = xp(local);
else
xp = zeros(n,1);
xp([local frontier]) = x;
x = xp(local);
end;
Lp = A(local,active);
outdegree = full(sum(Lp,2));
outdegree = [outdegree; zeros(length(frontier),1)];
siter = 0;
L = [Lp; sparse(length(frontier),length(active))];
x2 = [x; xp(frontier)];
while (siter < subiter)
y = full(c*L'*(invzero(outdegree).*x2));
omega = 1 - norm(y,1);
% the ordering of local is preseved, so these are always the
% correct vertices
y(1:length(p)) = y(1:length(p)) + omega*v;
x2 = y;
siter = siter+1;
end;
x2 = [x; xp(frontier)];
delta = norm(y-x2,1);
hist(iter+1) = delta;
if (options.verbose > 0)
fprintf('iter=%02i; delta=%0.03e expand=%i\n', iter, delta, length(toexpand));
end
x = y;
iter = iter + 1;
end;
dt = toc;
% resize hist
hist = hist(1:iter);
xpartial = x;
x = zeros(n,1);
x([local frontier]) = xpartial;
flag = 0;
if (delta > tol && iter == maxiter)
warning('pagerank:didNotConverge', ...
'The PageRank algorithm did not converge after %i iterations', ...
maxiter);
flag = 1;
end;
no% ===================================
% pagerank_eval
% ===================================
function [x flag hist dt] = pagerank_eval(P,options)
algs = {'power', 'gs', 'arnoldi', 'linsys', 'linsys'};
extra_opts = {struct(''), struct(''), struct(''), struct(''), ...
struct('linsys_solver',@(f,v,tol,its) gmres(f,v,8,tol, its))};
names = {'power', 'gs', 'arnoldi8', 'bicgstab', 'gmres8'};
v = options.v;
c = options.c;
x = cell(5,1);
flag = cell(5,1);
hist = cell(5,1);
dt = cell(5,1);
web('text://<html><body>Generating PageRank report...</body></html>','-noaddressbox');
htmlend = '</body></html>';
s = {};
s{1} = '<html><head><title>PageRank runtime report</title></head><body><h1>PageRank Report</h1>';
stemp = s;
stemp{end+1} = '<p>Generating graph statistics...</p>';
stemp{end+1} = htmlend;
A = spones(P);
d = dangling(P);
npages = size(P,1);
nedges = nnz(P);
ndangling = sum(d);
maxindeg = full(max(sum(A,1)));
maxoutdeg = full(max(sum(A,2)));
ncomp = components(A);
s{end+1} = '<h2>Graph statistics</h2>';
s{end+1} = '<table border="0" cellspacing="4">';
s{end+1} = sprintf('<tr><td style="font-weight: bold">%s:</td><td>%i</td></tr>', ...
'Number of pages', npages);
s{end+1} = sprintf('<tr><td style="font-weight: bold">%s:</td><td>%i</td></tr>', ...
'Number of edges', nedges);
s{end+1} = sprintf('<tr><td style="font-weight: bold">%s:</td><td>%i</td></tr>', ...
'Number of dangling nodes', ndangling);
s{end+1} = sprintf('<tr><td style="font-weight: bold">%s:</td><td>%i</td></tr>', ...
'Max in-degree', maxindeg);
s{end+1} = sprintf('<tr><td style="font-weight: bold">%s:</td><td>%i</td></tr>', ...
'Max out-degree', maxoutdeg);
s{end+1} = sprintf('<tr><td style="font-weight: bold">%s:</td><td>%i</td></tr>', ...
'Number of strong components:', ncomp);
s{end+1} = '</table>';
sOut = [stemp{:}];
web(['text://' sOut],'-noaddressbox');
s{end+1} = '<h2>Algorithm performance</h2>';
s{end+1} = '<table border="0">';
s{end+1} = sprintf('<tr><td style="text-align: right">%s</td><td>%0.3f</td></tr>', ...
'c = ', c);
s{end+1} = sprintf('<tr><td style="text-align: right">%s</td><td>%2.2e</td></tr>', ...
'tol = ', options.tol);
s{end+1} = sprintf('<tr><td style="text-align: right">%s</td><td>%i</td></tr>', ...
'maxiter = ', options.maxiter);
s{end+1} = '</table>';
s{end+1} = '<table border="0">';
s{end+1} = ['<tr style="text-align: left">' ...
'<th style="border-bottom:solid 1px">Algorithm</th>' ...
'<th style="border-bottom:solid 1px">Time</th>' ...
'<th style="border-bottom:solid 1px">Iterations</th>' ...
'<th style="border-bottom:solid 1px">Error</th></tr>'];
for (ii=1:length(algs))
alg = algs{ii};
extra_opt = extra_opts{ii};
name = names{ii};
stemp = s;
stemp{end+1} = '</table>';
stemp{end+1} = sprintf('<p>Solving for PageRank with %s...</p>', char(name));
stemp{end+1} = htmlend;
sOut = [stemp{:}];
web(['text://' sOut],'-noaddressbox');
extra_opt = merge_structs(struct('alg',char(alg)),extra_opt);
[pi flagi histi dti] = pagerank(P, merge_structs(extra_opt,options));
p{ii} = pi;
flag{ii} = flagi;
hist{ii} = histi;
dt{ii} = dti;
err = norm(pi - c*(pi'*P)' - c*(d'*pi)*v - (1-c)*v,1);
if (mod(ii,2) == 0)
s{end+1} = sprintf('<tr style="background-color: #cccccc"><td>%s</td><td>%.2f</td><td>%i</td><td>%2.2e</td></tr>',...
char(name), dti, length(histi), err);
else
s{end+1} = sprintf('<tr><td>%s</td><td>%.2f</td><td>%i</td><td>%2.2e</td></tr>',...
char(name), dti, length(histi), err);
end;
end
s{end+1} = '</table>';
s{end+1} = htmlend;
sOut = [s{:}];
web(['text://' sOut],'-noaddressbox');
s{end+1} = sprintf('<tr><td>%s</td><td></td><td></td><td></td></tr>',char(name));
%
% plot the time histogram
%
figure(1);
close(1);
figure(1);
dts = cell2mat(dt);
flags = cell2mat(flag);
h2 = bar(dts.*(flags==0));
set(h2,'FaceColor',[1 1 1]);
set(h2,'LineWidth',2.0);
set(gca,'XTick', 1:length(algs));
set(gca,'XTickLabel',names);
ylabel('time (sec)');
%
% plot the history results
%
figure(2);
close(2);
figure(2);
lso = get(0,'DefaultAxesLineStyleOrder');
lsc = get(0,'DefaultAxesColorOrder');
lso = {'o-', 'x:', '+-.', 's--', 'd-'};
nlso = length(lso);
curlso = 0;
nlsc = length(lsc);
curlsc = 0;
for ii=1:length(algs)
histi = hist{ii};
%legendname = fn{ii};
%line(1:length(mrval.hist), mrval.hist);
semilogy(1:length(histi),histi,...
lso{mod(curlso,nlso)+1}, ...
'Color',lsc(mod(curlsc,nlsc)+1,:),...
'MarkerSize',3);
hold on;
curlso = curlso+1;
curlsc = curlsc+1;
end;
title('PageRank algorithm convergence (WARNING: DIFFERENT Y-SCALES)');
xlabel('iteration')';
ylabel('convergence measure');
legend(names{:});
function S = merge_structs(A, B)
% MERGE_STRUCTS Merge two structures.
%
% S = merge_structs(A, B) makes the structure S have all the fields from A
% and B. Conflicts are resolved by using the value in A.
%
%
% merge_structs.m
% David Gleich
%
% Revision 1.00
% 19 Octoboer 2005
%
S = A;
fn = fieldnames(B);
for ii = 1:length(fn)
if (~isfield(A, fn{ii}))
S.(fn{ii}) = B.(fn{ii});
end;
end;
function P = normout(A)
% NORMOUT Normalize the outdegrees of the matrix A.
%
% P = normout(A)
%
% P has the same non-zero structure as A, but is normalized such that the
% sum of each row is 1, assuming that A has non-negative entries.
%
%
% normout.m
% David Gleich
%
% Revision 1.00
% 19 Octoboer 2005
%
% compute the row-sums/degrees
d = full(sum(A,2));
% invert the non-zeros in the data
id = invzero(d);
% scale the rows of the matrix
P = diag(sparse(id))*A;
function v = invzero(v)
% INVZERO Compute the inverse elements of a vector with zero entries.
%
% iv = invzero(v)
%
% iv is 1./v except where v = 0, in which case it is 0.
%
%
% invzero.m
% David Gleich
%
% Revision 1.00
% 19 Octoboer 2005
%
% sparse input are easy to handle
if (issparse(v))
[m n] = size(v);
[i j val] = find(v);
val = 1./val;
v = sparse(i,j,val,m,n);
return;
end;
% so are dense input
% compute the 0 mask
zm = abs(v) > eps(1);
% invert all non-zeros
v(zm) = 1./v(zm);
function dmask = dangling(A)
% DANGLING Compute the indicator vector for dangling links for webgraph A
%
% d = dangling(A)
%
d = full(sum(A,2));
dmask = d == 0;
function [k,sizes]=components(A)
% based on components.m from (MESHPART Toolkit)
% which had
% John Gilbert, Xerox PARC, 8 June 1992.
% Copyright (c) 1990-1996 by Xerox Corporation. All rights reserved.
% HELP COPYRIGHT for complete copyright and licensing notice
[p,p,r,r] = dmperm(A|speye(size(A)));
sizes = diff(r);
k = length(sizes);
