1 /* -----------------------------------------------------------------------------
2 * Programmer(s): David J. Gardner @ LLNL
3 * -----------------------------------------------------------------------------
4 * SUNDIALS Copyright Start
5 * Copyright (c) 2002-2020, Lawrence Livermore National Security
6 * and Southern Methodist University.
7 * All rights reserved.
8 *
9 * See the top-level LICENSE and NOTICE files for details.
10 *
11 * SPDX-License-Identifier: BSD-3-Clause
12 * SUNDIALS Copyright End
13 * -----------------------------------------------------------------------------
14 * This is the implementation file for a vector wrapper for an array of NVECTORS
15 * ---------------------------------------------------------------------------*/
16
17 #include <stdio.h>
18 #include <stdlib.h>
19
20 #include <stdarg.h>
21 #include <string.h>
22
23 #include <sundials/sundials_nvector.h>
24 #include <sundials/sundials_nvector_senswrapper.h>
25
26 #define ZERO RCONST(0.0)
27
28 /*==============================================================================
29 Constructors
30 ============================================================================*/
31
32 /*------------------------------------------------------------------------------
33 create a new empty vector wrapper with space for <nvecs> vectors
34 ----------------------------------------------------------------------------*/
N_VNewEmpty_SensWrapper(int nvecs)35 N_Vector N_VNewEmpty_SensWrapper(int nvecs)
36 {
37 int i;
38 N_Vector v;
39 N_Vector_Ops ops;
40 N_VectorContent_SensWrapper content;
41
42 /* return if wrapper is empty */
43 if (nvecs < 1) return(NULL);
44
45 /* create vector */
46 v = NULL;
47 v = (N_Vector) malloc(sizeof *v);
48 if (v == NULL) return(NULL);
49
50 /* create vector operation structure */
51 ops = NULL;
52 ops = (N_Vector_Ops) malloc(sizeof *ops);
53 if (ops == NULL) {free(v); return(NULL);}
54
55 ops->nvgetvectorid = NULL;
56 ops->nvclone = N_VClone_SensWrapper;
57 ops->nvcloneempty = N_VCloneEmpty_SensWrapper;
58 ops->nvdestroy = N_VDestroy_SensWrapper;
59 ops->nvspace = NULL;
60 ops->nvgetarraypointer = NULL;
61 ops->nvsetarraypointer = NULL;
62
63 /* standard vector operations */
64 ops->nvlinearsum = N_VLinearSum_SensWrapper;
65 ops->nvconst = N_VConst_SensWrapper;
66 ops->nvprod = N_VProd_SensWrapper;
67 ops->nvdiv = N_VDiv_SensWrapper;
68 ops->nvscale = N_VScale_SensWrapper;
69 ops->nvabs = N_VAbs_SensWrapper;
70 ops->nvinv = N_VInv_SensWrapper;
71 ops->nvaddconst = N_VAddConst_SensWrapper;
72 ops->nvdotprod = N_VDotProd_SensWrapper;
73 ops->nvmaxnorm = N_VMaxNorm_SensWrapper;
74 ops->nvwrmsnormmask = N_VWrmsNormMask_SensWrapper;
75 ops->nvwrmsnorm = N_VWrmsNorm_SensWrapper;
76 ops->nvmin = N_VMin_SensWrapper;
77 ops->nvwl2norm = N_VWL2Norm_SensWrapper;
78 ops->nvl1norm = N_VL1Norm_SensWrapper;
79 ops->nvcompare = N_VCompare_SensWrapper;
80 ops->nvinvtest = N_VInvTest_SensWrapper;
81 ops->nvconstrmask = N_VConstrMask_SensWrapper;
82 ops->nvminquotient = N_VMinQuotient_SensWrapper;
83
84 /* fused vector operations */
85 ops->nvlinearcombination = NULL;
86 ops->nvscaleaddmulti = NULL;
87 ops->nvdotprodmulti = NULL;
88
89 /* vector array operations */
90 ops->nvlinearsumvectorarray = NULL;
91 ops->nvscalevectorarray = NULL;
92 ops->nvconstvectorarray = NULL;
93 ops->nvwrmsnormvectorarray = NULL;
94 ops->nvwrmsnormmaskvectorarray = NULL;
95 ops->nvscaleaddmultivectorarray = NULL;
96 ops->nvlinearcombinationvectorarray = NULL;
97
98 /* create content */
99 content = NULL;
100 content = (N_VectorContent_SensWrapper) malloc(sizeof *content);
101 if (content == NULL) {free(ops); free(v); return(NULL);}
102
103 content->nvecs = nvecs;
104 content->own_vecs = SUNFALSE;
105 content->vecs = NULL;
106 content->vecs = (N_Vector*) malloc(nvecs * sizeof(N_Vector));
107 if (content->vecs == NULL) {free(ops); free(v); free(content); return(NULL);}
108
109 /* initialize vector array to null */
110 for (i=0; i < nvecs; i++)
111 content->vecs[i] = NULL;
112
113 /* attach content and ops */
114 v->content = content;
115 v->ops = ops;
116
117 return(v);
118 }
119
120
N_VNew_SensWrapper(int count,N_Vector w)121 N_Vector N_VNew_SensWrapper(int count, N_Vector w)
122 {
123 N_Vector v;
124 int i;
125
126 v = NULL;
127 v = N_VNewEmpty_SensWrapper(count);
128 if (v == NULL) return(NULL);
129
130 for (i=0; i < NV_NVECS_SW(v); i++) {
131 NV_VEC_SW(v,i) = N_VClone(w);
132 if (NV_VEC_SW(v,i) == NULL) { N_VDestroy(v); return(NULL); }
133 }
134
135 /* update own vectors status */
136 NV_OWN_VECS_SW(v) = SUNTRUE;
137
138 return(v);
139 }
140
141
142 /*==============================================================================
143 Clone operations
144 ============================================================================*/
145
146 /*------------------------------------------------------------------------------
147 create an empty clone of the vector wrapper w
148 ----------------------------------------------------------------------------*/
N_VCloneEmpty_SensWrapper(N_Vector w)149 N_Vector N_VCloneEmpty_SensWrapper(N_Vector w)
150 {
151 int i;
152 N_Vector v;
153 N_Vector_Ops ops;
154 N_VectorContent_SensWrapper content;
155
156 if (w == NULL) return(NULL);
157
158 if (NV_NVECS_SW(w) < 1) return(NULL);
159
160 /* create vector */
161 v = NULL;
162 v = (N_Vector) malloc(sizeof *v);
163 if (v == NULL) return(NULL);
164
165 /* create vector operation structure */
166 ops = NULL;
167 ops = (N_Vector_Ops) malloc(sizeof *ops);
168 if (ops == NULL) { free(v); return(NULL); }
169
170 ops->nvgetvectorid = w->ops->nvgetvectorid;
171 ops->nvclone = w->ops->nvclone;
172 ops->nvcloneempty = w->ops->nvcloneempty;
173 ops->nvdestroy = w->ops->nvdestroy;
174 ops->nvspace = w->ops->nvspace;
175 ops->nvgetarraypointer = w->ops->nvgetarraypointer;
176 ops->nvsetarraypointer = w->ops->nvsetarraypointer;
177
178 /* standard vector operations */
179 ops->nvlinearsum = w->ops->nvlinearsum;
180 ops->nvconst = w->ops->nvconst;
181 ops->nvprod = w->ops->nvprod;
182 ops->nvdiv = w->ops->nvdiv;
183 ops->nvscale = w->ops->nvscale;
184 ops->nvabs = w->ops->nvabs;
185 ops->nvinv = w->ops->nvinv;
186 ops->nvaddconst = w->ops->nvaddconst;
187 ops->nvdotprod = w->ops->nvdotprod;
188 ops->nvmaxnorm = w->ops->nvmaxnorm;
189 ops->nvwrmsnormmask = w->ops->nvwrmsnormmask;
190 ops->nvwrmsnorm = w->ops->nvwrmsnorm;
191 ops->nvmin = w->ops->nvmin;
192 ops->nvwl2norm = w->ops->nvwl2norm;
193 ops->nvl1norm = w->ops->nvl1norm;
194 ops->nvcompare = w->ops->nvcompare;
195 ops->nvinvtest = w->ops->nvinvtest;
196 ops->nvconstrmask = w->ops->nvconstrmask;
197 ops->nvminquotient = w->ops->nvminquotient;
198
199 /* fused vector operations */
200 ops->nvlinearcombination = w->ops->nvlinearcombination;
201 ops->nvscaleaddmulti = w->ops->nvscaleaddmulti;
202 ops->nvdotprodmulti = w->ops->nvdotprodmulti;
203
204 /* vector array operations */
205 ops->nvlinearsumvectorarray = w->ops->nvlinearsumvectorarray;
206 ops->nvscalevectorarray = w->ops->nvscalevectorarray;
207 ops->nvconstvectorarray = w->ops->nvconstvectorarray;
208 ops->nvwrmsnormvectorarray = w->ops->nvwrmsnormvectorarray;
209 ops->nvwrmsnormmaskvectorarray = w->ops->nvwrmsnormmaskvectorarray;
210 ops->nvscaleaddmultivectorarray = w->ops->nvscaleaddmultivectorarray;
211 ops->nvlinearcombinationvectorarray = w->ops->nvlinearcombinationvectorarray;
212
213 /* Create content */
214 content = NULL;
215 content = (N_VectorContent_SensWrapper) malloc(sizeof *content);
216 if (content == NULL) { free(ops); free(v); return(NULL); }
217
218 content->nvecs = NV_NVECS_SW(w);
219 content->own_vecs = SUNFALSE;
220 content->vecs = NULL;
221 content->vecs = (N_Vector*) malloc(NV_NVECS_SW(w) * sizeof(N_Vector));
222 if (content->vecs == NULL) {free(ops); free(v); free(content); return(NULL);}
223
224 /* initialize vector array to null */
225 for (i=0; i < NV_NVECS_SW(w); i++)
226 content->vecs[i] = NULL;
227
228 /* Attach content and ops */
229 v->content = content;
230 v->ops = ops;
231
232 return(v);
233 }
234
235
236 /*------------------------------------------------------------------------------
237 create a clone of the vector wrapper w
238 ----------------------------------------------------------------------------*/
N_VClone_SensWrapper(N_Vector w)239 N_Vector N_VClone_SensWrapper(N_Vector w)
240 {
241 N_Vector v;
242 int i;
243
244 /* create empty wrapper */
245 v = NULL;
246 v = N_VCloneEmpty_SensWrapper(w);
247 if (v == NULL) return(NULL);
248
249 /* update own vectors status */
250 NV_OWN_VECS_SW(v) = SUNTRUE;
251
252 /* allocate arrays */
253 for (i=0; i < NV_NVECS_SW(v); i++) {
254 NV_VEC_SW(v,i) = N_VClone(NV_VEC_SW(w,i));
255 if (NV_VEC_SW(v,i) == NULL) { N_VDestroy(v); return(NULL); }
256 }
257
258 return(v);
259 }
260
261
262 /*==============================================================================
263 Destructor
264 ============================================================================*/
265
N_VDestroy_SensWrapper(N_Vector v)266 void N_VDestroy_SensWrapper(N_Vector v)
267 {
268 int i;
269
270 if (NV_OWN_VECS_SW(v) == SUNTRUE) {
271 for (i=0; i < NV_NVECS_SW(v); i++) {
272 if (NV_VEC_SW(v,i)) N_VDestroy(NV_VEC_SW(v,i));
273 NV_VEC_SW(v,i) = NULL;
274 }
275 }
276
277 free(NV_VECS_SW(v)); NV_VECS_SW(v) = NULL;
278 free(v->content); v->content = NULL;
279 free(v->ops); v->ops = NULL;
280 free(v); v = NULL;
281
282 return;
283 }
284
285
286 /*==============================================================================
287 Standard vector operations
288 ============================================================================*/
289
N_VLinearSum_SensWrapper(realtype a,N_Vector x,realtype b,N_Vector y,N_Vector z)290 void N_VLinearSum_SensWrapper(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z)
291 {
292 int i;
293
294 for (i=0; i < NV_NVECS_SW(x); i++)
295 N_VLinearSum(a, NV_VEC_SW(x,i), b, NV_VEC_SW(y,i), NV_VEC_SW(z,i));
296
297 return;
298 }
299
300
N_VConst_SensWrapper(realtype c,N_Vector z)301 void N_VConst_SensWrapper(realtype c, N_Vector z)
302 {
303 int i;
304
305 for (i=0; i < NV_NVECS_SW(z); i++)
306 N_VConst(c, NV_VEC_SW(z,i));
307
308 return;
309 }
310
311
N_VProd_SensWrapper(N_Vector x,N_Vector y,N_Vector z)312 void N_VProd_SensWrapper(N_Vector x, N_Vector y, N_Vector z)
313 {
314 int i;
315
316 for (i=0; i < NV_NVECS_SW(x); i++)
317 N_VProd(NV_VEC_SW(x,i), NV_VEC_SW(y,i), NV_VEC_SW(z,i));
318
319 return;
320 }
321
322
N_VDiv_SensWrapper(N_Vector x,N_Vector y,N_Vector z)323 void N_VDiv_SensWrapper(N_Vector x, N_Vector y, N_Vector z)
324 {
325 int i;
326
327 for (i=0; i < NV_NVECS_SW(x); i++)
328 N_VDiv(NV_VEC_SW(x,i), NV_VEC_SW(y,i), NV_VEC_SW(z,i));
329
330 return;
331 }
332
333
N_VScale_SensWrapper(realtype c,N_Vector x,N_Vector z)334 void N_VScale_SensWrapper(realtype c, N_Vector x, N_Vector z)
335 {
336 int i;
337
338 for (i=0; i < NV_NVECS_SW(x); i++)
339 N_VScale(c, NV_VEC_SW(x,i), NV_VEC_SW(z,i));
340
341 return;
342 }
343
344
N_VAbs_SensWrapper(N_Vector x,N_Vector z)345 void N_VAbs_SensWrapper(N_Vector x, N_Vector z)
346 {
347 int i;
348
349 for (i=0; i < NV_NVECS_SW(x); i++)
350 N_VAbs(NV_VEC_SW(x,i), NV_VEC_SW(z,i));
351
352 return;
353 }
354
355
N_VInv_SensWrapper(N_Vector x,N_Vector z)356 void N_VInv_SensWrapper(N_Vector x, N_Vector z)
357 {
358 int i;
359
360 for (i=0; i < NV_NVECS_SW(x); i++)
361 N_VInv(NV_VEC_SW(x,i), NV_VEC_SW(z,i));
362
363 return;
364 }
365
366
N_VAddConst_SensWrapper(N_Vector x,realtype b,N_Vector z)367 void N_VAddConst_SensWrapper(N_Vector x, realtype b, N_Vector z)
368 {
369 int i;
370
371 for (i=0; i < NV_NVECS_SW(x); i++)
372 N_VAddConst(NV_VEC_SW(x,i), b, NV_VEC_SW(z,i));
373
374 return;
375 }
376
377
N_VDotProd_SensWrapper(N_Vector x,N_Vector y)378 realtype N_VDotProd_SensWrapper(N_Vector x, N_Vector y)
379 {
380 int i;
381 realtype sum;
382
383 sum = ZERO;
384
385 for (i=0; i < NV_NVECS_SW(x); i++)
386 sum += N_VDotProd(NV_VEC_SW(x,i), NV_VEC_SW(y,i));
387
388 return(sum);
389 }
390
391
N_VMaxNorm_SensWrapper(N_Vector x)392 realtype N_VMaxNorm_SensWrapper(N_Vector x)
393 {
394 int i;
395 realtype max, tmp;
396
397 max = ZERO;
398
399 for (i=0; i < NV_NVECS_SW(x); i++) {
400 tmp = N_VMaxNorm(NV_VEC_SW(x,i));
401 if (tmp > max) max = tmp;
402 }
403
404 return(max);
405 }
406
407
N_VWrmsNorm_SensWrapper(N_Vector x,N_Vector w)408 realtype N_VWrmsNorm_SensWrapper(N_Vector x, N_Vector w)
409 {
410 int i;
411 realtype nrm, tmp;
412
413 nrm = ZERO;
414
415 for (i=0; i < NV_NVECS_SW(x); i++) {
416 tmp = N_VWrmsNorm(NV_VEC_SW(x,i), NV_VEC_SW(w,i));
417 if (tmp > nrm) nrm = tmp;
418 }
419
420 return(nrm);
421 }
422
423
N_VWrmsNormMask_SensWrapper(N_Vector x,N_Vector w,N_Vector id)424 realtype N_VWrmsNormMask_SensWrapper(N_Vector x, N_Vector w, N_Vector id)
425 {
426 int i;
427 realtype nrm, tmp;
428
429 nrm = ZERO;
430
431 for (i=0; i < NV_NVECS_SW(x); i++) {
432 tmp = N_VWrmsNormMask(NV_VEC_SW(x,i), NV_VEC_SW(w,i), NV_VEC_SW(id,i));
433 if (tmp > nrm) nrm = tmp;
434 }
435
436 return(nrm);
437 }
438
439
N_VMin_SensWrapper(N_Vector x)440 realtype N_VMin_SensWrapper(N_Vector x)
441 {
442 int i;
443 realtype min, tmp;
444
445 min = N_VMin(NV_VEC_SW(x,0));
446
447 for (i=1; i < NV_NVECS_SW(x); i++) {
448 tmp = N_VMin(NV_VEC_SW(x,i));
449 if (tmp < min) min = tmp;
450 }
451
452 return(min);
453 }
454
455
N_VWL2Norm_SensWrapper(N_Vector x,N_Vector w)456 realtype N_VWL2Norm_SensWrapper(N_Vector x, N_Vector w)
457 {
458 int i;
459 realtype nrm, tmp;
460
461 nrm = ZERO;
462
463 for (i=0; i < NV_NVECS_SW(x); i++) {
464 tmp = N_VWL2Norm(NV_VEC_SW(x,i), NV_VEC_SW(w,i));
465 if (tmp > nrm) nrm = tmp;
466 }
467
468 return(nrm);
469 }
470
471
N_VL1Norm_SensWrapper(N_Vector x)472 realtype N_VL1Norm_SensWrapper(N_Vector x)
473 {
474 int i;
475 realtype nrm, tmp;
476
477 nrm = ZERO;
478
479 for (i=0; i < NV_NVECS_SW(x); i++) {
480 tmp = N_VL1Norm(NV_VEC_SW(x,i));
481 if (tmp > nrm) nrm = tmp;
482 }
483
484 return(nrm);
485 }
486
487
N_VCompare_SensWrapper(realtype c,N_Vector x,N_Vector z)488 void N_VCompare_SensWrapper(realtype c, N_Vector x, N_Vector z)
489 {
490 int i;
491
492 for (i=0; i < NV_NVECS_SW(x); i++)
493 N_VCompare(c, NV_VEC_SW(x,i), NV_VEC_SW(z,i));
494
495 return;
496 }
497
498
N_VInvTest_SensWrapper(N_Vector x,N_Vector z)499 booleantype N_VInvTest_SensWrapper(N_Vector x, N_Vector z)
500 {
501 int i;
502 booleantype no_zero_found, tmp;
503
504 no_zero_found = SUNTRUE;
505
506 for (i=0; i < NV_NVECS_SW(x); i++) {
507 tmp = N_VInvTest(NV_VEC_SW(x,i), NV_VEC_SW(z,i));
508 if (tmp != SUNTRUE) no_zero_found = SUNFALSE;
509 }
510
511 return(no_zero_found);
512 }
513
514
N_VConstrMask_SensWrapper(N_Vector c,N_Vector x,N_Vector m)515 booleantype N_VConstrMask_SensWrapper(N_Vector c, N_Vector x, N_Vector m)
516 {
517 int i;
518 booleantype test, tmp;
519
520 test = SUNTRUE;
521
522 for (i=0; i < NV_NVECS_SW(x); i++) {
523 tmp = N_VConstrMask(c, NV_VEC_SW(x,i), NV_VEC_SW(m,i));
524 if (tmp != SUNTRUE) test = SUNFALSE;
525 }
526
527 return(test);
528 }
529
530
N_VMinQuotient_SensWrapper(N_Vector num,N_Vector denom)531 realtype N_VMinQuotient_SensWrapper(N_Vector num, N_Vector denom)
532 {
533 int i;
534 realtype min, tmp;
535
536 min = N_VMinQuotient(NV_VEC_SW(num,0), NV_VEC_SW(denom,0));
537
538 for (i=1; i < NV_NVECS_SW(num); i++) {
539 tmp = N_VMinQuotient(NV_VEC_SW(num,i), NV_VEC_SW(denom,i));
540 if (tmp < min) min = tmp;
541 }
542
543 return(min);
544 }
545