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:   PetscErrorCode    ierr;
 23:   PetscInt          i,n_same=0;
 24:   PetscInt          n,low,high;
 25:   PetscInt          *same=NULL;
 26:   const PetscScalar *v1,*v2;

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

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

 45:     PetscMalloc1(n,&same);

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

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

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

 66:   Collective on S

 68:   Input Parameters:
 69: . Vec1, Vec2 - the two vectors to compare

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

 74:   Notes:
 75:   The two vectors must have the same parallel layout

 77:   For complex numbers this only compares the real part

 79:   Level: advanced
 80: @*/
 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;

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

 95:   VecGetOwnershipRange(Vec1,&low,&high);
 96:   VecGetLocalSize(Vec1,&n);
 97:   if (n>0){
 98:     if (Vec1 == Vec2){
 99:       VecGetArrayRead(Vec1,&v1);
100:       v2=v1;
101:     } else {
102:       VecGetArrayRead(Vec1,&v1);
103:       VecGetArrayRead(Vec2,&v2);
104:     }

106:     PetscMalloc1(n,&lt);

108:     for (i=0; i<n; ++i){
109:       if (PetscRealPart(v1[i]) < PetscRealPart(v2[i])) {lt[n_lt]=low+i; ++n_lt;}
110:     }

112:     if (Vec1 == Vec2){
113:       VecRestoreArrayRead(Vec1,&v1);
114:     } else {
115:       VecRestoreArrayRead(Vec1,&v1);
116:       VecRestoreArrayRead(Vec2,&v2);
117:     }
118:   }
119:   ISCreateGeneral(PetscObjectComm((PetscObject)Vec1),n_lt,lt,PETSC_OWN_POINTER,S);
120:   return(0);
121: }

123: /*@
124:   VecWhichGreaterThan - Creates an index set containing the indices
125:   where the vectors Vec1 > Vec2

127:   Collective on S

129:   Input Parameters:
130: . Vec1, Vec2 - the two vectors to compare

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

135:   Notes:
136:   The two vectors must have the same parallel layout

138:   For complex numbers this only compares the real part

140:   Level: advanced
141: @*/
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;

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

156:   VecGetOwnershipRange(Vec1,&low,&high);
157:   VecGetLocalSize(Vec1,&n);
158:   if (n>0){
159:     if (Vec1 == Vec2){
160:       VecGetArrayRead(Vec1,&v1);
161:       v2=v1;
162:     } else {
163:       VecGetArrayRead(Vec1,&v1);
164:       VecGetArrayRead(Vec2,&v2);
165:     }

167:     PetscMalloc1(n,&gt);

169:     for (i=0; i<n; ++i){
170:       if (PetscRealPart(v1[i]) > PetscRealPart(v2[i])) {gt[n_gt]=low+i; ++n_gt;}
171:     }

173:     if (Vec1 == Vec2){
174:       VecRestoreArrayRead(Vec1,&v1);
175:     } else {
176:       VecRestoreArrayRead(Vec1,&v1);
177:       VecRestoreArrayRead(Vec2,&v2);
178:     }
179:   }
180:   ISCreateGeneral(PetscObjectComm((PetscObject)Vec1),n_gt,gt,PETSC_OWN_POINTER,S);
181:   return(0);
182: }

184: /*@
185:   VecWhichBetween - Creates an index set containing the indices
186:                where  VecLow < V < VecHigh

188:   Collective on S

190:   Input Parameters:
191: + VecLow - lower bound
192: . V - Vector to compare
193: - VecHigh - higher bound

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

198:   Notes:
199:   The vectors must have the same parallel layout

201:   For complex numbers this only compares the real part

203:   Level: advanced
204: @*/
205: PetscErrorCode VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS *S)
206: {

208:   PetscErrorCode    ierr;
209:   PetscInt          i,n_vm=0;
210:   PetscInt          n,low,high;
211:   PetscInt          *vm=NULL;
212:   const PetscScalar *v1,*v2,*vmiddle;

220:   VecCheckSameSize(V,2,VecLow,1);
221:   VecCheckSameSize(V,2,VecHigh,3);

223:   VecGetOwnershipRange(VecLow,&low,&high);
224:   VecGetLocalSize(VecLow,&n);
225:   if (n>0){
226:     VecGetArrayRead(VecLow,&v1);
227:     if (VecLow != VecHigh){
228:       VecGetArrayRead(VecHigh,&v2);
229:     } else {
230:       v2=v1;
231:     }
232:     if (V != VecLow && V != VecHigh){
233:       VecGetArrayRead(V,&vmiddle);
234:     } else if (V==VecLow){
235:       vmiddle=v1;
236:     } else {
237:       vmiddle=v2;
238:     }

240:     PetscMalloc1(n,&vm);

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:     }

246:     VecRestoreArrayRead(VecLow,&v1);
247:     if (VecLow != VecHigh){
248:       VecRestoreArrayRead(VecHigh,&v2);
249:     }
250:     if (V != VecLow && V != VecHigh){
251:       VecRestoreArrayRead(V,&vmiddle);
252:     }
253:   }
254:   ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
255:   return(0);
256: }


259: /*@
260:   VecWhichBetweenOrEqual - Creates an index set containing the indices
261:   where  VecLow <= V <= VecHigh

263:   Collective on S

265:   Input Parameters:
266: + VecLow - lower bound
267: . V - Vector to compare
268: - VecHigh - higher bound

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

273:   Level: advanced
274: @*/

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;

290:   VecCheckSameSize(V,2,VecLow,1);
291:   VecCheckSameSize(V,2,VecHigh,3);

293:   VecGetOwnershipRange(VecLow,&low,&high);
294:   VecGetLocalSize(VecLow,&n);
295:   if (n>0){
296:     VecGetArrayRead(VecLow,&v1);
297:     if (VecLow != VecHigh){
298:       VecGetArrayRead(VecHigh,&v2);
299:     } else {
300:       v2=v1;
301:     }
302:     if (V != VecLow && V != VecHigh){
303:       VecGetArrayRead(V,&vmiddle);
304:     } else if (V==VecLow){
305:       vmiddle=v1;
306:     } else {
307:       vmiddle =v2;
308:     }

310:     PetscMalloc1(n,&vm);

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:     }

316:     VecRestoreArrayRead(VecLow,&v1);
317:     if (VecLow != VecHigh){
318:       VecRestoreArrayRead(VecHigh,&v2);
319:     }
320:     if (V != VecLow && V != VecHigh){
321:       VecRestoreArrayRead(V,&vmiddle);
322:     }
323:   }
324:   ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
325:   return(0);
326: }

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)

335:   Collective on S

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

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

347:   Level: advanced
348: @*/

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;

366:   VecCheckSameSize(V,2,VecLow,1);
367:   VecCheckSameSize(V,2,D,3);
368:   VecCheckSameSize(V,2,VecHigh,4);

370:   VecGetOwnershipRange(VecLow,&low,&high);
371:   VecGetLocalSize(VecLow,&n);
372:   if (n>0){
373:     VecGetArrayRead(VecLow,&v1);
374:     if (VecLow != VecHigh){
375:       VecGetArrayRead(VecHigh,&v2);
376:     } else {
377:       v2=v1;
378:     }
379:     if (V != VecLow && V != VecHigh){
380:       VecGetArrayRead(V,&v);
381:     } else if (V==VecLow){
382:       v = v1;
383:     } else {
384:       v = v2;
385:     }
386:     if (D != VecLow && D != VecHigh && D != V){
387:       VecGetArrayRead(D,&d);
388:     } else if (D==VecLow){
389:       d = v1;
390:     } else if (D==VecHigh){
391:       d = v2;
392:     } else {
393:       d = v;
394:     }

396:     PetscMalloc1(n,&vm);

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:     }

420:     VecRestoreArrayRead(VecLow,&v1);
421:     if (VecLow != VecHigh){
422:       VecRestoreArrayRead(VecHigh,&v2);
423:     }
424:     if (V != VecLow && V != VecHigh){
425:       VecRestoreArrayRead(V,&v);
426:     }
427:     if (D != VecLow && D != VecHigh && D != V){
428:       VecRestoreArrayRead(D,&d);
429:     }
430:   }
431:   ISCreateGeneral(PetscObjectComm((PetscObject)V),n_vm,vm,PETSC_OWN_POINTER,S);
432:   return(0);
433: }


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

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

446:   Output Parameters:
447: . vfull    - the sum of the full-space vector and reduced-space vector

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.

453:   Level: advanced

455: .seealso:  VecAXPY(), VecGetOwnershipRange()
456: @*/
457: PetscErrorCode VecISAXPY(Vec vfull, IS is, PetscScalar alpha, Vec vreduced)
458: {
459:   PetscInt       nfull,nreduced;

466:   VecGetSize(vfull,&nfull);
467:   VecGetSize(vreduced,&nreduced);

469:   if (nfull == nreduced) { /* Also takes care of masked vectors */
470:     VecAXPY(vfull,alpha,vreduced);
471:   } else {
472:     PetscScalar      *y;
473:     const PetscScalar *x;
474:     PetscInt          i,n,m,rstart,rend;
475:     const PetscInt    *id;

477:     VecGetArray(vfull,&y);
478:     VecGetArrayRead(vreduced,&x);
479:     ISGetIndices(is,&id);
480:     ISGetLocalSize(is,&n);
481:     VecGetLocalSize(vreduced,&m);
482:     if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"IS local length not equal to Vec local length");
483:     VecGetOwnershipRange(vfull,&rstart,&rend);
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:     ISRestoreIndices(is,&id);
500:     VecRestoreArray(vfull,&y);
501:     VecRestoreArrayRead(vreduced,&x);
502:   }
503:   return(0);
504: }

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

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

515:   Output Parameters:
516: . vfull    - the sum of the full-space vector and reduced-space vector

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.

522:     mode == SCATTER_FORWARD: vfull[is[i]] = vreduced[i]
523:     mode == SCATTER_REVERSE: vreduced[i] = vfull[is[i]]

525:   Level: advanced

527: .seealso:  VecAXPY(), VecGetOwnershipRange()
528: @*/
529: PetscErrorCode VecISCopy(Vec vfull, IS is, ScatterMode mode, Vec vreduced)
530: {
531:   PetscInt       nfull, nreduced;

538:   VecGetSize(vfull, &nfull);
539:   VecGetSize(vreduced, &nreduced);

541:   if (nfull == nreduced) { /* Also takes care of masked vectors */
542:     if (mode == SCATTER_FORWARD) {
543:       VecCopy(vreduced, vfull);
544:     } else {
545:       VecCopy(vfull, vreduced);
546:     }
547:   } else {
548:     const PetscInt *id;
549:     PetscInt        i, n, m, rstart, rend;

551:     ISGetIndices(is, &id);
552:     ISGetLocalSize(is, &n);
553:     VecGetLocalSize(vreduced, &m);
554:     VecGetOwnershipRange(vfull, &rstart, &rend);
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;

560:       VecGetArray(vfull, &y);
561:       VecGetArrayRead(vreduced, &x);
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:       VecRestoreArrayRead(vreduced, &x);
570:       VecRestoreArray(vfull, &y);
571:     } else if (mode == SCATTER_REVERSE) {
572:       PetscScalar       *x;
573:       const PetscScalar *y;

575:       VecGetArrayRead(vfull, &y);
576:       VecGetArray(vreduced, &x);
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:       VecRestoreArray(vreduced, &x);
583:       VecRestoreArrayRead(vfull, &y);
584:     } else SETERRQ(PetscObjectComm((PetscObject) vfull), PETSC_ERR_ARG_WRONG, "Only forward or reverse modes are legal");
585:     ISRestoreIndices(is, &id);
586:   }
587:   return(0);
588: }

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

593:    Collective on IS

595:    Input Parameter:
596: +  S -  a PETSc IS
597: -  V - the reference vector space

599:    Output Parameter:
600: .  T -  the complement of S

602:    Level: advanced

604: .seealso: ISCreateGeneral()
605: @*/
606: PetscErrorCode ISComplementVec(IS S, Vec V, IS *T)
607: {
609:   PetscInt       start, end;

612:   VecGetOwnershipRange(V,&start,&end);
613:   ISComplement(S,start,end,T);
614:   return(0);
615: }

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

620:    Input Parameter:
621: +  V - the vector
622: .  S - index set for the locations in the vector
623: -  c - the constant

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.

629:    Level: advanced

631: .seealso: VecSet(), VecGetOwnershipRange()
632: @*/
633: PetscErrorCode VecISSet(Vec V,IS S, PetscScalar c)
634: {
636:   PetscInt       nloc,low,high,i;
637:   const PetscInt *s;
638:   PetscScalar    *v;

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

646:   VecGetOwnershipRange(V,&low,&high);
647:   ISGetLocalSize(S,&nloc);
648:   ISGetIndices(S,&s);
649:   VecGetArray(V,&v);
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:   ISRestoreIndices(S,&s);
656:   VecRestoreArray(V,&v);
657:   return(0);
658: }

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);

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

673:   Output Parameter:
674: . GP - gradient projection vector

676:   Notes:
677:     GP may be the same vector as G

679:   Level: advanced
680: @*/
681: PetscErrorCode VecBoundGradientProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP)
682: {

684:   PetscErrorCode  ierr;
685:   PetscInt        n,i;
686:   const PetscReal *xptr,*xlptr,*xuptr;
687:   PetscReal       *gptr,*gpptr;
688:   PetscReal       xval,gpval;

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

698:   VecGetLocalSize(X,&n);

700:   VecGetArrayRead(X,&xptr);
701:   VecGetArrayRead(XL,&xlptr);
702:   VecGetArrayRead(XU,&xuptr);
703:   VecGetArrayPair(G,GP,&gptr,&gpptr);

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:   }

715:   VecRestoreArrayRead(X,&xptr);
716:   VecRestoreArrayRead(XL,&xlptr);
717:   VecRestoreArrayRead(XU,&xuptr);
718:   VecRestoreArrayPair(G,GP,&gptr,&gpptr);
719:   return(0);
720: }
721: #endif

723: /*@
724:      VecStepMaxBounded - See below

726:      Collective on Vec

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

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

737:   Level: intermediate

739: @*/
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;


753:   VecGetArrayRead(X,&xx);
754:   VecGetArrayRead(XL,&xl);
755:   VecGetArrayRead(XU,&xu);
756:   VecGetArrayRead(DX,&dx);
757:   VecGetLocalSize(X,&nn);
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:   VecRestoreArrayRead(X,&xx);
766:   VecRestoreArrayRead(XL,&xl);
767:   VecRestoreArrayRead(XU,&xu);
768:   VecRestoreArrayRead(DX,&dx);
769:   MPIU_Allreduce(&localmax,stepmax,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)X));
770:   return(0);
771: }

773: /*@
774:      VecStepBoundInfo - See below

776:      Collective on Vec

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

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]

789:      Notes:
790:     For complex numbers only compares the real part

792:   Level: advanced
793: @*/
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;


809:   VecGetArrayRead(X,&x);
810:   VecGetArrayRead(XL,&xl);
811:   VecGetArrayRead(XU,&xu);
812:   VecGetArrayRead(DX,&dx);
813:   VecGetLocalSize(X,&n);
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:   }

832:   VecRestoreArrayRead(X,&x);
833:   VecRestoreArrayRead(XL,&xl);
834:   VecRestoreArrayRead(XU,&xu);
835:   VecRestoreArrayRead(DX,&dx);
836:   PetscObjectGetComm((PetscObject)X,&comm);

838:   if (boundmin){
839:     MPIU_Allreduce(&localmin,boundmin,1,MPIU_REAL,MPIU_MIN,comm);
840:     PetscInfo1(X,"Step Bound Info: Closest Bound: %20.19e\n",(double)*boundmin);
841:   }
842:   if (wolfemin){
843:     MPIU_Allreduce(&localwolfemin,wolfemin,1,MPIU_REAL,MPIU_MIN,comm);
844:     PetscInfo1(X,"Step Bound Info: Wolfe: %20.19e\n",(double)*wolfemin);
845:   }
846:   if (boundmax) {
847:     MPIU_Allreduce(&localmax,boundmax,1,MPIU_REAL,MPIU_MAX,comm);
848:     if (*boundmax < 0) *boundmax=PETSC_INFINITY;
849:     PetscInfo1(X,"Step Bound Info: Max: %20.19e\n",(double)*boundmax);
850:   }
851:   return(0);
852: }

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

857:      Collective on Vec

859:      Input Parameters:
860: +      X  - vector with no negative entries
861: -      DX  - a step direction, can have negative, positive or zero entries

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

866:      Notes:
867:     For complex numbers only compares the real part

869:   Level: advanced
870:  @*/
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;


882:   VecGetLocalSize(X,&nn);
883:   VecGetArrayRead(X,&xx);
884:   VecGetArrayRead(DX,&dx);
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:   VecRestoreArrayRead(X,&xx);
890:   VecRestoreArrayRead(DX,&dx);
891:   MPIU_Allreduce(&stepmax,step,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)X));
892:   return(0);
893: }

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

898:   Logically Collective on v

900:   Input Parameter:
901: + v - the vector
902: - p - the exponent to use on each element

904:   Output Parameter:
905: . v - the vector

907:   Level: intermediate

909: @*/
910: PetscErrorCode VecPow(Vec v, PetscScalar p)
911: {
913:   PetscInt       n,i;
914:   PetscScalar    *v1;


919:   VecGetArray(v,&v1);
920:   VecGetLocalSize(v,&n);

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:   VecRestoreArray(v,&v1);
969:   return(0);
970: }

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.

977:   Logically Collective

979:   Input Parameters:
980: . Vec1, Vec2, Vec3 - The three vectors

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

985:   Level: advanced
986: @*/
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;


1000:   if (Vec1==Vec2 || Vec1==Vec3){
1001:     VecCopy(Vec1,VMedian);
1002:     return(0);
1003:   }
1004:   if (Vec2==Vec3){
1005:     VecCopy(Vec2,VMedian);
1006:     return(0);
1007:   }

1009:   /* Assert that Vec1 != Vec2 and Vec2 != Vec3 */
1020:   VecCheckSameSize(Vec1,1,Vec2,2);
1021:   VecCheckSameSize(Vec1,1,Vec3,3);
1022:   VecCheckSameSize(Vec1,1,VMedian,4);

1024:   VecGetOwnershipRange(Vec1,&low1,&high1);
1025:   VecGetLocalSize(Vec1,&n);
1026:   if (n>0){
1027:     VecGetArray(VMedian,&vmed);
1028:     if (Vec1 != VMedian){
1029:       VecGetArrayRead(Vec1,&v1);
1030:     } else {
1031:       v1=vmed;
1032:     }
1033:     if (Vec2 != VMedian){
1034:       VecGetArrayRead(Vec2,&v2);
1035:     } else {
1036:       v2=vmed;
1037:     }
1038:     if (Vec3 != VMedian){
1039:       VecGetArrayRead(Vec3,&v3);
1040:     } else {
1041:       v3=vmed;
1042:     }

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:     }

1048:     VecRestoreArray(VMedian,&vmed);
1049:     if (VMedian != Vec1){
1050:       VecRestoreArrayRead(Vec1,&v1);
1051:     }
1052:     if (VMedian != Vec2){
1053:       VecRestoreArrayRead(Vec2,&v2);
1054:     }
1055:     if (VMedian != Vec3){
1056:       VecRestoreArrayRead(Vec3,&v3);
1057:     }
1058:   }
1059:   return(0);
1060: }