Actual source code: ex1.c
1: /* Program usage: mpiexec ex1 [-help] [all PETSc options] */
3: static char help[] = "Basic vector routines.\n\n";
5: /*T
6: Concepts: vectors^basic routines;
7: Processors: n
8: T*/
10: /*
11: Include "petscvec.h" so that we can use vectors. Note that this file
12: automatically includes:
13: petscsys.h - base PETSc routines petscis.h - index sets
14: petscviewer.h - viewers
15: */
17: #include <petscvec.h>
21: int main(int argc,char **argv)
22: {
23: Vec x,y,w; /* vectors */
24: Vec *z; /* array of vectors */
25: PetscReal norm,v,v1,v2,maxval;
26: PetscInt n = 20,maxind;
28: PetscScalar one = 1.0,two = 2.0,three = 3.0,dots[3],dot;
30: PetscInitialize(&argc,&argv,(char*)0,help);
31: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
33: /*
34: Create a vector, specifying only its global dimension.
35: When using VecCreate(), VecSetSizes() and VecSetFromOptions(), the vector format
36: (currently parallel, shared, or sequential) is determined at runtime. Also, the
37: parallel partitioning of the vector is determined by PETSc at runtime.
39: Routines for creating particular vector types directly are:
40: VecCreateSeq() - uniprocessor vector
41: VecCreateMPI() - distributed vector, where the user can
42: determine the parallel partitioning
43: VecCreateShared() - parallel vector that uses shared memory
44: (available only on the SGI); otherwise,
45: is the same as VecCreateMPI()
47: With VecCreate(), VecSetSizes() and VecSetFromOptions() the option -vec_type mpi or
48: -vec_type shared causes the particular type of vector to be formed.
49: y
51: */
53: VecCreate(PETSC_COMM_WORLD,&x);
54: VecSetSizes(x,PETSC_DECIDE,n);
55: VecSetFromOptions(x);
56: /*
57: Duplicate some work vectors (of the same format and
58: partitioning as the initial vector).
59: */
60: VecDuplicate(x,&y);
61: VecDuplicate(x,&w);
62: VecNorm(w,NORM_2,&norm);
63: /*
64: Duplicate more work vectors (of the same format and
65: partitioning as the initial vector). Here we duplicate
66: an array of vectors, which is often more convenient than
67: duplicating individual ones.
68: */
69: VecDuplicateVecs(x,3,&z);
70: /*
71: Set the vectors to entries to a constant value.
72: */
73: VecSet(x,one);
74: VecSet(y,two);
75: VecSet(z[0],one);
76: VecSet(z[1],two);
77: VecSet(z[2],three);
78: /*
79: Demonstrate various basic vector routines.
80: */
81: VecDot(x,x,&dot);
82: VecMDot(x,3,z,dots);
84: /*
85: Note: If using a complex numbers version of PETSc, then
86: PETSC_USE_COMPLEX is defined in the makefiles; otherwise,
87: (when using real numbers) it is undefined.
88: */
89: //printf("Size Of Void* = %ld, Size Of PetscErrorCode = %ld\n",sizeof(void*),sizeof(PetscErrorCode));
90: //printf("Size of PetscBool = %ld\n",sizeof(PetscBool));
92: #if defined(PETSC_USE_COMPLEX)
93: PetscPrintf(PETSC_COMM_WORLD,"Vector length %D\n",(PetscInt) (PetscRealPart(dot)));
94: PetscPrintf(PETSC_COMM_WORLD,"Vector length %D %D %D\n",(PetscInt)PetscRealPart(dots[0]),
95: (PetscInt)PetscRealPart(dots[1]),(PetscInt)PetscRealPart(dots[2]));
96: #else
97: PetscPrintf(PETSC_COMM_WORLD,"Vector length %D\n",(PetscInt)dot);
98: PetscPrintf(PETSC_COMM_WORLD,"Vector length %D %D %D\n",(PetscInt)dots[0],
99: (PetscInt)dots[1],(PetscInt)dots[2]);
100: #endif
101: //struct timeval t1, t2, result;
102: //gettimeofday(&t1,NULL);
103: VecMax(x,&maxind,&maxval);
104: //gettimeofday(&t2,NULL);
105: //timersub(&t2,&t1,&result);
106: //printf("Time for VecMax = %lf\n",result.tv_sec+result.tv_usec/1000000.0);
107: PetscPrintf(PETSC_COMM_WORLD,"VecMax %g, VecInd %D\n",(double)maxval,maxind);
109: VecMin(x,&maxind,&maxval);
110: PetscPrintf(PETSC_COMM_WORLD,"VecMin %g, VecInd %D\n",(double)maxval,maxind);
111: /*
112: const PetscInt ix_kds[] = {3,7,14};
113: const PetscScalar y_kds[] = {13.2,69.3,-8.7};
114: VecSetValues(x,3,ix_kds,y_kds,INSERT_VALUES);
115: VecMax(x,&maxind,&maxval);
116: PetscPrintf(PETSC_COMM_WORLD,"VecMax %g, VecInd %D\n",maxval,maxind);
117: VecMin(x,&maxind,&maxval);
118: PetscPrintf(PETSC_COMM_WORLD,"VecMax %g, VecInd %D\n",maxval,maxind);
119: */
120: PetscPrintf(PETSC_COMM_WORLD,"All other values should be near zero\n");
123: VecScale(x,two);
124: VecNorm(x,NORM_2,&norm);
125: v = norm-2.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
126: PetscPrintf(PETSC_COMM_WORLD,"VecScale %g\n",(double)v);
129: VecCopy(x,w);
130: VecNorm(w,NORM_2,&norm);
131: v = norm-2.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
132: PetscPrintf(PETSC_COMM_WORLD,"VecCopy %g\n",(double)v);
134: VecAXPY(y,three,x);
135: VecNorm(y,NORM_2,&norm);
136: v = norm-8.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
137: PetscPrintf(PETSC_COMM_WORLD,"VecAXPY %g\n",(double)v);
139: VecAYPX(y,two,x);
140: VecNorm(y,NORM_2,&norm);
141: v = norm-18.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
142: PetscPrintf(PETSC_COMM_WORLD,"VecAYPX %g\n",(double)v);
144: VecSwap(x,y);
145: VecNorm(y,NORM_2,&norm);
146: v = norm-2.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
147: PetscPrintf(PETSC_COMM_WORLD,"VecSwap %g\n",(double)v);
148: VecNorm(x,NORM_2,&norm);
149: v = norm-18.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
150: PetscPrintf(PETSC_COMM_WORLD,"VecSwap %g\n",(double)v);
152: VecWAXPY(w,two,x,y);
153: VecNorm(w,NORM_2,&norm);
154: v = norm-38.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
155: PetscPrintf(PETSC_COMM_WORLD,"VecWAXPY %g\n",(double)v);
157: VecPointwiseMult(w,y,x);
158: VecNorm(w,NORM_2,&norm);
159: v = norm-36.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
160: PetscPrintf(PETSC_COMM_WORLD,"VecPointwiseMult %g\n",(double)v);
162: VecPointwiseDivide(w,x,y);
163: VecNorm(w,NORM_2,&norm);
164: v = norm-9.0*sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
165: PetscPrintf(PETSC_COMM_WORLD,"VecPointwiseDivide %g\n",(double)v);
167: dots[0] = one;
168: dots[1] = three;
169: dots[2] = two;
170: VecSet(x,one);
171: VecMAXPY(x,3,dots,z);
172: VecNorm(z[0],NORM_2,&norm);
173: v = norm-sqrt((double)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
174: VecNorm(z[1],NORM_2,&norm);
175: v1 = norm-2.0*sqrt((double)n); if (v1 > -PETSC_SMALL && v1 < PETSC_SMALL) v1 = 0.0;
176: VecNorm(z[2],NORM_2,&norm);
177: v2 = norm-3.0*sqrt((double)n); if (v2 > -PETSC_SMALL && v2 < PETSC_SMALL) v2 = 0.0;
178: PetscPrintf(PETSC_COMM_WORLD,"VecMAXPY %g %g %g \n",(double)v,(double)v1,(double)v2);
179:
181: /*
182: Free work space. All PETSc objects should be destroyed when they
183: are no longer needed.
184: */
185: VecDestroy(&x);
186: VecDestroy(&y);
187: VecDestroy(&w);
188: VecDestroyVecs(3,&z);
189: PetscFinalize();
190: return 0;
191: }
192: