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: /*@
      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
     13: - maskvalue    - the value to set the unused vector elements to (for `TAO_SUBSET_MASK` or `TAO_SUBSET_MATRIXFREE`)

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

     18:   Level: developer

     20:   Notes:
     21:   `maskvalue` should usually be `0.0`, unless a pointwise divide will be used.

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

     36:   PetscFunctionBegin;

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

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

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

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

     74:       PetscCall(VecSet(*vreduced, maskvalue));
     75:       PetscCall(ISGetLocalSize(is, &nlocal));
     76:       PetscCall(VecGetOwnershipRange(vfull, &flow, &fhigh));
     77:       PetscCall(VecGetArray(vfull, &fv));
     78:       PetscCall(VecGetArray(*vreduced, &rv));
     79:       PetscCall(ISGetIndices(is, &s));
     80:       PetscCheck(nlocal <= (fhigh - flow), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "IS local size %" PetscInt_FMT " > Vec local size %" PetscInt_FMT, nlocal, fhigh - flow);
     81:       for (i = 0; i < nlocal; ++i) rv[s[i] - flow] = fv[s[i] - flow];
     82:       PetscCall(ISRestoreIndices(is, &s));
     83:       PetscCall(VecRestoreArray(vfull, &fv));
     84:       PetscCall(VecRestoreArray(*vreduced, &rv));
     85:       break;
     86:     }
     87:   }
     88:   PetscFunctionReturn(PETSC_SUCCESS);
     89: }

     91: /*@
     92:   TaoMatGetSubMat - Gets a submatrix using the `IS`

     94:   Input Parameters:
     95: + M           - the full matrix (`n x n`)
     96: . is          - the index set for the submatrix (both row and column index sets need to be the same)
     97: . v1          - work vector of dimension n, needed for `TAO_SUBSET_MASK` option
     98: - subset_type - the method `Tao` is using for subsetting

    100:   Output Parameter:
    101: . Msub - the submatrix

    103:   Level: developer

    105: .seealso: `TaoVecGetSubVec()`, `TaoSubsetType`
    106: @*/
    107: PetscErrorCode TaoMatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub)
    108: {
    109:   IS        iscomp;
    110:   PetscBool flg = PETSC_TRUE;

    112:   PetscFunctionBegin;
    115:   PetscCall(MatDestroy(Msub));
    116:   switch (subset_type) {
    117:   case TAO_SUBSET_SUBVEC:
    118:     PetscCall(MatCreateSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub));
    119:     break;

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

    139:     /* Zero out rows and columns */
    140:     PetscCall(ISComplementVec(is, v1, &iscomp));

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

    145:     PetscCall(ISDestroy(&iscomp));
    146:     break;
    147:   case TAO_SUBSET_MATRIXFREE:
    148:     PetscCall(ISComplementVec(is, v1, &iscomp));
    149:     PetscCall(MatCreateSubMatrixFree(M, iscomp, iscomp, Msub));
    150:     PetscCall(ISDestroy(&iscomp));
    151:     break;
    152:   }
    153:   PetscFunctionReturn(PETSC_SUCCESS);
    154: }

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

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

    169:   Output Parameters:
    170: + bound_tol    - tolerance for the bound estimation
    171: . active_lower - index set for active variables at the lower bound
    172: . active_upper - index set for active variables at the upper bound
    173: . active_fixed - index set for fixed variables
    174: . active       - index set for all active variables
    175: - inactive     - complementary index set for inactive variables

    177:   Level: developer

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

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

    195:   PetscFunctionBegin;

    203:   if (XL) PetscCheckSameType(X, 1, XL, 2);
    204:   if (XU) PetscCheckSameType(X, 1, XU, 3);
    205:   PetscCheckSameType(X, 1, G, 4);
    206:   PetscCheckSameType(X, 1, S, 5);
    207:   PetscCheckSameType(X, 1, W, 6);
    208:   if (XL) PetscCheckSameComm(X, 1, XL, 2);
    209:   if (XU) PetscCheckSameComm(X, 1, XU, 3);
    210:   PetscCheckSameComm(X, 1, G, 4);
    211:   PetscCheckSameComm(X, 1, S, 5);
    212:   PetscCheckSameComm(X, 1, W, 6);
    213:   if (XL) VecCheckSameSize(X, 1, XL, 2);
    214:   if (XU) VecCheckSameSize(X, 1, XU, 3);
    215:   VecCheckSameSize(X, 1, G, 4);
    216:   VecCheckSameSize(X, 1, S, 5);
    217:   VecCheckSameSize(X, 1, W, 6);

    219:   /* Update the tolerance for bound detection (this is based on Bertsekas' method) */
    220:   PetscCall(VecCopy(X, W));
    221:   PetscCall(VecAXPBY(W, steplen, 1.0, S));
    222:   PetscCall(TaoBoundSolution(W, XL, XU, 0.0, &nDiff, W));
    223:   PetscCall(VecAXPBY(W, 1.0, -1.0, X));
    224:   PetscCall(VecNorm(W, NORM_2, &wnorm));
    225:   *bound_tol = PetscMin(*bound_tol, wnorm);

    227:   /* Clear all index sets */
    228:   PetscCall(ISDestroy(active_lower));
    229:   PetscCall(ISDestroy(active_upper));
    230:   PetscCall(ISDestroy(active_fixed));
    231:   PetscCall(ISDestroy(active));
    232:   PetscCall(ISDestroy(inactive));

    234:   PetscCall(VecGetOwnershipRange(X, &low, &high));
    235:   PetscCall(VecGetLocalSize(X, &n));
    236:   if (!XL && !XU) {
    237:     PetscCall(ISCreateStride(comm, n, low, 1, inactive));
    238:     PetscFunctionReturn(PETSC_SUCCESS);
    239:   }
    240:   if (n > 0) {
    241:     PetscCall(VecGetArrayRead(X, &x));
    242:     PetscCall(VecGetArrayRead(XL, &xl));
    243:     PetscCall(VecGetArrayRead(XU, &xu));
    244:     PetscCall(VecGetArrayRead(G, &g));

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

    278:     PetscCall(VecRestoreArrayRead(X, &x));
    279:     PetscCall(VecRestoreArrayRead(XL, &xl));
    280:     PetscCall(VecRestoreArrayRead(XU, &xu));
    281:     PetscCall(VecRestoreArrayRead(G, &g));
    282:   }

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

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

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

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

    337:   Output Parameter:
    338: . S - step direction to be modified

    340:   Level: developer

    342: .seealso: `TAOBNCG`, `TAOBNTL`, `TAOBNTR`, `TaoBoundSolution()`
    343: @*/
    344: PetscErrorCode TaoBoundStep(Vec X, Vec XL, Vec XU, IS active_lower, IS active_upper, IS active_fixed, PetscReal scale, Vec S)
    345: {
    346:   Vec step_lower, step_upper, step_fixed;
    347:   Vec x_lower, x_upper;
    348:   Vec bound_lower, bound_upper;

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

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

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

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

    389:   Collective

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

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

    401:   Level: developer

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

    411:   PetscFunctionBegin;
    416:   if (!XL && !XU) {
    417:     PetscCall(VecCopy(X, Xout));
    418:     *nDiff = 0.0;
    419:     PetscFunctionReturn(PETSC_SUCCESS);
    420:   }
    421:   PetscCheckSameType(X, 1, XL, 2);
    422:   PetscCheckSameType(X, 1, XU, 3);
    423:   PetscCheckSameType(X, 1, Xout, 6);
    424:   PetscCheckSameComm(X, 1, XL, 2);
    425:   PetscCheckSameComm(X, 1, XU, 3);
    426:   PetscCheckSameComm(X, 1, Xout, 6);
    427:   VecCheckSameSize(X, 1, XL, 2);
    428:   VecCheckSameSize(X, 1, XU, 3);
    429:   VecCheckSameSize(X, 1, Xout, 4);

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

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

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