Actual source code: ex39.c

  1: static char help[] = "This example is intended for showing how subvectors can\n\
  2:                       share the pointer with the main vector using VecGetArray()\n\
  3:                       and VecPlaceArray() routines so that vector operations done\n\
  4:                       on the subvectors automatically modify the values in the main vector.\n\n";

  6: #include <petscvec.h>

  8: /* This example shares the array pointers of vectors X,Y,and F with subvectors
  9:    X1,X2,Y1,Y2,F1,F2 and does vector addition on the subvectors F1 = X1 + Y1, F2 = X2 + Y2 so
 10:    that F gets updated as a result of sharing the pointers.
 11:  */

 13: int main(int argc, char **argv)
 14: {
 15:   PetscInt           N = 10, i;
 16:   Vec                X, Y, F, X1, Y1, X2, Y2, F1, F2;
 17:   PetscScalar        value, zero = 0.0;
 18:   const PetscScalar *x, *y;
 19:   PetscScalar       *f;

 21:   PetscFunctionBeginUser;
 22:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

 24:   /* create vectors X,Y and F and set values in it*/
 25:   PetscCall(VecCreate(PETSC_COMM_SELF, &X));
 26:   PetscCall(VecSetSizes(X, N, N));
 27:   PetscCall(VecSetFromOptions(X));
 28:   PetscCall(VecDuplicate(X, &Y));
 29:   PetscCall(VecDuplicate(X, &F));
 30:   PetscCall(PetscObjectSetName((PetscObject)F, "F"));
 31:   for (i = 0; i < N; i++) {
 32:     value = i;
 33:     PetscCall(VecSetValues(X, 1, &i, &value, INSERT_VALUES));
 34:     value = 100 + i;
 35:     PetscCall(VecSetValues(Y, 1, &i, &value, INSERT_VALUES));
 36:   }
 37:   PetscCall(VecSet(F, zero));

 39:   /* Create subvectors X1,X2,Y1,Y2,F1,F2 */
 40:   PetscCall(VecCreate(PETSC_COMM_SELF, &X1));
 41:   PetscCall(VecSetSizes(X1, N / 2, N / 2));
 42:   PetscCall(VecSetFromOptions(X1));
 43:   PetscCall(VecDuplicate(X1, &X2));
 44:   PetscCall(VecDuplicate(X1, &Y1));
 45:   PetscCall(VecDuplicate(X1, &Y2));
 46:   PetscCall(VecDuplicate(X1, &F1));
 47:   PetscCall(VecDuplicate(X1, &F2));

 49:   /* Get array pointers for X,Y,F */
 50:   PetscCall(VecGetArrayRead(X, &x));
 51:   PetscCall(VecGetArrayRead(Y, &y));
 52:   PetscCall(VecGetArray(F, &f));
 53:   /* Share X,Y,F array pointers with subvectors */
 54:   PetscCall(VecPlaceArray(X1, x));
 55:   PetscCall(VecPlaceArray(X2, x + N / 2));
 56:   PetscCall(VecPlaceArray(Y1, y));
 57:   PetscCall(VecPlaceArray(Y2, y + N / 2));
 58:   PetscCall(VecPlaceArray(F1, f));
 59:   PetscCall(VecPlaceArray(F2, f + N / 2));

 61:   /* Do subvector addition */
 62:   PetscCall(VecWAXPY(F1, 1.0, X1, Y1));
 63:   PetscCall(VecWAXPY(F2, 1.0, X2, Y2));

 65:   /* Reset subvectors */
 66:   PetscCall(VecResetArray(X1));
 67:   PetscCall(VecResetArray(X2));
 68:   PetscCall(VecResetArray(Y1));
 69:   PetscCall(VecResetArray(Y2));
 70:   PetscCall(VecResetArray(F1));
 71:   PetscCall(VecResetArray(F2));

 73:   /* Restore X,Y,and F */
 74:   PetscCall(VecRestoreArrayRead(X, &x));
 75:   PetscCall(VecRestoreArrayRead(Y, &y));
 76:   PetscCall(VecRestoreArray(F, &f));

 78:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "F = X + Y\n"));
 79:   PetscCall(VecView(F, 0));
 80:   /* Destroy vectors */
 81:   PetscCall(VecDestroy(&X));
 82:   PetscCall(VecDestroy(&Y));
 83:   PetscCall(VecDestroy(&F));
 84:   PetscCall(VecDestroy(&X1));
 85:   PetscCall(VecDestroy(&Y1));
 86:   PetscCall(VecDestroy(&F1));
 87:   PetscCall(VecDestroy(&X2));
 88:   PetscCall(VecDestroy(&Y2));
 89:   PetscCall(VecDestroy(&F2));

 91:   PetscCall(PetscFinalize());
 92:   return 0;
 93: }