Actual source code: ex44.c


  2: static char help[] = "Test VecConcatenate both in serial and parallel.\n";

  4: #include <petscvec.h>

  6: int main(int argc, char **args)
  7: {
  8:   Vec         *x, x_test, y, y_test;
  9:   IS          *x_is;
 10:   VecScatter   y_to_x, x_to_y;
 11:   PetscInt     i, j, nx, shift, x_size, y_size, *idx;
 12:   PetscScalar *vals;
 13:   PetscBool    flg, x_equal, y_equal;

 15:   PetscFunctionBeginUser;
 16:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 17:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nx", &nx, &flg));
 18:   if (!flg) nx = 3;

 20:   y_size = 0;
 21:   shift  = 0;
 22:   PetscCall(PetscMalloc1(nx, &x));
 23:   for (i = 0; i < nx; i++) {
 24:     x_size = 2 * (i + 1);
 25:     y_size += x_size;
 26:     PetscCall(VecCreate(PETSC_COMM_WORLD, &x[i]));
 27:     PetscCall(VecSetSizes(x[i], PETSC_DECIDE, x_size));
 28:     PetscCall(VecSetFromOptions(x[i]));
 29:     PetscCall(VecSetUp(x[i]));
 30:     PetscCall(PetscMalloc2(x_size, &idx, x_size, &vals));
 31:     for (j = 0; j < x_size; j++) {
 32:       idx[j]  = j;
 33:       vals[j] = (PetscScalar)(shift + j + 1);
 34:     }
 35:     shift += x_size;
 36:     PetscCall(VecSetValues(x[i], x_size, (const PetscInt *)idx, (const PetscScalar *)vals, INSERT_VALUES));
 37:     PetscCall(VecAssemblyBegin(x[i]));
 38:     PetscCall(VecAssemblyEnd(x[i]));
 39:     PetscCall(PetscFree2(idx, vals));
 40:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Original X[%" PetscInt_FMT "] vector\n", i));
 41:     PetscCall(VecView(x[i], PETSC_VIEWER_STDOUT_WORLD));
 42:   }
 43:   PetscCall(VecCreate(PETSC_COMM_WORLD, &y));
 44:   PetscCall(VecSetSizes(y, PETSC_DECIDE, y_size));
 45:   PetscCall(VecSetFromOptions(y));
 46:   PetscCall(VecSetUp(y));
 47:   PetscCall(PetscMalloc2(y_size, &idx, y_size, &vals));
 48:   for (j = 0; j < y_size; j++) {
 49:     idx[j]  = j;
 50:     vals[j] = (PetscScalar)(j + 1);
 51:   }
 52:   PetscCall(VecSetValues(y, y_size, (const PetscInt *)idx, (const PetscScalar *)vals, INSERT_VALUES));
 53:   PetscCall(VecAssemblyBegin(y));
 54:   PetscCall(VecAssemblyEnd(y));
 55:   PetscCall(PetscFree2(idx, vals));
 56:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected Y vector\n"));
 57:   PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD));
 58:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "------------------------------------------------------------\n"));

 60:   /* ---------- base VecConcatenate() test ----------- */
 61:   PetscCall(VecConcatenate(nx, (const Vec *)x, &y_test, &x_is));
 62:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing VecConcatenate() for Y = [X[1], X[2], ...]\n"));
 63:   PetscCall(VecView(y_test, PETSC_VIEWER_STDOUT_WORLD));
 64:   y_equal = PETSC_FALSE;
 65:   PetscCall(VecEqual(y_test, y, &y_equal));
 66:   if (!y_equal) {
 67:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  FAIL\n"));
 68:   } else {
 69:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  PASS\n"));
 70:   }
 71:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "------------------------------------------------------------\n"));

 73:   /* ---------- test VecConcatenate() without IS (checks for dangling memory from IS) ----------- */
 74:   PetscCall(VecDestroy(&y_test));
 75:   PetscCall(VecConcatenate(nx, (const Vec *)x, &y_test, NULL));
 76:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing VecConcatenate() for Y = [X[1], X[2], ...] w/o IS\n"));
 77:   PetscCall(VecView(y_test, PETSC_VIEWER_STDOUT_WORLD));
 78:   y_equal = PETSC_FALSE;
 79:   PetscCall(VecEqual(y_test, y, &y_equal));
 80:   if (!y_equal) {
 81:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  FAIL\n"));
 82:   } else {
 83:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  PASS\n"));
 84:   }
 85:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "------------------------------------------------------------\n"));

 87:   /* ---------- using index sets on expected Y instead of concatenated Y ----------- */
 88:   for (i = 0; i < nx; i++) {
 89:     PetscCall(VecGetSubVector(y, x_is[i], &x_test));
 90:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing index set for X[%" PetscInt_FMT "] component\n", i));
 91:     PetscCall(VecView(x_test, PETSC_VIEWER_STDOUT_WORLD));
 92:     x_equal = PETSC_FALSE;
 93:     PetscCall(VecEqual(x_test, x[i], &x_equal));
 94:     if (!x_equal) {
 95:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  FAIL\n"));
 96:     } else {
 97:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  PASS\n"));
 98:     }
 99:     PetscCall(VecRestoreSubVector(y, x_is[i], &x_test));
100:   }
101:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "------------------------------------------------------------\n"));

103:   /* ---------- using VecScatter to communicate data from Y to X[i] for all i ----------- */
104:   for (i = 0; i < nx; i++) {
105:     PetscCall(VecDuplicate(x[i], &x_test));
106:     PetscCall(VecZeroEntries(x_test));
107:     PetscCall(VecScatterCreate(y, x_is[i], x[i], NULL, &y_to_x));
108:     PetscCall(VecScatterBegin(y_to_x, y, x_test, INSERT_VALUES, SCATTER_FORWARD));
109:     PetscCall(VecScatterEnd(y_to_x, y, x_test, INSERT_VALUES, SCATTER_FORWARD));
110:     PetscCall(VecScatterDestroy(&y_to_x));
111:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing VecScatter for Y -> X[%" PetscInt_FMT "]\n", i));
112:     PetscCall(VecView(x_test, PETSC_VIEWER_STDOUT_WORLD));
113:     x_equal = PETSC_FALSE;
114:     PetscCall(VecEqual(x_test, x[i], &x_equal));
115:     if (!x_equal) {
116:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  FAIL\n"));
117:     } else {
118:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  PASS\n"));
119:     }
120:     PetscCall(VecDestroy(&x_test));
121:   }
122:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "------------------------------------------------------------\n"));

124:   /* ---------- using VecScatter to communicate data from X[i] to Y for all i ----------- */
125:   PetscCall(VecZeroEntries(y_test));
126:   for (i = 0; i < nx; i++) {
127:     PetscCall(VecScatterCreate(x[i], NULL, y, x_is[i], &x_to_y));
128:     PetscCall(VecScatterBegin(x_to_y, x[i], y_test, INSERT_VALUES, SCATTER_FORWARD));
129:     PetscCall(VecScatterEnd(x_to_y, x[i], y_test, INSERT_VALUES, SCATTER_FORWARD));
130:     PetscCall(VecScatterDestroy(&x_to_y));
131:   }
132:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing VecScatter for X[:] -> Y\n"));
133:   PetscCall(VecView(y_test, PETSC_VIEWER_STDOUT_WORLD));
134:   y_equal = PETSC_FALSE;
135:   PetscCall(VecEqual(y_test, y, &y_equal));
136:   if (!y_equal) {
137:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  FAIL\n"));
138:   } else {
139:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  PASS\n"));
140:   }
141:   PetscCall(VecDestroy(&y_test));

143:   PetscCall(VecDestroy(&y));
144:   for (i = 0; i < nx; i++) {
145:     PetscCall(VecDestroy(&x[i]));
146:     PetscCall(ISDestroy(&x_is[i]));
147:   }
148:   PetscCall(PetscFree(x));
149:   PetscCall(PetscFree(x_is));
150:   PetscCall(PetscFinalize());
151:   return 0;
152: }

154: /*TEST

156:     test:
157:         suffix: serial

159:     test:
160:         suffix: parallel
161:         nsize: 2

163:     test:
164:         suffix: cuda
165:         nsize: 2
166:         args: -vec_type cuda
167:         requires: cuda

169:     test:
170:         suffix: uneven
171:         nsize: 5

173: TEST*/