Actual source code: ex18.c

  1: static char help[] = "Computes the integral of 2*x/(1+x^2) from x=0..1 \nThis is equal to the ln(2).\n\n";

  3: /*
  4:    Contributed by Mike McCourt <mccomic@iit.edu> and Nathan Johnston <johnnat@iit.edu>
  5: */

  7: /*
  8:   Include "petscvec.h" so that we can use vectors.  Note that this file
  9:   automatically includes:
 10:      petscsys.h       - base PETSc routines   petscis.h     - index sets
 11:      petscviewer.h - viewers
 12: */
 13: #include <petscvec.h>

 15: PetscScalar func(PetscScalar a)
 16: {
 17:   return (PetscScalar)2. * a / ((PetscScalar)1. + a * a);
 18: }

 20: int main(int argc, char **argv)
 21: {
 22:   PetscMPIInt rank, size;
 23:   PetscInt    rstart, rend, i, k, N, numPoints = 1000000;
 24:   PetscScalar dummy, result = 0, h = 1.0 / numPoints, *xarray;
 25:   Vec         x, xend;

 27:   PetscFunctionBeginUser;
 28:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 29:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 30:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

 32:   /*
 33:      Create a parallel vector.
 34:        Here we set up our x vector which will be given values below.
 35:        The xend vector is a dummy vector to find the value of the
 36:          elements at the endpoints for use in the trapezoid rule.
 37:   */
 38:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 39:   PetscCall(VecSetSizes(x, PETSC_DECIDE, numPoints));
 40:   PetscCall(VecSetFromOptions(x));
 41:   PetscCall(VecGetSize(x, &N));
 42:   PetscCall(VecSet(x, result));
 43:   PetscCall(VecDuplicate(x, &xend));
 44:   result = 0.5;
 45:   if (rank == 0) {
 46:     i = 0;
 47:     PetscCall(VecSetValues(xend, 1, &i, &result, INSERT_VALUES));
 48:   }
 49:   if (rank == size - 1) {
 50:     i = N - 1;
 51:     PetscCall(VecSetValues(xend, 1, &i, &result, INSERT_VALUES));
 52:   }
 53:   /*
 54:      Assemble vector, using the 2-step process:
 55:        VecAssemblyBegin(), VecAssemblyEnd()
 56:      Computations can be done while messages are in transition
 57:      by placing code between these two statements.
 58:   */
 59:   PetscCall(VecAssemblyBegin(xend));
 60:   PetscCall(VecAssemblyEnd(xend));

 62:   /*
 63:      Set the x vector elements.
 64:       i*h will return 0 for i=0 and 1 for i=N-1.
 65:       The function evaluated (2x/(1+x^2)) is defined above.
 66:       Each evaluation is put into the local array of the vector without message passing.
 67:   */
 68:   PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
 69:   PetscCall(VecGetArray(x, &xarray));
 70:   k = 0;
 71:   for (i = rstart; i < rend; i++) {
 72:     xarray[k] = (PetscScalar)i * h;
 73:     xarray[k] = func(xarray[k]);
 74:     k++;
 75:   }
 76:   PetscCall(VecRestoreArray(x, &xarray));

 78:   /*
 79:      Evaluates the integral.  First the sum of all the points is taken.
 80:      That result is multiplied by the step size for the trapezoid rule.
 81:      Then half the value at each endpoint is subtracted,
 82:      this is part of the composite trapezoid rule.
 83:   */
 84:   PetscCall(VecSum(x, &result));
 85:   result = result * h;
 86:   PetscCall(VecDot(x, xend, &dummy));
 87:   result = result - h * dummy;

 89:   /*
 90:       Return the value of the integral.
 91:   */
 92:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ln(2) is %g\n", (double)PetscRealPart(result)));
 93:   PetscCall(VecDestroy(&x));
 94:   PetscCall(VecDestroy(&xend));

 96:   PetscCall(PetscFinalize());
 97:   return 0;
 98: }

100: /*TEST

102:      test:
103:        nsize: 1

105:      test:
106:        nsize: 2
107:        suffix: 2
108:        output_file: output/ex18_1.out

110: TEST*/