Actual source code: pvec.c

  1: #include "tao_general.h"
  2: #include "private/vecimpl.h"    /*I "tao_solver.h" I*/
  3: #include "petscksp.h"



  9: int  VecCompare(Vec V1,Vec V2, PetscTruth *flg){
 10:   int info;
 11:   PetscInt n1,n2,N1,N2;

 16:   info = VecGetSize(V1,&N1);CHKERRQ(info);
 17:   info = VecGetSize(V2,&N2);CHKERRQ(info);
 18:   info = VecGetLocalSize(V1,&n1);CHKERRQ(info);
 19:   info = VecGetLocalSize(V2,&n2);CHKERRQ(info);
 20:   if (N1==N2 && n1==n2) 
 21:     *flg=PETSC_TRUE;
 22:   else
 23:     *flg=PETSC_FALSE;

 25:   return(0);
 26: }

 28: /* ---------------------------------------------------------- */
 31: int VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
 32: {
 33:   int i,info;
 34:   PetscInt n,low1,low2,low3,low4,high1,high2,high3,high4;
 35:   PetscScalar *v1,*v2,*v3,*vmed;


 43:   if (Vec1==Vec2 || Vec1==Vec3){
 44:     info=VecCopy(Vec1,VMedian); CHKERRQ(info); 
 45:     return(0);
 46:   }
 47:   if (Vec2==Vec3){
 48:     info=VecCopy(Vec2,VMedian); CHKERRQ(info); 
 49:     return(0);
 50:   }


 58:   info = VecGetOwnershipRange(Vec1, &low1, &high1); CHKERRQ(info);
 59:   info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);
 60:   info = VecGetOwnershipRange(Vec3, &low3, &high3); CHKERRQ(info);
 61:   info = VecGetOwnershipRange(VMedian, &low4, &high4); CHKERRQ(info);
 62:   if ( low1!= low2 || low1!= low3 || low1!= low4 ||
 63:        high1!= high2 || high1!= high3 || high1!= high4)
 64:     SETERRQ(1,"InCompatible vector local lengths");

 66:   info = VecGetArray(Vec1,&v1); CHKERRQ(info);
 67:   info = VecGetArray(Vec2,&v2); CHKERRQ(info);
 68:   info = VecGetArray(Vec3,&v3); CHKERRQ(info);

 70:   if ( VMedian != Vec1 && VMedian != Vec2 && VMedian != Vec3){
 71:     info = VecGetArray(VMedian,&vmed); CHKERRQ(info);
 72:   } else if ( VMedian==Vec1 ){
 73:     vmed=v1;
 74:   } else if ( VMedian==Vec2 ){
 75:     vmed=v2;
 76:   } else {
 77:     vmed=v3;
 78:   }

 80:   info=VecGetLocalSize(Vec1,&n); CHKERRQ(info);

 82:   for (i=0;i<n;i++){
 83:     vmed[i]=PetscMax(PetscMax(PetscMin(v1[i],v2[i]),PetscMin(v1[i],v3[i])),PetscMin(v2[i],v3[i]));
 84:   }

 86:   info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
 87:   info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);
 88:   info = VecRestoreArray(Vec3,&v2); CHKERRQ(info);
 89:   
 90:   if (VMedian!=Vec1 && VMedian != Vec2 && VMedian != Vec3){
 91:     info = VecRestoreArray(VMedian,&vmed); CHKERRQ(info);
 92:   }

 94:   return(0);
 95: }

 97: /* ---------------------------------------------------------- 
100: int VecPointwiseMin(Vec Vec1, Vec Vec2, Vec VMin)
101: {
102:   int i,n,low1,low2,low3,high1,high2,high3,info;
103:   PetscScalar *v1,*v2,*vmin;


110:   if (Vec1==Vec2){
111:     info=VecCopy(Vec1,VMin); CHKERRQ(info); 
112:     return(0);
113:   }


118:   info = VecGetOwnershipRange(Vec1, &low1, &high1); CHKERRQ(info);
119:   info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);
120:   info = VecGetOwnershipRange(VMin, &low3, &high3); CHKERRQ(info);

122:   if ( low1!= low2 || low1!= low3 || high1!= high2 || high1!= high3)
123:     SETERRQ(1,"InCompatible vector local lengths");

125:   info = VecGetArray(Vec1,&v1); CHKERRQ(info);
126:   info = VecGetArray(Vec2,&v2); CHKERRQ(info);

128:   if ( VMin != Vec1 && VMin != Vec2){
129:     info = VecGetArray(VMin,&vmin); CHKERRQ(info);
130:   } else if ( VMin==Vec1 ){
131:     vmin=v1;
132:   } else {
133:     vmin =v2;
134:   }

136:   info=VecGetLocalSize(Vec1,&n); CHKERRQ(info);

138:   for (i=0; i<n; i++){
139:     vmin[i] = PetscMin(v1[i],v2[i]);
140:   }

142:   info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
143:   info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);

145:   if (VMin!=Vec1 && VMin != Vec2){
146:     info = VecRestoreArray(VMin,&vmin); CHKERRQ(info);
147:   }

149:   return(0);
150: }
151: */
152: /* ---------------------------------------------------------- 
155: int VecPointwiseMax(Vec Vec1, Vec Vec2, Vec VMax)
156: {
157:   int i,n,low1,low2,low3,high1,high2,high3,info;
158:   PetscScalar *v1,*v2,*vmax;


162:   if (Vec1==Vec2){
163:     info=VecCopy(Vec1,VMax); CHKERRQ(info); 
164:     return(0);
165:   }


174:   info = VecGetOwnershipRange(Vec1, &low1, &high1); CHKERRQ(info);
175:   info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);
176:   info = VecGetOwnershipRange(VMax, &low3, &high3); CHKERRQ(info);

178:   if ( low1!= low2 || low1!= low3 || high1!= high2 || high1!= high3)
179:     SETERRQ(1,"Vectors must be identically loaded over processors");

181:   info = VecGetArray(Vec1,&v1); CHKERRQ(info);
182:   info = VecGetArray(Vec2,&v2); CHKERRQ(info);

184:   if ( VMax != Vec1 && VMax != Vec2){
185:     info = VecGetArray(VMax,&vmax); CHKERRQ(info);
186:   } else if ( VMax==Vec1 ){ vmax=v1;
187:   } else { vmax =v2;
188:   }

190:   info=VecGetLocalSize(Vec1,&n); CHKERRQ(info);

192:   for (i=0; i<n; i++){
193:     vmax[i] = PetscMax(v1[i],v2[i]);
194:   }

196:   info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
197:   info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);

199:   if ( VMax != Vec1 && VMax != Vec2){
200:     info = VecRestoreArray(VMax,&vmax); CHKERRQ(info);
201:   }

203:   return(0);
204: }*/

208: inline static PetscScalar Fischer(PetscScalar a, PetscScalar b)
209: {
210:    // Method suggested by Bob Vanderbei
211:    if (a + b <= 0) {
212:      return sqrt(a*a + b*b) - (a + b);
213:    }
214:    return -2.0*a*b / (sqrt(a*a + b*b) + (a + b));
215: }

219: int VecFischer(Vec X, Vec F, Vec L, Vec U, Vec FF)
220: {
221:   PetscScalar *x, *f, *l, *u, *ff;
222:   PetscScalar xval, fval, lval, uval;
223:   int info, i;
224:   PetscInt low[5], high[5], n;


233:   info = VecGetOwnershipRange(X, low, high); CHKERRQ(info);
234:   info = VecGetOwnershipRange(F, low + 1, high + 1); CHKERRQ(info);
235:   info = VecGetOwnershipRange(L, low + 2, high + 2); CHKERRQ(info);
236:   info = VecGetOwnershipRange(U, low + 3, high + 3); CHKERRQ(info);
237:   info = VecGetOwnershipRange(FF, low + 4, high + 4); CHKERRQ(info);

239:   for (i = 1; i < 4; ++i) {
240:     if (low[0] != low[i] || high[0] != high[i])
241:       SETERRQ(1,"Vectors must be identically loaded over processors");
242:   }

244:   info = VecGetArray(X, &x); CHKERRQ(info);
245:   info = VecGetArray(F, &f); CHKERRQ(info);
246:   info = VecGetArray(L, &l); CHKERRQ(info);
247:   info = VecGetArray(U, &u); CHKERRQ(info);
248:   info = VecGetArray(FF, &ff); CHKERRQ(info);

250:   info = VecGetLocalSize(X, &n); CHKERRQ(info);

252:   for (i = 0; i < n; ++i) {
253:     xval = x[i]; fval = f[i];
254:     lval = l[i]; uval = u[i];

256:     if ((lval <= -TAO_INFINITY) && (uval >= TAO_INFINITY)) {
257:       ff[i] = -fval;
258:     } 
259:     else if (lval <= -TAO_INFINITY) {
260:       ff[i] = -Fischer(uval - xval, -fval);
261:     } 
262:     else if (uval >=  TAO_INFINITY) {
263:       ff[i] =  Fischer(xval - lval,  fval);
264:     } 
265:     else if (lval == uval) {
266:       ff[i] = lval - xval;
267:     }
268:     else {
269:       fval  =  Fischer(uval - xval, -fval);
270:       ff[i] =  Fischer(xval - lval,  fval);
271:     }
272:   }
273:   
274:   info = VecRestoreArray(X, &x); CHKERRQ(info);
275:   info = VecRestoreArray(F, &f); CHKERRQ(info);
276:   info = VecRestoreArray(L, &l); CHKERRQ(info);
277:   info = VecRestoreArray(U, &u); CHKERRQ(info);
278:   info = VecRestoreArray(FF, &ff); CHKERRQ(info);

280:   return(0);
281: }

285: inline static PetscScalar SFischer(PetscScalar a, PetscScalar b, PetscScalar c)
286: {
287:    // Method suggested by Bob Vanderbei
288:    if (a + b <= 0) {
289:      return sqrt(a*a + b*b + 2.0*c*c) - (a + b);
290:    }
291:    return 2.0*(c*c - a*b) / (sqrt(a*a + b*b + 2.0*c*c) + (a + b));
292: }

296: int VecSFischer(Vec X, Vec F, Vec L, Vec U, double mu, Vec FF)
297: {
298:   PetscScalar *x, *f, *l, *u, *ff;
299:   PetscScalar xval, fval, lval, uval;
300:   int info, i;
301:   PetscInt low[5], high[5], n;


310:   info = VecGetOwnershipRange(X, low, high); CHKERRQ(info);
311:   info = VecGetOwnershipRange(F, low + 1, high + 1); CHKERRQ(info);
312:   info = VecGetOwnershipRange(L, low + 2, high + 2); CHKERRQ(info);
313:   info = VecGetOwnershipRange(U, low + 3, high + 3); CHKERRQ(info);
314:   info = VecGetOwnershipRange(FF, low + 4, high + 4); CHKERRQ(info);

316:   for (i = 1; i < 4; ++i) {
317:     if (low[0] != low[i] || high[0] != high[i])
318:       SETERRQ(1,"Vectors must be identically loaded over processors");
319:   }

321:   info = VecGetArray(X, &x); CHKERRQ(info);
322:   info = VecGetArray(F, &f); CHKERRQ(info);
323:   info = VecGetArray(L, &l); CHKERRQ(info);
324:   info = VecGetArray(U, &u); CHKERRQ(info);
325:   info = VecGetArray(FF, &ff); CHKERRQ(info);

327:   info = VecGetLocalSize(X, &n); CHKERRQ(info);

329:   for (i = 0; i < n; ++i) {
330:     xval = (*x++); fval = (*f++);
331:     lval = (*l++); uval = (*u++);

333:     if ((lval <= -TAO_INFINITY) && (uval >= TAO_INFINITY)) {
334:       (*ff++) = -fval - mu*xval;
335:     } 
336:     else if (lval <= -TAO_INFINITY) {
337:       (*ff++) = -SFischer(uval - xval, -fval, mu);
338:     } 
339:     else if (uval >=  TAO_INFINITY) {
340:       (*ff++) =  SFischer(xval - lval,  fval, mu);
341:     } 
342:     else if (lval == uval) {
343:       (*ff++) = lval - xval;
344:     } 
345:     else {
346:       fval    =  SFischer(uval - xval, -fval, mu);
347:       (*ff++) =  SFischer(xval - lval,  fval, mu);
348:     }
349:   }
350:   x -= n; f -= n; l -=n; u -= n; ff -= n;

352:   info = VecRestoreArray(X, &x); CHKERRQ(info);
353:   info = VecRestoreArray(F, &f); CHKERRQ(info);
354:   info = VecRestoreArray(L, &l); CHKERRQ(info);
355:   info = VecRestoreArray(U, &u); CHKERRQ(info);
356:   info = VecRestoreArray(FF, &ff); CHKERRQ(info);
357:   return(0);
358: }

362: int VecBoundProjection(Vec G, Vec X, Vec XL, Vec XU, Vec GP){

364:   int i,info;
365:   PetscInt n;
366:   PetscScalar *xptr,*xlptr,*xuptr,*gptr,*gpptr;
367:   PetscScalar xval,gpval;

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

372:   info = VecGetLocalSize(X,&n); CHKERRQ(info);

374:   info=VecGetArray(X,&xptr); CHKERRQ(info);
375:   info=VecGetArray(XL,&xlptr); CHKERRQ(info);
376:   info=VecGetArray(XU,&xuptr); CHKERRQ(info);
377:   info=VecGetArray(G,&gptr); CHKERRQ(info);
378:   if (G!=GP){
379:     info=VecGetArray(GP,&gpptr); CHKERRQ(info);
380:   } else { gpptr=gptr; }

382:   for (i=0; i<n; ++i){
383:     gpval = gptr[i]; xval = xptr[i]; 

385:     if (gpval>0 && xval<=xlptr[i]){
386:       gpval = 0;
387:     } else if (gpval<0 && xval>=xuptr[i]){
388:       gpval = 0;
389:     }
390:     gpptr[i] = gpval;

392:     /*
393:     if (xlptr[i] >= xuptr[i]){
394:       gpptr[i]=0.0;
395:     } else if (xptr[i] <= xlptr[i]+eps){
396:       gpptr[i] = PetscMin(gptr[i],zero);
397:     } else if (xptr[i] >= xuptr[i]-eps){
398:       gpptr[i] = PetscMax(gptr[i],zero);
399:     } else {
400:       gpptr[i] = gptr[i];
401:     }
402:     */
403:   }

405:   info=VecRestoreArray(X,&xptr); CHKERRQ(info);
406:   info=VecRestoreArray(XL,&xlptr); CHKERRQ(info);
407:   info=VecRestoreArray(XU,&xuptr); CHKERRQ(info);
408:   info=VecRestoreArray(G,&gptr); CHKERRQ(info);
409:   if (G!=GP){
410:     info=VecRestoreArray(GP,&gpptr); CHKERRQ(info);
411:   }
412:   return(0);
413: }

417: int VecISAXPY(Vec Y, PetscScalar alpha, Vec X, IS YIS){

419:   int i,info;
420:   PetscInt in1,in2,xlow,xhigh,ylow,yhigh;
421:   const PetscInt *yis;
422:   PetscScalar *x,*y;

425:   info=ISGetLocalSize(YIS,&in1); CHKERRQ(info);
426:   info=VecGetLocalSize(X,&in2); CHKERRQ(info);

428:   if ( in1 != in2 )
429:     SETERRQ(1,"Index set and X vector must be identically loaded over processors");

431:   info=VecGetOwnershipRange(X, &xlow, &xhigh); CHKERRQ(info);
432:   info=VecGetOwnershipRange(Y, &ylow, &yhigh); CHKERRQ(info);

434:   info=ISGetIndices(YIS, &yis); CHKERRQ(info);

436:   info=VecGetArray(X,&x); CHKERRQ(info);
437:   if (X != Y){
438:     info=VecGetArray(Y,&y); CHKERRQ(info);
439:   } else {
440:     y=x;
441:   }

443:   for (i=0; i<in1; i++){
444:     if (yis[i]>=ylow && yis[i]<yhigh && i+xlow < xhigh){
445:       y[yis[i]-ylow] = y[yis[i]-ylow] + (alpha)*x[i];
446:     } else {
447:       SETERRQ(1,"IS index out of range");
448:     }
449:   }
450:   
451:   info=ISRestoreIndices(YIS, &yis); CHKERRQ(info);

453:   info=VecRestoreArray(X,&x); CHKERRQ(info);
454:   if (X != Y){
455:     info=VecRestoreArray(Y,&y); CHKERRQ(info);
456:   }

458:   info = PetscLogFlops(2*in1); CHKERRQ(info);

460:   return(0);
461: }

465: int VecCreateSubVec(Vec A,IS S,Vec *B)
466: {
467:   int     i,info;
468:   PetscInt nloc,n,low,high,nnloc;
469:   const PetscInt *ptrS;
470:   PetscScalar  zero=0.0;
471:   PetscScalar *ptrB,*ptrA;
472:   const VecType type_name;
473:   MPI_Comm comm;

478:   info=ISGetLocalSize(S,&nloc);CHKERRQ(info);
479:   info=ISGetSize(S,&n);CHKERRQ(info);
480:   info=PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(info);
481:   info=VecCreate(comm,B);CHKERRQ(info);
482:   info=VecSetSizes(*B,nloc,n);CHKERRQ(info);
483:   info=VecGetType(A,&type_name);CHKERRQ(info);
484:   info=VecSetType(*B,type_name);CHKERRQ(info);
485:   info=VecSet(*B, zero);CHKERRQ(info);
486:   
487:   info=VecGetOwnershipRange(A,&low,&high);CHKERRQ(info);
488:   info=VecGetLocalSize(A,&nnloc);CHKERRQ(info);
489:   info=VecGetArray(A,&ptrA);CHKERRQ(info);
490:   info=VecGetArray(*B,&ptrB);CHKERRQ(info);
491:   info=ISGetIndices(S,&ptrS);CHKERRQ(info);
492:   for (i=0;i<nloc;i++){ ptrB[i]=ptrA[ptrS[i]-low]; }
493:   info=VecRestoreArray(A,&ptrA);CHKERRQ(info);
494:   info=VecRestoreArray(*B,&ptrB);CHKERRQ(info);
495:   info=ISRestoreIndices(S,&ptrS);CHKERRQ(info);

497:   return(0);
498: }

502: int VecStepMax(Vec X, Vec DX, PetscReal*step){
503:   int i,info;
504:   PetscInt nn;
505:   PetscScalar stepmax=1.0e300;
506:   PetscScalar *xx,*dx;
507:   double t1,t2;

510:   info = VecGetLocalSize(X,&nn);CHKERRQ(info);
511:   info = VecGetArray(X,&xx);CHKERRQ(info);
512:   info = VecGetArray(DX,&dx);CHKERRQ(info);
513:   for (i=0;i<nn;i++){
514:     if (xx[i] < 0){
515:       SETERRQ(1,"Vector must be positive");
516:     } else if (dx[i]<0){ stepmax=PetscMin(stepmax,-xx[i]/dx[i]);
517:     }
518:   }
519:   info = VecRestoreArray(X,&xx);CHKERRQ(info);
520:   info = VecRestoreArray(DX,&dx);CHKERRQ(info);
521:   t1=(double)stepmax;
522:   info = MPI_Allreduce(&t1,&t2,1,MPI_DOUBLE,MPI_MIN,((PetscObject)X)->comm);
523:   CHKERRQ(info);
524:   *step=(PetscReal)t2;
525:   return(0);
526: }

530: int MatCreateProjectionOperator(Vec A, IS S , Vec B, Mat *MM){

532:   PetscInt m,n,M,N,low,high;
533:   int info,i;
534:   const PetscInt *iptr;
535:   PetscScalar one=1.0;
536:   Mat MMM;


544:   info = VecGetSize(B,&M);CHKERRQ(info);
545:   info = VecGetLocalSize(B,&m);CHKERRQ(info);
546:   info = VecGetSize(A,&N);CHKERRQ(info);
547:   info = VecGetLocalSize(A,&n);CHKERRQ(info);

549:   info = MatCreateMPIAIJ(((PetscObject)A)->comm,m,n,M,N,1,PETSC_NULL,1,PETSC_NULL,&MMM);
550:   CHKERRQ(info);
551:   /*
552:   info = MatCreateMPIBDiag(((PetscObject)pv)->comm,m,M,N,M,1,PETSC_NULL,PETSC_NULL,&A);
553:   */
554:   
555:   info = ISGetSize(S,&n);CHKERRQ(info);
556:   if (n!=M){    
557:     SETERRQ(1,"Size of index set does not match the second vector.");
558:   }
559:   
560:   info = ISGetLocalSize(S,&n);CHKERRQ(info);
561:   if (n!=m){    
562:     SETERRQ(1,"Local size of index set does not match the second vector.");
563:   }
564:   info=VecGetOwnershipRange(B,&low,&high);CHKERRQ(info);
565:   info = ISGetIndices(S,&iptr);CHKERRQ(info);
566:   for (i=0; i<n; i++){
567:     info=MatSetValue(MMM,low+i,iptr[i],one,INSERT_VALUES);CHKERRQ(info);
568:   }
569:   info = ISRestoreIndices(S,&iptr);CHKERRQ(info);

571:   info = MatAssemblyBegin(MMM,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
572:   info = MatAssemblyEnd(MMM,MAT_FINAL_ASSEMBLY);CHKERRQ(info);



576:   *MM = MMM;

578:   return(0);  
579: }

583: int VecStepMax2(Vec X, Vec DX, Vec XL, Vec XU, PetscReal *step2){

585:   int i,info;
586:   PetscInt nn1,nn2;
587:   PetscScalar *xx,*dx,*xl,*xu;
588:   PetscScalar stepmax2=0;
589:   double t1,t2;

591:   TaoFunctionBegin;
592:   info = VecGetArray(X,&xx);CHKERRQ(info);
593:   info = VecGetArray(XL,&xl);CHKERRQ(info);
594:   info = VecGetArray(XU,&xu);CHKERRQ(info);
595:   info = VecGetArray(DX,&dx);CHKERRQ(info);
596:   info = VecGetLocalSize(X,&nn1);CHKERRQ(info);
597:   info = VecGetLocalSize(DX,&nn2);CHKERRQ(info);

599:   for (i=0;i<nn1;i++){
600:     if (dx[i] > 0){
601:       stepmax2=PetscMax(stepmax2,(xu[i]-xx[i])/dx[i]);      
602:     } else if (dx[i]<0){ 
603:       stepmax2=PetscMax(stepmax2,(xl[i]-xx[i])/dx[i]);
604:     }
605:   }
606:   info = VecRestoreArray(X,&xx);CHKERRQ(info);
607:   info = VecRestoreArray(XL,&xl);CHKERRQ(info);
608:   info = VecRestoreArray(XU,&xu);CHKERRQ(info);
609:   info = VecRestoreArray(DX,&dx);CHKERRQ(info);

611:   t1=(double)stepmax2;
612:   info = MPI_Allreduce(&t1,&t2,1,MPI_DOUBLE,MPI_MAX,((PetscObject)X)->comm);
613:   CHKERRQ(info);
614:   *step2=(PetscReal)t2;

616:   TaoFunctionReturn(0);
617: }


622: int VecBoundStepInfo(Vec X, Vec XL, Vec XU, Vec DX, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax){
623:   int info,i;
624:   PetscInt n;
625:   PetscScalar *x,*xl,*xu,*dx;
626:   double tt,t1=1.0e+20, t2=1.0e+20, t3=0;
627:   double tt1,tt2,tt3;
628:   MPI_Comm comm;
629:   
631:   info=VecGetArray(X,&x);CHKERRQ(info);
632:   info=VecGetArray(XL,&xl);CHKERRQ(info);
633:   info=VecGetArray(XU,&xu);CHKERRQ(info);
634:   info=VecGetArray(DX,&dx);CHKERRQ(info);
635:   info = VecGetLocalSize(X,&n);CHKERRQ(info);
636:   for (i=0;i<n;i++){
637:     if (dx[i]>0){
638:       tt=(xu[i]-x[i])/dx[i];
639:       t1=PetscMin(t1,tt);
640:       if (t1>0){
641:         t2=PetscMin(t2,tt);
642:       }
643:       t3=PetscMax(t3,tt);
644:     } else if (dx[i]<0){
645:       tt=(xl[i]-x[i])/dx[i];
646:       t1=PetscMin(t1,tt);
647:       if (t1>0){
648:         t2=PetscMin(t2,tt);
649:       }
650:       t3=PetscMax(t3,tt);
651:     }
652:   }
653:   info=VecRestoreArray(X,&x);CHKERRQ(info);
654:   info=VecRestoreArray(XL,&xl);CHKERRQ(info);
655:   info=VecRestoreArray(XU,&xu);CHKERRQ(info);
656:   info=VecRestoreArray(DX,&dx);CHKERRQ(info);
657:   info=PetscObjectGetComm((PetscObject)X,&comm);CHKERRQ(info);
658:   
659:   if (boundmin){ info = MPI_Allreduce(&t1,&tt1,1,MPI_DOUBLE,MPI_MIN,comm);CHKERRQ(info);}
660:   if (wolfemin){ info = MPI_Allreduce(&t2,&tt2,1,MPI_DOUBLE,MPI_MIN,comm);CHKERRQ(info);}
661:   if (boundmax) { info = MPI_Allreduce(&t3,&tt3,1,MPI_DOUBLE,MPI_MAX,comm);CHKERRQ(info);}

663:   *boundmin=(PetscReal)tt1;
664:   *wolfemin=(PetscReal)tt2;
665:   *boundmax=(PetscReal)tt3;
666:   info = PetscInfo3(X,"Step Bound Info: Closest Bound: %6.4e, Wolfe: %6.4e, Max: %6.4e \n",*boundmin,*wolfemin,*boundmax); CHKERRQ(info);
667:   return(0);  
668: }

672: /*@C
673:    VecISSetToConstant - Sets the elements of a vector, specified by an index set, to a constant

675:    Input Parameter:
676: +  S -  a PETSc IS
677: .  c - the constant
678: -  V - a Vec

680: .seealso VecSet()

682:    Level: advanced
683: @*/
684: int VecISSetToConstant(IS S, PetscScalar c, Vec V){
685:   int info,i;
686:   PetscInt nloc,low,high;
687:   const PetscInt *s;
688:   PetscScalar *v;

695:   info = VecGetOwnershipRange(V, &low, &high); CHKERRQ(info);
696:   info = ISGetLocalSize(S,&nloc);CHKERRQ(info);

698:   info = ISGetIndices(S, &s); CHKERRQ(info);
699:   info = VecGetArray(V,&v); CHKERRQ(info);
700:   for (i=0; i<nloc; i++){
701:     v[s[i]-low] = c;
702:     /*
703:     if (s[i]>=low && s[i]<high){
704:     v[s[i]-low] = c;
705:     } else {
706:       SETERRQ(1,"IS index out of range");
707:     }
708:     */
709:   }
710:   
711:   info = ISRestoreIndices(S, &s); CHKERRQ(info);
712:   info = VecRestoreArray(V,&v); CHKERRQ(info);

714:   return(0);

716: }

720: int VecPow(Vec Vec1, double p)
721: {
722:   int i,info;
723:   PetscInt n;
724:   PetscScalar *v1;


729:   info = VecGetArray(Vec1, &v1); CHKERRQ(info);
730:   info = VecGetLocalSize(Vec1, &n); CHKERRQ(info);

732:   if (1.0 == p) {
733:   }
734:   else if (-1.0 == p) {
735:     for (i = 0; i < n; ++i){
736:       v1[i] = 1.0 / v1[i];
737:     }
738:   }
739:   else if (0.0 == p) {
740:     for (i = 0; i < n; ++i){
741:       // Not-a-number left alone
742:       // Infinity set to one 
743:       if (v1[i] == v1[i]) {
744:         v1[i] = 1.0;
745:       }
746:     }
747:   }
748:   else if (0.5 == p) {
749:     for (i = 0; i < n; ++i) {
750:       if (v1[i] >= 0) {
751:         v1[i] = sqrt(v1[i]);
752:       }
753:       else {
754:         v1[i] = TAO_INFINITY;
755:       }
756:     }
757:   }
758:   else if (-0.5 == p) {
759:     for (i = 0; i < n; ++i) {
760:       if (v1[i] >= 0) {
761:         v1[i] = 1.0 / sqrt(v1[i]);
762:       }
763:       else {
764:         v1[i] = TAO_INFINITY;
765:       }
766:     }
767:   }
768:   else if (2.0 == p) {
769:     for (i = 0; i < n; ++i){
770:       v1[i] *= v1[i];
771:     }
772:   }
773:   else if (-2.0 == p) {
774:     for (i = 0; i < n; ++i){
775:       v1[i] = 1.0 / (v1[i] * v1[i]);
776:     }
777:   }
778:   else {
779:     for (i = 0; i < n; ++i) {
780:       if (v1[i] >= 0) {
781:         v1[i] = pow(v1[i], p);
782:       }
783:       else {
784:         v1[i] = TAO_INFINITY;
785:       }
786:     }
787:   }

789:   info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
790:   return(0);
791: }