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: