1 #include <petsc/private/vecimpl.h> /*I "petscvec.h" I*/
2
3 /*@
4 VecWhichEqual - Creates an index set containing the indices
5 where the vectors Vec1 and Vec2 have identical elements.
6
7 Collective on Vec
8
9 Input Parameters:
10 . Vec1, Vec2 - the two vectors to compare
11
12 OutputParameter:
13 . S - The index set containing the indices i where vec1[i] == vec2[i]
14
15 Notes:
16 the two vectors must have the same parallel layout
17
18 Level: advanced
19 @*/
VecWhichEqual(Vec Vec1,Vec Vec2,IS * S)20 PetscErrorCode VecWhichEqual(Vec Vec1, Vec Vec2, IS *S)
21 {
22 PetscErrorCode ierr;
23 PetscInt i,n_same=0;
24 PetscInt n,low,high;
25 PetscInt *same=NULL;
26 const PetscScalar *v1,*v2;
27
28 PetscFunctionBegin;
29 PetscValidHeaderSpecific(Vec1,VEC_CLASSID,1);
30 PetscValidHeaderSpecific(Vec2,VEC_CLASSID,2);
31 PetscCheckSameComm(Vec1,1,Vec2,2);
32 VecCheckSameSize(Vec1,1,Vec2,2);
33
34 ierr = VecGetOwnershipRange(Vec1,&low,&high);CHKERRQ(ierr);
35 ierr = VecGetLocalSize(Vec1,&n);CHKERRQ(ierr);
36 if (n>0){
37 if (Vec1 == Vec2){
38 ierr = VecGetArrayRead(Vec1,&v1);CHKERRQ(ierr);
39 v2=v1;
40 } else {
41 ierr = VecGetArrayRead(Vec1,&v1);CHKERRQ(ierr);
42 ierr = VecGetArrayRead(Vec2,&v2);CHKERRQ(ierr);
43 }
44
45 ierr = PetscMalloc1(n,&same);CHKERRQ(ierr);
46
47 for (i=0; i<n; ++i){
48 if (v1[i] == v2[i]) {same[n_same]=low+i; ++n_same;}
49 }
50
51 if (Vec1 == Vec2){
52 ierr = VecRestoreArrayRead(Vec1,&v1);CHKERRQ(ierr);
53 } else {
54 ierr = VecRestoreArrayRead(Vec1,&v1);CHKERRQ(ierr);
55 ierr = VecRestoreArrayRead(Vec2,&v2);CHKERRQ(ierr);
56 }
57 }
58 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)Vec1),n_same,same,PETSC_OWN_POINTER,S);CHKERRQ(ierr);
59 PetscFunctionReturn(0);
60 }
61
62 /*@
63 VecWhichLessThan - Creates an index set containing the indices
64 where the vectors Vec1 < Vec2
65
66 Collective on S
67
68 Input Parameters:
69 . Vec1, Vec2 - the two vectors to compare
70
71 OutputParameter:
72 . S - The index set containing the indices i where vec1[i] < vec2[i]
73
74 Notes:
75 The two vectors must have the same parallel layout
76
77 For complex numbers this only compares the real part
78
79 Level: advanced
80 @*/
VecWhichLessThan(Vec Vec1,Vec Vec2,IS * S)81 PetscErrorCode VecWhichLessThan(Vec Vec1, Vec Vec2, IS *S)
82 {
83 PetscErrorCode ierr;
84 PetscInt i,n_lt=0;
85 PetscInt n,low,high;
86 PetscInt *lt=NULL;
87 const PetscScalar *v1,*v2;
88
89 PetscFunctionBegin;
90 PetscValidHeaderSpecific(Vec1,VEC_CLASSID,1);
91 PetscValidHeaderSpecific(Vec2,VEC_CLASSID,2);
92 PetscCheckSameComm(Vec1,1,Vec2,2);
93 VecCheckSameSize(Vec1,1,Vec2,2);
94
95 ierr = VecGetOwnershipRange(Vec1,&low,&high);CHKERRQ(ierr);
96 ierr = VecGetLocalSize(Vec1,&n);CHKERRQ(ierr);
97 if (n>0){
98 if (Vec1 == Vec2){
99 ierr = VecGetArrayRead(Vec1,&v1);CHKERRQ(ierr);
100 v2=v1;
101 } else {
102 ierr = VecGetArrayRead(Vec1,&v1);CHKERRQ(ierr);
103 ierr = VecGetArrayRead(Vec2,&v2);CHKERRQ(ierr);
104 }
105
106 ierr = PetscMalloc1(n,<);CHKERRQ(ierr);
107
108 for (i=0; i<n; ++i){
109 if (PetscRealPart(v1[i]) < PetscRealPart(v2[i])) {lt[n_lt]=low+i; ++n_lt;}
110 }
111
112 if (Vec1 == Vec2){
113 ierr = VecRestoreArrayRead(Vec1,&v1);CHKERRQ(ierr);
114 } else {
115 ierr = VecRestoreArrayRead(Vec1,&v1);CHKERRQ(ierr);
116 ierr = VecRestoreArrayRead(Vec2,&v2);CHKERRQ(ierr);
117 }
118 }
119 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)Vec1),n_lt,lt,PETSC_OWN_POINTER,S);CHKERRQ(ierr);
120 PetscFunctionReturn(0);
121 }
122
123 /*@
124 VecWhichGreaterThan - Creates an index set containing the indices
125 where the vectors Vec1 > Vec2
126
127 Collective on S
128
129 Input Parameters:
130 . Vec1, Vec2 - the two vectors to compare
131
132 OutputParameter:
133 . S - The index set containing the indices i where vec1[i] > vec2[i]
134
135 Notes:
136 The two vectors must have the same parallel layout
137
138 For complex numbers this only compares the real part
139
140 Level: advanced
141 @*/
VecWhichGreaterThan(Vec Vec1,Vec Vec2,IS * S)142 PetscErrorCode VecWhichGreaterThan(Vec Vec1, Vec Vec2, IS *S)
143 {
144 PetscErrorCode ierr;
145 PetscInt i,n_gt=0;
146 PetscInt n,low,high;
147 PetscInt *gt=NULL;
148 const PetscScalar *v1,*v2;
149
150 PetscFunctionBegin;
151 PetscValidHeaderSpecific(Vec1,VEC_CLASSID,1);
152 PetscValidHeaderSpecific(Vec2,VEC_CLASSID,2);
153 PetscCheckSameComm(Vec1,1,Vec2,2);
154 VecCheckSameSize(Vec1,1,Vec2,2);
155
156 ierr = VecGetOwnershipRange(Vec1,&low,&high);CHKERRQ(ierr);
157 ierr = VecGetLocalSize(Vec1,&n);CHKERRQ(ierr);
158 if (n>0){
159 if (Vec1 == Vec2){
160 ierr = VecGetArrayRead(Vec1,&v1);CHKERRQ(ierr);
161 v2=v1;
162 } else {
163 ierr = VecGetArrayRead(Vec1,&v1);CHKERRQ(ierr);
164 ierr = VecGetArrayRead(Vec2,&v2);CHKERRQ(ierr);
165 }
166
167 ierr = PetscMalloc1(n,>);CHKERRQ(ierr);
168
169 for (i=0; i<n; ++i){
170 if (PetscRealPart(v1[i]) > PetscRealPart(v2[i])) {gt[n_gt]=low+i; ++n_gt;}
171 }
172
173 if (Vec1 == Vec2){
174 ierr = VecRestoreArrayRead(Vec1,&v1);CHKERRQ(ierr);
175 } else {
176 ierr = VecRestoreArrayRead(Vec1,&v1);CHKERRQ(ierr);
177 ierr = VecRestoreArrayRead(Vec2,&v2);CHKERRQ(ierr);
178 }
179 }
180 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)Vec1),n_gt,gt,PETSC_OWN_POINTER,S);CHKERRQ(ierr);
181 PetscFunctionReturn(0);
182 }
183
184 /*@
185 VecWhichBetween - Creates an index set containing the indices
186 where VecLow < V < VecHigh
187
188 Collective on S
189
190 Input Parameters:
191 + VecLow - lower bound
192 . V - Vector to compare
193 - VecHigh - higher bound
194
195 OutputParameter:
196 . S - The index set containing the indices i where veclow[i] < v[i] < vechigh[i]
197
198 Notes:
199 The vectors must have the same parallel layout
200
201 For complex numbers this only compares the real part
202
203 Level: advanced
204 @*/
VecWhichBetween(Vec VecLow,Vec V,Vec VecHigh,IS * S)205 PetscErrorCode VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS *S)
206 {
207
208 PetscErrorCode ierr;
209 PetscInt i,n_vm=0;
210 PetscInt n,low,high;
211 PetscInt *vm=NULL;
212 const PetscScalar *v1,*v2,*vmiddle;
213
214 PetscFunctionBegin;
215 PetscValidHeaderSpecific(VecLow,VEC_CLASSID,1);
216 PetscValidHeaderSpecific(V,VEC_CLASSID,2);
217 PetscValidHeaderSpecific(VecHigh,VEC_CLASSID,3);
218 PetscCheckSameComm(V,2,VecLow,1);
219 PetscCheckSameComm(V,2,VecHigh,3);
220 VecCheckSameSize(V,2,VecLow,1);
221 VecCheckSameSize(V,2,VecHigh,3);
222
223 ierr = VecGetOwnershipRange(VecLow,&low,&high);CHKERRQ(ierr);
224 ierr = VecGetLocalSize(VecLow,&n);CHKERRQ(ierr);
225 if (n>0){
226 ierr = VecGetArrayRead(VecLow,&v1);CHKERRQ(ierr);
227 if (VecLow != VecHigh){
228 ierr = VecGetArrayRead(VecHigh,&v2);CHKERRQ(ierr);
229 } else {
230 v2=v1;
231 }
232 if (V != VecLow && V != VecHigh){
233 ierr = VecGetArrayRead(V,&vmiddle);CHKERRQ(ierr);
234 } else if (V==VecLow){
235 vmiddle=v1;
236 } else {
237 vmiddle=v2;
238 }
239
240 ierr = PetscMalloc1(n,&vm);CHKERRQ(ierr);
241
242 for (i=0; i<n; ++i){
243 if (PetscRealPart(v1[i]) < PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) < PetscRealPart(v2[i])) {vm[n_vm]=low+i; ++n_vm;}
244 }
245
246 ierr = VecRestoreArrayRead(VecLow,&v1);CHKERRQ(ierr);
247 if (VecLow != VecHigh){
248 ierr = VecRestoreArrayRead(VecHigh,&v2);CHKERRQ(ierr);
249 }
250 if (V != VecLow && V != VecHigh){
251 ierr = VecRestoreArrayRead(V,&vmiddle);CHKERRQ(ierr);
252 }
253 }
254 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);CHKERRQ(ierr);
255 PetscFunctionReturn(0);
256 }
257
258
259 /*@
260 VecWhichBetweenOrEqual - Creates an index set containing the indices
261 where VecLow <= V <= VecHigh
262
263 Collective on S
264
265 Input Parameters:
266 + VecLow - lower bound
267 . V - Vector to compare
268 - VecHigh - higher bound
269
270 OutputParameter:
271 . S - The index set containing the indices i where veclow[i] <= v[i] <= vechigh[i]
272
273 Level: advanced
274 @*/
275
VecWhichBetweenOrEqual(Vec VecLow,Vec V,Vec VecHigh,IS * S)276 PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS * S)
277 {
278 PetscErrorCode ierr;
279 PetscInt i,n_vm=0;
280 PetscInt n,low,high;
281 PetscInt *vm=NULL;
282 const PetscScalar *v1,*v2,*vmiddle;
283
284 PetscFunctionBegin;
285 PetscValidHeaderSpecific(VecLow,VEC_CLASSID,1);
286 PetscValidHeaderSpecific(V,VEC_CLASSID,2);
287 PetscValidHeaderSpecific(VecHigh,VEC_CLASSID,3);
288 PetscCheckSameComm(V,2,VecLow,1);
289 PetscCheckSameComm(V,2,VecHigh,3);
290 VecCheckSameSize(V,2,VecLow,1);
291 VecCheckSameSize(V,2,VecHigh,3);
292
293 ierr = VecGetOwnershipRange(VecLow,&low,&high);CHKERRQ(ierr);
294 ierr = VecGetLocalSize(VecLow,&n);CHKERRQ(ierr);
295 if (n>0){
296 ierr = VecGetArrayRead(VecLow,&v1);CHKERRQ(ierr);
297 if (VecLow != VecHigh){
298 ierr = VecGetArrayRead(VecHigh,&v2);CHKERRQ(ierr);
299 } else {
300 v2=v1;
301 }
302 if (V != VecLow && V != VecHigh){
303 ierr = VecGetArrayRead(V,&vmiddle);CHKERRQ(ierr);
304 } else if (V==VecLow){
305 vmiddle=v1;
306 } else {
307 vmiddle =v2;
308 }
309
310 ierr = PetscMalloc1(n,&vm);CHKERRQ(ierr);
311
312 for (i=0; i<n; ++i){
313 if (PetscRealPart(v1[i]) <= PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) <= PetscRealPart(v2[i])) {vm[n_vm]=low+i; ++n_vm;}
314 }
315
316 ierr = VecRestoreArrayRead(VecLow,&v1);CHKERRQ(ierr);
317 if (VecLow != VecHigh){
318 ierr = VecRestoreArrayRead(VecHigh,&v2);CHKERRQ(ierr);
319 }
320 if (V != VecLow && V != VecHigh){
321 ierr = VecRestoreArrayRead(V,&vmiddle);CHKERRQ(ierr);
322 }
323 }
324 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);CHKERRQ(ierr);
325 PetscFunctionReturn(0);
326 }
327
328 /*@
329 VecWhichInactive - Creates an index set containing the indices
330 where one of the following holds:
331 a) VecLow(i) < V(i) < VecHigh(i)
332 b) VecLow(i) = V(i) and D(i) <= 0 (< 0 when Strong is true)
333 c) VecHigh(i) = V(i) and D(i) >= 0 (> 0 when Strong is true)
334
335 Collective on S
336
337 Input Parameters:
338 + VecLow - lower bound
339 . V - Vector to compare
340 . D - Direction to compare
341 . VecHigh - higher bound
342 - Strong - indicator for applying strongly inactive test
343
344 OutputParameter:
345 . S - The index set containing the indices i where the bound is inactive
346
347 Level: advanced
348 @*/
349
VecWhichInactive(Vec VecLow,Vec V,Vec D,Vec VecHigh,PetscBool Strong,IS * S)350 PetscErrorCode VecWhichInactive(Vec VecLow, Vec V, Vec D, Vec VecHigh, PetscBool Strong, IS * S)
351 {
352 PetscErrorCode ierr;
353 PetscInt i,n_vm=0;
354 PetscInt n,low,high;
355 PetscInt *vm=NULL;
356 const PetscScalar *v1,*v2,*v,*d;
357
358 PetscFunctionBegin;
359 PetscValidHeaderSpecific(VecLow,VEC_CLASSID,1);
360 PetscValidHeaderSpecific(V,VEC_CLASSID,2);
361 PetscValidHeaderSpecific(D,VEC_CLASSID,3);
362 PetscValidHeaderSpecific(VecHigh,VEC_CLASSID,4);
363 PetscCheckSameComm(V,2,VecLow,1);
364 PetscCheckSameComm(V,2,D,3);
365 PetscCheckSameComm(V,2,VecHigh,4);
366 VecCheckSameSize(V,2,VecLow,1);
367 VecCheckSameSize(V,2,D,3);
368 VecCheckSameSize(V,2,VecHigh,4);
369
370 ierr = VecGetOwnershipRange(VecLow,&low,&high);CHKERRQ(ierr);
371 ierr = VecGetLocalSize(VecLow,&n);CHKERRQ(ierr);
372 if (n>0){
373 ierr = VecGetArrayRead(VecLow,&v1);CHKERRQ(ierr);
374 if (VecLow != VecHigh){
375 ierr = VecGetArrayRead(VecHigh,&v2);CHKERRQ(ierr);
376 } else {
377 v2=v1;
378 }
379 if (V != VecLow && V != VecHigh){
380 ierr = VecGetArrayRead(V,&v);CHKERRQ(ierr);
381 } else if (V==VecLow){
382 v = v1;
383 } else {
384 v = v2;
385 }
386 if (D != VecLow && D != VecHigh && D != V){
387 ierr = VecGetArrayRead(D,&d);CHKERRQ(ierr);
388 } else if (D==VecLow){
389 d = v1;
390 } else if (D==VecHigh){
391 d = v2;
392 } else {
393 d = v;
394 }
395
396 ierr = PetscMalloc1(n,&vm);CHKERRQ(ierr);
397
398 if (Strong){
399 for (i=0; i<n; ++i) {
400 if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])){
401 vm[n_vm]=low+i; ++n_vm;
402 } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) < 0){
403 vm[n_vm]=low+i; ++n_vm;
404 } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) > 0){
405 vm[n_vm]=low+i; ++n_vm;
406 }
407 }
408 } else {
409 for (i=0; i<n; ++i) {
410 if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])){
411 vm[n_vm]=low+i; ++n_vm;
412 } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) <= 0){
413 vm[n_vm]=low+i; ++n_vm;
414 } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) >= 0){
415 vm[n_vm]=low+i; ++n_vm;
416 }
417 }
418 }
419
420 ierr = VecRestoreArrayRead(VecLow,&v1);CHKERRQ(ierr);
421 if (VecLow != VecHigh){
422 ierr = VecRestoreArrayRead(VecHigh,&v2);CHKERRQ(ierr);
423 }
424 if (V != VecLow && V != VecHigh){
425 ierr = VecRestoreArrayRead(V,&v);CHKERRQ(ierr);
426 }
427 if (D != VecLow && D != VecHigh && D != V){
428 ierr = VecRestoreArrayRead(D,&d);CHKERRQ(ierr);
429 }
430 }
431 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);CHKERRQ(ierr);
432 PetscFunctionReturn(0);
433 }
434
435
436 /*@
437 VecISAXPY - Adds a reduced vector to the appropriate elements of a full-space vector.
438 vfull[is[i]] += alpha*vreduced[i]
439
440 Input Parameters:
441 + vfull - the full-space vector
442 . is - the index set for the reduced space
443 . alpha - the scalar coefficient
444 - vreduced - the reduced-space vector
445
446 Output Parameters:
447 . vfull - the sum of the full-space vector and reduced-space vector
448
449 Notes:
450 The index set identifies entries in the global vector.
451 Negative indices are skipped; indices outside the ownership range of vfull will raise an error.
452
453 Level: advanced
454
455 .seealso: VecAXPY(), VecGetOwnershipRange()
456 @*/
VecISAXPY(Vec vfull,IS is,PetscScalar alpha,Vec vreduced)457 PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha, Vec vreduced)
458 {
459 PetscInt nfull,nreduced;
460 PetscErrorCode ierr;
461
462 PetscFunctionBegin;
463 PetscValidHeaderSpecific(vfull,VEC_CLASSID,1);
464 PetscValidHeaderSpecific(is,IS_CLASSID,2);
465 PetscValidHeaderSpecific(vreduced,VEC_CLASSID,4);
466 ierr = VecGetSize(vfull,&nfull);CHKERRQ(ierr);
467 ierr = VecGetSize(vreduced,&nreduced);CHKERRQ(ierr);
468
469 if (nfull == nreduced) { /* Also takes care of masked vectors */
470 ierr = VecAXPY(vfull,alpha,vreduced);CHKERRQ(ierr);
471 } else {
472 PetscScalar *y;
473 const PetscScalar *x;
474 PetscInt i,n,m,rstart,rend;
475 const PetscInt *id;
476
477 ierr = VecGetArray(vfull,&y);CHKERRQ(ierr);
478 ierr = VecGetArrayRead(vreduced,&x);CHKERRQ(ierr);
479 ierr = ISGetIndices(is,&id);CHKERRQ(ierr);
480 ierr = ISGetLocalSize(is,&n);CHKERRQ(ierr);
481 ierr = VecGetLocalSize(vreduced,&m);CHKERRQ(ierr);
482 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IS local length not equal to Vec local length");
483 ierr = VecGetOwnershipRange(vfull,&rstart,&rend);CHKERRQ(ierr);
484 y -= rstart;
485 if (alpha == 1.0) {
486 for (i=0; i<n; ++i) {
487 if (id[i] < 0) continue;
488 if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
489 y[id[i]] += x[i];
490 }
491 } else {
492 for (i=0; i<n; ++i) {
493 if (id[i] < 0) continue;
494 if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
495 y[id[i]] += alpha*x[i];
496 }
497 }
498 y += rstart;
499 ierr = ISRestoreIndices(is,&id);CHKERRQ(ierr);
500 ierr = VecRestoreArray(vfull,&y);CHKERRQ(ierr);
501 ierr = VecRestoreArrayRead(vreduced,&x);CHKERRQ(ierr);
502 }
503 PetscFunctionReturn(0);
504 }
505
506 /*@
507 VecISCopy - Copies between a reduced vector and the appropriate elements of a full-space vector.
508
509 Input Parameters:
510 + vfull - the full-space vector
511 . is - the index set for the reduced space
512 . mode - the direction of copying, SCATTER_FORWARD or SCATTER_REVERSE
513 - vreduced - the reduced-space vector
514
515 Output Parameters:
516 . vfull - the sum of the full-space vector and reduced-space vector
517
518 Notes:
519 The index set identifies entries in the global vector.
520 Negative indices are skipped; indices outside the ownership range of vfull will raise an error.
521
522 mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
523 mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]
524
525 Level: advanced
526
527 .seealso: VecAXPY(), VecGetOwnershipRange()
528 @*/
VecISCopy(Vec vfull,IS is,ScatterMode mode,Vec vreduced)529 PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
530 {
531 PetscInt nfull, nreduced;
532 PetscErrorCode ierr;
533
534 PetscFunctionBegin;
535 PetscValidHeaderSpecific(vfull,VEC_CLASSID,1);
536 PetscValidHeaderSpecific(is,IS_CLASSID,2);
537 PetscValidHeaderSpecific(vreduced,VEC_CLASSID,3);
538 ierr = VecGetSize(vfull, &nfull);CHKERRQ(ierr);
539 ierr = VecGetSize(vreduced, &nreduced);CHKERRQ(ierr);
540
541 if (nfull == nreduced) { /* Also takes care of masked vectors */
542 if (mode == SCATTER_FORWARD) {
543 ierr = VecCopy(vreduced, vfull);CHKERRQ(ierr);
544 } else {
545 ierr = VecCopy(vfull, vreduced);CHKERRQ(ierr);
546 }
547 } else {
548 const PetscInt *id;
549 PetscInt i, n, m, rstart, rend;
550
551 ierr = ISGetIndices(is, &id);CHKERRQ(ierr);
552 ierr = ISGetLocalSize(is, &n);CHKERRQ(ierr);
553 ierr = VecGetLocalSize(vreduced, &m);CHKERRQ(ierr);
554 ierr = VecGetOwnershipRange(vfull, &rstart, &rend);CHKERRQ(ierr);
555 if (m != n) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "IS local length %D not equal to Vec local length %D", n, m);
556 if (mode == SCATTER_FORWARD) {
557 PetscScalar *y;
558 const PetscScalar *x;
559
560 ierr = VecGetArray(vfull, &y);CHKERRQ(ierr);
561 ierr = VecGetArrayRead(vreduced, &x);CHKERRQ(ierr);
562 y -= rstart;
563 for (i = 0; i < n; ++i) {
564 if (id[i] < 0) continue;
565 if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
566 y[id[i]] = x[i];
567 }
568 y += rstart;
569 ierr = VecRestoreArrayRead(vreduced, &x);CHKERRQ(ierr);
570 ierr = VecRestoreArray(vfull, &y);CHKERRQ(ierr);
571 } else if (mode == SCATTER_REVERSE) {
572 PetscScalar *x;
573 const PetscScalar *y;
574
575 ierr = VecGetArrayRead(vfull, &y);CHKERRQ(ierr);
576 ierr = VecGetArray(vreduced, &x);CHKERRQ(ierr);
577 for (i = 0; i < n; ++i) {
578 if (id[i] < 0) continue;
579 if (id[i] < rstart || id[i] >= rend) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
580 x[i] = y[id[i]-rstart];
581 }
582 ierr = VecRestoreArray(vreduced, &x);CHKERRQ(ierr);
583 ierr = VecRestoreArrayRead(vfull, &y);CHKERRQ(ierr);
584 } else SETERRQ(PetscObjectComm((PetscObject) vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
585 ierr = ISRestoreIndices(is, &id);CHKERRQ(ierr);
586 }
587 PetscFunctionReturn(0);
588 }
589
590 /*@
591 ISComplementVec - Creates the complement of the index set relative to a layout defined by a Vec
592
593 Collective on IS
594
595 Input Parameter:
596 + S - a PETSc IS
597 - V - the reference vector space
598
599 Output Parameter:
600 . T - the complement of S
601
602 Level: advanced
603
604 .seealso: ISCreateGeneral()
605 @*/
ISComplementVec(IS S,Vec V,IS * T)606 PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
607 {
608 PetscErrorCode ierr;
609 PetscInt start, end;
610
611 PetscFunctionBegin;
612 ierr = VecGetOwnershipRange(V,&start,&end);CHKERRQ(ierr);
613 ierr = ISComplement(S,start,end,T);CHKERRQ(ierr);
614 PetscFunctionReturn(0);
615 }
616
617 /*@
618 VecISSet - Sets the elements of a vector, specified by an index set, to a constant
619
620 Input Parameter:
621 + V - the vector
622 . S - index set for the locations in the vector
623 - c - the constant
624
625 Notes:
626 The index set identifies entries in the global vector.
627 Negative indices are skipped; indices outside the ownership range of V will raise an error.
628
629 Level: advanced
630
631 .seealso: VecSet(), VecGetOwnershipRange()
632 @*/
VecISSet(Vec V,IS S,PetscScalar c)633 PetscErrorCode VecISSet(Vec V,IS S, PetscScalar c)
634 {
635 PetscErrorCode ierr;
636 PetscInt nloc,low,high,i;
637 const PetscInt *s;
638 PetscScalar *v;
639
640 PetscFunctionBegin;
641 if (!S) PetscFunctionReturn(0); /* simply return with no-op if the index set is NULL */
642 PetscValidHeaderSpecific(V,VEC_CLASSID,1);
643 PetscValidHeaderSpecific(S,IS_CLASSID,2);
644 PetscValidType(V,3);
645
646 ierr = VecGetOwnershipRange(V,&low,&high);CHKERRQ(ierr);
647 ierr = ISGetLocalSize(S,&nloc);CHKERRQ(ierr);
648 ierr = ISGetIndices(S,&s);CHKERRQ(ierr);
649 ierr = VecGetArray(V,&v);CHKERRQ(ierr);
650 for (i=0; i<nloc; ++i){
651 if (s[i] < 0) continue;
652 if (s[i] < low || s[i] >= high) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only owned values supported");
653 v[s[i]-low] = c;
654 }
655 ierr = ISRestoreIndices(S,&s);CHKERRQ(ierr);
656 ierr = VecRestoreArray(V,&v);CHKERRQ(ierr);
657 PetscFunctionReturn(0);
658 }
659
660 #if !defined(PETSC_USE_COMPLEX)
661 /*@C
662 VecBoundGradientProjection - Projects vector according to this definition.
663 If XL[i] < X[i] < XU[i], then GP[i] = G[i];
664 If X[i] <= XL[i], then GP[i] = min(G[i],0);
665 If X[i] >= XU[i], then GP[i] = max(G[i],0);
666
667 Input Parameters:
668 + G - current gradient vector
669 . X - current solution vector with XL[i] <= X[i] <= XU[i]
670 . XL - lower bounds
671 - XU - upper bounds
672
673 Output Parameter:
674 . GP - gradient projection vector
675
676 Notes:
677 GP may be the same vector as G
678
679 Level: advanced
680 @*/
VecBoundGradientProjection(Vec G,Vec X,Vec XL,Vec XU,Vec GP)681 PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
682 {
683
684 PetscErrorCode ierr;
685 PetscInt n,i;
686 const PetscReal *xptr,*xlptr,*xuptr;
687 PetscReal *gptr,*gpptr;
688 PetscReal xval,gpval;
689
690 /* Project variables at the lower and upper bound */
691 PetscFunctionBegin;
692 PetscValidHeaderSpecific(G,VEC_CLASSID,1);
693 PetscValidHeaderSpecific(X,VEC_CLASSID,2);
694 PetscValidHeaderSpecific(XL,VEC_CLASSID,3);
695 PetscValidHeaderSpecific(XU,VEC_CLASSID,4);
696 PetscValidHeaderSpecific(GP,VEC_CLASSID,5);
697
698 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr);
699
700 ierr = VecGetArrayRead(X,&xptr);CHKERRQ(ierr);
701 ierr = VecGetArrayRead(XL,&xlptr);CHKERRQ(ierr);
702 ierr = VecGetArrayRead(XU,&xuptr);CHKERRQ(ierr);
703 ierr = VecGetArrayPair(G,GP,&gptr,&gpptr);CHKERRQ(ierr);
704
705 for (i=0; i<n; ++i){
706 gpval = gptr[i]; xval = xptr[i];
707 if (gpval>0.0 && xval<=xlptr[i]){
708 gpval = 0.0;
709 } else if (gpval<0.0 && xval>=xuptr[i]){
710 gpval = 0.0;
711 }
712 gpptr[i] = gpval;
713 }
714
715 ierr = VecRestoreArrayRead(X,&xptr);CHKERRQ(ierr);
716 ierr = VecRestoreArrayRead(XL,&xlptr);CHKERRQ(ierr);
717 ierr = VecRestoreArrayRead(XU,&xuptr);CHKERRQ(ierr);
718 ierr = VecRestoreArrayPair(G,GP,&gptr,&gpptr);CHKERRQ(ierr);
719 PetscFunctionReturn(0);
720 }
721 #endif
722
723 /*@
724 VecStepMaxBounded - See below
725
726 Collective on Vec
727
728 Input Parameters:
729 + X - vector with no negative entries
730 . XL - lower bounds
731 . XU - upper bounds
732 - DX - step direction, can have negative, positive or zero entries
733
734 Output Parameter:
735 . stepmax - minimum value so that X[i] + stepmax*DX[i] <= XL[i] or XU[i] <= X[i] + stepmax*DX[i]
736
737 Level: intermediate
738
739 @*/
VecStepMaxBounded(Vec X,Vec DX,Vec XL,Vec XU,PetscReal * stepmax)740 PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
741 {
742 PetscErrorCode ierr;
743 PetscInt i,nn;
744 const PetscScalar *xx,*dx,*xl,*xu;
745 PetscReal localmax=0;
746
747 PetscFunctionBegin;
748 PetscValidHeaderSpecific(X,VEC_CLASSID,2);
749 PetscValidHeaderSpecific(DX,VEC_CLASSID,5);
750 PetscValidHeaderSpecific(XL,VEC_CLASSID,3);
751 PetscValidHeaderSpecific(XU,VEC_CLASSID,4);
752
753 ierr = VecGetArrayRead(X,&xx);CHKERRQ(ierr);
754 ierr = VecGetArrayRead(XL,&xl);CHKERRQ(ierr);
755 ierr = VecGetArrayRead(XU,&xu);CHKERRQ(ierr);
756 ierr = VecGetArrayRead(DX,&dx);CHKERRQ(ierr);
757 ierr = VecGetLocalSize(X,&nn);CHKERRQ(ierr);
758 for (i=0;i<nn;i++){
759 if (PetscRealPart(dx[i]) > 0){
760 localmax=PetscMax(localmax,PetscRealPart((xu[i]-xx[i])/dx[i]));
761 } else if (PetscRealPart(dx[i])<0){
762 localmax=PetscMax(localmax,PetscRealPart((xl[i]-xx[i])/dx[i]));
763 }
764 }
765 ierr = VecRestoreArrayRead(X,&xx);CHKERRQ(ierr);
766 ierr = VecRestoreArrayRead(XL,&xl);CHKERRQ(ierr);
767 ierr = VecRestoreArrayRead(XU,&xu);CHKERRQ(ierr);
768 ierr = VecRestoreArrayRead(DX,&dx);CHKERRQ(ierr);
769 ierr = MPIU_Allreduce(&localmax,stepmax,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)X));CHKERRQ(ierr);
770 PetscFunctionReturn(0);
771 }
772
773 /*@
774 VecStepBoundInfo - See below
775
776 Collective on Vec
777
778 Input Parameters:
779 + X - vector with no negative entries
780 . XL - lower bounds
781 . XU - upper bounds
782 - DX - step direction, can have negative, positive or zero entries
783
784 Output Parameter:
785 + boundmin - (may be NULL this it is not computed) maximum value so that XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
786 . wolfemin - (may be NULL this it is not computed)
787 - boundmax - (may be NULL this it is not computed) minimum value so that X[i] + boundmax*DX[i] <= XL[i] or XU[i] <= X[i] + boundmax*DX[i]
788
789 Notes:
790 For complex numbers only compares the real part
791
792 Level: advanced
793 @*/
VecStepBoundInfo(Vec X,Vec DX,Vec XL,Vec XU,PetscReal * boundmin,PetscReal * wolfemin,PetscReal * boundmax)794 PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
795 {
796 PetscErrorCode ierr;
797 PetscInt n,i;
798 const PetscScalar *x,*xl,*xu,*dx;
799 PetscReal t;
800 PetscReal localmin=PETSC_INFINITY,localwolfemin=PETSC_INFINITY,localmax=-1;
801 MPI_Comm comm;
802
803 PetscFunctionBegin;
804 PetscValidHeaderSpecific(X,VEC_CLASSID,1);
805 PetscValidHeaderSpecific(XL,VEC_CLASSID,2);
806 PetscValidHeaderSpecific(XU,VEC_CLASSID,3);
807 PetscValidHeaderSpecific(DX,VEC_CLASSID,4);
808
809 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
810 ierr = VecGetArrayRead(XL,&xl);CHKERRQ(ierr);
811 ierr = VecGetArrayRead(XU,&xu);CHKERRQ(ierr);
812 ierr = VecGetArrayRead(DX,&dx);CHKERRQ(ierr);
813 ierr = VecGetLocalSize(X,&n);CHKERRQ(ierr);
814 for (i=0; i<n; ++i){
815 if (PetscRealPart(dx[i])>0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
816 t=PetscRealPart((xu[i]-x[i])/dx[i]);
817 localmin=PetscMin(t,localmin);
818 if (localmin>0){
819 localwolfemin = PetscMin(t,localwolfemin);
820 }
821 localmax = PetscMax(t,localmax);
822 } else if (PetscRealPart(dx[i])<0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
823 t=PetscRealPart((xl[i]-x[i])/dx[i]);
824 localmin = PetscMin(t,localmin);
825 if (localmin>0){
826 localwolfemin = PetscMin(t,localwolfemin);
827 }
828 localmax = PetscMax(t,localmax);
829 }
830 }
831
832 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
833 ierr = VecRestoreArrayRead(XL,&xl);CHKERRQ(ierr);
834 ierr = VecRestoreArrayRead(XU,&xu);CHKERRQ(ierr);
835 ierr = VecRestoreArrayRead(DX,&dx);CHKERRQ(ierr);
836 ierr = PetscObjectGetComm((PetscObject)X,&comm);CHKERRQ(ierr);
837
838 if (boundmin){
839 ierr = MPIU_Allreduce(&localmin,boundmin,1,MPIU_REAL,MPIU_MIN,comm);CHKERRQ(ierr);
840 ierr = PetscInfo1(X,"Step Bound Info: Closest Bound: %20.19e\n",(double)*boundmin);CHKERRQ(ierr);
841 }
842 if (wolfemin){
843 ierr = MPIU_Allreduce(&localwolfemin,wolfemin,1,MPIU_REAL,MPIU_MIN,comm);CHKERRQ(ierr);
844 ierr = PetscInfo1(X,"Step Bound Info: Wolfe: %20.19e\n",(double)*wolfemin);CHKERRQ(ierr);
845 }
846 if (boundmax) {
847 ierr = MPIU_Allreduce(&localmax,boundmax,1,MPIU_REAL,MPIU_MAX,comm);CHKERRQ(ierr);
848 if (*boundmax < 0) *boundmax=PETSC_INFINITY;
849 ierr = PetscInfo1(X,"Step Bound Info: Max: %20.19e\n",(double)*boundmax);CHKERRQ(ierr);
850 }
851 PetscFunctionReturn(0);
852 }
853
854 /*@
855 VecStepMax - Returns the largest value so that x[i] + step*DX[i] >= 0 for all i
856
857 Collective on Vec
858
859 Input Parameters:
860 + X - vector with no negative entries
861 - DX - a step direction, can have negative, positive or zero entries
862
863 Output Parameter:
864 . step - largest value such that x[i] + step*DX[i] >= 0 for all i
865
866 Notes:
867 For complex numbers only compares the real part
868
869 Level: advanced
870 @*/
VecStepMax(Vec X,Vec DX,PetscReal * step)871 PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
872 {
873 PetscErrorCode ierr;
874 PetscInt i,nn;
875 PetscReal stepmax=PETSC_INFINITY;
876 const PetscScalar *xx,*dx;
877
878 PetscFunctionBegin;
879 PetscValidHeaderSpecific(X,VEC_CLASSID,1);
880 PetscValidHeaderSpecific(DX,VEC_CLASSID,2);
881
882 ierr = VecGetLocalSize(X,&nn);CHKERRQ(ierr);
883 ierr = VecGetArrayRead(X,&xx);CHKERRQ(ierr);
884 ierr = VecGetArrayRead(DX,&dx);CHKERRQ(ierr);
885 for (i=0;i<nn;++i){
886 if (PetscRealPart(xx[i]) < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Vector must be positive");
887 else if (PetscRealPart(dx[i])<0) stepmax=PetscMin(stepmax,PetscRealPart(-xx[i]/dx[i]));
888 }
889 ierr = VecRestoreArrayRead(X,&xx);CHKERRQ(ierr);
890 ierr = VecRestoreArrayRead(DX,&dx);CHKERRQ(ierr);
891 ierr = MPIU_Allreduce(&stepmax,step,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)X));CHKERRQ(ierr);
892 PetscFunctionReturn(0);
893 }
894
895 /*@
896 VecPow - Replaces each component of a vector by x_i^p
897
898 Logically Collective on v
899
900 Input Parameter:
901 + v - the vector
902 - p - the exponent to use on each element
903
904 Output Parameter:
905 . v - the vector
906
907 Level: intermediate
908
909 @*/
VecPow(Vec v,PetscScalar p)910 PetscErrorCode VecPow(Vec v, PetscScalar p)
911 {
912 PetscErrorCode ierr;
913 PetscInt n,i;
914 PetscScalar *v1;
915
916 PetscFunctionBegin;
917 PetscValidHeaderSpecific(v,VEC_CLASSID,1);
918
919 ierr = VecGetArray(v,&v1);CHKERRQ(ierr);
920 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
921
922 if (1.0 == p) {
923 } else if (-1.0 == p) {
924 for (i = 0; i < n; ++i){
925 v1[i] = 1.0 / v1[i];
926 }
927 } else if (0.0 == p) {
928 for (i = 0; i < n; ++i){
929 /* Not-a-number left alone
930 Infinity set to one */
931 if (v1[i] == v1[i]) {
932 v1[i] = 1.0;
933 }
934 }
935 } else if (0.5 == p) {
936 for (i = 0; i < n; ++i) {
937 if (PetscRealPart(v1[i]) >= 0) {
938 v1[i] = PetscSqrtScalar(v1[i]);
939 } else {
940 v1[i] = PETSC_INFINITY;
941 }
942 }
943 } else if (-0.5 == p) {
944 for (i = 0; i < n; ++i) {
945 if (PetscRealPart(v1[i]) >= 0) {
946 v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
947 } else {
948 v1[i] = PETSC_INFINITY;
949 }
950 }
951 } else if (2.0 == p) {
952 for (i = 0; i < n; ++i){
953 v1[i] *= v1[i];
954 }
955 } else if (-2.0 == p) {
956 for (i = 0; i < n; ++i){
957 v1[i] = 1.0 / (v1[i] * v1[i]);
958 }
959 } else {
960 for (i = 0; i < n; ++i) {
961 if (PetscRealPart(v1[i]) >= 0) {
962 v1[i] = PetscPowScalar(v1[i],p);
963 } else {
964 v1[i] = PETSC_INFINITY;
965 }
966 }
967 }
968 ierr = VecRestoreArray(v,&v1);CHKERRQ(ierr);
969 PetscFunctionReturn(0);
970 }
971
972 /*@
973 VecMedian - Computes the componentwise median of three vectors
974 and stores the result in this vector. Used primarily for projecting
975 a vector within upper and lower bounds.
976
977 Logically Collective
978
979 Input Parameters:
980 . Vec1, Vec2, Vec3 - The three vectors
981
982 Output Parameter:
983 . VMedian - The median vector (this can be any one of the input vectors)
984
985 Level: advanced
986 @*/
VecMedian(Vec Vec1,Vec Vec2,Vec Vec3,Vec VMedian)987 PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
988 {
989 PetscErrorCode ierr;
990 PetscInt i,n,low1,high1;
991 const PetscScalar *v1,*v2,*v3;
992 PetscScalar *vmed;
993
994 PetscFunctionBegin;
995 PetscValidHeaderSpecific(Vec1,VEC_CLASSID,1);
996 PetscValidHeaderSpecific(Vec2,VEC_CLASSID,2);
997 PetscValidHeaderSpecific(Vec3,VEC_CLASSID,3);
998 PetscValidHeaderSpecific(VMedian,VEC_CLASSID,4);
999
1000 if (Vec1==Vec2 || Vec1==Vec3){
1001 ierr = VecCopy(Vec1,VMedian);CHKERRQ(ierr);
1002 PetscFunctionReturn(0);
1003 }
1004 if (Vec2==Vec3){
1005 ierr = VecCopy(Vec2,VMedian);CHKERRQ(ierr);
1006 PetscFunctionReturn(0);
1007 }
1008
1009 /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
1010 PetscValidType(Vec1,1);
1011 PetscValidType(Vec2,2);
1012 PetscValidType(Vec3,3);
1013 PetscValidType(VMedian,4);
1014 PetscCheckSameType(Vec1,1,Vec2,2);
1015 PetscCheckSameType(Vec1,1,Vec3,3);
1016 PetscCheckSameType(Vec1,1,VMedian,4);
1017 PetscCheckSameComm(Vec1,1,Vec2,2);
1018 PetscCheckSameComm(Vec1,1,Vec3,3);
1019 PetscCheckSameComm(Vec1,1,VMedian,4);
1020 VecCheckSameSize(Vec1,1,Vec2,2);
1021 VecCheckSameSize(Vec1,1,Vec3,3);
1022 VecCheckSameSize(Vec1,1,VMedian,4);
1023
1024 ierr = VecGetOwnershipRange(Vec1,&low1,&high1);CHKERRQ(ierr);
1025 ierr = VecGetLocalSize(Vec1,&n);CHKERRQ(ierr);
1026 if (n>0){
1027 ierr = VecGetArray(VMedian,&vmed);CHKERRQ(ierr);
1028 if (Vec1 != VMedian){
1029 ierr = VecGetArrayRead(Vec1,&v1);CHKERRQ(ierr);
1030 } else {
1031 v1=vmed;
1032 }
1033 if (Vec2 != VMedian){
1034 ierr = VecGetArrayRead(Vec2,&v2);CHKERRQ(ierr);
1035 } else {
1036 v2=vmed;
1037 }
1038 if (Vec3 != VMedian){
1039 ierr = VecGetArrayRead(Vec3,&v3);CHKERRQ(ierr);
1040 } else {
1041 v3=vmed;
1042 }
1043
1044 for (i=0;i<n;++i){
1045 vmed[i]=PetscMax(PetscMax(PetscMin(PetscRealPart(v1[i]),PetscRealPart(v2[i])),PetscMin(PetscRealPart(v1[i]),PetscRealPart(v3[i]))),PetscMin(PetscRealPart(v2[i]),PetscRealPart(v3[i])));
1046 }
1047
1048 ierr = VecRestoreArray(VMedian,&vmed);CHKERRQ(ierr);
1049 if (VMedian != Vec1){
1050 ierr = VecRestoreArrayRead(Vec1,&v1);CHKERRQ(ierr);
1051 }
1052 if (VMedian != Vec2){
1053 ierr = VecRestoreArrayRead(Vec2,&v2);CHKERRQ(ierr);
1054 }
1055 if (VMedian != Vec3){
1056 ierr = VecRestoreArrayRead(Vec3,&v3);CHKERRQ(ierr);
1057 }
1058 }
1059 PetscFunctionReturn(0);
1060 }
1061