1 /*
2
3 Copyright (C) 2008-2016 Michele Martone
4
5 This file is part of librsb.
6
7 librsb is free software; you can redistribute it and/or modify it
8 under the terms of the GNU Lesser General Public License as published
9 by the Free Software Foundation; either version 3 of the License, or
10 (at your option) any later version.
11
12 librsb is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 License for more details.
16
17 You should have received a copy of the GNU Lesser General Public
18 License along with librsb; see the file COPYING.
19 If not, see <http://www.gnu.org/licenses/>.
20
21 */
22 /* @cond INNERDOC */
23 /*!
24 * @file
25 * @author Michele Martone
26 * @brief
27 * This source file contains matrix setter functions.
28 * */
29
30 #include "rsb_internals.h"
31
32 RSB_INTERNALS_COMMON_HEAD_DECLS
33
34 #define RSB_ALLOW_MTX_UPD 1
35
rsb_do_has_coo_element_inner(const struct rsb_mtx_t * mtxAp,rsb_coo_idx_t i,rsb_coo_idx_t j)36 static const void * rsb_do_has_coo_element_inner(const struct rsb_mtx_t * mtxAp, rsb_coo_idx_t i, rsb_coo_idx_t j)
37 {
38 /*!
39 * \ingroup gr_internals
40 *
41 * FIXME: unfinished
42 *
43 * Should return a pointer to the value at (i,j), if present, and NULL if not present.
44 * Only for CSR/CSC.
45 *
46 * */
47 RSB_DEBUG_ASSERT(mtxAp);
48 RSB_DEBUG_ASSERT(!RSB_INVALID_COO_INDEX(i));
49 RSB_DEBUG_ASSERT(!RSB_INVALID_COO_INDEX(j));
50 RSB_DEBUG_ASSERT(rsb__is_css_matrix(mtxAp));
51
52 if(rsb__is_recursive_matrix(mtxAp->flags))
53 {
54 const struct rsb_mtx_t * submatrix = RSB_FIND_SUBMATRIX_CONTAINING(mtxAp,i+mtxAp->roff,j+mtxAp->coff);
55 rsb_coo_idx_t moff;
56 rsb_coo_idx_t koff;
57
58 if(!submatrix)
59 return NULL;
60 moff = submatrix->roff-mtxAp->roff;
61 koff = submatrix->coff-mtxAp->coff;
62 return rsb_do_has_coo_element_inner(submatrix,i-moff,j-koff);
63 }
64 else
65 {
66 rsb_nnz_idx_t si;
67 rsb_coo_idx_t Mi,mi;
68 const rsb_nnz_idx_t * bpntr = mtxAp->bpntr;
69 const rsb_coo_idx_t * bindx = mtxAp->bindx;
70
71 if( mtxAp->flags & RSB_FLAG_WANT_COLUMN_MAJOR_ORDER )
72 Mi = j, mi = i;
73 else
74 Mi = i, mi = j;
75
76 if(rsb__is_coo_matrix(mtxAp))
77 {
78 rsb_nnz_idx_t nnz1,nnz0,nnz = mtxAp->nnz;
79 if( mtxAp->flags & RSB_FLAG_USE_HALFWORD_INDICES)
80 {
81 // delimit the current row
82 const rsb_half_idx_t *IA = (const rsb_half_idx_t *)mtxAp->bpntr;
83 const rsb_half_idx_t *JA = (const rsb_half_idx_t *)mtxAp->bindx;
84 // we search the beginning of line Mi
85 nnz0 = rsb__nnz_split_hcoo_bsearch(IA,Mi,nnz);
86 // we search the end of line Mi
87 nnz1 = nnz0+rsb__nnz_split_hcoo_bsearch(IA+nnz0,Mi+1,nnz-nnz0);
88 if(nnz1-nnz0<1)
89 return NULL;// no row Mi
90 // in line Mi, we search the index of mi, if any
91 nnz0 += rsb__nnz_split_hcoo_bsearch(JA+nnz0,mi,nnz1-nnz0);
92 if(nnz1-nnz0<1 || JA[nnz0]!=mi)
93 return NULL;// no element mi
94 }
95 else
96 {
97 // delimit the current row
98 const rsb_coo_idx_t *IA = mtxAp->bpntr;
99 const rsb_coo_idx_t *JA = mtxAp->bindx;
100 // we search the beginning of line Mi
101 nnz0 = rsb__nnz_split_coo_bsearch(IA,Mi,nnz);
102 // we search the end of line Mi
103 nnz1 = nnz0+rsb__nnz_split_coo_bsearch(IA+nnz0,Mi+1,nnz-nnz0);
104 if(nnz1-nnz0<1)
105 return NULL;// no row Mi
106 // in line Mi, we search the index of mi, if any
107 nnz0 += rsb__nnz_split_coo_bsearch(JA+nnz0,mi,nnz1-nnz0);
108 if(nnz1-nnz0<1 || JA[nnz0]!=mi)
109 return NULL;// no element mi
110 }
111 return ((const rsb_char_t*)(mtxAp->VA))+mtxAp->el_size*(nnz0);
112 }
113 else
114 {
115 if(mtxAp->flags & RSB_FLAG_USE_HALFWORD_INDICES_CSR)
116 si = rsb__seek_half_idx_t(((rsb_half_idx_t*)bindx)+bpntr[Mi],mi,bpntr[Mi+1]-bpntr[Mi]);
117 else
118 si = rsb__seek_nnz_idx_t(bindx+bpntr[Mi],mi,bpntr[Mi+1]-bpntr[Mi]);
119 }
120
121 if(si == RSB_MARKER_NNZ_VALUE)
122 return NULL;
123 else
124 return ((const rsb_char_t*)(mtxAp->VA))+mtxAp->el_size*(bpntr[Mi]+si);
125 }
126 }
127
rsb__do_coo_element_inner_address(const struct rsb_mtx_t * mtxAp,rsb_coo_idx_t i,rsb_coo_idx_t j)128 const void * rsb__do_coo_element_inner_address(const struct rsb_mtx_t * mtxAp, rsb_coo_idx_t i, rsb_coo_idx_t j)
129 {
130 return rsb_do_has_coo_element_inner(mtxAp,i,j);
131 }
132
rsb__do_set_coo_elements(struct rsb_mtx_t * mtxAp,const void * VA,const rsb_coo_idx_t * IA,const rsb_coo_idx_t * JA,rsb_nnz_idx_t nnz)133 rsb_err_t rsb__do_set_coo_elements(struct rsb_mtx_t * mtxAp, const void * VA, const rsb_coo_idx_t *IA, const rsb_coo_idx_t *JA, rsb_nnz_idx_t nnz)
134 {
135 /*!
136 * \ingroup gr_internals
137 * FIXME: undocumented
138 * FIXME: should parallelize meaningfully
139 * */
140 rsb_err_t errval = RSB_ERR_NO_ERROR;
141 rsb_nnz_idx_t n;
142
143 if(nnz<rsb_global_session_handle.rsb_want_threads)
144 {
145 for(n=0;RSB_LIKELY(n<nnz);++n)
146 {
147 RSB_DO_ERROR_CUMULATE(errval,rsb__do_set_coo_element(mtxAp,((rsb_char_t*)VA)+mtxAp->el_size*n,IA[n],JA[n]));
148 if(RSB_SOME_ERROR(errval))
149 {
150 RSB_ERROR("error updating %dth element of %d: %d %d\n",n,nnz,IA[n],JA[n]);
151 RSB_PERR_GOTO(err,RSB_ERRM_ES)
152 }
153 }
154 }
155 else
156 {
157 #pragma omp parallel for schedule(static,1) reduction(|:errval) RSB_NTC
158 for(n=0;n<nnz;++n)
159 {
160 RSB_DO_ERROR_CUMULATE(errval,rsb__do_set_coo_element(mtxAp,((rsb_char_t*)VA)+mtxAp->el_size*n,IA[n],JA[n]));
161 if(RSB_SOME_ERROR(errval))
162 {
163 RSB_ERROR("error updating %dth element of %d: %d %d\n",n,nnz,IA[n],JA[n]);
164 // RSB_PERR_GOTO(err,RSB_ERRM_ES)
165 }
166 }
167 }
168 err:
169 RSB_DO_ERR_RETURN(errval)
170 }
171
rsb__do_set_coo_element(struct rsb_mtx_t * mtxAp,const void * vp,const rsb_coo_idx_t i,const rsb_coo_idx_t j)172 rsb_err_t rsb__do_set_coo_element(struct rsb_mtx_t * mtxAp, const void * vp, const rsb_coo_idx_t i, const rsb_coo_idx_t j)
173 {
174 return rsb__do_upd_coo_element(mtxAp, vp, i, j, RSB_FLAG_DUPLICATES_DEFAULT_HANDLE);
175 }
176
rsb__do_upd_coo_element(struct rsb_mtx_t * mtxAp,const void * vp,const rsb_coo_idx_t i,const rsb_coo_idx_t j,rsb_flags_t flags)177 rsb_err_t rsb__do_upd_coo_element(struct rsb_mtx_t * mtxAp, const void * vp, const rsb_coo_idx_t i, const rsb_coo_idx_t j, rsb_flags_t flags)
178 {
179 /*!
180 * \ingroup gr_internals
181 *
182 * Overwrites the element at (i,j), if present.
183 * If not present, returns RSB_ERR_GENERIC_ERROR.
184 * */
185 rsb_err_t errval = RSB_ERR_NO_ERROR;
186 void * OV = NULL;
187
188 if(!mtxAp || !vp || RSB_INVALID_COO_INDEX(i) || RSB_INVALID_COO_INDEX(j))
189 {
190 errval = RSB_ERR_BADARGS;
191 RSB_PERR_GOTO(err,RSB_ERRM_ES)
192 }
193
194 if(!rsb__is_css_matrix(mtxAp))
195 {
196 errval = RSB_ERR_UNIMPLEMENTED_YET;
197 RSB_PERR_GOTO(err,RSB_ERRM_ES)
198 }
199
200 if(!RSB_MATRIX_CONTAINS(mtxAp,i,j))
201 {
202 errval = RSB_ERR_BADARGS;
203 RSB_PERR_GOTO(err,RSB_ERRM_ES)
204 }
205
206 if( i == j && rsb__get_diagonal_type_flag(mtxAp)==RSB_DIAGONAL_I )
207 {
208 errval = RSB_ERR_BADARGS;
209 RSB_PERR_GOTO(err,RSB_ERRM_ES)
210 }
211
212 OV = (void*) rsb_do_has_coo_element_inner(mtxAp,i,j);
213
214 if(OV)
215 {
216 #if RSB_ALLOW_MTX_UPD
217 rsb_flags_t cflag = ( mtxAp->flags & RSB_FLAG_ALL_DUPLICATE_FLAGS )
218 | ( flags & RSB_FLAG_ALL_DUPLICATE_FLAGS );
219
220 if(RSB_DO_FLAG_HAS(cflag,RSB_FLAG_DUPLICATES_SUM))
221 { RSB_NUMERICAL_TYPE_SUM_AND_STORE_ELEMENTS(OV,vp,mtxAp->typecode);}
222 else
223 #endif
224 { RSB_NUMERICAL_TYPE_SET_ELEMENT(OV,vp,mtxAp->typecode); }
225 }
226 else
227 errval = RSB_ERR_GENERIC_ERROR;
228
229 err:
230 RSB_DO_ERR_RETURN(errval)
231 }
232
233 #if 1
234 /* rsb_do_locate_nnz_element and rsb_do_get_nnz_element are used (experimentally) by sparsersb */
rsb_do_locate_nnz_element(const struct rsb_mtx_t * mtxAp,void ** vpp,rsb_coo_idx_t * ip,rsb_coo_idx_t * jp,rsb_nnz_idx_t nzi)235 static rsb_err_t rsb_do_locate_nnz_element(const struct rsb_mtx_t * mtxAp, void ** vpp, rsb_coo_idx_t*ip, rsb_coo_idx_t*jp, rsb_nnz_idx_t nzi)
236 {
237 /* FIXME: * new, unfinished, untested */
238 rsb_err_t errval = RSB_ERR_BADARGS;
239 rsb_submatrix_idx_t i,j;
240 struct rsb_mtx_t * submatrix = NULL;
241
242 if(!mtxAp)
243 {
244 RSB_PERR_GOTO(err,RSB_ERRM_ES)
245 }
246
247 if(rsb__is_terminal_recursive_matrix(mtxAp)
248 && ( nzi >= mtxAp->nzoff )
249 && ( nzi < mtxAp->nzoff+mtxAp->nnz )
250 )
251 {
252 rsb_nnz_idx_t lnz = nzi-mtxAp->nzoff;
253 rsb_byte_t*OV = mtxAp->VA;
254 OV += mtxAp->el_size * lnz;
255
256 if( (!ip) || (!jp) )
257 goto noij;
258
259 if(rsb__is_coo_matrix(mtxAp))
260 {
261 if(RSB_DO_FLAG_HAS(mtxAp->flags,RSB_FLAG_USE_HALFWORD_INDICES))
262 {
263 RSB_DECLARE_CONST_HALFCOO_ARRAYS_FROM_MATRIX(mIA,mJA,mtxAp)
264 i = mIA[lnz];
265 j = mJA[lnz];
266 }
267 else
268 {
269 RSB_DECLARE_CONST_FULLCOO_ARRAYS_FROM_MATRIX(mIA,mJA,mtxAp)
270 i = mIA[lnz];
271 j = mJA[lnz];
272 }
273 }
274 else
275 {
276 if(mtxAp->flags & RSB_FLAG_USE_HALFWORD_INDICES_CSR)
277 {
278 RSB_DECLARE_CONST_HALFCSR_ARRAYS_FROM_MATRIX(mPA,mJA,mtxAp);
279 j = mJA[lnz];
280 i = rsb__nnz_split_coo_bsearch(mPA,lnz,mtxAp->nnz);
281 }
282 else
283 {
284 RSB_DECLARE_CONST_FULLCSR_ARRAYS_FROM_MATRIX(mPA,mJA,mtxAp);
285 j = mJA[lnz];
286 i = rsb__nnz_split_coo_bsearch(mPA,lnz,mtxAp->nnz);
287 }
288 }
289 if(ip)*ip = i;
290 if(jp)*jp = j;
291 noij:
292 if(vpp)*vpp = OV;
293 errval = RSB_ERR_NO_ERROR;
294 }
295 else
296 {
297 RSB_SUBMATRIX_FOREACH(mtxAp,submatrix,i,j)
298 if(submatrix
299 && ( nzi >= submatrix->nzoff )
300 && ( nzi < submatrix->nnz+submatrix->nzoff )
301 )
302 {
303 errval = rsb_do_locate_nnz_element(submatrix,vpp,ip,jp,nzi);
304 RSB_PERR_GOTO(err,RSB_ERRM_ES)
305 }
306 }
307 err:
308 RSB_DO_ERR_RETURN(errval)
309 }
310
rsb_do_get_nnz_element(const struct rsb_mtx_t * mtxAp,void * vp,rsb_coo_idx_t * ip,rsb_coo_idx_t * jp,rsb_nnz_idx_t nzi)311 rsb_err_t rsb_do_get_nnz_element(const struct rsb_mtx_t * mtxAp, void * vp, rsb_coo_idx_t*ip, rsb_coo_idx_t*jp, rsb_nnz_idx_t nzi)
312 {
313 /* FIXME: * new, unfinished (20130331) */
314 void * OV = NULL;
315 rsb_err_t errval = RSB_ERR_NO_ERROR;
316 errval = rsb_do_locate_nnz_element(mtxAp,&OV,ip,jp,nzi);
317 if((!RSB_SOME_ERROR(errval)) && OV)
318 {
319 RSB_NUMERICAL_TYPE_SET_ELEMENT(vp,OV,mtxAp->typecode);
320 }
321 RSB_DO_ERR_RETURN(errval)
322 }
323 #endif
324
325 #if 0
326 rsb_err_t rsb_do_set_nnz_element(const struct rsb_mtx_t * mtxAp, const void * vp, rsb_nnz_idx_t nzi)
327 {
328 /* FIXME: * new, unfinished (20130331) */
329 void * OV = NULL;
330 rsb_err_t errval = RSB_ERR_NO_ERROR;
331 errval = rsb_do_locate_nnz_element(mtxAp,&OV,NULL,NULL,nzi);
332 if((!RSB_SOME_ERROR(errval)) && OV)
333 {
334 RSB_NUMERICAL_TYPE_SET_ELEMENT(OV,vp,mtxAp->typecode);
335 }
336 RSB_DO_ERR_RETURN(errval)
337 }
338 #endif
339
rsb__do_get_coo_element(const struct rsb_mtx_t * mtxAp,void * vp,rsb_coo_idx_t i,rsb_coo_idx_t j)340 rsb_err_t rsb__do_get_coo_element(const struct rsb_mtx_t * mtxAp, void * vp, rsb_coo_idx_t i, rsb_coo_idx_t j)
341 {
342 /*!
343 * \ingroup gr_internals
344 *
345 * FIXME: undocumented
346 * Gets the element at (i,j), if present.
347 * If not present, returns RSB_ERR_GENERIC_ERROR and zeros the area.
348 * */
349 const void * OV = NULL;
350 rsb_err_t errval = RSB_ERR_NO_ERROR;
351
352 #if RSB_ALLOW_ZERO_DIM
353 if(RSB_ANY_MTX_DIM_ZERO(mtxAp))
354 goto err; /* FIXME: skipping further error checks */
355 #endif
356 if(!mtxAp || !vp || RSB_INVALID_COO_INDEX(i) || RSB_INVALID_COO_INDEX(j))
357 {
358 errval = RSB_ERR_BADARGS;
359 RSB_PERR_GOTO(err,RSB_ERRM_ES)
360 }
361
362 if(!rsb__is_css_matrix(mtxAp))
363 {
364 errval = RSB_ERR_UNIMPLEMENTED_YET;
365 RSB_PERR_GOTO(err,RSB_ERRM_ES)
366 }
367
368 if(!RSB_MATRIX_CONTAINS(mtxAp,i,j))
369 {
370 errval = RSB_ERR_BADARGS;
371 RSB_PERR_GOTO(err,RSB_ERRM_ES)
372 }
373
374 if( i == j && rsb__get_diagonal_type_flag(mtxAp)==RSB_DIAGONAL_I )
375 return rsb__fill_with_ones(vp,mtxAp->typecode,1,1);
376
377 OV = rsb_do_has_coo_element_inner(mtxAp,i,j);
378 if(!OV)
379 {
380 errval = RSB_ERR_GENERIC_ERROR;
381 rsb__cblas_Xscal(mtxAp->typecode,1,NULL,vp,1);
382 }
383 else
384 {RSB_NUMERICAL_TYPE_SET_ELEMENT(vp,OV,mtxAp->typecode);}
385
386 err:
387 RSB_DO_ERR_RETURN(errval)
388 }
389
rsb__do_reverse_odd_rows(struct rsb_mtx_t * mtxAp)390 rsb_err_t rsb__do_reverse_odd_rows(struct rsb_mtx_t * mtxAp)
391 {
392 /*!
393 * \ingroup gr_internals
394 * \rsb_warn_unoptimized_msg
395 * \note Works only for csr leaves.
396 * */
397 RSB_DEBUG_ASSERT(mtxAp);
398
399 if(rsb__is_recursive_matrix(mtxAp->flags))
400 {
401 struct rsb_mtx_t * submatrix;
402 rsb_submatrix_idx_t i,j;
403 RSB_SUBMATRIX_FOREACH(mtxAp,submatrix,i,j)
404 rsb__do_reverse_odd_rows(submatrix);
405 }
406 else
407 {
408 if(rsb__is_coo_matrix(mtxAp))
409 {
410 if( mtxAp->flags & RSB_FLAG_USE_HALFWORD_INDICES)
411 {
412 //RSB_DECLARE_CONST_HALFCOO_ARRAYS_FROM_MATRIX(IA,JA,mtxAp)
413 /* FIXME: not implemented yet for COO */
414 }
415 else
416 {
417 //RSB_DECLARE_CONST_FULLCOO_ARRAYS_FROM_MATRIX(IA,JA,mtxAp)
418 /* FIXME: not implemented yet for COO */
419 }
420 }
421 else
422 {
423 rsb_coo_idx_t i;
424 rsb_nnz_idx_t ib,ie;
425 if(mtxAp->flags & RSB_FLAG_USE_HALFWORD_INDICES)
426 {
427 RSB_DECLARE_HALFCSR_ARRAYS_FROM_MATRIX(IA,JA,mtxAp)
428 for(i=1;i<mtxAp->nr;i+=2)
429 ib = IA[i],ie = IA[i+1],
430 rsb__util_reverse_halfword_coo_array(JA+ib,ie-ib);
431 }
432 else
433 {
434 RSB_DECLARE_FULLCSR_ARRAYS_FROM_MATRIX(IA,JA,mtxAp)
435 for(i=1;i<mtxAp->nr;i+=2)
436 ib = IA[i],ie = IA[i+1],
437 rsb__util_reverse_fullword_coo_array(JA+ib,ie-ib);
438 }
439 }
440 }
441 return RSB_ERR_NO_ERROR;
442 }
443
rsb__do_zsort_coo_submatrices(struct rsb_mtx_t * mtxAp)444 rsb_err_t rsb__do_zsort_coo_submatrices(struct rsb_mtx_t * mtxAp)
445 {
446 /*!
447 * \ingroup gr_internals
448 * \rsb_warn_unoptimized_msg
449 * */
450 RSB_DEBUG_ASSERT(mtxAp);
451
452 if(rsb__is_recursive_matrix(mtxAp->flags))
453 {
454 struct rsb_mtx_t * submatrix;
455 rsb_submatrix_idx_t i,j;
456 RSB_SUBMATRIX_FOREACH(mtxAp,submatrix,i,j)
457 if(submatrix)
458 rsb__do_zsort_coo_submatrices(submatrix);
459 }
460 else
461 {
462 if(rsb__is_coo_matrix(mtxAp))
463 {
464 rsb_err_t errval = RSB_ERR_NO_ERROR;
465 if( mtxAp->flags & RSB_FLAG_USE_HALFWORD_INDICES)
466 {
467 RSB_DECLARE_FULLCOO_ARRAYS_FROM_MATRIX(IA,JA,mtxAp)
468 RSB_DO_ERROR_CUMULATE(errval,rsb__do_switch_to_fullword_coo(mtxAp));
469 RSB_DO_ERROR_CUMULATE(errval,rsb__do_index_based_z_morton_sort(NULL,NULL,NULL,IA,JA,mtxAp->VA,mtxAp->nr,mtxAp->nc,mtxAp->nnz,mtxAp->typecode,RSB_OP_FLAG_DEFAULT));
470 RSB_DO_ERROR_CUMULATE(errval,rsb__do_switch_to_halfword_coo(mtxAp));
471 }
472 else
473 {
474 RSB_DECLARE_FULLCOO_ARRAYS_FROM_MATRIX(IA,JA,mtxAp)
475 RSB_DO_ERROR_CUMULATE(errval,rsb__do_index_based_z_morton_sort(NULL,NULL,NULL,IA,JA,mtxAp->VA,mtxAp->nr,mtxAp->nc,mtxAp->nnz,mtxAp->typecode,RSB_OP_FLAG_DEFAULT));
476 }
477 return errval;
478 }
479 else
480 {
481 if(mtxAp->flags & RSB_FLAG_USE_HALFWORD_INDICES)
482 {
483 /* nothing to do */
484 }
485 else
486 {
487 /* nothing to do */
488 }
489 }
490 }
491 return RSB_ERR_NO_ERROR;
492 }
493
494 /* @endcond */
495