Actual source code: ex1e.c
2: /* Program usage: mpiexec ex1 [-help] [all PETSc options] */
4: static char help[] = "Demonstrates various vector routines.\n\n";
6: /*T
7: Concepts: vectors^basic routines;
8: Processors: n
9: T*/
11: /*
13: This uses the PETSc _ error checking routines. Put _ before the PETSc function call
14: and __ after the call (or ___ in a subroutine, not the main program). This is equivalent
15: to using the ... macros
18: Include "petscvec.h" so that we can use vectors. Note that this file
19: automatically includes:
20: petscsys.h - base PETSc routines petscis.h - index sets
21: petscviewer.h - viewers
22: */
23: #define PETSC_UNDERSCORE_CHKERR
25: #include petscvec.h
27: #if defined(PETSC_USE_SCALAR_SINGLE)
28: #define PETSC_EPS 1.e-5
29: #else
30: #define PETSC_EPS 1.e-10
31: #endif
35: int main(int argc,char **argv)
36: {
37: Vec x, y, w; /* vectors */
38: Vec *z; /* array of vectors */
39: PetscReal norm, v, v1, v2;
40: PetscInt n = 20;
41: PetscTruth flg;
42: PetscScalar one = 1.0, two = 2.0, three = 3.0, dots[3], dot;
44: _ PetscInitialize(&argc,&argv,(char*)0,help);___
45: _ PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);___
47: /*
48: Create a vector, specifying only its global dimension.
49: When using VecCreate(), VecSetSizes() and VecSetFromOptions(), the vector format
50: (currently parallel, shared, or sequential) is determined at runtime. Also, the
51: parallel partitioning of the vector is determined by PETSc at runtime.
53: Routines for creating particular vector types directly are:
54: VecCreateSeq() - uniprocessor vector
55: VecCreateMPI() - distributed vector, where the user can
56: determine the parallel partitioning
57: VecCreateShared() - parallel vector that uses shared memory
58: (available only on the SGI); otherwise,
59: is the same as VecCreateMPI()
61: With VecCreate(), VecSetSizes() and VecSetFromOptions() the option -vec_type mpi or
62: -vec_type shared causes the particular type of vector to be formed.
64: */
65: _ VecCreate(PETSC_COMM_WORLD,&x);___
66: _ VecSetSizes(x,PETSC_DECIDE,n);___
67: _ VecSetFromOptions(x);___
69: /*
70: Duplicate some work vectors (of the same format and
71: partitioning as the initial vector).
72: */
73: _ VecDuplicate(x,&y);___
74: _ VecDuplicate(x,&w);___
76: /*
77: Duplicate more work vectors (of the same format and
78: partitioning as the initial vector). Here we duplicate
79: an array of vectors, which is often more convenient than
80: duplicating individual ones.
81: */
82: _ VecDuplicateVecs(x,3,&z);___
84: /*
85: Set the vectors to entries to a constant value.
86: */
87: _ VecSet(x,one);___
88: _ VecSet(y,two);___
89: _ VecSet(z[0],one);___
90: _ VecSet(z[1],two);___
91: _ VecSet(z[2],three);___
93: /*
94: Demonstrate various basic vector routines.
95: */
96: _ VecDot(x,x,&dot);___
97: _ VecMDot(x,3,z,dots);___
99: /*
100: Note: If using a complex numbers version of PETSc, then
101: PETSC_USE_COMPLEX is defined in the makefiles; otherwise,
102: (when using real numbers) it is undefined.
103: */
104: #if defined(PETSC_USE_COMPLEX)
105: _ PetscPrintf(PETSC_COMM_WORLD,"Vector length %D\n", (int) (PetscRealPart(dot)));___
106: _ PetscPrintf(PETSC_COMM_WORLD,"Vector length %D %D %D\n",(PetscInt)PetscRealPart(dots[0]),
107: (PetscInt)PetscRealPart(dots[1]),(PetscInt)PetscRealPart(dots[2]));___
108: #else
109: _ PetscPrintf(PETSC_COMM_WORLD,"Vector length %D\n",(PetscInt) dot);___
110: _ PetscPrintf(PETSC_COMM_WORLD,"Vector length %D %D %D\n",(PetscInt)dots[0],
111: (PetscInt)dots[1],(PetscInt)dots[2]);___
112: #endif
114: _ PetscPrintf(PETSC_COMM_WORLD,"All other values should be near zero\n");___
116: _ VecScale(x,two);___
117: _ VecNorm(x,NORM_2,&norm);___
118: v = norm-2.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
119: _ PetscPrintf(PETSC_COMM_WORLD,"VecScale %G\n",v);___
121: _ VecCopy(x,w);___
122: _ VecNorm(w,NORM_2,&norm);___
123: v = norm-2.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
124: _ PetscPrintf(PETSC_COMM_WORLD,"VecCopy %G\n",v);___
126: _ VecAXPY(y,three,x);___
127: _ VecNorm(y,NORM_2,&norm);___
128: v = norm-8.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
129: _ PetscPrintf(PETSC_COMM_WORLD,"VecAXPY %G\n",v);___
131: _ VecAYPX(y,two,x);___
132: _ VecNorm(y,NORM_2,&norm);___
133: v = norm-18.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
134: _ PetscPrintf(PETSC_COMM_WORLD,"VecAYPX %G\n",v);___
136: _ VecSwap(x,y);___
137: _ VecNorm(y,NORM_2,&norm);___
138: v = norm-2.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
139: _ PetscPrintf(PETSC_COMM_WORLD,"VecSwap %G\n",v);___
140: _ VecNorm(x,NORM_2,&norm);___
141: v = norm-18.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
142: _ PetscPrintf(PETSC_COMM_WORLD,"VecSwap %G\n",v);___
144: _ VecWAXPY(w,two,x,y);___
145: _ VecNorm(w,NORM_2,&norm);___
146: v = norm-38.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
147: _ PetscPrintf(PETSC_COMM_WORLD,"VecWAXPY %G\n",v);___
149: _ VecPointwiseMult(w,y,x);___
150: _ VecNorm(w,NORM_2,&norm);___
151: v = norm-36.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
152: _ PetscPrintf(PETSC_COMM_WORLD,"VecPointwiseMult %G\n",v);___
154: _ VecPointwiseDivide(w,x,y);___
155: _ VecNorm(w,NORM_2,&norm);___
156: v = norm-9.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
157: _ PetscPrintf(PETSC_COMM_WORLD,"VecPointwiseDivide %G\n",v);___
159: dots[0] = one;
160: dots[1] = three;
161: dots[2] = two;
162: _ VecSet(x,one);___
163: _ VecMAXPY(x,3,dots,z);___
164: _ VecNorm(z[0],NORM_2,&norm);___
165: v = norm-sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
166: _ VecNorm(z[1],NORM_2,&norm);___
167: v1 = norm-2.0*sqrt((PetscReal) n); if (v1 > -PETSC_EPS && v1 < PETSC_EPS) v1 = 0.0;
168: _ VecNorm(z[2],NORM_2,&norm);___
169: v2 = norm-3.0*sqrt((PetscReal) n); if (v2 > -PETSC_EPS && v2 < PETSC_EPS) v2 = 0.0;
170: _ PetscPrintf(PETSC_COMM_WORLD,"VecMAXPY %G %G %G \n",v,v1,v2);___
172: /*
173: Test whether vector has been corrupted (just to demonstrate this
174: routine) not needed in most application codes.
175: */
176: _ VecValid(x,&flg);___
177: if (!flg) SETERRQ(1,"Corrupted vector.");
179: /*
180: Free work space. All PETSc objects should be destroyed when they
181: are no longer needed.
182: */
183: _ VecDestroy(x);___
184: _ VecDestroy(y);___
185: _ VecDestroy(w);___
186: _ VecDestroyVecs(z,3);___
187: _ PetscFinalize();___
188: return 0;
189: }
190: