Actual source code: ex1.c
2: /* Program usage: mpiexec ex1 [-help] [all PETSc options] */
4: static char help[] = "Basic vector routines.\n\n";
6: /*T
7: Concepts: vectors^basic routines;
8: Processors: n
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: */
18: #include petscvec.h
22: int main(int argc,char **argv)
23: {
24: Vec x,y,w; /* vectors */
25: Vec *z; /* array of vectors */
26: PetscReal norm,v,v1,v2;
27: PetscInt n = 20;
29: PetscTruth flg;
30: PetscScalar one = 1.0,two = 2.0,three = 3.0,dots[3],dot;
32: PetscInitialize(&argc,&argv,(char*)0,help);
33: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
35: /*
36: Create a vector, specifying only its global dimension.
37: When using VecCreate(), VecSetSizes() and VecSetFromOptions(), the vector format
38: (currently parallel, shared, or sequential) is determined at runtime. Also, the
39: parallel partitioning of the vector is determined by PETSc at runtime.
41: Routines for creating particular vector types directly are:
42: VecCreateSeq() - uniprocessor vector
43: VecCreateMPI() - distributed vector, where the user can
44: determine the parallel partitioning
45: VecCreateShared() - parallel vector that uses shared memory
46: (available only on the SGI); otherwise,
47: is the same as VecCreateMPI()
49: With VecCreate(), VecSetSizes() and VecSetFromOptions() the option -vec_type mpi or
50: -vec_type shared causes the particular type of vector to be formed.
52: */
54: VecCreate(PETSC_COMM_WORLD,&x);
55: VecSetSizes(x,PETSC_DECIDE,n);
56: VecSetFromOptions(x);
58: /*
59: Duplicate some work vectors (of the same format and
60: partitioning as the initial vector).
61: */
62: VecDuplicate(x,&y);
63: VecDuplicate(x,&w);
64: VecNorm(w,NORM_2,&norm);
67: /*
68: Duplicate more work vectors (of the same format and
69: partitioning as the initial vector). Here we duplicate
70: an array of vectors, which is often more convenient than
71: duplicating individual ones.
72: */
73: VecDuplicateVecs(x,3,&z);
75: /*
76: Set the vectors to entries to a constant value.
77: */
78: VecSet(x,one);
79: VecSet(y,two);
80: VecSet(z[0],one);
81: VecSet(z[1],two);
82: VecSet(z[2],three);
84: /*
85: Demonstrate various basic vector routines.
86: */
87: VecDot(x,x,&dot);
88: VecMDot(x,3,z,dots);
90: /*
91: Note: If using a complex numbers version of PETSc, then
92: PETSC_USE_COMPLEX is defined in the makefiles; otherwise,
93: (when using real numbers) it is undefined.
94: */
95: #if defined(PETSC_USE_COMPLEX)
96: PetscPrintf(PETSC_COMM_WORLD,"Vector length %D\n",(PetscInt) (PetscRealPart(dot)));
97: PetscPrintf(PETSC_COMM_WORLD,"Vector length %D %D %D\n",(PetscInt)PetscRealPart(dots[0]),
98: (PetscInt)PetscRealPart(dots[1]),(PetscInt)PetscRealPart(dots[2]));
99: #else
100: PetscPrintf(PETSC_COMM_WORLD,"Vector length %D\n",(PetscInt)dot);
101: PetscPrintf(PETSC_COMM_WORLD,"Vector length %D %D %D\n",(PetscInt)dots[0],
102: (PetscInt)dots[1],(PetscInt)dots[2]);
103: #endif
105: PetscPrintf(PETSC_COMM_WORLD,"All other values should be near zero\n");
107: VecScale(x,two);
108: VecNorm(x,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,"VecScale %G\n",v);
113: VecCopy(x,w);
114: VecNorm(w,NORM_2,&norm);
115: v = norm-2.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
116: PetscPrintf(PETSC_COMM_WORLD,"VecCopy %G\n",v);
118: VecAXPY(y,three,x);
119: VecNorm(y,NORM_2,&norm);
120: v = norm-8.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
121: PetscPrintf(PETSC_COMM_WORLD,"VecAXPY %G\n",v);
123: VecAYPX(y,two,x);
124: VecNorm(y,NORM_2,&norm);
125: v = norm-18.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
126: PetscPrintf(PETSC_COMM_WORLD,"VecAYPX %G\n",v);
128: VecSwap(x,y);
129: VecNorm(y,NORM_2,&norm);
130: v = norm-2.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
131: PetscPrintf(PETSC_COMM_WORLD,"VecSwap %G\n",v);
132: VecNorm(x,NORM_2,&norm);
133: v = norm-18.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
134: PetscPrintf(PETSC_COMM_WORLD,"VecSwap %G\n",v);
136: VecWAXPY(w,two,x,y);
137: VecNorm(w,NORM_2,&norm);
138: v = norm-38.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
139: PetscPrintf(PETSC_COMM_WORLD,"VecWAXPY %G\n",v);
141: VecPointwiseMult(w,y,x);
142: VecNorm(w,NORM_2,&norm);
143: v = norm-36.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
144: PetscPrintf(PETSC_COMM_WORLD,"VecPointwiseMult %G\n",v);
146: VecPointwiseDivide(w,x,y);
147: VecNorm(w,NORM_2,&norm);
148: v = norm-9.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
149: PetscPrintf(PETSC_COMM_WORLD,"VecPointwiseDivide %G\n",v);
151: dots[0] = one;
152: dots[1] = three;
153: dots[2] = two;
154: VecSet(x,one);
155: VecMAXPY(x,3,dots,z);
156: VecNorm(z[0],NORM_2,&norm);
157: v = norm-sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
158: VecNorm(z[1],NORM_2,&norm);
159: v1 = norm-2.0*sqrt((double)n); if (v1 > -PETSC_SMALL && v1 < PETSC_SMALL) v1 = 0.0;
160: VecNorm(z[2],NORM_2,&norm);
161: v2 = norm-3.0*sqrt((double)n); if (v2 > -PETSC_SMALL && v2 < PETSC_SMALL) v2 = 0.0;
162: PetscPrintf(PETSC_COMM_WORLD,"VecMAXPY %G %G %G \n",v,v1,v2);
164: /*
165: Test whether vector has been corrupted (just to demonstrate this
166: routine) not needed in most application codes.
167: */
168: VecValid(x,&flg);
169: if (!flg) SETERRQ(1,"Corrupted vector.");
171: /*
172: Free work space. All PETSc objects should be destroyed when they
173: are no longer needed.
174: */
175: VecDestroy(x);
176: VecDestroy(y);
177: VecDestroy(w);
178: VecDestroyVecs(z,3);
179: PetscFinalize();
180: return 0;
181: }
182: