1 /* -----------------------------------------------------------------
2 * Programmer(s): Radu Serban, Aaron Collier, and
3 * David J. Gardner @ LLNL
4 * Daniel R. Reynolds @ SMU
5 * -----------------------------------------------------------------
6 * SUNDIALS Copyright Start
7 * Copyright (c) 2002-2021, 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 implementation file for a generic NVECTOR package.
17 * It contains the implementation of the N_Vector operations listed
18 * in nvector.h.
19 * -----------------------------------------------------------------*/
20
21 #include <stdio.h>
22 #include <stdlib.h>
23
24 #include <sundials/sundials_nvector.h>
25
26 /* -----------------------------------------------------------------
27 * Create an empty NVector object
28 * -----------------------------------------------------------------*/
29
N_VNewEmpty()30 N_Vector N_VNewEmpty()
31 {
32 N_Vector v;
33 N_Vector_Ops ops;
34
35 /* create vector object */
36 v = NULL;
37 v = (N_Vector) malloc(sizeof *v);
38 if (v == NULL) return(NULL);
39
40 /* create vector ops structure */
41 ops = NULL;
42 ops = (N_Vector_Ops) malloc(sizeof *ops);
43 if (ops == NULL) { free(v); return(NULL); }
44
45 /* initialize operations to NULL */
46
47 /* constructors, destructors, and utility operations */
48 ops->nvgetvectorid = NULL;
49 ops->nvclone = NULL;
50 ops->nvcloneempty = NULL;
51 ops->nvdestroy = NULL;
52 ops->nvspace = NULL;
53 ops->nvgetarraypointer = NULL;
54 ops->nvgetdevicearraypointer = NULL;
55 ops->nvsetarraypointer = NULL;
56 ops->nvgetcommunicator = NULL;
57 ops->nvgetlength = NULL;
58
59 /* standard vector operations */
60 ops->nvlinearsum = NULL;
61 ops->nvconst = NULL;
62 ops->nvprod = NULL;
63 ops->nvdiv = NULL;
64 ops->nvscale = NULL;
65 ops->nvabs = NULL;
66 ops->nvinv = NULL;
67 ops->nvaddconst = NULL;
68 ops->nvdotprod = NULL;
69 ops->nvmaxnorm = NULL;
70 ops->nvwrmsnormmask = NULL;
71 ops->nvwrmsnorm = NULL;
72 ops->nvmin = NULL;
73 ops->nvwl2norm = NULL;
74 ops->nvl1norm = NULL;
75 ops->nvcompare = NULL;
76 ops->nvinvtest = NULL;
77 ops->nvconstrmask = NULL;
78 ops->nvminquotient = NULL;
79
80 /* fused vector operations (optional) */
81 ops->nvlinearcombination = NULL;
82 ops->nvscaleaddmulti = NULL;
83 ops->nvdotprodmulti = NULL;
84
85 /* vector array operations (optional) */
86 ops->nvlinearsumvectorarray = NULL;
87 ops->nvscalevectorarray = NULL;
88 ops->nvconstvectorarray = NULL;
89 ops->nvwrmsnormvectorarray = NULL;
90 ops->nvwrmsnormmaskvectorarray = NULL;
91 ops->nvscaleaddmultivectorarray = NULL;
92 ops->nvlinearcombinationvectorarray = NULL;
93
94 /* local reduction operations (optional) */
95 ops->nvdotprodlocal = NULL;
96 ops->nvmaxnormlocal = NULL;
97 ops->nvminlocal = NULL;
98 ops->nvl1normlocal = NULL;
99 ops->nvinvtestlocal = NULL;
100 ops->nvconstrmasklocal = NULL;
101 ops->nvminquotientlocal = NULL;
102 ops->nvwsqrsumlocal = NULL;
103 ops->nvwsqrsummasklocal = NULL;
104
105 /* XBraid interface operations */
106 ops->nvbufsize = NULL;
107 ops->nvbufpack = NULL;
108 ops->nvbufunpack = NULL;
109
110 /* debugging functions (called when SUNDIALS_DEBUG_PRINTVEC is defined) */
111 ops->nvprint = NULL;
112 ops->nvprintfile = NULL;
113
114 /* attach ops and initialize content to NULL */
115 v->ops = ops;
116 v->content = NULL;
117
118 return(v);
119 }
120
121
122 /* -----------------------------------------------------------------
123 * Free a generic N_Vector (assumes content is already empty)
124 * -----------------------------------------------------------------*/
125
N_VFreeEmpty(N_Vector v)126 void N_VFreeEmpty(N_Vector v)
127 {
128 if (v == NULL) return;
129
130 /* free non-NULL ops structure */
131 if (v->ops) free(v->ops);
132 v->ops = NULL;
133
134 /* free overall N_Vector object and return */
135 free(v);
136 return;
137 }
138
139
140 /* -----------------------------------------------------------------
141 * Copy a vector 'ops' structure
142 * -----------------------------------------------------------------*/
143
N_VCopyOps(N_Vector w,N_Vector v)144 int N_VCopyOps(N_Vector w, N_Vector v)
145 {
146 /* Check that ops structures exist */
147 if (w == NULL || v == NULL) return(-1);
148 if (w->ops == NULL || v->ops == NULL) return(-1);
149
150 /* Copy ops from w to v */
151
152 /* constructors, destructors, and utility operations */
153 v->ops->nvgetvectorid = w->ops->nvgetvectorid;
154 v->ops->nvclone = w->ops->nvclone;
155 v->ops->nvcloneempty = w->ops->nvcloneempty;
156 v->ops->nvdestroy = w->ops->nvdestroy;
157 v->ops->nvspace = w->ops->nvspace;
158 v->ops->nvgetarraypointer = w->ops->nvgetarraypointer;
159 v->ops->nvgetdevicearraypointer = w->ops->nvgetdevicearraypointer;
160 v->ops->nvsetarraypointer = w->ops->nvsetarraypointer;
161 v->ops->nvgetcommunicator = w->ops->nvgetcommunicator;
162 v->ops->nvgetlength = w->ops->nvgetlength;
163
164 /* standard vector operations */
165 v->ops->nvlinearsum = w->ops->nvlinearsum;
166 v->ops->nvconst = w->ops->nvconst;
167 v->ops->nvprod = w->ops->nvprod;
168 v->ops->nvdiv = w->ops->nvdiv;
169 v->ops->nvscale = w->ops->nvscale;
170 v->ops->nvabs = w->ops->nvabs;
171 v->ops->nvinv = w->ops->nvinv;
172 v->ops->nvaddconst = w->ops->nvaddconst;
173 v->ops->nvdotprod = w->ops->nvdotprod;
174 v->ops->nvmaxnorm = w->ops->nvmaxnorm;
175 v->ops->nvwrmsnormmask = w->ops->nvwrmsnormmask;
176 v->ops->nvwrmsnorm = w->ops->nvwrmsnorm;
177 v->ops->nvmin = w->ops->nvmin;
178 v->ops->nvwl2norm = w->ops->nvwl2norm;
179 v->ops->nvl1norm = w->ops->nvl1norm;
180 v->ops->nvcompare = w->ops->nvcompare;
181 v->ops->nvinvtest = w->ops->nvinvtest;
182 v->ops->nvconstrmask = w->ops->nvconstrmask;
183 v->ops->nvminquotient = w->ops->nvminquotient;
184
185 /* fused vector operations */
186 v->ops->nvlinearcombination = w->ops->nvlinearcombination;
187 v->ops->nvscaleaddmulti = w->ops->nvscaleaddmulti;
188 v->ops->nvdotprodmulti = w->ops->nvdotprodmulti;
189
190 /* vector array operations */
191 v->ops->nvlinearsumvectorarray = w->ops->nvlinearsumvectorarray;
192 v->ops->nvscalevectorarray = w->ops->nvscalevectorarray;
193 v->ops->nvconstvectorarray = w->ops->nvconstvectorarray;
194 v->ops->nvwrmsnormvectorarray = w->ops->nvwrmsnormvectorarray;
195 v->ops->nvwrmsnormmaskvectorarray = w->ops->nvwrmsnormmaskvectorarray;
196 v->ops->nvscaleaddmultivectorarray = w->ops->nvscaleaddmultivectorarray;
197 v->ops->nvlinearcombinationvectorarray = w->ops->nvlinearcombinationvectorarray;
198
199 /* local reduction operations */
200 v->ops->nvdotprodlocal = w->ops->nvdotprodlocal;
201 v->ops->nvmaxnormlocal = w->ops->nvmaxnormlocal;
202 v->ops->nvminlocal = w->ops->nvminlocal;
203 v->ops->nvl1normlocal = w->ops->nvl1normlocal;
204 v->ops->nvinvtestlocal = w->ops->nvinvtestlocal;
205 v->ops->nvconstrmasklocal = w->ops->nvconstrmasklocal;
206 v->ops->nvminquotientlocal = w->ops->nvminquotientlocal;
207 v->ops->nvwsqrsumlocal = w->ops->nvwsqrsumlocal;
208 v->ops->nvwsqrsummasklocal = w->ops->nvwsqrsummasklocal;
209
210 /* XBraid interface operations */
211 v->ops->nvbufsize = w->ops->nvbufsize;
212 v->ops->nvbufpack = w->ops->nvbufpack;
213 v->ops->nvbufunpack = w->ops->nvbufunpack;
214
215 /* debugging functions (called when SUNDIALS_DEBUG_PRINTVEC is defined) */
216 v->ops->nvprint = w->ops->nvprint;
217 v->ops->nvprintfile = w->ops->nvprintfile;
218
219 return(0);
220 }
221
222 /* -----------------------------------------------------------------
223 * Functions in the 'ops' structure
224 * -----------------------------------------------------------------*/
225
N_VGetVectorID(N_Vector w)226 N_Vector_ID N_VGetVectorID(N_Vector w)
227 {
228 return(w->ops->nvgetvectorid(w));
229 }
230
N_VClone(N_Vector w)231 N_Vector N_VClone(N_Vector w)
232 {
233 return(w->ops->nvclone(w));
234 }
235
N_VCloneEmpty(N_Vector w)236 N_Vector N_VCloneEmpty(N_Vector w)
237 {
238 return(w->ops->nvcloneempty(w));
239 }
240
N_VDestroy(N_Vector v)241 void N_VDestroy(N_Vector v)
242 {
243 if (v == NULL) return;
244
245 /* if the destroy operation exists use it */
246 if (v->ops)
247 if (v->ops->nvdestroy) { v->ops->nvdestroy(v); return; }
248
249 /* if we reach this point, either ops == NULL or nvdestroy == NULL,
250 try to cleanup by freeing the content, ops, and vector */
251 if (v->content) { free(v->content); v->content = NULL; }
252 if (v->ops) { free(v->ops); v->ops = NULL; }
253 free(v); v = NULL;
254
255 return;
256 }
257
N_VSpace(N_Vector v,sunindextype * lrw,sunindextype * liw)258 void N_VSpace(N_Vector v, sunindextype *lrw, sunindextype *liw)
259 {
260 v->ops->nvspace(v, lrw, liw);
261 return;
262 }
263
N_VGetArrayPointer(N_Vector v)264 realtype *N_VGetArrayPointer(N_Vector v)
265 {
266 return((realtype *) v->ops->nvgetarraypointer(v));
267 }
268
N_VGetDeviceArrayPointer(N_Vector v)269 realtype *N_VGetDeviceArrayPointer(N_Vector v)
270 {
271 if (v->ops->nvgetdevicearraypointer)
272 return((realtype *) v->ops->nvgetdevicearraypointer(v));
273 else
274 return(NULL);
275 }
276
N_VSetArrayPointer(realtype * v_data,N_Vector v)277 void N_VSetArrayPointer(realtype *v_data, N_Vector v)
278 {
279 v->ops->nvsetarraypointer(v_data, v);
280 return;
281 }
282
N_VGetCommunicator(N_Vector v)283 void *N_VGetCommunicator(N_Vector v)
284 {
285 if (v->ops->nvgetcommunicator)
286 return(v->ops->nvgetcommunicator(v));
287 else
288 return(NULL);
289 }
290
N_VGetLength(N_Vector v)291 sunindextype N_VGetLength(N_Vector v)
292 {
293 return((sunindextype) v->ops->nvgetlength(v));
294 }
295
296 /* -----------------------------------------------------------------
297 * standard vector operations
298 * -----------------------------------------------------------------*/
299
N_VLinearSum(realtype a,N_Vector x,realtype b,N_Vector y,N_Vector z)300 void N_VLinearSum(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z)
301 {
302 z->ops->nvlinearsum(a, x, b, y, z);
303 return;
304 }
305
N_VConst(realtype c,N_Vector z)306 void N_VConst(realtype c, N_Vector z)
307 {
308 z->ops->nvconst(c, z);
309 return;
310 }
311
N_VProd(N_Vector x,N_Vector y,N_Vector z)312 void N_VProd(N_Vector x, N_Vector y, N_Vector z)
313 {
314 z->ops->nvprod(x, y, z);
315 return;
316 }
317
N_VDiv(N_Vector x,N_Vector y,N_Vector z)318 void N_VDiv(N_Vector x, N_Vector y, N_Vector z)
319 {
320 z->ops->nvdiv(x, y, z);
321 return;
322 }
323
N_VScale(realtype c,N_Vector x,N_Vector z)324 void N_VScale(realtype c, N_Vector x, N_Vector z)
325 {
326 z->ops->nvscale(c, x, z);
327 return;
328 }
329
N_VAbs(N_Vector x,N_Vector z)330 void N_VAbs(N_Vector x, N_Vector z)
331 {
332 z->ops->nvabs(x, z);
333 return;
334 }
335
N_VInv(N_Vector x,N_Vector z)336 void N_VInv(N_Vector x, N_Vector z)
337 {
338 z->ops->nvinv(x, z);
339 return;
340 }
341
N_VAddConst(N_Vector x,realtype b,N_Vector z)342 void N_VAddConst(N_Vector x, realtype b, N_Vector z)
343 {
344 z->ops->nvaddconst(x, b, z);
345 return;
346 }
347
N_VDotProd(N_Vector x,N_Vector y)348 realtype N_VDotProd(N_Vector x, N_Vector y)
349 {
350 return((realtype) y->ops->nvdotprod(x, y));
351 }
352
N_VMaxNorm(N_Vector x)353 realtype N_VMaxNorm(N_Vector x)
354 {
355 return((realtype) x->ops->nvmaxnorm(x));
356 }
357
N_VWrmsNorm(N_Vector x,N_Vector w)358 realtype N_VWrmsNorm(N_Vector x, N_Vector w)
359 {
360 return((realtype) x->ops->nvwrmsnorm(x, w));
361 }
362
N_VWrmsNormMask(N_Vector x,N_Vector w,N_Vector id)363 realtype N_VWrmsNormMask(N_Vector x, N_Vector w, N_Vector id)
364 {
365 return((realtype) x->ops->nvwrmsnormmask(x, w, id));
366 }
367
N_VMin(N_Vector x)368 realtype N_VMin(N_Vector x)
369 {
370 return((realtype) x->ops->nvmin(x));
371 }
372
N_VWL2Norm(N_Vector x,N_Vector w)373 realtype N_VWL2Norm(N_Vector x, N_Vector w)
374 {
375 return((realtype) x->ops->nvwl2norm(x, w));
376 }
377
N_VL1Norm(N_Vector x)378 realtype N_VL1Norm(N_Vector x)
379 {
380 return((realtype) x->ops->nvl1norm(x));
381 }
382
N_VCompare(realtype c,N_Vector x,N_Vector z)383 void N_VCompare(realtype c, N_Vector x, N_Vector z)
384 {
385 z->ops->nvcompare(c, x, z);
386 return;
387 }
388
N_VInvTest(N_Vector x,N_Vector z)389 booleantype N_VInvTest(N_Vector x, N_Vector z)
390 {
391 return((booleantype) z->ops->nvinvtest(x, z));
392 }
393
N_VConstrMask(N_Vector c,N_Vector x,N_Vector m)394 booleantype N_VConstrMask(N_Vector c, N_Vector x, N_Vector m)
395 {
396 return((booleantype) x->ops->nvconstrmask(c, x, m));
397 }
398
N_VMinQuotient(N_Vector num,N_Vector denom)399 realtype N_VMinQuotient(N_Vector num, N_Vector denom)
400 {
401 return((realtype) num->ops->nvminquotient(num, denom));
402 }
403
404 /* -----------------------------------------------------------------
405 * OPTIONAL fused vector operations
406 * -----------------------------------------------------------------*/
407
N_VLinearCombination(int nvec,realtype * c,N_Vector * X,N_Vector z)408 int N_VLinearCombination(int nvec, realtype* c, N_Vector* X, N_Vector z)
409 {
410 int i;
411 realtype ONE=RCONST(1.0);
412
413 if (z->ops->nvlinearcombination != NULL) {
414
415 return(z->ops->nvlinearcombination(nvec, c, X, z));
416
417 } else {
418
419 z->ops->nvscale(c[0], X[0], z);
420 for (i=1; i<nvec; i++) {
421 z->ops->nvlinearsum(c[i], X[i], ONE, z, z);
422 }
423 return(0);
424
425 }
426 }
427
N_VScaleAddMulti(int nvec,realtype * a,N_Vector x,N_Vector * Y,N_Vector * Z)428 int N_VScaleAddMulti(int nvec, realtype* a, N_Vector x, N_Vector* Y, N_Vector* Z)
429 {
430 int i;
431 realtype ONE=RCONST(1.0);
432
433 if (x->ops->nvscaleaddmulti != NULL) {
434
435 return(x->ops->nvscaleaddmulti(nvec, a, x, Y, Z));
436
437 } else {
438
439 for (i=0; i<nvec; i++) {
440 x->ops->nvlinearsum(a[i], x, ONE, Y[i], Z[i]);
441 }
442 return(0);
443
444 }
445 }
446
N_VDotProdMulti(int nvec,N_Vector x,N_Vector * Y,realtype * dotprods)447 int N_VDotProdMulti(int nvec, N_Vector x, N_Vector* Y, realtype* dotprods)
448 {
449 int i;
450
451 if (x->ops->nvdotprodmulti != NULL) {
452
453 return(x->ops->nvdotprodmulti(nvec, x, Y, dotprods));
454
455 } else {
456
457 for (i=0; i<nvec; i++) {
458 dotprods[i] = x->ops->nvdotprod(x, Y[i]);
459 }
460 return(0);
461
462 }
463 }
464
465 /* -----------------------------------------------------------------
466 * OPTIONAL vector array operations
467 * -----------------------------------------------------------------*/
468
N_VLinearSumVectorArray(int nvec,realtype a,N_Vector * X,realtype b,N_Vector * Y,N_Vector * Z)469 int N_VLinearSumVectorArray(int nvec, realtype a, N_Vector* X,
470 realtype b, N_Vector* Y, N_Vector* Z)
471 {
472 int i;
473
474 if (Z[0]->ops->nvlinearsumvectorarray != NULL) {
475
476 return(Z[0]->ops->nvlinearsumvectorarray(nvec, a, X, b, Y, Z));
477
478 } else {
479
480 for (i=0; i<nvec; i++) {
481 Z[0]->ops->nvlinearsum(a, X[i], b, Y[i], Z[i]);
482 }
483 return(0);
484
485 }
486 }
487
N_VScaleVectorArray(int nvec,realtype * c,N_Vector * X,N_Vector * Z)488 int N_VScaleVectorArray(int nvec, realtype* c, N_Vector* X, N_Vector* Z)
489 {
490 int i;
491
492 if (Z[0]->ops->nvscalevectorarray != NULL) {
493
494 return(Z[0]->ops->nvscalevectorarray(nvec, c, X, Z));
495
496 } else {
497
498 for (i=0; i<nvec; i++) {
499 Z[0]->ops->nvscale(c[i], X[i], Z[i]);
500 }
501 return(0);
502
503 }
504 }
505
N_VConstVectorArray(int nvec,realtype c,N_Vector * Z)506 int N_VConstVectorArray(int nvec, realtype c, N_Vector* Z)
507 {
508 int i, ier;
509
510 if (Z[0]->ops->nvconstvectorarray != NULL) {
511
512 ier = Z[0]->ops->nvconstvectorarray(nvec, c, Z);
513 return(ier);
514
515 } else {
516
517 for (i=0; i<nvec; i++) {
518 Z[0]->ops->nvconst(c, Z[i]);
519 }
520 return(0);
521
522 }
523 }
524
N_VWrmsNormVectorArray(int nvec,N_Vector * X,N_Vector * W,realtype * nrm)525 int N_VWrmsNormVectorArray(int nvec, N_Vector* X, N_Vector* W, realtype* nrm)
526 {
527 int i, ier;
528
529 if (X[0]->ops->nvwrmsnormvectorarray != NULL) {
530
531 ier = X[0]->ops->nvwrmsnormvectorarray(nvec, X, W, nrm);
532 return(ier);
533
534 } else {
535
536 for (i=0; i<nvec; i++) {
537 nrm[i] = X[0]->ops->nvwrmsnorm(X[i], W[i]);
538 }
539 return(0);
540
541 }
542 }
543
N_VWrmsNormMaskVectorArray(int nvec,N_Vector * X,N_Vector * W,N_Vector id,realtype * nrm)544 int N_VWrmsNormMaskVectorArray(int nvec, N_Vector* X, N_Vector* W, N_Vector id,
545 realtype* nrm)
546 {
547 int i, ier;
548
549 if (id->ops->nvwrmsnormmaskvectorarray != NULL) {
550
551 ier = id->ops->nvwrmsnormmaskvectorarray(nvec, X, W, id, nrm);
552 return(ier);
553
554 } else {
555
556 for (i=0; i<nvec; i++) {
557 nrm[i] = id->ops->nvwrmsnormmask(X[i], W[i], id);
558 }
559 return(0);
560
561 }
562 }
563
N_VScaleAddMultiVectorArray(int nvec,int nsum,realtype * a,N_Vector * X,N_Vector ** Y,N_Vector ** Z)564 int N_VScaleAddMultiVectorArray(int nvec, int nsum, realtype* a, N_Vector* X,
565 N_Vector** Y, N_Vector** Z)
566 {
567 int i, j;
568 int ier=0;
569 realtype ONE=RCONST(1.0);
570 N_Vector* YY=NULL;
571 N_Vector* ZZ=NULL;
572
573 if (X[0]->ops->nvscaleaddmultivectorarray != NULL) {
574
575 ier = X[0]->ops->nvscaleaddmultivectorarray(nvec, nsum, a, X, Y, Z);
576 return(ier);
577
578 } else if (X[0]->ops->nvscaleaddmulti != NULL ) {
579
580 /* allocate arrays of vectors */
581 YY = (N_Vector*) malloc(nsum * sizeof(N_Vector));
582 ZZ = (N_Vector*) malloc(nsum * sizeof(N_Vector));
583
584 for (i=0; i<nvec; i++) {
585
586 for (j=0; j<nsum; j++) {
587 YY[j] = Y[j][i];
588 ZZ[j] = Z[j][i];
589 }
590
591 ier = X[0]->ops->nvscaleaddmulti(nsum, a, X[i], YY, ZZ);
592 if (ier != 0) break;
593 }
594
595 /* free array of vectors */
596 free(YY);
597 free(ZZ);
598
599 return(ier);
600
601 } else {
602
603 for (i=0; i<nvec; i++) {
604 for (j=0; j<nsum; j++) {
605 X[0]->ops->nvlinearsum(a[j], X[i], ONE, Y[j][i], Z[j][i]);
606 }
607 }
608 return(0);
609 }
610 }
611
N_VLinearCombinationVectorArray(int nvec,int nsum,realtype * c,N_Vector ** X,N_Vector * Z)612 int N_VLinearCombinationVectorArray(int nvec, int nsum, realtype* c,
613 N_Vector** X, N_Vector* Z)
614 {
615 int i, j;
616 int ier=0;
617 realtype ONE=RCONST(1.0);
618 N_Vector* Y=NULL;
619
620 if (Z[0]->ops->nvlinearcombinationvectorarray != NULL) {
621
622 ier = Z[0]->ops->nvlinearcombinationvectorarray(nvec, nsum, c, X, Z);
623 return(ier);
624
625 } else if (Z[0]->ops->nvlinearcombination != NULL ) {
626
627 /* allocate array of vectors */
628 Y = (N_Vector* ) malloc(nsum * sizeof(N_Vector));
629
630 for (i=0; i<nvec; i++) {
631
632 for (j=0; j<nsum; j++) {
633 Y[j] = X[j][i];
634 }
635
636 ier = Z[0]->ops->nvlinearcombination(nsum, c, Y, Z[i]);
637 if (ier != 0) break;
638 }
639
640 /* free array of vectors */
641 free(Y);
642
643 return(ier);
644
645 } else {
646
647 for (i=0; i<nvec; i++) {
648 Z[0]->ops->nvscale(c[0], X[0][i], Z[i]);
649 for (j=1; j<nsum; j++) {
650 Z[0]->ops->nvlinearsum(c[j], X[j][i], ONE, Z[i], Z[i]);
651 }
652 }
653 return(0);
654 }
655 }
656
657 /* -----------------------------------------------------------------
658 * OPTIONAL local reduction kernels (no parallel communication)
659 * -----------------------------------------------------------------*/
660
N_VDotProdLocal(N_Vector x,N_Vector y)661 realtype N_VDotProdLocal(N_Vector x, N_Vector y)
662 {
663 return((realtype) y->ops->nvdotprodlocal(x, y));
664 }
665
N_VMaxNormLocal(N_Vector x)666 realtype N_VMaxNormLocal(N_Vector x)
667 {
668 return((realtype) x->ops->nvmaxnormlocal(x));
669 }
670
N_VMinLocal(N_Vector x)671 realtype N_VMinLocal(N_Vector x)
672 {
673 return((realtype) x->ops->nvminlocal(x));
674 }
675
N_VL1NormLocal(N_Vector x)676 realtype N_VL1NormLocal(N_Vector x)
677 {
678 return((realtype) x->ops->nvl1normlocal(x));
679 }
680
N_VWSqrSumLocal(N_Vector x,N_Vector w)681 realtype N_VWSqrSumLocal(N_Vector x, N_Vector w)
682 {
683 return((realtype) x->ops->nvwsqrsumlocal(x,w));
684 }
685
N_VWSqrSumMaskLocal(N_Vector x,N_Vector w,N_Vector id)686 realtype N_VWSqrSumMaskLocal(N_Vector x, N_Vector w, N_Vector id)
687 {
688 return((realtype) x->ops->nvwsqrsummasklocal(x,w,id));
689 }
690
N_VInvTestLocal(N_Vector x,N_Vector z)691 booleantype N_VInvTestLocal(N_Vector x, N_Vector z)
692 {
693 return((booleantype) z->ops->nvinvtestlocal(x,z));
694 }
695
N_VConstrMaskLocal(N_Vector c,N_Vector x,N_Vector m)696 booleantype N_VConstrMaskLocal(N_Vector c, N_Vector x, N_Vector m)
697 {
698 return((booleantype) x->ops->nvconstrmasklocal(c,x,m));
699 }
700
N_VMinQuotientLocal(N_Vector num,N_Vector denom)701 realtype N_VMinQuotientLocal(N_Vector num, N_Vector denom)
702 {
703 return((realtype) num->ops->nvminquotientlocal(num,denom));
704 }
705
706 /* ------------------------------------
707 * OPTIONAL XBraid interface operations
708 * ------------------------------------*/
709
N_VBufSize(N_Vector x,sunindextype * size)710 int N_VBufSize(N_Vector x, sunindextype *size)
711 {
712 if (x->ops->nvbufsize == NULL) return(-1);
713 return(x->ops->nvbufsize(x, size));
714 }
715
N_VBufPack(N_Vector x,void * buf)716 int N_VBufPack(N_Vector x, void *buf)
717 {
718 if (x->ops->nvbufpack == NULL) return(-1);
719 return(x->ops->nvbufpack(x, buf));
720 }
721
N_VBufUnpack(N_Vector x,void * buf)722 int N_VBufUnpack(N_Vector x, void *buf)
723 {
724 if (x->ops->nvbufunpack == NULL) return(-1);
725 return(x->ops->nvbufunpack(x, buf));
726 }
727
728 /* -----------------------------------------------------------------
729 * Additional functions exported by the generic NVECTOR:
730 * N_VNewVectorArray
731 * N_VCloneEmptyVectorArray
732 * N_VCloneVectorArray
733 * N_VDestroyVectorArray
734 * -----------------------------------------------------------------*/
N_VNewVectorArray(int count)735 N_Vector* N_VNewVectorArray(int count)
736 {
737 N_Vector* vs = NULL;
738
739 if (count <= 0) return(NULL);
740
741 vs = (N_Vector* ) malloc(count * sizeof(N_Vector));
742 if(vs == NULL) return(NULL);
743
744 return(vs);
745 }
746
N_VCloneEmptyVectorArray(int count,N_Vector w)747 N_Vector* N_VCloneEmptyVectorArray(int count, N_Vector w)
748 {
749 N_Vector* vs = NULL;
750 int j;
751
752 if (count <= 0) return(NULL);
753
754 vs = (N_Vector* ) malloc(count * sizeof(N_Vector));
755 if(vs == NULL) return(NULL);
756
757 for (j = 0; j < count; j++) {
758 vs[j] = N_VCloneEmpty(w);
759 if (vs[j] == NULL) {
760 N_VDestroyVectorArray(vs, j-1);
761 return(NULL);
762 }
763 }
764
765 return(vs);
766 }
767
N_VCloneVectorArray(int count,N_Vector w)768 N_Vector* N_VCloneVectorArray(int count, N_Vector w)
769 {
770 N_Vector* vs = NULL;
771 int j;
772
773 if (count <= 0) return(NULL);
774
775 vs = (N_Vector* ) malloc(count * sizeof(N_Vector));
776 if(vs == NULL) return(NULL);
777
778 for (j = 0; j < count; j++) {
779 vs[j] = N_VClone(w);
780 if (vs[j] == NULL) {
781 N_VDestroyVectorArray(vs, j-1);
782 return(NULL);
783 }
784 }
785
786 return(vs);
787 }
788
N_VDestroyVectorArray(N_Vector * vs,int count)789 void N_VDestroyVectorArray(N_Vector* vs, int count)
790 {
791 int j;
792
793 if (vs==NULL) return;
794
795 for (j = 0; j < count; j++) N_VDestroy(vs[j]);
796
797 free(vs); vs = NULL;
798
799 return;
800 }
801
802 /* These function are really only for users of the Fortran interface */
N_VGetVecAtIndexVectorArray(N_Vector * vs,int index)803 N_Vector N_VGetVecAtIndexVectorArray(N_Vector* vs, int index)
804 {
805 if (vs==NULL) return NULL;
806 else if (index < 0) return NULL;
807 else return vs[index];
808 }
809
N_VSetVecAtIndexVectorArray(N_Vector * vs,int index,N_Vector w)810 void N_VSetVecAtIndexVectorArray(N_Vector* vs, int index, N_Vector w)
811 {
812 if (vs==NULL) return;
813 else if (index < 0) return;
814 else vs[index] = w;
815 }
816
817
818 /* -----------------------------------------------------------------
819 * Debugging functions
820 * ----------------------------------------------------------------- */
821
N_VPrint(N_Vector v)822 void N_VPrint(N_Vector v)
823 {
824 if (v == NULL) {
825 printf("NULL Vector\n");
826 return;
827 }
828 if (v->ops->nvprint == NULL) {
829 printf("NULL Print Op\n");
830 return;
831 }
832 v->ops->nvprint(v);
833 return;
834 }
835
836
N_VPrintFile(N_Vector v,FILE * outfile)837 void N_VPrintFile(N_Vector v, FILE* outfile)
838 {
839 if (v == NULL) {
840 fprintf(outfile, "NULL Vector\n");
841 return;
842 }
843 if (v->ops->nvprintfile == NULL) {
844 fprintf(outfile, "NULL PrintFile Op\n");
845 return;
846 }
847 v->ops->nvprintfile(v, outfile);
848 return;
849 }
850