Actual source code: projection.c

  1: #include <petsc/private/vecimpl.h>

  3: /*@
  4:   VecWhichEqual - Creates an index set containing the indices
  5:              where the vectors Vec1 and Vec2 have identical elements.

  7:   Collective on Vec

  9:   Input Parameters:
 10: . Vec1, Vec2 - the two vectors to compare

 12:   OutputParameter:
 13: . S - The index set containing the indices i where vec1[i] == vec2[i]

 15:   Notes:
 16:     the two vectors must have the same parallel layout

 18:   Level: advanced
 19: @*/
 20: PetscErrorCode VecWhichEqual(Vec Vec1, Vec Vec2, IS *S)
 21: {
 22:   PetscInt          i,n_same=0;
 23:   PetscInt          n,low,high;
 24:   PetscInt          *same=NULL;
 25:   const PetscScalar *v1,*v2;

 30:   VecCheckSameSize(Vec1,1,Vec2,2);

 32:   VecGetOwnershipRange(Vec1,&low,&high);
 33:   VecGetLocalSize(Vec1,&n);
 34:   if (n>0) {
 35:     if (Vec1 == Vec2) {
 36:       VecGetArrayRead(Vec1,&v1);
 37:       v2=v1;
 38:     } else {
 39:       VecGetArrayRead(Vec1,&v1);
 40:       VecGetArrayRead(Vec2,&v2);
 41:     }

 43:     PetscMalloc1(n,&same);

 45:     for (i=0; i<n; ++i) {
 46:       if (v1[i] == v2[i]) {same[n_same]=low+i; ++n_same;}
 47:     }

 49:     if (Vec1 == Vec2) {
 50:       VecRestoreArrayRead(Vec1,&v1);
 51:     } else {
 52:       VecRestoreArrayRead(Vec1,&v1);
 53:       VecRestoreArrayRead(Vec2,&v2);
 54:     }
 55:   }
 56:   ISCreateGeneral(PetscObjectComm((PetscObject)Vec1),n_same,same,PETSC_OWN_POINTER,S);
 57:   return 0;
 58: }

 60: /*@
 61:   VecWhichLessThan - Creates an index set containing the indices
 62:   where the vectors Vec1 < Vec2

 64:   Collective on S

 66:   Input Parameters:
 67: . Vec1, Vec2 - the two vectors to compare

 69:   OutputParameter:
 70: . S - The index set containing the indices i where vec1[i] < vec2[i]

 72:   Notes:
 73:   The two vectors must have the same parallel layout

 75:   For complex numbers this only compares the real part

 77:   Level: advanced
 78: @*/
 79: PetscErrorCode VecWhichLessThan(Vec Vec1, Vec Vec2, IS *S)
 80: {
 81:   PetscInt          i,n_lt=0;
 82:   PetscInt          n,low,high;
 83:   PetscInt          *lt=NULL;
 84:   const PetscScalar *v1,*v2;

 89:   VecCheckSameSize(Vec1,1,Vec2,2);

 91:   VecGetOwnershipRange(Vec1,&low,&high);
 92:   VecGetLocalSize(Vec1,&n);
 93:   if (n>0) {
 94:     if (Vec1 == Vec2) {
 95:       VecGetArrayRead(Vec1,&v1);
 96:       v2=v1;
 97:     } else {
 98:       VecGetArrayRead(Vec1,&v1);
 99:       VecGetArrayRead(Vec2,&v2);
100:     }

102:     PetscMalloc1(n,&lt);

104:     for (i=0; i<n; ++i) {
105:       if (PetscRealPart(v1[i]) < PetscRealPart(v2[i])) {lt[n_lt]=low+i; ++n_lt;}
106:     }

108:     if (Vec1 == Vec2) {
109:       VecRestoreArrayRead(Vec1,&v1);
110:     } else {
111:       VecRestoreArrayRead(Vec1,&v1);
112:       VecRestoreArrayRead(Vec2,&v2);
113:     }
114:   }
115:   ISCreateGeneral(PetscObjectComm((PetscObject)Vec1),n_lt,lt,PETSC_OWN_POINTER,S);
116:   return 0;
117: }

119: /*@
120:   VecWhichGreaterThan - Creates an index set containing the indices
121:   where the vectors Vec1 > Vec2

123:   Collective on S

125:   Input Parameters:
126: . Vec1, Vec2 - the two vectors to compare

128:   OutputParameter:
129: . S - The index set containing the indices i where vec1[i] > vec2[i]

131:   Notes:
132:   The two vectors must have the same parallel layout

134:   For complex numbers this only compares the real part

136:   Level: advanced
137: @*/
138: PetscErrorCode VecWhichGreaterThan(Vec Vec1, Vec Vec2, IS *S)
139: {
140:   PetscInt          i,n_gt=0;
141:   PetscInt          n,low,high;
142:   PetscInt          *gt=NULL;
143:   const PetscScalar *v1,*v2;

148:   VecCheckSameSize(Vec1,1,Vec2,2);

150:   VecGetOwnershipRange(Vec1,&low,&high);
151:   VecGetLocalSize(Vec1,&n);
152:   if (n>0) {
153:     if (Vec1 == Vec2) {
154:       VecGetArrayRead(Vec1,&v1);
155:       v2=v1;
156:     } else {
157:       VecGetArrayRead(Vec1,&v1);
158:       VecGetArrayRead(Vec2,&v2);
159:     }

161:     PetscMalloc1(n,&gt);

163:     for (i=0; i<n; ++i) {
164:       if (PetscRealPart(v1[i]) > PetscRealPart(v2[i])) {gt[n_gt]=low+i; ++n_gt;}
165:     }

167:     if (Vec1 == Vec2) {
168:       VecRestoreArrayRead(Vec1,&v1);
169:     } else {
170:       VecRestoreArrayRead(Vec1,&v1);
171:       VecRestoreArrayRead(Vec2,&v2);
172:     }
173:   }
174:   ISCreateGeneral(PetscObjectComm((PetscObject)Vec1),n_gt,gt,PETSC_OWN_POINTER,S);
175:   return 0;
176: }

178: /*@
179:   VecWhichBetween - Creates an index set containing the indices
180:                where  VecLow < V < VecHigh

182:   Collective on S

184:   Input Parameters:
185: + VecLow - lower bound
186: . V - Vector to compare
187: - VecHigh - higher bound

189:   OutputParameter:
190: . S - The index set containing the indices i where veclow[i] < v[i] < vechigh[i]

192:   Notes:
193:   The vectors must have the same parallel layout

195:   For complex numbers this only compares the real part

197:   Level: advanced
198: @*/
199: PetscErrorCode VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS *S)
200: {

202:   PetscInt          i,n_vm=0;
203:   PetscInt          n,low,high;
204:   PetscInt          *vm=NULL;
205:   const PetscScalar *v1,*v2,*vmiddle;

212:   VecCheckSameSize(V,2,VecLow,1);
213:   VecCheckSameSize(V,2,VecHigh,3);

215:   VecGetOwnershipRange(VecLow,&low,&high);
216:   VecGetLocalSize(VecLow,&n);
217:   if (n>0) {
218:     VecGetArrayRead(VecLow,&v1);
219:     if (VecLow != VecHigh) {
220:       VecGetArrayRead(VecHigh,&v2);
221:     } else {
222:       v2=v1;
223:     }
224:     if (V != VecLow && V != VecHigh) {
225:       VecGetArrayRead(V,&vmiddle);
226:     } else if (V==VecLow) {
227:       vmiddle=v1;
228:     } else {
229:       vmiddle=v2;
230:     }

232:     PetscMalloc1(n,&vm);

234:     for (i=0; i<n; ++i) {
235:       if (PetscRealPart(v1[i]) < PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) < PetscRealPart(v2[i])) {vm[n_vm]=low+i; ++n_vm;}
236:     }

238:     VecRestoreArrayRead(VecLow,&v1);
239:     if (VecLow != VecHigh) {
240:       VecRestoreArrayRead(VecHigh,&v2);
241:     }
242:     if (V != VecLow && V != VecHigh) {
243:       VecRestoreArrayRead(V,&vmiddle);
244:     }
245:   }
246:   ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
247:   return 0;
248: }

250: /*@
251:   VecWhichBetweenOrEqual - Creates an index set containing the indices
252:   where  VecLow <= V <= VecHigh

254:   Collective on S

256:   Input Parameters:
257: + VecLow - lower bound
258: . V - Vector to compare
259: - VecHigh - higher bound

261:   OutputParameter:
262: . S - The index set containing the indices i where veclow[i] <= v[i] <= vechigh[i]

264:   Level: advanced
265: @*/

267: PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS * S)
268: {
269:   PetscInt          i,n_vm=0;
270:   PetscInt          n,low,high;
271:   PetscInt          *vm=NULL;
272:   const PetscScalar *v1,*v2,*vmiddle;

279:   VecCheckSameSize(V,2,VecLow,1);
280:   VecCheckSameSize(V,2,VecHigh,3);

282:   VecGetOwnershipRange(VecLow,&low,&high);
283:   VecGetLocalSize(VecLow,&n);
284:   if (n>0) {
285:     VecGetArrayRead(VecLow,&v1);
286:     if (VecLow != VecHigh) {
287:       VecGetArrayRead(VecHigh,&v2);
288:     } else {
289:       v2=v1;
290:     }
291:     if (V != VecLow && V != VecHigh) {
292:       VecGetArrayRead(V,&vmiddle);
293:     } else if (V==VecLow) {
294:       vmiddle=v1;
295:     } else {
296:       vmiddle =v2;
297:     }

299:     PetscMalloc1(n,&vm);

301:     for (i=0; i<n; ++i) {
302:       if (PetscRealPart(v1[i]) <= PetscRealPart(vmiddle[i]) && PetscRealPart(vmiddle[i]) <= PetscRealPart(v2[i])) {vm[n_vm]=low+i; ++n_vm;}
303:     }

305:     VecRestoreArrayRead(VecLow,&v1);
306:     if (VecLow != VecHigh) {
307:       VecRestoreArrayRead(VecHigh,&v2);
308:     }
309:     if (V != VecLow && V != VecHigh) {
310:       VecRestoreArrayRead(V,&vmiddle);
311:     }
312:   }
313:   ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
314:   return 0;
315: }

317: /*@
318:    VecWhichInactive - Creates an index set containing the indices
319:   where one of the following holds:
320:     a) VecLow(i)  < V(i) < VecHigh(i)
321:     b) VecLow(i)  = V(i) and D(i) <= 0 (< 0 when Strong is true)
322:     c) VecHigh(i) = V(i) and D(i) >= 0 (> 0 when Strong is true)

324:   Collective on S

326:   Input Parameters:
327: + VecLow - lower bound
328: . V - Vector to compare
329: . D - Direction to compare
330: . VecHigh - higher bound
331: - Strong - indicator for applying strongly inactive test

333:   OutputParameter:
334: . S - The index set containing the indices i where the bound is inactive

336:   Level: advanced
337: @*/

339: PetscErrorCode VecWhichInactive(Vec VecLow, Vec V, Vec D, Vec VecHigh, PetscBool Strong, IS * S)
340: {
341:   PetscInt          i,n_vm=0;
342:   PetscInt          n,low,high;
343:   PetscInt          *vm=NULL;
344:   const PetscScalar *v1,*v2,*v,*d;

353:   VecCheckSameSize(V,2,VecLow,1);
354:   VecCheckSameSize(V,2,D,3);
355:   VecCheckSameSize(V,2,VecHigh,4);

357:   VecGetOwnershipRange(VecLow,&low,&high);
358:   VecGetLocalSize(VecLow,&n);
359:   if (n>0) {
360:     VecGetArrayRead(VecLow,&v1);
361:     if (VecLow != VecHigh) {
362:       VecGetArrayRead(VecHigh,&v2);
363:     } else {
364:       v2=v1;
365:     }
366:     if (V != VecLow && V != VecHigh) {
367:       VecGetArrayRead(V,&v);
368:     } else if (V==VecLow) {
369:       v = v1;
370:     } else {
371:       v = v2;
372:     }
373:     if (D != VecLow && D != VecHigh && D != V) {
374:       VecGetArrayRead(D,&d);
375:     } else if (D==VecLow) {
376:       d = v1;
377:     } else if (D==VecHigh) {
378:       d = v2;
379:     } else {
380:       d = v;
381:     }

383:     PetscMalloc1(n,&vm);

385:     if (Strong) {
386:       for (i=0; i<n; ++i) {
387:         if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
388:           vm[n_vm]=low+i; ++n_vm;
389:         } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) < 0) {
390:           vm[n_vm]=low+i; ++n_vm;
391:         } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) > 0) {
392:           vm[n_vm]=low+i; ++n_vm;
393:         }
394:       }
395:     } else {
396:       for (i=0; i<n; ++i) {
397:         if (PetscRealPart(v1[i]) < PetscRealPart(v[i]) && PetscRealPart(v[i]) < PetscRealPart(v2[i])) {
398:           vm[n_vm]=low+i; ++n_vm;
399:         } else if (PetscRealPart(v1[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) <= 0) {
400:           vm[n_vm]=low+i; ++n_vm;
401:         } else if (PetscRealPart(v2[i]) == PetscRealPart(v[i]) && PetscRealPart(d[i]) >= 0) {
402:           vm[n_vm]=low+i; ++n_vm;
403:         }
404:       }
405:     }

407:     VecRestoreArrayRead(VecLow,&v1);
408:     if (VecLow != VecHigh) {
409:       VecRestoreArrayRead(VecHigh,&v2);
410:     }
411:     if (V != VecLow && V != VecHigh) {
412:       VecRestoreArrayRead(V,&v);
413:     }
414:     if (D != VecLow && D != VecHigh && D != V) {
415:       VecRestoreArrayRead(D,&d);
416:     }
417:   }
418:   ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
419:   return 0;
420: }

422: /*@
423:   VecISAXPY - Adds a reduced vector to the appropriate elements of a full-space vector.
424:                   vfull[is[i]] += alpha*vreduced[i]

426:   Input Parameters:
427: + vfull    - the full-space vector
428: . is       - the index set for the reduced space
429: . alpha    - the scalar coefficient
430: - vreduced - the reduced-space vector

432:   Output Parameters:
433: . vfull    - the sum of the full-space vector and reduced-space vector

435:   Notes:
436:     The index set identifies entries in the global vector.
437:     Negative indices are skipped; indices outside the ownership range of vfull will raise an error.

439:   Level: advanced

441: .seealso:  VecAXPY(), VecGetOwnershipRange()
442: @*/
443: PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha, Vec vreduced)
444: {
445:   PetscInt       nfull,nreduced;

450:   VecGetSize(vfull,&nfull);
451:   VecGetSize(vreduced,&nreduced);

453:   if (nfull == nreduced) { /* Also takes care of masked vectors */
454:     VecAXPY(vfull,alpha,vreduced);
455:   } else {
456:     PetscScalar      *y;
457:     const PetscScalar *x;
458:     PetscInt          i,n,m,rstart,rend;
459:     const PetscInt    *id;

461:     VecGetArray(vfull,&y);
462:     VecGetArrayRead(vreduced,&x);
463:     ISGetIndices(is,&id);
464:     ISGetLocalSize(is,&n);
465:     VecGetLocalSize(vreduced,&m);
467:     VecGetOwnershipRange(vfull,&rstart,&rend);
468:     y   -= rstart;
469:     if (alpha == 1.0) {
470:       for (i=0; i<n; ++i) {
471:         if (id[i] < 0) continue;
473:         y[id[i]] += x[i];
474:       }
475:     } else {
476:       for (i=0; i<n; ++i) {
477:         if (id[i] < 0) continue;
479:         y[id[i]] += alpha*x[i];
480:       }
481:     }
482:     y += rstart;
483:     ISRestoreIndices(is,&id);
484:     VecRestoreArray(vfull,&y);
485:     VecRestoreArrayRead(vreduced,&x);
486:   }
487:   return 0;
488: }

490: /*@
491:   VecISCopy - Copies between a reduced vector and the appropriate elements of a full-space vector.

493:   Input Parameters:
494: + vfull    - the full-space vector
495: . is       - the index set for the reduced space
496: . mode     - the direction of copying, SCATTER_FORWARD or SCATTER_REVERSE
497: - vreduced - the reduced-space vector

499:   Output Parameters:
500: . vfull    - the sum of the full-space vector and reduced-space vector

502:   Notes:
503:     The index set identifies entries in the global vector.
504:     Negative indices are skipped; indices outside the ownership range of vfull will raise an error.

506:     mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
507:     mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]

509:   Level: advanced

511: .seealso:  VecAXPY(), VecGetOwnershipRange()
512: @*/
513: PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
514: {
515:   PetscInt       nfull, nreduced;

520:   VecGetSize(vfull, &nfull);
521:   VecGetSize(vreduced, &nreduced);

523:   if (nfull == nreduced) { /* Also takes care of masked vectors */
524:     if (mode == SCATTER_FORWARD) {
525:       VecCopy(vreduced, vfull);
526:     } else {
527:       VecCopy(vfull, vreduced);
528:     }
529:   } else {
530:     const PetscInt *id;
531:     PetscInt        i, n, m, rstart, rend;

533:     ISGetIndices(is, &id);
534:     ISGetLocalSize(is, &n);
535:     VecGetLocalSize(vreduced, &m);
536:     VecGetOwnershipRange(vfull, &rstart, &rend);
538:     if (mode == SCATTER_FORWARD) {
539:       PetscScalar       *y;
540:       const PetscScalar *x;

542:       VecGetArray(vfull, &y);
543:       VecGetArrayRead(vreduced, &x);
544:       y   -= rstart;
545:       for (i = 0; i < n; ++i) {
546:         if (id[i] < 0) continue;
548:         y[id[i]] = x[i];
549:       }
550:       y   += rstart;
551:       VecRestoreArrayRead(vreduced, &x);
552:       VecRestoreArray(vfull, &y);
553:     } else if (mode == SCATTER_REVERSE) {
554:       PetscScalar       *x;
555:       const PetscScalar *y;

557:       VecGetArrayRead(vfull, &y);
558:       VecGetArray(vreduced, &x);
559:       for (i = 0; i < n; ++i) {
560:         if (id[i] < 0) continue;
562:         x[i] = y[id[i]-rstart];
563:       }
564:       VecRestoreArray(vreduced, &x);
565:       VecRestoreArrayRead(vfull, &y);
566:     } else SETERRQ(PetscObjectComm((PetscObject) vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
567:     ISRestoreIndices(is, &id);
568:   }
569:   return 0;
570: }

572: /*@
573:    ISComplementVec - Creates the complement of the index set relative to a layout defined by a Vec

575:    Collective on IS

577:    Input Parameters:
578: +  S -  a PETSc IS
579: -  V - the reference vector space

581:    Output Parameter:
582: .  T -  the complement of S

584:    Level: advanced

586: .seealso: ISCreateGeneral()
587: @*/
588: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
589: {
590:   PetscInt       start, end;

592:   VecGetOwnershipRange(V,&start,&end);
593:   ISComplement(S,start,end,T);
594:   return 0;
595: }

597: /*@
598:    VecISSet - Sets the elements of a vector, specified by an index set, to a constant

600:    Input Parameters:
601: +  V - the vector
602: .  S - index set for the locations in the vector
603: -  c - the constant

605:   Notes:
606:     The index set identifies entries in the global vector.
607:     Negative indices are skipped; indices outside the ownership range of V will raise an error.

609:    Level: advanced

611: .seealso: VecSet(), VecGetOwnershipRange()
612: @*/
613: PetscErrorCode VecISSet(Vec V,IS S, PetscScalar c)
614: {
615:   PetscInt       nloc,low,high,i;
616:   const PetscInt *s;
617:   PetscScalar    *v;

619:   if (!S) return 0; /* simply return with no-op if the index set is NULL */

624:   VecGetOwnershipRange(V,&low,&high);
625:   ISGetLocalSize(S,&nloc);
626:   ISGetIndices(S,&s);
627:   VecGetArray(V,&v);
628:   for (i=0; i<nloc; ++i) {
629:     if (s[i] < 0) continue;
631:     v[s[i]-low] = c;
632:   }
633:   ISRestoreIndices(S,&s);
634:   VecRestoreArray(V,&v);
635:   return 0;
636: }

638: #if !defined(PETSC_USE_COMPLEX)
639: /*@C
640:   VecBoundGradientProjection - Projects  vector according to this definition.
641:   If XL[i] < X[i] < XU[i], then GP[i] = G[i];
642:   If X[i] <= XL[i], then GP[i] = min(G[i],0);
643:   If X[i] >= XU[i], then GP[i] = max(G[i],0);

645:   Input Parameters:
646: + G - current gradient vector
647: . X - current solution vector with XL[i] <= X[i] <= XU[i]
648: . XL - lower bounds
649: - XU - upper bounds

651:   Output Parameter:
652: . GP - gradient projection vector

654:   Notes:
655:     GP may be the same vector as G

657:   Level: advanced
658: @*/
659: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
660: {

662:   PetscInt        n,i;
663:   const PetscReal *xptr,*xlptr,*xuptr;
664:   PetscReal       *gptr,*gpptr;
665:   PetscReal       xval,gpval;

667:   /* Project variables at the lower and upper bound */

674:   VecGetLocalSize(X,&n);

676:   VecGetArrayRead(X,&xptr);
677:   VecGetArrayRead(XL,&xlptr);
678:   VecGetArrayRead(XU,&xuptr);
679:   VecGetArrayPair(G,GP,&gptr,&gpptr);

681:   for (i=0; i<n; ++i) {
682:     gpval = gptr[i]; xval = xptr[i];
683:     if (gpval>0.0 && xval<=xlptr[i]) {
684:       gpval = 0.0;
685:     } else if (gpval<0.0 && xval>=xuptr[i]) {
686:       gpval = 0.0;
687:     }
688:     gpptr[i] = gpval;
689:   }

691:   VecRestoreArrayRead(X,&xptr);
692:   VecRestoreArrayRead(XL,&xlptr);
693:   VecRestoreArrayRead(XU,&xuptr);
694:   VecRestoreArrayPair(G,GP,&gptr,&gpptr);
695:   return 0;
696: }
697: #endif

699: /*@
700:      VecStepMaxBounded - See below

702:      Collective on Vec

704:      Input Parameters:
705: +      X  - vector with no negative entries
706: .      XL - lower bounds
707: .      XU - upper bounds
708: -      DX  - step direction, can have negative, positive or zero entries

710:      Output Parameter:
711: .     stepmax -   minimum value so that X[i] + stepmax*DX[i] <= XL[i]  or  XU[i] <= X[i] + stepmax*DX[i]

713:   Level: intermediate

715: @*/
716: PetscErrorCode VecStepMaxBounded(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *stepmax)
717: {
718:   PetscInt          i,nn;
719:   const PetscScalar *xx,*dx,*xl,*xu;
720:   PetscReal         localmax=0;


727:   VecGetArrayRead(X,&xx);
728:   VecGetArrayRead(XL,&xl);
729:   VecGetArrayRead(XU,&xu);
730:   VecGetArrayRead(DX,&dx);
731:   VecGetLocalSize(X,&nn);
732:   for (i=0;i<nn;i++) {
733:     if (PetscRealPart(dx[i]) > 0) {
734:       localmax=PetscMax(localmax,PetscRealPart((xu[i]-xx[i])/dx[i]));
735:     } else if (PetscRealPart(dx[i])<0) {
736:       localmax=PetscMax(localmax,PetscRealPart((xl[i]-xx[i])/dx[i]));
737:     }
738:   }
739:   VecRestoreArrayRead(X,&xx);
740:   VecRestoreArrayRead(XL,&xl);
741:   VecRestoreArrayRead(XU,&xu);
742:   VecRestoreArrayRead(DX,&dx);
743:   MPIU_Allreduce(&localmax,stepmax,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)X));
744:   return 0;
745: }

747: /*@
748:      VecStepBoundInfo - See below

750:      Collective on Vec

752:      Input Parameters:
753: +      X  - vector with no negative entries
754: .      XL - lower bounds
755: .      XU - upper bounds
756: -      DX  - step direction, can have negative, positive or zero entries

758:      Output Parameters:
759: +     boundmin -  (may be NULL this it is not computed) maximum value so that   XL[i] <= X[i] + boundmax*DX[i] <= XU[i]
760: .     wolfemin -  (may be NULL this it is not computed)
761: -     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]

763:      Notes:
764:     For complex numbers only compares the real part

766:   Level: advanced
767: @*/
768: PetscErrorCode VecStepBoundInfo(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax)
769: {
770:   PetscInt          n,i;
771:   const PetscScalar *x,*xl,*xu,*dx;
772:   PetscReal         t;
773:   PetscReal         localmin=PETSC_INFINITY,localwolfemin=PETSC_INFINITY,localmax=-1;
774:   MPI_Comm          comm;


781:   VecGetArrayRead(X,&x);
782:   VecGetArrayRead(XL,&xl);
783:   VecGetArrayRead(XU,&xu);
784:   VecGetArrayRead(DX,&dx);
785:   VecGetLocalSize(X,&n);
786:   for (i=0; i<n; ++i) {
787:     if (PetscRealPart(dx[i])>0 && PetscRealPart(xu[i]) < PETSC_INFINITY) {
788:       t=PetscRealPart((xu[i]-x[i])/dx[i]);
789:       localmin=PetscMin(t,localmin);
790:       if (localmin>0) {
791:         localwolfemin = PetscMin(t,localwolfemin);
792:       }
793:       localmax = PetscMax(t,localmax);
794:     } else if (PetscRealPart(dx[i])<0 && PetscRealPart(xl[i]) > PETSC_NINFINITY) {
795:       t=PetscRealPart((xl[i]-x[i])/dx[i]);
796:       localmin = PetscMin(t,localmin);
797:       if (localmin>0) {
798:         localwolfemin = PetscMin(t,localwolfemin);
799:       }
800:       localmax = PetscMax(t,localmax);
801:     }
802:   }

804:   VecRestoreArrayRead(X,&x);
805:   VecRestoreArrayRead(XL,&xl);
806:   VecRestoreArrayRead(XU,&xu);
807:   VecRestoreArrayRead(DX,&dx);
808:   PetscObjectGetComm((PetscObject)X,&comm);

810:   if (boundmin) {
811:     MPIU_Allreduce(&localmin,boundmin,1,MPIU_REAL,MPIU_MIN,comm);
812:     PetscInfo(X,"Step Bound Info: Closest Bound: %20.19e\n",(double)*boundmin);
813:   }
814:   if (wolfemin) {
815:     MPIU_Allreduce(&localwolfemin,wolfemin,1,MPIU_REAL,MPIU_MIN,comm);
816:     PetscInfo(X,"Step Bound Info: Wolfe: %20.19e\n",(double)*wolfemin);
817:   }
818:   if (boundmax) {
819:     MPIU_Allreduce(&localmax,boundmax,1,MPIU_REAL,MPIU_MAX,comm);
820:     if (*boundmax < 0) *boundmax=PETSC_INFINITY;
821:     PetscInfo(X,"Step Bound Info: Max: %20.19e\n",(double)*boundmax);
822:   }
823:   return 0;
824: }

826: /*@
827:      VecStepMax - Returns the largest value so that x[i] + step*DX[i] >= 0 for all i

829:      Collective on Vec

831:      Input Parameters:
832: +      X  - vector with no negative entries
833: -      DX  - a step direction, can have negative, positive or zero entries

835:      Output Parameter:
836: .    step - largest value such that x[i] + step*DX[i] >= 0 for all i

838:      Notes:
839:     For complex numbers only compares the real part

841:   Level: advanced
842:  @*/
843: PetscErrorCode VecStepMax(Vec X, Vec DX, PetscReal *step)
844: {
845:   PetscInt          i,nn;
846:   PetscReal         stepmax=PETSC_INFINITY;
847:   const PetscScalar *xx,*dx;


852:   VecGetLocalSize(X,&nn);
853:   VecGetArrayRead(X,&xx);
854:   VecGetArrayRead(DX,&dx);
855:   for (i=0;i<nn;++i) {
857:     else if (PetscRealPart(dx[i])<0) stepmax=PetscMin(stepmax,PetscRealPart(-xx[i]/dx[i]));
858:   }
859:   VecRestoreArrayRead(X,&xx);
860:   VecRestoreArrayRead(DX,&dx);
861:   MPIU_Allreduce(&stepmax,step,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)X));
862:   return 0;
863: }

865: /*@
866:   VecPow - Replaces each component of a vector by x_i^p

868:   Logically Collective on v

870:   Input Parameters:
871: + v - the vector
872: - p - the exponent to use on each element

874:   Level: intermediate

876: @*/
877: PetscErrorCode VecPow(Vec v, PetscScalar p)
878: {
879:   PetscInt       n,i;
880:   PetscScalar    *v1;


884:   VecGetArray(v,&v1);
885:   VecGetLocalSize(v,&n);

887:   if (1.0 == p) {
888:   } else if (-1.0 == p) {
889:     for (i = 0; i < n; ++i) {
890:       v1[i] = 1.0 / v1[i];
891:     }
892:   } else if (0.0 == p) {
893:     for (i = 0; i < n; ++i) {
894:       /*  Not-a-number left alone
895:           Infinity set to one  */
896:       if (v1[i] == v1[i]) {
897:         v1[i] = 1.0;
898:       }
899:     }
900:   } else if (0.5 == p) {
901:     for (i = 0; i < n; ++i) {
902:       if (PetscRealPart(v1[i]) >= 0) {
903:         v1[i] = PetscSqrtScalar(v1[i]);
904:       } else {
905:         v1[i] = PETSC_INFINITY;
906:       }
907:     }
908:   } else if (-0.5 == p) {
909:     for (i = 0; i < n; ++i) {
910:       if (PetscRealPart(v1[i]) >= 0) {
911:         v1[i] = 1.0 / PetscSqrtScalar(v1[i]);
912:       } else {
913:         v1[i] = PETSC_INFINITY;
914:       }
915:     }
916:   } else if (2.0 == p) {
917:     for (i = 0; i < n; ++i) {
918:       v1[i] *= v1[i];
919:     }
920:   } else if (-2.0 == p) {
921:     for (i = 0; i < n; ++i) {
922:       v1[i] = 1.0 / (v1[i] * v1[i]);
923:     }
924:   } else {
925:     for (i = 0; i < n; ++i) {
926:       if (PetscRealPart(v1[i]) >= 0) {
927:         v1[i] = PetscPowScalar(v1[i],p);
928:       } else {
929:         v1[i] = PETSC_INFINITY;
930:       }
931:     }
932:   }
933:   VecRestoreArray(v,&v1);
934:   return 0;
935: }

937: /*@
938:   VecMedian - Computes the componentwise median of three vectors
939:   and stores the result in this vector.  Used primarily for projecting
940:   a vector within upper and lower bounds.

942:   Logically Collective

944:   Input Parameters:
945: + Vec1 - The first vector
946: . Vec2 - The second vector
947: - Vec3 - The third vector

949:   Output Parameter:
950: . VMedian - The median vector (this can be any one of the input vectors)

952:   Level: advanced
953: @*/
954: PetscErrorCode VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
955: {
956:   PetscInt          i,n,low1,high1;
957:   const PetscScalar *v1,*v2,*v3;
958:   PetscScalar       *vmed;


965:   if (Vec1==Vec2 || Vec1==Vec3) {
966:     VecCopy(Vec1,VMedian);
967:     return 0;
968:   }
969:   if (Vec2==Vec3) {
970:     VecCopy(Vec2,VMedian);
971:     return 0;
972:   }

974:   /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
985:   VecCheckSameSize(Vec1,1,Vec2,2);
986:   VecCheckSameSize(Vec1,1,Vec3,3);
987:   VecCheckSameSize(Vec1,1,VMedian,4);

989:   VecGetOwnershipRange(Vec1,&low1,&high1);
990:   VecGetLocalSize(Vec1,&n);
991:   if (n>0) {
992:     VecGetArray(VMedian,&vmed);
993:     if (Vec1 != VMedian) {
994:       VecGetArrayRead(Vec1,&v1);
995:     } else {
996:       v1=vmed;
997:     }
998:     if (Vec2 != VMedian) {
999:       VecGetArrayRead(Vec2,&v2);
1000:     } else {
1001:       v2=vmed;
1002:     }
1003:     if (Vec3 != VMedian) {
1004:       VecGetArrayRead(Vec3,&v3);
1005:     } else {
1006:       v3=vmed;
1007:     }

1009:     for (i=0;i<n;++i) {
1010:       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])));
1011:     }

1013:     VecRestoreArray(VMedian,&vmed);
1014:     if (VMedian != Vec1) {
1015:       VecRestoreArrayRead(Vec1,&v1);
1016:     }
1017:     if (VMedian != Vec2) {
1018:       VecRestoreArrayRead(Vec2,&v2);
1019:     }
1020:     if (VMedian != Vec3) {
1021:       VecRestoreArrayRead(Vec3,&v3);
1022:     }
1023:   }
1024:   return 0;
1025: }