Actual source code: vecnest.c

      1: #include <../src/vec/vec/impls/nest/vecnestimpl.h>

      3: /* check all blocks are filled */
      4: static PetscErrorCode VecAssemblyBegin_Nest(Vec v)
      5: {
      6:   Vec_Nest *vs = (Vec_Nest *)v->data;
      7:   PetscInt  i;

      9:   PetscFunctionBegin;
     10:   for (i = 0; i < vs->nb; i++) {
     11:     PetscCheck(vs->v[i], PetscObjectComm((PetscObject)v), PETSC_ERR_SUP, "Nest  vector cannot contain NULL blocks");
     12:     PetscCall(VecAssemblyBegin(vs->v[i]));
     13:   }
     14:   PetscFunctionReturn(PETSC_SUCCESS);
     15: }

     17: static PetscErrorCode VecAssemblyEnd_Nest(Vec v)
     18: {
     19:   Vec_Nest *vs = (Vec_Nest *)v->data;
     20:   PetscInt  i;

     22:   PetscFunctionBegin;
     23:   for (i = 0; i < vs->nb; i++) PetscCall(VecAssemblyEnd(vs->v[i]));
     24:   PetscFunctionReturn(PETSC_SUCCESS);
     25: }

     27: static PetscErrorCode VecDestroy_Nest(Vec v)
     28: {
     29:   Vec_Nest *vs = (Vec_Nest *)v->data;
     30:   PetscInt  i;

     32:   PetscFunctionBegin;
     33:   if (vs->v) {
     34:     for (i = 0; i < vs->nb; i++) PetscCall(VecDestroy(&vs->v[i]));
     35:     PetscCall(PetscFree(vs->v));
     36:   }
     37:   for (i = 0; i < vs->nb; i++) PetscCall(ISDestroy(&vs->is[i]));
     38:   PetscCall(PetscFree(vs->is));
     39:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSubVec_C", NULL));
     40:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSubVecs_C", NULL));
     41:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestSetSubVec_C", NULL));
     42:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestSetSubVecs_C", NULL));
     43:   PetscCall(PetscObjectComposeFunction((PetscObject)v, "VecNestGetSize_C", NULL));

     45:   PetscCall(PetscFree(vs));
     46:   PetscFunctionReturn(PETSC_SUCCESS);
     47: }

     49: static PetscErrorCode VecCopy_Nest(Vec x, Vec y)
     50: {
     51:   Vec_Nest *bx = (Vec_Nest *)x->data;
     52:   Vec_Nest *by = (Vec_Nest *)y->data;
     53:   PetscInt  i;

     55:   PetscFunctionBegin;
     56:   VecNestCheckCompatible2(x, 1, y, 2);
     57:   for (i = 0; i < bx->nb; i++) PetscCall(VecCopy(bx->v[i], by->v[i]));
     58:   PetscFunctionReturn(PETSC_SUCCESS);
     59: }

     61: static PetscErrorCode VecDuplicate_Nest(Vec x, Vec *y)
     62: {
     63:   Vec_Nest *bx = (Vec_Nest *)x->data;
     64:   Vec       Y;
     65:   Vec      *sub;
     66:   PetscInt  i;

     68:   PetscFunctionBegin;
     69:   PetscCall(PetscMalloc1(bx->nb, &sub));
     70:   for (i = 0; i < bx->nb; i++) PetscCall(VecDuplicate(bx->v[i], &sub[i]));
     71:   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)x), bx->nb, bx->is, sub, &Y));
     72:   for (i = 0; i < bx->nb; i++) PetscCall(VecDestroy(&sub[i]));
     73:   PetscCall(PetscFree(sub));
     74:   *y = Y;
     75:   PetscFunctionReturn(PETSC_SUCCESS);
     76: }

     78: static PetscErrorCode VecDot_Nest(Vec x, Vec y, PetscScalar *val)
     79: {
     80:   Vec_Nest   *bx = (Vec_Nest *)x->data;
     81:   Vec_Nest   *by = (Vec_Nest *)y->data;
     82:   PetscInt    i, nr;
     83:   PetscScalar x_dot_y, _val;

     85:   PetscFunctionBegin;
     86:   VecNestCheckCompatible2(x, 1, y, 2);
     87:   nr   = bx->nb;
     88:   _val = 0.0;
     89:   for (i = 0; i < nr; i++) {
     90:     PetscCall(VecDot(bx->v[i], by->v[i], &x_dot_y));
     91:     _val = _val + x_dot_y;
     92:   }
     93:   *val = _val;
     94:   PetscFunctionReturn(PETSC_SUCCESS);
     95: }

     97: static PetscErrorCode VecTDot_Nest(Vec x, Vec y, PetscScalar *val)
     98: {
     99:   Vec_Nest   *bx = (Vec_Nest *)x->data;
    100:   Vec_Nest   *by = (Vec_Nest *)y->data;
    101:   PetscInt    i, nr;
    102:   PetscScalar x_dot_y, _val;

    104:   PetscFunctionBegin;
    105:   VecNestCheckCompatible2(x, 1, y, 2);
    106:   nr   = bx->nb;
    107:   _val = 0.0;
    108:   for (i = 0; i < nr; i++) {
    109:     PetscCall(VecTDot(bx->v[i], by->v[i], &x_dot_y));
    110:     _val = _val + x_dot_y;
    111:   }
    112:   *val = _val;
    113:   PetscFunctionReturn(PETSC_SUCCESS);
    114: }

    116: static PetscErrorCode VecDotNorm2_Nest(Vec x, Vec y, PetscScalar *dp, PetscScalar *nm)
    117: {
    118:   Vec_Nest   *bx = (Vec_Nest *)x->data;
    119:   Vec_Nest   *by = (Vec_Nest *)y->data;
    120:   PetscInt    i, nr;
    121:   PetscScalar x_dot_y, _dp, _nm;
    122:   PetscReal   norm2_y;

    124:   PetscFunctionBegin;
    125:   VecNestCheckCompatible2(x, 1, y, 2);
    126:   nr  = bx->nb;
    127:   _dp = 0.0;
    128:   _nm = 0.0;
    129:   for (i = 0; i < nr; i++) {
    130:     PetscCall(VecDotNorm2(bx->v[i], by->v[i], &x_dot_y, &norm2_y));
    131:     _dp += x_dot_y;
    132:     _nm += norm2_y;
    133:   }
    134:   *dp = _dp;
    135:   *nm = _nm;
    136:   PetscFunctionReturn(PETSC_SUCCESS);
    137: }

    139: static PetscErrorCode VecAXPY_Nest(Vec y, PetscScalar alpha, Vec x)
    140: {
    141:   Vec_Nest *bx = (Vec_Nest *)x->data;
    142:   Vec_Nest *by = (Vec_Nest *)y->data;
    143:   PetscInt  i, nr;

    145:   PetscFunctionBegin;
    146:   VecNestCheckCompatible2(y, 1, x, 3);
    147:   nr = bx->nb;
    148:   for (i = 0; i < nr; i++) PetscCall(VecAXPY(by->v[i], alpha, bx->v[i]));
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }

    152: static PetscErrorCode VecAYPX_Nest(Vec y, PetscScalar alpha, Vec x)
    153: {
    154:   Vec_Nest *bx = (Vec_Nest *)x->data;
    155:   Vec_Nest *by = (Vec_Nest *)y->data;
    156:   PetscInt  i, nr;

    158:   PetscFunctionBegin;
    159:   VecNestCheckCompatible2(y, 1, x, 3);
    160:   nr = bx->nb;
    161:   for (i = 0; i < nr; i++) PetscCall(VecAYPX(by->v[i], alpha, bx->v[i]));
    162:   PetscFunctionReturn(PETSC_SUCCESS);
    163: }

    165: static PetscErrorCode VecAXPBY_Nest(Vec y, PetscScalar alpha, PetscScalar beta, Vec x)
    166: {
    167:   Vec_Nest *bx = (Vec_Nest *)x->data;
    168:   Vec_Nest *by = (Vec_Nest *)y->data;
    169:   PetscInt  i, nr;

    171:   PetscFunctionBegin;
    172:   VecNestCheckCompatible2(y, 1, x, 4);
    173:   nr = bx->nb;
    174:   for (i = 0; i < nr; i++) PetscCall(VecAXPBY(by->v[i], alpha, beta, bx->v[i]));
    175:   PetscFunctionReturn(PETSC_SUCCESS);
    176: }

    178: static PetscErrorCode VecAXPBYPCZ_Nest(Vec z, PetscScalar alpha, PetscScalar beta, PetscScalar gamma, Vec x, Vec y)
    179: {
    180:   Vec_Nest *bx = (Vec_Nest *)x->data;
    181:   Vec_Nest *by = (Vec_Nest *)y->data;
    182:   Vec_Nest *bz = (Vec_Nest *)z->data;
    183:   PetscInt  i, nr;

    185:   PetscFunctionBegin;
    186:   VecNestCheckCompatible3(z, 1, x, 5, y, 6);
    187:   nr = bx->nb;
    188:   for (i = 0; i < nr; i++) PetscCall(VecAXPBYPCZ(bz->v[i], alpha, beta, gamma, bx->v[i], by->v[i]));
    189:   PetscFunctionReturn(PETSC_SUCCESS);
    190: }

    192: static PetscErrorCode VecScale_Nest(Vec x, PetscScalar alpha)
    193: {
    194:   Vec_Nest *bx = (Vec_Nest *)x->data;
    195:   PetscInt  i, nr;

    197:   PetscFunctionBegin;
    198:   nr = bx->nb;
    199:   for (i = 0; i < nr; i++) PetscCall(VecScale(bx->v[i], alpha));
    200:   PetscFunctionReturn(PETSC_SUCCESS);
    201: }

    203: static PetscErrorCode VecPointwiseMult_Nest(Vec w, Vec x, Vec y)
    204: {
    205:   Vec_Nest *bx = (Vec_Nest *)x->data;
    206:   Vec_Nest *by = (Vec_Nest *)y->data;
    207:   Vec_Nest *bw = (Vec_Nest *)w->data;
    208:   PetscInt  i, nr;

    210:   PetscFunctionBegin;
    211:   VecNestCheckCompatible3(w, 1, x, 2, y, 3);
    212:   nr = bx->nb;
    213:   for (i = 0; i < nr; i++) PetscCall(VecPointwiseMult(bw->v[i], bx->v[i], by->v[i]));
    214:   PetscFunctionReturn(PETSC_SUCCESS);
    215: }

    217: static PetscErrorCode VecPointwiseDivide_Nest(Vec w, Vec x, Vec y)
    218: {
    219:   Vec_Nest *bx = (Vec_Nest *)x->data;
    220:   Vec_Nest *by = (Vec_Nest *)y->data;
    221:   Vec_Nest *bw = (Vec_Nest *)w->data;
    222:   PetscInt  i, nr;

    224:   PetscFunctionBegin;
    225:   VecNestCheckCompatible3(w, 1, x, 2, y, 3);
    226:   nr = bx->nb;
    227:   for (i = 0; i < nr; i++) PetscCall(VecPointwiseDivide(bw->v[i], bx->v[i], by->v[i]));
    228:   PetscFunctionReturn(PETSC_SUCCESS);
    229: }

    231: static PetscErrorCode VecReciprocal_Nest(Vec x)
    232: {
    233:   Vec_Nest *bx = (Vec_Nest *)x->data;
    234:   PetscInt  i, nr;

    236:   PetscFunctionBegin;
    237:   nr = bx->nb;
    238:   for (i = 0; i < nr; i++) PetscCall(VecReciprocal(bx->v[i]));
    239:   PetscFunctionReturn(PETSC_SUCCESS);
    240: }

    242: static PetscErrorCode VecNorm_Nest(Vec xin, NormType type, PetscReal *z)
    243: {
    244:   Vec_Nest *bx = (Vec_Nest *)xin->data;
    245:   PetscInt  i, nr;
    246:   PetscReal z_i;
    247:   PetscReal _z;

    249:   PetscFunctionBegin;
    250:   nr = bx->nb;
    251:   _z = 0.0;

    253:   if (type == NORM_2) {
    254:     PetscScalar dot;
    255:     PetscCall(VecDot(xin, xin, &dot));
    256:     _z = PetscAbsScalar(PetscSqrtScalar(dot));
    257:   } else if (type == NORM_1) {
    258:     for (i = 0; i < nr; i++) {
    259:       PetscCall(VecNorm(bx->v[i], type, &z_i));
    260:       _z = _z + z_i;
    261:     }
    262:   } else if (type == NORM_INFINITY) {
    263:     for (i = 0; i < nr; i++) {
    264:       PetscCall(VecNorm(bx->v[i], type, &z_i));
    265:       if (z_i > _z) _z = z_i;
    266:     }
    267:   }

    269:   *z = _z;
    270:   PetscFunctionReturn(PETSC_SUCCESS);
    271: }

    273: static PetscErrorCode VecMAXPY_Nest(Vec y, PetscInt nv, const PetscScalar alpha[], Vec *x)
    274: {
    275:   PetscFunctionBegin;
    276:   /* TODO: implement proper MAXPY. For now, do axpy on each vector */
    277:   for (PetscInt v = 0; v < nv; v++) PetscCall(VecAXPY(y, alpha[v], x[v]));
    278:   PetscFunctionReturn(PETSC_SUCCESS);
    279: }

    281: static PetscErrorCode VecMDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
    282: {
    283:   PetscFunctionBegin;
    284:   /* TODO: implement proper MDOT. For now, do dot on each vector */
    285:   for (PetscInt j = 0; j < nv; j++) PetscCall(VecDot(x, y[j], &val[j]));
    286:   PetscFunctionReturn(PETSC_SUCCESS);
    287: }

    289: static PetscErrorCode VecMTDot_Nest(Vec x, PetscInt nv, const Vec y[], PetscScalar *val)
    290: {
    291:   PetscFunctionBegin;
    292:   /* TODO: implement proper MTDOT. For now, do tdot on each vector */
    293:   for (PetscInt j = 0; j < nv; j++) PetscCall(VecTDot(x, y[j], &val[j]));
    294:   PetscFunctionReturn(PETSC_SUCCESS);
    295: }

    297: static PetscErrorCode VecSet_Nest(Vec x, PetscScalar alpha)
    298: {
    299:   Vec_Nest *bx = (Vec_Nest *)x->data;
    300:   PetscInt  i, nr;

    302:   PetscFunctionBegin;
    303:   nr = bx->nb;
    304:   for (i = 0; i < nr; i++) PetscCall(VecSet(bx->v[i], alpha));
    305:   PetscFunctionReturn(PETSC_SUCCESS);
    306: }

    308: static PetscErrorCode VecConjugate_Nest(Vec x)
    309: {
    310:   Vec_Nest *bx = (Vec_Nest *)x->data;
    311:   PetscInt  j, nr;

    313:   PetscFunctionBegin;
    314:   nr = bx->nb;
    315:   for (j = 0; j < nr; j++) PetscCall(VecConjugate(bx->v[j]));
    316:   PetscFunctionReturn(PETSC_SUCCESS);
    317: }

    319: static PetscErrorCode VecSwap_Nest(Vec x, Vec y)
    320: {
    321:   Vec_Nest *bx = (Vec_Nest *)x->data;
    322:   Vec_Nest *by = (Vec_Nest *)y->data;
    323:   PetscInt  i, nr;

    325:   PetscFunctionBegin;
    326:   VecNestCheckCompatible2(x, 1, y, 2);
    327:   nr = bx->nb;
    328:   for (i = 0; i < nr; i++) PetscCall(VecSwap(bx->v[i], by->v[i]));
    329:   PetscFunctionReturn(PETSC_SUCCESS);
    330: }

    332: static PetscErrorCode VecWAXPY_Nest(Vec w, PetscScalar alpha, Vec x, Vec y)
    333: {
    334:   Vec_Nest *bx = (Vec_Nest *)x->data;
    335:   Vec_Nest *by = (Vec_Nest *)y->data;
    336:   Vec_Nest *bw = (Vec_Nest *)w->data;
    337:   PetscInt  i, nr;

    339:   PetscFunctionBegin;
    340:   VecNestCheckCompatible3(w, 1, x, 3, y, 4);
    341:   nr = bx->nb;
    342:   for (i = 0; i < nr; i++) PetscCall(VecWAXPY(bw->v[i], alpha, bx->v[i], by->v[i]));
    343:   PetscFunctionReturn(PETSC_SUCCESS);
    344: }

    346: static PetscErrorCode VecMin_Nest(Vec x, PetscInt *p, PetscReal *v)
    347: {
    348:   PetscInt  i, lp = -1, lw = -1;
    349:   PetscReal lv;
    350:   Vec_Nest *bx = (Vec_Nest *)x->data;

    352:   PetscFunctionBegin;
    353:   *v = PETSC_MAX_REAL;
    354:   for (i = 0; i < bx->nb; i++) {
    355:     PetscInt tp;
    356:     PetscCall(VecMin(bx->v[i], &tp, &lv));
    357:     if (lv < *v) {
    358:       *v = lv;
    359:       lw = i;
    360:       lp = tp;
    361:     }
    362:   }
    363:   if (p && lw > -1) {
    364:     PetscInt        st, en;
    365:     const PetscInt *idxs;

    367:     *p = -1;
    368:     PetscCall(VecGetOwnershipRange(bx->v[lw], &st, &en));
    369:     if (lp >= st && lp < en) {
    370:       PetscCall(ISGetIndices(bx->is[lw], &idxs));
    371:       *p = idxs[lp - st];
    372:       PetscCall(ISRestoreIndices(bx->is[lw], &idxs));
    373:     }
    374:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
    375:   }
    376:   PetscFunctionReturn(PETSC_SUCCESS);
    377: }

    379: static PetscErrorCode VecMax_Nest(Vec x, PetscInt *p, PetscReal *v)
    380: {
    381:   PetscInt  i, lp = -1, lw = -1;
    382:   PetscReal lv;
    383:   Vec_Nest *bx = (Vec_Nest *)x->data;

    385:   PetscFunctionBegin;
    386:   *v = PETSC_MIN_REAL;
    387:   for (i = 0; i < bx->nb; i++) {
    388:     PetscInt tp;
    389:     PetscCall(VecMax(bx->v[i], &tp, &lv));
    390:     if (lv > *v) {
    391:       *v = lv;
    392:       lw = i;
    393:       lp = tp;
    394:     }
    395:   }
    396:   if (p && lw > -1) {
    397:     PetscInt        st, en;
    398:     const PetscInt *idxs;

    400:     *p = -1;
    401:     PetscCall(VecGetOwnershipRange(bx->v[lw], &st, &en));
    402:     if (lp >= st && lp < en) {
    403:       PetscCall(ISGetIndices(bx->is[lw], &idxs));
    404:       *p = idxs[lp - st];
    405:       PetscCall(ISRestoreIndices(bx->is[lw], &idxs));
    406:     }
    407:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, p, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)x)));
    408:   }
    409:   PetscFunctionReturn(PETSC_SUCCESS);
    410: }

    412: static PetscErrorCode VecView_Nest(Vec x, PetscViewer viewer)
    413: {
    414:   Vec_Nest *bx = (Vec_Nest *)x->data;
    415:   PetscBool isascii;
    416:   PetscInt  i;

    418:   PetscFunctionBegin;
    419:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    420:   if (isascii) {
    421:     PetscCall(PetscViewerASCIIPushTab(viewer));
    422:     PetscCall(PetscViewerASCIIPrintf(viewer, "VecNest, rows=%" PetscInt_FMT ",  structure: \n", bx->nb));
    423:     for (i = 0; i < bx->nb; i++) {
    424:       VecType  type;
    425:       char     name[256] = "", prefix[256] = "";
    426:       PetscInt NR;

    428:       PetscCall(VecGetSize(bx->v[i], &NR));
    429:       PetscCall(VecGetType(bx->v[i], &type));
    430:       if (((PetscObject)bx->v[i])->name) PetscCall(PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bx->v[i])->name));
    431:       if (((PetscObject)bx->v[i])->prefix) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bx->v[i])->prefix));

    433:       PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ") : %s%stype=%s, rows=%" PetscInt_FMT " \n", i, name, prefix, type, NR));

    435:       PetscCall(PetscViewerASCIIPushTab(viewer)); /* push1 */
    436:       PetscCall(VecView(bx->v[i], viewer));
    437:       PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop1 */
    438:     }
    439:     PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop0 */
    440:   }
    441:   PetscFunctionReturn(PETSC_SUCCESS);
    442: }

    444: static PetscErrorCode VecSize_Nest_Recursive(Vec x, PetscBool globalsize, PetscInt *L)
    445: {
    446:   Vec_Nest *bx;
    447:   PetscInt  size, i, nr;
    448:   PetscBool isnest;

    450:   PetscFunctionBegin;
    451:   PetscCall(PetscObjectTypeCompare((PetscObject)x, VECNEST, &isnest));
    452:   if (!isnest) {
    453:     /* Not nest */
    454:     if (globalsize) PetscCall(VecGetSize(x, &size));
    455:     else PetscCall(VecGetLocalSize(x, &size));
    456:     *L = *L + size;
    457:     PetscFunctionReturn(PETSC_SUCCESS);
    458:   }

    460:   /* Otherwise we have a nest */
    461:   bx = (Vec_Nest *)x->data;
    462:   nr = bx->nb;

    464:   /* now descend recursively */
    465:   for (i = 0; i < nr; i++) PetscCall(VecSize_Nest_Recursive(bx->v[i], globalsize, L));
    466:   PetscFunctionReturn(PETSC_SUCCESS);
    467: }

    469: /* Returns the sum of the global size of all the constituent vectors in the nest */
    470: static PetscErrorCode VecGetSize_Nest(Vec x, PetscInt *N)
    471: {
    472:   PetscFunctionBegin;
    473:   *N = x->map->N;
    474:   PetscFunctionReturn(PETSC_SUCCESS);
    475: }

    477: /* Returns the sum of the local size of all the constituent vectors in the nest */
    478: static PetscErrorCode VecGetLocalSize_Nest(Vec x, PetscInt *n)
    479: {
    480:   PetscFunctionBegin;
    481:   *n = x->map->n;
    482:   PetscFunctionReturn(PETSC_SUCCESS);
    483: }

    485: static PetscErrorCode VecMaxPointwiseDivide_Nest(Vec x, Vec y, PetscReal *max)
    486: {
    487:   Vec_Nest *bx = (Vec_Nest *)x->data;
    488:   Vec_Nest *by = (Vec_Nest *)y->data;
    489:   PetscInt  i, nr;
    490:   PetscReal local_max, m;

    492:   PetscFunctionBegin;
    493:   VecNestCheckCompatible2(x, 1, y, 2);
    494:   nr = bx->nb;
    495:   m  = 0.0;
    496:   for (i = 0; i < nr; i++) {
    497:     PetscCall(VecMaxPointwiseDivide(bx->v[i], by->v[i], &local_max));
    498:     if (local_max > m) m = local_max;
    499:   }
    500:   *max = m;
    501:   PetscFunctionReturn(PETSC_SUCCESS);
    502: }

    504: static PetscErrorCode VecGetSubVector_Nest(Vec X, IS is, Vec *x)
    505: {
    506:   Vec_Nest *bx = (Vec_Nest *)X->data;
    507:   PetscInt  i;

    509:   PetscFunctionBegin;
    510:   *x = NULL;
    511:   for (i = 0; i < bx->nb; i++) {
    512:     PetscBool issame = PETSC_FALSE;
    513:     PetscCall(ISEqual(is, bx->is[i], &issame));
    514:     if (issame) {
    515:       *x = bx->v[i];
    516:       PetscCall(PetscObjectReference((PetscObject)*x));
    517:       break;
    518:     }
    519:   }
    520:   PetscCheck(*x, PetscObjectComm((PetscObject)is), PETSC_ERR_ARG_OUTOFRANGE, "Index set not found in nested Vec");
    521:   PetscFunctionReturn(PETSC_SUCCESS);
    522: }

    524: static PetscErrorCode VecRestoreSubVector_Nest(Vec X, IS is, Vec *x)
    525: {
    526:   PetscFunctionBegin;
    527:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
    528:   PetscCall(VecDestroy(x));
    529:   PetscFunctionReturn(PETSC_SUCCESS);
    530: }

    532: static PetscErrorCode VecGetArray_Nest(Vec X, PetscScalar **x)
    533: {
    534:   Vec_Nest *bx = (Vec_Nest *)X->data;
    535:   PetscInt  i, m, rstart, rend;

    537:   PetscFunctionBegin;
    538:   PetscCall(VecGetOwnershipRange(X, &rstart, &rend));
    539:   PetscCall(VecGetLocalSize(X, &m));
    540:   PetscCall(PetscMalloc1(m, x));
    541:   for (i = 0; i < bx->nb; i++) {
    542:     Vec                subvec = bx->v[i];
    543:     IS                 isy    = bx->is[i];
    544:     PetscInt           j, sm;
    545:     const PetscInt    *ixy;
    546:     const PetscScalar *y;
    547:     PetscCall(VecGetLocalSize(subvec, &sm));
    548:     PetscCall(VecGetArrayRead(subvec, &y));
    549:     PetscCall(ISGetIndices(isy, &ixy));
    550:     for (j = 0; j < sm; j++) {
    551:       PetscInt ix = ixy[j];
    552:       PetscCheck(ix >= rstart && rend > ix, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for getting array from nonlocal subvec");
    553:       (*x)[ix - rstart] = y[j];
    554:     }
    555:     PetscCall(ISRestoreIndices(isy, &ixy));
    556:     PetscCall(VecRestoreArrayRead(subvec, &y));
    557:   }
    558:   PetscFunctionReturn(PETSC_SUCCESS);
    559: }

    561: static PetscErrorCode VecRestoreArray_Nest(Vec X, PetscScalar **x)
    562: {
    563:   Vec_Nest *bx = (Vec_Nest *)X->data;
    564:   PetscInt  i, m, rstart, rend;

    566:   PetscFunctionBegin;
    567:   PetscCall(VecGetOwnershipRange(X, &rstart, &rend));
    568:   PetscCall(VecGetLocalSize(X, &m));
    569:   for (i = 0; i < bx->nb; i++) {
    570:     Vec             subvec = bx->v[i];
    571:     IS              isy    = bx->is[i];
    572:     PetscInt        j, sm;
    573:     const PetscInt *ixy;
    574:     PetscScalar    *y;
    575:     PetscCall(VecGetLocalSize(subvec, &sm));
    576:     PetscCall(VecGetArray(subvec, &y));
    577:     PetscCall(ISGetIndices(isy, &ixy));
    578:     for (j = 0; j < sm; j++) {
    579:       PetscInt ix = ixy[j];
    580:       PetscCheck(ix >= rstart && rend > ix, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for getting array from nonlocal subvec");
    581:       y[j] = (*x)[ix - rstart];
    582:     }
    583:     PetscCall(ISRestoreIndices(isy, &ixy));
    584:     PetscCall(VecRestoreArray(subvec, &y));
    585:   }
    586:   PetscCall(PetscFree(*x));
    587:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
    588:   PetscFunctionReturn(PETSC_SUCCESS);
    589: }

    591: static PetscErrorCode VecRestoreArrayRead_Nest(Vec X, const PetscScalar **x)
    592: {
    593:   PetscFunctionBegin;
    594:   PetscCall(PetscFree(*x));
    595:   PetscFunctionReturn(PETSC_SUCCESS);
    596: }

    598: static PetscErrorCode VecConcatenate_Nest(PetscInt nx, const Vec X[], Vec *Y, IS *x_is[])
    599: {
    600:   PetscFunctionBegin;
    601:   PetscCheck(nx <= 0, PetscObjectComm((PetscObject)*X), PETSC_ERR_SUP, "VecConcatenate() is not supported for VecNest");
    602:   PetscFunctionReturn(PETSC_SUCCESS);
    603: }

    605: static PetscErrorCode VecCreateLocalVector_Nest(Vec v, Vec *w)
    606: {
    607:   Vec        *ww;
    608:   IS         *wis;
    609:   Vec_Nest   *bv = (Vec_Nest *)v->data;
    610:   PetscInt    i;
    611:   PetscMPIInt size;

    613:   PetscFunctionBegin;
    614:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
    615:   if (size == 1) {
    616:     PetscCall(VecDuplicate(v, w));
    617:     PetscFunctionReturn(PETSC_SUCCESS);
    618:   }
    619:   PetscCall(PetscMalloc2(bv->nb, &ww, bv->nb, &wis));
    620:   for (i = 0; i < bv->nb; i++) PetscCall(VecCreateLocalVector(bv->v[i], &ww[i]));
    621:   for (i = 0; i < bv->nb; i++) {
    622:     PetscInt off;

    624:     PetscCall(VecGetOwnershipRange(v, &off, NULL));
    625:     PetscCall(ISOnComm(bv->is[i], PetscObjectComm((PetscObject)ww[i]), PETSC_COPY_VALUES, &wis[i]));
    626:     PetscCall(ISShift(wis[i], -off, wis[i]));
    627:   }
    628:   PetscCall(VecCreateNest(PETSC_COMM_SELF, bv->nb, wis, ww, w));
    629:   for (i = 0; i < bv->nb; i++) {
    630:     PetscCall(VecDestroy(&ww[i]));
    631:     PetscCall(ISDestroy(&wis[i]));
    632:   }
    633:   PetscCall(PetscFree2(ww, wis));
    634:   PetscFunctionReturn(PETSC_SUCCESS);
    635: }

    637: static PetscErrorCode VecGetLocalVector_Nest(Vec v, Vec w)
    638: {
    639:   Vec_Nest *bv = (Vec_Nest *)v->data;
    640:   Vec_Nest *bw = (Vec_Nest *)w->data;
    641:   PetscInt  i;

    643:   PetscFunctionBegin;
    644:   PetscCheckSameType(v, 1, w, 2);
    645:   PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
    646:   for (i = 0; i < bv->nb; i++) PetscCall(VecGetLocalVector(bv->v[i], bw->v[i]));
    647:   PetscFunctionReturn(PETSC_SUCCESS);
    648: }

    650: static PetscErrorCode VecRestoreLocalVector_Nest(Vec v, Vec w)
    651: {
    652:   Vec_Nest *bv = (Vec_Nest *)v->data;
    653:   Vec_Nest *bw = (Vec_Nest *)w->data;
    654:   PetscInt  i;

    656:   PetscFunctionBegin;
    657:   PetscCheckSameType(v, 1, w, 2);
    658:   PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
    659:   for (i = 0; i < bv->nb; i++) PetscCall(VecRestoreLocalVector(bv->v[i], bw->v[i]));
    660:   PetscCall(PetscObjectStateIncrease((PetscObject)v));
    661:   PetscFunctionReturn(PETSC_SUCCESS);
    662: }

    664: static PetscErrorCode VecGetLocalVectorRead_Nest(Vec v, Vec w)
    665: {
    666:   Vec_Nest *bv = (Vec_Nest *)v->data;
    667:   Vec_Nest *bw = (Vec_Nest *)w->data;
    668:   PetscInt  i;

    670:   PetscFunctionBegin;
    671:   PetscCheckSameType(v, 1, w, 2);
    672:   PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
    673:   for (i = 0; i < bv->nb; i++) PetscCall(VecGetLocalVectorRead(bv->v[i], bw->v[i]));
    674:   PetscFunctionReturn(PETSC_SUCCESS);
    675: }

    677: static PetscErrorCode VecRestoreLocalVectorRead_Nest(Vec v, Vec w)
    678: {
    679:   Vec_Nest *bv = (Vec_Nest *)v->data;
    680:   Vec_Nest *bw = (Vec_Nest *)w->data;
    681:   PetscInt  i;

    683:   PetscFunctionBegin;
    684:   PetscCheckSameType(v, 1, w, 2);
    685:   PetscCheck(bv->nb = bw->nb, PetscObjectComm((PetscObject)w), PETSC_ERR_ARG_WRONG, "Invalid local vector");
    686:   for (i = 0; i < bv->nb; i++) PetscCall(VecRestoreLocalVectorRead(bv->v[i], bw->v[i]));
    687:   PetscFunctionReturn(PETSC_SUCCESS);
    688: }

    690: static PetscErrorCode VecSetRandom_Nest(Vec v, PetscRandom r)
    691: {
    692:   Vec_Nest *bv = (Vec_Nest *)v->data;

    694:   PetscFunctionBegin;
    695:   for (PetscInt i = 0; i < bv->nb; i++) PetscCall(VecSetRandom(bv->v[i], r));
    696:   PetscFunctionReturn(PETSC_SUCCESS);
    697: }

    699: static PetscErrorCode VecErrorWeightedNorms_Nest(Vec U, Vec Y, Vec E, NormType wnormtype, PetscReal atol, Vec vatol, PetscReal rtol, Vec vrtol, PetscReal ignore_max, PetscReal *norm, PetscInt *norm_loc, PetscReal *norma, PetscInt *norma_loc, PetscReal *normr, PetscInt *normr_loc)
    700: {
    701:   Vec_Nest *bu = (Vec_Nest *)U->data;
    702:   Vec_Nest *by = (Vec_Nest *)Y->data;
    703:   Vec_Nest *be = E ? (Vec_Nest *)E->data : NULL;
    704:   Vec_Nest *ba = vatol ? (Vec_Nest *)vatol->data : NULL;
    705:   Vec_Nest *br = vrtol ? (Vec_Nest *)vrtol->data : NULL;

    707:   PetscFunctionBegin;
    708:   VecNestCheckCompatible2(U, 1, Y, 2);
    709:   if (E) VecNestCheckCompatible2(U, 1, E, 3);
    710:   if (vatol) VecNestCheckCompatible2(U, 1, vatol, 6);
    711:   if (vrtol) VecNestCheckCompatible2(U, 1, vrtol, 8);
    712:   *norm      = 0.0;
    713:   *norma     = 0.0;
    714:   *normr     = 0.0;
    715:   *norm_loc  = 0;
    716:   *norma_loc = 0;
    717:   *normr_loc = 0;
    718:   for (PetscInt i = 0; i < bu->nb; i++) {
    719:     PetscReal n, na, nr;
    720:     PetscInt  n_loc, na_loc, nr_loc;

    722:     PetscCall(VecErrorWeightedNorms(bu->v[i], by->v[i], be ? be->v[i] : NULL, wnormtype, atol, ba ? ba->v[i] : NULL, rtol, br ? br->v[i] : NULL, ignore_max, &n, &n_loc, &na, &na_loc, &nr, &nr_loc));
    723:     if (wnormtype == NORM_INFINITY) {
    724:       *norm  = PetscMax(*norm, n);
    725:       *norma = PetscMax(*norma, na);
    726:       *normr = PetscMax(*normr, nr);
    727:     } else {
    728:       *norm += PetscSqr(n);
    729:       *norma += PetscSqr(na);
    730:       *normr += PetscSqr(nr);
    731:     }
    732:     *norm_loc += n_loc;
    733:     *norma_loc += na_loc;
    734:     *normr_loc += nr_loc;
    735:   }
    736:   if (wnormtype == NORM_2) {
    737:     *norm  = PetscSqrtReal(*norm);
    738:     *norma = PetscSqrtReal(*norma);
    739:     *normr = PetscSqrtReal(*normr);
    740:   }
    741:   PetscFunctionReturn(PETSC_SUCCESS);
    742: }

    744: static PetscErrorCode VecNestSetOps_Private(struct _VecOps *ops)
    745: {
    746:   PetscFunctionBegin;
    747:   ops->duplicate               = VecDuplicate_Nest;
    748:   ops->duplicatevecs           = VecDuplicateVecs_Default;
    749:   ops->destroyvecs             = VecDestroyVecs_Default;
    750:   ops->dot                     = VecDot_Nest;
    751:   ops->mdot                    = VecMDot_Nest;
    752:   ops->norm                    = VecNorm_Nest;
    753:   ops->tdot                    = VecTDot_Nest;
    754:   ops->mtdot                   = VecMTDot_Nest;
    755:   ops->scale                   = VecScale_Nest;
    756:   ops->copy                    = VecCopy_Nest;
    757:   ops->set                     = VecSet_Nest;
    758:   ops->swap                    = VecSwap_Nest;
    759:   ops->axpy                    = VecAXPY_Nest;
    760:   ops->axpby                   = VecAXPBY_Nest;
    761:   ops->maxpy                   = VecMAXPY_Nest;
    762:   ops->aypx                    = VecAYPX_Nest;
    763:   ops->waxpy                   = VecWAXPY_Nest;
    764:   ops->axpbypcz                = NULL;
    765:   ops->pointwisemult           = VecPointwiseMult_Nest;
    766:   ops->pointwisedivide         = VecPointwiseDivide_Nest;
    767:   ops->setvalues               = NULL;
    768:   ops->assemblybegin           = VecAssemblyBegin_Nest;
    769:   ops->assemblyend             = VecAssemblyEnd_Nest;
    770:   ops->getarray                = VecGetArray_Nest;
    771:   ops->getsize                 = VecGetSize_Nest;
    772:   ops->getlocalsize            = VecGetLocalSize_Nest;
    773:   ops->restorearray            = VecRestoreArray_Nest;
    774:   ops->restorearrayread        = VecRestoreArrayRead_Nest;
    775:   ops->max                     = VecMax_Nest;
    776:   ops->min                     = VecMin_Nest;
    777:   ops->setrandom               = NULL;
    778:   ops->setoption               = NULL;
    779:   ops->setvaluesblocked        = NULL;
    780:   ops->destroy                 = VecDestroy_Nest;
    781:   ops->view                    = VecView_Nest;
    782:   ops->placearray              = NULL;
    783:   ops->replacearray            = NULL;
    784:   ops->dot_local               = VecDot_Nest;
    785:   ops->tdot_local              = VecTDot_Nest;
    786:   ops->norm_local              = VecNorm_Nest;
    787:   ops->mdot_local              = VecMDot_Nest;
    788:   ops->mtdot_local             = VecMTDot_Nest;
    789:   ops->load                    = NULL;
    790:   ops->reciprocal              = VecReciprocal_Nest;
    791:   ops->conjugate               = VecConjugate_Nest;
    792:   ops->setlocaltoglobalmapping = NULL;
    793:   ops->resetarray              = NULL;
    794:   ops->setfromoptions          = NULL;
    795:   ops->maxpointwisedivide      = VecMaxPointwiseDivide_Nest;
    796:   ops->load                    = NULL;
    797:   ops->pointwisemax            = NULL;
    798:   ops->pointwisemaxabs         = NULL;
    799:   ops->pointwisemin            = NULL;
    800:   ops->getvalues               = NULL;
    801:   ops->sqrt                    = NULL;
    802:   ops->abs                     = NULL;
    803:   ops->exp                     = NULL;
    804:   ops->shift                   = NULL;
    805:   ops->create                  = NULL;
    806:   ops->stridegather            = NULL;
    807:   ops->stridescatter           = NULL;
    808:   ops->dotnorm2                = VecDotNorm2_Nest;
    809:   ops->getsubvector            = VecGetSubVector_Nest;
    810:   ops->restoresubvector        = VecRestoreSubVector_Nest;
    811:   ops->axpbypcz                = VecAXPBYPCZ_Nest;
    812:   ops->concatenate             = VecConcatenate_Nest;
    813:   ops->createlocalvector       = VecCreateLocalVector_Nest;
    814:   ops->getlocalvector          = VecGetLocalVector_Nest;
    815:   ops->getlocalvectorread      = VecGetLocalVectorRead_Nest;
    816:   ops->restorelocalvector      = VecRestoreLocalVector_Nest;
    817:   ops->restorelocalvectorread  = VecRestoreLocalVectorRead_Nest;
    818:   ops->setrandom               = VecSetRandom_Nest;
    819:   ops->errorwnorm              = VecErrorWeightedNorms_Nest;
    820:   PetscFunctionReturn(PETSC_SUCCESS);
    821: }

    823: static PetscErrorCode VecNestGetSubVecs_Private(Vec x, PetscInt m, const PetscInt idxm[], Vec vec[])
    824: {
    825:   Vec_Nest *b = (Vec_Nest *)x->data;
    826:   PetscInt  i;
    827:   PetscInt  row;

    829:   PetscFunctionBegin;
    830:   if (!m) PetscFunctionReturn(PETSC_SUCCESS);
    831:   for (i = 0; i < m; i++) {
    832:     row = idxm[i];
    833:     PetscCheck(row < b->nb, PetscObjectComm((PetscObject)x), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, b->nb - 1);
    834:     vec[i] = b->v[row];
    835:   }
    836:   PetscFunctionReturn(PETSC_SUCCESS);
    837: }

    839: static PetscErrorCode VecNestGetSubVec_Nest(Vec X, PetscInt idxm, Vec *sx)
    840: {
    841:   PetscFunctionBegin;
    842:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
    843:   PetscCall(VecNestGetSubVecs_Private(X, 1, &idxm, sx));
    844:   PetscFunctionReturn(PETSC_SUCCESS);
    845: }

    847: /*@
    848:   VecNestGetSubVec - Returns a single, sub-vector from a nest vector.

    850:   Not Collective

    852:   Input Parameters:
    853: + X    - nest vector
    854: - idxm - index of the vector within the nest

    856:   Output Parameter:
    857: . sx - vector at index `idxm` within the nest

    859:   Level: developer

    861: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVecs()`
    862: @*/
    863: PetscErrorCode VecNestGetSubVec(Vec X, PetscInt idxm, Vec *sx)
    864: {
    865:   PetscFunctionBegin;
    866:   PetscUseMethod(X, "VecNestGetSubVec_C", (Vec, PetscInt, Vec *), (X, idxm, sx));
    867:   PetscFunctionReturn(PETSC_SUCCESS);
    868: }

    870: static PetscErrorCode VecNestGetSubVecs_Nest(Vec X, PetscInt *N, Vec **sx)
    871: {
    872:   Vec_Nest *b = (Vec_Nest *)X->data;

    874:   PetscFunctionBegin;
    875:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
    876:   if (N) *N = b->nb;
    877:   if (sx) *sx = b->v;
    878:   PetscFunctionReturn(PETSC_SUCCESS);
    879: }

    881: /*@C
    882:   VecNestGetSubVecs - Returns the entire array of vectors defining a nest vector.

    884:   Not Collective

    886:   Input Parameter:
    887: . X - nest vector

    889:   Output Parameters:
    890: + N  - number of nested vecs
    891: - sx - array of vectors, can pass in `NULL`

    893:   Level: developer

    895:   Note:
    896:   The user should not free the array `sx`.

    898:   Fortran Notes:
    899:   The caller must allocate the array to hold the subvectors and pass it in.

    901: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
    902: @*/
    903: PetscErrorCode VecNestGetSubVecs(Vec X, PetscInt *N, Vec *sx[])
    904: {
    905:   PetscFunctionBegin;
    906:   PetscUseMethod(X, "VecNestGetSubVecs_C", (Vec, PetscInt *, Vec **), (X, N, sx));
    907:   PetscFunctionReturn(PETSC_SUCCESS);
    908: }

    910: static PetscErrorCode VecNestSetSubVec_Private(Vec X, PetscInt idxm, Vec x)
    911: {
    912:   Vec_Nest *bx = (Vec_Nest *)X->data;
    913:   PetscInt  i, offset = 0, n = 0, bs;
    914:   IS        is;
    915:   PetscBool issame = PETSC_FALSE;
    916:   PetscInt  N      = 0;

    918:   PetscFunctionBegin;
    919:   /* check if idxm < bx->nb */
    920:   PetscCheck(idxm < bx->nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Out of range index value %" PetscInt_FMT ", should be less than %" PetscInt_FMT, idxm, bx->nb);

    922:   PetscCall(PetscObjectReference((PetscObject)x));
    923:   PetscCall(VecDestroy(&bx->v[idxm]));
    924:   bx->v[idxm] = x;

    926:   /* check if we need to update the IS for the block */
    927:   offset = X->map->rstart;
    928:   for (i = 0; i < idxm; i++) {
    929:     n = 0;
    930:     PetscCall(VecGetLocalSize(bx->v[i], &n));
    931:     offset += n;
    932:   }

    934:   /* get the local size and block size */
    935:   PetscCall(VecGetLocalSize(x, &n));
    936:   PetscCall(VecGetBlockSize(x, &bs));

    938:   /* create the new IS */
    939:   PetscCall(ISCreateStride(PetscObjectComm((PetscObject)x), n, offset, 1, &is));
    940:   PetscCall(ISSetBlockSize(is, bs));

    942:   /* check if they are equal */
    943:   PetscCall(ISEqual(is, bx->is[idxm], &issame));

    945:   if (!issame) {
    946:     /* The IS of given vector has a different layout compared to the existing block vector.
    947:      Destroy the existing reference and update the IS. */
    948:     PetscCall(ISDestroy(&bx->is[idxm]));
    949:     PetscCall(ISDuplicate(is, &bx->is[idxm]));
    950:     PetscCall(ISCopy(is, bx->is[idxm]));

    952:     offset += n;
    953:     /* Since the current IS[idxm] changed, we need to update all the subsequent IS */
    954:     for (i = idxm + 1; i < bx->nb; i++) {
    955:       /* get the local size and block size */
    956:       PetscCall(VecGetLocalSize(bx->v[i], &n));
    957:       PetscCall(VecGetBlockSize(bx->v[i], &bs));

    959:       /* destroy the old and create the new IS */
    960:       PetscCall(ISDestroy(&bx->is[i]));
    961:       PetscCall(ISCreateStride(((PetscObject)bx->v[i])->comm, n, offset, 1, &bx->is[i]));
    962:       PetscCall(ISSetBlockSize(bx->is[i], bs));

    964:       offset += n;
    965:     }

    967:     n = 0;
    968:     PetscCall(VecSize_Nest_Recursive(X, PETSC_TRUE, &N));
    969:     PetscCall(VecSize_Nest_Recursive(X, PETSC_FALSE, &n));
    970:     PetscCall(PetscLayoutSetSize(X->map, N));
    971:     PetscCall(PetscLayoutSetLocalSize(X->map, n));
    972:   }

    974:   PetscCall(ISDestroy(&is));
    975:   PetscFunctionReturn(PETSC_SUCCESS);
    976: }

    978: static PetscErrorCode VecNestSetSubVec_Nest(Vec X, PetscInt idxm, Vec sx)
    979: {
    980:   PetscFunctionBegin;
    981:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
    982:   PetscCall(VecNestSetSubVec_Private(X, idxm, sx));
    983:   PetscFunctionReturn(PETSC_SUCCESS);
    984: }

    986: /*@
    987:   VecNestSetSubVec - Set a single component vector in a nest vector at specified index.

    989:   Not Collective

    991:   Input Parameters:
    992: + X    - nest vector
    993: . idxm - index of the vector within the nest vector
    994: - sx   - vector at index `idxm` within the nest vector

    996:   Level: developer

    998:   Notes:
    999:   The new vector `sx` does not have to be of same size as X[idxm]. Arbitrary vector layouts are allowed.

    1001:   The nest vector `X` keeps a reference to `sx` rather than creating a duplicate.

    1003: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestSetSubVecs()`, `VecNestGetSubVec()`
    1004: @*/
    1005: PetscErrorCode VecNestSetSubVec(Vec X, PetscInt idxm, Vec sx)
    1006: {
    1007:   PetscFunctionBegin;
    1008:   PetscUseMethod(X, "VecNestSetSubVec_C", (Vec, PetscInt, Vec), (X, idxm, sx));
    1009:   PetscFunctionReturn(PETSC_SUCCESS);
    1010: }

    1012: static PetscErrorCode VecNestSetSubVecs_Nest(Vec X, PetscInt N, PetscInt *idxm, Vec *sx)
    1013: {
    1014:   PetscInt i;

    1016:   PetscFunctionBegin;
    1017:   PetscCall(PetscObjectStateIncrease((PetscObject)X));
    1018:   for (i = 0; i < N; i++) PetscCall(VecNestSetSubVec_Private(X, idxm[i], sx[i]));
    1019:   PetscFunctionReturn(PETSC_SUCCESS);
    1020: }

    1022: /*@
    1023:   VecNestSetSubVecs - Sets the component vectors at the specified indices in a nest vector.

    1025:   Not Collective

    1027:   Input Parameters:
    1028: + X    - nest vector
    1029: . N    - number of component vecs in `sx`
    1030: . idxm - indices of component vectors that are to be replaced
    1031: - sx   - array of vectors

    1033:   Level: developer

    1035:   Notes:
    1036:   The components in the vector array `sx` do not have to be of the same size as corresponding
    1037:   components in `X`. The user can also free the array `sx` after the call.

    1039:   The nest vector `X` keeps references to `sx` vectors rather than creating duplicates.

    1041: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestGetSize()`, `VecNestGetSubVec()`
    1042: @*/
    1043: PetscErrorCode VecNestSetSubVecs(Vec X, PetscInt N, PetscInt idxm[], Vec sx[])
    1044: {
    1045:   PetscFunctionBegin;
    1046:   PetscUseMethod(X, "VecNestSetSubVecs_C", (Vec, PetscInt, PetscInt *, Vec *), (X, N, idxm, sx));
    1047:   PetscFunctionReturn(PETSC_SUCCESS);
    1048: }

    1050: static PetscErrorCode VecNestGetSize_Nest(Vec X, PetscInt *N)
    1051: {
    1052:   Vec_Nest *b = (Vec_Nest *)X->data;

    1054:   PetscFunctionBegin;
    1055:   *N = b->nb;
    1056:   PetscFunctionReturn(PETSC_SUCCESS);
    1057: }

    1059: /*@
    1060:   VecNestGetSize - Returns the size of the nest vector.

    1062:   Not Collective

    1064:   Input Parameter:
    1065: . X - nest vector

    1067:   Output Parameter:
    1068: . N - number of nested vecs

    1070:   Level: developer

    1072: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecNestGetSubVec()`, `VecNestGetSubVecs()`
    1073: @*/
    1074: PetscErrorCode VecNestGetSize(Vec X, PetscInt *N)
    1075: {
    1076:   PetscFunctionBegin;
    1078:   PetscAssertPointer(N, 2);
    1079:   PetscUseMethod(X, "VecNestGetSize_C", (Vec, PetscInt *), (X, N));
    1080:   PetscFunctionReturn(PETSC_SUCCESS);
    1081: }

    1083: static PetscErrorCode VecSetUp_Nest_Private(Vec V, PetscInt nb, Vec x[])
    1084: {
    1085:   Vec_Nest *ctx = (Vec_Nest *)V->data;
    1086:   PetscInt  i;

    1088:   PetscFunctionBegin;
    1089:   if (ctx->setup_called) PetscFunctionReturn(PETSC_SUCCESS);

    1091:   ctx->nb = nb;
    1092:   PetscCheck(ctx->nb >= 0, PetscObjectComm((PetscObject)V), PETSC_ERR_ARG_WRONG, "Cannot create VECNEST with < 0 blocks.");

    1094:   /* Create space */
    1095:   PetscCall(PetscMalloc1(ctx->nb, &ctx->v));
    1096:   for (i = 0; i < ctx->nb; i++) {
    1097:     ctx->v[i] = x[i];
    1098:     PetscCall(PetscObjectReference((PetscObject)x[i]));
    1099:     /* Do not allocate memory for internal sub blocks */
    1100:   }

    1102:   PetscCall(PetscMalloc1(ctx->nb, &ctx->is));

    1104:   ctx->setup_called = PETSC_TRUE;
    1105:   PetscFunctionReturn(PETSC_SUCCESS);
    1106: }

    1108: static PetscErrorCode VecSetUp_NestIS_Private(Vec V, PetscInt nb, IS is[])
    1109: {
    1110:   Vec_Nest *ctx = (Vec_Nest *)V->data;
    1111:   PetscInt  i, offset, m, n, M, N;

    1113:   PetscFunctionBegin;
    1114:   (void)nb;
    1115:   if (is) { /* Do some consistency checks and reference the is */
    1116:     offset = V->map->rstart;
    1117:     for (i = 0; i < ctx->nb; i++) {
    1118:       PetscCall(ISGetSize(is[i], &M));
    1119:       PetscCall(VecGetSize(ctx->v[i], &N));
    1120:       PetscCheck(M == N, PetscObjectComm((PetscObject)V), PETSC_ERR_ARG_INCOMP, "In slot %" PetscInt_FMT ", IS of size %" PetscInt_FMT " is not compatible with Vec of size %" PetscInt_FMT, i, M, N);
    1121:       PetscCall(ISGetLocalSize(is[i], &m));
    1122:       PetscCall(VecGetLocalSize(ctx->v[i], &n));
    1123:       PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "In slot %" PetscInt_FMT ", IS of local size %" PetscInt_FMT " is not compatible with Vec of local size %" PetscInt_FMT, i, m, n);
    1124:       PetscCall(PetscObjectReference((PetscObject)is[i]));
    1125:       ctx->is[i] = is[i];
    1126:       offset += n;
    1127:     }
    1128:   } else { /* Create a contiguous ISStride for each entry */
    1129:     offset = V->map->rstart;
    1130:     for (i = 0; i < ctx->nb; i++) {
    1131:       PetscInt bs;
    1132:       PetscCall(VecGetLocalSize(ctx->v[i], &n));
    1133:       PetscCall(VecGetBlockSize(ctx->v[i], &bs));
    1134:       PetscCall(ISCreateStride(((PetscObject)ctx->v[i])->comm, n, offset, 1, &ctx->is[i]));
    1135:       PetscCall(ISSetBlockSize(ctx->is[i], bs));
    1136:       offset += n;
    1137:     }
    1138:   }
    1139:   PetscFunctionReturn(PETSC_SUCCESS);
    1140: }

    1142: /*MC
    1143:   VECNEST - VECNEST = "nest" - Vector type consisting of nested subvectors, each stored separately.

    1145:   Level: intermediate

    1147:   Notes:
    1148:   This vector type reduces the number of copies for certain solvers applied to multi-physics problems.
    1149:   It is usually used with `MATNEST` and `DMCOMPOSITE` via `DMSetVecType()`.

    1151: .seealso: [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `VecType`, `VecCreateNest()`, `MatCreateNest()`
    1152: M*/

    1154: /*@C
    1155:   VecCreateNest - Creates a new vector containing several nested subvectors, each stored separately

    1157:   Collective

    1159:   Input Parameters:
    1160: + comm - Communicator for the new `Vec`
    1161: . nb   - number of nested blocks
    1162: . is   - array of `nb` index sets describing each nested block, or `NULL` to pack subvectors contiguously
    1163: - x    - array of `nb` sub-vectors

    1165:   Output Parameter:
    1166: . Y - new vector

    1168:   Level: advanced

    1170: .seealso: `VECNEST`,  [](ch_vectors), `Vec`, `VecType`, `VecCreate()`, `MatCreateNest()`, `DMSetVecType()`
    1171: @*/
    1172: PetscErrorCode VecCreateNest(MPI_Comm comm, PetscInt nb, IS is[], Vec x[], Vec *Y)
    1173: {
    1174:   Vec       V;
    1175:   Vec_Nest *s;
    1176:   PetscInt  n, N;

    1178:   PetscFunctionBegin;
    1179:   PetscCall(VecCreate(comm, &V));
    1180:   PetscCall(PetscObjectChangeTypeName((PetscObject)V, VECNEST));

    1182:   /* allocate and set pointer for implementation data */
    1183:   PetscCall(PetscNew(&s));
    1184:   V->data         = (void *)s;
    1185:   s->setup_called = PETSC_FALSE;
    1186:   s->nb           = -1;
    1187:   s->v            = NULL;

    1189:   PetscCall(VecSetUp_Nest_Private(V, nb, x));

    1191:   n = N = 0;
    1192:   PetscCall(VecSize_Nest_Recursive(V, PETSC_TRUE, &N));
    1193:   PetscCall(VecSize_Nest_Recursive(V, PETSC_FALSE, &n));
    1194:   PetscCall(PetscLayoutSetSize(V->map, N));
    1195:   PetscCall(PetscLayoutSetLocalSize(V->map, n));
    1196:   PetscCall(PetscLayoutSetBlockSize(V->map, 1));
    1197:   PetscCall(PetscLayoutSetUp(V->map));

    1199:   PetscCall(VecSetUp_NestIS_Private(V, nb, is));

    1201:   PetscCall(VecNestSetOps_Private(V->ops));
    1202:   V->petscnative = PETSC_FALSE;

    1204:   /* expose block api's */
    1205:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVec_C", VecNestGetSubVec_Nest));
    1206:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSubVecs_C", VecNestGetSubVecs_Nest));
    1207:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVec_C", VecNestSetSubVec_Nest));
    1208:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestSetSubVecs_C", VecNestSetSubVecs_Nest));
    1209:   PetscCall(PetscObjectComposeFunction((PetscObject)V, "VecNestGetSize_C", VecNestGetSize_Nest));

    1211:   *Y = V;
    1212:   PetscFunctionReturn(PETSC_SUCCESS);
    1213: }