Actual source code: ex1.c

      1: static char help[] = "VecTagger interface routines.\n\n";

      3: #include <petscvec.h>

      5: static PetscErrorCode ISGetBlockGlobalIS(IS is, Vec vec, PetscInt bs, IS *isBlockGlobal)
      6: {
      7:   const PetscInt *idxin;
      8:   PetscInt       *idxout, i, n, rstart;
      9:   PetscLayout     map;

     11:   PetscFunctionBegin;
     12:   PetscCall(VecGetLayout(vec, &map));
     13:   rstart = map->rstart / bs;
     14:   PetscCall(ISGetLocalSize(is, &n));
     15:   PetscCall(PetscMalloc1(n, &idxout));
     16:   PetscCall(ISGetIndices(is, &idxin));
     17:   for (i = 0; i < n; i++) idxout[i] = rstart + idxin[i];
     18:   PetscCall(ISRestoreIndices(is, &idxin));
     19:   PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)vec), bs, n, idxout, PETSC_OWN_POINTER, isBlockGlobal));
     20:   PetscFunctionReturn(PETSC_SUCCESS);
     21: }

     23: int main(int argc, char **argv)
     24: {
     25:   Vec           vec, tagged, untagged;
     26:   VecScatter    taggedScatter, untaggedScatter;
     27:   PetscInt      bs;
     28:   PetscInt      n, nloc, nint, i, j, k, localStart, localEnd, ntagged, nuntagged;
     29:   MPI_Comm      comm;
     30:   VecTagger     tagger;
     31:   PetscScalar  *array;
     32:   PetscRandom   rand;
     33:   VecTaggerBox *defaultBox;
     34:   VecTaggerBox *boxes;
     35:   IS            is, isBlockGlobal, isComp;
     36:   PetscBool     listed;

     38:   PetscFunctionBeginUser;
     39:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     40:   n    = 10.;
     41:   bs   = 1;
     42:   comm = PETSC_COMM_WORLD;
     43:   PetscOptionsBegin(comm, "", "VecTagger Test Options", "Vec");
     44:   PetscCall(PetscOptionsInt("-bs", "The block size of the vector", "ex1.c", bs, &bs, NULL));
     45:   PetscCall(PetscOptionsInt("-n", "The size of the vector (in blocks)", "ex1.c", n, &n, NULL));
     46:   PetscOptionsEnd();

     48:   PetscCall(PetscRandomCreate(comm, &rand));
     49:   PetscCall(PetscRandomSetFromOptions(rand));

     51:   PetscCall(VecCreate(comm, &vec));
     52:   PetscCall(PetscObjectSetName((PetscObject)vec, "Vec to Tag"));
     53:   PetscCall(VecSetBlockSize(vec, bs));
     54:   PetscCall(VecSetSizes(vec, PETSC_DECIDE, n));
     55:   PetscCall(VecSetUp(vec));
     56:   PetscCall(VecGetLocalSize(vec, &nloc));
     57:   PetscCall(VecGetArray(vec, &array));
     58:   for (i = 0; i < nloc; i++) {
     59:     PetscScalar val;

     61:     PetscCall(PetscRandomGetValue(rand, &val));
     62:     array[i] = val;
     63:   }
     64:   PetscCall(VecRestoreArray(vec, &array));
     65:   PetscCall(PetscRandomDestroy(&rand));
     66:   PetscCall(VecViewFromOptions(vec, NULL, "-vec_view"));

     68:   PetscCall(VecTaggerCreate(comm, &tagger));
     69:   PetscCall(VecTaggerSetBlockSize(tagger, bs));
     70:   PetscCall(VecTaggerSetType(tagger, VECTAGGERABSOLUTE));
     71:   PetscCall(PetscMalloc1(bs, &defaultBox));
     72:   for (i = 0; i < bs; i++) {
     73: #if !defined(PETSC_USE_COMPLEX)
     74:     defaultBox[i].min = 0.1;
     75:     defaultBox[i].max = 1.5;
     76: #else
     77:     defaultBox[i].min = PetscCMPLX(0.1, 0.1);
     78:     defaultBox[i].max = PetscCMPLX(1.5, 1.5);
     79: #endif
     80:   }
     81:   PetscCall(VecTaggerAbsoluteSetBox(tagger, defaultBox));
     82:   PetscCall(PetscFree(defaultBox));
     83:   PetscCall(VecTaggerSetFromOptions(tagger));
     84:   PetscCall(VecTaggerSetUp(tagger));
     85:   PetscCall(PetscObjectViewFromOptions((PetscObject)tagger, NULL, "-vec_tagger_view"));
     86:   PetscCall(VecTaggerGetBlockSize(tagger, &bs));

     88:   PetscCall(VecTaggerComputeBoxes(tagger, vec, &nint, &boxes, &listed));
     89:   if (listed) {
     90:     PetscViewer viewer = NULL;

     92:     PetscCall(PetscOptionsCreateViewer(comm, NULL, NULL, "-vec_tagger_boxes_view", &viewer, NULL, NULL));
     93:     if (viewer) {
     94:       PetscBool isascii;

     96:       PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
     97:       if (isascii) {
     98:         PetscCall(PetscViewerASCIIPrintf(viewer, "Num boxes: %" PetscInt_FMT "\n", nint));
     99:         PetscCall(PetscViewerASCIIPushTab(viewer));
    100:         for (i = 0, k = 0; i < nint; i++) {
    101:           PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT ": ", i));
    102:           for (j = 0; j < bs; j++, k++) {
    103:             if (j) PetscCall(PetscViewerASCIIPrintf(viewer, " x "));
    104: #if !defined(PETSC_USE_COMPLEX)
    105:             PetscCall(PetscViewerASCIIPrintf(viewer, "[%g,%g]", (double)boxes[k].min, (double)boxes[k].max));
    106: #else
    107:             PetscCall(PetscViewerASCIIPrintf(viewer, "[%g+%gi,%g+%gi]", (double)PetscRealPart(boxes[k].min), (double)PetscImaginaryPart(boxes[k].min), (double)PetscRealPart(boxes[k].max), (double)PetscImaginaryPart(boxes[k].max)));
    108: #endif
    109:           }
    110:           PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
    111:         }
    112:         PetscCall(PetscViewerASCIIPopTab(viewer));
    113:       }
    114:     }
    115:     PetscCall(PetscViewerDestroy(&viewer));
    116:     PetscCall(PetscFree(boxes));
    117:   }

    119:   PetscCall(VecTaggerComputeIS(tagger, vec, &is, &listed));
    120:   PetscCall(ISGetBlockGlobalIS(is, vec, bs, &isBlockGlobal));
    121:   PetscCall(PetscObjectSetName((PetscObject)isBlockGlobal, "Tagged IS (block global)"));
    122:   PetscCall(ISViewFromOptions(isBlockGlobal, NULL, "-tagged_is_view"));

    124:   PetscCall(VecGetOwnershipRange(vec, &localStart, &localEnd));
    125:   PetscCall(ISComplement(isBlockGlobal, localStart, localEnd, &isComp));
    126:   PetscCall(PetscObjectSetName((PetscObject)isComp, "Untagged IS (global)"));
    127:   PetscCall(ISViewFromOptions(isComp, NULL, "-untagged_is_view"));

    129:   PetscCall(ISGetLocalSize(isBlockGlobal, &ntagged));
    130:   PetscCall(ISGetLocalSize(isComp, &nuntagged));

    132:   PetscCall(VecCreate(comm, &tagged));
    133:   PetscCall(PetscObjectSetName((PetscObject)tagged, "Tagged selection"));
    134:   PetscCall(VecSetSizes(tagged, ntagged, PETSC_DETERMINE));
    135:   PetscCall(VecSetUp(tagged));

    137:   PetscCall(VecCreate(comm, &untagged));
    138:   PetscCall(PetscObjectSetName((PetscObject)untagged, "Untagged selection"));
    139:   PetscCall(VecSetSizes(untagged, nuntagged, PETSC_DETERMINE));
    140:   PetscCall(VecSetUp(untagged));

    142:   PetscCall(VecScatterCreate(vec, isBlockGlobal, tagged, NULL, &taggedScatter));
    143:   PetscCall(VecScatterCreate(vec, isComp, untagged, NULL, &untaggedScatter));

    145:   PetscCall(VecScatterBegin(taggedScatter, vec, tagged, INSERT_VALUES, SCATTER_FORWARD));
    146:   PetscCall(VecScatterEnd(taggedScatter, vec, tagged, INSERT_VALUES, SCATTER_FORWARD));
    147:   PetscCall(VecScatterBegin(untaggedScatter, vec, untagged, INSERT_VALUES, SCATTER_FORWARD));
    148:   PetscCall(VecScatterEnd(untaggedScatter, vec, untagged, INSERT_VALUES, SCATTER_FORWARD));

    150:   PetscCall(VecViewFromOptions(tagged, NULL, "-tagged_vec_view"));
    151:   PetscCall(VecViewFromOptions(untagged, NULL, "-untagged_vec_view"));

    153:   PetscCall(VecScatterDestroy(&untaggedScatter));
    154:   PetscCall(VecScatterDestroy(&taggedScatter));

    156:   PetscCall(VecDestroy(&untagged));
    157:   PetscCall(VecDestroy(&tagged));
    158:   PetscCall(ISDestroy(&isComp));
    159:   PetscCall(ISDestroy(&isBlockGlobal));
    160:   PetscCall(ISDestroy(&is));

    162:   PetscCall(VecTaggerDestroy(&tagger));
    163:   PetscCall(VecDestroy(&vec));
    164:   PetscCall(PetscFinalize());
    165:   return 0;
    166: }

    168: /*TEST

    170:   test:
    171:     suffix: 0
    172:     requires: !complex
    173:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view

    175:   test:
    176:     suffix: 1
    177:     requires: !complex
    178:     nsize: 3
    179:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view

    181:   test:
    182:     suffix: 2
    183:     requires: !complex
    184:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -bs 2

    186:   test:
    187:     suffix: 3
    188:     requires: !complex
    189:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_block_size 2 -vec_tagger_box 0.1,1.5,0.1,1.5

    191:   test:
    192:     suffix: 4
    193:     requires: !complex
    194:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_invert

    196:   test:
    197:     suffix: 5
    198:     requires: !complex
    199:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type relative -vec_tagger_box 0.25,0.75

    201:   test:
    202:     suffix: 6
    203:     requires: !complex
    204:     nsize: 3
    205:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type cdf -vec_tagger_box 0.25,0.75

    207:   test:
    208:     suffix: 7
    209:     requires: !complex
    210:     nsize: 3
    211:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type cdf -vec_tagger_box 0.25,0.75 -vec_tagger_cdf_method iterative -vec_tagger_cdf_max_it 10

    213:   test:
    214:     suffix: 8
    215:     requires: !complex
    216:     nsize: 3
    217:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type or -vec_tagger_num_subs 2 -sub_0_vec_tagger_type absolute -sub_0_vec_tagger_box 0.0,0.25 -sub_1_vec_tagger_type relative -sub_1_vec_tagger_box 0.75,inf
    218:     filter: sed -e s~Inf~inf~g

    220:   test:
    221:     suffix: 9
    222:     requires: !complex
    223:     nsize: 3
    224:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type and -vec_tagger_num_subs 2 -sub_0_vec_tagger_type absolute -sub_0_vec_tagger_box -inf,0.5 -sub_1_vec_tagger_type relative -sub_1_vec_tagger_box 0.25,0.75
    225:     filter: sed -e s~Inf~inf~g

    227:   test:
    228:     suffix: 10
    229:     requires: complex
    230:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view

    232:   test:
    233:     suffix: 11
    234:     requires: complex
    235:     nsize: 3
    236:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view

    238:   test:
    239:     suffix: 12
    240:     requires: complex
    241:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -bs 2

    243:   test:
    244:     suffix: 13
    245:     requires: complex
    246:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_block_size 2 -vec_tagger_box 0.1+0.1i,1.5+1.5i,0.1+0.1i,1.5+1.5i

    248:   test:
    249:     suffix: 14
    250:     requires: complex
    251:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_invert

    253:   test:
    254:     suffix: 15
    255:     requires: complex
    256:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type relative -vec_tagger_box 0.25+0.25i,0.75+0.75i

    258:   test:
    259:     suffix: 16
    260:     requires: complex
    261:     nsize: 3
    262:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type cdf -vec_tagger_box 0.25+0.25i,0.75+0.75i

    264:   test:
    265:     suffix: 17
    266:     requires: complex
    267:     nsize: 3
    268:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type cdf -vec_tagger_box 0.25+0.25i,0.75+0.75i -vec_tagger_cdf_method iterative -vec_tagger_cdf_max_it 10

    270:   test:
    271:     suffix: 18
    272:     requires: complex
    273:     nsize: 3
    274:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type or -vec_tagger_num_subs 2 -sub_0_vec_tagger_type absolute -sub_0_vec_tagger_box 0.0+0.0i,0.25+0.25i -sub_1_vec_tagger_type relative -sub_1_vec_tagger_box 0.75+0.75i,inf+infi
    275:     filter: sed -e s~Inf~inf~g

    277:   test:
    278:     suffix: 19
    279:     requires: complex
    280:     nsize: 3
    281:     args: -n 12 -vec_view -vec_tagger_view -vec_tagger_boxes_view -tagged_is_view -untagged_is_view -tagged_vec_view -untagged_vec_view -vec_tagger_type and -vec_tagger_num_subs 2 -sub_0_vec_tagger_type absolute -sub_0_vec_tagger_box -inf-infi,0.75+0.75i -sub_1_vec_tagger_type relative -sub_1_vec_tagger_box 0.25+0.25i,1.+1.i
    282:     filter: sed -e s~Inf~inf~g

    284: TEST*/