1function test51
2%TEST51 test GxB_subassign, multiply operations
3
4% SuiteSparse:GraphBLAS, Timothy A. Davis, (c) 2017-2021, All Rights Reserved.
5% SPDX-License-Identifier: Apache-2.0
6
7fprintf ('\n-----------performance test GB_mex_subassign, multiple ops\n') ;
8
9[save save_chunk] = nthreads_get ;
10chunk = 4096 ;
11nthreads = feature ('numcores') ;
12nthreads_set (nthreads, chunk) ;
13
14rng ('default')
15
16for problem = 1:5
17
18    fprintf ('\n----------------------------------------------------- %d\n',...
19    problem) ;
20
21    switch (problem)
22    case 1
23        Corig = sprand (1000, 1, 0.1) ;
24        nwork = 10 ;
25        m_max = 1000 ;
26        n_max = 1 ;
27        d = 0.1 ;
28    case 2
29        Prob = ssget ('HB/west0067') ; Corig = Prob.A ;
30        nwork = 1000 ;
31        m_max = 64 ;
32        n_max = 64 ;
33        d = 0.1 ;
34    case 3
35        Corig = abs (sprand (80, 80, 0.4)) ;
36        nwork = 1000 ;
37        m_max = 64 ;
38        n_max = 64 ;
39        d = 0.1 ;
40    case 4
41        Corig = sparse (rand (50000, 500)) ;
42        nwork = 50 ;
43        % [m n] = size (A) ;
44        [m n] = size (Corig) ;
45        m_max = m ;
46        n_max = n ;
47        d = 100 / (m*n) ;
48    case 5
49        Prob = ssget (2662)
50        Corig = Prob.A ;
51        nwork = 100 ;
52        m_max = 1024 ;
53        n_max = 1024 ;
54        d = 10000 / (ni*nj) ;
55    end
56
57    Corig = abs (Corig) ;
58    [m n] = size (Corig) ;
59
60    for opkind = 1 % 1:3
61
62        fprintf ('\n--------------------------\n') ;
63        switch opkind
64        case 1
65            fprintf ('C(I,J) = accum(C(I,J),A) one op, assemble at end\n') ;
66        case 2
67            fprintf ('C(I,J) = accum(C(I,J),A) two ops, change every 10th\n') ;
68        case 3
69            fprintf ('C(I,J) = A no accum\n') ;
70        end
71
72        nz = 0 ;
73        ni_min = inf ;
74        ni_max = -inf ;
75        nj_min = inf ;
76        nj_max = -inf ;
77        nz_min = inf ;
78        nz_max = -inf ;
79
80        clear Work
81
82        for k = 1:nwork
83            ni = double (irand (1, min (m,m_max))) ;
84            nj = double (irand (1, min (n,n_max))) ;
85            ni_min = min (ni_min, ni) ;
86            nj_min = min (nj_min, nj) ;
87            ni_max = max (ni_max, ni) ;
88            nj_max = max (nj_max, nj) ;
89
90            if (d == 1)
91                A = sparse (rand (ni, nj)) ;
92            else
93                A = sprand (ni, nj, d) ;
94            end
95            nz_min = min (nz_min, nnz (A)) ;
96            nz_max = max (nz_max, nnz (A)) ;
97            if (opkind == 1)
98                op = 'plus' ;
99            else
100                if (mod (k,5) == 0)
101                    op = 'max' ;
102                else
103                    op = 'plus' ;
104                end
105            end
106            nz = nz + nnz (A) ;
107            I = randperm (m,ni) ;
108            J = randperm (n,nj) ;
109            Work (k).A = A ;
110            Work (k).I = I ;
111            Work (k).J = J ;
112            Work (k).accum = op ;
113        end
114        fprintf ('number of C(I,J) = ... to do: %d\n', nwork) ;
115        fprintf ('ni: [ %d to %d]\n', ni_min, ni_max) ;
116        fprintf ('nj: [ %d to %d]\n', nj_min, nj_max) ;
117        fprintf ('nz: [ %d to %d]\n', nz_min, nz_max) ;
118        fprintf ('C is %d-by-%d nnz(A): %d  nz to add: %d  matrices: %d\n', ...
119            m, n, nnz (Corig), nz, nwork) ;
120
121        Work2 = Work ;
122        for k = 1:nwork
123            Work2 (k).I = uint64 (Work2 (k).I - 1) ;
124            Work2 (k).J = uint64 (Work2 (k).J - 1) ;
125        end
126        tic
127        C2 = GB_mex_subassign (Corig, Work2) ;
128        t1 = toc ;
129        fprintf ('GraphBLAS time: %g\n', t1) ;
130        fprintf ('final nnz: %d\n', nnz (C2.matrix)) ;
131
132        fprintf ('start MATLAB...\n') ;
133        tic
134        C = Corig ;
135        % full (C)
136        for k = 1:nwork
137            I = Work (k).I ;
138            J = Work (k).J ;
139            A = Work (k).A ;
140            % Afull = full (A)
141            if (isequal (Work (k).accum, 'plus'))
142                % fprintf (' %d : plus\n', k) ;
143                C (I,J) = C (I,J) + A ;
144            elseif (isequal (Work (k).accum, 'max'))
145                % fprintf (' %d : max\n', k) ;
146                C (I,J) = max (C (I,J),A) ;
147            else
148                C (I,J) = A ;
149            end
150            % full (C)
151        end
152        t2 = toc ;
153        fprintf ('MATLAB    time: %g\n', t2) ;
154        fprintf ('GraphBLAS speedup: %g\n', t2/t1) ;
155
156        % C2.matrix
157        % C - C2.matrix
158
159        % C2 = full (C2.matrix)
160        % C-C2
161        % assert (isequal (C, C2)) ;
162        assert (isequal (C, C2.matrix)) ;
163    end
164
165end
166
167fprintf ('\ntest51: all tests passed\n') ;
168
169nthreads_set (save, save_chunk) ;
170