Actual source code: ex4.c

petsc-3.13.4 2020-08-01
Report Typos and Errors
  1: static char help[]= "Test PetscSFFCompose when the ilocal array is not the identity\n\n";

  3:  #include <petsc.h>
  4:  #include <petscsf.h>

  6: int main(int argc, char **argv)
  7: {
  8:   PetscErrorCode     ierr;
  9:   PetscSF            sfA, sfB, sfBA;
 10:   PetscInt           nrootsA, nleavesA, nrootsB, nleavesB;
 11:   PetscInt          *ilocalA, *ilocalB;
 12:   PetscSFNode       *iremoteA, *iremoteB;
 13:   Vec                a, b, ba;
 14:   const PetscScalar *arrayR;
 15:   PetscScalar       *arrayW;
 16:   PetscMPIInt        size;
 17:   PetscInt           i;
 18:   PetscInt           maxleafB;
 19:   PetscBool          flag = PETSC_FALSE;

 21:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
 22:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 24:   if (size != 1) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER, "Only coded for one MPI process");

 26:   PetscSFCreate(PETSC_COMM_WORLD, &sfA); 
 27:   PetscSFCreate(PETSC_COMM_WORLD, &sfB); 
 28:   PetscSFSetFromOptions(sfA);
 29:   PetscSFSetFromOptions(sfB);

 31:   PetscOptionsGetBool(NULL,NULL,"-sparse_sfB",&flag,NULL);

 33:   if (flag) {
 34:     /* sfA permutes indices, sfB has sparse leaf space. */
 35:     nrootsA = 3;
 36:     nleavesA = 3;
 37:     nrootsB = 3;
 38:     nleavesB = 2;
 39:   } else {
 40:     /* sfA reverses indices, sfB is identity */
 41:     nrootsA = nrootsB = nleavesA = nleavesB = 4;
 42:   }
 43:   PetscMalloc1(nleavesA, &ilocalA); 
 44:   PetscMalloc1(nleavesA, &iremoteA); 
 45:   PetscMalloc1(nleavesB, &ilocalB); 
 46:   PetscMalloc1(nleavesB, &iremoteB); 

 48:   for (i = 0; i < nleavesA; i++) {
 49:     iremoteA[i].rank = 0;
 50:     iremoteA[i].index = i;
 51:     if (flag) {
 52:       ilocalA[i] = (i + 1) % nleavesA;
 53:     } else {
 54:       ilocalA[i] = nleavesA - i - 1;
 55:     }
 56:   }

 58:   for (i = 0; i < nleavesB; i++) {
 59:     iremoteB[i].rank = 0;
 60:     if (flag) {
 61:       ilocalB[i] = nleavesB - i;
 62:       iremoteB[i].index = nleavesB - i - 1;
 63:     } else {
 64:       ilocalB[i] = i;
 65:       iremoteB[i].index = i;
 66:     }
 67:   }

 69:   PetscSFSetGraph(sfA, nrootsA, nleavesA, ilocalA, PETSC_OWN_POINTER, iremoteA, PETSC_OWN_POINTER); 
 70:   PetscSFSetGraph(sfB, nrootsB, nleavesB, ilocalB, PETSC_OWN_POINTER, iremoteB, PETSC_OWN_POINTER); 
 71:   PetscSFSetUp(sfA); 
 72:   PetscSFSetUp(sfB); 
 73:   PetscObjectSetName((PetscObject)sfA, "sfA");
 74:   PetscObjectSetName((PetscObject)sfB, "sfB");

 76:   VecCreateSeq(PETSC_COMM_WORLD, nrootsA, &a); 
 77:   VecCreateSeq(PETSC_COMM_WORLD, nleavesA, &b); 
 78:   PetscSFGetLeafRange(sfB, NULL, &maxleafB); 
 79:   VecCreateSeq(PETSC_COMM_WORLD, maxleafB+1, &ba); 
 80:   VecGetArray(a, &arrayW); 
 81:   for ( i = 0; i < nrootsA; i++ ) {
 82:     arrayW[i] = (PetscScalar)i;
 83:   }
 84:   VecRestoreArray(a, &arrayW); 

 86:   PetscPrintf(PETSC_COMM_WORLD, "Initial Vec A\n"); 
 87:   VecView(a, NULL); 
 88:   VecGetArrayRead(a, &arrayR); 
 89:   VecGetArray(b, &arrayW); 

 91:   PetscSFBcastBegin(sfA, MPIU_SCALAR, arrayR, arrayW); 
 92:   PetscSFBcastEnd(sfA, MPIU_SCALAR, arrayR, arrayW); 
 93:   VecRestoreArray(b, &arrayW); 
 94:   VecRestoreArrayRead(a, &arrayR); 
 95:   PetscPrintf(PETSC_COMM_WORLD, "\nBroadcast A->B over sfA\n"); 
 96:   VecView(b, NULL); 

 98:   VecGetArrayRead(b, &arrayR); 
 99:   VecGetArray(ba, &arrayW); 
100:   arrayW[0] = 10.0;             /* Not touched by bcast */
101:   PetscSFBcastBegin(sfB, MPIU_SCALAR, arrayR, arrayW); 
102:   PetscSFBcastEnd(sfB, MPIU_SCALAR, arrayR, arrayW); 
103:   VecRestoreArrayRead(b, &arrayR); 
104:   VecRestoreArray(ba, &arrayW); 

106:   PetscPrintf(PETSC_COMM_WORLD, "\nBroadcast B->BA over sfB\n"); 
107:   VecView(ba, NULL); 

109:   PetscSFCompose(sfA, sfB, &sfBA); 
110:   PetscSFSetFromOptions(sfBA);
111:   PetscObjectSetName((PetscObject)sfBA, "(sfB o sfA)");
112:   VecGetArrayRead(a, &arrayR); 
113:   VecGetArray(ba, &arrayW); 
114:   arrayW[0] = 11.0;             /* Not touched by bcast */
115:   PetscSFBcastBegin(sfBA, MPIU_SCALAR, arrayR, arrayW); 
116:   PetscSFBcastEnd(sfBA, MPIU_SCALAR, arrayR, arrayW); 
117:   VecRestoreArray(ba, &arrayW); 
118:   VecRestoreArrayRead(a, &arrayR); 
119:   PetscPrintf(PETSC_COMM_WORLD, "\nBroadcast A->BA over sfBA (sfB o sfA)\n"); 
120:   VecView(ba, NULL); 

122:   VecDestroy(&ba); 
123:   VecDestroy(&b); 
124:   VecDestroy(&a); 

126:   PetscSFView(sfA, NULL); 
127:   PetscSFView(sfB, NULL); 
128:   PetscSFView(sfBA, NULL); 
129:   PetscSFDestroy(&sfA); 
130:   PetscSFDestroy(&sfB); 
131:   PetscSFDestroy(&sfBA); 

133:   PetscFinalize();

135:   return ierr;
136: }

138: /*TEST

140:    test:
141:      suffix: 1

143:    test:
144:      suffix: 2
145:      filter: grep -v "type" | grep -v "sort"
146:      args: -sparse_sfB

148:    test:
149:      suffix: 2_window
150:      filter: grep -v "type" | grep -v "sort"
151:      output_file: output/ex4_2.out
152:      args: -sparse_sfB -sf_type window -sf_window_sync {{fence active lock}} -sf_window_flavor {{create dynamic allocate}}
153:      requires: define(PETSC_HAVE_MPI_ONE_SIDED)

155:    # The nightly test suite with MPICH uses ch3:sock, which is broken when winsize == 0 in some of the processes
156:    test:
157:      suffix: 2_window_shared
158:      filter: grep -v "type" | grep -v "sort"
159:      output_file: output/ex4_2.out
160:      args: -sparse_sfB -sf_type window -sf_window_sync {{fence active lock}} -sf_window_flavor shared
161:      requires: define(PETSC_HAVE_MPI_PROCESS_SHARED_MEMORY) !define(PETSC_HAVE_MPICH_NUMVERSION) define(PETSC_HAVE_MPI_ONE_SIDED)

163: TEST*/