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, NULL, 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*/