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