1 /* -----------------------------------------------------------------
2 * Programmer(s): Slaven Peles @ LLNL
3 * -----------------------------------------------------------------
4 * Based on N_Vector_Parallel by Scott D. Cohen, Alan C. Hindmarsh,
5 * Radu Serban, and Aaron Collier @ LLNL
6 * -----------------------------------------------------------------
7 * SUNDIALS Copyright Start
8 * Copyright (c) 2002-2019, Lawrence Livermore National Security
9 * and Southern Methodist University.
10 * All rights reserved.
11 *
12 * See the top-level LICENSE and NOTICE files for details.
13 *
14 * SPDX-License-Identifier: BSD-3-Clause
15 * SUNDIALS Copyright End
16 * -----------------------------------------------------------------
17 * This is the implementation file for a PETSc implementation
18 * of the NVECTOR package.
19 * -----------------------------------------------------------------*/
20
21 #include <stdio.h>
22 #include <stdlib.h>
23
24 #include <nvector/nvector_petsc.h>
25 #include <sundials/sundials_math.h>
26 #include <sundials/sundials_mpi.h>
27
28 #define ZERO RCONST(0.0)
29 #define HALF RCONST(0.5)
30 #define ONE RCONST(1.0)
31 #define ONEPT5 RCONST(1.5)
32
33 /* Error Message */
34 #define BAD_N1 "N_VNewEmpty_Petsc -- Sum of local vector lengths differs from "
35 #define BAD_N2 "input global length. \n\n"
36 #define BAD_N BAD_N1 BAD_N2
37
38 /*
39 * -----------------------------------------------------------------
40 * Simplifying macros NV_CONTENT_PTC, NV_OWN_DATA_PTC,
41 * NV_LOCLENGTH_PTC, NV_GLOBLENGTH_PTC,
42 * NV_COMM_PTC
43 * -----------------------------------------------------------------
44 * In the descriptions below, the following user declarations
45 * are assumed:
46 *
47 * N_Vector v;
48 * sunindextype v_len, s_len, i;
49 *
50 * (1) NV_CONTENT_PTC
51 *
52 * This routines gives access to the contents of the PETSc
53 * vector wrapper N_Vector.
54 *
55 * The assignment v_cont = NV_CONTENT_PTC(v) sets v_cont to be
56 * a pointer to the N_Vector (PETSc wrapper) content structure.
57 *
58 * (2) NV_PVEC_PTC, NV_OWN_DATA_PTC, NV_LOCLENGTH_PTC, NV_GLOBLENGTH_PTC,
59 * and NV_COMM_PTC
60 *
61 * These routines give access to the individual parts of
62 * the content structure of a PETSc N_Vector wrapper.
63 *
64 * NV_PVEC_PTC(v) returns the PETSc vector (Vec) object.
65 *
66 * The assignment v_llen = NV_LOCLENGTH_PTC(v) sets v_llen to
67 * be the length of the local part of the vector v. The call
68 * NV_LOCLENGTH_PTC(v) = llen_v should NOT be used! It will
69 * change the value stored in the N_Vector content structure,
70 * but it will NOT change the length of the actual PETSc vector.
71 *
72 * The assignment v_glen = NV_GLOBLENGTH_PTC(v) sets v_glen to
73 * be the global length of the vector v. The call
74 * NV_GLOBLENGTH_PTC(v) = glen_v should NOT be used! It will
75 * change the value stored in the N_Vector content structure,
76 * but it will NOT change the length of the actual PETSc vector.
77 *
78 * The assignment v_comm = NV_COMM_PTC(v) sets v_comm to be the
79 * MPI communicator of the vector v. The assignment
80 * NV_COMM_PTC(v) = comm_v should NOT be used! It will change
81 * the value stored in the N_Vector content structure, but it
82 * will NOT change the MPI communicator of the actual PETSc
83 * vector.
84 *
85 * -----------------------------------------------------------------
86 */
87
88 #define NV_CONTENT_PTC(v) ( (N_VectorContent_Petsc)(v->content) )
89
90 #define NV_LOCLENGTH_PTC(v) ( NV_CONTENT_PTC(v)->local_length )
91
92 #define NV_GLOBLENGTH_PTC(v) ( NV_CONTENT_PTC(v)->global_length )
93
94 #define NV_OWN_DATA_PTC(v) ( NV_CONTENT_PTC(v)->own_data )
95
96 #define NV_PVEC_PTC(v) ( NV_CONTENT_PTC(v)->pvec )
97
98 #define NV_COMM_PTC(v) ( NV_CONTENT_PTC(v)->comm )
99
100
101
102 /* ----------------------------------------------------------------
103 * Returns vector type ID. Used to identify vector implementation
104 * from abstract N_Vector interface.
105 */
N_VGetVectorID_Petsc(N_Vector v)106 N_Vector_ID N_VGetVectorID_Petsc(N_Vector v)
107 {
108 return SUNDIALS_NVEC_PETSC;
109 }
110
111
112 /* ----------------------------------------------------------------
113 * Function to create a new N_Vector wrapper with an empty (NULL)
114 * PETSc vector.
115 */
116
N_VNewEmpty_Petsc(MPI_Comm comm,sunindextype local_length,sunindextype global_length)117 N_Vector N_VNewEmpty_Petsc(MPI_Comm comm,
118 sunindextype local_length,
119 sunindextype global_length)
120 {
121 N_Vector v;
122 N_Vector_Ops ops;
123 N_VectorContent_Petsc content;
124 sunindextype n, Nsum;
125 PetscErrorCode ierr;
126
127 /* Compute global length as sum of local lengths */
128 n = local_length;
129 ierr = MPI_Allreduce(&n, &Nsum, 1, PVEC_INTEGER_MPI_TYPE, MPI_SUM, comm);
130 CHKERRABORT(comm,ierr);
131 if (Nsum != global_length) {
132 STAN_SUNDIALS_FPRINTF(stderr, BAD_N);
133 return(NULL);
134 }
135
136 /* Create vector */
137 v = NULL;
138 v = (N_Vector) malloc(sizeof *v);
139 if (v == NULL) return(NULL);
140
141 /* Create vector operation structure */
142 ops = NULL;
143 ops = (N_Vector_Ops) malloc(sizeof(struct _generic_N_Vector_Ops));
144 if (ops == NULL) { free(v); return(NULL); }
145
146 ops->nvgetvectorid = N_VGetVectorID_Petsc;
147 ops->nvclone = N_VClone_Petsc;
148 ops->nvcloneempty = N_VCloneEmpty_Petsc;
149 ops->nvdestroy = N_VDestroy_Petsc;
150 ops->nvspace = N_VSpace_Petsc;
151 ops->nvgetarraypointer = N_VGetArrayPointer_Petsc;
152 ops->nvsetarraypointer = N_VSetArrayPointer_Petsc;
153
154 /* standard vector operations */
155 ops->nvlinearsum = N_VLinearSum_Petsc;
156 ops->nvconst = N_VConst_Petsc;
157 ops->nvprod = N_VProd_Petsc;
158 ops->nvdiv = N_VDiv_Petsc;
159 ops->nvscale = N_VScale_Petsc;
160 ops->nvabs = N_VAbs_Petsc;
161 ops->nvinv = N_VInv_Petsc;
162 ops->nvaddconst = N_VAddConst_Petsc;
163 ops->nvdotprod = N_VDotProd_Petsc;
164 ops->nvmaxnorm = N_VMaxNorm_Petsc;
165 ops->nvwrmsnormmask = N_VWrmsNormMask_Petsc;
166 ops->nvwrmsnorm = N_VWrmsNorm_Petsc;
167 ops->nvmin = N_VMin_Petsc;
168 ops->nvwl2norm = N_VWL2Norm_Petsc;
169 ops->nvl1norm = N_VL1Norm_Petsc;
170 ops->nvcompare = N_VCompare_Petsc;
171 ops->nvinvtest = N_VInvTest_Petsc;
172 ops->nvconstrmask = N_VConstrMask_Petsc;
173 ops->nvminquotient = N_VMinQuotient_Petsc;
174
175 /* fused vector operations (optional, NULL means disabled by default) */
176 ops->nvlinearcombination = NULL;
177 ops->nvscaleaddmulti = NULL;
178 ops->nvdotprodmulti = NULL;
179
180 /* vector array operations (optional, NULL means disabled by default) */
181 ops->nvlinearsumvectorarray = NULL;
182 ops->nvscalevectorarray = NULL;
183 ops->nvconstvectorarray = NULL;
184 ops->nvwrmsnormvectorarray = NULL;
185 ops->nvwrmsnormmaskvectorarray = NULL;
186 ops->nvscaleaddmultivectorarray = NULL;
187 ops->nvlinearcombinationvectorarray = NULL;
188
189 /* Create content */
190 content = NULL;
191 content = (N_VectorContent_Petsc) malloc(sizeof(struct _N_VectorContent_Petsc));
192 if (content == NULL) {
193 free(ops);
194 free(v);
195 return(NULL);
196 }
197
198 /* Attach lengths and communicator */
199 content->local_length = local_length;
200 content->global_length = global_length;
201 content->comm = comm;
202 content->own_data = SUNFALSE;
203 content->pvec = NULL;
204
205 /* Attach content and ops */
206 v->content = content;
207 v->ops = ops;
208
209 return(v);
210 }
211
212
213
214 /* ----------------------------------------------------------------
215 * Function to create an N_Vector wrapper for a PETSc vector.
216 */
217
N_VMake_Petsc(Vec pvec)218 N_Vector N_VMake_Petsc(Vec pvec)
219 {
220 N_Vector v = NULL;
221 MPI_Comm comm;
222 PetscInt local_length;
223 PetscInt global_length;
224
225 VecGetLocalSize(pvec, &local_length);
226 VecGetSize(pvec, &global_length);
227 PetscObjectGetComm((PetscObject) pvec, &comm);
228
229 v = N_VNewEmpty_Petsc(comm, local_length, global_length);
230 if (v == NULL)
231 return(NULL);
232
233 /* Attach data */
234 NV_OWN_DATA_PTC(v) = SUNFALSE;
235 NV_PVEC_PTC(v) = pvec;
236
237 return(v);
238 }
239
240 /* ----------------------------------------------------------------
241 * Function to create an array of new PETSc vector wrappers.
242 */
243
N_VCloneVectorArray_Petsc(int count,N_Vector w)244 N_Vector *N_VCloneVectorArray_Petsc(int count, N_Vector w)
245 {
246 N_Vector *vs;
247 int j;
248
249 if (count <= 0) return(NULL);
250
251 vs = NULL;
252 vs = (N_Vector *) malloc(count * sizeof(N_Vector));
253 if(vs == NULL) return(NULL);
254
255 for (j = 0; j < count; j++) {
256 vs[j] = NULL;
257 vs[j] = N_VClone_Petsc(w);
258 if (vs[j] == NULL) {
259 N_VDestroyVectorArray_Petsc(vs, j-1);
260 return(NULL);
261 }
262 }
263
264 return(vs);
265 }
266
267 /* ----------------------------------------------------------------
268 * Function to create an array of new PETSc vector wrappers with
269 * empty (NULL) PETSc vectors.
270 */
271
N_VCloneVectorArrayEmpty_Petsc(int count,N_Vector w)272 N_Vector *N_VCloneVectorArrayEmpty_Petsc(int count, N_Vector w)
273 {
274 N_Vector *vs;
275 int j;
276
277 if (count <= 0) return(NULL);
278
279 vs = NULL;
280 vs = (N_Vector *) malloc(count * sizeof(N_Vector));
281 if(vs == NULL) return(NULL);
282
283 for (j = 0; j < count; j++) {
284 vs[j] = NULL;
285 vs[j] = N_VCloneEmpty_Petsc(w);
286 if (vs[j] == NULL) {
287 N_VDestroyVectorArray_Petsc(vs, j-1);
288 return(NULL);
289 }
290 }
291
292 return(vs);
293 }
294
295 /* ----------------------------------------------------------------
296 * Function to free an array created with N_VCloneVectorArray_Petsc
297 */
298
N_VDestroyVectorArray_Petsc(N_Vector * vs,int count)299 void N_VDestroyVectorArray_Petsc(N_Vector *vs, int count)
300 {
301 int j;
302
303 for (j = 0; j < count; j++) N_VDestroy_Petsc(vs[j]);
304
305 free(vs);
306 vs = NULL;
307
308 return;
309 }
310
311 /* ----------------------------------------------------------------
312 * Function to extract PETSc vector
313 */
314
N_VGetVector_Petsc(N_Vector v)315 Vec N_VGetVector_Petsc(N_Vector v)
316 {
317 return NV_PVEC_PTC(v);
318 }
319
320 /* ----------------------------------------------------------------
321 * Function to print the global data in a PETSc vector to stdout
322 */
323
N_VPrint_Petsc(N_Vector x)324 void N_VPrint_Petsc(N_Vector x)
325 {
326 Vec xv = NV_PVEC_PTC(x);
327 MPI_Comm comm = NV_COMM_PTC(x);
328
329 VecView(xv, PETSC_VIEWER_STDOUT_(comm));
330
331 return;
332 }
333
334 /* ----------------------------------------------------------------
335 * Function to print the global data in a PETSc vector to fname
336 */
337
N_VPrintFile_Petsc(N_Vector x,const char fname[])338 void N_VPrintFile_Petsc(N_Vector x, const char fname[])
339 {
340 Vec xv = NV_PVEC_PTC(x);
341 MPI_Comm comm = NV_COMM_PTC(x);
342 PetscViewer viewer;
343
344 PetscViewerASCIIOpen(comm, fname, &viewer);
345
346 VecView(xv, viewer);
347
348 PetscViewerDestroy(&viewer);
349
350 return;
351 }
352
353 /*
354 * -----------------------------------------------------------------
355 * implementation of vector operations
356 * -----------------------------------------------------------------
357 */
358
N_VCloneEmpty_Petsc(N_Vector w)359 N_Vector N_VCloneEmpty_Petsc(N_Vector w)
360 {
361 N_Vector v;
362 N_Vector_Ops ops;
363 N_VectorContent_Petsc content;
364
365 if (w == NULL) return(NULL);
366
367 /* Create vector */
368 v = NULL;
369 v = (N_Vector) malloc(sizeof *v);
370 if (v == NULL) return(NULL);
371
372 /* Create vector operation structure */
373 ops = NULL;
374 ops = (N_Vector_Ops) malloc(sizeof(struct _generic_N_Vector_Ops));
375 if (ops == NULL) {
376 free(v);
377 return(NULL);
378 }
379
380 ops->nvgetvectorid = w->ops->nvgetvectorid;
381 ops->nvclone = w->ops->nvclone;
382 ops->nvcloneempty = w->ops->nvcloneempty;
383 ops->nvdestroy = w->ops->nvdestroy;
384 ops->nvspace = w->ops->nvspace;
385 ops->nvgetarraypointer = w->ops->nvgetarraypointer;
386 ops->nvsetarraypointer = w->ops->nvsetarraypointer;
387
388 /* standard vector operations */
389 ops->nvlinearsum = w->ops->nvlinearsum;
390 ops->nvconst = w->ops->nvconst;
391 ops->nvprod = w->ops->nvprod;
392 ops->nvdiv = w->ops->nvdiv;
393 ops->nvscale = w->ops->nvscale;
394 ops->nvabs = w->ops->nvabs;
395 ops->nvinv = w->ops->nvinv;
396 ops->nvaddconst = w->ops->nvaddconst;
397 ops->nvdotprod = w->ops->nvdotprod;
398 ops->nvmaxnorm = w->ops->nvmaxnorm;
399 ops->nvwrmsnormmask = w->ops->nvwrmsnormmask;
400 ops->nvwrmsnorm = w->ops->nvwrmsnorm;
401 ops->nvmin = w->ops->nvmin;
402 ops->nvwl2norm = w->ops->nvwl2norm;
403 ops->nvl1norm = w->ops->nvl1norm;
404 ops->nvcompare = w->ops->nvcompare;
405 ops->nvinvtest = w->ops->nvinvtest;
406 ops->nvconstrmask = w->ops->nvconstrmask;
407 ops->nvminquotient = w->ops->nvminquotient;
408
409 /* fused vector operations */
410 ops->nvlinearcombination = w->ops->nvlinearcombination;
411 ops->nvscaleaddmulti = w->ops->nvscaleaddmulti;
412 ops->nvdotprodmulti = w->ops->nvdotprodmulti;
413
414 /* vector array operations */
415 ops->nvlinearsumvectorarray = w->ops->nvlinearsumvectorarray;
416 ops->nvscalevectorarray = w->ops->nvscalevectorarray;
417 ops->nvconstvectorarray = w->ops->nvconstvectorarray;
418 ops->nvwrmsnormvectorarray = w->ops->nvwrmsnormvectorarray;
419 ops->nvwrmsnormmaskvectorarray = w->ops->nvwrmsnormmaskvectorarray;
420 ops->nvscaleaddmultivectorarray = w->ops->nvscaleaddmultivectorarray;
421 ops->nvlinearcombinationvectorarray = w->ops->nvlinearcombinationvectorarray;
422
423 /* Create content */
424 content = NULL;
425 content = (N_VectorContent_Petsc) malloc(sizeof(struct _N_VectorContent_Petsc));
426 if (content == NULL) {
427 free(ops);
428 free(v);
429 return(NULL);
430 }
431
432 /* Attach lengths and communicator */
433 content->local_length = NV_LOCLENGTH_PTC(w);
434 content->global_length = NV_GLOBLENGTH_PTC(w);
435 content->comm = NV_COMM_PTC(w);
436 content->own_data = SUNFALSE;
437 content->pvec = NULL;
438
439 /* Attach content and ops */
440 v->content = content;
441 v->ops = ops;
442
443 return(v);
444 }
445
N_VClone_Petsc(N_Vector w)446 N_Vector N_VClone_Petsc(N_Vector w)
447 {
448 N_Vector v = NULL;
449 Vec pvec = NULL;
450 Vec wvec = NV_PVEC_PTC(w);
451
452 /* PetscErrorCode ierr; */
453
454 v = N_VCloneEmpty_Petsc(w);
455 if (v == NULL)
456 return(NULL);
457
458 /* Create data */
459
460 /* Allocate empty PETSc vector */
461 pvec = (Vec) malloc(sizeof(Vec));
462 if(pvec == NULL) {
463 N_VDestroy_Petsc(v);
464 return(NULL);
465 }
466
467 /* ierr = */
468 VecDuplicate(wvec, &pvec);
469 if(pvec == NULL) {
470 N_VDestroy_Petsc(v);
471 return(NULL);
472 }
473
474 /* Attach data */
475 NV_OWN_DATA_PTC(v) = SUNTRUE;
476 NV_PVEC_PTC(v) = pvec;
477
478 return(v);
479 }
480
N_VDestroy_Petsc(N_Vector v)481 void N_VDestroy_Petsc(N_Vector v)
482 {
483 if (NV_OWN_DATA_PTC(v) == SUNTRUE) {
484 VecDestroy(&(NV_PVEC_PTC(v)));
485 NV_PVEC_PTC(v) = NULL;
486 }
487
488 free(v->content);
489 v->content = NULL;
490 free(v->ops);
491 v->ops = NULL;
492 free(v);
493 v = NULL;
494
495 return;
496 }
497
N_VSpace_Petsc(N_Vector v,sunindextype * lrw,sunindextype * liw)498 void N_VSpace_Petsc(N_Vector v, sunindextype *lrw, sunindextype *liw)
499 {
500 MPI_Comm comm;
501 int npes;
502
503 comm = NV_COMM_PTC(v);
504 MPI_Comm_size(comm, &npes);
505
506 *lrw = NV_GLOBLENGTH_PTC(v);
507 *liw = 2*npes;
508
509 return;
510 }
511
512 /*
513 * Not implemented for PETSc wrapper.
514 */
N_VGetArrayPointer_Petsc(N_Vector v)515 realtype *N_VGetArrayPointer_Petsc(N_Vector v)
516 {
517 return NULL;
518 }
519
520 /*
521 * Not implemented for PETSc wrapper.
522 */
N_VSetArrayPointer_Petsc(realtype * v_data,N_Vector v)523 void N_VSetArrayPointer_Petsc(realtype *v_data, N_Vector v)
524 {
525 return;
526 }
527
N_VLinearSum_Petsc(realtype a,N_Vector x,realtype b,N_Vector y,N_Vector z)528 void N_VLinearSum_Petsc(realtype a, N_Vector x, realtype b, N_Vector y, N_Vector z)
529 {
530 Vec xv = NV_PVEC_PTC(x);
531 Vec yv = NV_PVEC_PTC(y);
532 Vec zv = NV_PVEC_PTC(z);
533
534 if (x == y) {
535 N_VScale_Petsc(a + b, x, z); /* z <~ ax+bx */
536 return;
537 }
538
539 if (z == y) {
540 if (b == ONE) {
541 VecAXPY(yv, a, xv); /* BLAS usage: axpy y <- ax+y */
542 return;
543 }
544 VecAXPBY(yv, a, b, xv); /* BLAS usage: axpby y <- ax+by */
545 return;
546 }
547
548 if (z == x) {
549 if (a == ONE) {
550 VecAXPY(xv, b, yv); /* BLAS usage: axpy x <- by+x */
551 return;
552 }
553 VecAXPBY(xv, b, a, yv); /* BLAS usage: axpby x <- by+ax */
554 return;
555 }
556
557
558 /* Do all cases not handled above:
559 (1) a == other, b == 0.0 - user should have called N_VScale
560 (2) a == 0.0, b == other - user should have called N_VScale
561 (3) a,b == other, a !=b, a != -b */
562
563 VecAXPBYPCZ(zv, a, b, 0.0, xv, yv); /* PETSc, probably not optimal */
564
565 return;
566 }
567
N_VConst_Petsc(realtype c,N_Vector z)568 void N_VConst_Petsc(realtype c, N_Vector z)
569 {
570 Vec zv = NV_PVEC_PTC(z);
571
572 VecSet(zv, c);
573
574 return;
575 }
576
N_VProd_Petsc(N_Vector x,N_Vector y,N_Vector z)577 void N_VProd_Petsc(N_Vector x, N_Vector y, N_Vector z)
578 {
579 Vec xv = NV_PVEC_PTC(x);
580 Vec yv = NV_PVEC_PTC(y);
581 Vec zv = NV_PVEC_PTC(z);
582
583 VecPointwiseMult(zv, xv, yv);
584
585 return;
586 }
587
N_VDiv_Petsc(N_Vector x,N_Vector y,N_Vector z)588 void N_VDiv_Petsc(N_Vector x, N_Vector y, N_Vector z)
589 {
590 Vec xv = NV_PVEC_PTC(x);
591 Vec yv = NV_PVEC_PTC(y);
592 Vec zv = NV_PVEC_PTC(z);
593
594 VecPointwiseDivide(zv, xv, yv); /* z = x/y */
595
596 return;
597 }
598
N_VScale_Petsc(realtype c,N_Vector x,N_Vector z)599 void N_VScale_Petsc(realtype c, N_Vector x, N_Vector z)
600 {
601 Vec xv = NV_PVEC_PTC(x);
602 Vec zv = NV_PVEC_PTC(z);
603
604 if (z == x) { /* BLAS usage: scale x <- cx */
605 VecScale(xv, c);
606 return;
607 }
608
609 VecAXPBY(zv, c, 0.0, xv);
610
611 return;
612 }
613
N_VAbs_Petsc(N_Vector x,N_Vector z)614 void N_VAbs_Petsc(N_Vector x, N_Vector z)
615 {
616 Vec xv = NV_PVEC_PTC(x);
617 Vec zv = NV_PVEC_PTC(z);
618
619 if(z != x)
620 VecCopy(xv, zv); /* copy x~>z */
621 VecAbs(zv);
622
623 return;
624 }
625
N_VInv_Petsc(N_Vector x,N_Vector z)626 void N_VInv_Petsc(N_Vector x, N_Vector z)
627 {
628 Vec xv = NV_PVEC_PTC(x);
629 Vec zv = NV_PVEC_PTC(z);
630
631 if(z != x)
632 VecCopy(xv, zv); /* copy x~>z */
633 VecReciprocal(zv);
634
635 return;
636 }
637
N_VAddConst_Petsc(N_Vector x,realtype b,N_Vector z)638 void N_VAddConst_Petsc(N_Vector x, realtype b, N_Vector z)
639 {
640 Vec xv = NV_PVEC_PTC(x);
641 Vec zv = NV_PVEC_PTC(z);
642
643 if(z != x)
644 VecCopy(xv, zv); /* copy x~>z */
645 VecShift(zv, b);
646
647 return;
648 }
649
N_VDotProd_Petsc(N_Vector x,N_Vector y)650 realtype N_VDotProd_Petsc(N_Vector x, N_Vector y)
651 {
652 Vec xv = NV_PVEC_PTC(x);
653 Vec yv = NV_PVEC_PTC(y);
654 PetscScalar dotprod;
655
656 VecDot(xv, yv, &dotprod);
657
658 return dotprod;
659 }
660
N_VMaxNorm_Petsc(N_Vector x)661 realtype N_VMaxNorm_Petsc(N_Vector x)
662 {
663 Vec xv = NV_PVEC_PTC(x);
664 PetscReal norm;
665
666 VecNorm(xv, NORM_INFINITY, &norm);
667
668 return norm;
669 }
670
N_VWrmsNorm_Petsc(N_Vector x,N_Vector w)671 realtype N_VWrmsNorm_Petsc(N_Vector x, N_Vector w)
672 {
673 sunindextype i;
674 sunindextype N = NV_LOCLENGTH_PTC(x);
675 sunindextype N_global = NV_GLOBLENGTH_PTC(x);
676 MPI_Comm comm = NV_COMM_PTC(x);
677 Vec xv = NV_PVEC_PTC(x);
678 Vec wv = NV_PVEC_PTC(w);
679 PetscScalar *xd;
680 PetscScalar *wd;
681 PetscReal sum = ZERO;
682 realtype global_sum;
683
684 VecGetArray(xv, &xd);
685 VecGetArray(wv, &wd);
686 for (i = 0; i < N; i++) {
687 sum += PetscSqr(PetscAbsScalar(xd[i] * wd[i]));
688 }
689 VecRestoreArray(xv, &xd);
690 VecRestoreArray(wv, &wd);
691
692 global_sum = SUNMPI_Allreduce_scalar(sum, 1, comm);
693 return (SUNRsqrt(global_sum/N_global));
694 }
695
N_VWrmsNormMask_Petsc(N_Vector x,N_Vector w,N_Vector id)696 realtype N_VWrmsNormMask_Petsc(N_Vector x, N_Vector w, N_Vector id)
697 {
698 sunindextype i;
699 sunindextype N = NV_LOCLENGTH_PTC(x);
700 sunindextype N_global = NV_GLOBLENGTH_PTC(x);
701 MPI_Comm comm = NV_COMM_PTC(x);
702
703 Vec xv = NV_PVEC_PTC(x);
704 Vec wv = NV_PVEC_PTC(w);
705 Vec idv = NV_PVEC_PTC(id);
706 PetscScalar *xd;
707 PetscScalar *wd;
708 PetscScalar *idd;
709 PetscReal sum = ZERO;
710 realtype global_sum;
711
712 VecGetArray(xv, &xd);
713 VecGetArray(wv, &wd);
714 VecGetArray(idv, &idd);
715 for (i = 0; i < N; i++) {
716 PetscReal tag = (PetscReal) idd[i];
717 if (tag > ZERO) {
718 sum += PetscSqr(PetscAbsScalar(xd[i] * wd[i]));
719 }
720 }
721 VecRestoreArray(xv, &xd);
722 VecRestoreArray(wv, &wd);
723 VecRestoreArray(idv, &idd);
724
725 global_sum = SUNMPI_Allreduce_scalar(sum, 1, comm);
726 return (SUNRsqrt(global_sum/N_global));
727 }
728
N_VMin_Petsc(N_Vector x)729 realtype N_VMin_Petsc(N_Vector x)
730 {
731 Vec xv = NV_PVEC_PTC(x);
732 PetscReal minval;
733 PetscInt i;
734
735 VecMin(xv, &i, &minval);
736
737 return minval;
738 }
739
N_VWL2Norm_Petsc(N_Vector x,N_Vector w)740 realtype N_VWL2Norm_Petsc(N_Vector x, N_Vector w)
741 {
742 sunindextype i;
743 sunindextype N = NV_LOCLENGTH_PTC(x);
744 MPI_Comm comm = NV_COMM_PTC(x);
745
746 Vec xv = NV_PVEC_PTC(x);
747 Vec wv = NV_PVEC_PTC(w);
748 PetscScalar *xd;
749 PetscScalar *wd;
750 PetscReal sum = ZERO;
751 realtype global_sum;
752
753 VecGetArray(xv, &xd);
754 VecGetArray(wv, &wd);
755 for (i = 0; i < N; i++) {
756 sum += PetscSqr(PetscAbsScalar(xd[i] * wd[i]));
757 }
758 VecRestoreArray(xv, &xd);
759 VecRestoreArray(wv, &wd);
760
761 global_sum = SUNMPI_Allreduce_scalar(sum, 1, comm);
762 return (SUNRsqrt(global_sum));
763 }
764
N_VL1Norm_Petsc(N_Vector x)765 realtype N_VL1Norm_Petsc(N_Vector x)
766 {
767 Vec xv = NV_PVEC_PTC(x);
768 PetscReal norm;
769
770 VecNorm(xv, NORM_1, &norm);
771
772 return norm;
773 }
774
N_VCompare_Petsc(realtype c,N_Vector x,N_Vector z)775 void N_VCompare_Petsc(realtype c, N_Vector x, N_Vector z)
776 {
777 sunindextype i;
778 sunindextype N = NV_LOCLENGTH_PTC(x);
779 Vec xv = NV_PVEC_PTC(x);
780 Vec zv = NV_PVEC_PTC(z);
781 PetscReal cpet = c; /* <~ realtype should typedef to PETScReal */
782 PetscScalar *xdata;
783 PetscScalar *zdata;
784
785 VecGetArray(xv, &xdata);
786 VecGetArray(zv, &zdata);
787 for (i = 0; i < N; i++) {
788 zdata[i] = PetscAbsScalar(xdata[i]) >= cpet ? ONE : ZERO;
789 }
790 VecRestoreArray(xv, &xdata);
791 VecRestoreArray(zv, &zdata);
792
793 return;
794 }
795
N_VInvTest_Petsc(N_Vector x,N_Vector z)796 booleantype N_VInvTest_Petsc(N_Vector x, N_Vector z)
797 {
798 sunindextype i;
799 sunindextype N = NV_LOCLENGTH_PTC(x);
800 MPI_Comm comm = NV_COMM_PTC(x);
801 Vec xv = NV_PVEC_PTC(x);
802 Vec zv = NV_PVEC_PTC(z);
803 PetscScalar *xd;
804 PetscScalar *zd;
805 PetscReal val = ONE;
806
807 VecGetArray(xv, &xd);
808 VecGetArray(zv, &zd);
809 for (i = 0; i < N; i++) {
810 if (xd[i] == ZERO)
811 val = ZERO;
812 else
813 zd[i] = ONE/xd[i];
814 }
815 VecRestoreArray(xv, &xd);
816 VecRestoreArray(zv, &zd);
817
818 val = SUNMPI_Allreduce_scalar(val, 3, comm);
819
820 if (val == ZERO)
821 return(SUNFALSE);
822 else
823 return(SUNTRUE);
824 }
825
N_VConstrMask_Petsc(N_Vector c,N_Vector x,N_Vector m)826 booleantype N_VConstrMask_Petsc(N_Vector c, N_Vector x, N_Vector m)
827 {
828 sunindextype i;
829 sunindextype N = NV_LOCLENGTH_PTC(x);
830 MPI_Comm comm = NV_COMM_PTC(x);
831 realtype temp;
832 booleantype test;
833 Vec xv = NV_PVEC_PTC(x);
834 Vec cv = NV_PVEC_PTC(c);
835 Vec mv = NV_PVEC_PTC(m);
836 PetscScalar *xd;
837 PetscScalar *cd;
838 PetscScalar *md;
839
840 temp = ZERO;
841
842 VecGetArray(xv, &xd);
843 VecGetArray(cv, &cd);
844 VecGetArray(mv, &md);
845 for (i = 0; i < N; i++) {
846 PetscReal cc = (PetscReal) cd[i]; /* <~ Drop imaginary parts if any. */
847 PetscReal xx = (PetscReal) xd[i]; /* <~ Constraints defined on Re{x} */
848 md[i] = ZERO;
849
850 /* Continue if no constraints were set for the variable */
851 if (cc == ZERO)
852 continue;
853
854 /* Check if a set constraint has been violated */
855 test = (SUNRabs(cc) > ONEPT5 && xx*cc <= ZERO) ||
856 (SUNRabs(cc) > HALF && xx*cc < ZERO);
857 if (test) {
858 temp = md[i] = ONE;
859 }
860 }
861 VecRestoreArray(xv, &xd);
862 VecRestoreArray(cv, &cd);
863 VecRestoreArray(mv, &md);
864
865 /* Find max temp across all MPI ranks */
866 temp = SUNMPI_Allreduce_scalar(temp, 2, comm);
867
868 /* Return false if any constraint was violated */
869 return (temp == ONE) ? SUNFALSE : SUNTRUE;
870 }
871
N_VMinQuotient_Petsc(N_Vector num,N_Vector denom)872 realtype N_VMinQuotient_Petsc(N_Vector num, N_Vector denom)
873 {
874 booleantype notEvenOnce = SUNTRUE;
875 sunindextype i;
876 sunindextype N = NV_LOCLENGTH_PTC(num);
877 MPI_Comm comm = NV_COMM_PTC(num);
878
879 Vec nv = NV_PVEC_PTC(num);
880 Vec dv = NV_PVEC_PTC(denom);
881 PetscScalar *nd;
882 PetscScalar *dd;
883 PetscReal minval = BIG_REAL;
884
885 VecGetArray(nv, &nd);
886 VecGetArray(dv, &dd);
887 for (i = 0; i < N; i++) {
888 PetscReal nr = (PetscReal) nd[i];
889 PetscReal dr = (PetscReal) dd[i];
890 if (dr == ZERO)
891 continue;
892 else {
893 if (!notEvenOnce)
894 minval = SUNMIN(minval, nr/dr);
895 else {
896 minval = nr/dr;
897 notEvenOnce = SUNFALSE;
898 }
899 }
900 }
901 VecRestoreArray(nv, &nd);
902 VecRestoreArray(dv, &dd);
903
904 return(SUNMPI_Allreduce_scalar(minval, 3, comm));
905 }
906
907
908 /*
909 * -----------------------------------------------------------------
910 * fused vector operations
911 * -----------------------------------------------------------------
912 */
913
914
N_VLinearCombination_Petsc(int nvec,realtype * c,N_Vector * X,N_Vector z)915 int N_VLinearCombination_Petsc(int nvec, realtype* c, N_Vector* X, N_Vector z)
916 {
917 int i;
918 Vec* xv;
919 Vec zv;
920
921 /* invalid number of vectors */
922 if (nvec < 1) return(-1);
923
924 /* should have called N_VScale */
925 if (nvec == 1) {
926 N_VScale_Petsc(c[0], X[0], z);
927 return(0);
928 }
929
930 /* should have called N_VLinearSum */
931 if (nvec == 2) {
932 N_VLinearSum_Petsc(c[0], X[0], c[1], X[1], z);
933 return(0);
934 }
935
936 /* get petsc vectors */
937 xv = (Vec*) malloc(nvec * sizeof(Vec));
938 for (i=0; i<nvec; i++)
939 xv[i] = NV_PVEC_PTC(X[i]);
940
941 zv = NV_PVEC_PTC(z);
942
943 /*
944 * X[0] += c[i]*X[i], i = 1,...,nvec-1
945 */
946 if ((X[0] == z) && (c[0] == ONE)) {
947 VecMAXPY(zv, nvec-1, c+1, xv+1);
948 return(0);
949 }
950
951 /*
952 * X[0] = c[0] * X[0] + sum{ c[i] * X[i] }, i = 1,...,nvec-1
953 */
954 if (X[0] == z) {
955 VecScale(zv, c[0]);
956 VecMAXPY(zv, nvec-1, c+1, xv+1);
957 return(0);
958 }
959
960 /*
961 * z = sum{ c[i] * X[i] }, i = 0,...,nvec-1
962 */
963 VecAXPBY(zv, c[0], 0.0, xv[0]);
964 VecMAXPY(zv, nvec-1, c+1, xv+1);
965 return(0);
966 }
967
968
N_VScaleAddMulti_Petsc(int nvec,realtype * a,N_Vector x,N_Vector * Y,N_Vector * Z)969 int N_VScaleAddMulti_Petsc(int nvec, realtype* a, N_Vector x, N_Vector* Y, N_Vector* Z)
970 {
971 int i;
972 sunindextype j, N;
973 PetscScalar *xd, *yd, *zd;
974
975 /* invalid number of vectors */
976 if (nvec < 1) return(-1);
977
978 /* should have called N_VLinearSum */
979 if (nvec == 1) {
980 N_VLinearSum_Petsc(a[0], x, ONE, Y[0], Z[0]);
981 return(0);
982 }
983
984 /* get vector length and data array */
985 N = NV_LOCLENGTH_PTC(x);
986 VecGetArray(NV_PVEC_PTC(x), &xd);
987
988 /*
989 * Y[i][j] += a[i] * x[j]
990 */
991 if (Y == Z) {
992 for (i=0; i<nvec; i++) {
993 VecGetArray(NV_PVEC_PTC(Y[i]), &yd);
994 for (j=0; j<N; j++) {
995 yd[j] += a[i] * xd[j];
996 }
997 VecRestoreArray(NV_PVEC_PTC(Y[i]), &yd);
998 }
999 return(0);
1000 }
1001
1002 /*
1003 * Z[i][j] = Y[i][j] + a[i] * x[j]
1004 */
1005 for (i=0; i<nvec; i++) {
1006 VecGetArray(NV_PVEC_PTC(Y[i]), &yd);
1007 VecGetArray(NV_PVEC_PTC(Z[i]), &zd);
1008 for (j=0; j<N; j++) {
1009 zd[j] = a[i] * xd[j] + yd[j];
1010 }
1011 VecRestoreArray(NV_PVEC_PTC(Y[i]), &yd);
1012 VecRestoreArray(NV_PVEC_PTC(Z[i]), &zd);
1013 }
1014 return(0);
1015 }
1016
1017
N_VDotProdMulti_Petsc(int nvec,N_Vector x,N_Vector * Y,realtype * dotprods)1018 int N_VDotProdMulti_Petsc(int nvec, N_Vector x, N_Vector* Y, realtype* dotprods)
1019 {
1020 int i;
1021 Vec* yv;
1022 Vec xv;
1023
1024 /* invalid number of vectors */
1025 if (nvec < 1) return(-1);
1026
1027 /* should have called N_VDotProd */
1028 if (nvec == 1) {
1029 dotprods[0] = N_VDotProd_Petsc(x, Y[0]);
1030 return(0);
1031 }
1032
1033 /* get petsc vectors */
1034 yv = (Vec*) malloc(nvec * sizeof(Vec));
1035 for (i=0; i<nvec; i++)
1036 yv[i] = NV_PVEC_PTC(Y[i]);
1037
1038 xv = NV_PVEC_PTC(x);
1039
1040 VecMDot(xv, nvec, yv, dotprods);
1041 return(0);
1042 }
1043
1044
1045 /*
1046 * -----------------------------------------------------------------------------
1047 * vector array operations
1048 * -----------------------------------------------------------------------------
1049 */
1050
N_VLinearSumVectorArray_Petsc(int nvec,realtype a,N_Vector * X,realtype b,N_Vector * Y,N_Vector * Z)1051 int N_VLinearSumVectorArray_Petsc(int nvec,
1052 realtype a, N_Vector* X,
1053 realtype b, N_Vector* Y,
1054 N_Vector* Z)
1055 {
1056 int i;
1057 sunindextype j, N;
1058 PetscScalar *xd, *yd, *zd;
1059
1060 /* invalid number of vectors */
1061 if (nvec < 1) return(-1);
1062
1063 /* should have called N_VLinearSum */
1064 if (nvec == 1) {
1065 N_VLinearSum_Petsc(a, X[0], b, Y[0], Z[0]);
1066 return(0);
1067 }
1068
1069 /* get vector length */
1070 N = NV_LOCLENGTH_PTC(Z[0]);
1071
1072 /* compute linear sum for each vector pair in vector arrays */
1073 for (i=0; i<nvec; i++) {
1074 VecGetArray(NV_PVEC_PTC(X[i]), &xd);
1075 VecGetArray(NV_PVEC_PTC(Y[i]), &yd);
1076 VecGetArray(NV_PVEC_PTC(Z[i]), &zd);
1077 for (j=0; j<N; j++) {
1078 zd[j] = a * xd[j] + b * yd[j];
1079 }
1080 VecRestoreArray(NV_PVEC_PTC(X[i]), &xd);
1081 VecRestoreArray(NV_PVEC_PTC(Y[i]), &yd);
1082 VecRestoreArray(NV_PVEC_PTC(Z[i]), &zd);
1083 }
1084
1085 return(0);
1086 }
1087
1088
N_VScaleVectorArray_Petsc(int nvec,realtype * c,N_Vector * X,N_Vector * Z)1089 int N_VScaleVectorArray_Petsc(int nvec, realtype* c, N_Vector* X, N_Vector* Z)
1090 {
1091 int i;
1092 sunindextype j, N;
1093 PetscScalar *xd, *zd;
1094
1095 /* invalid number of vectors */
1096 if (nvec < 1) return(-1);
1097
1098 /* should have called N_VScale */
1099 if (nvec == 1) {
1100 N_VScale_Petsc(c[0], X[0], Z[0]);
1101 return(0);
1102 }
1103
1104 /* get vector length */
1105 N = NV_LOCLENGTH_PTC(Z[0]);
1106
1107 /*
1108 * X[i] *= c[i]
1109 */
1110 if (X == Z) {
1111 for (i=0; i<nvec; i++) {
1112 VecGetArray(NV_PVEC_PTC(X[i]), &xd);
1113 for (j=0; j<N; j++) {
1114 xd[j] *= c[i];
1115 }
1116 VecRestoreArray(NV_PVEC_PTC(X[i]), &xd);
1117 }
1118 return(0);
1119 }
1120
1121 /*
1122 * Z[i] = c[i] * X[i]
1123 */
1124 for (i=0; i<nvec; i++) {
1125 VecGetArray(NV_PVEC_PTC(X[i]), &xd);
1126 VecGetArray(NV_PVEC_PTC(Z[i]), &zd);
1127 for (j=0; j<N; j++) {
1128 zd[j] = c[i] * xd[j];
1129 }
1130 VecRestoreArray(NV_PVEC_PTC(X[i]), &xd);
1131 VecRestoreArray(NV_PVEC_PTC(Z[i]), &zd);
1132 }
1133 return(0);
1134 }
1135
1136
N_VConstVectorArray_Petsc(int nvec,realtype c,N_Vector * Z)1137 int N_VConstVectorArray_Petsc(int nvec, realtype c, N_Vector* Z)
1138 {
1139 int i;
1140 sunindextype j, N;
1141 PetscScalar *zd;
1142
1143 /* invalid number of vectors */
1144 if (nvec < 1) return(-1);
1145
1146 /* should have called N_VConst */
1147 if (nvec == 1) {
1148 N_VConst_Petsc(c, Z[0]);
1149 return(0);
1150 }
1151
1152 /* get vector length */
1153 N = NV_LOCLENGTH_PTC(Z[0]);
1154
1155 /* set each vector in the vector array to a constant */
1156 for (i=0; i<nvec; i++) {
1157 VecGetArray(NV_PVEC_PTC(Z[i]), &zd);
1158 for (j=0; j<N; j++) {
1159 zd[j] = c;
1160 }
1161 VecRestoreArray(NV_PVEC_PTC(Z[i]), &zd);
1162 }
1163
1164 return(0);
1165 }
1166
1167
N_VWrmsNormVectorArray_Petsc(int nvec,N_Vector * X,N_Vector * W,realtype * nrm)1168 int N_VWrmsNormVectorArray_Petsc(int nvec, N_Vector* X, N_Vector* W, realtype* nrm)
1169 {
1170 int i;
1171 sunindextype j, Nl, Ng;
1172 realtype* wd=NULL;
1173 realtype* xd=NULL;
1174 MPI_Comm comm;
1175
1176 /* invalid number of vectors */
1177 if (nvec < 1) return(-1);
1178
1179 /* should have called N_VWrmsNorm */
1180 if (nvec == 1) {
1181 nrm[0] = N_VWrmsNorm_Petsc(X[0], W[0]);
1182 return(0);
1183 }
1184
1185 /* get vector lengths and communicator */
1186 Nl = NV_LOCLENGTH_PTC(X[0]);
1187 Ng = NV_GLOBLENGTH_PTC(X[0]);
1188 comm = NV_COMM_PTC(X[0]);
1189
1190 /* compute the WRMS norm for each vector in the vector array */
1191 for (i=0; i<nvec; i++) {
1192 VecGetArray(NV_PVEC_PTC(X[i]), &xd);
1193 VecGetArray(NV_PVEC_PTC(W[i]), &wd);
1194 nrm[i] = ZERO;
1195 for (j=0; j<Nl; j++) {
1196 nrm[i] += PetscSqr(PetscAbsScalar(xd[j] * wd[j]));
1197 }
1198 VecRestoreArray(NV_PVEC_PTC(X[i]), &xd);
1199 VecRestoreArray(NV_PVEC_PTC(W[i]), &wd);
1200 }
1201 SUNMPI_Allreduce(nrm, nvec, 1, comm);
1202
1203 for (i=0; i<nvec; i++)
1204 nrm[i] = SUNRsqrt(nrm[i]/Ng);
1205
1206 return(0);
1207 }
1208
1209
N_VWrmsNormMaskVectorArray_Petsc(int nvec,N_Vector * X,N_Vector * W,N_Vector id,realtype * nrm)1210 int N_VWrmsNormMaskVectorArray_Petsc(int nvec, N_Vector* X, N_Vector* W,
1211 N_Vector id, realtype* nrm)
1212 {
1213 int i;
1214 sunindextype j, Nl, Ng;
1215 PetscScalar *wd, *xd, *idd;
1216 MPI_Comm comm;
1217
1218 /* invalid number of vectors */
1219 if (nvec < 1) return(-1);
1220
1221 /* should have called N_VWrmsNorm */
1222 if (nvec == 1) {
1223 nrm[0] = N_VWrmsNormMask_Petsc(X[0], W[0], id);
1224 return(0);
1225 }
1226
1227 /* get vector lengths and communicator */
1228 Nl = NV_LOCLENGTH_PTC(X[0]);
1229 Ng = NV_GLOBLENGTH_PTC(X[0]);
1230 comm = NV_COMM_PTC(X[0]);
1231
1232 /* compute the WRMS norm for each vector in the vector array */
1233 VecGetArray(NV_PVEC_PTC(id), &idd);
1234 for (i=0; i<nvec; i++) {
1235 VecGetArray(NV_PVEC_PTC(X[i]), &xd);
1236 VecGetArray(NV_PVEC_PTC(W[i]), &wd);
1237 nrm[i] = ZERO;
1238 for (j=0; j<Nl; j++) {
1239 if (idd[j] > ZERO)
1240 nrm[i] += SUNSQR(xd[j] * wd[j]);
1241 }
1242 VecRestoreArray(NV_PVEC_PTC(X[i]), &xd);
1243 VecRestoreArray(NV_PVEC_PTC(W[i]), &wd);
1244 }
1245 VecRestoreArray(NV_PVEC_PTC(id), &idd);
1246
1247 SUNMPI_Allreduce(nrm, nvec, 1, comm);
1248
1249 for (i=0; i<nvec; i++)
1250 nrm[i] = SUNRsqrt(nrm[i]/Ng);
1251
1252 return(0);
1253 }
1254
1255
N_VScaleAddMultiVectorArray_Petsc(int nvec,int nsum,realtype * a,N_Vector * X,N_Vector ** Y,N_Vector ** Z)1256 int N_VScaleAddMultiVectorArray_Petsc(int nvec, int nsum, realtype* a,
1257 N_Vector* X, N_Vector** Y, N_Vector** Z)
1258 {
1259 int i, j;
1260 sunindextype k, N;
1261 PetscScalar *xd, *yd, *zd;
1262
1263 int retval;
1264 N_Vector* YY;
1265 N_Vector* ZZ;
1266
1267 /* invalid number of vectors */
1268 if (nvec < 1) return(-1);
1269 if (nsum < 1) return(-1);
1270
1271 /* ---------------------------
1272 * Special cases for nvec == 1
1273 * --------------------------- */
1274
1275 if (nvec == 1) {
1276
1277 /* should have called N_VLinearSum */
1278 if (nsum == 1) {
1279 N_VLinearSum_Petsc(a[0], X[0], ONE, Y[0][0], Z[0][0]);
1280 return(0);
1281 }
1282
1283 /* should have called N_VScaleAddMulti */
1284 YY = (N_Vector *) malloc(nsum * sizeof(N_Vector));
1285 ZZ = (N_Vector *) malloc(nsum * sizeof(N_Vector));
1286
1287 for (j=0; j<nsum; j++) {
1288 YY[j] = Y[j][0];
1289 ZZ[j] = Z[j][0];
1290 }
1291
1292 retval = N_VScaleAddMulti_Petsc(nsum, a, X[0], YY, ZZ);
1293
1294 free(YY);
1295 free(ZZ);
1296 return(retval);
1297 }
1298
1299 /* --------------------------
1300 * Special cases for nvec > 1
1301 * -------------------------- */
1302
1303 /* should have called N_VLinearSumVectorArray */
1304 if (nsum == 1) {
1305 retval = N_VLinearSumVectorArray_Petsc(nvec, a[0], X, ONE, Y[0], Z[0]);
1306 return(retval);
1307 }
1308
1309 /* ----------------------------
1310 * Compute multiple linear sums
1311 * ---------------------------- */
1312
1313 /* get vector length */
1314 N = NV_LOCLENGTH_PTC(X[0]);
1315
1316 /*
1317 * Y[i][j] += a[i] * x[j]
1318 */
1319 if (Y == Z) {
1320 for (i=0; i<nvec; i++) {
1321 VecGetArray(NV_PVEC_PTC(X[i]), &xd);
1322 for (j=0; j<nsum; j++) {
1323 VecGetArray(NV_PVEC_PTC(Y[j][i]), &yd);
1324 for (k=0; k<N; k++) {
1325 yd[k] += a[j] * xd[k];
1326 }
1327 VecRestoreArray(NV_PVEC_PTC(Y[j][i]), &yd);
1328 }
1329 VecRestoreArray(NV_PVEC_PTC(X[i]), &xd);
1330 }
1331 return(0);
1332 }
1333
1334 /*
1335 * Z[i][j] = Y[i][j] + a[i] * x[j]
1336 */
1337 for (i=0; i<nvec; i++) {
1338 VecGetArray(NV_PVEC_PTC(X[i]), &xd);
1339 for (j=0; j<nsum; j++) {
1340 VecGetArray(NV_PVEC_PTC(Y[j][i]), &yd);
1341 VecGetArray(NV_PVEC_PTC(Z[j][i]), &zd);
1342 for (k=0; k<N; k++) {
1343 zd[k] = a[j] * xd[k] + yd[k];
1344 }
1345 VecRestoreArray(NV_PVEC_PTC(Y[j][i]), &yd);
1346 VecRestoreArray(NV_PVEC_PTC(Z[j][i]), &zd);
1347 }
1348 VecRestoreArray(NV_PVEC_PTC(X[i]), &xd);
1349 }
1350
1351 return(0);
1352 }
1353
1354
N_VLinearCombinationVectorArray_Petsc(int nvec,int nsum,realtype * c,N_Vector ** X,N_Vector * Z)1355 int N_VLinearCombinationVectorArray_Petsc(int nvec, int nsum,
1356 realtype* c,
1357 N_Vector** X,
1358 N_Vector* Z)
1359 {
1360 int i; /* vector arrays index in summation [0,nsum) */
1361 int j; /* vector index in vector array [0,nvec) */
1362 sunindextype k; /* element index in vector [0,N) */
1363 sunindextype N;
1364 PetscScalar *zd, *xd;
1365
1366 realtype* ctmp;
1367 N_Vector* Y;
1368
1369 /* invalid number of vectors */
1370 if (nvec < 1) return(-1);
1371 if (nsum < 1) return(-1);
1372
1373 /* ---------------------------
1374 * Special cases for nvec == 1
1375 * --------------------------- */
1376
1377 if (nvec == 1) {
1378
1379 /* should have called N_VScale */
1380 if (nsum == 1) {
1381 N_VScale_Petsc(c[0], X[0][0], Z[0]);
1382 return(0);
1383 }
1384
1385 /* should have called N_VLinearSum */
1386 if (nsum == 2) {
1387 N_VLinearSum_Petsc(c[0], X[0][0], c[1], X[1][0], Z[0]);
1388 return(0);
1389 }
1390
1391 /* should have called N_VLinearCombination */
1392 Y = (N_Vector *) malloc(nsum * sizeof(N_Vector));
1393
1394 for (i=0; i<nsum; i++) {
1395 Y[i] = X[i][0];
1396 }
1397
1398 N_VLinearCombination_Petsc(nsum, c, Y, Z[0]);
1399
1400 free(Y);
1401 return(0);
1402 }
1403
1404 /* --------------------------
1405 * Special cases for nvec > 1
1406 * -------------------------- */
1407
1408 /* should have called N_VScaleVectorArray */
1409 if (nsum == 1) {
1410
1411 ctmp = (realtype*) malloc(nvec * sizeof(realtype));
1412
1413 for (j=0; j<nvec; j++) {
1414 ctmp[j] = c[0];
1415 }
1416
1417 N_VScaleVectorArray_Petsc(nvec, ctmp, X[0], Z);
1418
1419 free(ctmp);
1420 return(0);
1421 }
1422
1423 /* should have called N_VLinearSumVectorArray */
1424 if (nsum == 2) {
1425 N_VLinearSumVectorArray_Petsc(nvec, c[0], X[0], c[1], X[1], Z);
1426 return(0);
1427 }
1428
1429 /* --------------------------
1430 * Compute linear combination
1431 * -------------------------- */
1432
1433 /* get vector length */
1434 N = NV_LOCLENGTH_PTC(Z[0]);
1435
1436 /*
1437 * X[0][j] += c[i]*X[i][j], i = 1,...,nvec-1
1438 */
1439 if ((X[0] == Z) && (c[0] == ONE)) {
1440 for (j=0; j<nvec; j++) {
1441 VecGetArray(NV_PVEC_PTC(Z[j]), &zd);
1442 for (i=1; i<nsum; i++) {
1443 VecGetArray(NV_PVEC_PTC(X[i][j]), &xd);
1444 for (k=0; k<N; k++) {
1445 zd[k] += c[i] * xd[k];
1446 }
1447 VecRestoreArray(NV_PVEC_PTC(X[i][j]), &xd);
1448 }
1449 VecRestoreArray(NV_PVEC_PTC(Z[j]), &zd);
1450 }
1451 return(0);
1452 }
1453
1454 /*
1455 * X[0][j] = c[0] * X[0][j] + sum{ c[i] * X[i][j] }, i = 1,...,nvec-1
1456 */
1457 if (X[0] == Z) {
1458 for (j=0; j<nvec; j++) {
1459 VecGetArray(NV_PVEC_PTC(Z[j]), &zd);
1460 for (k=0; k<N; k++) {
1461 zd[k] *= c[0];
1462 }
1463 for (i=1; i<nsum; i++) {
1464 VecGetArray(NV_PVEC_PTC(X[i][j]), &xd);
1465 for (k=0; k<N; k++) {
1466 zd[k] += c[i] * xd[k];
1467 }
1468 VecRestoreArray(NV_PVEC_PTC(X[i][j]), &xd);
1469 }
1470 VecRestoreArray(NV_PVEC_PTC(Z[j]), &zd);
1471 }
1472 return(0);
1473 }
1474
1475 /*
1476 * Z[j] = sum{ c[i] * X[i][j] }, i = 0,...,nvec-1
1477 */
1478 for (j=0; j<nvec; j++) {
1479 VecGetArray(NV_PVEC_PTC(X[0][j]), &xd);
1480 VecGetArray(NV_PVEC_PTC(Z[j]), &zd);
1481 for (k=0; k<N; k++) {
1482 zd[k] = c[0] * xd[k];
1483 }
1484 VecRestoreArray(NV_PVEC_PTC(X[0][j]), &xd);
1485 for (i=1; i<nsum; i++) {
1486 VecGetArray(NV_PVEC_PTC(X[i][j]), &xd);
1487 for (k=0; k<N; k++) {
1488 zd[k] += c[i] * xd[k];
1489 }
1490 VecRestoreArray(NV_PVEC_PTC(X[i][j]), &xd);
1491 }
1492 VecRestoreArray(NV_PVEC_PTC(Z[j]), &zd);
1493 }
1494 return(0);
1495 }
1496
1497 /*
1498 * -----------------------------------------------------------------
1499 * Enable / Disable fused and vector array operations
1500 * -----------------------------------------------------------------
1501 */
1502
N_VEnableFusedOps_Petsc(N_Vector v,booleantype tf)1503 int N_VEnableFusedOps_Petsc(N_Vector v, booleantype tf)
1504 {
1505 /* check that vector is non-NULL */
1506 if (v == NULL) return(-1);
1507
1508 /* check that ops structure is non-NULL */
1509 if (v->ops == NULL) return(-1);
1510
1511 if (tf) {
1512 /* enable all fused vector operations */
1513 v->ops->nvlinearcombination = N_VLinearCombination_Petsc;
1514 v->ops->nvscaleaddmulti = N_VScaleAddMulti_Petsc;
1515 v->ops->nvdotprodmulti = N_VDotProdMulti_Petsc;
1516 /* enable all vector array operations */
1517 v->ops->nvlinearsumvectorarray = N_VLinearSumVectorArray_Petsc;
1518 v->ops->nvscalevectorarray = N_VScaleVectorArray_Petsc;
1519 v->ops->nvconstvectorarray = N_VConstVectorArray_Petsc;
1520 v->ops->nvwrmsnormvectorarray = N_VWrmsNormVectorArray_Petsc;
1521 v->ops->nvwrmsnormmaskvectorarray = N_VWrmsNormMaskVectorArray_Petsc;
1522 v->ops->nvscaleaddmultivectorarray = N_VScaleAddMultiVectorArray_Petsc;
1523 v->ops->nvlinearcombinationvectorarray = N_VLinearCombinationVectorArray_Petsc;
1524 } else {
1525 /* disable all fused vector operations */
1526 v->ops->nvlinearcombination = NULL;
1527 v->ops->nvscaleaddmulti = NULL;
1528 v->ops->nvdotprodmulti = NULL;
1529 /* disable all vector array operations */
1530 v->ops->nvlinearsumvectorarray = NULL;
1531 v->ops->nvscalevectorarray = NULL;
1532 v->ops->nvconstvectorarray = NULL;
1533 v->ops->nvwrmsnormvectorarray = NULL;
1534 v->ops->nvwrmsnormmaskvectorarray = NULL;
1535 v->ops->nvscaleaddmultivectorarray = NULL;
1536 v->ops->nvlinearcombinationvectorarray = NULL;
1537 }
1538
1539 /* return success */
1540 return(0);
1541 }
1542
1543
N_VEnableLinearCombination_Petsc(N_Vector v,booleantype tf)1544 int N_VEnableLinearCombination_Petsc(N_Vector v, booleantype tf)
1545 {
1546 /* check that vector is non-NULL */
1547 if (v == NULL) return(-1);
1548
1549 /* check that ops structure is non-NULL */
1550 if (v->ops == NULL) return(-1);
1551
1552 /* enable/disable operation */
1553 if (tf)
1554 v->ops->nvlinearcombination = N_VLinearCombination_Petsc;
1555 else
1556 v->ops->nvlinearcombination = NULL;
1557
1558 /* return success */
1559 return(0);
1560 }
1561
N_VEnableScaleAddMulti_Petsc(N_Vector v,booleantype tf)1562 int N_VEnableScaleAddMulti_Petsc(N_Vector v, booleantype tf)
1563 {
1564 /* check that vector is non-NULL */
1565 if (v == NULL) return(-1);
1566
1567 /* check that ops structure is non-NULL */
1568 if (v->ops == NULL) return(-1);
1569
1570 /* enable/disable operation */
1571 if (tf)
1572 v->ops->nvscaleaddmulti = N_VScaleAddMulti_Petsc;
1573 else
1574 v->ops->nvscaleaddmulti = NULL;
1575
1576 /* return success */
1577 return(0);
1578 }
1579
N_VEnableDotProdMulti_Petsc(N_Vector v,booleantype tf)1580 int N_VEnableDotProdMulti_Petsc(N_Vector v, booleantype tf)
1581 {
1582 /* check that vector is non-NULL */
1583 if (v == NULL) return(-1);
1584
1585 /* check that ops structure is non-NULL */
1586 if (v->ops == NULL) return(-1);
1587
1588 /* enable/disable operation */
1589 if (tf)
1590 v->ops->nvdotprodmulti = N_VDotProdMulti_Petsc;
1591 else
1592 v->ops->nvdotprodmulti = NULL;
1593
1594 /* return success */
1595 return(0);
1596 }
1597
N_VEnableLinearSumVectorArray_Petsc(N_Vector v,booleantype tf)1598 int N_VEnableLinearSumVectorArray_Petsc(N_Vector v, booleantype tf)
1599 {
1600 /* check that vector is non-NULL */
1601 if (v == NULL) return(-1);
1602
1603 /* check that ops structure is non-NULL */
1604 if (v->ops == NULL) return(-1);
1605
1606 /* enable/disable operation */
1607 if (tf)
1608 v->ops->nvlinearsumvectorarray = N_VLinearSumVectorArray_Petsc;
1609 else
1610 v->ops->nvlinearsumvectorarray = NULL;
1611
1612 /* return success */
1613 return(0);
1614 }
1615
N_VEnableScaleVectorArray_Petsc(N_Vector v,booleantype tf)1616 int N_VEnableScaleVectorArray_Petsc(N_Vector v, booleantype tf)
1617 {
1618 /* check that vector is non-NULL */
1619 if (v == NULL) return(-1);
1620
1621 /* check that ops structure is non-NULL */
1622 if (v->ops == NULL) return(-1);
1623
1624 /* enable/disable operation */
1625 if (tf)
1626 v->ops->nvscalevectorarray = N_VScaleVectorArray_Petsc;
1627 else
1628 v->ops->nvscalevectorarray = NULL;
1629
1630 /* return success */
1631 return(0);
1632 }
1633
N_VEnableConstVectorArray_Petsc(N_Vector v,booleantype tf)1634 int N_VEnableConstVectorArray_Petsc(N_Vector v, booleantype tf)
1635 {
1636 /* check that vector is non-NULL */
1637 if (v == NULL) return(-1);
1638
1639 /* check that ops structure is non-NULL */
1640 if (v->ops == NULL) return(-1);
1641
1642 /* enable/disable operation */
1643 if (tf)
1644 v->ops->nvconstvectorarray = N_VConstVectorArray_Petsc;
1645 else
1646 v->ops->nvconstvectorarray = NULL;
1647
1648 /* return success */
1649 return(0);
1650 }
1651
N_VEnableWrmsNormVectorArray_Petsc(N_Vector v,booleantype tf)1652 int N_VEnableWrmsNormVectorArray_Petsc(N_Vector v, booleantype tf)
1653 {
1654 /* check that vector is non-NULL */
1655 if (v == NULL) return(-1);
1656
1657 /* check that ops structure is non-NULL */
1658 if (v->ops == NULL) return(-1);
1659
1660 /* enable/disable operation */
1661 if (tf)
1662 v->ops->nvwrmsnormvectorarray = N_VWrmsNormVectorArray_Petsc;
1663 else
1664 v->ops->nvwrmsnormvectorarray = NULL;
1665
1666 /* return success */
1667 return(0);
1668 }
1669
N_VEnableWrmsNormMaskVectorArray_Petsc(N_Vector v,booleantype tf)1670 int N_VEnableWrmsNormMaskVectorArray_Petsc(N_Vector v, booleantype tf)
1671 {
1672 /* check that vector is non-NULL */
1673 if (v == NULL) return(-1);
1674
1675 /* check that ops structure is non-NULL */
1676 if (v->ops == NULL) return(-1);
1677
1678 /* enable/disable operation */
1679 if (tf)
1680 v->ops->nvwrmsnormmaskvectorarray = N_VWrmsNormMaskVectorArray_Petsc;
1681 else
1682 v->ops->nvwrmsnormmaskvectorarray = NULL;
1683
1684 /* return success */
1685 return(0);
1686 }
1687
N_VEnableScaleAddMultiVectorArray_Petsc(N_Vector v,booleantype tf)1688 int N_VEnableScaleAddMultiVectorArray_Petsc(N_Vector v, booleantype tf)
1689 {
1690 /* check that vector is non-NULL */
1691 if (v == NULL) return(-1);
1692
1693 /* check that ops structure is non-NULL */
1694 if (v->ops == NULL) return(-1);
1695
1696 /* enable/disable operation */
1697 if (tf)
1698 v->ops->nvscaleaddmultivectorarray = N_VScaleAddMultiVectorArray_Petsc;
1699 else
1700 v->ops->nvscaleaddmultivectorarray = NULL;
1701
1702 /* return success */
1703 return(0);
1704 }
1705
N_VEnableLinearCombinationVectorArray_Petsc(N_Vector v,booleantype tf)1706 int N_VEnableLinearCombinationVectorArray_Petsc(N_Vector v, booleantype tf)
1707 {
1708 /* check that vector is non-NULL */
1709 if (v == NULL) return(-1);
1710
1711 /* check that ops structure is non-NULL */
1712 if (v->ops == NULL) return(-1);
1713
1714 /* enable/disable operation */
1715 if (tf)
1716 v->ops->nvlinearcombinationvectorarray = N_VLinearCombinationVectorArray_Petsc;
1717 else
1718 v->ops->nvlinearcombinationvectorarray = NULL;
1719
1720 /* return success */
1721 return(0);
1722 }
1723