Actual source code: ex1.c
1: /*$Id: ex1.c,v 1.74 2001/09/19 16:08:15 bsmith Exp $*/
3: /* Program usage: mpirun ex1 [-help] [all PETSc options] */
5: static char help[] = "Basic vector routines.nn";
7: /*T
8: Concepts: vectors^basic routines;
9: Processors: n
10: T*/
12: /*
13: Include "petscvec.h" so that we can use vectors. Note that this file
14: automatically includes:
15: petsc.h - base PETSc routines petscis.h - index sets
16: petscsys.h - system routines petscviewer.h - viewers
17: */
19: #include petscvec.h
21: #undef __FUNCT__
23: int main(int argc,char **argv)
24: {
25: Vec x,y,w; /* vectors */
26: Vec *z; /* array of vectors */
27: PetscReal norm,v,v1,v2;
28: int n = 20,ierr;
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(&one,x);
79: VecSet(&two,y);
80: VecSet(&one,z[0]);
81: VecSet(&two,z[1]);
82: VecSet(&three,z[2]);
84: /*
85: Demonstrate various basic vector routines.
86: */
87: VecDot(x,x,&dot);
88: VecMDot(3,x,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 %dn",int (PetscRealPart(dot)));
97: PetscPrintf(PETSC_COMM_WORLD,"Vector length %d %d %dn",(int)PetscRealPart(dots[0]),
98: (int)PetscRealPart(dots[1]),(int)PetscRealPart(dots[2]));
99: #else
100: PetscPrintf(PETSC_COMM_WORLD,"Vector length %dn",(int)dot);
101: PetscPrintf(PETSC_COMM_WORLD,"Vector length %d %d %dn",(int)dots[0],
102: (int)dots[1],(int)dots[2]);
103: #endif
105: PetscPrintf(PETSC_COMM_WORLD,"All other values should be near zeron");
107: VecScale(&two,x);
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 %gn",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 %gn",v);
118: VecAXPY(&three,x,y);
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 %gn",v);
123: VecAYPX(&two,x,y);
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 %gn",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 %gn",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 %gn",v);
136: VecWAXPY(&two,x,y,w);
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 %gn",v);
141: VecPointwiseMult(y,x,w);
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 %gn",v);
146: VecPointwiseDivide(x,y,w);
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 %gn",v);
151: dots[0] = one;
152: dots[1] = three;
153: dots[2] = two;
154: VecSet(&one,x);
155: VecMAXPY(3,dots,x,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: