Actual source code: isutil.c

      1: #include <petsctao.h>
      2: #include <petsc/private/vecimpl.h>
      3: #include <petsc/private/taoimpl.h>
      4: #include <../src/tao/matrix/submatfree.h>

      6: /*@C
      7:   TaoVecGetSubVec - Gets a subvector using the IS

      9:   Input Parameters:
     10: + vfull - the full matrix
     11: . is - the index set for the subvector
     12: . reduced_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK,  TAO_SUBSET_MATRIXFREE)
     13: - maskvalue - the value to set the unused vector elements to (for TAO_SUBSET_MASK or TAO_SUBSET_MATRIXFREE)

     15:   Output Parameters:
     16: . vreduced - the subvector

     18:   Notes:
     19:   maskvalue should usually be 0.0, unless a pointwise divide will be used.

     21:   Level: developer
     22: @*/
     23: PetscErrorCode TaoVecGetSubVec(Vec vfull, IS is, TaoSubsetType reduced_type, PetscReal maskvalue, Vec *vreduced)
     24: {
     26:   PetscInt       nfull,nreduced,nreduced_local,rlow,rhigh,flow,fhigh;
     27:   PetscInt       i,nlocal;
     28:   PetscReal      *fv,*rv;
     29:   const PetscInt *s;
     30:   IS             ident;
     31:   VecType        vtype;
     32:   VecScatter     scatter;
     33:   MPI_Comm       comm;

     39:   VecGetSize(vfull, &nfull);
     40:   ISGetSize(is, &nreduced);

     42:   if (nreduced == nfull) {
     43:     VecDestroy(vreduced);
     44:     VecDuplicate(vfull,vreduced);
     45:     VecCopy(vfull,*vreduced);
     46:   } else {
     47:     switch (reduced_type) {
     48:     case TAO_SUBSET_SUBVEC:
     49:       VecGetType(vfull,&vtype);
     50:       VecGetOwnershipRange(vfull,&flow,&fhigh);
     51:       ISGetLocalSize(is,&nreduced_local);
     52:       PetscObjectGetComm((PetscObject)vfull,&comm);
     53:       if (*vreduced) {
     54:         VecDestroy(vreduced);
     55:       }
     56:       VecCreate(comm,vreduced);
     57:       VecSetType(*vreduced,vtype);

     59:       VecSetSizes(*vreduced,nreduced_local,nreduced);
     60:       VecGetOwnershipRange(*vreduced,&rlow,&rhigh);
     61:       ISCreateStride(comm,nreduced_local,rlow,1,&ident);
     62:       VecScatterCreate(vfull,is,*vreduced,ident,&scatter);
     63:       VecScatterBegin(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);
     64:       VecScatterEnd(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);
     65:       VecScatterDestroy(&scatter);
     66:       ISDestroy(&ident);
     67:       break;

     69:     case TAO_SUBSET_MASK:
     70:     case TAO_SUBSET_MATRIXFREE:
     71:       /* vr[i] = vf[i]   if i in is
     72:        vr[i] = 0       otherwise */
     73:       if (!*vreduced) {
     74:         VecDuplicate(vfull,vreduced);
     75:       }

     77:       VecSet(*vreduced,maskvalue);
     78:       ISGetLocalSize(is,&nlocal);
     79:       VecGetOwnershipRange(vfull,&flow,&fhigh);
     80:       VecGetArray(vfull,&fv);
     81:       VecGetArray(*vreduced,&rv);
     82:       ISGetIndices(is,&s);
     83:       if (nlocal > (fhigh-flow)) SETERRQ2(PETSC_COMM_SELF,1,"IS local size %D > Vec local size %D",nlocal,fhigh-flow);
     84:       for (i=0;i<nlocal;++i) {
     85:         rv[s[i]-flow] = fv[s[i]-flow];
     86:       }
     87:       ISRestoreIndices(is,&s);
     88:       VecRestoreArray(vfull,&fv);
     89:       VecRestoreArray(*vreduced,&rv);
     90:       break;
     91:     }
     92:   }
     93:   return(0);
     94: }

     96: /*@C
     97:   TaoMatGetSubMat - Gets a submatrix using the IS

     99:   Input Parameters:
    100: + M - the full matrix (n x n)
    101: . is - the index set for the submatrix (both row and column index sets need to be the same)
    102: . v1 - work vector of dimension n, needed for TAO_SUBSET_MASK option
    103: - subset_type <TAO_SUBSET_SUBVEC,TAO_SUBSET_MASK,TAO_SUBSET_MATRIXFREE> - the method TAO is using for subsetting

    105:   Output Parameters:
    106: . Msub - the submatrix

    108:   Level: developer
    109: @*/
    110: PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
    111: {
    113:   IS             iscomp;
    114:   PetscBool      flg = PETSC_TRUE;

    119:   MatDestroy(Msub);
    120:   switch (subset_type) {
    121:   case TAO_SUBSET_SUBVEC:
    122:     MatCreateSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);
    123:     break;

    125:   case TAO_SUBSET_MASK:
    126:     /* Get Reduced Hessian
    127:      Msub[i,j] = M[i,j] if i,j in Free_Local or i==j
    128:      Msub[i,j] = 0      if i!=j and i or j not in Free_Local
    129:      */
    130:     PetscObjectOptionsBegin((PetscObject)M);
    131:     PetscOptionsBool("-overwrite_hessian","modify the existing hessian matrix when computing submatrices","TaoSubsetType",flg,&flg,NULL);
    132:     PetscOptionsEnd();
    133:     if (flg) {
    134:       MatDuplicate(M, MAT_COPY_VALUES, Msub);
    135:     } else {
    136:       /* Act on hessian directly (default) */
    137:       PetscObjectReference((PetscObject)M);
    138:       *Msub = M;
    139:     }
    140:     /* Save the diagonal to temporary vector */
    141:     MatGetDiagonal(*Msub,v1);

    143:     /* Zero out rows and columns */
    144:     ISComplementVec(is,v1,&iscomp);

    146:     /* Use v1 instead of 0 here because of PETSc bug */
    147:     MatZeroRowsColumnsIS(*Msub,iscomp,1.0,v1,v1);

    149:     ISDestroy(&iscomp);
    150:     break;
    151:   case TAO_SUBSET_MATRIXFREE:
    152:     ISComplementVec(is,v1,&iscomp);
    153:     MatCreateSubMatrixFree(M,iscomp,iscomp,Msub);
    154:     ISDestroy(&iscomp);
    155:     break;
    156:   }
    157:   return(0);
    158: }

    160: /*@C
    161:   TaoEstimateActiveBounds - Generates index sets for variables at the lower and upper
    162:   bounds, as well as fixed variables where lower and upper bounds equal each other.

    164:   Input Parameters:
    165: + X - solution vector
    166: . XL - lower bound vector
    167: . XU - upper bound vector
    168: . G - unprojected gradient
    169: . S - step direction with which the active bounds will be estimated
    170: . W - work vector of type and size of X
    171: - steplen - the step length at which the active bounds will be estimated (needs to be conservative)

    173:   Output Parameters:
    174: + bound_tol - tolerance for for the bound estimation
    175: . active_lower - index set for active variables at the lower bound
    176: . active_upper - index set for active variables at the upper bound
    177: . active_fixed - index set for fixed variables
    178: . active - index set for all active variables
    179: - inactive - complementary index set for inactive variables

    181:   Notes:
    182:   This estimation is based on Bertsekas' method, with a built in diagonal scaling value of 1.0e-3.

    184:   Level: developer
    185: @*/
    186: PetscErrorCode TaoEstimateActiveBounds(Vec X, Vec XL, Vec XU, Vec G, Vec S, Vec W, PetscReal steplen, PetscReal *bound_tol,
    187:                                        IS *active_lower, IS *active_upper, IS *active_fixed, IS *active, IS *inactive)
    188: {
    189:   PetscErrorCode               ierr;
    190:   PetscReal                    wnorm;
    191:   PetscReal                    zero = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0/3.0);
    192:   PetscInt                     i, n_isl=0, n_isu=0, n_isf=0, n_isa=0, n_isi=0;
    193:   PetscInt                     N_isl, N_isu, N_isf, N_isa, N_isi;
    194:   PetscInt                     n, low, high, nDiff;
    195:   PetscInt                     *isl=NULL, *isu=NULL, *isf=NULL, *isa=NULL, *isi=NULL;
    196:   const PetscScalar            *xl, *xu, *x, *g;
    197:   MPI_Comm                     comm = PetscObjectComm((PetscObject)X);

    223:   VecCheckSameSize(X,1,XL,2);
    224:   VecCheckSameSize(X,1,XU,3);
    225:   VecCheckSameSize(X,1,G,4);
    226:   VecCheckSameSize(X,1,S,5);
    227:   VecCheckSameSize(X,1,W,6);

    229:   /* Update the tolerance for bound detection (this is based on Bertsekas' method) */
    230:   VecCopy(X, W);
    231:   VecAXPBY(W, steplen, 1.0, S);
    232:   TaoBoundSolution(W, XL, XU, 0.0, &nDiff, W);
    233:   VecAXPBY(W, 1.0, -1.0, X);
    234:   VecNorm(W, NORM_2, &wnorm);
    235:   *bound_tol = PetscMin(*bound_tol, wnorm);

    237:   VecGetOwnershipRange(X, &low, &high);
    238:   VecGetLocalSize(X, &n);
    239:   if (n>0){
    240:     VecGetArrayRead(X, &x);
    241:     VecGetArrayRead(XL, &xl);
    242:     VecGetArrayRead(XU, &xu);
    243:     VecGetArrayRead(G, &g);

    245:     /* Loop over variables and categorize the indexes */
    246:     PetscMalloc1(n, &isl);
    247:     PetscMalloc1(n, &isu);
    248:     PetscMalloc1(n, &isf);
    249:     PetscMalloc1(n, &isa);
    250:     PetscMalloc1(n, &isi);
    251:     for (i=0; i<n; ++i) {
    252:       if (xl[i] == xu[i]) {
    253:         /* Fixed variables */
    254:         isf[n_isf]=low+i; ++n_isf;
    255:         isa[n_isa]=low+i; ++n_isa;
    256:       } else if ((xl[i] > PETSC_NINFINITY) && (x[i] <= xl[i] + *bound_tol) && (g[i] > zero)) {
    257:         /* Lower bounded variables */
    258:         isl[n_isl]=low+i; ++n_isl;
    259:         isa[n_isa]=low+i; ++n_isa;
    260:       } else if ((xu[i] < PETSC_INFINITY) && (x[i] >= xu[i] - *bound_tol) && (g[i] < zero)) {
    261:         /* Upper bounded variables */
    262:         isu[n_isu]=low+i; ++n_isu;
    263:         isa[n_isa]=low+i; ++n_isa;
    264:       } else {
    265:         /* Inactive variables */
    266:         isi[n_isi]=low+i; ++n_isi;
    267:       }
    268:     }

    270:     VecRestoreArrayRead(X, &x);
    271:     VecRestoreArrayRead(XL, &xl);
    272:     VecRestoreArrayRead(XU, &xu);
    273:     VecRestoreArrayRead(G, &g);
    274:   }

    276:   /* Clear all index sets */
    277:   ISDestroy(active_lower);
    278:   ISDestroy(active_upper);
    279:   ISDestroy(active_fixed);
    280:   ISDestroy(active);
    281:   ISDestroy(inactive);

    283:   /* Collect global sizes */
    284:   MPIU_Allreduce(&n_isl, &N_isl, 1, MPIU_INT, MPI_SUM, comm);
    285:   MPIU_Allreduce(&n_isu, &N_isu, 1, MPIU_INT, MPI_SUM, comm);
    286:   MPIU_Allreduce(&n_isf, &N_isf, 1, MPIU_INT, MPI_SUM, comm);
    287:   MPIU_Allreduce(&n_isa, &N_isa, 1, MPIU_INT, MPI_SUM, comm);
    288:   MPIU_Allreduce(&n_isi, &N_isi, 1, MPIU_INT, MPI_SUM, comm);

    290:   /* Create index set for lower bounded variables */
    291:   if (N_isl > 0) {
    292:     ISCreateGeneral(comm, n_isl, isl, PETSC_OWN_POINTER, active_lower);
    293:   } else {
    294:     PetscFree(isl);
    295:   }
    296:   /* Create index set for upper bounded variables */
    297:   if (N_isu > 0) {
    298:     ISCreateGeneral(comm, n_isu, isu, PETSC_OWN_POINTER, active_upper);
    299:   } else {
    300:     PetscFree(isu);
    301:   }
    302:   /* Create index set for fixed variables */
    303:   if (N_isf > 0) {
    304:     ISCreateGeneral(comm, n_isf, isf, PETSC_OWN_POINTER, active_fixed);
    305:   } else {
    306:     PetscFree(isf);
    307:   }
    308:   /* Create index set for all actively bounded variables */
    309:   if (N_isa > 0) {
    310:     ISCreateGeneral(comm, n_isa, isa, PETSC_OWN_POINTER, active);
    311:   } else {
    312:     PetscFree(isa);
    313:   }
    314:   /* Create index set for all inactive variables */
    315:   if (N_isi > 0) {
    316:     ISCreateGeneral(comm, n_isi, isi, PETSC_OWN_POINTER, inactive);
    317:   } else {
    318:     PetscFree(isi);
    319:   }

    321:   /* Clean up and exit */
    322:   return(0);
    323: }

    325: /*@C
    326:   TaoBoundStep - Ensures the correct zero or adjusted step direction
    327:   values for active variables.

    329:   Input Parameters:
    330: + X - solution vector
    331: . XL - lower bound vector
    332: . XU - upper bound vector
    333: . active_lower - index set for lower bounded active variables
    334: . active_upper - index set for lower bounded active variables
    335: . active_fixed - index set for fixed active variables
    336: - scale - amplification factor for the step that needs to be taken on actively bounded variables

    338:   Output Parameters:
    339: . S - step direction to be modified

    341:   Level: developer
    342: @*/
    343: PetscErrorCode TaoBoundStep(Vec X, Vec XL, Vec XU, IS active_lower, IS active_upper, IS active_fixed, PetscReal scale, Vec S)
    344: {
    345:   PetscErrorCode               ierr;

    347:   Vec                          step_lower, step_upper, step_fixed;
    348:   Vec                          x_lower, x_upper;
    349:   Vec                          bound_lower, bound_upper;

    352:   /* Adjust step for variables at the estimated lower bound */
    353:   if (active_lower) {
    354:     VecGetSubVector(S, active_lower, &step_lower);
    355:     VecGetSubVector(X, active_lower, &x_lower);
    356:     VecGetSubVector(XL, active_lower, &bound_lower);
    357:     VecCopy(bound_lower, step_lower);
    358:     VecAXPY(step_lower, -1.0, x_lower);
    359:     VecScale(step_lower, scale);
    360:     VecRestoreSubVector(S, active_lower, &step_lower);
    361:     VecRestoreSubVector(X, active_lower, &x_lower);
    362:     VecRestoreSubVector(XL, active_lower, &bound_lower);
    363:   }

    365:   /* Adjust step for the variables at the estimated upper bound */
    366:   if (active_upper) {
    367:     VecGetSubVector(S, active_upper, &step_upper);
    368:     VecGetSubVector(X, active_upper, &x_upper);
    369:     VecGetSubVector(XU, active_upper, &bound_upper);
    370:     VecCopy(bound_upper, step_upper);
    371:     VecAXPY(step_upper, -1.0, x_upper);
    372:     VecScale(step_upper, scale);
    373:     VecRestoreSubVector(S, active_upper, &step_upper);
    374:     VecRestoreSubVector(X, active_upper, &x_upper);
    375:     VecRestoreSubVector(XU, active_upper, &bound_upper);
    376:   }

    378:   /* Zero out step for fixed variables */
    379:   if (active_fixed) {
    380:     VecGetSubVector(S, active_fixed, &step_fixed);
    381:     VecSet(step_fixed, 0.0);
    382:     VecRestoreSubVector(S, active_fixed, &step_fixed);
    383:   }
    384:   return(0);
    385: }

    387: /*@C
    388:   TaoBoundSolution - Ensures that the solution vector is snapped into the bounds within a given tolerance.

    390:   Collective on Vec

    392:   Input Parameters:
    393: + X - solution vector
    394: . XL - lower bound vector
    395: . XU - upper bound vector
    396: - bound_tol - absolute tolerance in enforcing the bound

    398:   Output Parameters:
    399: + nDiff - total number of vector entries that have been bounded
    400: - Xout - modified solution vector satisfying bounds to bound_tol

    402:   Level: developer

    404: .seealso: TAOBNCG, TAOBNTL, TAOBNTR
    405: @*/
    406: PetscErrorCode TaoBoundSolution(Vec X, Vec XL, Vec XU, PetscReal bound_tol, PetscInt *nDiff, Vec Xout)
    407: {
    408:   PetscErrorCode    ierr;
    409:   PetscInt          i,n,low,high,nDiff_loc=0;
    410:   PetscScalar       *xout;
    411:   const PetscScalar *x,*xl,*xu;

    429:   VecCheckSameSize(X,1,XL,2);
    430:   VecCheckSameSize(X,1,XU,3);
    431:   VecCheckSameSize(X,1,Xout,4);

    433:   VecGetOwnershipRange(X,&low,&high);
    434:   VecGetLocalSize(X,&n);
    435:   if (n>0){
    436:     VecGetArrayRead(X, &x);
    437:     VecGetArrayRead(XL, &xl);
    438:     VecGetArrayRead(XU, &xu);
    439:     VecGetArray(Xout, &xout);

    441:     for (i=0;i<n;++i){
    442:       if ((xl[i] > PETSC_NINFINITY) && (x[i] <= xl[i] + bound_tol)) {
    443:         xout[i] = xl[i]; ++nDiff_loc;
    444:       } else if ((xu[i] < PETSC_INFINITY) && (x[i] >= xu[i] - bound_tol)) {
    445:         xout[i] = xu[i]; ++nDiff_loc;
    446:       }
    447:     }

    449:     VecRestoreArrayRead(X, &x);
    450:     VecRestoreArrayRead(XL, &xl);
    451:     VecRestoreArrayRead(XU, &xu);
    452:     VecRestoreArray(Xout, &xout);
    453:   }
    454:   MPIU_Allreduce(&nDiff_loc, nDiff, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)X));
    455:   return(0);
    456: }