Actual source code: ex1.c
petsc-3.4.5 2014-06-29
2: static char help[] = "Basic vector routines.\n\n";
4: /*T
5: Concepts: vectors^basic routines;
6: Processors: n
7: T*/
9: /*
10: Include "petscvec.h" so that we can use vectors. Note that this file
11: automatically includes:
12: petscsys.h - base PETSc routines petscis.h - index sets
13: petscviewer.h - viewers
14: */
16: #include <petscvec.h>
20: int main(int argc,char **argv)
21: {
22: Vec x,y,w; /* vectors */
23: Vec *z; /* array of vectors */
24: PetscReal norm,v,v1,v2,maxval;
25: PetscInt n = 20,maxind;
27: PetscScalar one = 1.0,two = 2.0,three = 3.0,dots[3],dot;
29: PetscInitialize(&argc,&argv,(char*)0,help);
32: PetscOptionsGetInt(NULL,"-n",&n,NULL);
34: /*
35: Create a vector, specifying only its global dimension.
36: When using VecCreate(), VecSetSizes() and VecSetFromOptions(), the vector format
37: (currently parallel, shared, or sequential) is determined at runtime. Also, the
38: parallel partitioning of the vector is determined by PETSc at runtime.
40: Routines for creating particular vector types directly are:
41: VecCreateSeq() - uniprocessor vector
42: VecCreateMPI() - distributed vector, where the user can
43: determine the parallel partitioning
44: VecCreateShared() - parallel vector that uses shared memory
45: (available only on the SGI); otherwise,
46: is the same as VecCreateMPI()
48: With VecCreate(), VecSetSizes() and VecSetFromOptions() the option -vec_type mpi or
49: -vec_type shared causes the particular type of vector to be formed.
50: y
52: */
54: VecCreate(PETSC_COMM_WORLD,&x);
55: VecSetSizes(x,PETSC_DECIDE,n);
56: VecSetFromOptions(x);
57: /*
58: Duplicate some work vectors (of the same format and
59: partitioning as the initial vector).
60: */
61: VecDuplicate(x,&y);
62: VecDuplicate(x,&w);
64: /*
65: Duplicate more work vectors (of the same format and
66: partitioning as the initial vector). Here we duplicate
67: an array of vectors, which is often more convenient than
68: duplicating individual ones.
69: */
70: VecDuplicateVecs(x,3,&z);
71: /*
72: Set the vectors to entries to a constant value.
73: */
74: VecSet(x,one);
75: VecSet(y,two);
76: VecSet(z[0],one);
77: VecSet(z[1],two);
78: VecSet(z[2],three);
79: /*
80: Demonstrate various basic vector routines.
81: */
82: MPI_Barrier(PETSC_COMM_WORLD);
83: VecDot(x,y,&dot);
84: VecMDot(x,3,z,dots);
86: /*
87: Note: If using a complex numbers version of PETSc, then
88: PETSC_USE_COMPLEX is defined in the makefiles; otherwise,
89: (when using real numbers) it is undefined.
90: */
92: PetscPrintf(PETSC_COMM_WORLD,"Vector length %D\n",n);
93: VecMax(x,&maxind,&maxval);
94: PetscPrintf(PETSC_COMM_WORLD,"VecMax %g, VecInd %D\n",(double)maxval,maxind);
96: VecMin(x,&maxind,&maxval);
97: PetscPrintf(PETSC_COMM_WORLD,"VecMin %g, VecInd %D\n",(double)maxval,maxind);
98: PetscPrintf(PETSC_COMM_WORLD,"All other values should be near zero\n");
101: VecScale(x,two);
102: VecNorm(x,NORM_2,&norm);
103: v = norm-2.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
104: PetscPrintf(PETSC_COMM_WORLD,"VecScale %g\n",(double)v);
107: VecCopy(x,w);
108: VecNorm(w,NORM_2,&norm);
109: v = norm-2.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
110: PetscPrintf(PETSC_COMM_WORLD,"VecCopy %g\n",(double)v);
112: VecAXPY(y,three,x);
113: VecNorm(y,NORM_2,&norm);
114: v = norm-8.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
115: PetscPrintf(PETSC_COMM_WORLD,"VecAXPY %g\n",(double)v);
117: VecAYPX(y,two,x);
118: VecNorm(y,NORM_2,&norm);
119: v = norm-18.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
120: PetscPrintf(PETSC_COMM_WORLD,"VecAYPX %g\n",(double)v);
122: VecSwap(x,y);
123: VecNorm(y,NORM_2,&norm);
124: v = norm-2.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
125: PetscPrintf(PETSC_COMM_WORLD,"VecSwap %g\n",(double)v);
126: VecNorm(x,NORM_2,&norm);
127: v = norm-18.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
128: PetscPrintf(PETSC_COMM_WORLD,"VecSwap %g\n",(double)v);
130: VecWAXPY(w,two,x,y);
131: VecNorm(w,NORM_2,&norm);
132: v = norm-38.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
133: PetscPrintf(PETSC_COMM_WORLD,"VecWAXPY %g\n",(double)v);
135: VecPointwiseMult(w,y,x);
136: VecNorm(w,NORM_2,&norm);
137: v = norm-36.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
138: PetscPrintf(PETSC_COMM_WORLD,"VecPointwiseMult %g\n",(double)v);
140: VecPointwiseDivide(w,x,y);
141: VecNorm(w,NORM_2,&norm);
142: v = norm-9.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
143: PetscPrintf(PETSC_COMM_WORLD,"VecPointwiseDivide %g\n",(double)v);
145: dots[0] = one;
146: dots[1] = three;
147: dots[2] = two;
149: VecSet(x,one);
150: VecMAXPY(x,3,dots,z);
151: VecNorm(z[0],NORM_2,&norm);
152: v = norm-sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
153: VecNorm(z[1],NORM_2,&norm);
154: v1 = norm-2.0*sqrt((double)n); if (v1 > -PETSC_SMALL && v1 < PETSC_SMALL) v1 = 0.0;
155: VecNorm(z[2],NORM_2,&norm);
156: v2 = norm-3.0*sqrt((double)n); if (v2 > -PETSC_SMALL && v2 < PETSC_SMALL) v2 = 0.0;
157: PetscPrintf(PETSC_COMM_WORLD,"VecMAXPY %g %g %g \n",(double)v,(double)v1,(double)v2);
159: /*
160: Free work space. All PETSc objects should be destroyed when they
161: are no longer needed.
162: */
163: VecDestroy(&x);
164: VecDestroy(&y);
165: VecDestroy(&w);
166: VecDestroyVecs(3,&z);
167: PetscFinalize();
168: return 0;
169: }