Actual source code: ex1e.c

      2: /* Program usage:  mpirun 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:      petsc.h       - base PETSc routines   petscis.h     - index sets
     21:      petscsys.h    - system routines       petscviewer.h - viewers
     22: */

     24:  #include petscvec.h

     26: #if defined(PETSC_USE_SINGLE)
     27: #define PETSC_EPS 1.e-5
     28: #else
     29: #define PETSC_EPS 1.e-10
     30: #endif

     34: int main(int argc,char **argv)
     35: {
     36:   Vec         x, y, w;               /* vectors */
     37:   Vec         *z;                    /* array of vectors */
     38:   PetscReal   norm, v, v1, v2;
     39:   PetscInt    n = 20;
     40:   PetscTruth  flg;
     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(3,x,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);___

    171:   /*
    172:      Test whether vector has been corrupted (just to demonstrate this
    173:      routine) not needed in most application codes.
    174:   */
    175: _ VecValid(x,&flg);___
    176:   if (!flg) SETERRQ(1,"Corrupted vector.");

    178:   /*
    179:      Free work space.  All PETSc objects should be destroyed when they
    180:      are no longer needed.
    181:   */
    182: _ VecDestroy(x);___
    183: _ VecDestroy(y);___
    184: _ VecDestroy(w);___
    185: _ VecDestroyVecs(z,3);___
    186: _ PetscFinalize();___
    187:   return 0;
    188: }
    189: