Actual source code: ex18.c
petsc-3.3-p7 2013-05-11
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){return 2*a/(1+a*a);}
23: int main(int argc,char **argv)
24: {
26: PetscMPIInt rank,nproc;
27: PetscInt rstart,rend,i,k,N,numPoints=1000000;
28: PetscScalar dummy,result=0,h=1.0/numPoints,*xarray;
29: Vec x,xend;
31: PetscInitialize(&argc,&argv,(char *)0,help);
32: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
33: MPI_Comm_size(PETSC_COMM_WORLD,&nproc);
35: /*
36: Create a parallel vector.
37: Here we set up our x vector which will be given values below.
38: The xend vector is a dummy vector to find the value of the
39: elements at the endpoints for use in the trapezoid rule.
40: */
41: VecCreate(PETSC_COMM_WORLD,&x);
42: VecSetSizes(x,PETSC_DECIDE,numPoints);
43: VecSetFromOptions(x);
44: VecGetSize(x,&N);
45: VecSet(x,result);
46: VecDuplicate(x,&xend);
47: result=0.5;
48: if (rank == 0){
49: i=0;
50: VecSetValues(xend,1,&i,&result,INSERT_VALUES);
51: } else if (rank == nproc){
52: i=N-1;
53: VecSetValues(xend,1,&i,&result,INSERT_VALUES);
54: }
55: /*
56: Assemble vector, using the 2-step process:
57: VecAssemblyBegin(), VecAssemblyEnd()
58: Computations can be done while messages are in transition
59: by placing code between these two statements.
60: */
61: VecAssemblyBegin(xend);
62: VecAssemblyEnd(xend);
64: /*
65: Set the x vector elements.
66: i*h will return 0 for i=0 and 1 for i=N-1.
67: The function evaluated (2x/(1+x^2)) is defined above.
68: Each evaluation is put into the local array of the vector without message passing.
69: */
70: VecGetOwnershipRange(x,&rstart,&rend);
71: VecGetArray(x,&xarray);
72: k = 0;
73: for (i=rstart; i<rend; i++){
74: xarray[k] = i*h;
75: xarray[k] = func(xarray[k]);
76: k++;
77: }
78: VecRestoreArray(x,&xarray);
80: /*
81: Evaluates the integral. First the sum of all the points is taken.
82: That result is multiplied by the step size for the trapezoid rule.
83: Then half the value at each endpoint is subtracted,
84: this is part of the composite trapezoid rule.
85: */
86: VecSum(x,&result);
87: result = result*h;
88: VecDot(x,xend,&dummy);
89: result = result-h*dummy;
91: /*
92: Return the value of the integral.
93: */
94: PetscPrintf(PETSC_COMM_WORLD,"ln(2) is %G\n",result);
95: VecDestroy(&x);
96: VecDestroy(&xend);
98: PetscFinalize();
99: return 0;
100: }
101: