Actual source code: pvec.c

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


  6: int PetscMapSame(PetscMap m1, PetscMap m2, PetscTruth *flag){
  7:   int i,info;
  8:   int *grange1,*grange2;
  9:   int size1,size2;
 10:   MPI_Comm comm1,comm2;

 13:   *flag=PETSC_TRUE;
 14:   info = PetscMapGetGlobalRange(m1,&grange1);CHKERRQ(info);
 15:   info = PetscMapGetGlobalRange(m2,&grange2);CHKERRQ(info);
 16:   info = PetscObjectGetComm((PetscObject)m1,&comm1);CHKERRQ(info);
 17:   info = PetscObjectGetComm((PetscObject)m2,&comm2);CHKERRQ(info);
 18:   info = MPI_Comm_size(comm1,&size1); CHKERRQ(info);
 19:   info = MPI_Comm_size(comm2,&size2); CHKERRQ(info);
 20:   if (size1!=size2){
 21:     *flag=PETSC_FALSE;
 22:   } else {
 23:     for (i=0;i<size1+1;i++){
 24:       if (grange1[i]!=grange2[i] ){
 25:         *flag=PETSC_FALSE;
 26:         break;
 27:       }
 28:     }
 29:   }
 30:   return(0);
 31: }

 35: int  VecCompare(Vec V1,Vec V2, PetscTruth *flg){
 36:   int info;
 37:   PetscMap m1,m2;
 41:   info=VecGetPetscMap(V1,&m1);CHKERRQ(info);
 42:   info=VecGetPetscMap(V2,&m2);CHKERRQ(info);
 43:   info=PetscMapSame(m1,m2,flg);CHKERRQ(info);
 44:   return(0);
 45: }

 47: /* ---------------------------------------------------------- */
 50: int VecMedian(Vec Vec1, Vec Vec2, Vec Vec3, Vec VMedian)
 51: {
 52:   int i,n,low1,low2,low3,low4,high1,high2,high3,high4,info;
 53:   PetscScalar *v1,*v2,*v3,*vmed;


 61:   if (Vec1==Vec2 || Vec1==Vec3){
 62:     info=VecCopy(Vec1,VMedian); CHKERRQ(info);
 63:     return(0);
 64:   }
 65:   if (Vec2==Vec3){
 66:     info=VecCopy(Vec2,VMedian); CHKERRQ(info);
 67:     return(0);
 68:   }


 76:   info = VecGetOwnershipRange(Vec1, &low1, &high1); CHKERRQ(info);
 77:   info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);
 78:   info = VecGetOwnershipRange(Vec3, &low3, &high3); CHKERRQ(info);
 79:   info = VecGetOwnershipRange(VMedian, &low4, &high4); CHKERRQ(info);
 80:   if ( low1!= low2 || low1!= low3 || low1!= low4 ||
 81:        high1!= high2 || high1!= high3 || high1!= high4)
 82:     SETERRQ(1,"InCompatible vector local lengths");

 84:   info = VecGetArray(Vec1,&v1); CHKERRQ(info);
 85:   info = VecGetArray(Vec2,&v2); CHKERRQ(info);
 86:   info = VecGetArray(Vec3,&v3); CHKERRQ(info);

 88:   if ( VMedian != Vec1 && VMedian != Vec2 && VMedian != Vec3){
 89:     info = VecGetArray(VMedian,&vmed); CHKERRQ(info);
 90:   } else if ( VMedian==Vec1 ){
 91:     vmed=v1;
 92:   } else if ( VMedian==Vec2 ){
 93:     vmed=v2;
 94:   } else {
 95:     vmed=v3;
 96:   }

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

100:   for (i=0;i<n;i++){
101:     vmed[i]=PetscMax(PetscMax(PetscMin(v1[i],v2[i]),PetscMin(v1[i],v3[i])),PetscMin(v2[i],v3[i]));
102:   }

104:   info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
105:   info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);
106:   info = VecRestoreArray(Vec3,&v2); CHKERRQ(info);
107: 
108:   if (VMedian!=Vec1 && VMedian != Vec2 && VMedian != Vec3){
109:     info = VecRestoreArray(VMedian,&vmed); CHKERRQ(info);
110:   }

112:   return(0);
113: }

115: /* ---------------------------------------------------------- 
118: int VecPointwiseMin(Vec Vec1, Vec Vec2, Vec VMin)
119: {
120:   int i,n,low1,low2,low3,high1,high2,high3,info;
121:   PetscScalar *v1,*v2,*vmin;


128:   if (Vec1==Vec2){
129:     info=VecCopy(Vec1,VMin); CHKERRQ(info);
130:     return(0);
131:   }


136:   info = VecGetOwnershipRange(Vec1, &low1, &high1); CHKERRQ(info);
137:   info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);
138:   info = VecGetOwnershipRange(VMin, &low3, &high3); CHKERRQ(info);

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

143:   info = VecGetArray(Vec1,&v1); CHKERRQ(info);
144:   info = VecGetArray(Vec2,&v2); CHKERRQ(info);

146:   if ( VMin != Vec1 && VMin != Vec2){
147:     info = VecGetArray(VMin,&vmin); CHKERRQ(info);
148:   } else if ( VMin==Vec1 ){
149:     vmin=v1;
150:   } else {
151:     vmin =v2;
152:   }

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

156:   for (i=0; i<n; i++){
157:     vmin[i] = PetscMin(v1[i],v2[i]);
158:   }

160:   info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
161:   info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);

163:   if (VMin!=Vec1 && VMin != Vec2){
164:     info = VecRestoreArray(VMin,&vmin); CHKERRQ(info);
165:   }

167:   return(0);
168: }
169: */
170: /* ---------------------------------------------------------- 
173: int VecPointwiseMax(Vec Vec1, Vec Vec2, Vec VMax)
174: {
175:   int i,n,low1,low2,low3,high1,high2,high3,info;
176:   PetscScalar *v1,*v2,*vmax;


180:   if (Vec1==Vec2){
181:     info=VecCopy(Vec1,VMax); CHKERRQ(info);
182:     return(0);
183:   }


192:   info = VecGetOwnershipRange(Vec1, &low1, &high1); CHKERRQ(info);
193:   info = VecGetOwnershipRange(Vec2, &low2, &high2); CHKERRQ(info);
194:   info = VecGetOwnershipRange(VMax, &low3, &high3); CHKERRQ(info);

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

199:   info = VecGetArray(Vec1,&v1); CHKERRQ(info);
200:   info = VecGetArray(Vec2,&v2); CHKERRQ(info);

202:   if ( VMax != Vec1 && VMax != Vec2){
203:     info = VecGetArray(VMax,&vmax); CHKERRQ(info);
204:   } else if ( VMax==Vec1 ){ vmax=v1;
205:   } else { vmax =v2;
206:   }

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

210:   for (i=0; i<n; i++){
211:     vmax[i] = PetscMax(v1[i],v2[i]);
212:   }

214:   info = VecRestoreArray(Vec1,&v1); CHKERRQ(info);
215:   info = VecRestoreArray(Vec2,&v2); CHKERRQ(info);

217:   if ( VMax != Vec1 && VMax != Vec2){
218:     info = VecRestoreArray(VMax,&vmax); CHKERRQ(info);
219:   }

221:   return(0);
222: }*/

226: inline static PetscScalar Fischer(PetscScalar a, PetscScalar b)
227: {

229: #ifdef TODD

231:    if (PetscAbsScalar(a) > PetscAbsScalar(b)) {
232:      return sqrt(a*a + b*b) - a - b;
233:    }
234:    return sqrt(a*a + b*b) - b - a;

236: #else

238:    // Method suggested by Bob Vanderbei

240:    if (a + b <= 0) {
241:      return sqrt(a*a + b*b) - (a + b);
242:    }
243:    return -2.0*a*b / (sqrt(a*a + b*b) + (a + b));

245: #endif

247: }

251: int VecFischer(Vec X, Vec F, Vec L, Vec U, Vec FF)
252: {
253:   PetscScalar *x, *f, *l, *u, *ff;
254:   PetscScalar xval, fval, lval, uval;
255:   int info, i, n;
256:   int low[5], high[5];


265:   info = VecGetOwnershipRange(X, low, high); CHKERRQ(info);
266:   info = VecGetOwnershipRange(F, low + 1, high + 1); CHKERRQ(info);
267:   info = VecGetOwnershipRange(L, low + 2, high + 2); CHKERRQ(info);
268:   info = VecGetOwnershipRange(U, low + 3, high + 3); CHKERRQ(info);
269:   info = VecGetOwnershipRange(FF, low + 4, high + 4); CHKERRQ(info);

271:   for (i = 1; i < 4; ++i) {
272:     if (low[0] != low[i] || high[0] != high[i])
273:       SETERRQ(1,"Vectors must be identically loaded over processors");
274:   }

276:   info = VecGetArray(X, &x); CHKERRQ(info);
277:   info = VecGetArray(F, &f); CHKERRQ(info);
278:   info = VecGetArray(L, &l); CHKERRQ(info);
279:   info = VecGetArray(U, &u); CHKERRQ(info);
280:   info = VecGetArray(FF, &ff); CHKERRQ(info);

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

284:   for (i = 0; i < n; ++i) {
285:     xval = (*x++); fval = (*f++);
286:     lval = (*l++); uval = (*u++);

288:     if ((lval <= -TAO_INFINITY) && (uval >= TAO_INFINITY)) {
289:       (*ff++) = -fval;
290:     }
291:     else if (lval <= -TAO_INFINITY) {
292:       (*ff++) = -Fischer(uval - xval, -fval);
293:     }
294:     else if (uval >=  TAO_INFINITY) {
295:       (*ff++) =  Fischer(xval - lval,  fval);
296:     }
297:     else if (lval == uval) {
298:       (*ff++) = lval - xval;
299:     }
300:     else {
301:       fval    =  Fischer(uval - xval, -fval);
302:       (*ff++) =  Fischer(xval - lval,  fval);
303:     }

305:   }
306:   x -= n; f -= n; l -=n; u -= n; ff -= n;
307: 
308:   info = VecRestoreArray(X, &x); CHKERRQ(info);
309:   info = VecRestoreArray(F, &f); CHKERRQ(info);
310:   info = VecRestoreArray(L, &l); CHKERRQ(info);
311:   info = VecRestoreArray(U, &u); CHKERRQ(info);
312:   info = VecRestoreArray(FF, &ff); CHKERRQ(info);

314:   return(0);
315: }

319: inline static PetscScalar SFischer(PetscScalar a, PetscScalar b, PetscScalar c)
320: {

322: #ifdef TODD

324:    if (PetscAbsScalar(a) > PetscAbsScalar(b)) {
325:      return sqrt(a*a + b*b + 2.0*c*c) - a - b;
326:    }
327:    return sqrt(a*a + b*b + 2.0*c*c) - b - a;

329: #else

331:    // Method suggested by Bob Vanderbei

333:    if (a + b <= 0) {
334:      return sqrt(a*a + b*b + 2.0*c*c) - (a + b);
335:    }
336:    return 2.0*(c*c - a*b) / (sqrt(a*a + b*b + 2.0*c*c) + (a + b));

338: #endif

340: }

344: int VecSFischer(Vec X, Vec F, Vec L, Vec U, double mu, Vec FF)
345: {
346:   PetscScalar *x, *f, *l, *u, *ff;
347:   PetscScalar xval, fval, lval, uval;
348:   int info, i, n;
349:   int low[5], high[5];


358:   info = VecGetOwnershipRange(X, low, high); CHKERRQ(info);
359:   info = VecGetOwnershipRange(F, low + 1, high + 1); CHKERRQ(info);
360:   info = VecGetOwnershipRange(L, low + 2, high + 2); CHKERRQ(info);
361:   info = VecGetOwnershipRange(U, low + 3, high + 3); CHKERRQ(info);
362:   info = VecGetOwnershipRange(FF, low + 4, high + 4); CHKERRQ(info);

364:   for (i = 1; i < 4; ++i) {
365:     if (low[0] != low[i] || high[0] != high[i])
366:       SETERRQ(1,"Vectors must be identically loaded over processors");
367:   }

369:   info = VecGetArray(X, &x); CHKERRQ(info);
370:   info = VecGetArray(F, &f); CHKERRQ(info);
371:   info = VecGetArray(L, &l); CHKERRQ(info);
372:   info = VecGetArray(U, &u); CHKERRQ(info);
373:   info = VecGetArray(FF, &ff); CHKERRQ(info);

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

377:   for (i = 0; i < n; ++i) {
378:     xval = (*x++); fval = (*f++);
379:     lval = (*l++); uval = (*u++);

381:     if ((lval <= -TAO_INFINITY) && (uval >= TAO_INFINITY)) {
382:       (*ff++) = -fval - mu*xval;
383:     }
384:     else if (lval <= -TAO_INFINITY) {
385:       (*ff++) = -SFischer(uval - xval, -fval, mu);
386:     }
387:     else if (uval >=  TAO_INFINITY) {
388:       (*ff++) =  SFischer(xval - lval,  fval, mu);
389:     }
390:     else if (lval == uval) {
391:       (*ff++) = lval - xval;
392:     }
393:     else {
394:       fval    =  SFischer(uval - xval, -fval, mu);
395:       (*ff++) =  SFischer(xval - lval,  fval, mu);
396:     }
397:   }
398:   x -= n; f -= n; l -=n; u -= n; ff -= n;

400:   info = VecRestoreArray(X, &x); CHKERRQ(info);
401:   info = VecRestoreArray(F, &f); CHKERRQ(info);
402:   info = VecRestoreArray(L, &l); CHKERRQ(info);
403:   info = VecRestoreArray(U, &u); CHKERRQ(info);
404:   info = VecRestoreArray(FF, &ff); CHKERRQ(info);
405:   return(0);
406: }

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

412:   int i,n,info;
413:   PetscScalar *xptr,*xlptr,*xuptr,*gptr,*gpptr;
414:   PetscScalar xval,gpval;

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

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

421:   info=VecGetArray(X,&xptr); CHKERRQ(info);
422:   info=VecGetArray(XL,&xlptr); CHKERRQ(info);
423:   info=VecGetArray(XU,&xuptr); CHKERRQ(info);
424:   info=VecGetArray(G,&gptr); CHKERRQ(info);
425:   if (G!=GP){
426:     info=VecGetArray(GP,&gpptr); CHKERRQ(info);
427:   } else { gpptr=gptr; }

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

432:     if (gpval>0 && xval<=xlptr[i]){
433:       gpval = 0;
434:     } else if (gpval<0 && xval>=xuptr[i]){
435:       gpval = 0;
436:     }
437:     gpptr[i] = gpval;

439:     /*
440:     if (xlptr[i] >= xuptr[i]){
441:       gpptr[i]=0.0;
442:     } else if (xptr[i] <= xlptr[i]+eps){
443:       gpptr[i] = PetscMin(gptr[i],zero);
444:     } else if (xptr[i] >= xuptr[i]-eps){
445:       gpptr[i] = PetscMax(gptr[i],zero);
446:     } else {
447:       gpptr[i] = gptr[i];
448:     }
449:     */
450:   }

452:   info=VecRestoreArray(X,&xptr); CHKERRQ(info);
453:   info=VecRestoreArray(XL,&xlptr); CHKERRQ(info);
454:   info=VecRestoreArray(XU,&xuptr); CHKERRQ(info);
455:   info=VecRestoreArray(G,&gptr); CHKERRQ(info);
456:   if (G!=GP){
457:     info=VecRestoreArray(GP,&gpptr); CHKERRQ(info);
458:   }
459:   return(0);
460: }

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

466:   int i,info,in1,in2,xlow,xhigh,ylow,yhigh,*yis;
467:   PetscScalar *x,*y;

470:   info=ISGetLocalSize(YIS,&in1); CHKERRQ(info);
471:   info=VecGetLocalSize(X,&in2); CHKERRQ(info);

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

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

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

481:   info=VecGetArray(X,&x); CHKERRQ(info);
482:   if (X != Y){
483:     info=VecGetArray(Y,&y); CHKERRQ(info);
484:   } else {
485:     y=x;
486:   }

488:   for (i=0; i<in1; i++){
489:     if (yis[i]>=ylow && yis[i]<yhigh && i+xlow < xhigh){
490:       y[yis[i]-ylow] = y[yis[i]-ylow] + (alpha)*x[i];
491:     } else {
492:       SETERRQ(1,"IS index out of range");
493:     }
494:   }
495: 
496:   info=ISRestoreIndices(YIS, &yis); CHKERRQ(info);

498:   info=VecRestoreArray(X,&x); CHKERRQ(info);
499:   if (X != Y){
500:     info=VecRestoreArray(Y,&y); CHKERRQ(info);
501:   }

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

505:   return(0);
506: }

510: int VecCreateSubVec(Vec A,IS S,Vec *B)
511: {
512:   int     info,nloc,n;
513:   int i,low,high,nnloc;
514:   int *ptrS;
515:   PetscScalar  zero=0.0;
516:   PetscScalar *ptrB,*ptrA;
517:   VecType type_name;
518:   MPI_Comm comm;

523:   info=ISGetLocalSize(S,&nloc);CHKERRQ(info);
524:   info=ISGetSize(S,&n);CHKERRQ(info);
525:   info=PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(info);
526:   info=VecCreate(comm,B);CHKERRQ(info);
527:   info=VecSetSizes(*B,nloc,n);CHKERRQ(info);
528:   info=VecGetType(A,&type_name);CHKERRQ(info);
529:   info=VecSetType(*B,type_name);CHKERRQ(info);
530:   info=VecSet(*B, zero);CHKERRQ(info);
531: 
532:   info=VecGetOwnershipRange(A,&low,&high);CHKERRQ(info);
533:   info=VecGetLocalSize(A,&nnloc);CHKERRQ(info);
534:   info=VecGetArray(A,&ptrA);CHKERRQ(info);
535:   info=VecGetArray(*B,&ptrB);CHKERRQ(info);
536:   info=ISGetIndices(S,&ptrS);CHKERRQ(info);
537:   for (i=0;i<nloc;i++){ ptrB[i]=ptrA[ptrS[i]-low]; }
538:   info=VecRestoreArray(A,&ptrA);CHKERRQ(info);
539:   info=VecRestoreArray(*B,&ptrB);CHKERRQ(info);
540:   info=ISRestoreIndices(S,&ptrS);CHKERRQ(info);

542:   return(0);
543: }

547: int VecStepMax(Vec X, Vec DX, PetscReal*step){
548:   int i,info,nn;
549:   PetscScalar stepmax=1.0e300;
550:   PetscScalar *xx,*dx;
551:   double t1,t2;

554:   info = VecGetLocalSize(X,&nn);CHKERRQ(info);
555:   info = VecGetArray(X,&xx);CHKERRQ(info);
556:   info = VecGetArray(DX,&dx);CHKERRQ(info);
557:   for (i=0;i<nn;i++){
558:     if (xx[i] < 0){
559:       SETERRQ(1,"Vector must be positive");
560:     } else if (dx[i]<0){ stepmax=PetscMin(stepmax,-xx[i]/dx[i]);
561:     }
562:   }
563:   info = VecRestoreArray(X,&xx);CHKERRQ(info);
564:   info = VecRestoreArray(DX,&dx);CHKERRQ(info);
565:   t1=(double)stepmax;
566:   info = MPI_Allreduce(&t1,&t2,1,MPI_DOUBLE,MPI_MIN,X->comm);
567:   CHKERRQ(info);
568:   *step=(PetscReal)t2;
569:   return(0);
570: }

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

576:   int m,n,M,N,low,high;
577:   int info,i,*iptr;
578:   PetscScalar one=1.0;
579:   Mat MMM;


587:   info = VecGetSize(B,&M);CHKERRQ(info);
588:   info = VecGetLocalSize(B,&m);CHKERRQ(info);
589:   info = VecGetSize(A,&N);CHKERRQ(info);
590:   info = VecGetLocalSize(A,&n);CHKERRQ(info);

592:   info = MatCreateMPIAIJ(A->comm,m,n,M,N,1,PETSC_NULL,1,PETSC_NULL,&MMM);
593:   CHKERRQ(info);
594:   /*
595:   info = MatCreateMPIBDiag(pv->comm,m,M,N,M,1,PETSC_NULL,PETSC_NULL,&A);
596:   */
597: 
598:   info = ISGetSize(S,&n);CHKERRQ(info);
599:   if (n!=M){
600:     SETERRQ(1,"Size of index set does not match the second vector.");
601:   }
602: 
603:   info = ISGetLocalSize(S,&n);CHKERRQ(info);
604:   if (n!=m){
605:     SETERRQ(1,"Local size of index set does not match the second vector.");
606:   }
607:   info=VecGetOwnershipRange(B,&low,&high);CHKERRQ(info);
608:   info = ISGetIndices(S,&iptr);CHKERRQ(info);
609:   for (i=0; i<n; i++){
610:     info=MatSetValue(MMM,low+i,iptr[i],one,INSERT_VALUES);CHKERRQ(info);
611:   }
612:   info = ISRestoreIndices(S,&iptr);CHKERRQ(info);

614:   info = MatAssemblyBegin(MMM,MAT_FINAL_ASSEMBLY);CHKERRQ(info);
615:   info = MatAssemblyEnd(MMM,MAT_FINAL_ASSEMBLY);CHKERRQ(info);

617:   info = MatCompress(MMM);CHKERRQ(info);

619:   *MM = MMM;

621:   return(0);
622: }

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

628:   int i,nn1,nn2,info;
629:   PetscScalar *xx,*dx,*xl,*xu;
630:   PetscScalar stepmax2=0;
631:   double t1,t2;

633:   TaoFunctionBegin;
634:   info = VecGetArray(X,&xx);CHKERRQ(info);
635:   info = VecGetArray(XL,&xl);CHKERRQ(info);
636:   info = VecGetArray(XU,&xu);CHKERRQ(info);
637:   info = VecGetArray(DX,&dx);CHKERRQ(info);
638:   info = VecGetLocalSize(X,&nn1);CHKERRQ(info);
639:   info = VecGetLocalSize(DX,&nn2);CHKERRQ(info);

641:   for (i=0;i<nn1;i++){
642:     if (dx[i] > 0){
643:       stepmax2=PetscMax(stepmax2,(xu[i]-xx[i])/dx[i]);
644:     } else if (dx[i]<0){
645:       stepmax2=PetscMax(stepmax2,(xl[i]-xx[i])/dx[i]);
646:     }
647:   }
648:   info = VecRestoreArray(X,&xx);CHKERRQ(info);
649:   info = VecRestoreArray(XL,&xl);CHKERRQ(info);
650:   info = VecRestoreArray(XU,&xu);CHKERRQ(info);
651:   info = VecRestoreArray(DX,&dx);CHKERRQ(info);

653:   t1=(double)stepmax2;
654:   info = MPI_Allreduce(&t1,&t2,1,MPI_DOUBLE,MPI_MAX,X->comm);
655:   CHKERRQ(info);
656:   *step2=(PetscReal)t2;

658:   TaoFunctionReturn(0);
659: }


664: int VecBoundStepInfo(Vec X, Vec XL, Vec XU, Vec DX, PetscReal *boundmin, PetscReal *wolfemin, PetscReal *boundmax){
665:   int info,i,n;
666:   PetscScalar *x,*xl,*xu,*dx;
667:   double tt,t1=1.0e+20, t2=1.0e+20, t3=0;
668:   double tt1,tt2,tt3;
669:   MPI_Comm comm;
670: 
672:   info=VecGetArray(X,&x);CHKERRQ(info);
673:   info=VecGetArray(XL,&xl);CHKERRQ(info);
674:   info=VecGetArray(XU,&xu);CHKERRQ(info);
675:   info=VecGetArray(DX,&dx);CHKERRQ(info);
676:   info = VecGetLocalSize(X,&n);CHKERRQ(info);
677:   for (i=0;i<n;i++){
678:     if (dx[i]>0){
679:       tt=(xu[i]-x[i])/dx[i];
680:       t1=PetscMin(t1,tt);
681:       if (t1>0){
682:         t2=PetscMin(t2,tt);
683:       }
684:       t3=PetscMax(t3,tt);
685:     } else if (dx[i]<0){
686:       tt=(xl[i]-x[i])/dx[i];
687:       t1=PetscMin(t1,tt);
688:       if (t1>0){
689:         t2=PetscMin(t2,tt);
690:       }
691:       t3=PetscMax(t3,tt);
692:     }
693:   }
694:   info=VecRestoreArray(X,&x);CHKERRQ(info);
695:   info=VecRestoreArray(XL,&xl);CHKERRQ(info);
696:   info=VecRestoreArray(XU,&xu);CHKERRQ(info);
697:   info=VecRestoreArray(DX,&dx);CHKERRQ(info);
698:   info=PetscObjectGetComm((PetscObject)X,&comm);CHKERRQ(info);
699: 
700:   if (boundmin){ info = MPI_Allreduce(&t1,&tt1,1,MPI_DOUBLE,MPI_MIN,comm);CHKERRQ(info);}
701:   if (wolfemin){ info = MPI_Allreduce(&t2,&tt2,1,MPI_DOUBLE,MPI_MIN,comm);CHKERRQ(info);}
702:   if (boundmax) { info = MPI_Allreduce(&t3,&tt3,1,MPI_DOUBLE,MPI_MAX,comm);CHKERRQ(info);}

704:   *boundmin=(PetscReal)tt1;
705:   *wolfemin=(PetscReal)tt2;
706:   *boundmax=(PetscReal)tt3;
707:   info = PetscLogInfo((X,"Step Bound Info: Closest Bound: %6.4e, Wolfe: %6.4e, Max: %6.4e \n",*boundmin,*wolfemin,*boundmax)); CHKERRQ(info);
708:   return(0);
709: }

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

716:    Input Parameter:
717: +  S -  a PETSc IS
718: .  c - the constant
719: -  V - a Vec

721: .seealso VecSet()

723:    Level: advanced
724: @*/
725: int VecISSetToConstant(IS S, PetscScalar c, Vec V){
726:   int info,i,nloc,low,high,*s;
727:   PetscScalar *v;

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

737:   info = ISGetIndices(S, &s); CHKERRQ(info);
738:   info = VecGetArray(V,&v); CHKERRQ(info);
739:   for (i=0; i<nloc; i++){
740:     v[s[i]-low] = c;
741:     /*
742:     if (s[i]>=low && s[i]<high){
743:     v[s[i]-low] = c;
744:     } else {
745:       SETERRQ(1,"IS index out of range");
746:     }
747:     */
748:   }
749: 
750:   info = ISRestoreIndices(S, &s); CHKERRQ(info);
751:   info = VecRestoreArray(V,&v); CHKERRQ(info);

753:   return(0);

755: }