33 changed files with 479087 additions and 2617 deletions
@ -1,8 +0,0 @@ |
|||||
The datasets on this folder where downloaded from the website of the Stanford Network Analysis Project (SNAP), found [here](https://snap.stanford.edu/data/). |
|
||||
|
|
||||
More details about the datasets can be found in the table bellow. |
|
||||
|
|
||||
| Dataset directory | Description | Nodes | Edges | URL link | |
|
||||
| ----------- | ----------- | ----------- | ----------- | ----------- | |
|
||||
| web-Google | "web-Google" | 875,713 | 5,105,039 | [link](https://snap.stanford.edu/data/web-Google.html) | |
|
||||
| wiki-Talk | "wiki-Talk" | 2,394,385 | 5,021,410 | [link](https://snap.stanford.edu/data/wiki-Talk.html) | |
|
||||
Binary file not shown.
Binary file not shown.
File diff suppressed because it is too large
File diff suppressed because it is too large
@ -1,162 +0,0 @@ |
|||||
|
|
||||
//ALGORITHM 6
|
|
||||
//p_1 is period 1, p_2 is period 2 in no. of iterations
|
|
||||
float** filterMARP(int** A, float* x, float e, int max_iterations, float a, int p_1, int p_2, int size){ |
|
||||
|
|
||||
float* x_new, x_c, y; |
|
||||
int count = 0; |
|
||||
float norm = 1.0; |
|
||||
int* N, C, N_new, C_new; |
|
||||
int** Ann, Acn; |
|
||||
int i; |
|
||||
/*********** ...malloc... *************/ |
|
||||
x_new = malloc(size*sizeof(float)); |
|
||||
|
|
||||
x_c = malloc(size*sizeof(float)); |
|
||||
y = malloc(size*sizeof(float)); |
|
||||
N = malloc(size*sizeof(int)); |
|
||||
C = malloc(size*sizeof(int)); |
|
||||
N_new = malloc(size*sizeof(int)); |
|
||||
C_new = malloc(size*sizeof(int)); |
|
||||
Ann = malloc(size*sizeof(int *)); |
|
||||
for(i=0; i<size; i++){ |
|
||||
Ann[i] = malloc(size*sizeof(int)); |
|
||||
} |
|
||||
Acn = malloc(size*sizeof(int *)); |
|
||||
for(i=0; i<size; i++){ |
|
||||
Acn[i] = malloc(size*sizeof(int)); |
|
||||
} |
|
||||
|
|
||||
/******************* start iterations ********************/ |
|
||||
while(count<max_iterations && norm<e){ |
|
||||
if(count == 0){ |
|
||||
gauss_Seidel(x, x_new, A, NULL, NULL, a, size); |
|
||||
} |
|
||||
else{ |
|
||||
gauss_Seidel(x, x_new, Ann, Acn, x_c, a, size); |
|
||||
} |
|
||||
|
|
||||
count++; |
|
||||
if(count%p_1==0){ |
|
||||
N = N_new; //might not be needed
|
|
||||
C = C_new; |
|
||||
detect_Converged(N_new, C_new, x_new, x, e, size); //sinthiki stin prwti selida toy prwtou kefalaiou
|
|
||||
filterA(A, Ann, Acn, C_new, N_new, size); //
|
|
||||
y = compy(Acn, x, size);//Acn*x ;
|
|
||||
x_c = filterx(x_new, C, size); |
|
||||
} |
|
||||
if(count%p_2==0){ |
|
||||
norm = comp_norm(x_new, A); //||A*x - x||
|
|
||||
} |
|
||||
x = x_new; |
|
||||
} |
|
||||
return x_new; |
|
||||
} |
|
||||
|
|
||||
void detect_Converged(int* N_new, int* C_new, float* x_new, float* x, float e, int size){ |
|
||||
int i; |
|
||||
for( i=0; i<size; i++){ |
|
||||
if((norm(x_new[i] - x[i])/norm(x[i]))<e){ // 10 ^ -3
|
|
||||
//converged
|
|
||||
N_new[i] = 0; |
|
||||
C_new[i] = 1; |
|
||||
} |
|
||||
else{ |
|
||||
N_new[i] = 1; |
|
||||
C_new[i] = 0; |
|
||||
} |
|
||||
} |
|
||||
return; |
|
||||
} |
|
||||
|
|
||||
float* compy(int** Acn, float* x, int size){ //Acn*x pollaplasiasmos
|
|
||||
float* y; |
|
||||
y = malloc(size*sizeof(float)); |
|
||||
int i,j; |
|
||||
for(i=0; i<size; i++){ |
|
||||
for(j=0; j<size; j++){ |
|
||||
y[i] += Acn[i][j]*x[j]; |
|
||||
} |
|
||||
} |
|
||||
return y; |
|
||||
} |
|
||||
|
|
||||
void filterA(int** A, int** Ann, int** Acn, int* C_new, int* N_new, int size){ |
|
||||
int i,j; |
|
||||
for(i=0; i<size; i++){ |
|
||||
if(C_new[i] == 1){ |
|
||||
for(j = 0; j<size;j++){ |
|
||||
Ann[i][j] = 0; //sxesi 10 sto pdf
|
|
||||
} |
|
||||
} |
|
||||
else{ |
|
||||
for(j = 0; j<size;j++){ |
|
||||
Ann[i][j] = A[i][j]; //sxesi 10 sto pdf
|
|
||||
} |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
|
|
||||
//Acn part: pages that have converged to pages that haven't
|
|
||||
// assuming A is column - based transition matrix (since it was transposed)
|
|
||||
|
|
||||
for(i=0; i<size; i++){ |
|
||||
for(j=0; j<size ; j++){ |
|
||||
if(C_new[j] == 1 && N_new[i] == 1){ |
|
||||
Acn[i][j] = A[i][j]; |
|
||||
} |
|
||||
else{ |
|
||||
Acn[i][j] = 0; |
|
||||
} |
|
||||
} |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
float* filterx(float* xnew, int* C_new, int size){ |
|
||||
int i,j; |
|
||||
float *x_c; |
|
||||
x_c = malloc(size*sizeof(float)); |
|
||||
for(i=0; i<size; i++){ |
|
||||
if(C_new[i] == 1){ |
|
||||
x_c[i] = x[i]; |
|
||||
} |
|
||||
else{ |
|
||||
x_c[i] = 0; |
|
||||
} |
|
||||
} |
|
||||
return x_c; |
|
||||
} |
|
||||
|
|
||||
void gauss_Seidel(float* x, float* x_new, int** Ann, int** Acn, float* x_c, float* y, float a, int size){ |
|
||||
int i, j; |
|
||||
float sum1 = 0; |
|
||||
float sum2 = 0; |
|
||||
x_new = x; //giati ston algorithmo leei most updated values na xrisimopoiountai gia to sum1
|
|
||||
if(Acn == NULL && x_c == NULL){ |
|
||||
|
|
||||
for(i=0; i<size; i++){ |
|
||||
for(j=0;j<i;j++){ |
|
||||
sum1 += Ann[i][j]*x_new[j]; //einai o pinakas A afou einai i prwti epanalipsi
|
|
||||
} |
|
||||
for(j=i;j<size;j++){ |
|
||||
sum2 += Ann[i][j]*x[j]; |
|
||||
} |
|
||||
x_new[i] = (1-a)+a*sum1+a*sum2; |
|
||||
} |
|
||||
} |
|
||||
else{ |
|
||||
|
|
||||
for(i=0; i<size; i++){ |
|
||||
for(j=0;j<i;j++){ |
|
||||
sum1 += Ann[i][j]*x_new[j]; |
|
||||
} |
|
||||
for(j=i;j<size;j++){ |
|
||||
sum2 += Ann[i][j]*x[j]; |
|
||||
} |
|
||||
|
|
||||
x_new[i] = (1-a)+a*sum1+a*sum2+y[i]+x_c[i]; |
|
||||
} |
|
||||
|
|
||||
} |
|
||||
|
|
||||
} |
|
||||
Binary file not shown.
@ -1,893 +0,0 @@ |
|||||
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); |
|
||||
@ -1,109 +0,0 @@ |
|||||
/*
|
|
||||
* ============================================================= |
|
||||
* pagerank_gs_mult.c Compute the matrix vector multiplication |
|
||||
* for the gauss seidel iteration in an efficient manner |
|
||||
* (that is, by overwriting the vector x in place.) |
|
||||
* |
|
||||
* David Gleich |
|
||||
* Stanford University |
|
||||
* 28 January 2006 |
|
||||
* ============================================================= |
|
||||
*/ |
|
||||
|
|
||||
#include "mex.h" |
|
||||
|
|
||||
/*
|
|
||||
* The mex function just computes one matrix-vector product. |
|
||||
*/ |
|
||||
void mexFunction(int nlhs, mxArray *plhs[], |
|
||||
int nrhs, const mxArray *prhs[]) |
|
||||
{ |
|
||||
int i, j, k; |
|
||||
int n, mrows, ncols; |
|
||||
|
|
||||
/* sparse matrix */ |
|
||||
int A_nz; |
|
||||
int *A_row, *A_col; |
|
||||
double *A_val; |
|
||||
|
|
||||
double *x, *b; |
|
||||
double *xold; |
|
||||
|
|
||||
if (nrhs != 3) |
|
||||
{ |
|
||||
mexErrMsgTxt("Three inputs required."); |
|
||||
} |
|
||||
else if (nlhs > 1) |
|
||||
{ |
|
||||
mexErrMsgTxt("Too many output arguments"); |
|
||||
} |
|
||||
|
|
||||
mrows = mxGetM(prhs[0]); |
|
||||
ncols = mxGetN(prhs[0]); |
|
||||
if (mrows != ncols || |
|
||||
!mxIsSparse(prhs[0]) || |
|
||||
!mxIsDouble(prhs[0]) || |
|
||||
mxIsComplex(prhs[0])) |
|
||||
{ |
|
||||
mexErrMsgTxt("Input must be a noncomplex square sparse matrix."); |
|
||||
} |
|
||||
|
|
||||
/* okay, the first input passes */ |
|
||||
n = mrows; |
|
||||
|
|
||||
/* The second input must be a vector. */ |
|
||||
if (mxGetM(prhs[1])*mxGetN(prhs[1]) != n || |
|
||||
mxIsSparse(prhs[1]) || !mxIsDouble(prhs[1])) |
|
||||
{ |
|
||||
mexErrMsgTxt("Invalid vector."); |
|
||||
} |
|
||||
|
|
||||
/* The third input must be a vector. */ |
|
||||
if (mxGetM(prhs[2])*mxGetN(prhs[2]) != n || |
|
||||
mxIsSparse(prhs[2]) || !mxIsDouble(prhs[2])) |
|
||||
{ |
|
||||
mexErrMsgTxt("Invalid vector."); |
|
||||
} |
|
||||
|
|
||||
/* Get the sparse matrix */ |
|
||||
A_nz = mxGetNzmax(prhs[0]); |
|
||||
A_val = mxGetPr(prhs[0]); |
|
||||
A_row = mxGetIr(prhs[0]); |
|
||||
A_col = mxGetJc(prhs[0]); |
|
||||
|
|
||||
/* Get the vector x */ |
|
||||
x = mxGetPr(prhs[1]); |
|
||||
|
|
||||
/* Get the vector b */ |
|
||||
b = mxGetPr(prhs[2]); |
|
||||
|
|
||||
/* if they request x old, then we need to copy x to xold */ |
|
||||
if (nlhs > 0) |
|
||||
{ |
|
||||
plhs[0] = mxDuplicateArray(prhs[1]); |
|
||||
} |
|
||||
|
|
||||
/* Update x in place, this means we have to iterate over columns
|
|
||||
* of the matrix A. */ |
|
||||
|
|
||||
for (i = 0; i < n; i++) |
|
||||
{ |
|
||||
/* we actually compute one iteration for the
|
|
||||
* system (I+A')x = b */ |
|
||||
double aself = 1.0; |
|
||||
double xnew = b[i]; |
|
||||
|
|
||||
|
|
||||
for (j = A_col[i]; j < A_col[i+1]; ++j) |
|
||||
{ |
|
||||
/* add to aself only if the row = i (the column) */ |
|
||||
aself += A_val[j]*(A_row[j] == i); |
|
||||
|
|
||||
/* add to xnew only if row != i */ |
|
||||
xnew -= A_val[j]*x[A_row[j]]*(A_row[j] != i); |
|
||||
} |
|
||||
|
|
||||
x[i] = xnew/aself; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
Binary file not shown.
Binary file not shown.
@ -1,128 +0,0 @@ |
|||||
/*
|
|
||||
* ============================================================= |
|
||||
* pagerank_mult.c Compute the matrix vector multiplication |
|
||||
* between the PageRank matrix and a vector |
|
||||
* |
|
||||
* David Gleich |
|
||||
* Stanford University |
|
||||
* 14 February 2006 |
|
||||
* ============================================================= |
|
||||
*/ |
|
||||
|
|
||||
#include "mex.h" |
|
||||
|
|
||||
/*
|
|
||||
* The mex function just computes one matrix-vector product. |
|
||||
* |
|
||||
* function y = pagerank_mult(x,Pt,c,d,v) |
|
||||
* y = c*Pt*x + (c*(d'*x))*v + (1-c)*sum(x)*v; |
|
||||
*/ |
|
||||
void mexFunction(int nlhs, mxArray *plhs[], |
|
||||
int nrhs, const mxArray *prhs[]) |
|
||||
{ |
|
||||
int i, j, k; |
|
||||
int n, mrows, ncols; |
|
||||
|
|
||||
/* sparse matrix */ |
|
||||
int *A_row, *A_col; |
|
||||
double *A_val; |
|
||||
|
|
||||
double *x, *d, *v; |
|
||||
double c; |
|
||||
double *y; |
|
||||
|
|
||||
double sum_x; |
|
||||
double dtx; |
|
||||
double xval; |
|
||||
|
|
||||
|
|
||||
if (nrhs != 5) |
|
||||
{ |
|
||||
mexErrMsgTxt("5 inputs required."); |
|
||||
} |
|
||||
else if (nlhs > 1) |
|
||||
{ |
|
||||
mexErrMsgTxt("Too many output arguments"); |
|
||||
} |
|
||||
|
|
||||
mrows = mxGetM(prhs[1]); |
|
||||
ncols = mxGetN(prhs[1]); |
|
||||
if (mrows != ncols || |
|
||||
!mxIsSparse(prhs[1]) || |
|
||||
!mxIsDouble(prhs[1]) || |
|
||||
mxIsComplex(prhs[1])) |
|
||||
{ |
|
||||
mexErrMsgTxt("Input must be a noncomplex square sparse matrix."); |
|
||||
} |
|
||||
|
|
||||
/* okay, the second input passes */ |
|
||||
n = mrows; |
|
||||
|
|
||||
/* The first input must be a vector. */ |
|
||||
if (mxGetM(prhs[0])*mxGetN(prhs[0]) != n || |
|
||||
mxIsSparse(prhs[0]) || !mxIsDouble(prhs[0])) |
|
||||
{ |
|
||||
mexErrMsgTxt("Invalid vector 1."); |
|
||||
} |
|
||||
|
|
||||
/* The third input must be a scalar. */ |
|
||||
if (mxGetM(prhs[2])*mxGetN(prhs[2]) != 1 || !mxIsDouble(prhs[0])) |
|
||||
{ |
|
||||
mxErrMsgTxt("Invalid scalar 3."); |
|
||||
} |
|
||||
|
|
||||
/* The fourth input must be a scalar. */ |
|
||||
if (mxGetM(prhs[3])*mxGetN(prhs[3]) != n || |
|
||||
mxIsSparse(prhs[3]) || !mxIsDouble(prhs[3])) |
|
||||
{ |
|
||||
mexErrMsgTxt("Invalid vector 4."); |
|
||||
} |
|
||||
|
|
||||
/* The fifth input must be a scalar. */ |
|
||||
if (mxGetM(prhs[4])*mxGetN(prhs[4]) != n || |
|
||||
mxIsSparse(prhs[4]) || !mxIsDouble(prhs[4])) |
|
||||
{ |
|
||||
mexErrMsgTxt("Invalid vector 5."); |
|
||||
} |
|
||||
|
|
||||
/* Get the sparse matrix */ |
|
||||
A_val = mxGetPr(prhs[1]); |
|
||||
A_row = mxGetIr(prhs[1]); |
|
||||
A_col = mxGetJc(prhs[1]); |
|
||||
|
|
||||
/* Get the vector x */ |
|
||||
x = mxGetPr(prhs[0]); |
|
||||
|
|
||||
/* Get the vector d */ |
|
||||
d = mxGetPr(prhs[3]); |
|
||||
|
|
||||
/* Get the vector v */ |
|
||||
v = mxGetPr(prhs[4]); |
|
||||
|
|
||||
c = mxGetScalar(prhs[2]); |
|
||||
|
|
||||
plhs[0] = mxCreateDoubleMatrix(n,1,mxREAL); |
|
||||
y = mxGetPr(plhs[0]); |
|
||||
|
|
||||
sum_x = 0.0; |
|
||||
dtx = 0.0; |
|
||||
|
|
||||
for (i = 0; i < n; i++) |
|
||||
{ |
|
||||
xval = x[i]; |
|
||||
sum_x += xval; |
|
||||
dtx += d[i]*xval; |
|
||||
|
|
||||
for (j = A_col[i]; j < A_col[i+1]; ++j) |
|
||||
{ |
|
||||
y[A_row[j]] += c*A_val[j]*xval; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
xval = c*dtx + (1-c)*sum_x; |
|
||||
for (i=0; i < n;i++) |
|
||||
{ |
|
||||
y[i] += xval*v[i]; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
Binary file not shown.
Binary file not shown.
@ -1,78 +0,0 @@ |
|||||
/*
|
|
||||
* ============================================================== |
|
||||
* spmatvec_mult.c Compute a sparse matrix vector multiplication |
|
||||
* |
|
||||
* David Gleich |
|
||||
* 14 February 2006 |
|
||||
* ============================================================= |
|
||||
*/ |
|
||||
|
|
||||
#include "mex.h" |
|
||||
|
|
||||
/*
|
|
||||
* The mex function just computes one matrix-vector product. |
|
||||
* |
|
||||
* function y = A*x; |
|
||||
*/ |
|
||||
void mexFunction(int nlhs, mxArray *plhs[], |
|
||||
int nrhs, const mxArray *prhs[]) |
|
||||
{ |
|
||||
int i, j, k; |
|
||||
int mrows, ncols; |
|
||||
|
|
||||
/* sparse matrix */ |
|
||||
int *A_row, *A_col; |
|
||||
double *A_val; |
|
||||
|
|
||||
double *x; |
|
||||
double *y; |
|
||||
|
|
||||
double xval; |
|
||||
|
|
||||
if (nrhs != 2) |
|
||||
{ |
|
||||
mexErrMsgTxt("2 inputs required."); |
|
||||
} |
|
||||
else if (nlhs > 1) |
|
||||
{ |
|
||||
mexErrMsgTxt("Too many output arguments"); |
|
||||
} |
|
||||
|
|
||||
mrows = mxGetM(prhs[0]); |
|
||||
ncols = mxGetN(prhs[0]); |
|
||||
if (!mxIsSparse(prhs[0]) || |
|
||||
!mxIsDouble(prhs[0]) || |
|
||||
mxIsComplex(prhs[0])) |
|
||||
{ |
|
||||
mexErrMsgTxt("Input must be a noncomplex sparse matrix."); |
|
||||
} |
|
||||
|
|
||||
/* The first input must be a vector. */ |
|
||||
if (mxGetM(prhs[1])*mxGetN(prhs[1]) != ncols || |
|
||||
mxIsSparse(prhs[1]) || !mxIsDouble(prhs[1])) |
|
||||
{ |
|
||||
mexErrMsgTxt("Invalid vector."); |
|
||||
} |
|
||||
|
|
||||
/* Get the sparse matrix */ |
|
||||
A_val = mxGetPr(prhs[0]); |
|
||||
A_row = mxGetIr(prhs[0]); |
|
||||
A_col = mxGetJc(prhs[0]); |
|
||||
|
|
||||
/* Get the vector x */ |
|
||||
x = mxGetPr(prhs[1]); |
|
||||
|
|
||||
plhs[0] = mxCreateDoubleMatrix(mrows,1,mxREAL); |
|
||||
y = mxGetPr(plhs[0]); |
|
||||
|
|
||||
for (i = 0; i < ncols; i++) |
|
||||
{ |
|
||||
xval = x[i]; |
|
||||
|
|
||||
for (j = A_col[i]; j < A_col[i+1]; ++j) |
|
||||
{ |
|
||||
y[A_row[j]] += A_val[j]*xval; |
|
||||
} |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
Binary file not shown.
Binary file not shown.
@ -1,82 +0,0 @@ |
|||||
/*
|
|
||||
* ============================================================= |
|
||||
* spmatvec_mult.c Compute a sparse matrix vector multiplication |
|
||||
* using a transposed matrix. |
|
||||
* |
|
||||
* David Gleich |
|
||||
* Stanford University |
|
||||
* 14 February 2006 |
|
||||
* ============================================================= |
|
||||
*/ |
|
||||
|
|
||||
#include "mex.h" |
|
||||
|
|
||||
/*
|
|
||||
* The mex function just computes one matrix-vector product. |
|
||||
* |
|
||||
* function y = A'*x |
|
||||
*/ |
|
||||
void mexFunction(int nlhs, mxArray *plhs[], |
|
||||
int nrhs, const mxArray *prhs[]) |
|
||||
{ |
|
||||
int i, j, k; |
|
||||
int mrows, ncols; |
|
||||
|
|
||||
/* sparse matrix */ |
|
||||
int *A_row, *A_col; |
|
||||
double *A_val; |
|
||||
|
|
||||
double *x; |
|
||||
double *y; |
|
||||
|
|
||||
double yval; |
|
||||
|
|
||||
if (nrhs != 2) |
|
||||
{ |
|
||||
mexErrMsgTxt("2 inputs required."); |
|
||||
} |
|
||||
else if (nlhs > 1) |
|
||||
{ |
|
||||
mexErrMsgTxt("Too many output arguments"); |
|
||||
} |
|
||||
|
|
||||
mrows = mxGetM(prhs[0]); |
|
||||
ncols = mxGetN(prhs[0]); |
|
||||
if (!mxIsSparse(prhs[0]) || |
|
||||
!mxIsDouble(prhs[0]) || |
|
||||
mxIsComplex(prhs[0])) |
|
||||
{ |
|
||||
mexErrMsgTxt("Input must be a noncomplex sparse matrix."); |
|
||||
} |
|
||||
|
|
||||
/* The second input must be a vector. */ |
|
||||
if (mxGetM(prhs[1])*mxGetN(prhs[1]) != mrows || |
|
||||
mxIsSparse(prhs[1]) || !mxIsDouble(prhs[1])) |
|
||||
{ |
|
||||
mexErrMsgTxt("Invalid vector 2."); |
|
||||
} |
|
||||
|
|
||||
/* Get the sparse matrix */ |
|
||||
A_val = mxGetPr(prhs[0]); |
|
||||
A_row = mxGetIr(prhs[0]); |
|
||||
A_col = mxGetJc(prhs[0]); |
|
||||
|
|
||||
/* Get the vector x */ |
|
||||
x = mxGetPr(prhs[1]); |
|
||||
|
|
||||
plhs[0] = mxCreateDoubleMatrix(ncols,1,mxREAL); |
|
||||
y = mxGetPr(plhs[0]); |
|
||||
|
|
||||
for (i = 0; i < ncols; i++) |
|
||||
{ |
|
||||
yval = 0.0; |
|
||||
|
|
||||
for (j = A_col[i]; j < A_col[i+1]; ++j) |
|
||||
{ |
|
||||
yval += A_val[j]*x[A_row[j]]; |
|
||||
} |
|
||||
|
|
||||
y[i] = yval; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
Binary file not shown.
Binary file not shown.
@ -1,37 +0,0 @@ |
|||||
SHELL := /bin/bash |
|
||||
|
|
||||
# ============================================
|
|
||||
# COMMANDS
|
|
||||
|
|
||||
CC = gcc -std=gnu99 -pthread |
|
||||
RM = rm -f |
|
||||
CFLAGS_DEBUG=-O0 -ggdb3 -Wall -I. |
|
||||
CFLAGS=-O3 -Wall -I. |
|
||||
OBJ=serial_gs_pagerank.o serial_gs_pagerank_functions.o coo_sparse_matrix.o csr_sparse_matrix.o |
|
||||
DEPS=serial_gs_pagerank_functions.h coo_sparse_matrix.h csr_sparse_matrix.h |
|
||||
|
|
||||
# ==========================================
|
|
||||
# TARGETS
|
|
||||
|
|
||||
EXECUTABLES = pagerank.out |
|
||||
|
|
||||
.PHONY: all clean |
|
||||
|
|
||||
all: $(EXECUTABLES) |
|
||||
|
|
||||
# ==========================================
|
|
||||
# DEPENDENCIES (HEADERS)
|
|
||||
|
|
||||
%.o: %.c $(DEPS) |
|
||||
$(CC) -c -o $@ $< $(CFLAGS) |
|
||||
|
|
||||
.PRECIOUS: $(EXECUTABLES) $(OBJ) |
|
||||
|
|
||||
# ==========================================
|
|
||||
# EXECUTABLE (MAIN)
|
|
||||
|
|
||||
$(EXECUTABLES): $(OBJ) |
|
||||
$(CC) -o $@ $^ $(CFLAGS) |
|
||||
|
|
||||
clean: |
|
||||
$(RM) *.o *~ $(EXECUTABLES) |
|
||||
@ -1,132 +0,0 @@ |
|||||
#include "coo_sparse_matrix.h" |
|
||||
|
|
||||
CooSparseMatrix initCooSparseMatrix() { |
|
||||
CooSparseMatrix sparseMatrix; |
|
||||
sparseMatrix.size = 0; |
|
||||
sparseMatrix.numberOfNonZeroElements = 0; |
|
||||
sparseMatrix.elements = NULL; |
|
||||
return sparseMatrix; |
|
||||
} |
|
||||
|
|
||||
void allocMemoryForCoo(CooSparseMatrix *sparseMatrix, int numberOfElements) { |
|
||||
sparseMatrix->elements = (CooSparseMatrixElement **) malloc( |
|
||||
numberOfElements * sizeof(CooSparseMatrixElement *)); |
|
||||
sparseMatrix->size = numberOfElements; |
|
||||
} |
|
||||
|
|
||||
void addElement(CooSparseMatrix *sparseMatrix, double value, int row, int column) { |
|
||||
// Checks if there is enough space allocated
|
|
||||
if (sparseMatrix->numberOfNonZeroElements == sparseMatrix->size) { |
|
||||
printf("Number of non zero elements exceeded size of matrix!\n"); |
|
||||
exit(EXIT_FAILURE); |
|
||||
} |
|
||||
|
|
||||
// Creates the new element
|
|
||||
CooSparseMatrixElement *newElement = (CooSparseMatrixElement *) malloc( |
|
||||
sizeof(CooSparseMatrixElement)); |
|
||||
newElement->value = value; |
|
||||
newElement->rowIndex = row; |
|
||||
newElement->columnIndex = column; |
|
||||
|
|
||||
// Adds the new element to the first empty (NULL) address of the matrix
|
|
||||
sparseMatrix->elements[sparseMatrix->numberOfNonZeroElements] = newElement; |
|
||||
sparseMatrix->numberOfNonZeroElements = sparseMatrix->numberOfNonZeroElements + 1; |
|
||||
} |
|
||||
|
|
||||
void transposeSparseMatrix(CooSparseMatrix *sparseMatrix) { |
|
||||
for (int i=0; i<sparseMatrix->numberOfNonZeroElements; ++i) { |
|
||||
CooSparseMatrixElement *element = sparseMatrix->elements[i]; |
|
||||
int tempRow = element->rowIndex; |
|
||||
element->rowIndex = element->columnIndex; |
|
||||
element->columnIndex = tempRow; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
/*
|
|
||||
* This function is a port of the one found here: |
|
||||
* https://github.com/scipy/scipy/blob/3b36a57/scipy/sparse/sparsetools/coo.h#L34
|
|
||||
*/ |
|
||||
void transformToCSR(CooSparseMatrix initialSparseMatrix, |
|
||||
CsrSparseMatrix *transformedSparseMatrix) { |
|
||||
// Checks if the sizes of the two matrices fit
|
|
||||
if (initialSparseMatrix.numberOfNonZeroElements > transformedSparseMatrix->size) { |
|
||||
printf("Transformed CSR matrix does not have enough space!\n"); |
|
||||
exit(EXIT_FAILURE); |
|
||||
} |
|
||||
|
|
||||
// Calculates the elements per row
|
|
||||
for (int i=0; i<initialSparseMatrix.numberOfNonZeroElements; ++i){ |
|
||||
int rowIndex = initialSparseMatrix.elements[i]->rowIndex; |
|
||||
transformedSparseMatrix->rowCumulativeIndexes[rowIndex] = |
|
||||
transformedSparseMatrix->rowCumulativeIndexes[rowIndex] + 1; |
|
||||
} |
|
||||
|
|
||||
// Cumulative sums the non zero elements per row
|
|
||||
for (int i=0, sum=0; i<transformedSparseMatrix->size+1; ++i){ |
|
||||
int temp = transformedSparseMatrix->rowCumulativeIndexes[i]; |
|
||||
transformedSparseMatrix->rowCumulativeIndexes[i] = sum; |
|
||||
sum += temp; |
|
||||
} |
|
||||
|
|
||||
// Copies the values and columns of the elements
|
|
||||
for (int i=0; i<initialSparseMatrix.numberOfNonZeroElements; ++i){ |
|
||||
int row = initialSparseMatrix.elements[i]->rowIndex; |
|
||||
int destinationIndex = transformedSparseMatrix->rowCumulativeIndexes[row]; |
|
||||
|
|
||||
transformedSparseMatrix->columnIndexes[destinationIndex] = initialSparseMatrix.elements[i]->columnIndex; |
|
||||
transformedSparseMatrix->values[destinationIndex] = initialSparseMatrix.elements[i]->value; |
|
||||
|
|
||||
transformedSparseMatrix->rowCumulativeIndexes[row]++; |
|
||||
} |
|
||||
|
|
||||
// Fixes the cumulative sum
|
|
||||
for (int i=0, last=0; i<=transformedSparseMatrix->size; i++){ |
|
||||
int temp = transformedSparseMatrix->rowCumulativeIndexes[i]; |
|
||||
transformedSparseMatrix->rowCumulativeIndexes[i] = last; |
|
||||
last = temp; |
|
||||
} |
|
||||
|
|
||||
transformedSparseMatrix->numberOfNonZeroElements = initialSparseMatrix.numberOfNonZeroElements; |
|
||||
} |
|
||||
|
|
||||
void cooSparseMatrixVectorMultiplication(CooSparseMatrix sparseMatrix, |
|
||||
double *vector, double **product, int vectorSize) { |
|
||||
// Initializes the elements of the product vector to zero
|
|
||||
for (int i=0; i<vectorSize; ++i) { |
|
||||
(*product)[i] = 0; |
|
||||
} |
|
||||
|
|
||||
CooSparseMatrixElement *element; |
|
||||
for (int i=0; i<sparseMatrix.numberOfNonZeroElements; ++i) { |
|
||||
element = sparseMatrix.elements[i]; |
|
||||
int row = element->rowIndex, column = element->columnIndex; |
|
||||
|
|
||||
if (row >= vectorSize) { |
|
||||
printf("Error at sparseMatrixVectorMultiplication. Matrix has more rows than vector!\n"); |
|
||||
printf("row = %d\n", row); |
|
||||
exit(EXIT_FAILURE); |
|
||||
} |
|
||||
|
|
||||
(*product)[row] = (*product)[row] + element->value * vector[column]; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
void destroyCooSparseMatrix(CooSparseMatrix *sparseMatrix) { |
|
||||
for (int i=0; i<sparseMatrix->numberOfNonZeroElements; ++i) { |
|
||||
free(sparseMatrix->elements[i]); |
|
||||
} |
|
||||
free(sparseMatrix->elements); |
|
||||
} |
|
||||
|
|
||||
void printCooSparseMatrix(CooSparseMatrix sparseMatrix) { |
|
||||
if (sparseMatrix.numberOfNonZeroElements == 0) { |
|
||||
return; |
|
||||
} |
|
||||
|
|
||||
CooSparseMatrixElement *element; |
|
||||
for (int i=0; i<sparseMatrix.numberOfNonZeroElements; ++i) { |
|
||||
element = sparseMatrix.elements[i]; |
|
||||
printf("[%d,%d] = %f\n", element->rowIndex, element->columnIndex, |
|
||||
element->value); |
|
||||
} |
|
||||
} |
|
||||
@ -1,60 +0,0 @@ |
|||||
#ifndef COO_SPARSE_MATRIX_H /* Include guard */ |
|
||||
#define COO_SPARSE_MATRIX_H |
|
||||
|
|
||||
/* ===== INCLUDES ===== */ |
|
||||
|
|
||||
#include <stdbool.h> |
|
||||
#include <stdlib.h> |
|
||||
#include <stdio.h> |
|
||||
#include <stdlib.h> |
|
||||
|
|
||||
#include "csr_sparse_matrix.h" |
|
||||
|
|
||||
/* ===== STRUCTURES ===== */ |
|
||||
|
|
||||
// One element of the coordinate formated sparse matrix.
|
|
||||
typedef struct cooSparseMatrixElement { |
|
||||
double value; |
|
||||
int rowIndex, columnIndex; |
|
||||
} CooSparseMatrixElement; |
|
||||
|
|
||||
// A sparse matrix in COOrdinate format (aka triplet format).
|
|
||||
typedef struct cooSparseMatrix { |
|
||||
int size, numberOfNonZeroElements; |
|
||||
CooSparseMatrixElement **elements; |
|
||||
} CooSparseMatrix; |
|
||||
|
|
||||
/* ===== FUNCTION DEFINITIONS ===== */ |
|
||||
|
|
||||
// initCooSparseMatrix creates and initializes the members of a CooSparseMatrix
|
|
||||
// structure instance.
|
|
||||
CooSparseMatrix initCooSparseMatrix(); |
|
||||
|
|
||||
//allocMemoryForCoo allocates memory for the elements of the matrix.
|
|
||||
void allocMemoryForCoo(CooSparseMatrix *sparseMatrix, int numberOfElements); |
|
||||
|
|
||||
// addElement adds an element representing the triplet passed in the arguments
|
|
||||
// to the first empty address of the space allocated for the elements.
|
|
||||
void addElement(CooSparseMatrix *sparseMatrix, double value, int row, |
|
||||
int column); |
|
||||
|
|
||||
// transposeSparseMatrix transposes the matrix.
|
|
||||
void transposeSparseMatrix(CooSparseMatrix *sparseMatrix); |
|
||||
|
|
||||
// transformToCSR transforms the sparse matrix representation format from COO
|
|
||||
// to CSR.
|
|
||||
void transformToCSR(CooSparseMatrix initialSparseMatrix, |
|
||||
CsrSparseMatrix *transformedSparseMatrix); |
|
||||
|
|
||||
// cooSparseMatrixVectorMultiplication calculates the product of a
|
|
||||
// CooSparseMatrix and a vector.
|
|
||||
void cooSparseMatrixVectorMultiplication(CooSparseMatrix sparseMatrix, |
|
||||
double *vector, double **product, int vectorSize); |
|
||||
|
|
||||
// destroyCooSparseMatrix frees all space used by the CooSparseMatrix.
|
|
||||
void destroyCooSparseMatrix(CooSparseMatrix *sparseMatrix); |
|
||||
|
|
||||
// printCooSparseMatrix prints the values of a CooSparseMatrix.
|
|
||||
void printCooSparseMatrix(CooSparseMatrix sparseMatrix); |
|
||||
|
|
||||
#endif // COO_SPARSE_MATRIX_H
|
|
||||
@ -1,92 +0,0 @@ |
|||||
#include "csr_sparse_matrix.h" |
|
||||
|
|
||||
CsrSparseMatrix initCsrSparseMatrix() { |
|
||||
CsrSparseMatrix sparseMatrix; |
|
||||
sparseMatrix.size = 0; |
|
||||
sparseMatrix.numberOfNonZeroElements = 0; |
|
||||
|
|
||||
sparseMatrix.values = NULL; |
|
||||
sparseMatrix.columnIndexes = NULL; |
|
||||
sparseMatrix.rowCumulativeIndexes = NULL; |
|
||||
return sparseMatrix; |
|
||||
} |
|
||||
|
|
||||
void allocMemoryForCsr(CsrSparseMatrix *sparseMatrix, int numberOfElements) { |
|
||||
sparseMatrix->values = (double *) malloc(numberOfElements * sizeof(double)); |
|
||||
sparseMatrix->columnIndexes = (int *) malloc( |
|
||||
numberOfElements * sizeof(int)); |
|
||||
sparseMatrix->rowCumulativeIndexes = (int *) malloc( |
|
||||
(numberOfElements + 1) * sizeof(int)); |
|
||||
|
|
||||
for (int i=0; i<numberOfElements+1; ++i) { |
|
||||
sparseMatrix->rowCumulativeIndexes[i] = 0; |
|
||||
} |
|
||||
sparseMatrix->size = numberOfElements; |
|
||||
} |
|
||||
|
|
||||
void zeroOutRow(CsrSparseMatrix *sparseMatrix, int row) { |
|
||||
// Gets start and end indexes of the row's elements
|
|
||||
int startIndex = sparseMatrix->rowCumulativeIndexes[row], |
|
||||
endIndex = sparseMatrix->rowCumulativeIndexes[row+1]; |
|
||||
for (int i=startIndex; i<endIndex; ++i) { |
|
||||
sparseMatrix->values[i] = 0; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
void zeroOutColumn(CsrSparseMatrix *sparseMatrix, int column) { |
|
||||
for (int i=0; i<sparseMatrix->numberOfNonZeroElements; ++i){ |
|
||||
if(sparseMatrix->columnIndexes[i] == column){ |
|
||||
sparseMatrix->values[i] = 0; |
|
||||
} |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
void csrSparseMatrixVectorMultiplication(CsrSparseMatrix sparseMatrix, |
|
||||
double *vector, double **product, int vectorSize) { |
|
||||
// Initializes the elements of the product vector to zero
|
|
||||
for (int i=0; i<vectorSize; ++i) { |
|
||||
(*product)[i] = 0; |
|
||||
} |
|
||||
|
|
||||
for (int i=0; i<sparseMatrix.size; ++i) { |
|
||||
// Gets start and end indexes of this row's elements
|
|
||||
int startIndex = sparseMatrix.rowCumulativeIndexes[i], |
|
||||
endIndex = sparseMatrix.rowCumulativeIndexes[i+1]; |
|
||||
|
|
||||
if (startIndex == endIndex) { |
|
||||
// This row has no elements
|
|
||||
continue; |
|
||||
} |
|
||||
|
|
||||
double sum = 0; |
|
||||
for(int j=startIndex; j<endIndex; ++j){ |
|
||||
int elementColumn = sparseMatrix.columnIndexes[j]; |
|
||||
sum += sparseMatrix.values[j] * vector[elementColumn]; |
|
||||
} |
|
||||
|
|
||||
(*product)[i] = sum; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
void destroyCsrSparseMatrix(CsrSparseMatrix *sparseMatrix) { |
|
||||
free(sparseMatrix->values); |
|
||||
free(sparseMatrix->rowCumulativeIndexes); |
|
||||
free(sparseMatrix->columnIndexes); |
|
||||
} |
|
||||
|
|
||||
void printCsrSparseMatrix(CsrSparseMatrix sparseMatrix) { |
|
||||
if (sparseMatrix.size == 0) { |
|
||||
return; |
|
||||
} |
|
||||
|
|
||||
for (int i=0; i<sparseMatrix.size; ++i){ |
|
||||
int startIndex = sparseMatrix.rowCumulativeIndexes[i], |
|
||||
endIndex = sparseMatrix.rowCumulativeIndexes[i+1]; |
|
||||
for(int j=startIndex; j<endIndex; ++j){ |
|
||||
printf("Row [%d] has [%d] nz elements: \n at column[%d] is value = %f \n", |
|
||||
i, endIndex-startIndex, |
|
||||
sparseMatrix.columnIndexes[j], |
|
||||
sparseMatrix.values[j]); |
|
||||
} |
|
||||
} |
|
||||
} |
|
||||
@ -1,47 +0,0 @@ |
|||||
#ifndef CSR_SPARSE_MATRIX_H /* Include guard */ |
|
||||
#define CSR_SPARSE_MATRIX_H |
|
||||
|
|
||||
/* ===== INCLUDES ===== */ |
|
||||
|
|
||||
#include <stdbool.h> |
|
||||
#include <stdlib.h> |
|
||||
#include <stdio.h> |
|
||||
#include <stdlib.h> |
|
||||
|
|
||||
/* ===== STRUCTURES ===== */ |
|
||||
|
|
||||
// A sparse matrix in compressed SparseRow format.
|
|
||||
typedef struct csrSparseMatrix { |
|
||||
int size, numberOfNonZeroElements; |
|
||||
int *rowCumulativeIndexes, *columnIndexes; |
|
||||
double *values; |
|
||||
} CsrSparseMatrix; |
|
||||
|
|
||||
/* ===== FUNCTION DEFINITIONS ===== */ |
|
||||
|
|
||||
// initCsrSparseMatrix creates and initializes the members of a CsrSparseMatrix
|
|
||||
// structure instance.
|
|
||||
CsrSparseMatrix initCsrSparseMatrix(); |
|
||||
|
|
||||
// allocMemoryForCsr allocates memory for the elements of the matrix.
|
|
||||
void allocMemoryForCsr(CsrSparseMatrix *sparseMatrix, int numberOfElements); |
|
||||
|
|
||||
// zeroOutRow assigns a zero value to all the elements of a row in the matrix.
|
|
||||
void zeroOutRow(CsrSparseMatrix *sparseMatrix, int row); |
|
||||
|
|
||||
// zeroOutColumn assigns a zero value to all the elements of a column in the
|
|
||||
// matrix.
|
|
||||
void zeroOutColumn(CsrSparseMatrix *sparseMatrix, int column); |
|
||||
|
|
||||
// csrSparseMatrixVectorMultiplication calculates the product of a
|
|
||||
// CsrSparseMatrix and a vector.
|
|
||||
void csrSparseMatrixVectorMultiplication(CsrSparseMatrix sparseMatrix, |
|
||||
double *vector, double **product, int vectorSize); |
|
||||
|
|
||||
// destroyCsrSparseMatrix frees all space used by the CsrSparseMatrix.
|
|
||||
void destroyCsrSparseMatrix(CsrSparseMatrix *sparseMatrix); |
|
||||
|
|
||||
// printCsrSparseMatrix prints the values of a CsrSparseMatrix.
|
|
||||
void printCsrSparseMatrix(CsrSparseMatrix sparseMatrix); |
|
||||
|
|
||||
#endif // CSR_SPARSE_MATRIX_H
|
|
||||
Binary file not shown.
@ -1,42 +0,0 @@ |
|||||
#include <sys/time.h> |
|
||||
|
|
||||
#include "serial_gs_pagerank_functions.h" |
|
||||
|
|
||||
struct timeval startwtime, endwtime; |
|
||||
|
|
||||
int main(int argc, char **argv) { |
|
||||
CsrSparseMatrix transitionMatrix = initCsrSparseMatrix(); |
|
||||
double *pagerankVector; |
|
||||
bool convergenceStatus; |
|
||||
Parameters parameters; |
|
||||
|
|
||||
parseArguments(argc, argv, ¶meters); |
|
||||
|
|
||||
initialize(&transitionMatrix, &pagerankVector, ¶meters); |
|
||||
|
|
||||
// Starts wall-clock timer
|
|
||||
gettimeofday (&startwtime, NULL); |
|
||||
|
|
||||
int iterations = pagerank(&transitionMatrix, &pagerankVector, |
|
||||
&convergenceStatus, parameters); |
|
||||
if (parameters.verbose) { |
|
||||
printf(ANSI_COLOR_YELLOW "\n----- RESULTS -----\n" ANSI_COLOR_RESET); |
|
||||
if (convergenceStatus) { |
|
||||
printf(ANSI_COLOR_GREEN "Pagerank converged after %d iterations!\n" \ |
|
||||
ANSI_COLOR_RESET, iterations); |
|
||||
} else { |
|
||||
printf(ANSI_COLOR_RED "Pagerank did not converge after max number of" \ |
|
||||
" iterations (%d) was reached!\n" ANSI_COLOR_RESET, iterations); |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
// Stops wall-clock timer
|
|
||||
gettimeofday (&endwtime, NULL); |
|
||||
double seq_time = (double)((endwtime.tv_usec - startwtime.tv_usec)/1.0e6 + |
|
||||
endwtime.tv_sec - startwtime.tv_sec); |
|
||||
printf("%s wall clock time = %f\n","Pagerank (Gauss-Seidel method), serial implementation", |
|
||||
seq_time); |
|
||||
|
|
||||
free(pagerankVector); |
|
||||
destroyCsrSparseMatrix(&transitionMatrix); |
|
||||
} |
|
||||
@ -1,610 +0,0 @@ |
|||||
/* ===== INCLUDES ===== */ |
|
||||
|
|
||||
#include "serial_gs_pagerank_functions.h" |
|
||||
#include <pthread.h> |
|
||||
|
|
||||
/* ===== CONSTANTS ===== */ |
|
||||
|
|
||||
const char *ARGUMENT_CONVERGENCE_TOLERANCE = "-c"; |
|
||||
const char *ARGUMENT_MAX_ITERATIONS = "-m"; |
|
||||
const char *ARGUMENT_DAMPING_FACTOR = "-a"; |
|
||||
const char *ARGUMENT_VERBAL_OUTPUT = "-v"; |
|
||||
const char *ARGUMENT_OUTPUT_HISTORY = "-h"; |
|
||||
const char *ARGUMENT_OUTPUT_FILENAME = "-o"; |
|
||||
|
|
||||
const int NUMERICAL_BASE = 10; |
|
||||
char *DEFAULT_OUTPUT_FILENAME = "pagerank_output"; |
|
||||
const int FILE_READ_BUFFER_SIZE = 4096; |
|
||||
|
|
||||
const int CONVERGENCE_CHECK_ITERATION_PERIOD = 2; |
|
||||
const int SPARSITY_INCREASE_ITERATION_PERIOD = 10; |
|
||||
|
|
||||
/* ===== THREAD STUFF ====== */ |
|
||||
|
|
||||
pthread_mutex_t Q = PTHREAD_MUTEX_INITIALIZER;; |
|
||||
|
|
||||
typedef struct threadArgs{ |
|
||||
CsrSparseMatrix* transitionMatrix; |
|
||||
double* pagerankVector; |
|
||||
double* previousPagerankVector; |
|
||||
int numberOfPages; |
|
||||
double webUniformProbability; |
|
||||
double *linksFromConvergedPagesPagerankVector; |
|
||||
double *convergedPagerankVector; |
|
||||
int position; |
|
||||
double dF; |
|
||||
}threadArgs; |
|
||||
|
|
||||
/* ===== FUNCTIONS ===== */ |
|
||||
|
|
||||
int pagerank(CsrSparseMatrix *transitionMatrix, double **pagerankVector, |
|
||||
bool *convergenceStatus, Parameters parameters) { |
|
||||
// Variables declaration
|
|
||||
int iterations = 0, numberOfPages = parameters.numberOfPages; |
|
||||
double delta, *pagerankDifference, *previousPagerankVector, |
|
||||
*convergedPagerankVector, *linksFromConvergedPagesPagerankVector; |
|
||||
CooSparseMatrix linksFromConvergedPages = initCooSparseMatrix(); |
|
||||
bool *convergenceMatrix; |
|
||||
|
|
||||
|
|
||||
int threadNum = parameters.numberOfPages; |
|
||||
pthread_t* threads = (pthread_t *)malloc(threadNum*sizeof(pthread_t)); |
|
||||
|
|
||||
// Space allocation
|
|
||||
{ |
|
||||
size_t sizeofDouble = sizeof(double); |
|
||||
// pagerankDifference used to calculate delta
|
|
||||
pagerankDifference = (double *) malloc(numberOfPages * sizeofDouble); |
|
||||
// previousPagerankVector holds last iteration's pagerank vector
|
|
||||
previousPagerankVector = (double *) malloc(numberOfPages * sizeofDouble); |
|
||||
// convergedPagerankVector is the pagerank vector of converged pages only
|
|
||||
convergedPagerankVector = (double *) malloc(numberOfPages * sizeofDouble); |
|
||||
// linksFromConvergedPagesPagerankVector holds the partial sum of the
|
|
||||
// pagerank vector, that describes effect of the links from converged
|
|
||||
// pages to non converged pages
|
|
||||
linksFromConvergedPagesPagerankVector = (double *) malloc(numberOfPages * sizeofDouble); |
|
||||
// convergenceMatrix indicates which pages have converged
|
|
||||
convergenceMatrix = (bool *) malloc(numberOfPages * sizeof(bool)); |
|
||||
*convergenceStatus = false; |
|
||||
|
|
||||
// Initialization
|
|
||||
allocMemoryForCoo(&linksFromConvergedPages, transitionMatrix->numberOfNonZeroElements); |
|
||||
for (int i=0; i<numberOfPages; ++i) { |
|
||||
convergedPagerankVector[i] = 0; |
|
||||
convergenceMatrix[i] = false; |
|
||||
linksFromConvergedPagesPagerankVector[i] = 0; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
if (parameters.verbose) { |
|
||||
printf(ANSI_COLOR_YELLOW "\n----- Starting iterations -----\n" ANSI_COLOR_RESET); |
|
||||
} |
|
||||
|
|
||||
do { |
|
||||
// Stores previous pagerank vector
|
|
||||
memcpy(previousPagerankVector, *pagerankVector, numberOfPages * sizeof(double)); |
|
||||
|
|
||||
// Calculates new pagerank vector
|
|
||||
calculateNextPagerank(transitionMatrix, previousPagerankVector, |
|
||||
pagerankVector, linksFromConvergedPagesPagerankVector, |
|
||||
convergedPagerankVector, numberOfPages, |
|
||||
parameters.dampingFactor, threads, threadNum); |
|
||||
|
|
||||
if (parameters.history) { |
|
||||
// Outputs pagerank vector to file
|
|
||||
savePagerankToFile(parameters.outputFilename, iterations != 0, |
|
||||
*pagerankVector, numberOfPages, iterations); |
|
||||
} |
|
||||
|
|
||||
// Periodically checks for convergence
|
|
||||
if (!(iterations % CONVERGENCE_CHECK_ITERATION_PERIOD)) { |
|
||||
// Builds pagerank vectors difference
|
|
||||
for (int i=0; i<numberOfPages; ++i) { |
|
||||
pagerankDifference[i] = (*pagerankVector)[i] - previousPagerankVector[i]; |
|
||||
} |
|
||||
// Calculates convergence
|
|
||||
delta = vectorNorm(pagerankDifference, numberOfPages); |
|
||||
|
|
||||
if (delta < parameters.convergenceCriterion) { |
|
||||
// Converged
|
|
||||
*convergenceStatus = true; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
// Periodically increases sparsity
|
|
||||
if (iterations && !(iterations % SPARSITY_INCREASE_ITERATION_PERIOD)) { |
|
||||
bool *newlyConvergedPages = (bool *) malloc(numberOfPages * sizeof(bool)); |
|
||||
// Checks each individual page for convergence
|
|
||||
for (int i=0; i<numberOfPages; ++i) { |
|
||||
double difference = fabs((*pagerankVector)[i] - |
|
||||
previousPagerankVector[i]) / fabs(previousPagerankVector[i]); |
|
||||
|
|
||||
newlyConvergedPages[i] = false; |
|
||||
if (!convergenceMatrix[i] && difference < parameters.convergenceCriterion){ |
|
||||
// Page converged
|
|
||||
newlyConvergedPages[i] = true; |
|
||||
convergenceMatrix[i] = true; |
|
||||
convergedPagerankVector[i] = (*pagerankVector)[i]; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
for (int i=0; i<numberOfPages; ++i) { |
|
||||
// Filters newly converged pages
|
|
||||
if (newlyConvergedPages[i] == true) { |
|
||||
// Checks if this converged page has an out-link to a non converged one
|
|
||||
int rowStartIndex = transitionMatrix->rowCumulativeIndexes[i], |
|
||||
rowEndIndex = transitionMatrix->rowCumulativeIndexes[i+1]; |
|
||||
if (rowEndIndex > rowStartIndex) { |
|
||||
// This row (page) has non zero elements (out-links)
|
|
||||
for (int j=rowStartIndex; j<rowEndIndex; ++j) { |
|
||||
// Checks for links from converged pages to non converged
|
|
||||
int pageLinksTo = transitionMatrix->columnIndexes[j]; |
|
||||
if (convergenceMatrix[pageLinksTo] == false){ |
|
||||
// Link exists, adds element to the vector
|
|
||||
addElement(&linksFromConvergedPages, |
|
||||
transitionMatrix->values[j], i, pageLinksTo); |
|
||||
} |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
// Increases sparsity of the transition matrix by zeroing
|
|
||||
// out elements that correspond to converged pages
|
|
||||
zeroOutRow(transitionMatrix, i); |
|
||||
zeroOutColumn(transitionMatrix, i); |
|
||||
|
|
||||
// Builds the new linksFromConvergedPagesPagerankVector
|
|
||||
cooSparseMatrixVectorMultiplication(linksFromConvergedPages, |
|
||||
*pagerankVector, &linksFromConvergedPagesPagerankVector, |
|
||||
numberOfPages); |
|
||||
} |
|
||||
} |
|
||||
free(newlyConvergedPages); |
|
||||
} |
|
||||
|
|
||||
++iterations; |
|
||||
// Outputs information about this iteration
|
|
||||
if (iterations%2) { |
|
||||
printf(ANSI_COLOR_BLUE "Iteration %d: delta = %f\n" ANSI_COLOR_RESET, iterations, delta); |
|
||||
} else { |
|
||||
printf(ANSI_COLOR_CYAN "Iteration %d: delta = %f\n" ANSI_COLOR_RESET, iterations, delta); |
|
||||
} |
|
||||
} while (!*convergenceStatus && (parameters.maxIterations == 0 || |
|
||||
iterations < parameters.maxIterations)); |
|
||||
|
|
||||
if (!parameters.history) { |
|
||||
// Always outputs last pagerank vector to file
|
|
||||
savePagerankToFile(parameters.outputFilename, false, *pagerankVector, |
|
||||
numberOfPages, iterations); |
|
||||
} |
|
||||
|
|
||||
// Frees memory
|
|
||||
free(pagerankDifference); |
|
||||
free(previousPagerankVector); |
|
||||
free(convergedPagerankVector); |
|
||||
free(linksFromConvergedPagesPagerankVector); |
|
||||
free(convergenceMatrix); |
|
||||
destroyCooSparseMatrix(&linksFromConvergedPages); |
|
||||
|
|
||||
return iterations; |
|
||||
} |
|
||||
|
|
||||
/*
|
|
||||
* initialize allocates required memory for arrays, reads the web graph from the |
|
||||
* from the file and creates the initial transition probability distribution |
|
||||
* matrix. |
|
||||
*/ |
|
||||
void initialize(CsrSparseMatrix *transitionMatrix, |
|
||||
double **pagerankVector, Parameters *parameters) { |
|
||||
|
|
||||
// Reads web graph from file
|
|
||||
if ((*parameters).verbose) { |
|
||||
printf(ANSI_COLOR_YELLOW "----- Reading graph from file -----\n" ANSI_COLOR_RESET); |
|
||||
} |
|
||||
generateNormalizedTransitionMatrixFromFile(transitionMatrix, parameters); |
|
||||
|
|
||||
// Outputs the algorithm parameters to the console
|
|
||||
if ((*parameters).verbose) { |
|
||||
printf(ANSI_COLOR_YELLOW "\n----- Running with parameters -----\n" ANSI_COLOR_RESET\ |
|
||||
"Number of pages: %d", (*parameters).numberOfPages); |
|
||||
if (!(*parameters).maxIterations) { |
|
||||
printf("\nMaximum number of iterations: inf"); |
|
||||
} else { |
|
||||
printf("\nMaximum number of iterations: %d", (*parameters).maxIterations); |
|
||||
} |
|
||||
printf("\nConvergence criterion: %f" \ |
|
||||
"\nDamping factor: %f" \ |
|
||||
"\nGraph filename: %s\n", (*parameters).convergenceCriterion, |
|
||||
(*parameters).dampingFactor, (*parameters).graphFilename); |
|
||||
} |
|
||||
|
|
||||
// Allocates memory for the pagerank vector
|
|
||||
(*pagerankVector) = (double *) malloc((*parameters).numberOfPages * sizeof(double)); |
|
||||
double webUniformProbability = 1. / (*parameters).numberOfPages; |
|
||||
for (int i=0; i<(*parameters).numberOfPages; ++i) { |
|
||||
(*pagerankVector)[i] = webUniformProbability; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
// ==================== MATH UTILS ====================
|
|
||||
|
|
||||
/*
|
|
||||
* calculateNextPagerank calculates the product of the multiplication |
|
||||
* between a matrix and the a vector in a cheap way. |
|
||||
*/ |
|
||||
void calculateNextPagerank(CsrSparseMatrix *transitionMatrix, |
|
||||
double *previousPagerankVector, double **pagerankVector, |
|
||||
double *linksFromConvergedPagesPagerankVector, |
|
||||
double *convergedPagerankVector, int vectorSize, double dampingFactor, pthread_t *threads, int threadNum) { |
|
||||
|
|
||||
pthread_attr_t attr; |
|
||||
pthread_attr_init(&attr); |
|
||||
pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE); |
|
||||
|
|
||||
// Calculates the web uniform probability once.
|
|
||||
double webUniformProbability = 1. / vectorSize; |
|
||||
int runningThreads = 0; |
|
||||
|
|
||||
for (int i=0; i<vectorSize; ++i) { |
|
||||
pthread_mutex_lock(&Q); |
|
||||
threadArgs arg; |
|
||||
arg.position = i; |
|
||||
arg.transitionMatrix = transitionMatrix; |
|
||||
arg.pagerankVector = *pagerankVector; |
|
||||
arg.previousPagerankVector = previousPagerankVector; |
|
||||
arg.linksFromConvergedPagesPagerankVector = linksFromConvergedPagesPagerankVector; |
|
||||
arg.convergedPagerankVector = convergedPagerankVector; |
|
||||
arg.webUniformProbability = webUniformProbability; |
|
||||
arg.dF = dampingFactor; |
|
||||
if(runningThreads < threadNum){ |
|
||||
if(pthread_create(&threads[i], &attr, csrSparseMatrixVectorMultiplication_threads, (void *) &arg)){ |
|
||||
printf("Error creating thread %d", i); |
|
||||
pthread_mutex_unlock(&Q); |
|
||||
exit(-1); |
|
||||
} |
|
||||
else{ |
|
||||
++runningThreads; |
|
||||
pthread_mutex_unlock(&Q); |
|
||||
} |
|
||||
} |
|
||||
else{ |
|
||||
pthread_mutex_unlock(&Q); |
|
||||
printf("not enough threads\n"); |
|
||||
exit(-1); |
|
||||
} |
|
||||
pthread_join(threads[i], NULL); |
|
||||
if(runningThreads < threadNum){ |
|
||||
if(pthread_create(&threads[i], &attr, compPagerankVector_threads, (void *) &arg)){ |
|
||||
printf("Error creating thread %d", i); |
|
||||
pthread_mutex_unlock(&Q); |
|
||||
exit(-1); |
|
||||
} |
|
||||
else{ |
|
||||
++runningThreads; |
|
||||
pthread_mutex_unlock(&Q); |
|
||||
} |
|
||||
} |
|
||||
else{ |
|
||||
printf("You have a problem \n"); |
|
||||
exit(-1); |
|
||||
pthread_mutex_unlock(&Q); |
|
||||
} |
|
||||
|
|
||||
} |
|
||||
for(int i=0; i<vectorSize; ++i){ |
|
||||
pthread_join(threads[i], NULL); |
|
||||
} |
|
||||
free(threads); |
|
||||
|
|
||||
} |
|
||||
|
|
||||
void compPagerankVector_threads(void* arg){ |
|
||||
threadArgs* co = (threadArgs *)arg; |
|
||||
co->pagerankVector[co->position] = co->dF * co->pagerankVector[co->position]; |
|
||||
|
|
||||
|
|
||||
double normDifference = vectorNorm(co->previousPagerankVector, co->numberOfPages) - |
|
||||
vectorNorm(co->pagerankVector, co->numberOfPages); |
|
||||
|
|
||||
|
|
||||
co->pagerankVector[co->position] += normDifference * co->webUniformProbability + |
|
||||
co->linksFromConvergedPagesPagerankVector[co->position] + co->convergedPagerankVector[co->position]; |
|
||||
|
|
||||
} |
|
||||
|
|
||||
void csrSparseMatrixVectorMultiplication_threads(void* arg){ |
|
||||
threadArgs* co = (threadArgs *)arg; |
|
||||
//(CsrSparseMatrix sparseMatrix,
|
|
||||
//double *vector, double **product, int vectorSize) {
|
|
||||
// Initializes the elements of the product vector to zero
|
|
||||
//for (int i=0; i<vectorSize; ++i) {
|
|
||||
co->pagerankVector[co->position] = 0; |
|
||||
//}
|
|
||||
|
|
||||
//for (int i=0; i<sparseMatrix.size; ++i) {
|
|
||||
// Gets start and end indexes of this row's elements
|
|
||||
int startIndex = co->transitionMatrix[co->position].rowCumulativeIndexes[co->position], |
|
||||
endIndex = co->transitionMatrix[co->position].rowCumulativeIndexes[co->position+1]; |
|
||||
|
|
||||
if (startIndex == endIndex) { |
|
||||
// This row has no elements
|
|
||||
return; |
|
||||
} |
|
||||
|
|
||||
double sum = 0; |
|
||||
for(int j=startIndex; j<endIndex; ++j){ |
|
||||
int elementColumn = co->transitionMatrix[co->position].columnIndexes[j]; |
|
||||
sum += co->transitionMatrix[co->position].values[j] * co->previousPagerankVector[elementColumn]; |
|
||||
} |
|
||||
|
|
||||
co->pagerankVector[co->position] = sum; |
|
||||
//}
|
|
||||
} |
|
||||
|
|
||||
|
|
||||
/*
|
|
||||
* vectorNorm calculates the first norm of a vector. |
|
||||
*/ |
|
||||
double vectorNorm(double *vector, int vectorSize) { |
|
||||
double norm = 0.; |
|
||||
|
|
||||
for (int i=0; i<vectorSize; ++i) { |
|
||||
norm += fabs(vector[i]); |
|
||||
} |
|
||||
|
|
||||
return norm; |
|
||||
} |
|
||||
|
|
||||
// ==================== PROGRAM INPUT AND OUTPUT UTILS ====================
|
|
||||
|
|
||||
/*
|
|
||||
* parseArguments parses the command line arguments given by the user. |
|
||||
*/ |
|
||||
void parseArguments(int argumentCount, char **argumentVector, Parameters *parameters) { |
|
||||
if (argumentCount < 2 || argumentCount > 10) { |
|
||||
validUsage(argumentVector[0]); |
|
||||
} |
|
||||
|
|
||||
(*parameters).numberOfPages = 0; |
|
||||
(*parameters).maxIterations = 0; |
|
||||
(*parameters).convergenceCriterion = 1; |
|
||||
(*parameters).dampingFactor = 0.85; |
|
||||
(*parameters).verbose = false; |
|
||||
(*parameters).history = false; |
|
||||
(*parameters).outputFilename = DEFAULT_OUTPUT_FILENAME; |
|
||||
|
|
||||
char *endPointer; |
|
||||
int argumentIndex = 1; |
|
||||
|
|
||||
while (argumentIndex < argumentCount) { |
|
||||
if (!strcmp(argumentVector[argumentIndex], ARGUMENT_CONVERGENCE_TOLERANCE)) { |
|
||||
argumentIndex = checkIncrement(argumentIndex, argumentCount, argumentVector[0]); |
|
||||
|
|
||||
double convergenceInput = strtod(argumentVector[argumentIndex], &endPointer); |
|
||||
if (convergenceInput == 0) { |
|
||||
printf("Invalid convergence argument\n"); |
|
||||
exit(EXIT_FAILURE); |
|
||||
} |
|
||||
(*parameters).convergenceCriterion = convergenceInput; |
|
||||
} else if (!strcmp(argumentVector[argumentIndex], ARGUMENT_MAX_ITERATIONS)) { |
|
||||
argumentIndex = checkIncrement(argumentIndex, argumentCount, argumentVector[0]); |
|
||||
|
|
||||
size_t iterationsInput = strtol(argumentVector[argumentIndex], &endPointer, NUMERICAL_BASE); |
|
||||
if (iterationsInput == 0 && endPointer) { |
|
||||
printf("Invalid iterations argument\n"); |
|
||||
exit(EXIT_FAILURE); |
|
||||
} |
|
||||
(*parameters).maxIterations = iterationsInput; |
|
||||
} else if (!strcmp(argumentVector[argumentIndex], ARGUMENT_DAMPING_FACTOR)) { |
|
||||
argumentIndex = checkIncrement(argumentIndex, argumentCount, argumentVector[0]); |
|
||||
|
|
||||
double alphaInput = strtod(argumentVector[argumentIndex], &endPointer); |
|
||||
if ((alphaInput == 0 || alphaInput > 1) && endPointer) { |
|
||||
printf("Invalid alpha argument\n"); |
|
||||
exit(EXIT_FAILURE); |
|
||||
} |
|
||||
(*parameters).dampingFactor = alphaInput; |
|
||||
} else if (!strcmp(argumentVector[argumentIndex], ARGUMENT_VERBAL_OUTPUT)) { |
|
||||
(*parameters).verbose = true; |
|
||||
} else if (!strcmp(argumentVector[argumentIndex], ARGUMENT_OUTPUT_HISTORY)) { |
|
||||
(*parameters).history = true; |
|
||||
} else if (!strcmp(argumentVector[argumentIndex], ARGUMENT_OUTPUT_FILENAME)) { |
|
||||
argumentIndex = checkIncrement(argumentIndex, argumentCount, argumentVector[0]); |
|
||||
|
|
||||
if (fopen(argumentVector[argumentIndex], "w") == NULL) { |
|
||||
printf("Invalid output filename. Reverting to default.\n"); |
|
||||
continue; |
|
||||
} |
|
||||
(*parameters).outputFilename = argumentVector[argumentIndex]; |
|
||||
} else if (argumentIndex == argumentCount - 1) { |
|
||||
(*parameters).graphFilename = argumentVector[argumentIndex]; |
|
||||
} else { |
|
||||
validUsage(argumentVector[0]); |
|
||||
exit(EXIT_FAILURE); |
|
||||
} |
|
||||
++argumentIndex; |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
/*
|
|
||||
* readGraphFromFile loads the file supplied in the command line arguments to an |
|
||||
* array (directedWebGraph) that represents the graph. |
|
||||
*/ |
|
||||
void generateNormalizedTransitionMatrixFromFile(CsrSparseMatrix *transitionMatrix, |
|
||||
Parameters *parameters){ |
|
||||
FILE *graphFile; |
|
||||
|
|
||||
// Opens the file for reading
|
|
||||
graphFile = fopen((*parameters).graphFilename, "r+"); |
|
||||
if (!graphFile) { |
|
||||
printf("Error opening file \n"); |
|
||||
exit(EXIT_FAILURE); |
|
||||
} |
|
||||
|
|
||||
char buffer[FILE_READ_BUFFER_SIZE]; |
|
||||
char *readResult; |
|
||||
// Skips the first two lines
|
|
||||
readResult = fgets(buffer, FILE_READ_BUFFER_SIZE, graphFile); |
|
||||
readResult = fgets(buffer, FILE_READ_BUFFER_SIZE, graphFile); |
|
||||
if (readResult == NULL) { |
|
||||
printf("Error while reading from the file. Does the file have the correct format?\n"); |
|
||||
exit(EXIT_FAILURE); |
|
||||
} |
|
||||
|
|
||||
// Third line contains the numbers of nodes and edges
|
|
||||
int numberOfNodes = 0, numberOfEdges = 0; |
|
||||
|
|
||||
readResult = fgets(buffer, FILE_READ_BUFFER_SIZE, graphFile); |
|
||||
if (readResult == NULL) { |
|
||||
printf("Error while reading from the file. Does the file have the correct format?\n"); |
|
||||
exit(EXIT_FAILURE); |
|
||||
} |
|
||||
|
|
||||
// Parses the number of nodes and number of edges
|
|
||||
{ |
|
||||
// Splits string to whitespace
|
|
||||
char *token = strtok(buffer, " "); |
|
||||
bool nextIsNodes = false, nextIsEdges = false; |
|
||||
|
|
||||
while (token != NULL) { |
|
||||
if (strcmp(token, "Nodes:") == 0) { |
|
||||
nextIsNodes = true; |
|
||||
} else if (nextIsNodes) { |
|
||||
numberOfNodes = atoi(token); |
|
||||
nextIsNodes = false; |
|
||||
} else if (strcmp(token, "Edges:") == 0) { |
|
||||
nextIsEdges = true; |
|
||||
} else if (nextIsEdges) { |
|
||||
numberOfEdges = atoi(token); |
|
||||
break; |
|
||||
} |
|
||||
|
|
||||
// Gets next string token
|
|
||||
token = strtok (NULL, " ,.-"); |
|
||||
} |
|
||||
} |
|
||||
|
|
||||
if ((*parameters).verbose) { |
|
||||
printf("File claims number of pages is: %d\nThe number of edges is: %d\n", |
|
||||
numberOfNodes, numberOfEdges); |
|
||||
} |
|
||||
|
|
||||
// Skips the fourth line
|
|
||||
readResult = fgets(buffer, 512, graphFile); |
|
||||
if (readResult == NULL) { |
|
||||
printf("Error while reading from the file. Does the file have the correct format?\n"); |
|
||||
exit(EXIT_FAILURE); |
|
||||
} |
|
||||
|
|
||||
|
|
||||
int maxPageIndex = 0; |
|
||||
CooSparseMatrix tempMatrix = initCooSparseMatrix(); |
|
||||
allocMemoryForCoo(&tempMatrix, numberOfEdges); |
|
||||
|
|
||||
for (int i=0; i<numberOfEdges; i++) { |
|
||||
int fileFrom = 0, fileTo = 0; |
|
||||
if (!fscanf(graphFile, "%d %d", &fileFrom, &fileTo)) { |
|
||||
break; |
|
||||
} |
|
||||
|
|
||||
if (fileFrom > maxPageIndex) { |
|
||||
maxPageIndex = fileFrom; |
|
||||
} |
|
||||
if (fileTo > maxPageIndex) { |
|
||||
maxPageIndex = fileTo; |
|
||||
} |
|
||||
addElement(&tempMatrix, 1, fileFrom, fileTo); |
|
||||
} |
|
||||
|
|
||||
if ((*parameters).verbose) { |
|
||||
printf("Max page index found is: %d\n", maxPageIndex); |
|
||||
} |
|
||||
(*parameters).numberOfPages = maxPageIndex + 1; |
|
||||
|
|
||||
// Calculates the outdegree of each page and assigns the uniform probability
|
|
||||
// of transition to the elements of the corresponding row
|
|
||||
int* pageOutdegree = malloc((*parameters).numberOfPages*sizeof(int)); |
|
||||
for (int i=0; i<(*parameters).numberOfPages; ++i){ |
|
||||
pageOutdegree[i] = 0; |
|
||||
} |
|
||||
|
|
||||
for (int i=0; i<numberOfEdges; ++i) { |
|
||||
int currentRow = tempMatrix.elements[i]->rowIndex; |
|
||||
++pageOutdegree[currentRow]; |
|
||||
} |
|
||||
|
|
||||
for (int i=0; i<tempMatrix.size; ++i) { |
|
||||
tempMatrix.elements[i]->value = 1./pageOutdegree[tempMatrix.elements[i]->rowIndex]; |
|
||||
} |
|
||||
free(pageOutdegree); |
|
||||
|
|
||||
// Transposes the temporary transition matrix (P^T).
|
|
||||
transposeSparseMatrix(&tempMatrix); |
|
||||
|
|
||||
allocMemoryForCsr(transitionMatrix, numberOfEdges); |
|
||||
// Transforms the temporary COO matrix to the desired CSR format
|
|
||||
transformToCSR(tempMatrix, transitionMatrix); |
|
||||
destroyCooSparseMatrix(&tempMatrix); |
|
||||
|
|
||||
fclose(graphFile); |
|
||||
} |
|
||||
|
|
||||
/*
|
|
||||
* validUsage outputs a message to the console that informs the user of the |
|
||||
* correct (valid) way to use the program. |
|
||||
*/ |
|
||||
void validUsage(char *programName) { |
|
||||
printf("%s [-c convergence_criterion] [-m max_iterations] [-a alpha] [-v] [-h] [-o output_filename] <graph_file>" \ |
|
||||
"\n-c convergence_criterion" \ |
|
||||
"\n\tthe convergence tolerance criterion" \ |
|
||||
"\n-m max_iterations" \ |
|
||||
"\n\tmaximum number of iterations to perform" \ |
|
||||
"\n-a alpha" \ |
|
||||
"\n\tthe damping factor" \ |
|
||||
"\n-v enable verbal output" \ |
|
||||
"\n-h enable history output to file" \ |
|
||||
"\n-o output_filename" \ |
|
||||
"\n\tfilename and path for the output" \ |
|
||||
"\n", programName); |
|
||||
exit(EXIT_FAILURE); |
|
||||
} |
|
||||
|
|
||||
/*
|
|
||||
* checkIncrement is a helper function for parseArguments function. |
|
||||
*/ |
|
||||
int checkIncrement(int previousIndex, int maxIndex, char *programName) { |
|
||||
if (previousIndex == maxIndex) { |
|
||||
validUsage(programName); |
|
||||
exit(EXIT_FAILURE); |
|
||||
} |
|
||||
return ++previousIndex; |
|
||||
} |
|
||||
|
|
||||
void savePagerankToFile(char *filename, bool append, double *pagerankVector, |
|
||||
int vectorSize, int iteration) { |
|
||||
FILE *outputFile; |
|
||||
|
|
||||
if (append) { |
|
||||
outputFile = fopen(filename, "a"); |
|
||||
} else { |
|
||||
outputFile = fopen(filename, "w"); |
|
||||
} |
|
||||
|
|
||||
if (outputFile == NULL) { |
|
||||
printf("Error while opening the output file.\n"); |
|
||||
return; |
|
||||
} |
|
||||
|
|
||||
// Saves the pagerank vector
|
|
||||
//fprintf(outputFile, "Iteration %d:\t", iteration);
|
|
||||
double sum = 0; |
|
||||
for (int i=0; i<vectorSize; ++i) { |
|
||||
sum += pagerankVector[i]; |
|
||||
} |
|
||||
//fprintf(outputFile, "%f\n", sum);
|
|
||||
|
|
||||
for (int i=0; i<vectorSize; ++i) { |
|
||||
fprintf(outputFile, "%d = %.10g\n", i, pagerankVector[i]/sum); |
|
||||
} |
|
||||
|
|
||||
fclose(outputFile); |
|
||||
} |
|
||||
@ -1,100 +0,0 @@ |
|||||
#ifndef SERIAL_GS_PAGERANK_FUNCTIONS_H /* Include guard */ |
|
||||
#define SERIAL_GS_PAGERANK_FUNCTIONS_H |
|
||||
|
|
||||
/* ===== INCLUDES ===== */ |
|
||||
|
|
||||
#include <stdbool.h> |
|
||||
#include <stdio.h> |
|
||||
#include <stdlib.h> |
|
||||
#include <string.h> |
|
||||
#include <math.h> |
|
||||
|
|
||||
#include "coo_sparse_matrix.h" |
|
||||
|
|
||||
/* ===== DEFINITIONS ===== */ |
|
||||
|
|
||||
//Colors used for better console output formating.
|
|
||||
#define ANSI_COLOR_RED "\x1B[31m" |
|
||||
#define ANSI_COLOR_GREEN "\x1B[32m" |
|
||||
#define ANSI_COLOR_YELLOW "\x1B[33m" |
|
||||
#define ANSI_COLOR_BLUE "\x1B[34m" |
|
||||
#define ANSI_COLOR_CYAN "\x1B[36m" |
|
||||
#define ANSI_COLOR_RESET "\x1B[0m" |
|
||||
|
|
||||
/* ===== CONSTANTS DEFINITION ===== */ |
|
||||
|
|
||||
// Constant strings that store the command line options available.
|
|
||||
extern const char *ARGUMENT_CONVERGENCE_TOLERANCE; |
|
||||
extern const char *ARGUMENT_MAX_ITERATIONS; |
|
||||
extern const char *ARGUMENT_DAMPING_FACTOR; |
|
||||
extern const char *ARGUMENT_VERBAL_OUTPUT; |
|
||||
extern const char *ARGUMENT_OUTPUT_HISTORY; |
|
||||
extern const char *ARGUMENT_OUTPUT_FILENAME; |
|
||||
// The numerical base used when parsing numerical command line arguments.
|
|
||||
extern const int NUMERICAL_BASE; |
|
||||
// Default filename used for the output.
|
|
||||
extern char *DEFAULT_OUTPUT_FILENAME; |
|
||||
// The size of the buffer used for reading the graph input file.
|
|
||||
extern const int FILE_READ_BUFFER_SIZE; |
|
||||
|
|
||||
/* ===== STRUCTURES ===== */ |
|
||||
|
|
||||
// A data structure to conveniently hold the algorithm's parameters.
|
|
||||
typedef struct parameters { |
|
||||
int numberOfPages, maxIterations; |
|
||||
double convergenceCriterion, dampingFactor; |
|
||||
bool verbose, history; |
|
||||
char *outputFilename, *graphFilename; |
|
||||
} Parameters; |
|
||||
|
|
||||
/* ===== FUNCTION DEFINITIONS ===== */ |
|
||||
|
|
||||
// Function validUsage outputs the correct way to use the program with command
|
|
||||
// line arguments.
|
|
||||
void validUsage(char *programName); |
|
||||
|
|
||||
// Function checkIncrement is a helper function used in parseArguments (see
|
|
||||
// bellow).
|
|
||||
int checkIncrement(int previousIndex, int maxIndex, char *programName); |
|
||||
|
|
||||
// Function parseArguments parses command line arguments.
|
|
||||
void parseArguments(int argumentCount, char **argumentVector, |
|
||||
Parameters *parameters); |
|
||||
|
|
||||
// Function generateNormalizedTransitionMatrixFromFile reads through the entries
|
|
||||
// of the file specified in the arguments (parameters->graphFilename), using
|
|
||||
// them to populate the sparse array (transitionMatrix). The entries of the file
|
|
||||
// represent the edges of the web transition graph. The entries are then
|
|
||||
// modified to become the rows of the transition matrix.
|
|
||||
void generateNormalizedTransitionMatrixFromFile(CsrSparseMatrix *transitionMatrix, |
|
||||
Parameters *parameters); |
|
||||
|
|
||||
// Function savePagerankToFile appends or overwrites the pagerank vector
|
|
||||
// "pagerankVector" to the file with the filename supplied in the arguments.
|
|
||||
void savePagerankToFile(char *filename, bool append, double *pagerankVector, |
|
||||
int vectorSize, int iteration); |
|
||||
|
|
||||
// Function initialize allocates memory for the pagerank vector, reads the
|
|
||||
// dataset from the file and creates the transition probability distribution
|
|
||||
// matrix.
|
|
||||
void initialize(CsrSparseMatrix *transitionMatrix, double **pagerankVector, |
|
||||
Parameters *parameters); |
|
||||
|
|
||||
// Function vectorNorm calculates the first norm of a vector.
|
|
||||
double vectorNorm(double *vector, int vectorSize); |
|
||||
|
|
||||
// Function calculateNextPagerank calculates the next pagerank vector.
|
|
||||
void calculateNextPagerank(CsrSparseMatrix *transitionMatrix, |
|
||||
double *previousPagerankVector, double **pagerankVector, |
|
||||
double *linksFromConvergedPagesPagerankVector, |
|
||||
double *convergedPagerankVector, int vectorSize, double dampingFactor, pthread_t *threads, int threadNum); |
|
||||
|
|
||||
// Function pagerank iteratively calculates the pagerank of each page until
|
|
||||
// either the convergence criterion is met or the maximum number of iterations
|
|
||||
// is reached.
|
|
||||
int pagerank(CsrSparseMatrix *transitionMatrix, double **pagerankVector, |
|
||||
bool *convergenceStatus, Parameters parameters); |
|
||||
|
|
||||
void csrSparseMatrixVectorMultiplication_threads(void* arg); |
|
||||
void compPagerankVector_threads(void* arg); |
|
||||
#endif // SERIAL_GS_PAGERANK_FUNCTIONS_H
|
|
||||
Loading…
Reference in new issue