1
2 /*
3 Some useful vector utility functions.
4 */
5 #include <../src/vec/vec/impls/mpi/pvecimpl.h> /*I "petscvec.h" I*/
6 extern MPI_Op MPIU_MAXINDEX_OP;
7 extern MPI_Op MPIU_MININDEX_OP;
8
9 /*@
10 VecStrideSet - Sets a subvector of a vector defined
11 by a starting point and a stride with a given value
12
13 Logically Collective on Vec
14
15 Input Parameter:
16 + v - the vector
17 . start - starting point of the subvector (defined by a stride)
18 - s - value to set for each entry in that subvector
19
20 Notes:
21 One must call VecSetBlockSize() before this routine to set the stride
22 information, or use a vector created from a multicomponent DMDA.
23
24 This will only work if the desire subvector is a stride subvector
25
26 Level: advanced
27
28
29 .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
30 @*/
VecStrideSet(Vec v,PetscInt start,PetscScalar s)31 PetscErrorCode VecStrideSet(Vec v,PetscInt start,PetscScalar s)
32 {
33 PetscErrorCode ierr;
34 PetscInt i,n,bs;
35 PetscScalar *x;
36
37 PetscFunctionBegin;
38 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
39 PetscValidLogicalCollectiveInt(v,start,2);
40 PetscValidLogicalCollectiveScalar(v,s,3);
41 ierr = VecSetErrorIfLocked(v,1);CHKERRQ(ierr);
42
43 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
44 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
45 ierr = VecGetBlockSize(v,&bs);CHKERRQ(ierr);
46 if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
47 else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
48 x += start;
49 for (i=0; i<n; i+=bs) x[i] = s;
50 x -= start;
51 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
52 PetscFunctionReturn(0);
53 }
54
55 /*@
56 VecStrideScale - Scales a subvector of a vector defined
57 by a starting point and a stride.
58
59 Logically Collective on Vec
60
61 Input Parameter:
62 + v - the vector
63 . start - starting point of the subvector (defined by a stride)
64 - scale - value to multiply each subvector entry by
65
66 Notes:
67 One must call VecSetBlockSize() before this routine to set the stride
68 information, or use a vector created from a multicomponent DMDA.
69
70 This will only work if the desire subvector is a stride subvector
71
72 Level: advanced
73
74
75 .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
76 @*/
VecStrideScale(Vec v,PetscInt start,PetscScalar scale)77 PetscErrorCode VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
78 {
79 PetscErrorCode ierr;
80 PetscInt i,n,bs;
81 PetscScalar *x;
82
83 PetscFunctionBegin;
84 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
85 PetscValidLogicalCollectiveInt(v,start,2);
86 PetscValidLogicalCollectiveScalar(v,scale,3);
87 ierr = VecSetErrorIfLocked(v,1);CHKERRQ(ierr);
88
89 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
90 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
91 ierr = VecGetBlockSize(v,&bs);CHKERRQ(ierr);
92 if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
93 else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
94 x += start;
95 for (i=0; i<n; i+=bs) x[i] *= scale;
96 x -= start;
97 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
98 PetscFunctionReturn(0);
99 }
100
101 /*@
102 VecStrideNorm - Computes the norm of subvector of a vector defined
103 by a starting point and a stride.
104
105 Collective on Vec
106
107 Input Parameter:
108 + v - the vector
109 . start - starting point of the subvector (defined by a stride)
110 - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
111
112 Output Parameter:
113 . norm - the norm
114
115 Notes:
116 One must call VecSetBlockSize() before this routine to set the stride
117 information, or use a vector created from a multicomponent DMDA.
118
119 If x is the array representing the vector x then this computes the norm
120 of the array (x[start],x[start+stride],x[start+2*stride], ....)
121
122 This is useful for computing, say the norm of the pressure variable when
123 the pressure is stored (interlaced) with other variables, say density etc.
124
125 This will only work if the desire subvector is a stride subvector
126
127 Level: advanced
128
129
130 .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
131 @*/
VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal * nrm)132 PetscErrorCode VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
133 {
134 PetscErrorCode ierr;
135 PetscInt i,n,bs;
136 const PetscScalar *x;
137 PetscReal tnorm;
138 MPI_Comm comm;
139
140 PetscFunctionBegin;
141 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
142 PetscValidRealPointer(nrm,3);
143 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
144 ierr = VecGetArrayRead(v,&x);CHKERRQ(ierr);
145 ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
146
147 ierr = VecGetBlockSize(v,&bs);CHKERRQ(ierr);
148 if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
149 else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
150 x += start;
151
152 if (ntype == NORM_2) {
153 PetscScalar sum = 0.0;
154 for (i=0; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
155 tnorm = PetscRealPart(sum);
156 ierr = MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);CHKERRQ(ierr);
157 *nrm = PetscSqrtReal(*nrm);
158 } else if (ntype == NORM_1) {
159 tnorm = 0.0;
160 for (i=0; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
161 ierr = MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);CHKERRQ(ierr);
162 } else if (ntype == NORM_INFINITY) {
163 PetscReal tmp;
164 tnorm = 0.0;
165
166 for (i=0; i<n; i+=bs) {
167 if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
168 /* check special case of tmp == NaN */
169 if (tmp != tmp) {tnorm = tmp; break;}
170 }
171 ierr = MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,comm);CHKERRQ(ierr);
172 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
173 ierr = VecRestoreArrayRead(v,&x);CHKERRQ(ierr);
174 PetscFunctionReturn(0);
175 }
176
177 /*@
178 VecStrideMax - Computes the maximum of subvector of a vector defined
179 by a starting point and a stride and optionally its location.
180
181 Collective on Vec
182
183 Input Parameter:
184 + v - the vector
185 - start - starting point of the subvector (defined by a stride)
186
187 Output Parameter:
188 + index - the location where the maximum occurred (pass NULL if not required)
189 - nrm - the maximum value in the subvector
190
191 Notes:
192 One must call VecSetBlockSize() before this routine to set the stride
193 information, or use a vector created from a multicomponent DMDA.
194
195 If xa is the array representing the vector x, then this computes the max
196 of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
197
198 This is useful for computing, say the maximum of the pressure variable when
199 the pressure is stored (interlaced) with other variables, e.g., density, etc.
200 This will only work if the desire subvector is a stride subvector.
201
202 Level: advanced
203
204
205 .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
206 @*/
VecStrideMax(Vec v,PetscInt start,PetscInt * idex,PetscReal * nrm)207 PetscErrorCode VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
208 {
209 PetscErrorCode ierr;
210 PetscInt i,n,bs,id;
211 const PetscScalar *x;
212 PetscReal max,tmp;
213 MPI_Comm comm;
214
215 PetscFunctionBegin;
216 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
217 PetscValidRealPointer(nrm,3);
218
219 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
220 ierr = VecGetArrayRead(v,&x);CHKERRQ(ierr);
221 ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
222
223 ierr = VecGetBlockSize(v,&bs);CHKERRQ(ierr);
224 if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
225 else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
226 x += start;
227
228 id = -1;
229 if (!n) max = PETSC_MIN_REAL;
230 else {
231 id = 0;
232 max = PetscRealPart(x[0]);
233 for (i=bs; i<n; i+=bs) {
234 if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
235 }
236 }
237 ierr = VecRestoreArrayRead(v,&x);CHKERRQ(ierr);
238
239 if (!idex) {
240 ierr = MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,comm);CHKERRQ(ierr);
241 } else {
242 PetscReal in[2],out[2];
243 PetscInt rstart;
244
245 ierr = VecGetOwnershipRange(v,&rstart,NULL);CHKERRQ(ierr);
246 in[0] = max;
247 in[1] = rstart+id+start;
248 ierr = MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MAXINDEX_OP,PetscObjectComm((PetscObject)v));CHKERRQ(ierr);
249 *nrm = out[0];
250 *idex = (PetscInt)out[1];
251 }
252 PetscFunctionReturn(0);
253 }
254
255 /*@
256 VecStrideMin - Computes the minimum of subvector of a vector defined
257 by a starting point and a stride and optionally its location.
258
259 Collective on Vec
260
261 Input Parameter:
262 + v - the vector
263 - start - starting point of the subvector (defined by a stride)
264
265 Output Parameter:
266 + idex - the location where the minimum occurred. (pass NULL if not required)
267 - nrm - the minimum value in the subvector
268
269 Level: advanced
270
271 Notes:
272 One must call VecSetBlockSize() before this routine to set the stride
273 information, or use a vector created from a multicomponent DMDA.
274
275 If xa is the array representing the vector x, then this computes the min
276 of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
277
278 This is useful for computing, say the minimum of the pressure variable when
279 the pressure is stored (interlaced) with other variables, e.g., density, etc.
280 This will only work if the desire subvector is a stride subvector.
281
282
283 .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
284 @*/
VecStrideMin(Vec v,PetscInt start,PetscInt * idex,PetscReal * nrm)285 PetscErrorCode VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
286 {
287 PetscErrorCode ierr;
288 PetscInt i,n,bs,id;
289 const PetscScalar *x;
290 PetscReal min,tmp;
291 MPI_Comm comm;
292
293 PetscFunctionBegin;
294 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
295 PetscValidRealPointer(nrm,4);
296
297 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
298 ierr = VecGetArrayRead(v,&x);CHKERRQ(ierr);
299 ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
300
301 ierr = VecGetBlockSize(v,&bs);CHKERRQ(ierr);
302 if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
303 else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\nHave you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
304 x += start;
305
306 id = -1;
307 if (!n) min = PETSC_MAX_REAL;
308 else {
309 id = 0;
310 min = PetscRealPart(x[0]);
311 for (i=bs; i<n; i+=bs) {
312 if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
313 }
314 }
315 ierr = VecRestoreArrayRead(v,&x);CHKERRQ(ierr);
316
317 if (!idex) {
318 ierr = MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);CHKERRQ(ierr);
319 } else {
320 PetscReal in[2],out[2];
321 PetscInt rstart;
322
323 ierr = VecGetOwnershipRange(v,&rstart,NULL);CHKERRQ(ierr);
324 in[0] = min;
325 in[1] = rstart+id;
326 ierr = MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MININDEX_OP,PetscObjectComm((PetscObject)v));CHKERRQ(ierr);
327 *nrm = out[0];
328 *idex = (PetscInt)out[1];
329 }
330 PetscFunctionReturn(0);
331 }
332
333 /*@
334 VecStrideScaleAll - Scales the subvectors of a vector defined
335 by a starting point and a stride.
336
337 Logically Collective on Vec
338
339 Input Parameter:
340 + v - the vector
341 - scales - values to multiply each subvector entry by
342
343 Notes:
344 One must call VecSetBlockSize() before this routine to set the stride
345 information, or use a vector created from a multicomponent DMDA.
346
347 The dimension of scales must be the same as the vector block size
348
349
350 Level: advanced
351
352
353 .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
354 @*/
VecStrideScaleAll(Vec v,const PetscScalar * scales)355 PetscErrorCode VecStrideScaleAll(Vec v,const PetscScalar *scales)
356 {
357 PetscErrorCode ierr;
358 PetscInt i,j,n,bs;
359 PetscScalar *x;
360
361 PetscFunctionBegin;
362 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
363 PetscValidScalarPointer(scales,2);
364 ierr = VecSetErrorIfLocked(v,1);CHKERRQ(ierr);
365
366 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
367 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
368 ierr = VecGetBlockSize(v,&bs);CHKERRQ(ierr);
369
370 /* need to provide optimized code for each bs */
371 for (i=0; i<n; i+=bs) {
372 for (j=0; j<bs; j++) x[i+j] *= scales[j];
373 }
374 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
375 PetscFunctionReturn(0);
376 }
377
378 /*@
379 VecStrideNormAll - Computes the norms of subvectors of a vector defined
380 by a starting point and a stride.
381
382 Collective on Vec
383
384 Input Parameter:
385 + v - the vector
386 - ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
387
388 Output Parameter:
389 . nrm - the norms
390
391 Notes:
392 One must call VecSetBlockSize() before this routine to set the stride
393 information, or use a vector created from a multicomponent DMDA.
394
395 If x is the array representing the vector x then this computes the norm
396 of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride
397
398 The dimension of nrm must be the same as the vector block size
399
400 This will only work if the desire subvector is a stride subvector
401
402 Level: advanced
403
404
405 .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
406 @*/
VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])407 PetscErrorCode VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
408 {
409 PetscErrorCode ierr;
410 PetscInt i,j,n,bs;
411 const PetscScalar *x;
412 PetscReal tnorm[128];
413 MPI_Comm comm;
414
415 PetscFunctionBegin;
416 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
417 PetscValidRealPointer(nrm,3);
418 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
419 ierr = VecGetArrayRead(v,&x);CHKERRQ(ierr);
420 ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
421
422 ierr = VecGetBlockSize(v,&bs);CHKERRQ(ierr);
423 if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
424
425 if (ntype == NORM_2) {
426 PetscScalar sum[128];
427 for (j=0; j<bs; j++) sum[j] = 0.0;
428 for (i=0; i<n; i+=bs) {
429 for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
430 }
431 for (j=0; j<bs; j++) tnorm[j] = PetscRealPart(sum[j]);
432
433 ierr = MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);CHKERRQ(ierr);
434 for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
435 } else if (ntype == NORM_1) {
436 for (j=0; j<bs; j++) tnorm[j] = 0.0;
437
438 for (i=0; i<n; i+=bs) {
439 for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
440 }
441
442 ierr = MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);CHKERRQ(ierr);
443 } else if (ntype == NORM_INFINITY) {
444 PetscReal tmp;
445 for (j=0; j<bs; j++) tnorm[j] = 0.0;
446
447 for (i=0; i<n; i+=bs) {
448 for (j=0; j<bs; j++) {
449 if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
450 /* check special case of tmp == NaN */
451 if (tmp != tmp) {tnorm[j] = tmp; break;}
452 }
453 }
454 ierr = MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);CHKERRQ(ierr);
455 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
456 ierr = VecRestoreArrayRead(v,&x);CHKERRQ(ierr);
457 PetscFunctionReturn(0);
458 }
459
460 /*@
461 VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
462 by a starting point and a stride and optionally its location.
463
464 Collective on Vec
465
466 Input Parameter:
467 . v - the vector
468
469 Output Parameter:
470 + index - the location where the maximum occurred (not supported, pass NULL,
471 if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
472 - nrm - the maximum values of each subvector
473
474 Notes:
475 One must call VecSetBlockSize() before this routine to set the stride
476 information, or use a vector created from a multicomponent DMDA.
477
478 The dimension of nrm must be the same as the vector block size
479
480 Level: advanced
481
482
483 .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
484 @*/
VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])485 PetscErrorCode VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
486 {
487 PetscErrorCode ierr;
488 PetscInt i,j,n,bs;
489 const PetscScalar *x;
490 PetscReal max[128],tmp;
491 MPI_Comm comm;
492
493 PetscFunctionBegin;
494 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
495 PetscValidRealPointer(nrm,3);
496 if (idex) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
497 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
498 ierr = VecGetArrayRead(v,&x);CHKERRQ(ierr);
499 ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
500
501 ierr = VecGetBlockSize(v,&bs);CHKERRQ(ierr);
502 if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
503
504 if (!n) {
505 for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
506 } else {
507 for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);
508
509 for (i=bs; i<n; i+=bs) {
510 for (j=0; j<bs; j++) {
511 if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
512 }
513 }
514 }
515 ierr = MPIU_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);CHKERRQ(ierr);
516
517 ierr = VecRestoreArrayRead(v,&x);CHKERRQ(ierr);
518 PetscFunctionReturn(0);
519 }
520
521 /*@
522 VecStrideMinAll - Computes the minimum of subvector of a vector defined
523 by a starting point and a stride and optionally its location.
524
525 Collective on Vec
526
527 Input Parameter:
528 . v - the vector
529
530 Output Parameter:
531 + idex - the location where the minimum occurred (not supported, pass NULL,
532 if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
533 - nrm - the minimums of each subvector
534
535 Level: advanced
536
537 Notes:
538 One must call VecSetBlockSize() before this routine to set the stride
539 information, or use a vector created from a multicomponent DMDA.
540
541 The dimension of nrm must be the same as the vector block size
542
543
544 .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
545 @*/
VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])546 PetscErrorCode VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
547 {
548 PetscErrorCode ierr;
549 PetscInt i,n,bs,j;
550 const PetscScalar *x;
551 PetscReal min[128],tmp;
552 MPI_Comm comm;
553
554 PetscFunctionBegin;
555 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
556 PetscValidRealPointer(nrm,3);
557 if (idex) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
558 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
559 ierr = VecGetArrayRead(v,&x);CHKERRQ(ierr);
560 ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
561
562 ierr = VecGetBlockSize(v,&bs);CHKERRQ(ierr);
563 if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
564
565 if (!n) {
566 for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
567 } else {
568 for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);
569
570 for (i=bs; i<n; i+=bs) {
571 for (j=0; j<bs; j++) {
572 if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
573 }
574 }
575 }
576 ierr = MPIU_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);CHKERRQ(ierr);
577
578 ierr = VecRestoreArrayRead(v,&x);CHKERRQ(ierr);
579 PetscFunctionReturn(0);
580 }
581
582 /*----------------------------------------------------------------------------------------------*/
583 /*@
584 VecStrideGatherAll - Gathers all the single components from a multi-component vector into
585 separate vectors.
586
587 Collective on Vec
588
589 Input Parameter:
590 + v - the vector
591 - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
592
593 Output Parameter:
594 . s - the location where the subvectors are stored
595
596 Notes:
597 One must call VecSetBlockSize() before this routine to set the stride
598 information, or use a vector created from a multicomponent DMDA.
599
600 If x is the array representing the vector x then this gathers
601 the arrays (x[start],x[start+stride],x[start+2*stride], ....)
602 for start=0,1,2,...bs-1
603
604 The parallel layout of the vector and the subvector must be the same;
605 i.e., nlocal of v = stride*(nlocal of s)
606
607 Not optimized; could be easily
608
609 Level: advanced
610
611 .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
612 VecStrideScatterAll()
613 @*/
VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)614 PetscErrorCode VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
615 {
616 PetscErrorCode ierr;
617 PetscInt i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
618 PetscScalar **y;
619 const PetscScalar *x;
620
621 PetscFunctionBegin;
622 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
623 PetscValidPointer(s,2);
624 PetscValidHeaderSpecific(*s,VEC_CLASSID,2);
625 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
626 ierr = VecGetLocalSize(s[0],&n2);CHKERRQ(ierr);
627 ierr = VecGetArrayRead(v,&x);CHKERRQ(ierr);
628 ierr = VecGetBlockSize(v,&bs);CHKERRQ(ierr);
629 if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
630
631 ierr = PetscMalloc2(bs,&y,bs,&bss);CHKERRQ(ierr);
632 nv = 0;
633 nvc = 0;
634 for (i=0; i<bs; i++) {
635 ierr = VecGetBlockSize(s[i],&bss[i]);CHKERRQ(ierr);
636 if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
637 ierr = VecGetArray(s[i],&y[i]);CHKERRQ(ierr);
638 nvc += bss[i];
639 nv++;
640 if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
641 if (nvc == bs) break;
642 }
643
644 n = n/bs;
645
646 jj = 0;
647 if (addv == INSERT_VALUES) {
648 for (j=0; j<nv; j++) {
649 for (k=0; k<bss[j]; k++) {
650 for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
651 }
652 jj += bss[j];
653 }
654 } else if (addv == ADD_VALUES) {
655 for (j=0; j<nv; j++) {
656 for (k=0; k<bss[j]; k++) {
657 for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
658 }
659 jj += bss[j];
660 }
661 #if !defined(PETSC_USE_COMPLEX)
662 } else if (addv == MAX_VALUES) {
663 for (j=0; j<nv; j++) {
664 for (k=0; k<bss[j]; k++) {
665 for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
666 }
667 jj += bss[j];
668 }
669 #endif
670 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
671
672 ierr = VecRestoreArrayRead(v,&x);CHKERRQ(ierr);
673 for (i=0; i<nv; i++) {
674 ierr = VecRestoreArray(s[i],&y[i]);CHKERRQ(ierr);
675 }
676
677 ierr = PetscFree2(y,bss);CHKERRQ(ierr);
678 PetscFunctionReturn(0);
679 }
680
681 /*@
682 VecStrideScatterAll - Scatters all the single components from separate vectors into
683 a multi-component vector.
684
685 Collective on Vec
686
687 Input Parameter:
688 + s - the location where the subvectors are stored
689 - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
690
691 Output Parameter:
692 . v - the multicomponent vector
693
694 Notes:
695 One must call VecSetBlockSize() before this routine to set the stride
696 information, or use a vector created from a multicomponent DMDA.
697
698 The parallel layout of the vector and the subvector must be the same;
699 i.e., nlocal of v = stride*(nlocal of s)
700
701 Not optimized; could be easily
702
703 Level: advanced
704
705 .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
706 VecStrideScatterAll()
707 @*/
VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)708 PetscErrorCode VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
709 {
710 PetscErrorCode ierr;
711 PetscInt i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
712 PetscScalar *x;
713 PetscScalar const **y;
714
715 PetscFunctionBegin;
716 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
717 PetscValidPointer(s,2);
718 PetscValidHeaderSpecific(*s,VEC_CLASSID,2);
719 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
720 ierr = VecGetLocalSize(s[0],&n2);CHKERRQ(ierr);
721 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
722 ierr = VecGetBlockSize(v,&bs);CHKERRQ(ierr);
723 if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
724
725 ierr = PetscMalloc2(bs,(PetscScalar***)&y,bs,&bss);CHKERRQ(ierr);
726 nv = 0;
727 nvc = 0;
728 for (i=0; i<bs; i++) {
729 ierr = VecGetBlockSize(s[i],&bss[i]);CHKERRQ(ierr);
730 if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1 Re: [PETSC #8241] VecStrideGatherAll */
731 ierr = VecGetArrayRead(s[i],&y[i]);CHKERRQ(ierr);
732 nvc += bss[i];
733 nv++;
734 if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
735 if (nvc == bs) break;
736 }
737
738 n = n/bs;
739
740 jj = 0;
741 if (addv == INSERT_VALUES) {
742 for (j=0; j<nv; j++) {
743 for (k=0; k<bss[j]; k++) {
744 for (i=0; i<n; i++) x[bs*i+jj+k] = y[j][i*bss[j] + k];
745 }
746 jj += bss[j];
747 }
748 } else if (addv == ADD_VALUES) {
749 for (j=0; j<nv; j++) {
750 for (k=0; k<bss[j]; k++) {
751 for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
752 }
753 jj += bss[j];
754 }
755 #if !defined(PETSC_USE_COMPLEX)
756 } else if (addv == MAX_VALUES) {
757 for (j=0; j<nv; j++) {
758 for (k=0; k<bss[j]; k++) {
759 for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
760 }
761 jj += bss[j];
762 }
763 #endif
764 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
765
766 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
767 for (i=0; i<nv; i++) {
768 ierr = VecRestoreArrayRead(s[i],&y[i]);CHKERRQ(ierr);
769 }
770 ierr = PetscFree2(*(PetscScalar***)&y,bss);CHKERRQ(ierr);
771 PetscFunctionReturn(0);
772 }
773
774 /*@
775 VecStrideGather - Gathers a single component from a multi-component vector into
776 another vector.
777
778 Collective on Vec
779
780 Input Parameter:
781 + v - the vector
782 . start - starting point of the subvector (defined by a stride)
783 - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
784
785 Output Parameter:
786 . s - the location where the subvector is stored
787
788 Notes:
789 One must call VecSetBlockSize() before this routine to set the stride
790 information, or use a vector created from a multicomponent DMDA.
791
792 If x is the array representing the vector x then this gathers
793 the array (x[start],x[start+stride],x[start+2*stride], ....)
794
795 The parallel layout of the vector and the subvector must be the same;
796 i.e., nlocal of v = stride*(nlocal of s)
797
798 Not optimized; could be easily
799
800 Level: advanced
801
802 .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
803 VecStrideScatterAll()
804 @*/
VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)805 PetscErrorCode VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
806 {
807 PetscErrorCode ierr;
808
809 PetscFunctionBegin;
810 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
811 PetscValidHeaderSpecific(s,VEC_CLASSID,3);
812 if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
813 if (start >= v->map->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,v->map->bs);
814 if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
815 ierr = (*v->ops->stridegather)(v,start,s,addv);CHKERRQ(ierr);
816 PetscFunctionReturn(0);
817 }
818
819 /*@
820 VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
821
822 Collective on Vec
823
824 Input Parameter:
825 + s - the single-component vector
826 . start - starting point of the subvector (defined by a stride)
827 - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
828
829 Output Parameter:
830 . v - the location where the subvector is scattered (the multi-component vector)
831
832 Notes:
833 One must call VecSetBlockSize() on the multi-component vector before this
834 routine to set the stride information, or use a vector created from a multicomponent DMDA.
835
836 The parallel layout of the vector and the subvector must be the same;
837 i.e., nlocal of v = stride*(nlocal of s)
838
839 Not optimized; could be easily
840
841 Level: advanced
842
843 .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
844 VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
845 @*/
VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)846 PetscErrorCode VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
847 {
848 PetscErrorCode ierr;
849
850 PetscFunctionBegin;
851 PetscValidHeaderSpecific(s,VEC_CLASSID,1);
852 PetscValidHeaderSpecific(v,VEC_CLASSID,3);
853 if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
854 if (start >= v->map->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,v->map->bs);
855 if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
856 ierr = (*v->ops->stridescatter)(s,start,v,addv);CHKERRQ(ierr);
857 PetscFunctionReturn(0);
858 }
859
860 /*@
861 VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
862 another vector.
863
864 Collective on Vec
865
866 Input Parameter:
867 + v - the vector
868 . nidx - the number of indices
869 . idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
870 . idxs - the indices of the components 0 <= idxs[0] ...idxs[nidx-1] < bs(s), they need not be sorted, may be null if nidx == bs(s) or is PETSC_DETERMINE
871 - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
872
873 Output Parameter:
874 . s - the location where the subvector is stored
875
876 Notes:
877 One must call VecSetBlockSize() on both vectors before this routine to set the stride
878 information, or use a vector created from a multicomponent DMDA.
879
880
881 The parallel layout of the vector and the subvector must be the same;
882
883 Not optimized; could be easily
884
885 Level: advanced
886
887 .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
888 VecStrideScatterAll()
889 @*/
VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)890 PetscErrorCode VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
891 {
892 PetscErrorCode ierr;
893
894 PetscFunctionBegin;
895 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
896 PetscValidHeaderSpecific(s,VEC_CLASSID,5);
897 if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
898 if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
899 ierr = (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);CHKERRQ(ierr);
900 PetscFunctionReturn(0);
901 }
902
903 /*@
904 VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.
905
906 Collective on Vec
907
908 Input Parameter:
909 + s - the smaller-component vector
910 . nidx - the number of indices in idx
911 . idxs - the indices of the components in the smaller-component vector, 0 <= idxs[0] ...idxs[nidx-1] < bs(s) they need not be sorted, may be null if nidx == bs(s) or is PETSC_DETERMINE
912 . idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
913 - addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
914
915 Output Parameter:
916 . v - the location where the subvector is into scattered (the multi-component vector)
917
918 Notes:
919 One must call VecSetBlockSize() on the vectors before this
920 routine to set the stride information, or use a vector created from a multicomponent DMDA.
921
922 The parallel layout of the vector and the subvector must be the same;
923
924 Not optimized; could be easily
925
926 Level: advanced
927
928 .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
929 VecStrideScatterAll()
930 @*/
VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)931 PetscErrorCode VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
932 {
933 PetscErrorCode ierr;
934
935 PetscFunctionBegin;
936 PetscValidHeaderSpecific(s,VEC_CLASSID,1);
937 PetscValidHeaderSpecific(v,VEC_CLASSID,5);
938 if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
939 if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
940 ierr = (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);CHKERRQ(ierr);
941 PetscFunctionReturn(0);
942 }
943
VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)944 PetscErrorCode VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
945 {
946 PetscErrorCode ierr;
947 PetscInt i,n,bs,ns;
948 const PetscScalar *x;
949 PetscScalar *y;
950
951 PetscFunctionBegin;
952 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
953 ierr = VecGetLocalSize(s,&ns);CHKERRQ(ierr);
954 ierr = VecGetArrayRead(v,&x);CHKERRQ(ierr);
955 ierr = VecGetArray(s,&y);CHKERRQ(ierr);
956
957 bs = v->map->bs;
958 if (n != ns*bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for gather from original vector %D",ns*bs,n);
959 x += start;
960 n = n/bs;
961
962 if (addv == INSERT_VALUES) {
963 for (i=0; i<n; i++) y[i] = x[bs*i];
964 } else if (addv == ADD_VALUES) {
965 for (i=0; i<n; i++) y[i] += x[bs*i];
966 #if !defined(PETSC_USE_COMPLEX)
967 } else if (addv == MAX_VALUES) {
968 for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
969 #endif
970 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
971
972 ierr = VecRestoreArrayRead(v,&x);CHKERRQ(ierr);
973 ierr = VecRestoreArray(s,&y);CHKERRQ(ierr);
974 PetscFunctionReturn(0);
975 }
976
VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)977 PetscErrorCode VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
978 {
979 PetscErrorCode ierr;
980 PetscInt i,n,bs,ns;
981 PetscScalar *x;
982 const PetscScalar *y;
983
984 PetscFunctionBegin;
985 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
986 ierr = VecGetLocalSize(s,&ns);CHKERRQ(ierr);
987 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
988 ierr = VecGetArrayRead(s,&y);CHKERRQ(ierr);
989
990 bs = v->map->bs;
991 if (n != ns*bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for scatter to multicomponent vector %D",ns*bs,n);
992 x += start;
993 n = n/bs;
994
995 if (addv == INSERT_VALUES) {
996 for (i=0; i<n; i++) x[bs*i] = y[i];
997 } else if (addv == ADD_VALUES) {
998 for (i=0; i<n; i++) x[bs*i] += y[i];
999 #if !defined(PETSC_USE_COMPLEX)
1000 } else if (addv == MAX_VALUES) {
1001 for (i=0; i<n; i++) x[bs*i] = PetscMax(y[i],x[bs*i]);
1002 #endif
1003 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1004
1005 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
1006 ierr = VecRestoreArrayRead(s,&y);CHKERRQ(ierr);
1007 PetscFunctionReturn(0);
1008 }
1009
VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)1010 PetscErrorCode VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
1011 {
1012 PetscErrorCode ierr;
1013 PetscInt i,j,n,bs,bss,ns;
1014 const PetscScalar *x;
1015 PetscScalar *y;
1016
1017 PetscFunctionBegin;
1018 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
1019 ierr = VecGetLocalSize(s,&ns);CHKERRQ(ierr);
1020 ierr = VecGetArrayRead(v,&x);CHKERRQ(ierr);
1021 ierr = VecGetArray(s,&y);CHKERRQ(ierr);
1022
1023 bs = v->map->bs;
1024 bss = s->map->bs;
1025 n = n/bs;
1026
1027 if (PetscDefined(USE_DEBUG)) {
1028 if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1029 for (j=0; j<nidx; j++) {
1030 if (idxv[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
1031 if (idxv[j] >= bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxv[j],bs);
1032 }
1033 if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
1034 }
1035
1036 if (addv == INSERT_VALUES) {
1037 if (!idxs) {
1038 for (i=0; i<n; i++) {
1039 for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
1040 }
1041 } else {
1042 for (i=0; i<n; i++) {
1043 for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
1044 }
1045 }
1046 } else if (addv == ADD_VALUES) {
1047 if (!idxs) {
1048 for (i=0; i<n; i++) {
1049 for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
1050 }
1051 } else {
1052 for (i=0; i<n; i++) {
1053 for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
1054 }
1055 }
1056 #if !defined(PETSC_USE_COMPLEX)
1057 } else if (addv == MAX_VALUES) {
1058 if (!idxs) {
1059 for (i=0; i<n; i++) {
1060 for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(y[bss*i+j],x[bs*i+idxv[j]]);
1061 }
1062 } else {
1063 for (i=0; i<n; i++) {
1064 for (j=0; j<bss; j++) y[bss*i+idxs[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
1065 }
1066 }
1067 #endif
1068 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1069
1070 ierr = VecRestoreArrayRead(v,&x);CHKERRQ(ierr);
1071 ierr = VecRestoreArray(s,&y);CHKERRQ(ierr);
1072 PetscFunctionReturn(0);
1073 }
1074
VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)1075 PetscErrorCode VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
1076 {
1077 PetscErrorCode ierr;
1078 PetscInt j,i,n,bs,ns,bss;
1079 PetscScalar *x;
1080 const PetscScalar *y;
1081
1082 PetscFunctionBegin;
1083 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
1084 ierr = VecGetLocalSize(s,&ns);CHKERRQ(ierr);
1085 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
1086 ierr = VecGetArrayRead(s,&y);CHKERRQ(ierr);
1087
1088 bs = v->map->bs;
1089 bss = s->map->bs;
1090 n = n/bs;
1091
1092 if (PetscDefined(USE_DEBUG)) {
1093 if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
1094 for (j=0; j<bss; j++) {
1095 if (idxs) {
1096 if (idxs[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxs[j]);
1097 if (idxs[j] >= bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxs[j],bs);
1098 }
1099 }
1100 if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
1101 }
1102
1103 if (addv == INSERT_VALUES) {
1104 if (!idxs) {
1105 for (i=0; i<n; i++) {
1106 for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
1107 }
1108 } else {
1109 for (i=0; i<n; i++) {
1110 for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
1111 }
1112 }
1113 } else if (addv == ADD_VALUES) {
1114 if (!idxs) {
1115 for (i=0; i<n; i++) {
1116 for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
1117 }
1118 } else {
1119 for (i=0; i<n; i++) {
1120 for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
1121 }
1122 }
1123 #if !defined(PETSC_USE_COMPLEX)
1124 } else if (addv == MAX_VALUES) {
1125 if (!idxs) {
1126 for (i=0; i<n; i++) {
1127 for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
1128 }
1129 } else {
1130 for (i=0; i<n; i++) {
1131 for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
1132 }
1133 }
1134 #endif
1135 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
1136
1137 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
1138 ierr = VecRestoreArrayRead(s,&y);CHKERRQ(ierr);
1139 PetscFunctionReturn(0);
1140 }
1141
VecReciprocal_Default(Vec v)1142 PetscErrorCode VecReciprocal_Default(Vec v)
1143 {
1144 PetscErrorCode ierr;
1145 PetscInt i,n;
1146 PetscScalar *x;
1147
1148 PetscFunctionBegin;
1149 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
1150 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
1151 for (i=0; i<n; i++) {
1152 if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
1153 }
1154 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
1155 PetscFunctionReturn(0);
1156 }
1157
1158 /*@
1159 VecExp - Replaces each component of a vector by e^x_i
1160
1161 Not collective
1162
1163 Input Parameter:
1164 . v - The vector
1165
1166 Output Parameter:
1167 . v - The vector of exponents
1168
1169 Level: beginner
1170
1171 .seealso: VecLog(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1172
1173 @*/
VecExp(Vec v)1174 PetscErrorCode VecExp(Vec v)
1175 {
1176 PetscScalar *x;
1177 PetscInt i, n;
1178 PetscErrorCode ierr;
1179
1180 PetscFunctionBegin;
1181 PetscValidHeaderSpecific(v, VEC_CLASSID,1);
1182 if (v->ops->exp) {
1183 ierr = (*v->ops->exp)(v);CHKERRQ(ierr);
1184 } else {
1185 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr);
1186 ierr = VecGetArray(v, &x);CHKERRQ(ierr);
1187 for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
1188 ierr = VecRestoreArray(v, &x);CHKERRQ(ierr);
1189 }
1190 PetscFunctionReturn(0);
1191 }
1192
1193 /*@
1194 VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
1195
1196 Not collective
1197
1198 Input Parameter:
1199 . v - The vector
1200
1201 Output Parameter:
1202 . v - The vector of logs
1203
1204 Level: beginner
1205
1206 .seealso: VecExp(), VecAbs(), VecSqrtAbs(), VecReciprocal()
1207
1208 @*/
VecLog(Vec v)1209 PetscErrorCode VecLog(Vec v)
1210 {
1211 PetscScalar *x;
1212 PetscInt i, n;
1213 PetscErrorCode ierr;
1214
1215 PetscFunctionBegin;
1216 PetscValidHeaderSpecific(v, VEC_CLASSID,1);
1217 if (v->ops->log) {
1218 ierr = (*v->ops->log)(v);CHKERRQ(ierr);
1219 } else {
1220 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr);
1221 ierr = VecGetArray(v, &x);CHKERRQ(ierr);
1222 for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
1223 ierr = VecRestoreArray(v, &x);CHKERRQ(ierr);
1224 }
1225 PetscFunctionReturn(0);
1226 }
1227
1228 /*@
1229 VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.
1230
1231 Not collective
1232
1233 Input Parameter:
1234 . v - The vector
1235
1236 Output Parameter:
1237 . v - The vector square root
1238
1239 Level: beginner
1240
1241 Note: The actual function is sqrt(|x_i|)
1242
1243 .seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()
1244
1245 @*/
VecSqrtAbs(Vec v)1246 PetscErrorCode VecSqrtAbs(Vec v)
1247 {
1248 PetscScalar *x;
1249 PetscInt i, n;
1250 PetscErrorCode ierr;
1251
1252 PetscFunctionBegin;
1253 PetscValidHeaderSpecific(v, VEC_CLASSID,1);
1254 if (v->ops->sqrt) {
1255 ierr = (*v->ops->sqrt)(v);CHKERRQ(ierr);
1256 } else {
1257 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr);
1258 ierr = VecGetArray(v, &x);CHKERRQ(ierr);
1259 for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
1260 ierr = VecRestoreArray(v, &x);CHKERRQ(ierr);
1261 }
1262 PetscFunctionReturn(0);
1263 }
1264
1265 /*@
1266 VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
1267
1268 Collective on Vec
1269
1270 Input Parameter:
1271 + s - first vector
1272 - t - second vector
1273
1274 Output Parameter:
1275 + dp - s'conj(t)
1276 - nm - t'conj(t)
1277
1278 Level: advanced
1279
1280 Notes:
1281 conj(x) is the complex conjugate of x when x is complex
1282
1283
1284 .seealso: VecDot(), VecNorm(), VecDotBegin(), VecNormBegin(), VecDotEnd(), VecNormEnd()
1285
1286 @*/
VecDotNorm2(Vec s,Vec t,PetscScalar * dp,PetscReal * nm)1287 PetscErrorCode VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
1288 {
1289 const PetscScalar *sx, *tx;
1290 PetscScalar dpx = 0.0, nmx = 0.0,work[2],sum[2];
1291 PetscInt i, n;
1292 PetscErrorCode ierr;
1293
1294 PetscFunctionBegin;
1295 PetscValidHeaderSpecific(s, VEC_CLASSID,1);
1296 PetscValidHeaderSpecific(t, VEC_CLASSID,2);
1297 PetscValidScalarPointer(dp,3);
1298 PetscValidScalarPointer(nm,4);
1299 PetscValidType(s,1);
1300 PetscValidType(t,2);
1301 PetscCheckSameTypeAndComm(s,1,t,2);
1302 if (s->map->N != t->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1303 if (s->map->n != t->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1304
1305 ierr = PetscLogEventBegin(VEC_DotNorm2,s,t,0,0);CHKERRQ(ierr);
1306 if (s->ops->dotnorm2) {
1307 ierr = (*s->ops->dotnorm2)(s,t,dp,&dpx);CHKERRQ(ierr);
1308 *nm = PetscRealPart(dpx);
1309 } else {
1310 ierr = VecGetLocalSize(s, &n);CHKERRQ(ierr);
1311 ierr = VecGetArrayRead(s, &sx);CHKERRQ(ierr);
1312 ierr = VecGetArrayRead(t, &tx);CHKERRQ(ierr);
1313
1314 for (i = 0; i<n; i++) {
1315 dpx += sx[i]*PetscConj(tx[i]);
1316 nmx += tx[i]*PetscConj(tx[i]);
1317 }
1318 work[0] = dpx;
1319 work[1] = nmx;
1320
1321 ierr = MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));CHKERRQ(ierr);
1322 *dp = sum[0];
1323 *nm = PetscRealPart(sum[1]);
1324
1325 ierr = VecRestoreArrayRead(t, &tx);CHKERRQ(ierr);
1326 ierr = VecRestoreArrayRead(s, &sx);CHKERRQ(ierr);
1327 ierr = PetscLogFlops(4.0*n);CHKERRQ(ierr);
1328 }
1329 ierr = PetscLogEventEnd(VEC_DotNorm2,s,t,0,0);CHKERRQ(ierr);
1330 PetscFunctionReturn(0);
1331 }
1332
1333 /*@
1334 VecSum - Computes the sum of all the components of a vector.
1335
1336 Collective on Vec
1337
1338 Input Parameter:
1339 . v - the vector
1340
1341 Output Parameter:
1342 . sum - the result
1343
1344 Level: beginner
1345
1346 .seealso: VecNorm()
1347 @*/
VecSum(Vec v,PetscScalar * sum)1348 PetscErrorCode VecSum(Vec v,PetscScalar *sum)
1349 {
1350 PetscErrorCode ierr;
1351 PetscInt i,n;
1352 const PetscScalar *x;
1353 PetscScalar lsum = 0.0;
1354
1355 PetscFunctionBegin;
1356 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
1357 PetscValidScalarPointer(sum,2);
1358 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
1359 ierr = VecGetArrayRead(v,&x);CHKERRQ(ierr);
1360 for (i=0; i<n; i++) lsum += x[i];
1361 ierr = MPIU_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));CHKERRQ(ierr);
1362 ierr = VecRestoreArrayRead(v,&x);CHKERRQ(ierr);
1363 PetscFunctionReturn(0);
1364 }
1365
1366 /*@
1367 VecImaginaryPart - Replaces a complex vector with its imginary part
1368
1369 Collective on Vec
1370
1371 Input Parameter:
1372 . v - the vector
1373
1374 Level: beginner
1375
1376 .seealso: VecNorm(), VecRealPart()
1377 @*/
VecImaginaryPart(Vec v)1378 PetscErrorCode VecImaginaryPart(Vec v)
1379 {
1380 PetscErrorCode ierr;
1381 PetscInt i,n;
1382 PetscScalar *x;
1383
1384 PetscFunctionBegin;
1385 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
1386 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
1387 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
1388 for (i=0; i<n; i++) x[i] = PetscImaginaryPart(x[i]);
1389 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
1390 PetscFunctionReturn(0);
1391 }
1392
1393 /*@
1394 VecRealPart - Replaces a complex vector with its real part
1395
1396 Collective on Vec
1397
1398 Input Parameter:
1399 . v - the vector
1400
1401 Level: beginner
1402
1403 .seealso: VecNorm(), VecImaginaryPart()
1404 @*/
VecRealPart(Vec v)1405 PetscErrorCode VecRealPart(Vec v)
1406 {
1407 PetscErrorCode ierr;
1408 PetscInt i,n;
1409 PetscScalar *x;
1410
1411 PetscFunctionBegin;
1412 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
1413 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
1414 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
1415 for (i=0; i<n; i++) x[i] = PetscRealPart(x[i]);
1416 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
1417 PetscFunctionReturn(0);
1418 }
1419
1420 /*@
1421 VecShift - Shifts all of the components of a vector by computing
1422 x[i] = x[i] + shift.
1423
1424 Logically Collective on Vec
1425
1426 Input Parameters:
1427 + v - the vector
1428 - shift - the shift
1429
1430 Output Parameter:
1431 . v - the shifted vector
1432
1433 Level: intermediate
1434
1435 @*/
VecShift(Vec v,PetscScalar shift)1436 PetscErrorCode VecShift(Vec v,PetscScalar shift)
1437 {
1438 PetscErrorCode ierr;
1439 PetscInt i,n;
1440 PetscScalar *x;
1441
1442 PetscFunctionBegin;
1443 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
1444 PetscValidLogicalCollectiveScalar(v,shift,2);
1445 ierr = VecSetErrorIfLocked(v,1);CHKERRQ(ierr);
1446 if (shift == 0.0) PetscFunctionReturn(0);
1447
1448 if (v->ops->shift) {
1449 ierr = (*v->ops->shift)(v,shift);CHKERRQ(ierr);
1450 } else {
1451 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
1452 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
1453 for (i=0; i<n; i++) x[i] += shift;
1454 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
1455 }
1456 PetscFunctionReturn(0);
1457 }
1458
1459 /*@
1460 VecAbs - Replaces every element in a vector with its absolute value.
1461
1462 Logically Collective on Vec
1463
1464 Input Parameters:
1465 . v - the vector
1466
1467 Level: intermediate
1468
1469 @*/
VecAbs(Vec v)1470 PetscErrorCode VecAbs(Vec v)
1471 {
1472 PetscErrorCode ierr;
1473 PetscInt i,n;
1474 PetscScalar *x;
1475
1476 PetscFunctionBegin;
1477 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
1478 ierr = VecSetErrorIfLocked(v,1);CHKERRQ(ierr);
1479
1480 if (v->ops->abs) {
1481 ierr = (*v->ops->abs)(v);CHKERRQ(ierr);
1482 } else {
1483 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
1484 ierr = VecGetArray(v,&x);CHKERRQ(ierr);
1485 for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
1486 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
1487 }
1488 PetscFunctionReturn(0);
1489 }
1490
1491 /*@
1492 VecPermute - Permutes a vector in place using the given ordering.
1493
1494 Input Parameters:
1495 + vec - The vector
1496 . order - The ordering
1497 - inv - The flag for inverting the permutation
1498
1499 Level: beginner
1500
1501 Note: This function does not yet support parallel Index Sets with non-local permutations
1502
1503 .seealso: MatPermute()
1504 @*/
VecPermute(Vec x,IS row,PetscBool inv)1505 PetscErrorCode VecPermute(Vec x, IS row, PetscBool inv)
1506 {
1507 const PetscScalar *array;
1508 PetscScalar *newArray;
1509 const PetscInt *idx;
1510 PetscInt i,rstart,rend;
1511 PetscErrorCode ierr;
1512
1513 PetscFunctionBegin;
1514 ierr = VecSetErrorIfLocked(x,1);CHKERRQ(ierr);
1515
1516 ierr = VecGetOwnershipRange(x,&rstart,&rend);CHKERRQ(ierr);
1517 ierr = ISGetIndices(row, &idx);CHKERRQ(ierr);
1518 ierr = VecGetArrayRead(x, &array);CHKERRQ(ierr);
1519 ierr = PetscMalloc1(x->map->n, &newArray);CHKERRQ(ierr);
1520 if (PetscDefined(USE_DEBUG)) {
1521 for (i = 0; i < x->map->n; i++) {
1522 if ((idx[i] < rstart) || (idx[i] >= rend)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT, "Permutation index %D is out of bounds: %D", i, idx[i]);
1523 }
1524 }
1525 if (!inv) {
1526 for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
1527 } else {
1528 for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
1529 }
1530 ierr = VecRestoreArrayRead(x, &array);CHKERRQ(ierr);
1531 ierr = ISRestoreIndices(row, &idx);CHKERRQ(ierr);
1532 ierr = VecReplaceArray(x, newArray);CHKERRQ(ierr);
1533 PetscFunctionReturn(0);
1534 }
1535
1536 /*@
1537 VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
1538 or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
1539 Does NOT take round-off errors into account.
1540
1541 Collective on Vec
1542
1543 Input Parameters:
1544 + vec1 - the first vector
1545 - vec2 - the second vector
1546
1547 Output Parameter:
1548 . flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
1549
1550 Level: intermediate
1551
1552
1553 @*/
VecEqual(Vec vec1,Vec vec2,PetscBool * flg)1554 PetscErrorCode VecEqual(Vec vec1,Vec vec2,PetscBool *flg)
1555 {
1556 const PetscScalar *v1,*v2;
1557 PetscErrorCode ierr;
1558 PetscInt n1,n2,N1,N2;
1559 PetscBool flg1;
1560
1561 PetscFunctionBegin;
1562 PetscValidHeaderSpecific(vec1,VEC_CLASSID,1);
1563 PetscValidHeaderSpecific(vec2,VEC_CLASSID,2);
1564 PetscValidBoolPointer(flg,3);
1565 if (vec1 == vec2) *flg = PETSC_TRUE;
1566 else {
1567 ierr = VecGetSize(vec1,&N1);CHKERRQ(ierr);
1568 ierr = VecGetSize(vec2,&N2);CHKERRQ(ierr);
1569 if (N1 != N2) flg1 = PETSC_FALSE;
1570 else {
1571 ierr = VecGetLocalSize(vec1,&n1);CHKERRQ(ierr);
1572 ierr = VecGetLocalSize(vec2,&n2);CHKERRQ(ierr);
1573 if (n1 != n2) flg1 = PETSC_FALSE;
1574 else {
1575 ierr = VecGetArrayRead(vec1,&v1);CHKERRQ(ierr);
1576 ierr = VecGetArrayRead(vec2,&v2);CHKERRQ(ierr);
1577 ierr = PetscArraycmp(v1,v2,n1,&flg1);CHKERRQ(ierr);
1578 ierr = VecRestoreArrayRead(vec1,&v1);CHKERRQ(ierr);
1579 ierr = VecRestoreArrayRead(vec2,&v2);CHKERRQ(ierr);
1580 }
1581 }
1582 /* combine results from all processors */
1583 ierr = MPIU_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));CHKERRQ(ierr);
1584 }
1585 PetscFunctionReturn(0);
1586 }
1587
1588 /*@
1589 VecUniqueEntries - Compute the number of unique entries, and those entries
1590
1591 Collective on Vec
1592
1593 Input Parameter:
1594 . vec - the vector
1595
1596 Output Parameters:
1597 + n - The number of unique entries
1598 - e - The entries
1599
1600 Level: intermediate
1601
1602 @*/
VecUniqueEntries(Vec vec,PetscInt * n,PetscScalar ** e)1603 PetscErrorCode VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
1604 {
1605 const PetscScalar *v;
1606 PetscScalar *tmp, *vals;
1607 PetscMPIInt *N, *displs, l;
1608 PetscInt ng, m, i, j, p;
1609 PetscMPIInt size;
1610 PetscErrorCode ierr;
1611
1612 PetscFunctionBegin;
1613 PetscValidHeaderSpecific(vec,VEC_CLASSID,1);
1614 PetscValidIntPointer(n,2);
1615 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) vec), &size);CHKERRQ(ierr);
1616 ierr = VecGetLocalSize(vec, &m);CHKERRQ(ierr);
1617 ierr = VecGetArrayRead(vec, &v);CHKERRQ(ierr);
1618 ierr = PetscMalloc2(m,&tmp,size,&N);CHKERRQ(ierr);
1619 for (i = 0, j = 0, l = 0; i < m; ++i) {
1620 /* Can speed this up with sorting */
1621 for (j = 0; j < l; ++j) {
1622 if (v[i] == tmp[j]) break;
1623 }
1624 if (j == l) {
1625 tmp[j] = v[i];
1626 ++l;
1627 }
1628 }
1629 ierr = VecRestoreArrayRead(vec, &v);CHKERRQ(ierr);
1630 /* Gather serial results */
1631 ierr = MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject) vec));CHKERRQ(ierr);
1632 for (p = 0, ng = 0; p < size; ++p) {
1633 ng += N[p];
1634 }
1635 ierr = PetscMalloc2(ng,&vals,size+1,&displs);CHKERRQ(ierr);
1636 for (p = 1, displs[0] = 0; p <= size; ++p) {
1637 displs[p] = displs[p-1] + N[p-1];
1638 }
1639 ierr = MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject) vec));CHKERRQ(ierr);
1640 /* Find unique entries */
1641 #ifdef PETSC_USE_COMPLEX
1642 SETERRQ(PetscObjectComm((PetscObject) vec), PETSC_ERR_SUP, "Does not work with complex numbers");
1643 #else
1644 *n = displs[size];
1645 ierr = PetscSortRemoveDupsReal(n, (PetscReal *) vals);CHKERRQ(ierr);
1646 if (e) {
1647 PetscValidPointer(e,3);
1648 ierr = PetscMalloc1(*n, e);CHKERRQ(ierr);
1649 for (i = 0; i < *n; ++i) {
1650 (*e)[i] = vals[i];
1651 }
1652 }
1653 ierr = PetscFree2(vals,displs);CHKERRQ(ierr);
1654 ierr = PetscFree2(tmp,N);CHKERRQ(ierr);
1655 PetscFunctionReturn(0);
1656 #endif
1657 }
1658
1659
1660
1661