1 /*
2
3 Copyright (C) 2008-2021 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 Sparse BLAS interface testing code
27 * */
28 #include "rsb_common.h"
29 #include "blas_sparse/blas_enum.h"
30 #include "rsb_libspblas.h"
31 /* #include "rsb_libspblas_handle.h" */
32 #include "rsb_psblas.h" /* (in rsb) header for rsb__do_psblas_trans_to_rsb_trans */
33 #include <stdio.h> /* fileno */
34 #include <unistd.h> /* isatty */
35 #include "rsb_libspblas_tests.h"
36 #define RSB_WANT_SPGEMM_TESTING_FOR_ONLY_FIRST_DIMI 3
37 #define RSB_WANT_VERBOSE_FAILURES 1
38 #define RSB_TESTER_ALLOW_TIMEOUT 1
39 #define RSB_BLAS_INVALID_MATRIX (-1)
40 #define RSB_INVALID_BLAS_INT_IDX_VAL -1
41 #define RSB_WANT_AUTOTUNING_TESTING 1
42
43 RSB_INTERNALS_COMMON_HEAD_DECLS
44 RSB_INTERNALS_RSBENCH_HEAD_DECLS
45
46 /* TODO: shall use the following throughout the tester routine */
47 #define RSB_LSTERR(MSG) {RSB_ERROR(MSG);goto err;}
48 #define RSB_LSTPROBE(EXP,MSG) if( RSB_SOME_ERROR(errval=(EXP))){RSB_ERROR(MSG);goto err;} /* error is not expected here */
49 #define RSB_LSTPROBI(EXP,MSG) if(!RSB_SOME_ERROR(errval=(EXP))){errval=RSB_ERR_NO_ERROR;RSB_ERROR(MSG);goto err;}else{errval = RSB_ERR_INTERNAL_ERROR;} /* error is expected here but ignore: better use it after a non-zealot-error-reporting function, but rather, an internal one */
50 #define RSB_EMPTY_STRING ""
51
52 #define RSB_BLAS_MT_STR(MT) ((MT)==blas_lower_triangular?"LT": \
53 ((MT)==blas_upper_triangular?"UT": \
54 ((MT)==blas_lower_symmetric?"LS": \
55 ((MT)==blas_upper_symmetric?"US": \
56 ((MT)==blas_lower_hermitian?"LH": \
57 ((MT)==blas_upper_hermitian?"UH": \
58 ((MT)==blas_general?"GE":"??") \
59 ))))))
60
rsb_blas_tester_options_init(struct rsb_tester_options_t * top)61 rsb_err_t rsb_blas_tester_options_init(struct rsb_tester_options_t * top)
62 {
63 /* This function shall not need any library initialization to work. */
64 rsb_err_t errval = RSB_ERR_NO_ERROR;
65 RSB_BZERO_P(top);
66 top->mtt = RSB_TIME_ZERO;
67 top->rrm = RSB_BOOL_FALSE;
68 top->tur = RSB_BOOL_FALSE;
69 top->wqt = RSB_BOOL_FALSE;
70 top->wqc = RSB_BOOL_FALSE;
71 top->wcs = RSB_BOOL_FALSE;
72 return errval;
73 }
74
rsb_blas_single_allocation_tester(void)75 static blas_sparse_matrix rsb_blas_single_allocation_tester(void)
76 {
77 /**
78 * \ingroup gr_internals
79 * */
80 blas_sparse_matrix A = RSB_BLAS_INVALID_MATRIX;
81 const rsb_coo_idx_t IA[]={0};
82 const rsb_coo_idx_t JA[]={0};
83 const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE VA[]={11};
84 const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE X[]={1};
85 RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE Y[]={0};
86 const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE alpha = 1.0, beta = 1.0;
87 const rsb_coo_idx_t m=1,n=1;
88 const int nz=1;
89 rsb_type_t typecode = RSB_NUMERICAL_TYPE_FIRST_BLAS;
90 if( RSB_NUMERICAL_TYPE_FIRST_BLAS == RSB_NUMERICAL_TYPE_INVALID_TYPE )
91 { RSB_INFO("SKIPPING A TEST (no BLAS types in)\n"); goto err; }
92 if( (A= rsb__BLAS_Xuscr_begin( m, n, typecode )) == RSB_BLAS_INVALID_MATRIX )
93 {RSB_ERROR("error calling BLAS_duscr_begin\n"); goto err;}
94 if( rsb__BLAS_Xuscr_insert_entries( A, nz, VA, IA, JA) == RSB_BLAS_ERROR )
95 {RSB_ERROR("error calling BLAS_duscr_insert_entries\n"); goto err;}
96 if( rsb__BLAS_Xuscr_end(A) == RSB_BLAS_ERROR )
97 {RSB_ERROR("error calling BLAS_duscr_end\n"); goto err;}
98 if( rsb__BLAS_Xusmv( blas_no_trans, &alpha, A, X, 1, &beta, Y, 1) == RSB_BLAS_ERROR)
99 {RSB_ERROR("error calling BLAS_dusmv\n"); goto err;}
100
101 return A;
102 err:
103 if(A != RSB_BLAS_INVALID_MATRIX)
104 {
105 if (rsb__BLAS_Xusds(A)!=RSB_BLAS_NO_ERROR)
106 {RSB_ERROR("error destroying the matrix after error\n"); goto err;}
107
108 }
109 return RSB_BLAS_INVALID_MATRIX;
110 }
111
112
rsb_blas_allocation_tester(void)113 static rsb_err_t rsb_blas_allocation_tester(void)
114 {
115 /**
116 * \ingroup gr_internals
117 * Descriptor handling machinery tester.
118 * */
119 rsb_err_t errval = RSB_ERR_NO_ERROR;
120 const rsb_submatrix_idx_t mcount = RSB_MIN(1024,RSB_BLAS_MATRICES_MAX);
121 rsb_submatrix_idx_t count=0;
122 blas_sparse_matrix bsms[mcount];
123 blas_sparse_matrix A = RSB_BLAS_INVALID_MATRIX;
124
125 if( RSB_NUMERICAL_TYPE_FIRST_BLAS == RSB_NUMERICAL_TYPE_INVALID_TYPE )
126 { RSB_INFO("SKIPPING A TEST (no BLAS types in)\n"); goto err; }
127
128 for(count=0;count<mcount;++count)
129 bsms[count] = RSB_BLAS_INVALID_MATRIX;
130
131 for(count=0;count<mcount;++count)
132 {
133 A = rsb_blas_single_allocation_tester();
134 if(A == RSB_BLAS_INVALID_MATRIX)
135 {
136 RSB_ERROR(RSB_ERRM_ES);
137 errval = RSB_ERR_INTERNAL_ERROR;
138 goto out;
139 }
140 else
141 bsms[count]=A;
142 }
143 out:
144 if(count<mcount)
145 RSB_ERROR("failed allocating %d matrices: only allocated %d!\n",mcount,count);
146
147 for(count=mcount-1;count+1>0;--count)
148 {
149 if((bsms[count] != RSB_BLAS_INVALID_MATRIX) && (rsb__BLAS_Xusds(bsms[count])==RSB_BLAS_ERROR))
150 {
151 RSB_ERROR(RSB_ERRM_ES);
152 RSB_DO_ERROR_CUMULATE(errval,RSB_ERR_INTERNAL_ERROR);
153 }
154 bsms[count] = RSB_BLAS_INVALID_MATRIX;
155 }
156 err:
157 RSB_DO_ERR_RETURN(errval)
158 }
159
rsb_blas_mini_tester(void)160 rsb_err_t rsb_blas_mini_tester(void)
161 {
162 /**
163 * \ingroup gr_internals
164 * */
165 blas_sparse_matrix A = RSB_BLAS_INVALID_MATRIX;
166 const rsb_coo_idx_t IA[]={0,1,2,3};
167 const rsb_coo_idx_t JA[]={0,1,2,3};
168 const rsb_coo_idx_t BR[]={2};
169 const rsb_coo_idx_t BC[]={2};
170 const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE VA[]={0,11,22,33};
171 const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE X[]={4,3,2,1};
172 RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE Y[]={0,0,0,0};
173 const RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE alpha = 1.0, beta = 1.0;
174 const rsb_coo_idx_t m=4,n=4;
175 const int nz=4;
176 rsb_char_t optstr[RSB_MAX_LINE_LENGTH];
177 rsb_err_t errval = RSB_ERR_NO_ERROR;
178 rsb_type_t typecode = RSB_NUMERICAL_TYPE_FIRST_BLAS;
179 /* const char*tsep="*\n"; */
180 const char*tsep="%s";
181
182 if( RSB_NUMERICAL_TYPE_FIRST_BLAS == RSB_NUMERICAL_TYPE_INVALID_TYPE )
183 { RSB_INFO("SKIPPING A TEST (no BLAS types in)\n"); goto ret; }
184
185 #ifndef RSB_NUMERICAL_TYPE_DOUBLE
186 RSB_INFO("SKIPPING BASIC SPARSE BLAS TEST (UNFINISHED TESTING SUITE)\n");
187 goto ret; /* FIXME: we are being overly tolerant here, because we don't want to break int-type-only cases */
188 #endif /* RSB_NUMERICAL_TYPE_DOUBLE */
189
190 #if 1
191 RSB_INFO("BASIC SPARSE BLAS TEST: BEGIN\n");
192 if((A= rsb__BLAS_Xuscr_begin( m, n, typecode )) == RSB_BLAS_INVALID_MATRIX )
193 {RSB_ERROR(RSB_ERRM_NL); goto err;}
194 RSB_INFO(tsep,"");
195 if( rsb__BLAS_Xuscr_insert_entries( A, nz, VA, IA, JA) == RSB_BLAS_INVALID_MATRIX )
196 {RSB_ERROR(RSB_ERRM_NL); goto err;}
197 RSB_INFO(tsep,"");
198 if( rsb__BLAS_Xuscr_end(A) == RSB_BLAS_INVALID_MATRIX )
199 {RSB_ERROR(RSB_ERRM_NL); goto err;}
200 if( rsb__BLAS_usgp( A, blas_num_rows ) != m ) {RSB_ERROR(RSB_ERRM_NL); goto err;}
201 if( rsb__BLAS_usgp( A, blas_num_cols ) != n ) {RSB_ERROR(RSB_ERRM_NL); goto err;}
202 if( rsb__BLAS_usgp( A, blas_num_nonzeros ) != nz ) {RSB_ERROR(RSB_ERRM_NL); goto err;}
203 RSB_INFO(tsep,"");
204 if( rsb__BLAS_Xusmv( blas_no_trans, &alpha, A, X, 1, &beta, Y, 1) != RSB_BLAS_NO_ERROR )
205 {RSB_ERROR(RSB_ERRM_NL); goto err;}
206 RSB_INFO(tsep,"");
207 #if 0
208 /* missing lower triangular mark */
209 if( BLAS_dussv( blas_no_trans, alpha, A, X, 1) != RSB_BLAS_NO_ERROR )
210 goto err;
211 RSB_INFO("*\n");
212 #endif
213
214
215 #if RSB_WANT_ALLOCATOR_LIMITS
216 if(1)
217 {
218 /* TODO: in the future, may constrain the whole test within memory limits */
219 size_t sval=0;
220 sval=0;
221 RSB_DO_REINIT_SINGLE_VALUE_SET(RSB_IO_WANT_MAX_MEMORY_ALLOCATIONS,&sval,errval); RSB_LSTPROBE(errval,"");
222 RSB_DO_REINIT_SINGLE_VALUE_SET(RSB_IO_WANT_MAX_MEMORY_ALLOCATED,&sval,errval); RSB_LSTPROBE(errval,"");
223 }
224 #endif /* RSB_WANT_ALLOCATOR_LIMITS */
225 if(1)
226 {
227 rsb_int val=0;
228 enum rsb_opt_t key=RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE;
229
230 RSB_INFO("INIT INTERFACE TEST: BEGIN\n");
231 RSB_DO_REINIT_SINGLE_VALUE_GET(key,&val,errval); RSB_LSTPROBE(errval,"");
232 if(val!=-1)
233 { RSB_DO_REINIT_SINGLE_VALUE_SET(key,&val,errval); }
234 RSB_LSTPROBE(errval,"");
235 rsb__sprintf(optstr,"got RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE: %d",val);
236 if(val!=-1)
237 { RSB_LSTPROBE(rsb__do_set_initopt_as_string("RSB_IO_WANT_EXTRA_VERBOSE_INTERFACE",optstr),""); }
238 RSB_INFO("%s\n",optstr);
239
240 key=RSB_IO_WANT_IS_INITIALIZED_MARKER;
241 RSB_DO_REINIT_SINGLE_VALUE_GET(key,&val,errval); RSB_LSTPROBE(errval,"");
242 RSB_DO_REINIT_SINGLE_VALUE_SET(key,&val,errval); RSB_LSTPROBE(errval,"");
243 rsb__sprintf(optstr,"%d",val);
244 RSB_LSTPROBE(rsb__do_set_initopt_as_string("RSB_IO_WANT_IS_INITIALIZED_MARKER",optstr),"");
245 RSB_INFO("got RSB_IO_WANT_IS_INITIALIZED_MARKER: %s\n",optstr);
246
247 RSB_INFO("INIT INTERFACE TEST: END (SUCCESS)\n");
248 }
249 RSB_INFO("PRINT TEST: BEGIN\n");
250 errval = rsb__do_file_mtx_save(rsb__BLAS_inner_matrix_retrieve(A),NULL);
251 /* rsb_mtx_file_render(rsb__BLAS_inner_matrix_retrieve(A)); */
252 errval = rsb__debug_print_vectors_diff(VA,VA+1,nz-1,typecode,1,1,RSB_VECTORS_DIFF_DISPLAY_N);
253 errval = rsb__do_mtx_render(NULL,rsb__BLAS_inner_matrix_retrieve(A), 100, 100, RSB_MARF_EPS | RSB_MARF_EPS_B);
254 /* rsb__do_file_mtx_rndr(void * pmp, const char * filename, rsb_coo_idx_t pmlWidth, rsb_coo_idx_t pmWidth, rsb_coo_idx_t pmHeight, rsb_marf_t rflags) */
255 RSB_INFO("PRINT TEST: END (SUCCESS)\n");
256
257 if(rsb__BLAS_Xusds( A ) == RSB_BLAS_ERROR)
258 {RSB_ERROR(RSB_ERRM_NL); goto err;}
259 RSB_INFO(tsep,"");
260 #endif
261 #if 1
262 if((A= rsb__BLAS_Xuscr_block_begin( 1, 1, 2, 2, typecode )) == RSB_BLAS_INVALID_MATRIX )
263 {RSB_ERROR(RSB_ERRM_NL); goto err;}
264 RSB_INFO(tsep,"");
265 if( rsb__BLAS_Xuscr_insert_block( A, VA, 1, 1, 0, 0) == RSB_BLAS_INVALID_MATRIX )
266 {RSB_ERROR(RSB_ERRM_NL); goto err;}
267 RSB_INFO(tsep,"");
268 if( rsb__BLAS_Xuscr_end(A) == RSB_BLAS_INVALID_MATRIX )
269 {RSB_ERROR(RSB_ERRM_NL); goto err;}
270
271 if( rsb__BLAS_usgp( A, blas_num_rows ) != 2 ) {RSB_ERROR(RSB_ERRM_NL); goto err;}
272 if( rsb__BLAS_usgp( A, blas_num_cols ) != 2 ) {RSB_ERROR(RSB_ERRM_NL); goto err;}
273 if( rsb__BLAS_usgp( A, blas_num_nonzeros ) != 4 ) {RSB_ERROR(RSB_ERRM_NL); goto err;}
274
275 RSB_INFO(tsep,"");
276 if( rsb__BLAS_Xusmv( blas_no_trans, &alpha, A, X, 1, &beta, Y, 1) != RSB_BLAS_NO_ERROR )
277 {RSB_ERROR(RSB_ERRM_NL); goto err;}
278 RSB_INFO(tsep,"");
279 if(rsb__BLAS_Xusds( A ) == RSB_BLAS_ERROR)
280 {RSB_ERROR(RSB_ERRM_NL); goto err;}
281 RSB_INFO(tsep,"");
282 #endif
283 #if 1
284 if((A= rsb__BLAS_Xuscr_variable_block_begin( 1, 1, &BR[0], &BC[0], typecode )) == RSB_BLAS_INVALID_MATRIX )
285 {RSB_ERROR(RSB_ERRM_NL); goto err;}
286 RSB_INFO(tsep,"");
287 if( rsb__BLAS_Xuscr_insert_block( A, VA, 1, 1, 0, 0) == RSB_BLAS_INVALID_MATRIX )
288 {RSB_ERROR(RSB_ERRM_NL); goto err;}
289 RSB_INFO(tsep,"");
290 if( rsb__BLAS_Xuscr_end(A) == RSB_BLAS_INVALID_MATRIX )
291 {RSB_ERROR(RSB_ERRM_NL); goto err;}
292 if( rsb__BLAS_usgp( A, blas_num_rows ) != 2 ) {RSB_ERROR(RSB_ERRM_NL); goto err;}
293 if( rsb__BLAS_usgp( A, blas_num_cols ) != 2 ) {RSB_ERROR(RSB_ERRM_NL); goto err;}
294 if( rsb__BLAS_usgp( A, blas_num_nonzeros ) != 4 ) {RSB_ERROR(RSB_ERRM_NL); goto err;}
295 RSB_INFO(tsep,"");
296 if( rsb__BLAS_Xusmv( blas_no_trans, &alpha, A, X, 1, &beta, Y, 1) != RSB_BLAS_NO_ERROR )
297 {RSB_ERROR(RSB_ERRM_NL); goto err;}
298 RSB_INFO(tsep,"");
299 if(rsb__BLAS_Xusds( A ) == RSB_BLAS_ERROR)
300 {RSB_ERROR(RSB_ERRM_NL); goto err;}
301 RSB_INFO(tsep,"");
302 #endif
303 RSB_INFO("BASIC SPARSE BLAS TEST: END (SUCCESS)\n");
304
305 #if 0
306 RSB_INFO("BIGGER MATRICES SPARSE BLAS TEST: BEGIN\n");
307 if(rsb_blas_bigger_matrices_tester(NULL))
308 goto err;
309 RSB_INFO("BIGGER MATRICES SPARSE BLAS TEST: END\n");
310 #endif
311
312 RSB_INFO("STRESS SPARSE BLAS TEST: BEGIN\n");
313 if(RSB_SOME_ERROR(rsb_blas_allocation_tester()))
314 {RSB_ERROR(RSB_ERRM_NL); goto err;}
315 RSB_INFO("STRESS SPARSE BLAS TEST: END (SUCCESS)\n");
316
317 RSB_INFO("SPARSE BLAS TESTS: END (SUCCESS)\n");
318 goto ret;
319 err:
320 RSB_INFO("SPARSE BLAS TESTS: FAILURE!\n");
321 errval = RSB_ERR_GENERIC_ERROR;
322 ret:
323 return errval;
324 }
325
rsb_basic_primitives_tester(void)326 static rsb_err_t rsb_basic_primitives_tester(void)
327 {
328 /**
329 * \ingroup gr_internals
330 * */
331 rsb_err_t errval = RSB_ERR_NO_ERROR;
332 const size_t n=1024; // FIXME
333 rsb_nnz_idx_t i=0;
334 rsb_coo_idx_t *cp = rsb__calloc(sizeof(rsb_coo_idx_t)*n);
335 rsb_half_idx_t*hp=(rsb_half_idx_t*)cp;
336 RSB_INFO("BASIC PRIMITIVES TEST: BEGIN\n");
337 if(cp==NULL){ RSB_ERROR(RSB_ERRM_ES); errval = RSB_ERR_ENOMEM; goto err; }
338 // RSB_XCOO_ISET(hp,0,n);
339 for(i=0;i<n;++i) hp[i]=i;
340 for(i=0;i<n;++i)if(hp[i]!=i){ RSB_ERROR("half word assignment is broken");errval = RSB_ERR_INTERNAL_ERROR; goto err;}
341 rsb__do_switch_array_to_fullword_coo(hp,n,0);
342 for(i=0;i<n;++i)if(cp[i]!=i){ RSB_ERROR("half to full word conversion is broken (has %d instead of %d)",cp[i],i);errval = RSB_ERR_INTERNAL_ERROR; }
343 if(RSB_SOME_ERROR(errval))
344 goto err;
345 rsb__do_switch_array_to_halfword_coo(cp,n,0);
346 for(i=0;i<n;++i)if(hp[i]!=i){ RSB_ERROR("full to half word conversion is broken");errval = RSB_ERR_INTERNAL_ERROR; goto err;}
347
348 RSB_CONDITIONAL_FREE(cp);
349
350 err:
351 if(RSB_SOME_ERROR(errval))
352 {
353 rsb__do_perror(NULL,errval);
354 RSB_INFO("BASIC PRIMITIVES TEST: END (FAILURE)\n");
355 }
356 else
357 {
358 RSB_INFO("BASIC PRIMITIVES TEST: END (SUCCESS)\n");
359 }
360 return errval;
361 }
362
rsb_blas_limit_mul_tester(const rsb_coo_idx_t * aIA,const rsb_coo_idx_t * aJA,const void * aVA,const rsb_coo_idx_t * bIA,const rsb_coo_idx_t * bJA,const void * bVA,const rsb_coo_idx_t m,const rsb_coo_idx_t k,const rsb_coo_idx_t n,const rsb_nnz_idx_t annz,const rsb_nnz_idx_t bnnz,rsb_type_t typecode)363 static rsb_err_t rsb_blas_limit_mul_tester(
364 const rsb_coo_idx_t*aIA, const rsb_coo_idx_t*aJA, const void* aVA,
365 const rsb_coo_idx_t*bIA, const rsb_coo_idx_t*bJA, const void* bVA,
366 const rsb_coo_idx_t m, const rsb_coo_idx_t k, const rsb_coo_idx_t n,
367 const rsb_nnz_idx_t annz, const rsb_nnz_idx_t bnnz, rsb_type_t typecode)
368 {
369 /* FIXME: need a complete checking suite, here. */
370 blas_sparse_matrix A = RSB_BLAS_INVALID_MATRIX;
371 blas_sparse_matrix B = RSB_BLAS_INVALID_MATRIX;
372 //blas_sparse_matrix C = RSB_BLAS_INVALID_MATRIX;
373 rsb_err_t errval = RSB_ERR_NO_ERROR;
374 struct rsb_mtx_t * mtxAp=NULL;
375 struct rsb_mtx_t * mtxBp=NULL;
376 struct rsb_mtx_t * mtxCp=NULL;
377 rsb_trans_t trans = RSB_TRANSPOSITION_N;
378 if((A = rsb__BLAS_Xuscr_begin( m, k, typecode )) == RSB_BLAS_INVALID_MATRIX ) goto err;
379 if( rsb__BLAS_Xuscr_insert_entries( A, annz, aVA, aIA, aJA) == RSB_BLAS_INVALID_MATRIX ) goto err;
380 if( rsb__BLAS_Xuscr_end(A) == RSB_BLAS_INVALID_MATRIX ) goto err;
381 if((B= rsb__BLAS_Xuscr_begin( k, n, typecode )) == RSB_BLAS_INVALID_MATRIX ) goto err;
382 if( rsb__BLAS_Xuscr_insert_entries( B, bnnz, bVA, bIA, bJA) == RSB_BLAS_INVALID_MATRIX ) goto err;
383 if( rsb__BLAS_Xuscr_end(B) == RSB_BLAS_INVALID_MATRIX ) goto err;
384 #if 1
385 mtxAp = rsb__BLAS_inner_matrix_retrieve(A);
386 mtxBp = rsb__BLAS_inner_matrix_retrieve(B);
387 if(!mtxAp || !mtxBp)
388 {
389 RSB_ERROR(RSB_ERRM_EM);
390 goto err;// it's not ok.
391 }
392
393
394 #endif
395 /* TODO: need a complete checking suite, here. */
396 if((mtxCp = rsb__do_matrix_mul(typecode,RSB_TRANSPOSITION_N,NULL,mtxAp,trans,NULL,mtxBp,&errval))==NULL)
397 {
398 if(errval == RSB_ERR_LIMITS)
399 {
400 RSB_INFO("failed computing a dense %d x %d matrix (for internal memory limits reasons--it's ok)!\n",m,n);
401 errval = RSB_ERR_NO_ERROR;
402 }
403 else
404 {
405 RSB_INFO("failed computing a dense %d x %d matrix (unknown reasons--it's not ok)!\n",m,n);
406 }
407 }
408 else
409 {
410 if(!rsb__mtx_chk(mtxCp))
411 {
412 RSB_ERROR("matrix does not seem to be correctly built\n");
413 RSB_DO_ERROR_CUMULATE(errval,RSB_ERR_INTERNAL_ERROR);
414 }
415 RSB_MTX_FREE(mtxCp);
416 }
417
418 if(rsb__BLAS_Xusds( A ) == RSB_BLAS_ERROR) goto err;
419 if(rsb__BLAS_Xusds( B ) == RSB_BLAS_ERROR) goto err;
420 goto ret;
421 err:
422 errval = RSB_ERR_INTERNAL_ERROR;
423 ret:
424 return errval;
425 }
426
rsb_blas_limit_instancing_tester(const rsb_coo_idx_t * IA,const rsb_coo_idx_t * JA,const void * VA,const rsb_coo_idx_t m,const rsb_coo_idx_t k,const rsb_nnz_idx_t nnz,rsb_type_t typecode)427 static rsb_err_t rsb_blas_limit_instancing_tester(const rsb_coo_idx_t*IA, const rsb_coo_idx_t*JA, const void* VA, const rsb_coo_idx_t m, const rsb_coo_idx_t k, const rsb_nnz_idx_t nnz, rsb_type_t typecode)
428 {
429 /* FIXME: need a complete checking suite, here. */
430 rsb_err_t errval = RSB_ERR_NO_ERROR;
431 struct rsb_mtx_t * mtxAp=NULL;
432 #if 0
433 blas_sparse_matrix A = RSB_BLAS_INVALID_MATRIX;
434 if((A=BLAS_duscr_begin( m, k )) == RSB_BLAS_INVALID_MATRIX )
435 goto err;
436 // RSB_INFO("*\n");
437 if( rsb__BLAS_Xuscr_insert_entries( A, nnz, VA, IA, JA) == RSB_BLAS_INVALID_MATRIX )
438 goto err;
439 // RSB_INFO("*\n");
440 if( rsb__BLAS_Xuscr_end(A) == RSB_BLAS_INVALID_MATRIX )
441 goto err;
442
443 mtxAp = rsb__BLAS_inner_matrix_retrieve(A);
444 #else
445 mtxAp = rsb__do_mtx_alloc_from_coo_const(VA,IA,JA,nnz,typecode,m,k,RSB_DEFAULT_ROW_BLOCKING,RSB_DEFAULT_COL_BLOCKING,RSB_FLAG_NOFLAGS,&errval);
446 #endif
447 if(!mtxAp)
448 {
449 //RSB_ERROR(RSB_ERRM_EM);
450 //goto err;// it's ok.
451 if(RSB_SOME_ERROR(errval)) { RSB_ERROR("failed allocating a %d x %d matrix !\n",m,k); }
452 }
453 else
454 {
455 if(!rsb__mtx_chk(mtxAp))
456 {
457 RSB_ERROR("matrix does not seem to be correctly built\n");
458 RSB_DO_ERROR_CUMULATE(errval,RSB_ERR_INTERNAL_ERROR);
459 }
460 }
461 /* FIXME: need a complete checking suite, here. */
462 #if 0
463 if(rsb__BLAS_Xusds( A ) == RSB_BLAS_INVALID_MATRIX)
464 goto err;
465 #else
466 RSB_MTX_FREE(mtxAp);
467 if(errval == RSB_ERR_LIMITS)
468 {
469 RSB_INFO("failed instancing of (dense?) %d x %d matrix (it's ok)!\n",m,k);
470 errval = RSB_ERR_NO_ERROR;
471 }
472 else
473 {
474 errval = RSB_ERR_INTERNAL_ERROR;// goto err;
475 goto err;
476 }
477 #endif
478 // RSB_INFO("*\n");
479 RSB_INFO("instancing %d x %d, %d nnz succeeded\n",m,k,nnz);
480 goto ret;
481 err:
482 errval = RSB_ERR_INTERNAL_ERROR;
483 ret:
484 return errval;
485 }
486
rsb_blas_limit_cases_tester(void)487 rsb_err_t rsb_blas_limit_cases_tester(void)
488 {
489 /**
490 * \ingroup gr_internals
491 *
492 * FIXME: shall perform some serious (more iterations) test, here.
493 * FIXME: shall test on limits nonzeroes for various operations (sort, dups, etc)
494 * */
495 rsb_err_t errval = RSB_ERR_NO_ERROR;
496 #ifndef RSB_NUMERICAL_TYPE_DOUBLE
497 RSB_INFO("SKIPPING BASIC LIMIT CASES TEST (UNFINISHED TESTING SUITE)\n");
498 #else /* RSB_NUMERICAL_TYPE_DOUBLE */
499 rsb_type_t typecode = RSB_NUMERICAL_TYPE_DOUBLE;
500 const rsb_nnz_idx_t nnz=4;
501 const rsb_coo_idx_t dima[]={
502 4,RSB_MAX_SHORTIDX_MATRIX_DIM
503 //,RSB_MAX_MATRIX_DIM-1
504 ,2<<18 ,2<<20
505 //,RSB_MAX_MATRIX_DIM
506 //RSB_MAX_MATRIX_DIM+1000000+RSB_NNZ_BLK_MAX-3
507 ,2<<22 ,2<<24
508 //,2<<26 ,2<<28
509 //,RSB_MAX_MATRIX_DIM-1000000
510 , RSB_MAX_MATRIX_DIM
511 };
512 //const rsb_coo_idx_t dima[]={2};
513 rsb_int_t dimi;
514 RSB_INFO("BASIC LIMIT CASES TEST: BEGIN\n");
515 RSB_INFO("BASIC LIMIT CASES TEST: BEGIN\n");
516 RSB_INFO("(please do not worry if some tests fail due to insufficient memory)\n");
517 #if 1
518 RSB_INFO("(forcing allocations to be memory resident)\n");
519 rsb__lock_as_memory_resident(RSB_BOOL_TRUE); /* TODO: check return value here! */
520 #else
521 #endif
522
523 if(1)
524 for(dimi=0;dimi<sizeof(dima)/sizeof(dima[0]);++dimi)
525 {
526 const rsb_coo_idx_t dim=dima[dimi];
527 rsb_coo_idx_t IA[nnz];
528 rsb_coo_idx_t JA[nnz];
529 const double VA[]={11,22,33,44};
530 //const double X[dim],Y[dim]; const double alpha = 1.0;
531 //rsb__util_set_array_to_converted_integer(X,typecode,m,1,1);
532 //rsb__util_set_array_to_converted_integer(Y,typecode,k,1,0);
533 rsb_coo_idx_t m,k;
534 RSB_INFO("testing instantiation %d-sized, %d nnz\n",dim,nnz);
535 /* * * * * * * * * * * * * * * * * * * * * * * */
536 /* FIXME: need a `rotation' routine, here */
537 /* * * * * * * * * * * * * * * * * * * * * * * */
538 m=dim,k=dim;
539 IA[0]=0; IA[1]=0; IA[2]=dim-1; IA[3]=dim-1;
540 JA[0]=0; JA[1]=dim-1; JA[2]=0; JA[3]=dim-1;
541 errval = rsb_blas_limit_instancing_tester(IA, JA, VA, m, k, nnz, typecode);
542 if(RSB_SOME_ERROR(errval)) { }
543 /* * * * * * * * * * * * * * * * * * * * * * * */
544 m=1,k=dim;
545 IA[0]=0; IA[1]=0; IA[2]=0; IA[3]=0;
546 JA[0]=0; JA[1]=1; JA[2]=dim-2; JA[3]=dim-1;
547 errval = rsb_blas_limit_instancing_tester(IA, JA, VA, m, k, nnz, typecode);
548 if(RSB_SOME_ERROR(errval)) { }
549 /* * * * * * * * * * * * * * * * * * * * * * * */
550 m=dim,k=1;
551 IA[0]=0; IA[1]=1; IA[2]=dim-2; IA[3]=dim-1;
552 JA[0]=0; JA[1]=0; JA[2]=0; JA[3]=0;
553 errval = rsb_blas_limit_instancing_tester(IA, JA, VA, m, k, nnz, typecode);
554 if(RSB_SOME_ERROR(errval)) { }
555 /* * * * * * * * * * * * * * * * * * * * * * * */
556 m=dim,k=dim; IA[0]=0; JA[0]=0;
557 errval = rsb_blas_limit_instancing_tester(IA, JA, VA, m, k, 1, typecode);
558 if(RSB_SOME_ERROR(errval)) { }
559 /* * * * * * * * * * * * * * * * * * * * * * * */
560 m=dim,k=dim; IA[0]=dim-1; JA[0]=0;
561 errval = rsb_blas_limit_instancing_tester(IA, JA, VA, m, k, 1, typecode);
562 if(RSB_SOME_ERROR(errval)) { }
563 /* * * * * * * * * * * * * * * * * * * * * * * */
564 m=dim,k=dim; IA[0]=0; JA[0]=dim-1;
565 errval = rsb_blas_limit_instancing_tester(IA, JA, VA, m, k, 1, typecode);
566 if(RSB_SOME_ERROR(errval)) { }
567 /* * * * * * * * * * * * * * * * * * * * * * * */
568 m=dim,k=dim; IA[0]=dim-1; JA[0]=dim-1;
569 errval = rsb_blas_limit_instancing_tester(IA, JA, VA, m, k, 1, typecode);
570 if(RSB_SOME_ERROR(errval)) { }
571 /* * * * * * * * * * * * * * * * * * * * * * * */
572 }
573
574 if(1)
575 {
576 const rsb_nnz_idx_t dim = RSB_MAX_SHORTIDX_MATRIX_DIM+1;
577 //const rsb_nnz_idx_t dim=10;
578 rsb_coo_idx_t*aIA=NULL, *aJA=NULL, *bIA=NULL, *bJA=NULL;
579 void* aVA=NULL, * bVA=NULL;
580 const rsb_coo_idx_t m=dim, k=dim, n=dim;
581 const rsb_nnz_idx_t annz=dim+1, bnnz=dim+1;
582 rsb_type_t typecode = RSB_NUMERICAL_TYPE_DOUBLE;
583 /* size_t el_size = RSB_NUMERICAL_TYPE_SIZE(typecode); */
584 RSB_INFO("testing spmult for %d-sized, %d nnz\n",dim,nnz);
585 if(RSB_SOME_ERROR(errval = rsb__util_coo_alloc(&aVA,&aIA,&aJA,annz,typecode,RSB_BOOL_TRUE))){goto erra;}
586 if(RSB_SOME_ERROR(errval = rsb__util_coo_alloc(&bVA,&bIA,&bJA,bnnz,typecode,RSB_BOOL_TRUE))){goto erra;}
587 rsb__util_coo_array_set(aJA,annz,0);
588 rsb__util_coo_array_set_sequence(aIA,annz,0,1);
589 rsb__util_coo_array_set(bIA,bnnz,0);
590 rsb__util_coo_array_set_sequence(bJA,bnnz,0,1);
591 aIA[annz-1]=dim/2; aJA[bnnz-1]=dim/2;
592 bIA[annz-1]=dim/2; bJA[bnnz-1]=dim/2;
593 if(RSB_SOME_ERROR(rsb__fill_with_ones (aVA,typecode,dim,1))){ errval = RSB_ERR_INTERNAL_ERROR; goto erra; }
594 if(RSB_SOME_ERROR(rsb__fill_with_ones (bVA,typecode,dim,1))){ errval = RSB_ERR_INTERNAL_ERROR; goto erra; }
595 errval = rsb_blas_limit_mul_tester( aIA, aJA, aVA, bIA, bJA, bVA, m, k, n, annz, bnnz, typecode);
596 erra:
597 RSB_CONDITIONAL_FREE(aIA);
598 RSB_CONDITIONAL_FREE(aJA);
599 RSB_CONDITIONAL_FREE(aVA);
600 RSB_CONDITIONAL_FREE(bIA);
601 RSB_CONDITIONAL_FREE(bJA);
602 RSB_CONDITIONAL_FREE(bVA);
603 if(RSB_SOME_ERROR(errval)) {RSB_ERROR("!\n"); goto err;}
604 }
605
606 RSB_INFO("BASIC LIMIT CASES TEST: END\n");
607 goto ret;
608 err:
609 RSB_INFO("BASIC LIMIT CASES TEST: END : FAILURE\n");
610 ret:
611 #endif /* RSB_NUMERICAL_TYPE_DOUBLE */
612 return errval;
613 }
614
615
616 #if RSB_WANT_COO_BEGIN
rsb_mtx_alloc_from_coo_test(void)617 static rsb_err_t rsb_mtx_alloc_from_coo_test(void)
618 {
619 rsb_err_t errval = RSB_ERR_NO_ERROR;
620 const rsb_nnz_idx_t nnzA=4; /* matrix nonzeroes count */
621 const rsb_coo_idx_t nrA=3; /* matrix rows count */
622 const rsb_coo_idx_t ncA=3; /* matrix columns count */
623 rsb_coo_idx_t IA[]={0,1,2,2};
624 rsb_coo_idx_t JA[]={0,1,2,2};
625 RSB_DEFAULT_POSSIBLY_FIRST_BLAS_TYPE VA[]={11,22,32,1};/* values of nonzeroes */
626 rsb_type_t typecode = RSB_NUMERICAL_TYPE_FIRST_BLAS;
627 struct rsb_mtx_t * mtxAp = NULL;
628
629 if( RSB_NUMERICAL_TYPE_FIRST_BLAS == RSB_NUMERICAL_TYPE_INVALID_TYPE )
630 { RSB_INFO("SKIPPING A TEST (no BLAS types in)\n"); goto ret; }
631
632 mtxAp = rsb__do_mtx_alloc_from_coo_begin(nnzA,typecode,nrA,ncA,RSB_FLAG_NOFLAGS,&errval);
633 if(RSB_SOME_ERROR(errval))goto err;
634 if( mtxAp == NULL ){ errval = RSB_ERR_INTERNAL_ERROR;goto err; }
635 if(RSB_SOME_ERROR(errval = rsb__do_set_elements(mtxAp,VA,IA,JA,nnzA,RSB_FLAG_NOFLAGS)))goto err;
636 if(RSB_SOME_ERROR(errval = rsb__do_mtx_alloc_from_coo_end(&mtxAp)))goto err;
637 RSB_MTX_FREE(mtxAp);
638 goto ret;
639 err:
640 RSB_INFO("!\n");
641 ret:
642 return errval;
643 }
644 #endif /* RSB_WANT_COO_BEGIN */
645
646 #if RSB_ALLOW_INTERNAL_GETENVS
rsb__txt_ar(const char * c,int * ap,int * lp)647 static void rsb__txt_ar(const char*c, int* ap, int*lp)
648 {
649 int nul = 0,ci,l=0;
650
651 if(!c)
652 goto err;
653
654 do
655 {
656 while(*c!=nul && !isdigit(*c))++c;
657 ci = rsb__util_atoi(c);/* Flawfinder: ignore */
658
659 if(isdigit(*c))
660 ap[l++] = ci;
661 while(*c && isdigit(*c))++c;
662 }
663 while(*c);
664
665 RSB_ASSIGN_IF(lp,l)
666 err:
667 return;
668 }
669 #endif /* RSB_ALLOW_INTERNAL_GETENVS */
670
rsb_blas_bigger_matrices_tester(struct rsb_tester_options_t * top)671 rsb_err_t rsb_blas_bigger_matrices_tester(struct rsb_tester_options_t * top)
672 {
673 /**
674 * \ingroup gr_internals
675 * */
676 rsb_err_t errvalf = RSB_ERR_NO_ERROR;
677 rsb_err_t errval = RSB_ERR_NO_ERROR;
678 // full blas compliance:
679 enum blas_trans_type transTa[]={blas_no_trans,blas_trans,blas_conj_trans};
680 enum blas_symmetry_type stypea[]={
681 blas_lower_triangular
682 ,blas_upper_triangular
683 ,blas_lower_symmetric
684 ,blas_general
685 /* ,blas_upper_symmetric*/ /* one symmetry is enough for testing purposes ... */
686 ,blas_lower_hermitian
687 //,blas_upper_hermitian
688 };
689 rsb_blas_int_t incXa[]={1,2};
690 rsb_blas_int_t incBa[]={1,2};
691 rsb_blas_int_t alphaa[]={-2,-1,1,2};
692 #if (RSB_IMPLEMENTED_SOME_BLAS_TYPES>0)
693 rsb_type_t typecodea[]=RSB_MATRIX_SPBLAS_TYPE_CODES_ARRAY;
694 #else /* RSB_IMPLEMENTED_SOME_BLAS_TYPES */
695 rsb_type_t typecodea[]={RSB_NUMERICAL_TYPE_INVALID_TYPE};/* bogus definition */
696 #endif /* RSB_IMPLEMENTED_SOME_BLAS_TYPES */
697 enum blas_diag_type diaga[]={blas_non_unit_diag,blas_unit_diag};
698
699 // FIXME: should implement a routine to conjugate complex matrices before testing !
700
701 //enum blas_trans_type transTa[]={blas_no_trans};
702 //enum blas_trans_type transTa[]={blas_trans};
703 //enum blas_trans_type transTa[]={blas_conj_trans};
704 //enum blas_trans_type transTa[]={blas_no_trans,blas_trans};
705 //enum blas_symmetry_type stypea[]={blas_lower_triangular};
706 //enum blas_symmetry_type stypea[]={blas_upper_triangular};
707 //rsb_blas_int_t incXa[]={1};
708 //rsb_blas_int_t incBa[]={2};
709 //rsb_blas_int_t alphaa[]={1};
710 //rsb_blas_int_t alphaa[]={-1};
711 //rsb_blas_int_t incXa[]={2};
712 //rsb_blas_int_t alphaa[]={-1,1};
713 //rsb_blas_int_t alphaa[]={-2};
714 rsb_blas_int_t betaa[]={1,0};// FIXME: there the Sparse BLAS interface works only with beta=1
715 //rsb_type_t typecodea[]={RSB_NUMERICAL_TYPE_DOUBLE};
716 //rsb_type_t typecodea[]={RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX};
717 //rsb_type_t typecodea[]={RSB_NUMERICAL_TYPE_FLOAT};
718 //rsb_type_t typecodea[]=RSB_MATRIX_TYPE_CODES_ARRAY;
719
720 //const rsb_blas_int_t dima_lcl=1;
721 const rsb_blas_int_t dima_lcl=2;
722 const rsb_blas_int_t dima_pcl=dima_lcl;
723 const rsb_blas_int_t dima_tcl=dima_lcl;
724 //enum blas_diag_type diaga[]={blas_unit_diag};
725 //enum blas_diag_type diaga[]={blas_non_unit_diag};
726 //enum blas_diag_type diaga[]={blas_unit_diag};
727 const rsb_blas_int_t dimas=dima_pcl+dima_tcl*RSB_MAX_SUPPORTED_CACHE_LEVELS+dima_lcl;
728 rsb_blas_int_t dima[dimas];
729 rsb_blas_int_t dims=0;
730 rsb_blas_int_t diagi=RSB_INVALID_BLAS_INT_IDX_VAL;
731 rsb_blas_int_t transTi=RSB_INVALID_BLAS_INT_IDX_VAL;
732 rsb_blas_int_t alphai=RSB_INVALID_BLAS_INT_IDX_VAL;
733 rsb_blas_int_t betai=RSB_INVALID_BLAS_INT_IDX_VAL;
734 rsb_blas_int_t stypei=RSB_INVALID_BLAS_INT_IDX_VAL;
735 rsb_blas_int_t incXi=RSB_INVALID_BLAS_INT_IDX_VAL;
736 rsb_blas_int_t dimi=RSB_INVALID_BLAS_INT_IDX_VAL;
737 rsb_blas_int_t typecodei=RSB_INVALID_BLAS_INT_IDX_VAL;
738 rsb_blas_int_t incBi=RSB_INVALID_BLAS_INT_IDX_VAL;
739 rsb_blas_int_t cl=RSB_INVALID_BLAS_INT_IDX_VAL,cln = rsb__get_cache_levels_num();
740 rsb_blas_int_t passed=0,failed=0;
741 rsb_blas_int_t instantiated_some_recursive=0;
742 #if RSB_TESTER_ALLOW_TIMEOUT
743 rsb_time_t tt = RSB_TIME_ZERO,tt0=rsb_time();
744 struct rsb_tester_options_t to;
745 #endif /* RSB_TESTER_ALLOW_TIMEOUT */
746 rsb_blas_int_t isempty=0,isinvertible=1;
747 /* FIXME: in the future, may use these indices (isemptym) to fill the matrix with a particular value */
748 #if RSB_ALLOW_EMPTY_MATRICES
749 #if RSB_ALLOW_ZERO_DIM
750 rsb_blas_int_t isemptym=3;
751 #else
752 rsb_blas_int_t isemptym=2;
753 #endif
754 #else /* RSB_ALLOW_EMPTY_MATRICES */
755 rsb_blas_int_t isemptym=1;
756 #endif /* RSB_ALLOW_EMPTY_MATRICES */
757 const rsb_char_t*btps=RSB_EMPTY_STRING;
758 rsb_int_t iat=1;
759 #if RSB_ALLOW_INTERNAL_GETENVS
760 rsb_blas_int_t maxtc = 0; /* max tests count, current tests count */
761 #endif /* RSB_ALLOW_INTERNAL_GETENVS */
762
763 if( (sizeof(typecodea)==0)
764 #if (RSB_IMPLEMENTED_SOME_BLAS_TYPES==0)
765 || 1
766 #endif /* RSB_IMPLEMENTED_SOME_BLAS_TYPES */
767 )
768 {
769 // FIXME: new
770 RSB_INFO("Did not configure any BLAS-standard type: thus skipping BLAS-based testing.\n");
771 goto ret;
772 }
773 #if RSB_TESTER_ALLOW_TIMEOUT
774 if(!top)
775 errval = rsb_blas_tester_options_init(&to);
776 else
777 rsb__memcpy(&to,top,sizeof(to));
778 #endif /* RSB_TESTER_ALLOW_TIMEOUT */
779 if(RSB_SOME_ERROR(errval))
780 {
781 RSB_ERROR("!\n");
782 goto ret;
783 }
784 errval = rsb_basic_primitives_tester();
785 if(RSB_SOME_ERROR(errval))
786 {
787 RSB_ERROR("!\n");
788 goto ret;
789 }
790
791 #if RSB_WANT_COO_BEGIN
792 /* FIXME: this mini-test is insufficient! */
793 if(RSB_SOME_ERROR(errval=rsb_mtx_alloc_from_coo_test()))
794 { RSB_ERROR("!\n"); goto ret; }
795 #endif /* RSB_WANT_COO_BEGIN */
796
797
798 #if RSB_HAVE_ISATTY
799 #if RSB_HAVE_STREAMS
800 if( rsb_global_session_handle.out_stream )
801 iat=( isatty(rsb__fileno(rsb_global_session_handle.out_stream)) );
802 else
803 iat=0;
804 #endif /* RSB_HAVE_STREAMS */
805 #endif /* RSB_HAVE_ISATTY */
806
807 if(to.wcs==RSB_BOOL_TRUE)
808 btps=RSB_CLEARTERM_STRING;
809 if((to.wqc==RSB_BOOL_TRUE) && (!iat))
810 to.wqt=RSB_BOOL_TRUE;
811 RSB_INFO("ADVANCED SPARSE BLAS TEST: BEGIN\n");
812 #if 1
813 for(cl=0;cl<dima_pcl;++cl)
814 {
815 // 1,2,..
816 dima[dims++]=1<<cl;
817 }
818 #else
819 // if(dims<dimas)dima[dims++]=39;
820 // if(dims<dimas)dima[dims++]=724;
821 // if(dims<dimas)dima[dims++]=362;
822 // if(dims<dimas)dima[dims++]=1;
823 if(dims<dimas)dima[dims++]=2;
824 // if(dims<dimas)dima[dims++]=3;
825 #endif
826 #if 1
827 for(cl=1;cl<=cln && dims<dimas;++cl)
828 {
829 // around cache size
830 long cs = rsb__get_lnc_size(cl);
831 rsb_blas_int_t i;
832 for(i=1;i<=dima_tcl;++i)
833 dima[dims++]=((1<<i)*2*sqrt(cs))/(4*sizeof(rsb_coo_idx_t));
834 }
835 if((cl=cln)>0)
836 {
837 // more than outermost cache size
838 rsb_blas_int_t i;
839 long cs = rsb__get_lnc_size(cl);
840 for(i=1;i<=dima_lcl && dims<dimas;++i)
841 {
842 long ndim=(((i)*(1<<dima_tcl)*2*sqrt(cs))/(4*sizeof(rsb_coo_idx_t)));
843 if(ndim > dima[dims]) // to avoid duplicates
844 dima[dims++]=ndim;
845 }
846 }
847 #endif
848
849 #if RSB_ALLOW_INTERNAL_GETENVS
850 rsb__txt_ar(getenv("RSB_BMT_ALPHA"), & alphaa[0], NULL);
851 rsb__txt_ar(getenv("RSB_BMT_INCXA"), & incXa[0], NULL);
852 rsb__txt_ar(getenv("RSB_BMT_INCBA"), & incBa[0], NULL);
853 rsb__txt_ar(getenv("RSB_BMT_SYMMA"), & stypea[0], NULL);
854 rsb__txt_ar(getenv("RSB_BMT_STYPA"), &typecodea[0], NULL);
855 rsb__txt_ar(getenv("RSB_BMT_DIAGA"), & diaga[0], NULL);
856 rsb__txt_ar(getenv("RSB_BMT_TRANSA"), & transTa[0], NULL);
857 rsb__txt_ar(getenv("RSB_BMT_DIMA"), & dima[0],&dims);
858 #if RSB_ALLOW_EMPTY_MATRICES
859 if(getenv("RSB_BMT_ISEMPTYM")) isemptym = rsb__util_atoi(getenv("RSB_BMT_ISEMPTYM"));
860 #endif /* RSB_ALLOW_EMPTY_MATRICES */
861 #endif /* RSB_ALLOW_INTERNAL_GETENVS */
862
863 #if 1
864 if(
865 (rsb__do_psblas_trans_to_rsb_trans(RSB_PSBLAS_TRANS_N) != RSB_TRANSPOSITION_N) ||
866 (rsb__do_psblas_trans_to_rsb_trans(RSB_PSBLAS_TRANS_T) != RSB_TRANSPOSITION_T) ||
867 (rsb__do_psblas_trans_to_rsb_trans(RSB_PSBLAS_TRANS_C) != RSB_TRANSPOSITION_C)
868 )
869 {RSB_ERROR("!\n"); goto ret;}
870 #endif
871
872 //dims=0;
873 //dima[dims++]=45;
874 //dima[dims++]=362;
875 //dima[dims++]=500;
876 //dima[dims++]=499;
877 //dima[dims++]=724;
878 //dima[dims++]=1448;
879 // typecodei=3;
880 for(dimi=0;dimi<dims;++dimi)
881 for(stypei=0;stypei<sizeof(stypea)/sizeof(enum blas_symmetry_type);++stypei)
882 for(typecodei=0;typecodei<sizeof(typecodea)/sizeof(rsb_type_t);++typecodei)
883 for(incXi=0;incXi<sizeof(incXa)/sizeof(rsb_blas_int_t);++incXi)
884 for(incBi=0;incBi<sizeof(incBa)/sizeof(rsb_blas_int_t);++incBi)
885 for(alphai=0;alphai<sizeof(alphaa)/sizeof(rsb_blas_int_t);++alphai)
886 for(betai=0;betai<sizeof(betaa)/sizeof(rsb_blas_int_t);++betai)
887 for(diagi=0;diagi<sizeof(diaga)/sizeof(enum blas_diag_type);++diagi)
888 for(transTi=0;transTi<sizeof(transTa)/sizeof(enum blas_trans_type);++transTi)
889 #if RSB_ALLOW_EMPTY_MATRICES
890 for(isempty=0;isempty<isemptym;++isempty)
891 #endif /* RSB_ALLOW_EMPTY_MATRICES */
892 {
893 const rsb_blas_int_t is_really_empty=isempty && (diaga[diagi]!=blas_unit_diag);
894 blas_sparse_matrix T = RSB_BLAS_INVALID_MATRIX;
895 void *B=NULL,*X=NULL,*D=NULL;
896 rsb_blas_int_t dim=dima[dimi];
897 rsb_coo_idx_t *IA=NULL,*JA=NULL;
898 void * VA=NULL;
899 rsb_blas_int_t nnz = RSB_BLAS_INT_MAX;
900 rsb_type_t typecode=typecodea[typecodei];
901 enum blas_trans_type transT=transTa[transTi];
902 rsb_blas_int_t incX=incXa[incXi],incB=incBa[incBi],incD=1;
903 size_t el_size = RSB_NUMERICAL_TYPE_SIZE(typecode);
904 enum blas_symmetry_type stype=stypea[stypei];
905 rsb_aligned_t alpha_inv[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
906 rsb_aligned_t inrm[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
907 rsb_aligned_t inrm_[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
908 rsb_aligned_t alpha[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
909 rsb_aligned_t beta[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
910 struct rsb_mtx_t * mtxAp=NULL;
911 struct rsb_mtx_t *cmatrix=NULL;
912 struct rsb_mtx_t *kmatrix=NULL;
913 size_t extra_vels=0;
914 rsb_nnz_idx_t rnnz=0,ndnnz,rnz=0;
915 rsb_submatrix_idx_t submatrices=0;
916 rsb_trans_t trans = rsb__blas_trans_to_rsb_trans(transT);
917 rsb_char_t tc = RSB_TRANSPOSITION_AS_CHAR(trans);
918 rsb_blas_int_t mmi,msmd=100;
919 rsb_coo_idx_t coov;
920 rsb_nnz_idx_t nnzv;
921 rsb_aligned_t zero[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
922 rsb_aligned_t one[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
923 rsb_aligned_t two[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
924 rsb_aligned_t three[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
925 rsb_int do_tune_test = 0;
926 #if RSB_ALLOW_ZERO_DIM
927 if(isempty>=2 && dim > 1) continue;
928 if(isempty>=2) dim=0;
929 #endif
930 rsb__util_set_area_to_converted_integer(one,typecode,1);
931 rsb__util_set_area_to_converted_integer(two,typecode,2);
932 rsb__util_set_area_to_converted_integer(three,typecode,3);
933
934 rsb__util_set_area_to_fraction_of_integer(alpha_inv,alphaa[alphai],typecode);
935
936 // ... need asserts ...
937 rsb__util_set_area_to_converted_integer(alpha,typecode,alphaa[alphai]);
938 rsb__util_set_area_to_converted_integer(beta,typecode,betaa[betai]);
939 rsb__util_set_area_to_converted_integer(zero,typecode,0);
940 #if RSB_ALLOW_ZERO_DIM
941 if(isempty>=2) extra_vels=1;
942 #endif
943 X = rsb__calloc(el_size*(dim*incX+extra_vels));
944 B = rsb__calloc(el_size*(dim*incB+extra_vels));
945 D = rsb__calloc(el_size*(dim*incD+extra_vels));
946 if(!X || !B || !D)
947 {
948 RSB_ERROR("failed allocating a vector!\n"); goto err;
949 }
950
951 /* generate a triangular matrix */
952 /* FIXME: should make sure that the matrix is recursive, somehow. */
953 errval = rsb__generate_dense_lower_triangular_coo(dim,1,&IA,&JA,&VA,&nnz,typecode);
954 if(RSB_SOME_ERROR(errval))
955 {
956 RSB_ERROR("!\n"); goto err;
957 }
958 #if RSB_ALLOW_EMPTY_MATRICES
959 if(isempty)
960 {
961 RSB_DO_ERROR_CUMULATE(errval,rsb__cblas_Xscal(typecode,nnz,&zero,VA,1));
962 if(RSB_SOME_ERROR(errval))
963 {RSB_ERROR("!\n"); goto err;}
964 }
965 #endif /* RSB_ALLOW_EMPTY_MATRICES */
966 isinvertible=(diaga[diagi]==blas_unit_diag||!isempty);
967 isinvertible&=(stype != blas_general);
968 /* isinvertible&=(stype != blas_lower_symmetric);
969 isinvertible&=(stype != blas_upper_symmetric);
970 */
971 isinvertible&=(stype != blas_upper_hermitian);
972 isinvertible&=(stype != blas_lower_hermitian);
973 ndnnz=nnz-(diaga[diagi]==blas_unit_diag?dim:0);
974 if(ndnnz>nnz)
975 {
976 RSB_ERROR("!\n"); goto err;
977 }
978 if(nnz > 0 && (!VA || !IA || !JA))
979 {
980 RSB_ERROR("!\n");
981 goto err;
982 }
983
984 if(stype==blas_upper_triangular)
985 RSB_SWAP(rsb_coo_idx_t*,IA,JA);
986
987 #if 1 /* 20110425 */
988 if(incX==1)
989 if(incB==1)/* FIXME: shall propagate incX to the test routine, someday */
990 if(nnz>0)/* empty matrices are not supported for now */
991 if(!(diaga[diagi]==blas_unit_diag))/* FIXME: the accuracy test needs cleaned up input (i.e.: won't remove the diagonal) */
992
993 if(!RSB_IS_MATRIX_TYPE_COMPLEX(typecode))/* FIXME: shall fix many vector-operating routines, first */
994 {
995 /* FIXME: to be complete, shall implement symmetry/lower/upper/diagonal flags */
996 rsb_flags_t aflags = RSB_FLAG_NOFLAGS;
997 struct rsb_coo_matrix_t coo;
998 if(stype==blas_lower_symmetric) RSB_DO_FLAG_ADD(aflags,RSB_FLAG_SYMMETRIC);
999 if(stype==blas_upper_symmetric) RSB_DO_FLAG_ADD(aflags,RSB_FLAG_SYMMETRIC);
1000 if(stype==blas_upper_hermitian) RSB_DO_FLAG_ADD(aflags,RSB_FLAG_UPPER_HERMITIAN);
1001 if(stype==blas_lower_hermitian) RSB_DO_FLAG_ADD(aflags,RSB_FLAG_LOWER_HERMITIAN);
1002 if(diaga[diagi]==blas_unit_diag)RSB_DO_FLAG_ADD(aflags,RSB_FLAG_UNIT_DIAG_IMPLICIT);
1003 rsb__fill_coo_struct(&coo,VA,IA,JA,dim,dim,nnz,typecode);
1004 RSB_DO_ERROR_CUMULATE(errval,rsb__do_spmv_accuracy_test(&coo,NULL,0,aflags));
1005 if(RSB_SOME_ERROR(errval))
1006 {
1007 RSB_ERROR("!\n");
1008 RSB_DO_ERROR_CUMULATE(errval,rsb__do_spmv_accuracy_test(&coo,NULL,0,aflags));
1009 goto err;
1010 }
1011 }
1012 #endif
1013 T = rsb__BLAS_Xuscr_begin(dim,dim,typecode);
1014 if( T == RSB_BLAS_INVALID_MATRIX )
1015 {errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("error while calling uscr_begin\n"); goto err;}
1016 if( BLAS_ussp(T,stype) != RSB_BLAS_NO_ERROR )
1017 {errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("error while calling ussp(%d)\n",stype); goto err;}
1018 if( BLAS_ussp(T,diaga[diagi]) != RSB_BLAS_NO_ERROR )
1019 {errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("error while calling ussp(%d)\n",diaga[diagi]); goto err;}
1020 if( rsb__BLAS_Xuscr_insert_entries(T,nnz,VA,IA,JA) != RSB_BLAS_NO_ERROR)
1021 {errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("error while calling cr_insert_entries\n"); goto err;}
1022 if( rsb__BLAS_Xuscr_end(T) != RSB_BLAS_NO_ERROR )
1023 {errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("error finalizing matrix!\n"); goto err;}
1024 mtxAp = rsb__BLAS_inner_matrix_retrieve(T);
1025 if(!mtxAp)
1026 {
1027 RSB_ERROR(RSB_ERRM_NL);
1028 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1029 }
1030
1031 if( mtxAp->nnz>0 && rsb__get_index_storage_amount(mtxAp)==0 )
1032 {
1033 RSB_ERROR(RSB_ERRM_NL);
1034 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1035 }
1036
1037 if(diaga[diagi]==blas_non_unit_diag) /* FIXME */
1038 if(nnz > 0)
1039 {
1040 if( rsb__BLAS_Xusset_elements(T, IA, JA, VA, mtxAp->nnz) != RSB_BLAS_NO_ERROR )
1041 {
1042 RSB_ERROR(RSB_ERRM_NL);
1043 errval = RSB_ERR_INTERNAL_ERROR;
1044 goto err;
1045 }
1046 }
1047
1048 rnz=mtxAp->nnz;
1049 if(!rsb__mtx_chk(mtxAp))
1050 {
1051 RSB_ERROR(RSB_ERRM_NL);
1052 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1053 }
1054
1055 coov=0;
1056 rsb__do_mtx_get_info(mtxAp, RSB_MIF_MATRIX_COLS__TO__RSB_COO_INDEX_T ,&coov);
1057 if(coov!=dim){errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL);goto err;}
1058 coov=0;
1059 rsb__do_get_matrix_info_from_string(mtxAp,"RSB_MIF_MATRIX_COLS__TO__RSB_COO_INDEX_T",&coov,0);
1060 if(coov!=dim){errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL);goto err;}
1061
1062 coov=0;
1063 rsb__do_mtx_get_info(mtxAp, RSB_MIF_MATRIX_ROWS__TO__RSB_COO_INDEX_T ,&coov);
1064 if(coov!=dim){errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL);goto err;}
1065 coov=0;
1066 rsb__do_get_matrix_info_from_string(mtxAp,"RSB_MIF_MATRIX_ROWS__TO__RSB_COO_INDEX_T",&coov,0);
1067 if(coov!=dim){errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL);goto err;}
1068
1069 nnzv=0;
1070 rsb__do_mtx_get_info (mtxAp, RSB_MIF_MATRIX_NNZ__TO__RSB_NNZ_INDEX_T ,&nnzv);
1071 if(nnzv!=mtxAp->nnz){errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL);goto err;}
1072 nnzv=0;
1073 rsb__do_get_matrix_info_from_string(mtxAp,"RSB_MIF_MATRIX_NNZ__TO__RSB_NNZ_INDEX_T",&nnzv,0);
1074 if(nnzv!=mtxAp->nnz){errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL);goto err;}
1075
1076 {
1077 /* TODO: need a systematic tester */
1078 const int errstrlen=128;
1079 char errstr[errstrlen];
1080 strcpy(errstr,"");
1081 rsb__do_strerror_r(RSB_ERR_BADARGS,errstr,errstrlen);
1082 if(strlen(errstr)<1)
1083 RSB_LSTERR(RSB_ERRM_ES);
1084
1085 }
1086
1087 /* TODO: extract && test other values as well... */
1088
1089 // submatrices = rsb__submatrices(mtxAp);
1090 submatrices = rsb__terminal_recursive_matrix_count(mtxAp);
1091 instantiated_some_recursive+=(rsb__submatrices(mtxAp)>1?1:0);
1092 if(rsb__get_sizeof(mtxAp)<(
1093 (mtxAp->el_size*mtxAp->nnz)+ (sizeof(rsb_half_idx_t)*mtxAp->nnz)+0))
1094 {
1095 errval = RSB_ERR_INTERNAL_ERROR;
1096 RSB_ERROR("!\n");
1097 goto err;
1098 }
1099
1100 #if 1
1101 #if RSB_ALLOW_INTERNAL_GETENVS
1102 if(! ( getenv("RSB_BMT_SPMV") && rsb__util_atoi(getenv("RSB_BMT_SPMV")) == 0 ) )
1103 #endif /* RSB_ALLOW_INTERNAL_GETENVS */
1104 if(RSB_DO_FLAG_HAS(mtxAp->flags,RSB_FLAG_SYMMETRIC) || RSB_DO_FLAG_HAS(mtxAp->flags,RSB_FLAG_HERMITIAN))
1105 {
1106 /* FIXME: this gets NOT covered, it seems */
1107 for(mmi=0;mmi< (dim<msmd?3:2) ;++mmi)
1108 if(! (mmi==1 && ((incX!= 1) || (incB!=1) ) ))
1109 {
1110 const int nrhs=1;/* TODO: need more ... */
1111 /* TODO : should fill X and B with sentinel values ! */
1112 if(RSB_SOME_ERROR(rsb__cblas_Xscal(typecode,dim,&zero,X,incX)) || RSB_SOME_ERROR(rsb__fill_with_ones (B,typecode,dim,incB)))
1113 { errval = RSB_ERR_INTERNAL_ERROR; goto err; }
1114
1115
1116 if(mmi==0)
1117 if( rsb__BLAS_Xusmv(transT,alpha,T,B,incB,beta,X,incX) != RSB_BLAS_NO_ERROR )
1118 {
1119 errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("Symmetric USMV failed!\n"); goto err;
1120 }
1121
1122 if(mmi==1)
1123 if(RSB_SOME_ERROR( rsb__do_spmm(trans,alpha,mtxAp,nrhs,RSB_FLAG_WANT_COLUMN_MAJOR_ORDER,B,dim,beta,X,dim,RSB_OP_FLAG_DEFAULT)) )
1124 {
1125 errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("Symmetric USMV failed!\n"); goto err;
1126 }
1127
1128 if(mmi==2)
1129 if(RSB_SOME_ERROR( rsb__do_spmv_general(trans,alpha,mtxAp,B,incB,beta,X,incX,RSB_OP_FLAG_WANT_SERIAL RSB_DEFAULT_OUTER_NRHS_SPMV_ARGS)))
1130 {
1131 errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("Symmetric USMV failed!\n"); goto err;
1132 }
1133
1134 /* if(!isinvertible) */
1135 if(is_really_empty)
1136 rsb__util_set_array_to_converted_integer(B,typecode,dim,incB,0 );
1137 else
1138 {
1139 if(isempty)
1140 rsb__util_set_array_to_converted_integer(B,typecode,dim,incB, alphaa[alphai]);
1141 else
1142 rsb__util_set_array_to_converted_integer(B,typecode,dim,incB,dim*alphaa[alphai]);
1143 }
1144 if( RSB_SOME_ERROR(rsb__do_are_same(B,X,dim,typecode,incB,incX)) )
1145 {
1146 errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("Symmetric USMV computed wrong results!\n"); goto err;
1147 }
1148 goto err; /* ok. skip the remaining tests FIXME */
1149 }
1150 }
1151 #endif
1152 #if 1
1153 #if RSB_ALLOW_INTERNAL_GETENVS
1154 if(! ( getenv("RSB_BMT_SPGEMM") && rsb__util_atoi(getenv("RSB_BMT_SPGEMM")) == 0 ) )
1155 #endif /* RSB_ALLOW_INTERNAL_GETENVS */
1156 // FIXME: this is quite slow : SPGEMM is O(dim^3, and so shall be limited down to a certain threshold)
1157 if(dimi <= RSB_WANT_SPGEMM_TESTING_FOR_ONLY_FIRST_DIMI)// FIXME: spgemm cost increases quadratically..
1158 if(incX==1 && incB==1)
1159 if(stype==blas_lower_triangular || stype==blas_upper_triangular)
1160 if(!RSB_DO_TOOFEWNNZFORCSR(nnz,dim) /*&& typecode == RSB_NUMERICAL_TYPE_DOUBLE*/ )
1161 if(diaga[diagi]==blas_non_unit_diag)
1162 if( trans == RSB_TRANSPOSITION_N ) /* FIXME: this is just a workaround, 20140324 */
1163 {
1164 rsb_nnz_idx_t ldd=2*dim;
1165 rsb_bool_t rowmajor=(stype==blas_lower_triangular)?RSB_BOOL_TRUE:RSB_BOOL_FALSE;
1166 //rsb_nnz_idx_t ldd=1*dim;
1167 void *dVA = rsb__calloc(el_size*dim*ldd);
1168 rsb_aligned_t sum1[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
1169 rsb_aligned_t sum2[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
1170 rsb_aligned_t two[RSB_CONST_ENOUGH_ALIGNED_FOR_ANY_TYPE];
1171
1172 if(dVA)
1173 {
1174 /* rsb_trans_t ttrans = rsb__do_transpose_transposition(trans); */
1175 rsb_trans_t ttrans = (trans); /* FIXME: this is just a workaround, 20140324 */
1176
1177 RSB_DO_ERROR_CUMULATE(errval,rsb__do_spgemm_to_dense(typecode,trans,alpha,mtxAp,ttrans,beta,mtxAp,ldd,dim,dim,rowmajor,dVA,NULL,NULL));
1178 if(RSB_SOME_ERROR(errval)) { RSB_ERROR("!\n"); goto err; }
1179 RSB_DO_ERROR_CUMULATE(errval,rsb__util_vector_sum(sum1,dVA,typecode,ldd*dim));
1180 if(RSB_SOME_ERROR(errval)) { RSB_ERROR("!\n"); goto err; }
1181 // now: c <- alpha a ^ trans * beta a ^ trans
1182 RSB_DO_ERROR_CUMULATE(errval,rsb__do_spgemm_to_dense(typecode,trans,alpha,mtxAp,ttrans,beta,mtxAp,ldd,dim,dim,rowmajor,dVA,NULL,NULL));
1183 if(RSB_SOME_ERROR(errval)) { RSB_ERROR("!\n"); goto err; }
1184 RSB_DO_ERROR_CUMULATE(errval,rsb__util_vector_sum(sum2,dVA,typecode,ldd*dim));
1185 if(RSB_SOME_ERROR(errval)) { RSB_ERROR("!\n"); goto err; }
1186 // now: c <- 2 ( alpha a ^ trans * beta a ^ trans )
1187 if(RSB_SOME_ERROR(errval))
1188 {
1189 RSB_CONDITIONAL_FREE(dVA);
1190 RSB_ERROR(RSB_ERRM_FMMTDT);
1191 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1192 }
1193 rsb__util_set_area_to_converted_integer(two,typecode,2);
1194 rsb__util_vector_div(sum2,two,typecode,1);
1195 // TODO: there is risk of overflow, though..
1196 if( RSB_SOME_ERROR(rsb__do_are_same(sum1,sum2,1,typecode,1,1) ))
1197 {
1198 RSB_CONDITIONAL_FREE(dVA);
1199 RSB_ERROR(RSB_ERRM_FMMTDT);
1200 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1201 }
1202 // since a is lower triangular full with ones,
1203 // c is full, with values (given s=2(alpha+beta))
1204 // s s s s ...
1205 // s 2s 2s 2s ...
1206 // s 2s 3s 3s ...
1207 // s 2s 3s 4s ...
1208 // ...
1209 // (transposed, in the case trans is)
1210 // now: c <- 2 ( alpha a ^ trans * beta ^ trans ) + alpha a ^ trans
1211
1212 RSB_DO_ERROR_CUMULATE(errval,rsb__cblas_Xscal(typecode,ldd*dim,&zero,dVA,1));
1213 RSB_DO_ERROR_CUMULATE(errval,rsb__do_matrix_add_to_dense(alpha,mtxAp,ldd,dim,dim,rowmajor,dVA));
1214 RSB_DO_ERROR_CUMULATE(errval,rsb__util_vector_sum(sum1,dVA,typecode,ldd*dim));
1215 RSB_DO_ERROR_CUMULATE(errval,rsb__do_matrix_add_to_dense(alpha,mtxAp,ldd,dim,dim,rowmajor,dVA));
1216 RSB_DO_ERROR_CUMULATE(errval,rsb__util_vector_sum(sum2,dVA,typecode,ldd*dim));
1217 rsb__util_vector_div(sum2,two,typecode,1);
1218 // TODO: there is risk of overflow, though..
1219 if( RSB_SOME_ERROR(rsb__do_are_same(sum1,sum2,1,typecode,1,1) ))
1220 {
1221 RSB_CONDITIONAL_FREE(dVA);
1222 RSB_ERROR(RSB_ERRM_FMATDBC);
1223 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1224 }
1225 RSB_CONDITIONAL_FREE(dVA);
1226 if(RSB_SOME_ERROR(errval))
1227 {
1228 RSB_ERROR(RSB_ERRM_FMATD);
1229 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1230 }
1231 }
1232 //if((cmatrix = rsb__do_matrix_mul(RSB_TRANSPOSITION_N,NULL,mtxAp,trans,NULL,mtxAp,&errval))!=RSB_ERR_NO_ERROR)
1233 if((cmatrix = rsb__do_matrix_mul(typecode,RSB_TRANSPOSITION_N,NULL,mtxAp,trans,NULL,mtxAp,&errval))!=NULL)
1234 {
1235 // FIXME: ignoring scaling values!
1236 RSB_MTX_FREE(cmatrix);
1237 }
1238 else
1239 {
1240 RSB_ERROR(RSB_ERRM_FMM);
1241 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1242 }
1243 }
1244 #endif
1245 #if 1
1246 #if RSB_ALLOW_INTERNAL_GETENVS
1247 if(! ( getenv("RSB_BMT_SUM") && rsb__util_atoi(getenv("RSB_BMT_SUM")) == 0 ) )
1248 #endif /* RSB_ALLOW_INTERNAL_GETENVS */
1249 {
1250 if(diaga[diagi]!=blas_unit_diag && nnz>0)
1251 {
1252 cmatrix = NULL;
1253 errval = rsb__mtx_clone(&cmatrix,RSB_NUMERICAL_TYPE_SAME_TYPE,RSB_TRANSPOSITION_N,NULL,mtxAp,RSB_FLAG_IDENTICAL_FLAGS);
1254 // cmatrix = rsb__mtx_clone_simple(mtxAp);
1255 if( (cmatrix == NULL) || (!rsb__mtx_chk(cmatrix) ) )
1256 {
1257 if(cmatrix==NULL)
1258 {
1259 RSB_ERROR(RSB_ERRM_FMC);
1260 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1261 }
1262 else
1263 {
1264 RSB_ERROR(RSB_ERRM_CMINBC);
1265 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1266 }
1267 }
1268 else
1269 {
1270 /* very, very slow sparse matrices sum testing */
1271 struct rsb_mtx_t *smatrix=NULL;
1272 struct rsb_coo_matrix_t coo,coc;
1273 // s = m^trans + 2 * m^trans = 3 m^trans
1274 //smatrix = rsb__do_matrix_sum(typecode,mtxAp,one,trans,cmatrix,two,trans,&errval);
1275 RSB_BZERO_P(&coo);
1276 RSB_BZERO_P(&coc);
1277
1278 smatrix = rsb__do_matrix_sum(typecode,RSB_TRANSPOSITION_N,one,mtxAp,RSB_TRANSPOSITION_N,two,cmatrix,&errval);
1279 if(!smatrix ) {RSB_ERROR(RSB_ERRM_FCMS);errval = RSB_ERR_INTERNAL_ERROR;}
1280 if( RSB_SOME_ERROR(errval)) { RSB_ERROR(RSB_ERRM_FCMS); goto smerr; }
1281 #if 0
1282 if( smatrix->nnz > 1)
1283 {
1284 rsb_nnz_idx_t i;
1285 for(i=0;i<cmatrix->nnz;++i)RSB_STDOUT("%d \n",((rsb_half_idx_t*)(mtxAp->bindx))[i]);
1286 RSB_INFO("cmatrix:\n");
1287 rsb_print_matrix_t(cmatrix);
1288 RSB_INFO(" mtxAp:\n");
1289 rsb_print_matrix_t( mtxAp);
1290 RSB_STDOUT_MATRIX_SUMMARY(mtxAp), RSB_INFO("\n");
1291 RSB_INFO("smatrix:\n");
1292 rsb_print_matrix_t(smatrix);
1293 RSB_INFO("\n");
1294 }
1295 #endif
1296 if( smatrix->nnz != mtxAp->nnz)
1297 #if RSB_ALLOW_EMPTY_MATRICES
1298 if( !isempty )
1299 #endif /* RSB_ALLOW_EMPTY_MATRICES */
1300 {
1301 RSB_ERROR("seems like matrix sum failed (same pattern, no cancellation possible): %d + %d to %d nnz)\n",mtxAp->nnz,cmatrix->nnz,smatrix->nnz);
1302 errval = RSB_ERR_INTERNAL_ERROR;
1303 goto smerr;
1304 }
1305 //if(RSB_SOME_ERROR(errval = rsb_mtx_elemental_scale(cmatrix,&three)))
1306 if(RSB_SOME_ERROR(errval = rsb__do_upd_vals(cmatrix,RSB_ELOPF_MUL,&three)))
1307 { RSB_ERROR(RSB_ERRM_ES); goto smerr; }
1308 coo.typecode=smatrix->typecode; coo.nnz=smatrix->nnz;
1309 RSB_DO_FLAG_ADD(smatrix->flags,RSB_FLAG_EXTERNALLY_ALLOCATED_ARRAYS);
1310 errval = rsb__do_switch_rsb_mtx_to_coo(smatrix,&coo.VA,&coo.IA,&coo.JA,RSB_FLAG_SORTED_INPUT);
1311 if(RSB_SOME_ERROR(errval)){ RSB_ERROR(RSB_ERRM_ES); goto smerr; }
1312 smatrix=NULL;
1313 coc.typecode=cmatrix->typecode; coc.nnz=cmatrix->nnz;
1314 RSB_DO_FLAG_ADD(cmatrix->flags,RSB_FLAG_EXTERNALLY_ALLOCATED_ARRAYS);
1315 errval = rsb__do_switch_rsb_mtx_to_coo(cmatrix,&coc.VA,&coc.IA,&coc.JA,RSB_FLAG_SORTED_INPUT);
1316 if(RSB_SOME_ERROR(errval)){ RSB_ERROR(RSB_ERRM_ES); goto smerr; }
1317 cmatrix=NULL;
1318 if(RSB_SOME_ERROR(errval)) { RSB_ERROR(RSB_ERRM_ES); goto smerr; }
1319 //if(trans == RSB_TRANSPOSITION_T)RSB_SWAP(rsb_coo_idx_t*,coo.IA,coo.JA);
1320 //if(trans == RSB_TRANSPOSITION_C)errval = rsb__util_do_conjugate(coo.VA,coo.typecode,coo.nnz);
1321 //coo.VA=coc.VA=NULL;
1322 #if RSB_ALLOW_EMPTY_MATRICES
1323 if((!isempty) || (!rsb__are_coo_matrices_both_empty(&coo,RSB_FLAG_NOFLAGS,&coc,RSB_FLAG_NOFLAGS)))
1324 #endif /* RSB_ALLOW_EMPTY_MATRICES */
1325 if(!rsb__are_coo_matrices_equal(&coo,&coc))
1326 { errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR("matrices do not match!\n"); goto smerr; }
1327 //
1328 smerr:
1329 rsb__destroy_coo_matrix_t(&coo);
1330 rsb__destroy_coo_matrix_t(&coc);
1331 RSB_MTX_FREE(cmatrix);
1332 RSB_MTX_FREE(smatrix);
1333 if(RSB_SOME_ERROR(errval)) { RSB_ERROR("some error occurred while testing matrices sum functionality\n"); goto err; }
1334 #if RSB_WANT_AUTOTUNING_TESTING
1335 #if RSB_ALLOW_INTERNAL_GETENVS
1336 if( getenv("RSB_BMT_AUTOTUNE") )
1337 do_tune_test = rsb__util_atoi(getenv("RSB_BMT_AUTOTUNE"));
1338 #endif /* RSB_ALLOW_INTERNAL_GETENVS */
1339 if( do_tune_test > 0 )
1340 if( nnz > 0 && nnz < 100 ) /* FIXME: these limits here are only for time reasons */
1341 {
1342 struct rsb_mtx_t * mtxOp = NULL;
1343 #if !RSB_AT_DESTROYS_MTX
1344 struct rsb_mtx_t * mtxQp = NULL;
1345 #endif /* RSB_AT_DESTROYS_MTX */
1346 rsb_real_t *sfp = NULL;
1347 rsb_int_t *tnp = NULL;
1348 rsb_int_t oitmax = 0;
1349 rsb_time_t tmax=/*RSB_TIME_ZERO*/0.003;
1350 rsb_trans_t transA = trans;
1351 const void * alphap = NULL;
1352 // struct rsb_mtx_t * mtxAp = ;
1353 rsb_coo_idx_t nrhs = 1;
1354 rsb_flags_t order = RSB_FLAG_WANT_COLUMN_MAJOR_ORDER;
1355 const void * Bp = NULL;
1356 rsb_nnz_idx_t ldB = 0;
1357 const void * betap = NULL;
1358 void * Cp = NULL;
1359 rsb_nnz_idx_t ldC = 0;
1360
1361 // mtxOp = rsb__mtx_clone_simple(mtxAp);
1362 errval = rsb__mtx_clone(&mtxOp,RSB_NUMERICAL_TYPE_SAME_TYPE,RSB_TRANSPOSITION_N,NULL,mtxAp,RSB_FLAG_IDENTICAL_FLAGS);
1363
1364 ldB = rsb__do_get_rows_of(mtxOp,transA);
1365 ldC = rsb__do_get_columns_of(mtxOp,transA);
1366 #if !RSB_AT_DESTROYS_MTX
1367 mtxQp = mtxOp;
1368 #endif /* RSB_AT_DESTROYS_MTX */
1369 if( RSB_SOME_ERROR( errval = rsb__do_tune_spmm(&mtxOp, sfp, tnp, oitmax, tmax, transA, alphap, NULL, nrhs, order, Bp, ldB, betap, Cp, ldC)))
1370 {
1371 errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR("rsb_tune_spmm failed !\n"); goto sterr;
1372 }
1373 #if !RSB_AT_DESTROYS_MTX
1374 if(mtxQp != mtxOp) RSB_MTX_FREE(mtxQp);
1375 #endif /* RSB_AT_DESTROYS_MTX */
1376 #if 0
1377 RSB_MTX_FREE(mtxOp);
1378 mtxOp = rsb__mtx_clone_simple(mtxAp);
1379 #else
1380 errval = rsb__mtx_clone(&mtxOp,RSB_NUMERICAL_TYPE_SAME_TYPE,RSB_TRANSPOSITION_N,NULL,mtxAp,RSB_FLAG_IDENTICAL_FLAGS);
1381 #endif
1382 #if !RSB_AT_DESTROYS_MTX
1383 mtxQp = mtxOp;
1384 #endif /* RSB_AT_DESTROYS_MTX */
1385 #if RSB_ALLOW_EMPTY_MATRICES
1386 if(!isempty)
1387 #endif
1388 if(stype==blas_upper_triangular || stype==blas_lower_triangular )
1389 if( RSB_SOME_ERROR( errval = rsb__do_tune_spsm(&mtxOp, sfp, tnp, oitmax, tmax, transA, alphap, NULL, nrhs, order, Bp, ldB, betap, Cp, ldC)))
1390 {
1391 errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR("rsb_tune_spsm failed !\n"); goto sterr;
1392 }
1393 #if !RSB_AT_DESTROYS_MTX
1394 if(mtxQp != mtxOp) RSB_MTX_FREE(mtxQp);
1395 #endif /* RSB_AT_DESTROYS_MTX */
1396 RSB_MTX_FREE(mtxOp);
1397 sterr:
1398 if(RSB_SOME_ERROR(errval)) { RSB_ERROR("...\n"); goto err; }
1399 }
1400 #endif /* RSB_WANT_AUTOTUNING_TESTING */
1401 }
1402 }
1403 }
1404 #endif
1405 #if 1
1406 cmatrix = NULL;
1407 errval = rsb__mtx_clone(&cmatrix,RSB_NUMERICAL_TYPE_SAME_TYPE,RSB_TRANSPOSITION_N,NULL,mtxAp,RSB_FLAG_IDENTICAL_FLAGS);
1408 // cmatrix = rsb__mtx_clone_simple(mtxAp);
1409
1410 if( cmatrix == NULL )
1411 {
1412 RSB_ERROR("failed matrix cloning\n");
1413 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1414 }
1415 else
1416 {
1417 struct rsb_coo_matrix_t coo,csr;
1418 rsb_flags_t cflags = RSB_DO_FLAG_FILTEROUT(cmatrix->flags,RSB_FLAG_FORTRAN_INDICES_INTERFACE);
1419
1420 RSB_BZERO_P(&coo);
1421 RSB_BZERO_P(&csr);
1422
1423 RSB_DO_FLAG_ADD(cflags,RSB_FLAG_SORTED_INPUT); // NEW, TO SPEEDUP THIS CODE (WEAKENS THE TESTING EFFECTIVENESS, THOUGH)
1424 if(!rsb__mtx_chk(cmatrix))
1425 {
1426 RSB_ERROR("cloned matrix is not built correctly\n");
1427 errval = RSB_ERR_INTERNAL_ERROR; goto converr;
1428 }
1429 if(!cmatrix->nnz)
1430 goto cmedone;
1431 RSB_INIT_CXX_FROM_MTX(&coo,cmatrix);
1432 csr=coo;
1433 if((rsb__allocate_coo_matrix_t(&coo)!=&coo) || (rsb__allocate_coo_matrix_t(&csr)!=&csr))
1434 {
1435 RSB_ERROR("allocaton problem\n");
1436 errval = RSB_ERR_INTERNAL_ERROR; goto converr;
1437 }
1438 coo.nnz=csr.nnz=mtxAp->nnz;
1439 getcsrcooagain:
1440 // RSB -> COO
1441 errval = rsb__do_get_coo_noalloc(mtxAp,coo.VA,coo.IA,coo.JA,NULL,cflags);
1442 if(RSB_SOME_ERROR(errval))
1443 {
1444 RSB_ERROR("coo extraction problems\n");
1445 errval = RSB_ERR_INTERNAL_ERROR; goto converr;
1446 }
1447 // RSB -> CSR
1448 errval = rsb__do_get_csr(typecode,mtxAp,csr.VA,csr.IA,csr.JA,cflags);
1449 //errval = rsb__do_get_csr(typecode,mtxAp,csr.VA,csr.IA,csr.JA,cflags|RSB_FLAG_DEFAULT_CSR_MATRIX_FLAGS);
1450 //errval = rsb__do_get_csr(typecode,mtxAp,csr.VA,csr.IA,csr.JA,RSB_FLAG_DEFAULT_CSR_MATRIX_FLAGS);
1451 if(RSB_SOME_ERROR(errval))
1452 {
1453 RSB_ERROR("csr extraction problems\n");
1454 errval = RSB_ERR_INTERNAL_ERROR; goto converr;
1455 }
1456 // CSR -> COO
1457 errval = rsb__util_uncompress_row_pointers_array(csr.IA,csr.nr,cflags,cflags,csr.IA);
1458 if(RSB_SOME_ERROR(errval))
1459 {
1460 RSB_ERROR("coo->csr conversion problems\n");
1461 errval = RSB_ERR_INTERNAL_ERROR; goto converr;
1462 }
1463 // let's check if 'csr' was converted in coo
1464 if(!rsb__are_coo_matrices_equal(&coo,&csr))
1465 {
1466 RSB_ERROR("no match in coo/csr extractors\n");
1467 errval = RSB_ERR_INTERNAL_ERROR; goto converr;
1468 }
1469 // COO -> CSR
1470 errval = rsb__util_compress_to_row_pointers_array(NULL,csr.nnz,csr.nr,cflags,cflags,csr.IA);
1471 if(RSB_SOME_ERROR(errval))
1472 {
1473 RSB_ERROR("csr->coo conversion failed!\n");
1474 errval = RSB_ERR_INTERNAL_ERROR; goto converr;
1475 }
1476 else
1477 /* FIXME: checks are missing for the following ! */
1478 if(!RSB_DO_FLAG_HAS(cflags,RSB_FLAG_USE_HALFWORD_INDICES))
1479 if(RSB_SOME_ERROR(errval = rsb__csr_chk(csr.IA,csr.JA,csr.nr,csr.nc,csr.nnz,0)))
1480 {
1481 RSB_ERROR("csr->coo conversion produced corrupt results!\n");
1482 errval = RSB_ERR_INTERNAL_ERROR; goto converr;
1483 }
1484 // RSB -> COO
1485 // FIXME:
1486 // errval = rsb__do_switch_recursive_in_place_matrix_to_in_place_coo_unsorteda..
1487 // if ..
1488 // rsb__destroy_coo_matrix_t(&icoo);
1489 RSB_MTX_FREE(cmatrix);
1490 // CSR -> RSB
1491
1492 if(incX==1 && incB==1) if(alphai==0 && betai==0) /* agnostic to these parameters */
1493 {
1494 kmatrix = rsb__do_mtx_alloc_from_csr_const(csr.VA,csr.IA,csr.JA,csr.nnz,csr.typecode,csr.nr,csr.nc,RSB_DEFAULT_ROW_BLOCKING,RSB_DEFAULT_COL_BLOCKING,cflags,&errval);
1495 if((RSB_SOME_ERROR(errval)) || (!kmatrix) || (!rsb__mtx_chk(kmatrix)))
1496 { RSB_ERROR("csr->rsb construction problems\n"); goto err;}
1497 RSB_MTX_FREE(kmatrix);
1498 }
1499
1500 cmatrix = rsb__do_mtx_alloc_from_csr_inplace(csr.VA,csr.IA,csr.JA,csr.nnz,csr.typecode,csr.nr,csr.nc,RSB_DEFAULT_ROW_BLOCKING,RSB_DEFAULT_COL_BLOCKING,cflags,&errval);
1501 if((RSB_SOME_ERROR(errval)) || (!cmatrix) || (!rsb__mtx_chk(cmatrix)))
1502 {
1503 if(RSB_SOME_ERROR(errval))
1504 RSB_ERROR("csr->rsb construction problems\n");
1505 else
1506 if(!cmatrix)
1507 {
1508 errval = RSB_ERR_INTERNAL_ERROR;
1509 RSB_ERROR("csr->rsb construction problems: did not succeed\n");
1510 }
1511 else
1512 {
1513 RSB_ERROR("csr->rsb construction problems: built a corrupted matrix\n");
1514 errval = RSB_ERR_INTERNAL_ERROR;
1515 }
1516 goto err;
1517 }
1518 RSB_DO_FLAG_ADD(cmatrix->flags,RSB_FLAG_EXTERNALLY_ALLOCATED_ARRAYS);
1519 // RSB -> COO
1520 #if 0
1521 errval = rsb__do_switch_rsb_mtx_to_coo(cmatrix,&csr.VA,&csr.IA,&csr.JA,cflags|RSB_FLAG_SORTED_INPUT);
1522 #else
1523 // still broken!
1524 if(nnz<42) /* coverage testing purpose :P */
1525 {
1526 const rsb_nnz_idx_t nnz = cmatrix->nnz; // cmatrix may be freed in-block with csr.VA &co
1527 errval = rsb__do_switch_rsb_mtx_to_coo(cmatrix,&csr.VA,&csr.IA,&csr.JA,RSB_DO_FLAG_FILTEROUT(cflags,RSB_FLAG_SORTED_INPUT));
1528 errval |= rsb__util_sort_row_major_inner(csr.VA,csr.IA,csr.JA,nnz,dim,dim,typecode,RSB_DO_FLAG_FILTEROUT(cflags,RSB_FLAG_SORTED_INPUT));
1529 }
1530 else
1531 errval = rsb__do_switch_rsb_mtx_to_coo(cmatrix,&csr.VA,&csr.IA,&csr.JA,cflags|RSB_FLAG_SORTED_INPUT);
1532 #endif
1533 if((RSB_SOME_ERROR(errval)) || !rsb__are_coo_matrices_equal(&coo,&csr))
1534
1535 {
1536 RSB_ERROR("rsb->coo conversion problems\n");
1537 errval = RSB_ERR_INTERNAL_ERROR; goto converr;
1538 }
1539 // COO -> RSB
1540 cmatrix = rsb__do_mtx_alloc_from_coo_inplace(csr.VA,csr.IA,csr.JA,csr.nnz,csr.typecode,csr.nr,csr.nc,RSB_DEFAULT_ROW_BLOCKING,RSB_DEFAULT_COL_BLOCKING,cflags,&errval);
1541 if((RSB_SOME_ERROR(errval)) || (!cmatrix) || (!rsb__mtx_chk(cmatrix)))
1542 {
1543 RSB_ERROR("csr->coo conversion problems\n");
1544 errval = RSB_ERR_INTERNAL_ERROR; goto converr;
1545 }
1546 // RSB -> CSR
1547 RSB_DO_FLAG_ADD(cmatrix->flags,RSB_FLAG_EXTERNALLY_ALLOCATED_ARRAYS);
1548 errval = rsb__do_switch_rsb_mtx_to_csr_sorted(cmatrix,&csr.VA,&csr.IA,&csr.JA,cflags);
1549 cmatrix=NULL;
1550 if(RSB_SOME_ERROR(errval))
1551 {
1552 RSB_ERROR("coo->csr conversion problems\n");
1553 errval = RSB_ERR_INTERNAL_ERROR; goto converr;
1554 }
1555 errval = rsb__util_uncompress_row_pointers_array(csr.IA,csr.nr,cflags,cflags,csr.IA);
1556 if((RSB_SOME_ERROR(errval)) || !rsb__are_coo_matrices_equal(&coo,&csr))
1557 {
1558 RSB_ERROR("csr->coo conversion problems\n");
1559 // rsb__debug_print_index_vectors_diff(coo.IA,csr.IA,csr.nnz,RSB_VECTORS_DIFF_DISPLAY_N_SMALL);
1560 // rsb__debug_print_index_vectors_diff(coo.JA,csr.JA,csr.nnz,RSB_VECTORS_DIFF_DISPLAY_N_SMALL);
1561 // rsb__debug_print_vectors_diff(coo.VA,csr.VA,csr.nnz,typecode,1,1,RSB_VECTORS_DIFF_DISPLAY_N_SMALL);
1562 errval = RSB_ERR_INTERNAL_ERROR; goto converr;
1563 }
1564
1565 if(!RSB_DO_FLAG_HAS(cflags,RSB_FLAG_FORTRAN_INDICES_INTERFACE))
1566 {
1567 RSB_DO_FLAG_ADD(cflags,RSB_FLAG_FORTRAN_INDICES_INTERFACE);
1568 goto getcsrcooagain;
1569 }
1570 else
1571 RSB_DO_FLAG_DEL(cflags,RSB_FLAG_FORTRAN_INDICES_INTERFACE);
1572 converr:
1573 rsb__destroy_coo_matrix_t(&coo);
1574 rsb__destroy_coo_matrix_t(&csr);
1575 cmedone:
1576 RSB_MTX_FREE(cmatrix);
1577 if(RSB_SOME_ERROR(errval))
1578 RSB_PERR_GOTO(err,RSB_ERRM_ES);
1579 }
1580 #endif
1581 #if RSB_ALLOW_INTERNAL_GETENVS
1582 if(! ( getenv("RSB_BMT_SPSV") && rsb__util_atoi(getenv("RSB_BMT_SPSV")) == 0 ) )
1583 #endif /* RSB_ALLOW_INTERNAL_GETENVS */
1584 if(!RSB_DO_FLAG_HAS(mtxAp->flags,RSB_FLAG_SYMMETRIC))
1585 {
1586 for(mmi=0;mmi< (dim<msmd?3:2) ;++mmi)
1587 if(! (mmi==1 && ((incX!= 1) || (incB!=1) ) ))
1588 {
1589 const rsb_int nrhs=1;
1590
1591 /* TODO : should fill X and B with sentinel values ! */
1592 if(RSB_SOME_ERROR(rsb__cblas_Xscal(typecode,dim,&zero,X,incX)) || RSB_SOME_ERROR(rsb__fill_with_ones (B,typecode,dim,incB)))
1593 { errval = RSB_ERR_INTERNAL_ERROR; goto err; }
1594
1595 if(mmi==0)
1596 if( rsb__BLAS_Xusmv(transT,alpha,T,B,incB,beta,X,incX) != RSB_BLAS_NO_ERROR )
1597 {
1598 errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("error while performing Unsymmetric USMV\n"); goto err;
1599 }
1600
1601 if(mmi==1)
1602 if(RSB_SOME_ERROR( rsb__do_spmm(trans,alpha,mtxAp,nrhs,RSB_FLAG_WANT_COLUMN_MAJOR_ORDER,B,dim,beta,X,dim,RSB_OP_FLAG_DEFAULT)) )
1603 {
1604 errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("Unsymmetric USMV failed!\n"); goto err;
1605 }
1606
1607 if(mmi==2)
1608 if(RSB_SOME_ERROR( rsb__do_spmv_general(trans,alpha,mtxAp,B,incB,beta,X,incX,RSB_OP_FLAG_WANT_SERIAL RSB_DEFAULT_OUTER_NRHS_SPMV_ARGS)))
1609 {
1610 errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("Unsymmetric USMV failed!\n"); goto err;
1611 }
1612
1613 if(isinvertible)
1614 {
1615 if(mmi==0)
1616 if( rsb__BLAS_Xussv(transT,alpha_inv,T,X,incX) != RSB_BLAS_NO_ERROR )
1617 {
1618 errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("error while performing USSV\n"); goto err;
1619 }
1620
1621 if(mmi==1)
1622 if(RSB_SOME_ERROR( rsb_spsm(trans,alpha_inv,mtxAp,nrhs,RSB_FLAG_WANT_COLUMN_MAJOR_ORDER,alpha_inv,X,dim,X,dim)) )
1623 {
1624 errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("error while performing USSV\n"); goto err;
1625 }
1626
1627 if(mmi==2)
1628 if( rsb__do_spsv(trans,alpha_inv,mtxAp,X,incX,X,incX) != RSB_BLAS_NO_ERROR )
1629 {
1630 errval = RSB_ERR_INTERNAL_ERROR;RSB_ERROR("error while performing USSV\n"); goto err;
1631 }
1632 }
1633 if(!isinvertible)
1634 rsb__cblas_Xscal(typecode,dim,&zero,B,incB);
1635 if(stype != blas_general)
1636 if(stype != blas_lower_hermitian) /* FIXME: complete this case */
1637 if(stype != blas_upper_hermitian) /* FIXME: complete this case */
1638 if( RSB_SOME_ERROR(rsb__do_are_same(B,X,dim,typecode,incB,incX)) )
1639 {
1640 RSB_ERROR("failed post combined USMV-USSV check!\n");
1641 rsb__debug_print_vectors_diff(B,X,dim,typecode,incB,incX,RSB_VECTORS_DIFF_DISPLAY_N);
1642 errval = RSB_ERR_INTERNAL_ERROR;
1643 goto err;
1644 }
1645 }
1646 }
1647
1648 if(betai > 0 || alphai > 0)
1649 goto err; /* only the previous tests were affected by alpha and beta */
1650 /*
1651 TODO: complete the following ...
1652
1653 if( rsb__BLAS_Xusget_rows_sums(T,rs,transT) != RSB_BLAS_NO_ERROR )
1654 {
1655 RSB_ERROR("error getting rows sum!\n");
1656 errval = RSB_ERR_INTERNAL_ERROR;
1657 goto err;
1658 }
1659 else
1660 {
1661 */
1662
1663 /* TODO: need parameters scan here: */
1664 if( RSB_SOME_ERROR(errval = rsb__do_upd_vals(mtxAp,RSB_ELOPF_NEG,NULL))) { RSB_ERROR("Failed negating.\n"); goto smerr; }
1665 if( RSB_SOME_ERROR(errval = rsb__do_upd_vals(mtxAp,RSB_ELOPF_NEG,NULL))) { RSB_ERROR("Failed negating.\n"); goto smerr; }
1666 if(RSB_SOME_ERROR(rsb__cblas_Xscal(typecode,1,&zero,inrm,1))){ errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR("!\n"); goto err; }
1667
1668 rsb__util_set_area_to_converted_integer(inrm,typecode,0);
1669
1670 if( rsb__BLAS_Xusget_infinity_norm(T,inrm,transT) != RSB_BLAS_NO_ERROR )
1671 {
1672 RSB_ERROR("error getting infinity norm!\n");
1673 errval = RSB_ERR_INTERNAL_ERROR;
1674 goto err;
1675 }
1676 else
1677 {
1678 if(is_really_empty)
1679 {
1680 rsb__util_set_area_to_converted_integer(D,typecode,0);
1681 }
1682 else
1683 {
1684 if(isempty)
1685 rsb__util_set_area_to_converted_integer(D,typecode,1 );
1686 else
1687 rsb__util_set_area_to_converted_integer(D,typecode,dim);
1688 }
1689 if( RSB_SOME_ERROR(rsb__do_are_same(inrm,D,1,typecode,1,1)) )
1690 {
1691 RSB_ERROR("matrix norm is not what was expected!\n");
1692 rsb__debug_print_vectors_diff(inrm,D,1,typecode,1,1,RSB_VECTORS_DIFF_DISPLAY_N);
1693 errval = RSB_ERR_INTERNAL_ERROR;
1694 goto err;
1695 }
1696
1697 mtxAp = rsb__BLAS_inner_matrix_retrieve(T);
1698 if(incX==1 && incB==1) if(alphai==0 && betai==0) /* agnostic to these parameters */
1699 if((!isempty) && !(mtxAp->nnz == 0 && diaga[diagi]==blas_unit_diag))
1700 {
1701 /* FIXME: Frobenius norm (RSB_EXTF_NORM_TWO) is untested! */
1702 rsb__util_set_area_to_converted_integer(D,typecode,1);
1703 rsb__do_upd_vals(mtxAp,RSB_ELOPF_MUL,D);
1704 rsb__util_set_area_to_converted_integer(inrm_,typecode,0);
1705 RSB_LSTPROBE(rsb__BLAS_Xusget_infinity_norm(T,inrm_,transT),"");
1706 RSB_LSTPROBE(rsb__do_are_same(inrm_,inrm,1,typecode,1,1),"");
1707 rsb__util_set_area_to_converted_integer(D,typecode,2);
1708
1709 rsb__do_upd_vals(mtxAp,RSB_ELOPF_MUL,D);
1710 rsb__util_set_area_to_converted_integer(inrm_,typecode,0);
1711 RSB_LSTPROBE(rsb__BLAS_Xusget_infinity_norm(T,inrm_,transT),"");
1712 if(mtxAp->nnz && diaga[diagi]==blas_unit_diag)
1713 rsb__util_increase_by_one(inrm_,0,mtxAp->typecode);
1714 rsb__util_vector_div(inrm_,D,typecode,1);
1715 RSB_LSTPROBE(rsb__do_are_same(inrm_,inrm,1,typecode,1,1),"");
1716
1717 rsb__do_upd_vals(mtxAp,RSB_ELOPF_DIV,D);
1718 rsb__util_set_area_to_converted_integer(inrm_,typecode,0);
1719 RSB_LSTPROBE(rsb__BLAS_Xusget_infinity_norm(T,inrm_,transT),"");
1720 RSB_LSTPROBE(rsb__do_are_same(inrm_,inrm,1,typecode,1,1),"");
1721 /* TODO: there are many more subcases for rsb__mtx_clone! */
1722 // if( ((cmatrix = rsb__mtx_clone_simple(mtxAp))!=NULL))
1723 cmatrix = NULL;
1724 errval = rsb__mtx_clone(&cmatrix,RSB_NUMERICAL_TYPE_SAME_TYPE,RSB_TRANSPOSITION_N,NULL,mtxAp,RSB_FLAG_IDENTICAL_FLAGS);
1725 if( (cmatrix !=NULL))
1726 {
1727 RSB_LSTPROBE(rsb__mtx_clone(&cmatrix,RSB_NUMERICAL_TYPE_SAME_TYPE,RSB_TRANSPOSITION_T,NULL,cmatrix,RSB_FLAG_IDENTICAL_FLAGS),"");
1728 RSB_LSTPROBE(rsb__mtx_clone(&cmatrix,RSB_NUMERICAL_TYPE_SAME_TYPE,RSB_TRANSPOSITION_T,NULL,cmatrix,RSB_FLAG_IDENTICAL_FLAGS),"");
1729 if(RSB_IS_MATRIX_TYPE_COMPLEX(typecode))/* FIXME: shall fix many vector-operating routines, first */
1730 {
1731 /* TODO: more checking */
1732 RSB_LSTPROBE(rsb__mtx_clone(&cmatrix,RSB_NUMERICAL_TYPE_SAME_TYPE,RSB_TRANSPOSITION_C,NULL,cmatrix,RSB_FLAG_IDENTICAL_FLAGS),"");
1733 RSB_LSTPROBE(rsb__mtx_clone(&cmatrix,RSB_NUMERICAL_TYPE_SAME_TYPE,RSB_TRANSPOSITION_C,NULL,cmatrix,RSB_FLAG_IDENTICAL_FLAGS),"");
1734 }
1735 RSB_MTX_FREE(cmatrix);
1736 }
1737
1738 if(dim>0)
1739 {
1740 IA[0]=dim-1; JA[0]=dim-1;
1741 RSB_LSTPROBE(rsb__do_get_elements(mtxAp,VA,IA,JA,1,mtxAp->flags),"");
1742 if(!(diaga[diagi]==blas_unit_diag))
1743 RSB_LSTPROBE(rsb__do_set_elements(mtxAp,VA,IA,JA,1,mtxAp->flags),"");
1744 }
1745 if(dim>1)
1746 {
1747 const rsb_int mmudim=100;
1748 IA[0]=dim-1; JA[0]=0;
1749 /* TODO: shall check value! */
1750 if(stype==blas_upper_triangular)
1751 {
1752 RSB_LSTPROBE(rsb__do_get_elements(mtxAp,VA,JA,IA,1,mtxAp->flags),"");
1753 if(!(diaga[diagi]==blas_unit_diag))
1754 RSB_LSTPROBE(rsb__do_set_elements(mtxAp,VA,JA,IA,1,mtxAp->flags),"");
1755 RSB_LSTPROBI(rsb__do_get_elements/*rsb__do_get_elements*/(mtxAp,VA,IA,JA,1,mtxAp->flags),"");
1756 RSB_LSTPROBI(rsb__do_set_elements/*rsb__do_set_elements*/(mtxAp,VA,IA,JA,1,mtxAp->flags),"");
1757 }
1758 else
1759 {
1760 RSB_LSTPROBE(rsb__do_get_elements(mtxAp,VA,IA,JA,1,mtxAp->flags),"");
1761 if(!(diaga[diagi]==blas_unit_diag))
1762 RSB_LSTPROBE(rsb__do_set_elements(mtxAp,VA,IA,JA,1,mtxAp->flags),"");
1763 RSB_LSTPROBI(rsb__do_get_elements/*rsb__do_get_elements*/(mtxAp,VA,JA,IA,1,mtxAp->flags),"");
1764 RSB_LSTPROBI(rsb__do_set_elements/*rsb__do_set_elements*/(mtxAp,VA,JA,IA,1,mtxAp->flags),"");
1765 }
1766
1767 if(dim>1 && mtxAp->nnz>0 )
1768 if(dim < mmudim )
1769 {
1770 struct rsb_mtx_t*LU[]={NULL,NULL};
1771 RSB_LSTPROBE(rsb__do_get_preconditioner(LU,mtxAp,RSB_PRECF_ILU0,NULL),"");
1772 RSB_MTX_FREE(LU[0]);
1773 RSB_MTX_FREE(LU[1]);
1774 }
1775 #if 0
1776 if(dim < mmudim && ( cmatrix = rsb__mtx_clone_simple(mtxAp)) !=NULL)
1777 {
1778 if(stype==blas_upper_triangular)
1779 {
1780 RSB_LSTPROBE(rsb_mtx_set_values_pattern_changing(&cmatrix,VA,IA,JA,1,mtxAp->flags),"");
1781 }
1782 else
1783 {
1784 RSB_LSTPROBE(rsb_mtx_set_values_pattern_changing(&cmatrix,VA,JA,IA,1,mtxAp->flags),"");
1785 }
1786 RSB_MTX_FREE(cmatrix);
1787 }
1788 #endif
1789 }
1790
1791 }
1792 }
1793
1794 RSB_LSTPROBE(rsb__do_elemental_binop(mtxAp, RSB_ELOPF_POW, &three),""); /* FIXME: shall test systematically all the others as well !*/
1795
1796 #if RSB_ALLOW_INTERNAL_GETENVS
1797 if(! ( getenv("RSB_BMT_GET") && rsb__util_atoi(getenv("RSB_BMT_GET")) == 0 ) )
1798 #endif /* RSB_ALLOW_INTERNAL_GETENVS */
1799 if(1)
1800 {
1801 rsb_coo_idx_t m=dim,k=dim;
1802 rsb_coo_idx_t rc=m/3,fr = rc,lr = RSB_MIN(m-1,2*rc),ri,
1803 cc=k/3,fc=cc,lc = RSB_MIN(k-1,2*cc),ci;
1804 rsb_nnz_idx_t bnnz=0,cnnz=0,off=0;
1805 const void*vp=NULL;
1806 bnnz = rsb__do_get_block_nnz(mtxAp,fr,lr,fc,lc,RSB_FLAG_C_INDICES_INTERFACE,&errval);
1807 // FIXME: TODO: should also test rsb__do_get_block_sparse()
1808 if(RSB_SOME_ERROR(errval))
1809 {
1810 RSB_ERROR("sparse subblocks nnz count mechanisms seem broken.\n");
1811 rsb__do_perror(NULL,errval);
1812 errval = RSB_ERR_INTERNAL_ERROR;
1813 goto err;
1814 }
1815 for(ri=fr;ri<=lr;++ri)
1816 for(ci=fc;ci<=lc;++ci)
1817 cnnz+=(rsb__do_coo_element_inner_address(mtxAp,ri,ci)!=NULL);
1818 if(bnnz!=cnnz)
1819 {
1820 RSB_ERROR("sparse subblocks nnz count mechanisms seem broken (%d vs %d counted in (%d,%d)..(%d,%d)).\n",bnnz,cnnz,fr,fc,lr,lc);
1821 errval = RSB_ERR_INTERNAL_ERROR;
1822 goto err;
1823 }
1824
1825 rsb__util_coo_array_set(IA,nnz,RSB_MARKER_COO_VALUE);
1826 rsb__util_coo_array_set(JA,nnz,RSB_MARKER_COO_VALUE);
1827 errval = rsb__do_get_block_sparse(mtxAp,VA,IA,JA,fr,lr,fc,lc,NULL,NULL,&cnnz,RSB_FLAG_C_INDICES_INTERFACE);
1828 if(RSB_SOME_ERROR(errval))
1829 {
1830 RSB_ERROR("sparse subblocks nnz get mechanisms seem broken.\n");
1831 rsb__do_perror(NULL,errval);
1832 errval = RSB_ERR_INTERNAL_ERROR;
1833 goto err;
1834 }
1835
1836 if(bnnz!=cnnz)
1837 {
1838 RSB_ERROR("sparse subblocks nnz get mechanisms seem broken (%d vs %d counted in (%d,%d)..(%d,%d)).\n",bnnz,cnnz,fr,fc,lr,lc);
1839 errval = RSB_ERR_INTERNAL_ERROR;
1840 goto err;
1841 }
1842
1843 for(off=0;off<bnnz;++off)
1844 if((vp = rsb__do_coo_element_inner_address(mtxAp,IA[off],JA[off]))!=NULL)
1845 {
1846 if(RSB_VA_MEMCMP(vp,0,VA,off,mtxAp->el_size))
1847 {
1848 RSB_ERROR("address of (%d,%d)@%d extracted from sparse seems not the right one\n",IA[off],JA[off],off);
1849 errval = RSB_ERR_INTERNAL_ERROR;
1850 goto err;
1851 }
1852 }
1853 else
1854 {
1855 RSB_ERROR("an element (%d,%d)@%d extracted from sparse seems not present\n",IA[off],JA[off],off);
1856 errval = RSB_ERR_INTERNAL_ERROR;
1857 goto err;
1858 }
1859 }
1860
1861 #if RSB_ALLOW_INTERNAL_GETENVS
1862 if(! ( getenv("RSB_BMT_SCALE") && rsb__util_atoi(getenv("RSB_BMT_SCALE")) == 0 ) )
1863 #endif /* RSB_ALLOW_INTERNAL_GETENVS */
1864 {
1865 if(RSB_SOME_ERROR(rsb__cblas_Xscal(typecode,dim,&zero,D,incD))){ errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR("!\n"); goto err; }
1866 if( rsb__BLAS_Xusget_diag(T,D) != RSB_BLAS_NO_ERROR )
1867 {
1868 RSB_ERROR("!\n");
1869 errval = RSB_ERR_INTERNAL_ERROR;
1870 goto err;
1871 }
1872 else
1873 {
1874 if(is_really_empty)
1875 {if(RSB_SOME_ERROR(rsb__cblas_Xscal(typecode,dim,&zero,B,incB))){ errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL); goto err; }}
1876 else
1877 {if(RSB_SOME_ERROR(rsb__fill_with_ones(B,typecode,dim,incB))){ errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL); goto err; }}
1878 if(RSB_SOME_ERROR(rsb__do_are_same(D,B,dim,typecode,incD,incB)))
1879 {
1880 RSB_ERROR("diagonal vector not what expected!\n");
1881 rsb__debug_print_vectors_diff(D,B,dim,typecode,incD,incB,RSB_VECTORS_DIFF_DISPLAY_N);
1882 errval = RSB_ERR_INTERNAL_ERROR;
1883 goto err;
1884 }
1885 }
1886
1887 incB=1; /* FIXME: now on, no stride */
1888 incD=1; /* FIXME: now on, no stride */
1889 if(RSB_SOME_ERROR(rsb__fill_with_increasing_values(B,typecode,dim))){ errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL); goto err; } /* rsb__fill_with_increasing_values..) -> ,every) */
1890 if(rsb__BLAS_Xusrows_scale(T,B,transT) != RSB_BLAS_NO_ERROR)
1891 {
1892 RSB_ERROR("!\n");
1893 errval = RSB_ERR_INTERNAL_ERROR;
1894 goto err;
1895 }
1896
1897 if(RSB_SOME_ERROR(rsb__cblas_Xscal(typecode,dim,&zero,D,incD))){ errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL); goto err; }
1898 if( rsb__BLAS_Xusget_diag(T,D) != RSB_BLAS_NO_ERROR )
1899 {
1900 RSB_ERROR("!\n");
1901 errval = RSB_ERR_INTERNAL_ERROR;
1902 goto err;
1903 }
1904 else
1905 if(diaga[diagi]==blas_non_unit_diag && !isempty) // diagonal implicit won't be scaled :)
1906 {
1907 rsb_nnz_idx_t n;
1908 if( RSB_SOME_ERROR(rsb__do_are_same(D,B,dim,typecode,incD,incB)) )
1909 {
1910 RSB_ERROR("!\n");
1911 rsb__debug_print_vectors_diff(D,B,dim,typecode,incD,incB,RSB_VECTORS_DIFF_DISPLAY_N);
1912 errval = RSB_ERR_INTERNAL_ERROR;
1913 goto err;
1914 }
1915 if(RSB_SOME_ERROR(rsb__cblas_Xscal(typecode,dim,&zero,B,incB))){ errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL); goto err; }
1916 for(n=0;n<dim;++n) rsb__BLAS_Xusget_element(T,n,n,((rsb_char_t*)B)+el_size*n*incB);
1917 if( RSB_SOME_ERROR(rsb__do_are_same(D,B,dim,typecode,incD,incB)) )
1918 {
1919 RSB_ERROR("!\n");
1920 rsb__debug_print_vectors_diff(D,B,dim,typecode,incD,incB,RSB_VECTORS_DIFF_DISPLAY_N);
1921 errval = RSB_ERR_INTERNAL_ERROR;
1922 goto err;
1923 }
1924 if(RSB_SOME_ERROR(rsb__fill_with_increasing_values(B,typecode,dim))){errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL); goto err; }
1925 for(n=0;n<dim;++n) rsb__BLAS_Xusset_element(T,n,n,((rsb_char_t*)B)+el_size*n*incB);
1926 if(RSB_SOME_ERROR(rsb__cblas_Xscal(typecode,dim,&zero,D,incD))){ errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL); goto err; }
1927 for(n=0;n<dim;++n) rsb__BLAS_Xusget_element(T,n,n,((rsb_char_t*)D)+el_size*n*incB);
1928 if( RSB_SOME_ERROR(rsb__do_are_same(D,B,dim,typecode,incD,incB)) )
1929 {
1930 RSB_ERROR("!\n");
1931 rsb__debug_print_vectors_diff(D,B,dim,typecode,incD,incB,RSB_VECTORS_DIFF_DISPLAY_N);
1932 errval = RSB_ERR_INTERNAL_ERROR;
1933 goto err;
1934 }
1935 }
1936
1937 }
1938 //if((rnnz = rsb__dodo_get_rows_nnz(mtxAp,0,dim-1,RSB_FLAG_C_INDICES_INTERFACE,NULL))!=ndnnz)
1939 rnnz=10;
1940 if(!isempty)/* FIXME */
1941 if((rsb__BLAS_Xusget_rows_nnz(T,0,dim-1,&rnnz)!=RSB_BLAS_NO_ERROR) || (rnnz!=ndnnz))
1942 {
1943 RSB_ERROR("Mismatch in the extracted rows nonzeroes counts vs non diagonal nonzero count: %d != %d\n",(int)rnnz,(int)ndnnz);
1944 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1945 }
1946
1947 rnnz=0;
1948 if(!isempty)/* FIXME */
1949 if(rsb__BLAS_Xusget_matrix_nnz(T,&rnnz)!=RSB_BLAS_NO_ERROR || rnnz!=ndnnz)
1950 {
1951 RSB_ERROR("Mismatch in the wffective matrix counts vs input nonzero count: %d != %d\n",(int)rnnz,(int)ndnnz);
1952 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1953 }
1954
1955 if(RSB_SOME_ERROR(rsb__cblas_Xscal(typecode,nnz,&zero,VA,1))){ errval = RSB_ERR_INTERNAL_ERROR; RSB_ERROR(RSB_ERRM_NL); goto err; }
1956 rsb__util_coo_array_set(IA,nnz,RSB_MARKER_COO_VALUE);
1957 rsb__util_coo_array_set(JA,nnz,RSB_MARKER_COO_VALUE);
1958 rnnz=0;
1959 //if(rsb__do_get_rows_sparse(RSB_TRANSPOSITION_N,NULL,mtxAp,VA,IA,JA,0,dim-1,&rnnz,RSB_FLAG_NOFLAGS|RSB_FLAG_SORT_INPUT))
1960 if(rsb__BLAS_Xusget_rows_sparse(T,VA,IA,JA,&rnnz,0,dim-1)!=RSB_BLAS_NO_ERROR)
1961 {
1962 RSB_ERROR(RSB_ERRM_NL);
1963 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1964 }
1965 else
1966 if(diaga[diagi]==blas_non_unit_diag)
1967 {
1968 rsb_nnz_idx_t n;
1969 for(n=0;n<nnz;++n)
1970 if(IA[n] != JA[n])
1971 if(( rsb__do_coo_element_inner_address(mtxAp,IA[n],JA[n] ) == NULL) ||
1972 (RSB_SOME_ERROR(rsb__do_are_same(rsb__do_coo_element_inner_address(mtxAp,IA[n],JA[n]),
1973 ((rsb_char_t*)VA)+el_size*n,1,typecode,1,1) ) ))
1974 {
1975 RSB_ERROR("@%d (%d,%d) : 0x%x\n",n,IA[n],JA[n],rsb__do_coo_element_inner_address(mtxAp,IA[n],JA[n] ));
1976 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1977 }
1978 }
1979 if(!RSB_DO_TOOFEWNNZFORCSR(nnz,dim))/* we don't want IA overwrite */
1980 {
1981 if(RSB_SOME_ERROR(rsb__do_get_csr(typecode,mtxAp,(void*)(VA),IA,JA,RSB_FLAG_DEFAULT_CSR_MATRIX_FLAGS)))
1982 {
1983 RSB_ERROR(RSB_ERRM_NL);
1984 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1985 }
1986 else
1987 if(diaga[diagi]==blas_non_unit_diag)
1988 {
1989 // TODO: is_csr_sorted ?
1990 rsb_nnz_idx_t n;
1991 rsb_coo_idx_t i;
1992
1993 for(i=0;i<dim;++i)
1994 if(!rsb__util_is_nnz_array_sorted_up(JA+IA[i],IA[i+1]-IA[i]))
1995 {
1996 RSB_ERROR(RSB_ERRM_NL);
1997 errval = RSB_ERR_INTERNAL_ERROR; goto err;
1998 }
1999
2000 for(i=0;i<dim;++i)
2001 for(n=IA[i];n<IA[i+1];++n)
2002 if(JA[n]<0 || JA[n]>=dim)
2003 {
2004 RSB_ERROR(RSB_ERRM_NL);
2005 errval = RSB_ERR_INTERNAL_ERROR; goto err;
2006 }
2007
2008 for(i=0;i<dim;++i)
2009 for(n=IA[i];n<IA[i+1];++n)
2010 if(i != JA[n])
2011 if(( rsb__do_coo_element_inner_address(mtxAp,i,JA[n] ) == NULL) ||
2012 (RSB_SOME_ERROR(rsb__do_are_same(rsb__do_coo_element_inner_address(mtxAp,i,JA[n]),
2013 ((rsb_char_t*)VA)+el_size*n,1,typecode,1,1) ) ))
2014 {
2015 RSB_ERROR("@%d, %d %d : 0x%x\n",n,i,JA[n],rsb__do_coo_element_inner_address(mtxAp,i,JA[n] ));
2016 errval = RSB_ERR_INTERNAL_ERROR; goto err;
2017 }
2018 }
2019 }
2020 if(!RSB_DO_TOOFEWNNZFORCSR(nnz,dim))/* we don't want JA overwrite */
2021 {
2022 if(RSB_SOME_ERROR(rsb__do_get_csc(mtxAp,(void*)(&VA),&JA,&IA)))
2023 {
2024 RSB_ERROR(RSB_ERRM_NL);
2025 errval = RSB_ERR_INTERNAL_ERROR; goto err;
2026 }
2027 else
2028 if(diaga[diagi]==blas_non_unit_diag)
2029 {
2030 rsb_nnz_idx_t n;
2031 rsb_coo_idx_t j;
2032
2033 if( RSB_SOME_ERROR(rsb__csc_chk(JA,IA,dim,dim,JA[dim],0) ) )
2034 {
2035 RSB_ERROR(RSB_ERRM_NL);
2036 errval = RSB_ERR_INTERNAL_ERROR; goto err;
2037 }
2038
2039 for(j=0;j<dim;++j)
2040 for(n=JA[j];n<JA[j+1];++n)
2041 if(j != IA[n])
2042 if(( rsb__do_coo_element_inner_address(mtxAp,IA[n],j) == NULL) ||
2043 ( RSB_SOME_ERROR(rsb__do_are_same(rsb__do_coo_element_inner_address(mtxAp,IA[n],j),
2044 ((rsb_char_t*)VA)+el_size*n,1,typecode,1,1))) )
2045 {
2046 RSB_ERROR("@%d, %d %d : 0x%x\n",n,IA[n],j,rsb__do_coo_element_inner_address(mtxAp,IA[n],j ));
2047 errval = RSB_ERR_INTERNAL_ERROR; goto err;
2048 }
2049 }
2050
2051 if(incX==1 && incB==1) if(alphai==0 && betai==0) /* agnostic to these parameters */
2052 {
2053 rsb_flags_t cflags = RSB_DO_FLAG_FILTEROUT(mtxAp->flags,RSB_FLAG_FORTRAN_INDICES_INTERFACE);
2054 kmatrix = rsb__do_mtx_alloc_from_csc_const(VA,IA,JA,/*nnz*/JA[dim],typecode,dim,dim,RSB_DEFAULT_ROW_BLOCKING,RSB_DEFAULT_COL_BLOCKING,cflags,&errval);
2055 if(RSB_SOME_ERROR(errval) || (!kmatrix) || (!rsb__mtx_chk(kmatrix)))
2056 { RSB_ERROR("csc->rsb construction problems\n"); goto err;}
2057 RSB_MTX_FREE(kmatrix);
2058 }
2059
2060
2061 }
2062 err:
2063 if(errval == RSB_ERR_NO_ERROR)
2064 {
2065 }
2066 else
2067 if(mtxAp && X && B)
2068 {
2069 if(mtxAp->nnz<20)
2070 rsb__do_file_mtx_save(mtxAp,NULL),
2071 RSB_INFO("\n"),
2072 RSB_INFO("actual results vs correct results:\n"),
2073 rsb__debug_print_vectors(X,B,dim,incX,incB,typecode);
2074 else
2075 if( RSB_SOME_ERROR(rsb__do_are_same(B,X,dim,typecode,incB,incX) ))
2076 {
2077 RSB_INFO("actual results vs correct results:\n"),
2078 rsb__debug_print_vectors_diff(X,B,dim,typecode,incX,incB,RSB_VECTORS_DIFF_DISPLAY_N);
2079 errval = RSB_ERR_INTERNAL_ERROR;
2080 }
2081 #if RSB_WANT_VERBOSE_FAILURES
2082 RSB_INFO("Matrix summary:\n");
2083 RSB_INFO_MATRIX_SUMMARY(mtxAp);
2084 RSB_INFO("\n");
2085 rsb__do_print_matrix_stats(mtxAp,RSB_CONST_DUMP_RECURSION_BRIEF,NULL);
2086 RSB_INFO("\n");
2087 #endif /* RSB_WANT_VERBOSE_FAILURES */
2088 }
2089 if( T != RSB_BLAS_INVALID_MATRIX && rsb__BLAS_Xusds(T) != RSB_BLAS_NO_ERROR )
2090 errval = RSB_ERR_INTERNAL_ERROR;
2091 T = RSB_BLAS_INVALID_MATRIX;
2092 RSB_MTX_FREE(cmatrix);
2093
2094 RSB_CONDITIONAL_FREE(X);
2095 RSB_CONDITIONAL_FREE(B);
2096 RSB_CONDITIONAL_FREE(D);
2097 RSB_CONDITIONAL_FREE(IA);
2098 RSB_CONDITIONAL_FREE(JA);
2099 RSB_CONDITIONAL_FREE(VA);
2100 if(to.wqt!=RSB_BOOL_TRUE)
2101 RSB_INFO("%s%7d: type:%c sym:%s incX:%d incB:%d dim:%10d transT:%c alpha:%+2d beta:%+2d diag:%c subms:%5d nz:%d",btps,passed,(char)typecode,RSB_BLAS_MT_STR(stype),incX,incB,dim,tc,alphaa[alphai],betaa[betai],RSB_BLAS_DIAG_CHAR(diaga[diagi]),submatrices,rnz);
2102 errvalf|=errval;
2103
2104 if(errval == RSB_ERR_NO_ERROR)
2105 {
2106 if(to.wqt!=RSB_BOOL_TRUE)
2107 RSB_INFO(" is ok\n");
2108 ++passed;
2109 }
2110 else
2111 {
2112 if(to.wqt!=RSB_BOOL_TRUE)
2113 RSB_INFO(" is not ok\n");
2114 ++failed;
2115 RSB_INFO("Terminating testing due to errors.\n");
2116 goto done;
2117 }
2118
2119 if(RSB_SHALL_QUIT)
2120 {
2121 RSB_INFO("Terminating testing earlier due to interactive user request: test took %lf s, max allowed was %lf.\n",tt,to.mtt);
2122 goto done;
2123 }
2124
2125 #if RSB_TESTER_ALLOW_TIMEOUT
2126 if(to.mtt != RSB_TIME_ZERO && (tt = rsb_time()-tt0)>to.mtt)
2127 {
2128 RSB_INFO("Terminating testing earlier due to user timeout request: test took %lf s, max allowed was %lf.\n",tt,to.mtt);
2129 goto done;
2130 }
2131 #endif /* RSB_TESTER_ALLOW_TIMEOUT */
2132 #if RSB_ALLOW_INTERNAL_GETENVS
2133 if( maxtc != 0 && passed + failed >= maxtc )
2134 {
2135 RSB_INFO("Terminating testing earlier due to user limit request to %d tests.\n",maxtc);
2136 goto done;
2137 }
2138 #endif /* RSB_ALLOW_INTERNAL_GETENVS */
2139 if(to.tur==RSB_BOOL_TRUE && instantiated_some_recursive==1 && failed==0)
2140 {
2141 RSB_INFO("ALL TESTS PASSED SO FAR, AND ALSO INSTANTIATED ONE \"RECURSIVE\" MATRIX... THIS IS ENOUGH\n");
2142 errval = RSB_ERR_NO_ERROR;
2143 goto done;
2144 }
2145 }
2146 done:
2147 if(to.rrm==RSB_BOOL_TRUE && instantiated_some_recursive==0 && failed==0)
2148 {
2149 RSB_INFO("STRANGE: TESTS PASSED, BUT DID NOT INSTANTIATE ANY \"RECURSIVE\" MATRIX... RAISING AN ERROR FOR THIS\n");
2150 errvalf |= RSB_ERR_INTERNAL_ERROR;
2151 rsb__do_perror(NULL,RSB_ERR_INTERNAL_ERROR);
2152 }
2153 RSB_INFO(" PASSED:%d\n FAILED:%d\n",passed,failed);
2154 RSB_INFO("ADVANCED SPARSE BLAS TEST: END (SUCCESS)\n");
2155 // return RSB_ERR_NO_ERROR;
2156 errval=errvalf;
2157 ret:
2158 RSB_DO_ERR_RETURN(errval)
2159 }
2160
rsb_blas_runtime_limits_tester(void)2161 rsb_err_t rsb_blas_runtime_limits_tester(void)
2162 {
2163 /**
2164 \ingroup gr_internals
2165
2166 TODO: INCOMPLETE
2167 */
2168 rsb_err_t errval = RSB_ERR_NO_ERROR;
2169 rsb_nnz_idx_t maxcoo = RSB_MAX_MATRIX_DIM;
2170 rsb_coo_idx_t * IA=NULL;
2171 size_t maxcoo_bytes=0;
2172 const size_t minmem=1;
2173 size_t free_mem=0;
2174 size_t tot_mem=0;
2175 rsb_nnz_idx_t i;
2176 rsb_nnz_idx_t fel;
2177
2178 // FIXME: should fix the code revolving around the following:
2179 // RSB_STDOUT("%u>=%u : %d\n", rsb__nearest_power_of_two(maxcoo), maxcoo , rsb__nearest_power_of_two(maxcoo)>= maxcoo );
2180 RSB_INFO("Beginning large binary search test.\n");
2181 maxcoo_bytes=((size_t)maxcoo)*sizeof(rsb_coo_idx_t);
2182
2183 free_mem = rsb__sys_free_system_memory();
2184 tot_mem = rsb__sys_total_system_memory();
2185
2186 RSB_INFO("Detected %zu bytes of memory, comprehensive of %zu of free memory.\n",tot_mem,free_mem);
2187
2188 if(tot_mem<minmem || free_mem<minmem)
2189 {
2190 RSB_INFO("Too little memory detected: seems like your system is not well supported or not standards compliant.\n");
2191 tot_mem=free_mem=1024*1024*16;
2192 maxcoo=tot_mem/sizeof(rsb_coo_idx_t);
2193 RSB_INFO("Will try setting a reasonably small value: %zu for detected free memory.\n",free_mem);
2194 //goto skip_max_coo_test;
2195 }
2196 RSB_INFO("On this system, maximal array of coordinates can have %zu elements and occupy %zu bytes.\n",((size_t)maxcoo),maxcoo_bytes);
2197
2198 if(tot_mem<maxcoo_bytes || RSB_MUL_OVERFLOW(maxcoo,sizeof(rsb_coo_idx_t),rsb_nnz_idx_t,size_t))
2199 {
2200 /* FIXME: overflow cases shall be handled better */
2201 maxcoo=(3*free_mem/4)/sizeof(rsb_coo_idx_t);
2202 maxcoo_bytes=sizeof(rsb_coo_idx_t)*maxcoo;
2203 RSB_INFO("Will perform the test using less memory (%zu MB) than on the maximal coordinate indices array (%zu) allows.\n",maxcoo_bytes/(1024*1024),maxcoo_bytes);
2204 }
2205
2206 if(tot_mem<maxcoo_bytes)
2207 {
2208 RSB_INFO("Skipping test: too little memory.\n");
2209 goto skip_max_coo_test;
2210 }
2211
2212 if(free_mem<maxcoo_bytes)
2213 {
2214 RSB_INFO("Detected %zd bytes of free memory, needed %zd\nlet's see if test succeed ..\n",free_mem,maxcoo_bytes);
2215 //RSB_STDOUT("detected %zd bytes of free memory, needed %zd.\n",free_mem,maxcoo_bytes);
2216 //RSB_STDOUT("detected %zd bytes of free memory, needed %zd.\n");
2217 //RSB_STDOUT("detected %zd bytes of free memory, needed %zd. skipping test.\n",free_mem,maxcoo_bytes);
2218 //goto skip_max_coo_test;
2219 }
2220 IA = rsb__calloc(sizeof(rsb_coo_idx_t)*maxcoo);
2221 for ( i=99; i>0 && !IA; --i )
2222 {
2223 const size_t nmaxcoo = ( maxcoo / 100 ) * i;
2224 IA = rsb__calloc((sizeof(rsb_coo_idx_t)*(nmaxcoo)));
2225 if(IA)
2226 {
2227 RSB_INFO("WARNING: Failed (c)allocating of %zd nnz (%zd bytes)\n",(size_t)maxcoo,maxcoo_bytes);
2228 maxcoo = nmaxcoo;
2229 maxcoo_bytes = sizeof(rsb_coo_idx_t)*maxcoo;
2230 RSB_INFO("But made it with %zd nnz (%zd bytes, %d%% of that). Are you running in a containerized environment?\n",(size_t)maxcoo,maxcoo_bytes,i);
2231 }
2232 }
2233 if(!IA)
2234 {
2235 RSB_INFO("Failed (c)allocating of %zd nnz (%zd bytes)\n",(size_t)maxcoo,maxcoo_bytes);
2236 if(free_mem>maxcoo_bytes)
2237 {
2238 errval = RSB_ERR_ENOMEM;
2239 goto err;
2240 }
2241 else
2242 goto skip_max_coo_test;
2243 }
2244 else
2245 {
2246 RSB_INFO("(c)allocated %zd nnz (%zd bytes)\n",(size_t)maxcoo,maxcoo_bytes);
2247 }
2248
2249 for(i=0;i<maxcoo;++i)IA[i]=i;
2250
2251 fel = rsb__seek_coo_idx_t(IA,maxcoo-1,maxcoo);
2252 if(fel == RSB_MARKER_NNZ_VALUE || IA[fel]!=fel)
2253 {
2254 RSB_INFO("Failed retrieving array last element!\n");
2255 errval = RSB_ERR_INTERNAL_ERROR;
2256 goto err;
2257 }
2258 else
2259 RSB_INFO("Succeeded retrieving array last element.\n");
2260
2261 RSB_INFO("Successfully performed large binary search test.\n");
2262 goto done;
2263 skip_max_coo_test:
2264 RSB_INFO("Skipping large binary search test.\n");
2265 done:
2266 err:
2267 rsb__do_perror(NULL,errval);
2268 RSB_CONDITIONAL_FREE(IA);
2269 RSB_DO_ERR_RETURN(errval)
2270 }
2271
2272 /* @endcond */
2273