Actual source code: ex4.c

petsc-3.3-p7 2013-05-11
  1: static char help[] = "Test Basic vector routines.\n\n";

  3: /*
  4:   Include "petscthreadcomm.h" so that we can use the PetscThreadComm interface.
  5: */
  6: #include <petscthreadcomm.h>
  7: #include <petscvec.h>

 11: int main(int argc,char **argv)
 12: {
 14:   PetscScalar    dot=0.0,v;
 15:   Vec            x,y;
 16:   PetscInt       N=8;
 17:   PetscScalar    one=1.0,two=2.0,alpha=2.0;

 19:   PetscInitialize(&argc,&argv,(char *)0,help);

 21: #if defined(PETSC_THREADCOMM_ACTIVE)
 22:   PetscThreadCommView(PETSC_COMM_WORLD,PETSC_VIEWER_STDOUT_WORLD);
 23: #endif
 24:   PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);

 26:   VecCreate(PETSC_COMM_WORLD,&x);
 27:   VecSetSizes(x,PETSC_DECIDE,N);
 28:   VecSetFromOptions(x);
 29:   VecSet(x,one);
 30:   PetscPrintf(PETSC_COMM_WORLD,"x = %lf\n",one);

 32:   VecCreate(PETSC_COMM_WORLD,&y);
 33:   VecSetSizes(y,PETSC_DECIDE,N);
 34:   VecSetFromOptions(y);
 35:   VecSet(y,two);
 36:   PetscPrintf(PETSC_COMM_WORLD,"y = %lf\n",two);

 38:   VecAXPY(y,alpha,x);
 39:   v = two+alpha*one;
 40:   PetscPrintf(PETSC_COMM_WORLD,"x+%lfy = %lf\n",alpha,v);
 41: 
 42:   VecDot(x,y,&dot);

 44: #if defined(PETSC_THREADCOMM_ACTIVE)
 45:   PetscThreadCommBarrier(PETSC_COMM_WORLD);
 46: #endif

 48:   PetscPrintf(PETSC_COMM_WORLD,"Dot product %d*(%lf*%lf) is %lf\n",N,one,v,dot);
 49:   VecDestroy(&x);
 50:   VecDestroy(&y);
 51:   PetscFinalize();
 52:   return 0;
 53: }