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: