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: /*T
  5:    Concepts: vectors^assembling vectors;
  6:    Processors: n

  8:    Contributed by Mike McCourt <mccomic@iit.edu> and Nathan Johnston <johnnat@iit.edu>
  9: T*/

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

 19: PetscScalar func(PetscScalar a)
 20: {
 21:   return 2*a/(1+a*a);
 22: }

 26: int main(int argc,char **argv)
 27: {
 29:   PetscMPIInt    rank,size;
 30:   PetscInt       rstart,rend,i,k,N,numPoints=1000000;
 31:   PetscScalar    dummy,result=0,h=1.0/numPoints,*xarray;
 32:   Vec            x,xend;

 34:   PetscInitialize(&argc,&argv,(char*)0,help);
 35:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 36:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

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

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

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

 95:   /*
 96:       Return the value of the integral.
 97:   */
 98:   PetscPrintf(PETSC_COMM_WORLD,"ln(2) is %g\n",(double)result);
 99:   VecDestroy(&x);
100:   VecDestroy(&xend);

102:   PetscFinalize();
103:   return 0;
104: }