Actual source code: ex44.c

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

  3: #include <petscvec.h>

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

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

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

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

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

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

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

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

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

153: /*TEST

155:     test:
156:         suffix: serial

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

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

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

172: TEST*/