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:   PetscInitialize(&argc,&args,(char*)0,help);
 16:   PetscOptionsGetInt(NULL,NULL,"-nx",&nx,&flg);
 17:   if (!flg) nx = 3;

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

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

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

 86:   /* ---------- using index sets on expected Y instead of concatenated Y ----------- */
 87:   for (i=0; i<nx; i++) {
 88:     VecGetSubVector(y, x_is[i], &x_test);
 89:     PetscPrintf(PETSC_COMM_WORLD, "Testing index set for X[%" PetscInt_FMT "] component\n", i);
 90:     VecView(x_test, PETSC_VIEWER_STDOUT_WORLD);
 91:     x_equal = PETSC_FALSE;
 92:     VecEqual(x_test, x[i], &x_equal);
 93:     if (!x_equal) {
 94:       PetscPrintf(PETSC_COMM_WORLD, "  FAIL\n");
 95:     } else {
 96:       PetscPrintf(PETSC_COMM_WORLD, "  PASS\n");
 97:     }
 98:     VecRestoreSubVector(y, x_is[i], &x_test);
 99:   }
100:   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:     VecDuplicate(x[i], &x_test);
105:     VecZeroEntries(x_test);
106:     VecScatterCreate(y, x_is[i], x[i], NULL, &y_to_x);
107:     VecScatterBegin(y_to_x, y, x_test, INSERT_VALUES, SCATTER_FORWARD);
108:     VecScatterEnd(y_to_x, y, x_test, INSERT_VALUES, SCATTER_FORWARD);
109:     VecScatterDestroy(&y_to_x);
110:     PetscPrintf(PETSC_COMM_WORLD, "Testing VecScatter for Y -> X[%" PetscInt_FMT "]\n", i);
111:     VecView(x_test, PETSC_VIEWER_STDOUT_WORLD);
112:     x_equal = PETSC_FALSE;
113:     VecEqual(x_test, x[i], &x_equal);
114:     if (!x_equal) {
115:       PetscPrintf(PETSC_COMM_WORLD, "  FAIL\n");
116:     } else {
117:       PetscPrintf(PETSC_COMM_WORLD, "  PASS\n");
118:     }
119:     VecDestroy(&x_test);
120:   }
121:   PetscPrintf(PETSC_COMM_WORLD, "------------------------------------------------------------\n");

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

142:   VecDestroy(&y);
143:   for (i=0; i<nx; i++) {
144:     VecDestroy(&x[i]);
145:     ISDestroy(&x_is[i]);
146:   }
147:   PetscFree(x);
148:   PetscFree(x_is);
149:   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*/