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