1 /*
2 * -----------------------------------------------------------------
3 * Programmer(s): David J. Gardner, Cody J. Balos @ LLNL
4 * Daniel R. Reynolds @ SMU
5 * -----------------------------------------------------------------
6 * SUNDIALS Copyright Start
7 * Copyright (c) 2002-2020, Lawrence Livermore National Security
8 * and Southern Methodist University.
9 * All rights reserved.
10 *
11 * See the top-level LICENSE and NOTICE files for details.
12 *
13 * SPDX-License-Identifier: BSD-3-Clause
14 * SUNDIALS Copyright End
15 * -----------------------------------------------------------------
16 * These test functions are designed to check a SUNMatrix module
17 * implementation.
18 * -----------------------------------------------------------------
19 */
20
21 #include <sundials/sundials_matrix.h>
22 #include <sundials/sundials_types.h>
23 #include <sundials/sundials_math.h>
24
25 #include <math.h> /* include isnan */
26 #include <stdio.h>
27 #include <stdlib.h>
28
29 #include "test_sunmatrix.h"
30
31 #if defined( SUNDIALS_HAVE_POSIX_TIMERS) && defined(_POSIX_TIMERS)
32 #include <time.h>
33 #include <unistd.h>
34 #endif
35
36 /* private functions */
37 static double get_time();
38
39 int print_time = 0;
40 int print_all_ranks = 0;
41
42 #define FMT "%s Time: %22.15e\n\n"
43 #define PRINT_TIME(format, time) if(print_time) printf(format, time)
44
45 /* ----------------------------------------------------------------------
46 * SUNMatGetID Test
47 * --------------------------------------------------------------------*/
Test_SUNMatGetID(SUNMatrix A,SUNMatrix_ID sunid,int myid)48 int Test_SUNMatGetID(SUNMatrix A, SUNMatrix_ID sunid, int myid)
49 {
50 double start_time, stop_time;
51 SUNMatrix_ID mysunid;
52
53 start_time = get_time();
54 mysunid = SUNMatGetID(A);
55 stop_time = get_time();
56
57 if (sunid != mysunid) {
58 TEST_STATUS(">>> FAILED test -- SUNMatGetID \n", myid);
59 PRINT_TIME(" SUNMatGetID Time: %22.15e \n \n", stop_time - start_time);
60 return(1);
61 } else {
62 TEST_STATUS(" PASSED test -- SUNMatGetID \n", myid);
63 }
64
65 if (myid == 0)
66 PRINT_TIME(" SUNMatGetID Time: %22.15e \n \n", stop_time - start_time);
67
68 return(0);
69 }
70
71
72 /* ----------------------------------------------------------------------
73 * SUNMatClone Test
74 *
75 * NOTE: This routine depends on SUNMatCopy to check matrix data.
76 * --------------------------------------------------------------------*/
Test_SUNMatClone(SUNMatrix A,int myid)77 int Test_SUNMatClone(SUNMatrix A, int myid)
78 {
79 int failure;
80 double start_time, stop_time;
81 realtype tol=10*UNIT_ROUNDOFF;
82 SUNMatrix B;
83
84 /* clone vector */
85 start_time = get_time();
86 B = SUNMatClone(A);
87 stop_time = get_time();
88
89 /* check cloned matrix */
90 if (B == NULL) {
91 TEST_STATUS(">>> FAILED test -- SUNMatClone \n", myid);
92 TEST_STATUS(" After SUNMatClone, B == NULL \n \n", myid);
93 return(1);
94 }
95
96 /* check cloned matrix data */
97 if (!has_data(B)) {
98 TEST_STATUS(">>> FAILED test -- SUNMatClone \n", myid);
99 TEST_STATUS(" Matrix data == NULL \n \n", myid);
100 SUNMatDestroy(B);
101 return(1);
102 }
103
104 failure = SUNMatCopy(A, B);
105 if (failure) {
106 TEST_STATUS(">>> FAILED test -- SUNMatClone (error using SUNMatCopy) \n", myid);
107 SUNMatDestroy(B);
108 return(1);
109 }
110
111 failure = check_matrix(B, A, tol);
112 if (failure) {
113 TEST_STATUS(">>> FAILED test -- SUNMatClone \n", myid);
114 TEST_STATUS(" Failed SUNMatClone check \n \n", myid);
115 SUNMatDestroy(B);
116 return(1);
117 } else {
118 TEST_STATUS(" PASSED test -- SUNMatClone \n", myid);
119 }
120
121 if (myid == 0)
122 PRINT_TIME(" SUNMatClone Time: %22.15e \n \n", stop_time - start_time);
123
124 SUNMatDestroy(B);
125 return(0);
126 }
127
128
129
130 /* ----------------------------------------------------------------------
131 * SUNMatZero Test
132 * --------------------------------------------------------------------*/
Test_SUNMatZero(SUNMatrix A,int myid)133 int Test_SUNMatZero(SUNMatrix A, int myid)
134 {
135 int failure;
136 double start_time, stop_time;
137 realtype tol=10*UNIT_ROUNDOFF;
138 SUNMatrix B;
139
140 /* protect A */
141 B = SUNMatClone(A);
142
143 /* set matrix data to zero */
144 start_time = get_time();
145 failure = SUNMatZero(B);
146 stop_time = get_time();
147
148 if (failure) {
149 TEST_STATUS2(">>> FAILED test -- SUNMatZero returned %d \n", failure, myid);
150 SUNMatDestroy(B);
151 return(1);
152 }
153
154 /* A data should be a vector of zeros */
155 failure = check_matrix_entry(B, ZERO, tol);
156
157 if (failure) {
158 TEST_STATUS(">>> FAILED test -- SUNMatZero check \n", myid);
159 PRINT_TIME(" SUNMatZero Time: %22.15e \n \n", stop_time - start_time);
160 SUNMatDestroy(B);
161 return(1);
162 }
163 else {
164 TEST_STATUS(" PASSED test -- SUNMatZero \n", myid);
165 }
166
167 if (myid == 0)
168 PRINT_TIME(" SUNMatZero Time: %22.15e \n \n", stop_time - start_time);
169
170 SUNMatDestroy(B);
171 return(0);
172 }
173
174
175 /* ----------------------------------------------------------------------
176 * SUNMatCopy Test
177 * --------------------------------------------------------------------*/
Test_SUNMatCopy(SUNMatrix A,int myid)178 int Test_SUNMatCopy(SUNMatrix A, int myid)
179 {
180 int failure;
181 double start_time, stop_time;
182 SUNMatrix B;
183 realtype tol=10*UNIT_ROUNDOFF;
184
185 B = SUNMatClone(A);
186
187 /* copy matrix data */
188 start_time = get_time();
189 failure = SUNMatCopy(A, B);
190 stop_time = get_time();
191
192 if (failure) {
193 TEST_STATUS2(">>> FAILED test -- SUNMatCopy returned %d \n", failure, myid);
194 SUNMatDestroy(B);
195 return(1);
196 }
197
198 /* check matrix entries */
199 failure = check_matrix(B, A, tol);
200
201 if (failure) {
202 TEST_STATUS(">>> FAILED test -- SUNMatCopy \n", myid);
203 PRINT_TIME(" SUNMatCopy Time: %22.15e \n \n", stop_time - start_time);
204 SUNMatDestroy(B);
205 return(1);
206 }
207 else {
208 TEST_STATUS(" PASSED test -- SUNMatCopy \n", myid);
209 }
210
211 if (myid == 0)
212 PRINT_TIME(" SUNMatCopy Time: %22.15e \n \n", stop_time - start_time);
213
214 SUNMatDestroy(B);
215 return(0);
216 }
217
218
219
220 /* ----------------------------------------------------------------------
221 * SUNMatScaleAdd Test: A = c * A + B
222 *
223 * NOTE: Sparse matrices will need additional testing for possibly
224 * different sparsity patterns
225 * --------------------------------------------------------------------*/
Test_SUNMatScaleAdd(SUNMatrix A,SUNMatrix I,int myid)226 int Test_SUNMatScaleAdd(SUNMatrix A, SUNMatrix I, int myid)
227 {
228 int failure;
229 double start_time, stop_time;
230 SUNMatrix B, C, D;
231 realtype tol=10*UNIT_ROUNDOFF;
232
233 /*
234 * Case 1: same sparsity/bandwith pattern
235 */
236
237 /* protect A */
238 B = SUNMatClone(A);
239 failure = SUNMatCopy(A, B);
240 if (failure) {
241 TEST_STATUS2(">>> FAILED test -- SUNMatCopy returned %d \n", failure, myid);
242 SUNMatDestroy(B);
243 return(1);
244 }
245
246 /* fill vector data */
247 start_time = get_time();
248 failure = SUNMatScaleAdd(NEG_ONE, B, B);
249 stop_time = get_time();
250
251 if (failure) {
252 TEST_STATUS2(">>> FAILED test -- SUNMatScaleAdd returned %d \n", failure, myid);
253 SUNMatDestroy(B);
254 return(1);
255 }
256
257 /* check matrix entries */
258 failure = check_matrix_entry(B, ZERO, tol);
259
260 if (failure) {
261 TEST_STATUS(">>> FAILED test -- SUNMatScaleAdd case 1 check \n", myid);
262 PRINT_TIME(" SUNMatScaleAdd Time: %22.15e \n \n", stop_time - start_time);
263 SUNMatDestroy(B);
264 return(1);
265 }
266 else {
267 TEST_STATUS(" PASSED test -- SUNMatScaleAdd case 1 \n", myid);
268 }
269
270 if (myid == 0) {
271 PRINT_TIME(" SUNMatScaleAdd Time: %22.15e \n \n", stop_time - start_time);
272 }
273
274 /*
275 * Case 2: different sparsity/bandwidth patterns
276 */
277 if (is_square(A) && SUNMatGetID(A) != SUNMATRIX_CUSPARSE) {
278
279 /* protect A and I */
280 D = SUNMatClone(A);
281 failure = SUNMatCopy(A, D);
282 if (failure) {
283 TEST_STATUS2(">>> FAILED test -- SUNMatCopy returned %d \n", failure, myid);
284 SUNMatDestroy(B);
285 SUNMatDestroy(D);
286 return(1);
287 }
288 C = SUNMatClone(I);
289 failure = SUNMatCopy(I, C);
290 if (failure) {
291 TEST_STATUS2(">>> FAILED test -- SUNMatCopy returned %d \n", failure, myid);
292 SUNMatDestroy(B);
293 SUNMatDestroy(C);
294 SUNMatDestroy(D);
295 return(1);
296 }
297
298 /* fill B and C */
299 start_time = get_time();
300 failure = SUNMatScaleAdd(ONE, D, I);
301 if (failure) {
302 TEST_STATUS2(">>> FAILED test -- SUNMatScaleAdd returned %d \n", failure, myid);
303 SUNMatDestroy(B);
304 SUNMatDestroy(C);
305 SUNMatDestroy(D);
306 return(1);
307 }
308 failure = SUNMatScaleAdd(ONE, C, A);
309 if (failure) {
310 TEST_STATUS2(">>> FAILED test -- SUNMatScaleAdd returned %d \n", failure, myid);
311 SUNMatDestroy(B);
312 SUNMatDestroy(C);
313 SUNMatDestroy(D);
314 return(1);
315 }
316 stop_time = get_time();
317
318 failure = check_matrix(D, C, tol);
319 if (failure) {
320 TEST_STATUS(">>> FAILED test -- SUNMatScaleAdd case 2 \n", myid);
321 PRINT_TIME(" SUNMatScaleAdd Time: %22.15e \n \n", stop_time - start_time);
322 SUNMatDestroy(B);
323 SUNMatDestroy(C);
324 SUNMatDestroy(D);
325 return(1);
326 }
327 else {
328 TEST_STATUS(" PASSED test -- SUNMatScaleAdd case 2 \n", myid);
329 }
330
331 if (myid == 0)
332 PRINT_TIME(" SUNMatScaleAdd Time: %22.15e \n \n", stop_time - start_time);
333
334 SUNMatDestroy(C);
335 SUNMatDestroy(D);
336 }
337
338 SUNMatDestroy(B);
339
340 return(0);
341 }
342
343
344 /* ----------------------------------------------------------------------
345 * SUNMatScaleAddI Tests
346 *
347 * NOTE: Sparse matrices will need additional testing for possibly
348 * different sparsity patterns
349 * --------------------------------------------------------------------*/
Test_SUNMatScaleAddI(SUNMatrix A,SUNMatrix I,int myid)350 int Test_SUNMatScaleAddI(SUNMatrix A, SUNMatrix I, int myid)
351 {
352 int failure;
353 double start_time, stop_time;
354 SUNMatrix B;
355 realtype tol=10*UNIT_ROUNDOFF;
356
357 /* protect A */
358 B = SUNMatClone(A);
359
360 failure = SUNMatCopy(I, B);
361 if (failure) {
362 TEST_STATUS2(">>> FAILED test -- SUNMatCopy returned %d \n", failure, myid);
363 SUNMatDestroy(B);
364 return(1);
365 }
366
367 /* perform operation */
368 start_time = get_time();
369 failure = SUNMatScaleAddI(NEG_ONE, B);
370 stop_time = get_time();
371
372 if (failure) {
373 TEST_STATUS2(">>> FAILED test -- SUNMatScaleAddI returned %d \n", failure, myid);
374 SUNMatDestroy(B);
375 return(1);
376 }
377
378 /* check matrix */
379 failure = check_matrix_entry(B, ZERO, tol);
380
381 if (failure) {
382 TEST_STATUS(">>> FAILED test -- SUNMatScaleAddI check \n", myid);
383 PRINT_TIME(" SUNMatScaleAddI Time: %22.15e \n \n", stop_time - start_time);
384 SUNMatDestroy(B);
385 return(1);
386 }
387 else {
388 TEST_STATUS(" PASSED test -- SUNMatScaleAddI \n", myid);
389 }
390
391 if (myid == 0)
392 PRINT_TIME(" SUNMatScaleAddI Time: %22.15e \n \n", stop_time - start_time);
393
394 SUNMatDestroy(B);
395 return(0);
396 }
397
398 /* ----------------------------------------------------------------------
399 * SUNMatMatvecSetup Test
400 * --------------------------------------------------------------------*/
Test_SUNMatMatvecSetup(SUNMatrix A,int myid)401 int Test_SUNMatMatvecSetup(SUNMatrix A, int myid)
402 {
403 int failure;
404 double start_time, stop_time;
405
406 if (A->ops->matvecsetup == NULL) {
407 TEST_STATUS(" PASSED test -- SUNMatMatvecSetup not implemented\n", myid);
408 return(0);
409 }
410
411 start_time = get_time();
412 failure = SUNMatMatvecSetup(A);
413 stop_time = get_time();
414
415 if (failure) {
416 TEST_STATUS(">>> FAILED test -- SUNMatMatvecSetup \n", myid);
417 return(1);
418 } else {
419 TEST_STATUS(" PASSED test -- SUNMatMatvecSetup\n", myid);
420 }
421
422 if (myid == 0)
423 PRINT_TIME(" SUNMatMatvecSetup Time: %22.15e \n \n", stop_time - start_time);
424
425 return(0);
426 }
427
428
429 /* ----------------------------------------------------------------------
430 * SUNMatMatvec Test (y should be correct A*x product)
431 * --------------------------------------------------------------------*/
Test_SUNMatMatvec(SUNMatrix A,N_Vector x,N_Vector y,int myid)432 int Test_SUNMatMatvec(SUNMatrix A, N_Vector x, N_Vector y, int myid)
433 {
434 int failure;
435 double start_time, stop_time;
436 SUNMatrix B;
437 N_Vector z, w;
438 realtype tol=100*UNIT_ROUNDOFF;
439
440 if (A->ops->matvec == NULL) {
441 TEST_STATUS(" PASSED test -- SUNMatMatvec not implemented\n", myid);
442 return(0);
443 }
444
445 /* harder tests for square matrices */
446 if (is_square(A)) {
447
448 /* protect A */
449 B = SUNMatClone(A);
450 failure = SUNMatCopy(A, B);
451 if (failure) {
452 TEST_STATUS2(">>> FAILED test -- SUNMatCopy returned %d \n", failure, myid);
453 SUNMatDestroy(B);
454 return(1);
455 }
456
457 /* compute matrix vector product */
458 failure = SUNMatScaleAddI(THREE,B);
459 if (failure) {
460 TEST_STATUS2(">>> FAILED test -- SUNMatScaleAddI returned %d \n", failure, myid);
461 SUNMatDestroy(B);
462 return(1);
463 }
464
465 z = N_VClone(y); /* will be computed with matvec */
466 w = N_VClone(y); /* will be the reference */
467
468 /* Call the Setup function before the Matvec if it exists */
469 if (B->ops->matvecsetup)
470 SUNMatMatvecSetup(B);
471
472 start_time = get_time();
473 failure = SUNMatMatvec(B,x,z); /* z = (3A+I)x */
474 stop_time = get_time();
475
476 if (failure) {
477 TEST_STATUS2(">>> FAILED test -- SUNMatMatvec returned %d \n", failure, myid);
478 SUNMatDestroy(B);
479 return(1);
480 }
481
482 N_VLinearSum(THREE,y,ONE,x,w); /* w = 3x + x */
483
484 failure = check_vector(w,z,tol);
485
486 SUNMatDestroy(B);
487 N_VDestroy(z);
488 N_VDestroy(w);
489
490 } else {
491
492 z = N_VClone(y); /* will be computed with matvec */
493
494 start_time = get_time();
495 failure = SUNMatMatvec(A,x,z); /* z = Ax */
496 stop_time = get_time();
497
498 if (failure) {
499 TEST_STATUS2(">>> FAILED test -- SUNMatMatvec returned %d \n", failure, myid);
500 return(1);
501 }
502
503 failure = check_vector(y,z,tol);
504
505 N_VDestroy(z);
506
507 }
508
509 if (failure) {
510 TEST_STATUS(">>> FAILED test -- SUNMatMatvec check \n", myid);
511 PRINT_TIME(" SUNMatMatvec Time: %22.15e \n \n", stop_time - start_time);
512 return(1);
513 }
514 else {
515 TEST_STATUS(" PASSED test -- SUNMatMatvec \n", myid);
516 }
517
518 if (myid == 0)
519 PRINT_TIME(" SUNMatMatvec Time: %22.15e \n \n", stop_time - start_time);
520
521 return(0);
522 }
523
524
525
526 /* ----------------------------------------------------------------------
527 * SUNMatSpace Test
528 * --------------------------------------------------------------------*/
Test_SUNMatSpace(SUNMatrix A,int myid)529 int Test_SUNMatSpace(SUNMatrix A, int myid)
530 {
531 int failure;
532 double start_time, stop_time;
533 long int lenrw, leniw;
534
535 start_time = get_time();
536 failure = SUNMatSpace(A, &lenrw, &leniw);
537 stop_time = get_time();
538
539 if (failure) {
540 TEST_STATUS(">>> FAILED test -- SUNMatSpace \n", myid);
541 PRINT_TIME(" SUNMatSpace Time: %22.15e \n \n", stop_time - start_time);
542 return(1);
543 } else {
544 TEST_STATUS2(" PASSED test -- SUNMatSpace lenrw=%li ", lenrw, myid);
545 TEST_STATUS2("leniw = %li\n", leniw, myid);
546 }
547
548 if (myid == 0)
549 PRINT_TIME(" SUNMatSpace Time: %22.15e \n \n", stop_time - start_time);
550
551 return(0);
552 }
553
554
555 /* ======================================================================
556 * Private functions
557 * ====================================================================*/
558
559 #if defined( SUNDIALS_HAVE_POSIX_TIMERS) && defined(_POSIX_TIMERS)
560 time_t base_time_tv_sec = 0; /* Base time; makes time values returned
561 by get_time easier to read when
562 printed since they will be zero
563 based.
564 */
565 #endif
566
SetTiming(int onoff)567 void SetTiming(int onoff)
568 {
569 print_time = onoff;
570
571 #if defined( SUNDIALS_HAVE_POSIX_TIMERS) && defined(_POSIX_TIMERS)
572 struct timespec spec;
573 clock_gettime( CLOCK_MONOTONIC_RAW, &spec );
574 base_time_tv_sec = spec.tv_sec;
575 #endif
576 }
577
SetPrintAllRanks(int onoff)578 void SetPrintAllRanks(int onoff)
579 {
580 print_all_ranks = onoff;
581 }
582
583 /* ----------------------------------------------------------------------
584 * Timer
585 * --------------------------------------------------------------------*/
get_time()586 static double get_time()
587 {
588 #if defined( SUNDIALS_HAVE_POSIX_TIMERS) && defined(_POSIX_TIMERS)
589 struct timespec spec;
590 clock_gettime( CLOCK_MONOTONIC_RAW, &spec );
591 double time = (double)(spec.tv_sec - base_time_tv_sec) + ((double)(spec.tv_nsec) / 1E9);
592 #else
593 double time = 0;
594 #endif
595 return time;
596 }
597
598
599