2: 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";

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

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

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

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

 29:   PetscInitialize(&argc, &argv, (char *)0, help);
 30:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 31:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

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

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

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

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

 97:   PetscFinalize();
 98:   return 0;
 99: }

101: /*TEST

103:      test:
104:        nsize: 1

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

111: TEST*/