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