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;

 28:   PetscInitialize(&argc,&argv,(char*)0,help);
 29:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 30:   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:   VecCreate(PETSC_COMM_WORLD,&x);
 39:   VecSetSizes(x,PETSC_DECIDE,numPoints);
 40:   VecSetFromOptions(x);
 41:   VecGetSize(x,&N);
 42:   VecSet(x,result);
 43:   VecDuplicate(x,&xend);
 44:   result = 0.5;
 45:   if (rank == 0) {
 46:     i    = 0;
 47:     VecSetValues(xend,1,&i,&result,INSERT_VALUES);
 48:   }
 49:   if (rank == size-1) {
 50:     i    = N-1;
 51:     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:   VecAssemblyBegin(xend);
 60:   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:   VecGetOwnershipRange(x,&rstart,&rend);
 69:   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:   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:   VecSum(x,&result);
 85:   result = result*h;
 86:   VecDot(x,xend,&dummy);
 87:   result = result-h*dummy;

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

 96:   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*/