1 /****************************************************************/
2 /* Parallel Combinatorial BLAS Library (for Graph Computations) */
3 /* version 1.6 -------------------------------------------------*/
4 /* date: 6/15/2017 ---------------------------------------------*/
5 /* authors: Ariful Azad, Aydin Buluc --------------------------*/
6 /****************************************************************/
7 /*
8 Copyright (c) 2010-2017, The Regents of the University of California
9
10 Permission is hereby granted, free of charge, to any person obtaining a copy
11 of this software and associated documentation files (the "Software"), to deal
12 in the Software without restriction, including without limitation the rights
13 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
14 copies of the Software, and to permit persons to whom the Software is
15 furnished to do so, subject to the following conditions:
16
17 The above copyright notice and this permission notice shall be included in
18 all copies or substantial portions of the Software.
19
20 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
21 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
23 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
25 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
26 THE SOFTWARE.
27 */
28
29
30 #include "SpCCols.h"
31 #include "Deleter.h"
32 #include <algorithm>
33 #include <functional>
34 #include <vector>
35 #include <climits>
36 #include <iomanip>
37 #include <cassert>
38
39 namespace combblas {
40
41 /****************************************************************************/
42 /********************* PUBLIC CONSTRUCTORS/DESTRUCTORS **********************/
43 /****************************************************************************/
44
45 template <class IT, class NT>
46 const IT SpCCols<IT,NT>::esscount = static_cast<IT>(3);
47
48
49 template <class IT, class NT>
SpCCols()50 SpCCols<IT,NT>::SpCCols():csc(NULL), m(0), n(0), nnz(0), splits(0){
51 }
52
53 // Allocate all the space necessary
54 template <class IT, class NT>
SpCCols(IT size,IT nRow,IT nCol)55 SpCCols<IT,NT>::SpCCols(IT size, IT nRow, IT nCol)
56 :m(nRow), n(nCol), nnz(size), splits(0)
57 {
58 if(nnz > 0)
59 csc = new Csc<IT,NT>(nnz, n);
60 else
61 csc = NULL;
62 }
63
64 template <class IT, class NT>
~SpCCols()65 SpCCols<IT,NT>::~SpCCols()
66 {
67 if(nnz > 0)
68 {
69 if(csc != NULL)
70 {
71 if(splits > 0)
72 {
73 for(int i=0; i<splits; ++i)
74 delete cscarr[i];
75 delete [] cscarr;
76 }
77 else
78 {
79 delete csc;
80 }
81 }
82 }
83 }
84
85 // Copy constructor (constructs a new object. i.e. this is NEVER called on an existing object)
86 // Derived's copy constructor can safely call Base's default constructor as base has no data members
87 template <class IT, class NT>
SpCCols(const SpCCols<IT,NT> & rhs)88 SpCCols<IT,NT>::SpCCols(const SpCCols<IT,NT> & rhs)
89 : m(rhs.m), n(rhs.n), nnz(rhs.nnz), splits(rhs.splits)
90 {
91 if(splits > 0)
92 {
93 for(int i=0; i<splits; ++i)
94 CopyCsc(rhs.cscarr[i]);
95 }
96 else
97 {
98 CopyCsc(rhs.csc);
99 }
100 }
101
102 /**
103 * Constructor for converting SpTuples matrix -> SpCCols
104 * @param[in] rhs if transpose=true,
105 * \n then rhs is assumed to be a row sorted SpTuples object
106 * \n else rhs is assumed to be a column sorted SpTuples object
107 **/
108 template <class IT, class NT>
SpCCols(const SpTuples<IT,NT> & rhs,bool transpose)109 SpCCols<IT,NT>::SpCCols(const SpTuples<IT, NT> & rhs, bool transpose)
110 : m(rhs.m), n(rhs.n), nnz(rhs.nnz), splits(0)
111 {
112
113 if(nnz == 0) // m by n matrix of complete zeros
114 {
115 if(transpose) std::swap(m,n);
116 csc = NULL;
117 }
118 else
119 {
120 if(transpose)
121 {
122 std::swap(m,n);
123 csc = new Csc<IT,NT>(nnz,n); // the swap is already done here
124 std::vector< std::pair<IT,NT> > tosort (nnz);
125 std::vector<IT> work(n+1, (IT) 0 ); // workspace, zero initialized, first entry stays zero
126 for (IT k = 0 ; k < nnz ; ++k)
127 {
128 IT tmp = rhs.rowindex(k);
129 work [ tmp+1 ]++ ; // column counts (i.e, w holds the "col difference array")
130 }
131 if(nnz > 0)
132 {
133 std::partial_sum(work.begin(), work.end(), work.begin());
134 std::copy(work.begin(), work.end(), csc->jc);
135 IT last;
136 for (IT k = 0 ; k < nnz ; ++k)
137 {
138 tosort[ work[ rhs.rowindex(k) ]++] = std::make_pair( rhs.colindex(k), rhs.numvalue(k));
139 }
140 #ifdef _OPENMP
141 #pragma omp parallel for
142 #endif
143 for(int i=0; i< n; ++i)
144 {
145 sort(tosort.begin() + csc->jc[i], tosort.begin() + csc->jc[i+1]);
146
147 IT ind;
148 typename std::vector<std::pair<IT,NT> >::iterator itr; // iterator is a dependent name
149 for(itr = tosort.begin() + csc->jc[i], ind = csc->jc[i]; itr != tosort.begin() + csc->jc[i+1]; ++itr, ++ind)
150 {
151 csc->ir[ind] = itr->first;
152 csc->num[ind] = itr->second;
153 }
154 }
155 }
156 }
157 else
158 {
159 csc = new Csc<IT,NT>(nnz,n); // the swap is already done here
160 std::vector< std::pair<IT,NT> > tosort (nnz);
161 std::vector<IT> work(n+1, (IT) 0 ); // workspace, zero initialized, first entry stays zero
162 for (IT k = 0 ; k < nnz ; ++k)
163 {
164 IT tmp = rhs.colindex(k);
165 work [ tmp+1 ]++ ; // column counts (i.e, w holds the "col difference array")
166 }
167 if(nnz > 0)
168 {
169 std::partial_sum(work.begin(), work.end(), work.begin());
170 std::copy(work.begin(), work.end(), csc->jc);
171 IT last;
172 for (IT k = 0 ; k < nnz ; ++k)
173 {
174 tosort[ work[ rhs.colindex(k) ]++] = std::make_pair( rhs.rowindex(k), rhs.numvalue(k));
175 }
176 #ifdef _OPENMP
177 #pragma omp parallel for
178 #endif
179 for(int i=0; i< n; ++i)
180 {
181 sort(tosort.begin() + csc->jc[i], tosort.begin() + csc->jc[i+1]);
182
183 IT ind;
184 typename std::vector<std::pair<IT,NT> >::iterator itr; // iterator is a dependent name
185 for(itr = tosort.begin() + csc->jc[i], ind = csc->jc[i]; itr != tosort.begin() + csc->jc[i+1]; ++itr, ++ind)
186 {
187 csc->ir[ind] = itr->first;
188 csc->num[ind] = itr->second;
189 }
190 }
191 }
192 }
193 }
194 }
195
196 /****************************************************************************/
197 /************************** PUBLIC OPERATORS ********************************/
198 /****************************************************************************/
199
200 /**
201 * The assignment operator operates on an existing object
202 * The assignment operator is the only operator that is not inherited.
203 * But there is no need to call base's assigment operator as it has no data members
204 */
205 template <class IT, class NT>
operator =(const SpCCols<IT,NT> & rhs)206 SpCCols<IT,NT> & SpCCols<IT,NT>::operator=(const SpCCols<IT,NT> & rhs)
207 {
208 // this pointer stores the address of the class instance
209 // check for self assignment using address comparison
210 if(this != &rhs)
211 {
212 if(csc != NULL && nnz > 0)
213 {
214 delete csc;
215 }
216 if(rhs.csc != NULL)
217 {
218 csc = new Csc<IT,NT>(*(rhs.csc));
219 nnz = rhs.nnz;
220 }
221 else
222 {
223 csc = NULL;
224 nnz = 0;
225 }
226
227 m = rhs.m;
228 n = rhs.n;
229 splits = rhs.splits;
230 }
231 return *this;
232 }
233
234
235 template <class IT, class NT>
RowSplit(int numsplits)236 void SpCCols<IT,NT>::RowSplit(int numsplits)
237 {
238 splits = numsplits;
239 IT perpiece = m / splits;
240 std::vector<IT> nnzs(splits, 0);
241 std::vector < std::vector < std::tuple<IT,IT,NT> > > colrowpairs(splits);
242 std::vector< std::vector<IT> > colcnts(splits);
243 for(int i=0; i< splits; ++i)
244 colcnts[i].resize(n, 0);
245
246 if(nnz > 0 && csc != NULL)
247 {
248 for(IT i=0; i< csc->n; ++i)
249 {
250 for(IT j = csc->jc[i]; j< csc->jc[i+1]; ++j)
251 {
252 IT rowid = csc->ir[j]; // colid=i
253 IT owner = std::min(rowid / perpiece, static_cast<IT>(splits-1));
254 colrowpairs[owner].push_back(std::make_tuple(i, rowid - owner*perpiece, csc->num[i]));
255
256 ++(colcnts[owner][i]);
257 ++(nnzs[owner]);
258 }
259 }
260 }
261 delete csc; // claim memory
262 cscarr = new Csc<IT,NT>*[splits];
263
264 #ifdef _OPENMP
265 #pragma omp parallel for
266 #endif
267 for(int i=0; i< splits; ++i) // i iterates over splits
268 {
269 cscarr[i] = new Csc<IT,NT>(nnzs[i],n);
270 sort(colrowpairs[i].begin(), colrowpairs[i].end()); // sort w.r.t. columns first and rows second
271 cscarr[i]->jc[0] = 0;
272 std::partial_sum(colcnts[i].begin(), colcnts[i].end(), cscarr[i]->jc+1);
273 std::copy(cscarr[i]->jc, cscarr[i]->jc+n, colcnts[i].begin()); // reuse the colcnts as "current column pointers"
274
275
276 for(IT k=0; k<nnzs[i]; ++k) // k iterates over all nonzeros
277 {
278 IT cindex = std::get<0>(colrowpairs[i][k]);
279 IT rindex = std::get<1>(colrowpairs[i][k]);
280 NT value = std::get<2>(colrowpairs[i][k]);
281
282 IT curcptr = (colcnts[i][cindex])++; // fetch the pointer and post increment
283 cscarr[i]->ir[curcptr] = rindex;
284 cscarr[i]->num[curcptr] = value;
285 }
286 }
287 }
288
289
290 template<class IT, class NT>
PrintInfo() const291 void SpCCols<IT,NT>::PrintInfo() const
292 {
293 std::cout << "m: " << m ;
294 std::cout << ", n: " << n ;
295 std::cout << ", nnz: "<< nnz ;
296
297 if(splits > 0)
298 {
299 std::cout << ", local splits: " << splits << std::endl;
300 #ifdef _OPENMP
301 if(omp_get_thread_num() == 0)
302 {
303 SubPrintInfo(cscarr[0]);
304 }
305 #endif
306 }
307 else
308 {
309 std::cout << std::endl;
310 SubPrintInfo(csc);
311 }
312 }
313
314
315
316
317 /****************************************************************************/
318 /************************* PRIVATE MEMBER FUNCTIONS *************************/
319 /****************************************************************************/
320
321
322 template <class IT, class NT>
SubPrintInfo(Csc<IT,NT> * mycsc) const323 void SpCCols<IT,NT>::SubPrintInfo(Csc<IT,NT> * mycsc) const
324 {
325 #ifdef _OPENMP
326 std::cout << "Printing for thread " << omp_get_thread_num() << std::endl;
327 #endif
328 if(m < PRINT_LIMIT && n < PRINT_LIMIT) // small enough to print
329 {
330 NT ** A = SpHelper::allocate2D<NT>(m,n);
331 for(IT i=0; i< m; ++i)
332 for(IT j=0; j<n; ++j)
333 A[i][j] = NT();
334 if(mycsc != NULL)
335 {
336 for(IT i=0; i< n; ++i)
337 {
338 for(IT j = mycsc->jc[i]; j< mycsc->jc[i+1]; ++j)
339 {
340 IT rowid = mycsc->ir[j];
341 A[rowid][i] = mycsc->num[j];
342 }
343 }
344 }
345 for(IT i=0; i< m; ++i)
346 {
347 for(IT j=0; j<n; ++j)
348 {
349 std::cout << std::setiosflags(std::ios::fixed) << std::setprecision(2) << A[i][j];
350 std::cout << " ";
351 }
352 std::cout << std::endl;
353 }
354 SpHelper::deallocate2D(A,m);
355 }
356 }
357
358
359 template <class IT, class NT>
CopyCsc(Csc<IT,NT> * source)360 inline void SpCCols<IT,NT>::CopyCsc(Csc<IT,NT> * source)
361 {
362 // source csc will be NULL if number of nonzeros is zero
363 if(source != NULL)
364 csc = new Csc<IT,NT>(*source);
365 else
366 csc = NULL;
367 }
368
369 }
370