1: static char help[] = "Test to demonstrate interface for thread reductions and passing scalar values.\n\n";

  3: /*T
  4:    Concepts: PetscThreadComm^basic example: Threaded reductions and passing scalar values
  5: T*/

  7: /*
  8:   Include "petscthreadcomm.h" so that we can use the PetscThreadComm interface.
  9: */
 10: #include <petscthreadcomm.h>

 12: PetscInt    *trstarts;

 14: PetscErrorCode set_kernel(PetscInt myrank,PetscScalar *a,PetscScalar *alphap)
 15: {
 16:   PetscScalar alpha=*alphap;
 17:   PetscInt    i;

 19:   for(i=trstarts[myrank];i < trstarts[myrank+1];i++) a[i] = alpha;

 21:   return 0;
 22: }

 24: PetscErrorCode reduce_kernel(PetscInt myrank,PetscScalar *a,PetscThreadCommRedCtx red)
 25: {
 26:   PetscScalar my_sum=0.0;
 27:   PetscInt    i;

 29:   for(i=trstarts[myrank];i < trstarts[myrank+1];i++) my_sum += a[i];

 31:   PetscThreadReductionKernelBegin(myrank,red,&my_sum);
 32: 
 33:   return 0;
 34: }

 38: int main(int argc,char **argv)
 39: {
 40:   PetscErrorCode         ierr;
 41:   PetscInt               N=64;
 42:   PetscScalar           *a,sum=0.0,alpha=1.0,*scalar;
 43:   PetscThreadCommRedCtx  red;

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

 47:   PetscThreadCommView(PETSC_COMM_WORLD,PETSC_VIEWER_STDOUT_WORLD);

 49:   PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
 50:   PetscMalloc(N*sizeof(PetscScalar),&a);
 51: 
 52:   /* Set thread ownership ranges for the array */
 53:   PetscThreadCommGetOwnershipRanges(PETSC_COMM_WORLD,N,&trstarts);

 55:   /* Set a[i] = 1.0 .. i = 1,N */
 56:   /* Get location to store the scalar value alpha from threadcomm */
 57:   PetscThreadCommGetScalars(PETSC_COMM_WORLD,&scalar,PETSC_NULL,PETSC_NULL);
 58:   *scalar = alpha;
 59:   PetscThreadCommRunKernel(PETSC_COMM_WORLD,(PetscThreadKernel)set_kernel,2,a,scalar);

 61:   PetscThreadReductionBegin(PETSC_COMM_WORLD,THREADCOMM_SUM,PETSC_SCALAR,&red);
 62:   PetscThreadCommRunKernel(PETSC_COMM_WORLD,(PetscThreadKernel)reduce_kernel,2,a,red);
 63:   PetscThreadReductionEnd(red,&sum);

 65:   PetscThreadCommBarrier(PETSC_COMM_WORLD);

 67:   PetscPrintf(PETSC_COMM_SELF,"Sum(x) = %f\n",sum);
 68:   PetscFree(a);
 69:   PetscFree(trstarts);
 70:   PetscFinalize();
 71:   return 0;
 72: }