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