1 /* spmatrix/test_source.c
2 *
3 * Copyright (C) 2018, 2019, 2020 Patrick Alken
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 /*
21 test_random()
22 Create a random sparse matrix with approximately
23 M*N*density non-zero entries in [lower,upper]
24
25 Inputs: M - number of rows
26 N - number of columns
27 density - sparse density \in [0,1]
28 0 = no non-zero entries
29 1 = all m*n entries are filled
30 lower - lower bound on entries
31 upper - upper bound on entries
32 r - random number generator
33
34 Return: pointer to sparse matrix in triplet format (must be freed by caller)
35
36 Notes:
37 1) non-zero matrix entries are uniformly distributed in [0,1]
38 */
39
TYPE(gsl_spmatrix)40 static TYPE (gsl_spmatrix) *
41 FUNCTION (test, random)(const size_t M, const size_t N, const double density,
42 const double lower, const double upper, const gsl_rng *r)
43 {
44 size_t nnzwanted = (size_t) floor(M * N * GSL_MIN(density, 1.0));
45 TYPE (gsl_spmatrix) * m = FUNCTION (gsl_spmatrix, alloc_nzmax) (M, N, nnzwanted, GSL_SPMATRIX_COO);
46
47 while (FUNCTION (gsl_spmatrix, nnz) (m) < nnzwanted)
48 {
49 /* generate a random row and column */
50 size_t i = gsl_rng_uniform(r) * M;
51 size_t j = gsl_rng_uniform(r) * N;
52
53 /* generate random m_{ij} and add it */
54 BASE x = (BASE) (gsl_rng_uniform(r) * (upper - lower) + lower);
55 FUNCTION (gsl_spmatrix, set) (m, i, j, x);
56 }
57
58 return m;
59 }
60
TYPE(gsl_spmatrix)61 static TYPE (gsl_spmatrix) *
62 FUNCTION (test, random_int)(const size_t M, const size_t N, const double density,
63 const double lower, const double upper, const gsl_rng *r)
64 {
65 size_t nnzwanted = (size_t) floor(M * N * GSL_MIN(density, 1.0));
66 TYPE (gsl_spmatrix) * m = FUNCTION (gsl_spmatrix, alloc_nzmax) (M, N, nnzwanted, GSL_SPMATRIX_COO);
67
68 while (FUNCTION (gsl_spmatrix, nnz) (m) < nnzwanted)
69 {
70 /* generate a random row and column */
71 size_t i = gsl_rng_uniform(r) * M;
72 size_t j = gsl_rng_uniform(r) * N;
73
74 /* generate random m_{ij} and add it */
75 int x = (int) (gsl_rng_uniform(r) * (upper - lower) + lower);
76 FUNCTION (gsl_spmatrix, set) (m, i, j, (BASE) x);
77 }
78
79 return m;
80 }
81
82 static void
FUNCTION(test,random_dense)83 FUNCTION (test, random_dense)(TYPE (gsl_matrix) * m, const double lower,
84 const double upper, const gsl_rng *r)
85 {
86 const size_t M = m->size1;
87 const size_t N = m->size2;
88 size_t i, j;
89
90 for (i = 0; i < M; ++i)
91 {
92 for (j = 0; j < N; ++j)
93 {
94 double mij = gsl_rng_uniform(r) * (upper - lower) + lower;
95 FUNCTION (gsl_matrix, set) (m, i, j, (BASE) mij);
96 }
97 }
98 }
99
100 static void
FUNCTION(test,random_vector)101 FUNCTION (test, random_vector)(TYPE (gsl_vector) * x, const double lower,
102 const double upper, const gsl_rng *r)
103 {
104 size_t i;
105
106 for (i = 0; i < x->size; ++i)
107 {
108 double xi = gsl_rng_uniform(r) * (upper - lower) + lower;
109 FUNCTION (gsl_vector, set) (x, i, (BASE) xi);
110 }
111 }
112
113 static void
FUNCTION(test,alloc)114 FUNCTION (test, alloc) (const size_t M, const size_t N, const int sptype)
115 {
116 TYPE (gsl_spmatrix) * m = FUNCTION (gsl_spmatrix, alloc_nzmax) (M, N, 10, sptype);
117 size_t i, j;
118
119 gsl_test (m->data == 0, NAME (gsl_spmatrix) "_alloc(%s) returns valid data pointer", FUNCTION (gsl_spmatrix, type) (m));
120 gsl_test (m->i == 0, NAME (gsl_spmatrix) "_alloc(%s) returns valid i pointer", FUNCTION (gsl_spmatrix, type) (m));
121 gsl_test (m->p == 0, NAME (gsl_spmatrix) "_alloc(%s) returns valid p pointer", FUNCTION (gsl_spmatrix, type) (m));
122 gsl_test (m->size1 != M, NAME (gsl_spmatrix) "_alloc(%s) returns valid size1", FUNCTION (gsl_spmatrix, type) (m));
123 gsl_test (m->size2 != N, NAME (gsl_spmatrix) "_alloc(%s) returns valid size2", FUNCTION (gsl_spmatrix, type) (m));
124 gsl_test (m->nz != 0, NAME (gsl_spmatrix) "_alloc(%s) returns valid nz", FUNCTION (gsl_spmatrix, type) (m));
125 gsl_test (m->sptype != sptype, NAME (gsl_spmatrix) "_alloc(%s) returns valid sptype", FUNCTION (gsl_spmatrix, type) (m));
126
127 /* test that a newly allocated matrix is initialized to zero */
128 status = 0;
129 for (i = 0; i < M; ++i)
130 {
131 for (j = 0; j < N; ++j)
132 {
133 BASE mij = FUNCTION (gsl_spmatrix, get) (m, i, j);
134
135 if (mij != (BASE) 0)
136 status = 1;
137 }
138 }
139
140 gsl_test (status, NAME (gsl_spmatrix) "_alloc(%s) initializes matrix to zero", FUNCTION (gsl_spmatrix, type) (m));
141
142 FUNCTION (gsl_spmatrix, free) (m);
143 }
144
145 static void
FUNCTION(test,realloc)146 FUNCTION (test, realloc) (const size_t M, const size_t N)
147 {
148 TYPE (gsl_spmatrix) * m = FUNCTION (gsl_spmatrix, alloc_nzmax) (M, N, 10, GSL_SPMATRIX_COO);
149 const BASE x = (BASE) 5;
150 const size_t nnz_total = 30;
151 size_t i, j, nnz;
152
153 /* add enough elements to m to trigger a realloc call, then check the result */
154
155 nnz = 0;
156 i = 0;
157 j = 0;
158 while (nnz++ < nnz_total)
159 {
160 if (i >= m->size1)
161 {
162 i = 0;
163 if (++j >= m->size2)
164 break;
165 }
166
167 FUNCTION (gsl_spmatrix, set) (m, i++, j, x);
168 }
169
170 status = 0;
171 nnz = 0;
172 i = 0;
173 j = 0;
174 while (nnz++ < nnz_total)
175 {
176 if (i >= m->size1)
177 {
178 i = 0;
179 if (++j >= m->size2)
180 break;
181 }
182
183 if (FUNCTION (gsl_spmatrix, get) (m, i++, j) != x)
184 status = 1;
185 }
186
187 gsl_test (status, NAME (gsl_spmatrix) "_realloc[%zu,%zu]", M, N);
188
189 FUNCTION (gsl_spmatrix, free) (m);
190 }
191
192 static void
FUNCTION(test,getset)193 FUNCTION (test, getset) (const size_t M, const size_t N, const int sptype)
194 {
195 TYPE (gsl_spmatrix) * A = FUNCTION (gsl_spmatrix, alloc_nzmax) (M, N, M * N, GSL_SPMATRIX_COO);
196 TYPE (gsl_spmatrix) * B = NULL;
197 size_t i, j, idx;
198 size_t k = 0;
199
200 for (i = 0; i < M; ++i)
201 {
202 for (j = 0; j < N; ++j)
203 {
204 k++;
205 FUNCTION (gsl_spmatrix, set) (A, i, j, (BASE) k);
206 }
207 }
208
209 status = 0;
210 k = 0;
211 idx = 0;
212 for (i = 0; i < M; ++i)
213 {
214 for (j = 0; j < N; ++j)
215 {
216 k++;
217 if ((A->i[idx] != (int) i) || (A->p[idx] != (int) j) ||
218 (A->data[idx] != (BASE) k))
219 status = 1;
220 idx++;
221 }
222 }
223
224 gsl_test (status, NAME (gsl_spmatrix) "_set[%zu,%zu] writes into arrays", M, N);
225
226 B = FUNCTION (gsl_spmatrix, compress) (A, sptype);
227
228 status = 0;
229 k = 0;
230 for (i = 0; i < M; ++i)
231 {
232 for (j = 0; j < N; ++j)
233 {
234 k++;
235 if (FUNCTION (gsl_spmatrix, get) (B, i, j) != (BASE) k)
236 status = 1;
237 }
238 }
239
240 gsl_test (status, NAME (gsl_spmatrix) "_get[%zu,%zu](%s) reads from arrays",
241 M, N, FUNCTION (gsl_spmatrix, type) (B));
242
243 status = FUNCTION (gsl_spmatrix, nnz) (B) != M * N;
244 gsl_test (status, NAME (gsl_spmatrix) "_nnz[%zu,%zu](%s) %zu/%zu",
245 M, N, FUNCTION (gsl_spmatrix, type) (B),
246 FUNCTION (gsl_spmatrix, nnz) (B), M * N);
247
248 /* test gsl_spmatrix_ptr */
249 status = 0;
250 for (i = 0; i < M; ++i)
251 {
252 for (j = 0; j < N; ++j)
253 {
254 BASE mij = FUNCTION (gsl_spmatrix, get) (B, i, j);
255 BASE * ptr = FUNCTION (gsl_spmatrix, ptr) (B, i, j);
256
257 if (ptr != NULL)
258 {
259 *ptr += (ATOMIC) 2;
260 if (FUNCTION (gsl_spmatrix, get) (B, i, j) != (BASE) (mij + (ATOMIC) 2))
261 status = 1;
262 }
263 }
264 }
265
266 gsl_test (status, NAME (gsl_spmatrix) "_ptr[%zu,%zu](%s)",
267 M, N, FUNCTION (gsl_spmatrix, type) (B));
268
269 /* test set_zero */
270 FUNCTION (gsl_spmatrix, set_zero) (B);
271 status = FUNCTION (gsl_spmatrix, nnz) (B) != 0;
272 gsl_test (status, NAME (gsl_spmatrix) "_set_zero[%zu,%zu](%s) nnz=%zu",
273 M, N, FUNCTION (gsl_spmatrix, type) (B), FUNCTION (gsl_spmatrix, nnz) (B));
274
275 if (B->sptype == GSL_SPMATRIX_COO)
276 {
277 const size_t min = GSL_MIN(M, N);
278 const size_t expected_nnz = min;
279 size_t nnz;
280
281 /* test setting duplicate values */
282 status = 0;
283 k = 0;
284 for (i = 0; i < min; ++i)
285 {
286 for (j = 0; j < 5; ++j)
287 {
288 BASE x = (BASE) ++k;
289 BASE y;
290
291 FUNCTION (gsl_spmatrix, set) (B, i, i, x);
292 y = FUNCTION (gsl_spmatrix, get) (B, i, i);
293 if (x != y)
294 status = 1;
295 }
296 }
297
298 gsl_test (status, NAME (gsl_spmatrix) " duplicates[%zu,%zu]", M, N);
299
300 nnz = FUNCTION (gsl_spmatrix, nnz) (B);
301 status = nnz != expected_nnz;
302 gsl_test (status, NAME (gsl_spmatrix) " duplicates nnz[%zu,%zu] %zu/%zu",
303 M, N, nnz, expected_nnz);
304
305 /* test setting an element to 0 */
306 FUNCTION (gsl_spmatrix, set) (B, 0, 0, (BASE) 1);
307 FUNCTION (gsl_spmatrix, set) (B, 0, 0, (BASE) 0);
308 status = FUNCTION (gsl_spmatrix, get) (B, 0, 0) != (BASE) 0;
309 gsl_test (status, NAME (gsl_spmatrix) " zero element[%zu,%zu] %f",
310 M, N, FUNCTION (gsl_spmatrix, get) (B, 0, 0));
311 }
312
313 FUNCTION (gsl_spmatrix, free) (A);
314 FUNCTION (gsl_spmatrix, free) (B);
315 }
316
317 static void
FUNCTION(test,memcpy)318 FUNCTION (test, memcpy) (const size_t M, const size_t N, const int sptype,
319 const double density, gsl_rng * r)
320 {
321 /* test memcpy */
322 {
323 TYPE (gsl_spmatrix) * m = FUNCTION (test, random) (M, N, density, 1.0, 20.0, r);
324 TYPE (gsl_spmatrix) * A = FUNCTION (gsl_spmatrix, compress) (m, sptype);
325 TYPE (gsl_spmatrix) * B = FUNCTION (gsl_spmatrix, alloc_nzmax) (M, N, m->nz, sptype);
326 TYPE (gsl_spmatrix) * C = FUNCTION (gsl_spmatrix, alloc_nzmax) (M, N, 1, sptype);
327
328 FUNCTION (gsl_spmatrix, memcpy) (B, A);
329 status = FUNCTION (gsl_spmatrix, equal) (A, B) != 1;
330 gsl_test (status, NAME (gsl_spmatrix) "_memcpy[%zu,%zu](%s) equality",
331 M, N, FUNCTION (gsl_spmatrix, type) (A));
332
333 /* this will force a realloc call */
334 FUNCTION (gsl_spmatrix, memcpy) (C, A);
335 status = FUNCTION (gsl_spmatrix, equal) (A, C) != 1;
336 gsl_test (status, NAME (gsl_spmatrix) "_memcpy[%zu,%zu](%s) realloc equality",
337 M, N, FUNCTION (gsl_spmatrix, type) (A));
338
339 FUNCTION (gsl_spmatrix, free) (m);
340 FUNCTION (gsl_spmatrix, free) (A);
341 FUNCTION (gsl_spmatrix, free) (B);
342 FUNCTION (gsl_spmatrix, free) (C);
343 }
344
345 /* test transpose memcpy */
346 {
347 TYPE (gsl_spmatrix) * m = FUNCTION (test, random) (M, N, density, 1.0, 20.0, r);
348 TYPE (gsl_spmatrix) * A = FUNCTION (gsl_spmatrix, compress) (m, sptype);
349 TYPE (gsl_spmatrix) * AT = FUNCTION (gsl_spmatrix, alloc_nzmax) (N, M, m->nz, sptype);
350 TYPE (gsl_spmatrix) * BT = FUNCTION (gsl_spmatrix, alloc_nzmax) (N, M, 1, sptype);
351 size_t i, j;
352
353 FUNCTION (gsl_spmatrix, transpose_memcpy) (AT, A);
354 FUNCTION (gsl_spmatrix, transpose_memcpy) (BT, A); /* force realloc call */
355
356 status = 0;
357 for (i = 0; i < M; ++i)
358 {
359 for (j = 0; j < N; ++j)
360 {
361 BASE Aij = FUNCTION (gsl_spmatrix, get) (A, i, j);
362 BASE ATji = FUNCTION (gsl_spmatrix, get) (AT, j, i);
363 BASE BTji = FUNCTION (gsl_spmatrix, get) (BT, j, i);
364
365 if (Aij != ATji)
366 status = 1;
367 if (Aij != BTji)
368 status = 2;
369 }
370 }
371
372 gsl_test(status == 1, NAME (gsl_spmatrix) "_transpose_memcpy[%zu,%zu](%s) AT",
373 M, N, FUNCTION (gsl_spmatrix, type) (A));
374 gsl_test(status == 2, NAME (gsl_spmatrix) "_transpose_memcpy[%zu,%zu](%s) BT",
375 M, N, FUNCTION (gsl_spmatrix, type) (A));
376
377 FUNCTION (gsl_spmatrix, free) (m);
378 FUNCTION (gsl_spmatrix, free) (A);
379 FUNCTION (gsl_spmatrix, free) (AT);
380 FUNCTION (gsl_spmatrix, free) (BT);
381 }
382 }
383
384 static void
FUNCTION(test,transpose)385 FUNCTION (test, transpose) (const size_t M, const size_t N, const int sptype,
386 const double density, gsl_rng * r)
387 {
388 TYPE (gsl_spmatrix) * A = FUNCTION (test, random) (M, N, density, 1.0, 20.0, r);
389 TYPE (gsl_spmatrix) * B = FUNCTION (gsl_spmatrix, compress) (A, sptype);
390 size_t i;
391
392 FUNCTION (gsl_spmatrix, transpose) (B);
393
394 status = 0;
395 for (i = 0; i < A->nz; ++i)
396 {
397 BASE Aij = A->data[i];
398 BASE Bji = FUNCTION (gsl_spmatrix, get) (B, A->p[i], A->i[i]);
399
400 if (Aij != Bji)
401 status = 1;
402 }
403
404 gsl_test(status, NAME (gsl_spmatrix) "_transpose[%zu,%zu](%s)",
405 M, N, FUNCTION (gsl_spmatrix, type) (B));
406
407 FUNCTION (gsl_spmatrix, free) (A);
408 FUNCTION (gsl_spmatrix, free) (B);
409 }
410
411 static void
FUNCTION(test,scale)412 FUNCTION (test, scale) (const size_t M, const size_t N, const int sptype,
413 const double density, gsl_rng * r)
414 {
415 TYPE (gsl_spmatrix) * A = FUNCTION (test, random) (M, N, density, 1.0, 20.0, r);
416 TYPE (gsl_spmatrix) * B = FUNCTION (gsl_spmatrix, compress) (A, sptype);
417 TYPE (gsl_vector) * x;
418 size_t i, j;
419
420 FUNCTION (gsl_spmatrix, scale) (B, 2.0);
421
422 status = 0;
423 for (i = 0; i < M; ++i)
424 {
425 for (j = 0; j < N; ++j)
426 {
427 BASE aij = FUNCTION (gsl_spmatrix, get) (A, i, j);
428 BASE bij = FUNCTION (gsl_spmatrix, get) (B, i, j);
429 if (bij != (ATOMIC) (2*aij))
430 status = 1;
431 }
432 }
433
434 gsl_test (status, NAME (gsl_spmatrix) "_scale[%zu,%zu](%s)",
435 M, N, FUNCTION (gsl_spmatrix, type) (B));
436
437 /* reset B = A */
438 FUNCTION (gsl_spmatrix, free) (B);
439 B = FUNCTION (gsl_spmatrix, compress) (A, sptype);
440
441 /* test column scaling */
442 x = FUNCTION (gsl_vector, alloc) (A->size2);
443 FUNCTION (test, random_vector) (x, 1.0, 20.0, r);
444
445 FUNCTION (gsl_spmatrix, scale_columns) (B, x);
446
447 status = 0;
448 for (i = 0; i < A->nz; ++i)
449 {
450 BASE aij = A->data[i];
451 BASE bij = FUNCTION (gsl_spmatrix, get) (B, A->i[i], A->p[i]);
452 BASE xj = FUNCTION (gsl_vector, get) (x, A->p[i]);
453
454 if (bij != (ATOMIC) (xj * aij))
455 status = 1;
456 }
457
458 gsl_test (status, NAME (gsl_spmatrix) "_scale_columns[%zu,%zu](%s)",
459 M, N, FUNCTION (gsl_spmatrix, type) (B));
460
461 FUNCTION (gsl_vector, free) (x);
462
463 /* reset B = A */
464 FUNCTION (gsl_spmatrix, free) (B);
465 B = FUNCTION (gsl_spmatrix, compress) (A, sptype);
466
467 /* test row scaling */
468 x = FUNCTION (gsl_vector, alloc) (A->size1);
469 FUNCTION (test, random_vector) (x, 1.0, 20.0, r);
470
471 FUNCTION (gsl_spmatrix, scale_rows) (B, x);
472
473 status = 0;
474 for (i = 0; i < A->nz; ++i)
475 {
476 BASE aij = A->data[i];
477 BASE bij = FUNCTION (gsl_spmatrix, get) (B, A->i[i], A->p[i]);
478 BASE xi = FUNCTION (gsl_vector, get) (x, A->i[i]);
479
480 if (bij != (ATOMIC) (xi * aij))
481 status = 1;
482 }
483
484 gsl_test (status, NAME (gsl_spmatrix) "_scale_rows[%zu,%zu](%s)",
485 M, N, FUNCTION (gsl_spmatrix, type) (B));
486
487 FUNCTION (gsl_vector, free) (x);
488 FUNCTION (gsl_spmatrix, free) (A);
489 FUNCTION (gsl_spmatrix, free) (B);
490 }
491
492 static void
FUNCTION(test,add)493 FUNCTION (test, add) (const size_t M, const size_t N, const int sptype,
494 const double density, gsl_rng * r)
495 {
496 TYPE (gsl_spmatrix) * m1 = FUNCTION (test, random) (M, N, density, 1.0, 20.0, r);
497 TYPE (gsl_spmatrix) * m2 = FUNCTION (test, random) (M, N, density, 1.0, 20.0, r);
498 TYPE (gsl_spmatrix) * A = FUNCTION (gsl_spmatrix, compress) (m1, sptype);
499 TYPE (gsl_spmatrix) * B = FUNCTION (gsl_spmatrix, compress) (m2, sptype);
500 TYPE (gsl_spmatrix) * C = FUNCTION (gsl_spmatrix, alloc_nzmax) (M, N, 1, sptype);
501 size_t i, j;
502
503 FUNCTION (gsl_spmatrix, add) (C, A, B);
504
505 status = 0;
506 for (i = 0; i < M; ++i)
507 {
508 for (j = 0; j < N; ++j)
509 {
510 BASE aij = FUNCTION (gsl_spmatrix, get) (A, i, j);
511 BASE bij = FUNCTION (gsl_spmatrix, get) (B, i, j);
512 BASE cij = FUNCTION (gsl_spmatrix, get) (C, i, j);
513
514 if (aij + bij != cij)
515 status = 1;
516 }
517 }
518
519 gsl_test (status, NAME (gsl_spmatrix) "_add[%zu,%zu](%s)",
520 M, N, FUNCTION (gsl_spmatrix, type) (A));
521
522 /* test again with C = 2*A */
523 FUNCTION (gsl_spmatrix, add) (C, A, A);
524
525 status = 0;
526 for (i = 0; i < M; ++i)
527 {
528 for (j = 0; j < N; ++j)
529 {
530 BASE aij = FUNCTION (gsl_spmatrix, get) (A, i, j);
531 BASE cij = FUNCTION (gsl_spmatrix, get) (C, i, j);
532
533 if (aij + aij != cij)
534 status = 1;
535 }
536 }
537
538 gsl_test (status, NAME (gsl_spmatrix) "_add[%zu,%zu](%s) duplicate",
539 M, N, FUNCTION (gsl_spmatrix, type) (A));
540
541 FUNCTION (gsl_spmatrix, free) (A);
542 FUNCTION (gsl_spmatrix, free) (B);
543 FUNCTION (gsl_spmatrix, free) (C);
544 FUNCTION (gsl_spmatrix, free) (m1);
545 FUNCTION (gsl_spmatrix, free) (m2);
546 }
547
548 static void
FUNCTION(test,dense_add)549 FUNCTION (test, dense_add) (const size_t M, const size_t N, const int sptype,
550 const double density, gsl_rng * r)
551 {
552 TYPE (gsl_spmatrix) * m = FUNCTION (test, random) (M, N, density, 1.0, 20.0, r);
553 TYPE (gsl_spmatrix) * B = FUNCTION (gsl_spmatrix, compress) (m, sptype);
554 TYPE (gsl_matrix) * A = FUNCTION (gsl_matrix, alloc) (M, N);
555 TYPE (gsl_matrix) * A_copy = FUNCTION (gsl_matrix, alloc) (M, N);
556 size_t i, j;
557
558 FUNCTION (test, random_dense) (A, 1.0, 20.0, r);
559 FUNCTION (gsl_matrix, memcpy) (A_copy, A);
560
561 FUNCTION (gsl_spmatrix, dense_add) (A, B);
562
563 status = 0;
564 for (i = 0; i < M; ++i)
565 {
566 for (j = 0; j < N; ++j)
567 {
568 BASE aij = FUNCTION (gsl_matrix, get) (A_copy, i, j);
569 BASE bij = FUNCTION (gsl_spmatrix, get) (B, i, j);
570 BASE cij = FUNCTION (gsl_matrix, get) (A, i, j);
571
572 if (aij + bij != cij)
573 status = 1;
574 }
575 }
576
577 gsl_test (status, NAME (gsl_spmatrix) "_dense_add[%zu,%zu](%s)",
578 M, N, FUNCTION (gsl_spmatrix, type) (B));
579
580 FUNCTION (gsl_matrix, free) (A);
581 FUNCTION (gsl_matrix, free) (A_copy);
582 FUNCTION (gsl_spmatrix, free) (B);
583 FUNCTION (gsl_spmatrix, free) (m);
584 }
585
586 static void
FUNCTION(test,dense_sub)587 FUNCTION (test, dense_sub) (const size_t M, const size_t N, const int sptype,
588 const double density, gsl_rng * r)
589 {
590 TYPE (gsl_spmatrix) * m = FUNCTION (test, random) (M, N, density, 1.0, 20.0, r);
591 TYPE (gsl_spmatrix) * B = FUNCTION (gsl_spmatrix, compress) (m, sptype);
592 TYPE (gsl_matrix) * A = FUNCTION (gsl_matrix, alloc) (M, N);
593 TYPE (gsl_matrix) * A_copy = FUNCTION (gsl_matrix, alloc) (M, N);
594 size_t i, j;
595
596 FUNCTION (test, random_dense) (A, 1.0, 20.0, r);
597 FUNCTION (gsl_matrix, memcpy) (A_copy, A);
598
599 FUNCTION (gsl_spmatrix, dense_sub) (A, B);
600
601 status = 0;
602 for (i = 0; i < M; ++i)
603 {
604 for (j = 0; j < N; ++j)
605 {
606 BASE aij = FUNCTION (gsl_matrix, get) (A_copy, i, j);
607 BASE bij = FUNCTION (gsl_spmatrix, get) (B, i, j);
608 BASE cij = FUNCTION (gsl_matrix, get) (A, i, j);
609
610 if ((BASE) (aij - bij) != cij)
611 status = 1;
612 }
613 }
614
615 gsl_test (status, NAME (gsl_spmatrix) "_dense_sub[%zu,%zu](%s)",
616 M, N, FUNCTION (gsl_spmatrix, type) (B));
617
618 FUNCTION (gsl_matrix, free) (A);
619 FUNCTION (gsl_matrix, free) (A_copy);
620 FUNCTION (gsl_spmatrix, free) (B);
621 FUNCTION (gsl_spmatrix, free) (m);
622 }
623
624 static void
FUNCTION(test,convert)625 FUNCTION (test, convert) (const size_t M, const size_t N, const int sptype,
626 const double density, gsl_rng * r)
627 {
628 TYPE (gsl_spmatrix) * A = FUNCTION (test, random) (M, N, density, 1.0, 20.0, r);
629 TYPE (gsl_spmatrix) * B = FUNCTION (gsl_spmatrix, compress) (A, sptype);
630 TYPE (gsl_spmatrix) * C = FUNCTION (gsl_spmatrix, alloc_nzmax) (M, N, 1, GSL_SPMATRIX_COO);
631 TYPE (gsl_matrix) * D = FUNCTION (gsl_matrix, alloc) (M, N);
632 size_t i, j;
633
634 /* convert D := B */
635 FUNCTION (gsl_spmatrix, sp2d) (D, B);
636
637 status = 0;
638 for (i = 0; i < M; ++i)
639 {
640 for (j = 0; j < N; ++j)
641 {
642 BASE bij = FUNCTION (gsl_spmatrix, get) (B, i, j);
643 BASE dij = FUNCTION (gsl_matrix, get) (D, i, j);
644
645 if (bij != dij)
646 status = 1;
647 }
648 }
649
650 gsl_test (status, NAME (gsl_spmatrix) "_sp2d[%zu,%zu](%s)",
651 M, N, FUNCTION (gsl_spmatrix, type) (B));
652
653 /* convert C := D */
654 FUNCTION (gsl_spmatrix, d2sp) (C, D);
655
656 status = FUNCTION (gsl_spmatrix, equal) (A, C) != 1;
657 gsl_test (status, NAME (gsl_spmatrix) "_d2sp[%zu,%zu](%s)",
658 M, N, FUNCTION (gsl_spmatrix, type) (B));
659
660 FUNCTION (gsl_spmatrix, free) (A);
661 FUNCTION (gsl_spmatrix, free) (B);
662 FUNCTION (gsl_spmatrix, free) (C);
663 FUNCTION (gsl_matrix, free) (D);
664 }
665
666 static void
FUNCTION(test,minmax)667 FUNCTION (test, minmax) (const size_t M, const size_t N, const int sptype,
668 const double density, gsl_rng * r)
669 {
670 TYPE (gsl_spmatrix) * A = FUNCTION (test, random) (M, N, density, 5.0, 20.0, r);
671 TYPE (gsl_spmatrix) * B;
672 BASE min, max;
673 size_t imin, jmin;
674
675 FUNCTION (gsl_spmatrix, set) (A, 4, 3, (BASE) 12);
676 FUNCTION (gsl_spmatrix, set) (A, 3, 5, (BASE) 5);
677 FUNCTION (gsl_spmatrix, set) (A, 1, 1, (BASE) 10);
678 FUNCTION (gsl_spmatrix, set) (A, 7, 3, (BASE) 2);
679 FUNCTION (gsl_spmatrix, set) (A, 1, 0, (BASE) 30);
680
681 B = FUNCTION (gsl_spmatrix, compress) (A, sptype);
682
683 FUNCTION (gsl_spmatrix, minmax) (B, &min, &max);
684
685 status = min != (BASE) 2;
686 gsl_test (status, NAME (gsl_spmatrix) "_minmax[%zu,%zu](%s) minimum",
687 M, N, FUNCTION (gsl_spmatrix, type) (B));
688
689 status = max != (BASE) 30;
690 gsl_test (status, NAME (gsl_spmatrix) "_minmax[%zu,%zu](%s) maximum",
691 M, N, FUNCTION (gsl_spmatrix, type) (B));
692
693 FUNCTION (gsl_spmatrix, min_index) (B, &imin, &jmin);
694
695 status = imin != 7;
696 gsl_test (status, NAME (gsl_spmatrix) "_min_index[%zu,%zu](%s) imin",
697 M, N, FUNCTION (gsl_spmatrix, type) (B));
698
699 status = jmin != 3;
700 gsl_test (status, NAME (gsl_spmatrix) "_min_index[%zu,%zu](%s) jmin",
701 M, N, FUNCTION (gsl_spmatrix, type) (B));
702
703 FUNCTION (gsl_spmatrix, free) (A);
704 FUNCTION (gsl_spmatrix, free) (B);
705 }
706
707 #if !defined(UNSIGNED) && !defined(BASE_CHAR)
708
709 static void
FUNCTION(test,norm)710 FUNCTION (test, norm) (const size_t M, const size_t N, const int sptype)
711 {
712 TYPE (gsl_spmatrix) * A = FUNCTION (gsl_spmatrix, alloc) (M, N);
713 TYPE (gsl_spmatrix) * B;
714 ATOMIC norm, norm_expected;
715 size_t i, j, k = 0;
716
717 for (i = 0; i < M; i++)
718 {
719 for (j = 0; j < N; j++)
720 {
721 k++;
722 FUNCTION (gsl_spmatrix, set) (A, i, j, (BASE) k);
723 }
724 }
725
726 B = FUNCTION (gsl_spmatrix, compress) (A, sptype);
727
728 norm = FUNCTION (gsl_spmatrix, norm1) (B);
729 norm_expected = N*M*(M+1)/2;
730
731 status = norm != norm_expected;
732 gsl_test (status, NAME (gsl_spmatrix) "_norm1[%zu,%zu](%s)",
733 M, N, FUNCTION (gsl_spmatrix, type) (B));
734
735 FUNCTION (gsl_spmatrix, free) (A);
736 FUNCTION (gsl_spmatrix, free) (B);
737 }
738
739 #endif
740
741 static void
FUNCTION(test,io_ascii)742 FUNCTION (test, io_ascii) (const size_t M, const size_t N, const int sptype,
743 const double density, gsl_rng * r)
744 {
745 TYPE (gsl_spmatrix) * A = FUNCTION (test, random_int) (M, N, density, 1.0, 20.0, r);
746 TYPE (gsl_spmatrix) * B = FUNCTION (gsl_spmatrix, compress) (A, sptype);
747 TYPE (gsl_spmatrix) * C;
748 char filename[] = "test.dat";
749 FILE *f;
750
751 f = fopen (filename, "w");
752 FUNCTION (gsl_spmatrix, fprintf) (f, B, OUT_FORMAT);
753 fclose (f);
754
755 f = fopen (filename, "r");
756 C = FUNCTION (gsl_spmatrix, fscanf) (f);
757 fclose (f);
758
759 status = FUNCTION (gsl_spmatrix, equal) (A, C) != 1;
760 gsl_test (status, NAME (gsl_spmatrix) "_fscanf[%zu,%zu](%s)",
761 M, N, FUNCTION (gsl_spmatrix, type) (B));
762
763 unlink (filename);
764
765 FUNCTION (gsl_spmatrix, free) (A);
766 FUNCTION (gsl_spmatrix, free) (B);
767 FUNCTION (gsl_spmatrix, free) (C);
768 }
769
770 static void
FUNCTION(test,io_binary)771 FUNCTION (test, io_binary) (const size_t M, const size_t N, const int sptype,
772 const double density, gsl_rng * r)
773 {
774 TYPE (gsl_spmatrix) * A = FUNCTION (test, random) (M, N, density, 1.0, 20.0, r);
775 TYPE (gsl_spmatrix) * B = FUNCTION (gsl_spmatrix, compress) (A, sptype);
776 TYPE (gsl_spmatrix) * C = FUNCTION (gsl_spmatrix, alloc_nzmax) (M, N, A->nz, sptype);
777 char filename[] = "test.dat";
778 FILE *f;
779
780 f = fopen (filename, "wb");
781 FUNCTION (gsl_spmatrix, fwrite) (f, B);
782 fclose (f);
783
784 f = fopen (filename, "rb");
785 FUNCTION (gsl_spmatrix, fread) (f, C);
786 fclose (f);
787
788 status = FUNCTION (gsl_spmatrix, equal) (B, C) != 1;
789 gsl_test (status, NAME (gsl_spmatrix) "_fread[%zu,%zu](%s)",
790 M, N, FUNCTION (gsl_spmatrix, type) (B));
791
792 unlink (filename);
793
794 FUNCTION (gsl_spmatrix, free) (A);
795 FUNCTION (gsl_spmatrix, free) (B);
796 FUNCTION (gsl_spmatrix, free) (C);
797 }
798
799 static void
FUNCTION(test,all)800 FUNCTION (test, all) (const size_t M, const size_t N, const double density, gsl_rng * r)
801 {
802 FUNCTION (test, alloc) (M, N, GSL_SPMATRIX_COO);
803 FUNCTION (test, alloc) (M, N, GSL_SPMATRIX_CSC);
804 FUNCTION (test, alloc) (M, N, GSL_SPMATRIX_CSR);
805
806 FUNCTION (test, realloc) (M, N);
807
808 FUNCTION (test, getset) (M, N, GSL_SPMATRIX_COO);
809 FUNCTION (test, getset) (M, N, GSL_SPMATRIX_CSC);
810 FUNCTION (test, getset) (M, N, GSL_SPMATRIX_CSR);
811
812 FUNCTION (test, memcpy) (M, N, GSL_SPMATRIX_COO, density, r);
813 FUNCTION (test, memcpy) (M, N, GSL_SPMATRIX_CSC, density, r);
814 FUNCTION (test, memcpy) (M, N, GSL_SPMATRIX_CSR, density, r);
815
816 FUNCTION (test, transpose) (M, N, GSL_SPMATRIX_COO, density, r);
817 FUNCTION (test, transpose) (M, N, GSL_SPMATRIX_CSC, density, r);
818 FUNCTION (test, transpose) (M, N, GSL_SPMATRIX_CSR, density, r);
819
820 FUNCTION (test, scale) (M, N, GSL_SPMATRIX_COO, density, r);
821 FUNCTION (test, scale) (M, N, GSL_SPMATRIX_CSC, density, r);
822 FUNCTION (test, scale) (M, N, GSL_SPMATRIX_CSR, density, r);
823
824 FUNCTION (test, add) (M, N, GSL_SPMATRIX_CSC, density, r);
825 FUNCTION (test, add) (M, N, GSL_SPMATRIX_CSR, density, r);
826
827 FUNCTION (test, dense_add) (M, N, GSL_SPMATRIX_COO, density, r);
828 FUNCTION (test, dense_add) (M, N, GSL_SPMATRIX_CSC, density, r);
829 FUNCTION (test, dense_add) (M, N, GSL_SPMATRIX_CSR, density, r);
830
831 FUNCTION (test, dense_sub) (M, N, GSL_SPMATRIX_COO, density, r);
832 FUNCTION (test, dense_sub) (M, N, GSL_SPMATRIX_CSC, density, r);
833 FUNCTION (test, dense_sub) (M, N, GSL_SPMATRIX_CSR, density, r);
834
835 FUNCTION (test, convert) (M, N, GSL_SPMATRIX_COO, density, r);
836 FUNCTION (test, convert) (M, N, GSL_SPMATRIX_CSC, density, r);
837 FUNCTION (test, convert) (M, N, GSL_SPMATRIX_CSR, density, r);
838
839 FUNCTION (test, minmax) (M, N, GSL_SPMATRIX_COO, density, r);
840 FUNCTION (test, minmax) (M, N, GSL_SPMATRIX_CSC, density, r);
841 FUNCTION (test, minmax) (M, N, GSL_SPMATRIX_CSR, density, r);
842
843 #if !defined(UNSIGNED) && !defined(BASE_CHAR)
844 FUNCTION (test, norm) (M, N, GSL_SPMATRIX_COO);
845 FUNCTION (test, norm) (M, N, GSL_SPMATRIX_CSC);
846 FUNCTION (test, norm) (M, N, GSL_SPMATRIX_CSR);
847 #endif
848
849 FUNCTION (test, io_ascii) (M, N, GSL_SPMATRIX_COO, density, r);
850 FUNCTION (test, io_ascii) (M, N, GSL_SPMATRIX_CSC, density, r);
851 FUNCTION (test, io_ascii) (M, N, GSL_SPMATRIX_CSR, density, r);
852
853 FUNCTION (test, io_binary) (M, N, GSL_SPMATRIX_COO, density, r);
854 FUNCTION (test, io_binary) (M, N, GSL_SPMATRIX_CSC, density, r);
855 FUNCTION (test, io_binary) (M, N, GSL_SPMATRIX_CSR, density, r);
856 }
857