1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156
|
//------------------------------------------------------------------------------
// irowscale: scale the rows of an adjacency matrix by out-degree
//------------------------------------------------------------------------------
// SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2020, All Rights Reserved.
// http://suitesparse.com See GraphBLAS/Doc/License.txt for license.
//------------------------------------------------------------------------------
// on input, A is a square unsymmetric binary matrix of size n-by-n, of any
// built-in type. On output, C is a rowscaled version of A, of type
// GrB_UINT64, with C = D*A + I. The diagonal matrix D has diagonal entries
// D(i,i)=(2^30)//sum(A(i,:)), or D(i,i) is not present if A(i,:) is empty.
// The matrix I has entries only on the diagonal, all equal to zero. This
// optional step ensures that C has no empty columns, which speeds up the
// subsequent PageRank computation.
/* MATLAB equivalent (excluding the addition of I):
function C = rowscale (A)
%ROWSCALE row scale an adjacency matrix by out-degree
% C = rowscale (A)
% scale the adjacency matrix by out-degree
dout = sum (A,2) ; % dout(i) is the out-degree of node i
is_nonempty = (dout > 0) ; % find vertices with outgoing edges
nonempty = find (is_nonempty) ; % list of vertices with outgoing edges
n = size (A,1) ;
% divide each non-empty row by its out-degree
dinv = 1 ./ dout (is_nonempty) ;
D = sparse (nonempty, nonempty, dinv, n, n) ;
C = D*A ;
C = floor ((2^30) * C) ; % scale the result to integer
*/
#include "GraphBLAS.h"
//------------------------------------------------------------------------------
// helper macros
//------------------------------------------------------------------------------
// free all workspace
#define FREEWORK \
{ \
GrB_Vector_free (&dout) ; \
GrB_Matrix_free (&D) ; \
if (I != NULL) free (I) ; \
if (X != NULL) free (X) ; \
}
// error handler: free workspace and the output matrix C
#define FREE_ALL \
{ \
GrB_Matrix_free (&C) ; \
FREEWORK ; \
}
#undef GB_PUBLIC
#define GB_LIBRARY
#include "graphblas_demos.h"
//------------------------------------------------------------------------------
// irowscale: C = D*A + I*0 where D(i,i) = ZSCALE/sum(A(i,:)
//------------------------------------------------------------------------------
GB_PUBLIC
GrB_Info irowscale // GrB_SUCCESS or error condition
(
GrB_Matrix *Chandle, // output matrix C = rowscale (A)
GrB_Matrix A // input matrix, not modified
)
{
//--------------------------------------------------------------------------
// intializations
//--------------------------------------------------------------------------
GrB_Info info ;
GrB_Vector dout = NULL ;
GrB_Matrix D = NULL, C = NULL ;
GrB_Index *I = NULL ;
uint64_t *X = NULL ;
// n = size (A,1) ;
GrB_Index n ;
OK (GrB_Matrix_nrows (&n, A)) ;
//--------------------------------------------------------------------------
// dout = sum (A,2) ; // dout(i) is the out-degree of node i
//--------------------------------------------------------------------------
OK (GrB_Vector_new (&dout, GrB_UINT64, n)) ;
OK (GrB_Matrix_reduce_BinaryOp (dout, NULL, NULL, GrB_PLUS_UINT64,
A, NULL)) ;
//--------------------------------------------------------------------------
// construct scaling matrix D
//--------------------------------------------------------------------------
// D is diagonal with D(i,i) = ZSCALE/sum(A(i,:)), or D(i,i) is not present
// if row i of A has no entries.
// [I,~,X] = find (dout) ;
I = (GrB_Index *) malloc ((n+1) * sizeof (GrB_Index)) ;
X = (uint64_t *) malloc ((n+1) * sizeof (uint64_t)) ;
CHECK (I != NULL && X != NULL, GrB_OUT_OF_MEMORY) ;
GrB_Index nvals = n ;
OK (GrB_Vector_extractTuples_UINT64 (I, X, &nvals, dout)) ;
// I and X exclude empty columns of A. This condition is always true.
CHECK (nvals <= n, GrB_PANIC) ;
// D = diag (ZSCALE./dout) ;
OK (GrB_Matrix_new (&D, GrB_UINT64, n, n)) ;
for (int64_t i = 0 ; i < nvals ; i++)
{
// X (i) = ZSCALE / X (i), but make sure it doesn't underflow to zero.
// Underflow would only occur if a node has degree higher than 2^30.
uint64_t y = ZSCALE / X [i] ;
if (y == 0) y = 1 ;
X [i] = y ;
}
OK (GrB_Matrix_build_UINT64 (D, I, I, X, nvals, GrB_PLUS_UINT64)) ;
//--------------------------------------------------------------------------
// C = diagonal matrix with explicit zeros on diagonal
//--------------------------------------------------------------------------
// This step is optional, but it ensures that C has no non-empty columns,
// which speeds up the pagerank computation by ensuring r*C remains a dense
// vector.
for (int64_t i = 0 ; i < n ; i++)
{
I [i] = i ;
X [i] = 0 ;
}
OK (GrB_Matrix_new (&C, GrB_UINT64, n, n)) ;
OK (GrB_Matrix_build_UINT64 (C, I, I, X, n, GrB_PLUS_UINT64)) ;
//--------------------------------------------------------------------------
// C += D*A
//--------------------------------------------------------------------------
OK (GrB_mxm (C, NULL, GrB_PLUS_UINT64, GxB_PLUS_TIMES_UINT64, D, A, NULL)) ;
//--------------------------------------------------------------------------
// free workspace and return result
//--------------------------------------------------------------------------
FREEWORK ;
(*Chandle) = C ;
return (GrB_SUCCESS) ;
}
|