Actual source code: ex41.c

  1: static char help[] = "Nest vector set subvector functionality.\n\n";

  3: #include <petscvec.h>

  5: PetscErrorCode test_vec_ops(void)
  6: {
  7:   Vec         X, Y, a, b;
  8:   Vec         c, d, e, f, g, h;
  9:   PetscScalar val;
 10:   PetscInt    tmp_ind[2];
 11:   Vec         tmp_buf[2];

 13:   PetscFunctionBegin;
 14:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "============== %s ==============\n", PETSC_FUNCTION_NAME));

 16:   /* create 4 worker vectors */
 17:   PetscCall(VecCreate(PETSC_COMM_WORLD, &c));
 18:   PetscCall(VecSetSizes(c, PETSC_DECIDE, 4));
 19:   PetscCall(VecSetType(c, VECMPI));
 20:   PetscCall(VecDuplicate(c, &d));
 21:   PetscCall(VecDuplicate(c, &e));
 22:   PetscCall(VecDuplicate(c, &f));

 24:   /* create two more workers of different sizes */
 25:   PetscCall(VecCreate(PETSC_COMM_WORLD, &g));
 26:   PetscCall(VecSetSizes(g, PETSC_DECIDE, 6));
 27:   PetscCall(VecSetType(g, VECMPI));
 28:   PetscCall(VecCreate(PETSC_COMM_WORLD, &h));
 29:   PetscCall(VecSetSizes(h, PETSC_DECIDE, 8));
 30:   PetscCall(VecSetType(h, VECMPI));

 32:   /* set the 6 vectors to some numbers */
 33:   PetscCall(VecSet(c, 1.0));
 34:   PetscCall(VecSet(d, 2.0));
 35:   PetscCall(VecSet(e, 3.0));
 36:   PetscCall(VecSet(f, 4.0));
 37:   PetscCall(VecSet(g, 5.0));
 38:   PetscCall(VecSet(h, 6.0));

 40:   /* assemble a */
 41:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "a = [c d] \n"));
 42:   tmp_buf[0] = c;
 43:   tmp_buf[1] = d;

 45:   PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &a));
 46:   PetscCall(VecView(a, PETSC_VIEWER_STDOUT_WORLD));
 47:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "a = [d c] \n"));
 48:   PetscCall(VecNestSetSubVec(a, 1, c));
 49:   PetscCall(VecNestSetSubVec(a, 0, d));
 50:   PetscCall(VecAssemblyBegin(a));
 51:   PetscCall(VecAssemblyEnd(a));
 52:   PetscCall(VecView(a, PETSC_VIEWER_STDOUT_WORLD));

 54:   /* assemble b */
 55:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "b = [e f] \n"));
 56:   tmp_buf[0] = e;
 57:   tmp_buf[1] = f;

 59:   PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &b));
 60:   PetscCall(VecView(b, PETSC_VIEWER_STDOUT_WORLD));
 61:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "b = [f e] \n"));
 62:   PetscCall(VecNestSetSubVec(b, 1, e));
 63:   PetscCall(VecNestSetSubVec(b, 0, f));
 64:   PetscCall(VecAssemblyBegin(b));
 65:   PetscCall(VecAssemblyEnd(b));
 66:   PetscCall(VecView(b, PETSC_VIEWER_STDOUT_WORLD));

 68:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X = [a b] \n"));
 69:   tmp_buf[0] = a;
 70:   tmp_buf[1] = b;

 72:   PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &X));
 73:   PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
 74:   PetscCall(VecDot(X, X, &val));
 75:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.X = %g \n", (double)PetscRealPart(val)));

 77:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X = [b a] \n"));
 78:   /* re-order components of X */
 79:   PetscCall(VecNestSetSubVec(X, 1, a));
 80:   PetscCall(VecNestSetSubVec(X, 0, b));
 81:   PetscCall(VecAssemblyBegin(X));
 82:   PetscCall(VecAssemblyEnd(X));
 83:   PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
 84:   PetscCall(VecDot(X, X, &val));
 85:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.X = %g \n", (double)PetscRealPart(val)));

 87:   /* re-assemble X */
 88:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X = [g h] \n"));
 89:   PetscCall(VecNestSetSubVec(X, 1, g));
 90:   PetscCall(VecNestSetSubVec(X, 0, h));
 91:   PetscCall(VecAssemblyBegin(X));
 92:   PetscCall(VecAssemblyEnd(X));
 93:   PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
 94:   PetscCall(VecDot(X, X, &val));
 95:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.X = %g \n", (double)PetscRealPart(val)));

 97:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Y = X \n"));
 98:   PetscCall(VecDuplicate(X, &Y));
 99:   PetscCall(VecCopy(X, Y));
100:   PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));
101:   PetscCall(VecDot(Y, Y, &val));
102:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Y.Y = %g \n", (double)PetscRealPart(val)));

104:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Y = [a b] \n"));
105:   tmp_buf[0] = a;
106:   tmp_buf[1] = b;
107:   tmp_ind[0] = 0;
108:   tmp_ind[1] = 1;

110:   PetscCall(VecNestSetSubVecs(Y, 2, tmp_ind, tmp_buf));
111:   PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));

113:   PetscCall(VecDestroy(&c));
114:   PetscCall(VecDestroy(&d));
115:   PetscCall(VecDestroy(&e));
116:   PetscCall(VecDestroy(&f));
117:   PetscCall(VecDestroy(&g));
118:   PetscCall(VecDestroy(&h));
119:   PetscCall(VecDestroy(&a));
120:   PetscCall(VecDestroy(&b));
121:   PetscCall(VecDestroy(&X));
122:   PetscCall(VecDestroy(&Y));
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

126: int main(int argc, char **args)
127: {
128:   PetscFunctionBeginUser;
129:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
130:   PetscCall(test_vec_ops());
131:   PetscCall(PetscFinalize());
132:   return 0;
133: }

135: /*TEST

137:    test:

139: TEST*/