1 /*
2 * -----------------------------------------------------------------
3 * Programmer(s): Daniel Reynolds @ SMU
4 * David Gardner @ LLNL
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 * This is the testing routine to check the SUNMatrix Dense module
17 * implementation.
18 * -----------------------------------------------------------------
19 */
20
21 #include <stdio.h>
22 #include <stdlib.h>
23
24 #include <sundials/sundials_types.h>
25 #include <sunmatrix/sunmatrix_dense.h>
26 #include <nvector/nvector_serial.h>
27 #include <sundials/sundials_math.h>
28 #include "test_sunmatrix.h"
29
30 #if defined(SUNDIALS_EXTENDED_PRECISION)
31 #define GSYM "Lg"
32 #define ESYM "Le"
33 #define FSYM "Lf"
34 #else
35 #define GSYM "g"
36 #define ESYM "e"
37 #define FSYM "f"
38 #endif
39
40 /* ----------------------------------------------------------------------
41 * Main SUNMatrix Testing Routine
42 * --------------------------------------------------------------------*/
main(int argc,char * argv[])43 int main(int argc, char *argv[])
44 {
45 int fails = 0; /* counter for test failures */
46 sunindextype matrows, matcols; /* vector length */
47 N_Vector x, y; /* test vectors */
48 realtype *xdata, *ydata; /* pointers to vector data */
49 SUNMatrix A, I; /* test matrices */
50 realtype *Adata, *Idata; /* pointers to matrix data */
51 int print_timing, square;
52 sunindextype i, j, m, n;
53
54 /* check input and set vector length */
55 if (argc < 4){
56 printf("ERROR: THREE (3) Input required: matrix rows, matrix cols, print timing \n");
57 return(-1);
58 }
59
60 matrows = (sunindextype) atol(argv[1]);
61 if (matrows <= 0) {
62 printf("ERROR: number of rows must be a positive integer \n");
63 return(-1);
64 }
65
66 matcols = (sunindextype) atol(argv[2]);
67 if (matcols <= 0) {
68 printf("ERROR: number of cols must be a positive integer \n");
69 return(-1);
70 }
71
72 print_timing = atoi(argv[3]);
73 SetTiming(print_timing);
74
75 square = (matrows == matcols) ? 1 : 0;
76 printf("\nDense matrix test: size %ld by %ld\n\n",
77 (long int) matrows, (long int) matcols);
78
79 /* Initialize vectors and matrices to NULL */
80 x = NULL;
81 y = NULL;
82 A = NULL;
83 I = NULL;
84
85 /* Create vectors and matrices */
86 x = N_VNew_Serial(matcols);
87 y = N_VNew_Serial(matrows);
88 A = SUNDenseMatrix(matrows, matcols);
89 I = NULL;
90 if (square)
91 I = SUNDenseMatrix(matrows, matcols);
92
93 /* Fill matrices and vectors */
94 Adata = SUNDenseMatrix_Data(A);
95 for(j=0; j < matcols; j++) {
96 for(i=0; i < matrows; i++) {
97 Adata[j*matrows + i] = (j+1)*(i+j);
98 }
99 }
100
101 if (square) {
102 Idata = SUNDenseMatrix_Data(I);
103 for(i=0, j=0; i < matrows; i++, j++) {
104 Idata[j*matrows + i] = ONE;
105 }
106 }
107
108 xdata = N_VGetArrayPointer(x);
109 for(i=0; i < matcols; i++) {
110 xdata[i] = ONE / (i+1);
111 }
112
113 ydata = N_VGetArrayPointer(y);
114 for(i=0; i < matrows; i++) {
115 m = i;
116 n = m + matcols - 1;
117 ydata[i] = HALF*(n+1-m)*(n+m);
118 }
119
120 /* SUNMatrix Tests */
121 fails += Test_SUNMatGetID(A, SUNMATRIX_DENSE, 0);
122 fails += Test_SUNMatClone(A, 0);
123 fails += Test_SUNMatCopy(A, 0);
124 fails += Test_SUNMatZero(A, 0);
125 if (square) {
126 fails += Test_SUNMatScaleAdd(A, I, 0);
127 fails += Test_SUNMatScaleAddI(A, I, 0);
128 }
129 fails += Test_SUNMatMatvec(A, x, y, 0);
130 fails += Test_SUNMatSpace(A, 0);
131
132 /* Print result */
133 if (fails) {
134 printf("FAIL: SUNMatrix module failed %i tests \n \n", fails);
135 printf("\nA =\n");
136 SUNDenseMatrix_Print(A,stdout);
137 if (square) {
138 printf("\nI =\n");
139 SUNDenseMatrix_Print(I,stdout);
140 }
141 printf("\nx =\n");
142 N_VPrint_Serial(x);
143 printf("\ny =\n");
144 N_VPrint_Serial(y);
145 } else {
146 printf("SUCCESS: SUNMatrix module passed all tests \n \n");
147 }
148
149 /* Free vectors and matrices */
150 N_VDestroy(x);
151 N_VDestroy(y);
152 SUNMatDestroy(A);
153 if (square)
154 SUNMatDestroy(I);
155
156 return(fails);
157 }
158
159 /* ----------------------------------------------------------------------
160 * Check matrix
161 * --------------------------------------------------------------------*/
check_matrix(SUNMatrix A,SUNMatrix B,realtype tol)162 int check_matrix(SUNMatrix A, SUNMatrix B, realtype tol)
163 {
164 int failure = 0;
165 realtype *Adata, *Bdata;
166 sunindextype Aldata, Bldata;
167 sunindextype i;
168
169 /* get data pointers */
170 Adata = SUNDenseMatrix_Data(A);
171 Bdata = SUNDenseMatrix_Data(B);
172
173 /* get and check data lengths */
174 Aldata = SUNDenseMatrix_LData(A);
175 Bldata = SUNDenseMatrix_LData(B);
176
177 if (Aldata != Bldata) {
178 printf(">>> ERROR: check_matrix: Different data array lengths \n");
179 return(1);
180 }
181
182 /* compare data */
183 for(i=0; i < Aldata; i++){
184 failure += FNEQ(Adata[i], Bdata[i], tol);
185 }
186
187 if (failure > ZERO)
188 return(1);
189 else
190 return(0);
191 }
192
check_matrix_entry(SUNMatrix A,realtype val,realtype tol)193 int check_matrix_entry(SUNMatrix A, realtype val, realtype tol)
194 {
195 int failure = 0;
196 realtype *Adata;
197 sunindextype Aldata;
198 sunindextype i;
199
200 /* get data pointer */
201 Adata = SUNDenseMatrix_Data(A);
202
203 /* compare data */
204 Aldata = SUNDenseMatrix_LData(A);
205 for(i=0; i < Aldata; i++){
206 failure += FNEQ(Adata[i], val, tol);
207 }
208
209 if (failure > ZERO) {
210 printf("Check_matrix_entry failures:\n");
211 for(i=0; i < Aldata; i++)
212 if (FNEQ(Adata[i], val, tol) != 0)
213 printf(" Adata[%ld] = %"GSYM" != %"GSYM" (err = %"GSYM")\n", (long int) i,
214 Adata[i], val, SUNRabs(Adata[i]-val));
215 }
216
217 if (failure > ZERO)
218 return(1);
219 else
220 return(0);
221 }
222
check_vector(N_Vector x,N_Vector y,realtype tol)223 int check_vector(N_Vector x, N_Vector y, realtype tol)
224 {
225 int failure = 0;
226 realtype *xdata, *ydata;
227 sunindextype xldata, yldata;
228 sunindextype i;
229
230 /* get vector data */
231 xdata = N_VGetArrayPointer(x);
232 ydata = N_VGetArrayPointer(y);
233
234 /* check data lengths */
235 xldata = N_VGetLength(x);
236 yldata = N_VGetLength(y);
237
238 if (xldata != yldata) {
239 printf(">>> ERROR: check_vector: Different data array lengths \n");
240 return(1);
241 }
242
243 /* check vector data */
244 for(i=0; i < xldata; i++)
245 failure += FNEQ(xdata[i], ydata[i], tol);
246
247 if (failure > ZERO) {
248 printf("Check_vector failures:\n");
249 for(i=0; i < xldata; i++)
250 if (FNEQ(xdata[i], ydata[i], tol) != 0)
251 printf(" xdata[%ld] = %"GSYM" != %"GSYM" (err = %"GSYM")\n", (long int) i,
252 xdata[i], ydata[i], SUNRabs(xdata[i]-ydata[i]));
253 }
254
255 if (failure > ZERO)
256 return(1);
257 else
258 return(0);
259 }
260
has_data(SUNMatrix A)261 booleantype has_data(SUNMatrix A)
262 {
263 realtype *Adata = SUNDenseMatrix_Data(A);
264 if (Adata == NULL)
265 return SUNFALSE;
266 else
267 return SUNTRUE;
268 }
269
is_square(SUNMatrix A)270 booleantype is_square(SUNMatrix A)
271 {
272 if (SUNDenseMatrix_Rows(A) == SUNDenseMatrix_Columns(A))
273 return SUNTRUE;
274 else
275 return SUNFALSE;
276 }
277