Actual source code: ex1.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  2: static char help[] = "Basic vector routines.\n\n";

  4: /*T
  5:    Concepts: vectors^basic routines;
  6:    Processors: n
  7: T*/

  9: /*
 10:   Include "petscvec.h" so that we can use vectors.  Note that this file
 11:   automatically includes:
 12:      petscsys.h       - base PETSc routines   petscis.h     - index sets
 13:      petscviewer.h - viewers
 14: */

 16: #include <petscvec.h>

 20: int main(int argc,char **argv)
 21: {
 22:   Vec            x,y,w;               /* vectors */
 23:   Vec            *z;                    /* array of vectors */
 24:   PetscReal      norm,v,v1,v2,maxval;
 25:   PetscInt       n = 20,maxind;
 27:   PetscScalar    one = 1.0,two = 2.0,three = 3.0,dots[3],dot;

 29:   PetscInitialize(&argc,&argv,(char*)0,help);
 30:   PetscOptionsGetInt(NULL,"-n",&n,NULL);

 32:   /*
 33:      Create a vector, specifying only its global dimension.
 34:      When using VecCreate(), VecSetSizes() and VecSetFromOptions(), the vector format
 35:      (currently parallel, shared, or sequential) is determined at runtime.  Also, the
 36:      parallel partitioning of the vector is determined by PETSc at runtime.

 38:      Routines for creating particular vector types directly are:
 39:         VecCreateSeq() - uniprocessor vector
 40:         VecCreateMPI() - distributed vector, where the user can
 41:                          determine the parallel partitioning
 42:         VecCreateShared() - parallel vector that uses shared memory
 43:                             (available only on the SGI); otherwise,
 44:                             is the same as VecCreateMPI()

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

 49:   */

 51:   VecCreate(PETSC_COMM_WORLD,&x);
 52:   VecSetSizes(x,PETSC_DECIDE,n);
 53:   VecSetFromOptions(x);
 54:   /*
 55:      Duplicate some work vectors (of the same format and
 56:      partitioning as the initial vector).
 57:   */
 58:   VecDuplicate(x,&y);
 59:   VecDuplicate(x,&w);

 61:   /*
 62:      Duplicate more work vectors (of the same format and
 63:      partitioning as the initial vector).  Here we duplicate
 64:      an array of vectors, which is often more convenient than
 65:      duplicating individual ones.
 66:   */
 67:   VecDuplicateVecs(x,3,&z);
 68:   /*
 69:      Set the vectors to entries to a constant value.
 70:   */
 71:   VecSet(x,one);
 72:   VecSet(y,two);
 73:   VecSet(z[0],one);
 74:   VecSet(z[1],two);
 75:   VecSet(z[2],three);
 76:   /*
 77:      Demonstrate various basic vector routines.
 78:   */
 79:   MPI_Barrier(PETSC_COMM_WORLD);
 80:   VecDot(x,y,&dot);
 81:   VecMDot(x,3,z,dots);

 83:   /*
 84:      Note: If using a complex numbers version of PETSc, then
 85:      PETSC_USE_COMPLEX is defined in the makefiles; otherwise,
 86:      (when using real numbers) it is undefined.
 87:   */

 89:   PetscPrintf(PETSC_COMM_WORLD,"Vector length %D\n",n);
 90:   VecMax(x,&maxind,&maxval);
 91:   PetscPrintf(PETSC_COMM_WORLD,"VecMax %g, VecInd %D\n",(double)maxval,maxind);

 93:   VecMin(x,&maxind,&maxval);
 94:   PetscPrintf(PETSC_COMM_WORLD,"VecMin %g, VecInd %D\n",(double)maxval,maxind);
 95:   PetscPrintf(PETSC_COMM_WORLD,"All other values should be near zero\n");


 98:   VecScale(x,two);
 99:   VecNorm(x,NORM_2,&norm);
100:   v    = norm-2.0*PetscSqrtReal((PetscReal)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
101:   PetscPrintf(PETSC_COMM_WORLD,"VecScale %g\n",(double)v);


104:   VecCopy(x,w);
105:   VecNorm(w,NORM_2,&norm);
106:   v    = norm-2.0*PetscSqrtReal((PetscReal)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
107:   PetscPrintf(PETSC_COMM_WORLD,"VecCopy  %g\n",(double)v);

109:   VecAXPY(y,three,x);
110:   VecNorm(y,NORM_2,&norm);
111:   v    = norm-8.0*PetscSqrtReal((PetscReal)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
112:   PetscPrintf(PETSC_COMM_WORLD,"VecAXPY %g\n",(double)v);

114:   VecAYPX(y,two,x);
115:   VecNorm(y,NORM_2,&norm);
116:   v    = norm-18.0*PetscSqrtReal((PetscReal)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
117:   PetscPrintf(PETSC_COMM_WORLD,"VecAYPX %g\n",(double)v);

119:   VecSwap(x,y);
120:   VecNorm(y,NORM_2,&norm);
121:   v    = norm-2.0*PetscSqrtReal((PetscReal)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
122:   PetscPrintf(PETSC_COMM_WORLD,"VecSwap  %g\n",(double)v);
123:   VecNorm(x,NORM_2,&norm);
124:   v = norm-18.0*PetscSqrtReal((PetscReal)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
125:   PetscPrintf(PETSC_COMM_WORLD,"VecSwap  %g\n",(double)v);

127:   VecWAXPY(w,two,x,y);
128:   VecNorm(w,NORM_2,&norm);
129:   v    = norm-38.0*PetscSqrtReal((PetscReal)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
130:   PetscPrintf(PETSC_COMM_WORLD,"VecWAXPY %g\n",(double)v);

132:   VecPointwiseMult(w,y,x);
133:   VecNorm(w,NORM_2,&norm);
134:   v    = norm-36.0*PetscSqrtReal((PetscReal)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
135:   PetscPrintf(PETSC_COMM_WORLD,"VecPointwiseMult %g\n",(double)v);

137:   VecPointwiseDivide(w,x,y);
138:   VecNorm(w,NORM_2,&norm);
139:   v    = norm-9.0*PetscSqrtReal((PetscReal)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
140:   PetscPrintf(PETSC_COMM_WORLD,"VecPointwiseDivide %g\n",(double)v);

142:   dots[0] = one;
143:   dots[1] = three;
144:   dots[2] = two;

146:   VecSet(x,one);
147:   VecMAXPY(x,3,dots,z);
148:   VecNorm(z[0],NORM_2,&norm);
149:   v    = norm-PetscSqrtReal((PetscReal)n); if (v > -PETSC_SMALL && v < PETSC_SMALL) v = 0.0;
150:   VecNorm(z[1],NORM_2,&norm);
151:   v1   = norm-2.0*PetscSqrtReal((PetscReal)n); if (v1 > -PETSC_SMALL && v1 < PETSC_SMALL) v1 = 0.0;
152:   VecNorm(z[2],NORM_2,&norm);
153:   v2   = norm-3.0*PetscSqrtReal((PetscReal)n); if (v2 > -PETSC_SMALL && v2 < PETSC_SMALL) v2 = 0.0;
154:   PetscPrintf(PETSC_COMM_WORLD,"VecMAXPY %g %g %g \n",(double)v,(double)v1,(double)v2);

156:   /*
157:      Free work space.  All PETSc objects should be destroyed when they
158:      are no longer needed.
159:   */
160:   VecDestroy(&x);
161:   VecDestroy(&y);
162:   VecDestroy(&w);
163:   VecDestroyVecs(3,&z);
164:   PetscFinalize();
165:   return 0;
166: }