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_REAL_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:   PetscScalar one = 1.0, two = 2.0, three = 3.0, dots[3], dot;

 43: _ PetscInitialize(&argc,&argv,(char*)0,help);___
 44: _ PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);___

 46:   /* 
 47:      Create a vector, specifying only its global dimension.
 48:      When using VecCreate(), VecSetSizes() and VecSetFromOptions(), the vector format 
 49:      (currently parallel, shared, or sequential) is determined at runtime.  Also, the 
 50:      parallel partitioning of the vector is determined by PETSc at runtime.

 52:      Routines for creating particular vector types directly are:
 53:         VecCreateSeq() - uniprocessor vector
 54:         VecCreateMPI() - distributed vector, where the user can
 55:                          determine the parallel partitioning
 56:         VecCreateShared() - parallel vector that uses shared memory
 57:                             (available only on the SGI); otherwise,
 58:                             is the same as VecCreateMPI()

 60:      With VecCreate(), VecSetSizes() and VecSetFromOptions() the option -vec_type mpi or 
 61:      -vec_type shared causes the particular type of vector to be formed.

 63:   */
 64: _ VecCreate(PETSC_COMM_WORLD,&x);___
 65: _ VecSetSizes(x,PETSC_DECIDE,n);___
 66: _ VecSetFromOptions(x);___

 68:   /*
 69:      Duplicate some work vectors (of the same format and
 70:      partitioning as the initial vector).
 71:   */
 72: _ VecDuplicate(x,&y);___
 73: _ VecDuplicate(x,&w);___

 75:   /*
 76:      Duplicate more work vectors (of the same format and
 77:      partitioning as the initial vector).  Here we duplicate
 78:      an array of vectors, which is often more convenient than
 79:      duplicating individual ones.
 80:   */
 81: _ VecDuplicateVecs(x,3,&z);___

 83:   /*
 84:      Set the vectors to entries to a constant value.
 85:   */
 86: _ VecSet(x,one);___
 87: _ VecSet(y,two);___
 88: _ VecSet(z[0],one);___
 89: _ VecSet(z[1],two);___
 90: _ VecSet(z[2],three);___

 92:   /*
 93:      Demonstrate various basic vector routines.
 94:   */
 95: _ VecDot(x,x,&dot);___
 96: _ VecMDot(x,3,z,dots);___

 98:   /* 
 99:      Note: If using a complex numbers version of PETSc, then
100:      PETSC_USE_COMPLEX is defined in the makefiles; otherwise,
101:      (when using real numbers) it is undefined.
102:   */
103: #if defined(PETSC_USE_COMPLEX)
104: _ PetscPrintf(PETSC_COMM_WORLD,"Vector length %D\n", (int) (PetscRealPart(dot)));___
105: _ PetscPrintf(PETSC_COMM_WORLD,"Vector length %D %D %D\n",(PetscInt)PetscRealPart(dots[0]),
106:                              (PetscInt)PetscRealPart(dots[1]),(PetscInt)PetscRealPart(dots[2]));___
107: #else
108: _ PetscPrintf(PETSC_COMM_WORLD,"Vector length %D\n",(PetscInt) dot);___
109: _ PetscPrintf(PETSC_COMM_WORLD,"Vector length %D %D %D\n",(PetscInt)dots[0],
110:                              (PetscInt)dots[1],(PetscInt)dots[2]);___
111: #endif

113: _ PetscPrintf(PETSC_COMM_WORLD,"All other values should be near zero\n");___

115: _ VecScale(x,two);___
116: _ VecNorm(x,NORM_2,&norm);___
117:   v = norm-2.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
118: _ PetscPrintf(PETSC_COMM_WORLD,"VecScale %G\n",v);___

120: _ VecCopy(x,w);___
121: _ VecNorm(w,NORM_2,&norm);___
122:   v = norm-2.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
123: _ PetscPrintf(PETSC_COMM_WORLD,"VecCopy  %G\n",v);___

125: _ VecAXPY(y,three,x);___
126: _ VecNorm(y,NORM_2,&norm);___
127:   v = norm-8.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
128: _ PetscPrintf(PETSC_COMM_WORLD,"VecAXPY %G\n",v);___

130: _ VecAYPX(y,two,x);___
131: _ VecNorm(y,NORM_2,&norm);___
132:   v = norm-18.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
133: _ PetscPrintf(PETSC_COMM_WORLD,"VecAYPX %G\n",v);___

135: _ VecSwap(x,y);___
136: _ VecNorm(y,NORM_2,&norm);___
137:   v = norm-2.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
138: _ PetscPrintf(PETSC_COMM_WORLD,"VecSwap  %G\n",v);___
139: _ VecNorm(x,NORM_2,&norm);___
140:   v = norm-18.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
141: _ PetscPrintf(PETSC_COMM_WORLD,"VecSwap  %G\n",v);___

143: _ VecWAXPY(w,two,x,y);___
144: _ VecNorm(w,NORM_2,&norm);___
145:   v = norm-38.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
146: _ PetscPrintf(PETSC_COMM_WORLD,"VecWAXPY %G\n",v);___

148: _ VecPointwiseMult(w,y,x);___
149: _ VecNorm(w,NORM_2,&norm);___
150:   v = norm-36.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
151: _ PetscPrintf(PETSC_COMM_WORLD,"VecPointwiseMult %G\n",v);___

153: _ VecPointwiseDivide(w,x,y);___
154: _ VecNorm(w,NORM_2,&norm);___
155:   v = norm-9.0*sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
156: _ PetscPrintf(PETSC_COMM_WORLD,"VecPointwiseDivide %G\n",v);___

158:   dots[0] = one;
159:   dots[1] = three;
160:   dots[2] = two;
161: _ VecSet(x,one);___
162: _ VecMAXPY(x,3,dots,z);___
163: _ VecNorm(z[0],NORM_2,&norm);___
164:   v = norm-sqrt((PetscReal) n); if (v > -PETSC_EPS && v < PETSC_EPS) v = 0.0;
165: _ VecNorm(z[1],NORM_2,&norm);___
166:   v1 = norm-2.0*sqrt((PetscReal) n); if (v1 > -PETSC_EPS && v1 < PETSC_EPS) v1 = 0.0;
167: _ VecNorm(z[2],NORM_2,&norm);___
168:   v2 = norm-3.0*sqrt((PetscReal) n); if (v2 > -PETSC_EPS && v2 < PETSC_EPS) v2 = 0.0;
169: _ PetscPrintf(PETSC_COMM_WORLD,"VecMAXPY %G %G %G \n",v,v1,v2);___


172:   /* 
173:      Free work space.  All PETSc objects should be destroyed when they
174:      are no longer needed.
175:   */
176: _ VecDestroy(&x);___
177: _ VecDestroy(&y);___
178: _ VecDestroy(&w);___
179: _ VecDestroyVecs(3,&z);___
180: _ PetscFinalize();___
181:   return 0;
182: }
183: