File: irowscale.c

package info (click to toggle)
suitesparse 1%3A5.8.1%2Bdfsg-2
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 152,716 kB
  • sloc: ansic: 774,385; cpp: 24,213; makefile: 6,310; fortran: 1,927; java: 1,826; csh: 1,686; ruby: 725; sh: 535; perl: 225; python: 209; sed: 164; awk: 60
file content (156 lines) | stat: -rw-r--r-- 5,830 bytes parent folder | download
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) ;
}