Actual source code: ex18.c
1: 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";
3: /*
4: Contributed by Mike McCourt <mccomic@iit.edu> and Nathan Johnston <johnnat@iit.edu>
5: */
7: /*
8: Include "petscvec.h" so that we can use vectors. Note that this file
9: automatically includes:
10: petscsys.h - base PETSc routines petscis.h - index sets
11: petscviewer.h - viewers
12: */
13: #include <petscvec.h>
15: PetscScalar func(PetscScalar a)
16: {
17: return (PetscScalar)2. * a / ((PetscScalar)1. + a * a);
18: }
20: int main(int argc, char **argv)
21: {
22: PetscMPIInt rank, size;
23: PetscInt rstart, rend, i, k, N, numPoints = 1000000;
24: PetscScalar dummy, result = 0, h = 1.0 / numPoints, *xarray;
25: Vec x, xend;
27: PetscFunctionBeginUser;
28: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
29: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
30: PetscCallMPI(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: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
39: PetscCall(VecSetSizes(x, PETSC_DECIDE, numPoints));
40: PetscCall(VecSetFromOptions(x));
41: PetscCall(VecGetSize(x, &N));
42: PetscCall(VecSet(x, result));
43: PetscCall(VecDuplicate(x, &xend));
44: result = 0.5;
45: if (rank == 0) {
46: i = 0;
47: PetscCall(VecSetValues(xend, 1, &i, &result, INSERT_VALUES));
48: }
49: if (rank == size - 1) {
50: i = N - 1;
51: PetscCall(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: PetscCall(VecAssemblyBegin(xend));
60: PetscCall(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: PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
69: PetscCall(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: PetscCall(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: PetscCall(VecSum(x, &result));
85: result = result * h;
86: PetscCall(VecDot(x, xend, &dummy));
87: result = result - h * dummy;
89: /*
90: Return the value of the integral.
91: */
92: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "ln(2) is %g\n", (double)PetscRealPart(result)));
93: PetscCall(VecDestroy(&x));
94: PetscCall(VecDestroy(&xend));
96: PetscCall(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*/