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"; /*T Concepts: vectors^assembling vectors; Processors: n Contributed by Mike McCourt and Nathan Johnston T*/ /* Include "petscvec.h" so that we can use vectors. Note that this file automatically includes: petscsys.h - base PETSc routines petscis.h - index sets petscviewer.h - viewers */ #include PetscScalar func(PetscScalar a) { return (PetscScalar)2.*a/((PetscScalar)1.+a*a); } int main(int argc,char **argv) { PetscErrorCode ierr; PetscMPIInt rank,size; PetscInt rstart,rend,i,k,N,numPoints=1000000; PetscScalar dummy,result=0,h=1.0/numPoints,*xarray; Vec x,xend; ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); /* Create a parallel vector. Here we set up our x vector which will be given values below. The xend vector is a dummy vector to find the value of the elements at the endpoints for use in the trapezoid rule. */ ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); ierr = VecSetSizes(x,PETSC_DECIDE,numPoints);CHKERRQ(ierr); ierr = VecSetFromOptions(x);CHKERRQ(ierr); ierr = VecGetSize(x,&N);CHKERRQ(ierr); ierr = VecSet(x,result);CHKERRQ(ierr); ierr = VecDuplicate(x,&xend);CHKERRQ(ierr); result = 0.5; if (!rank) { i = 0; ierr = VecSetValues(xend,1,&i,&result,INSERT_VALUES);CHKERRQ(ierr); } if (rank == size-1) { i = N-1; ierr = VecSetValues(xend,1,&i,&result,INSERT_VALUES);CHKERRQ(ierr); } /* Assemble vector, using the 2-step process: VecAssemblyBegin(), VecAssemblyEnd() Computations can be done while messages are in transition by placing code between these two statements. */ ierr = VecAssemblyBegin(xend);CHKERRQ(ierr); ierr = VecAssemblyEnd(xend);CHKERRQ(ierr); /* Set the x vector elements. i*h will return 0 for i=0 and 1 for i=N-1. The function evaluated (2x/(1+x^2)) is defined above. Each evaluation is put into the local array of the vector without message passing. */ ierr = VecGetOwnershipRange(x,&rstart,&rend);CHKERRQ(ierr); ierr = VecGetArray(x,&xarray);CHKERRQ(ierr); k = 0; for (i=rstart; i