1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */

 11: static char help[] = "Test split reductions in BV.\n\n";

 13: #include <slepcbv.h>

 15: int main(int argc,char **argv)
 16: {
 17:   Vec            t,v,w,y,z,zsplit;
 18:   BV             X;
 19:   PetscInt       i,j,n=10,k=5;
 20:   PetscScalar    *zarray;
 21:   PetscReal      nrm;

 23:   PetscFunctionBeginUser;
 24:   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
 25:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
 26:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
 27:   PetscCheck(k>2,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Should specify at least k=3 columns");
 28:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"BV split ops (%" PetscInt_FMT " columns of dimension %" PetscInt_FMT ").\n",k,n));

 30:   /* Create template vector */
 31:   PetscCall(VecCreate(PETSC_COMM_WORLD,&t));
 32:   PetscCall(VecSetSizes(t,PETSC_DECIDE,n));
 33:   PetscCall(VecSetFromOptions(t));
 34:   PetscCall(VecDuplicate(t,&v));
 35:   PetscCall(VecSet(v,1.0));

 37:   /* Create BV object X */
 38:   PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
 39:   PetscCall(PetscObjectSetName((PetscObject)X,"X"));
 40:   PetscCall(BVSetSizesFromVec(X,t,k));
 41:   PetscCall(BVSetFromOptions(X));

 43:   /* Fill X entries */
 44:   for (j=0;j<k;j++) {
 45:     PetscCall(BVGetColumn(X,j,&w));
 46:     PetscCall(VecSet(w,0.0));
 47:     for (i=0;i<4;i++) {
 48:       if (i+j<n) PetscCall(VecSetValue(w,i+j,(PetscScalar)(3*i+j-2),INSERT_VALUES));
 49:     }
 50:     PetscCall(VecAssemblyBegin(w));
 51:     PetscCall(VecAssemblyEnd(w));
 52:     PetscCall(BVRestoreColumn(X,j,&w));
 53:   }

 55:   /* Use regular operations */
 56:   PetscCall(VecCreateSeq(PETSC_COMM_SELF,k+6,&z));
 57:   PetscCall(VecGetArray(z,&zarray));
 58:   PetscCall(BVGetColumn(X,0,&w));
 59:   PetscCall(VecDot(w,v,zarray));
 60:   PetscCall(BVRestoreColumn(X,0,&w));
 61:   PetscCall(BVDotVec(X,v,zarray+1));
 62:   PetscCall(BVDotColumn(X,2,zarray+1+k));

 64:   PetscCall(BVGetColumn(X,1,&y));
 65:   PetscCall(VecNorm(y,NORM_2,&nrm));
 66:   PetscCall(BVRestoreColumn(X,1,&y));
 67:   zarray[k+3] = nrm;
 68:   PetscCall(BVNormVec(X,v,NORM_2,&nrm));
 69:   zarray[k+4] = nrm;
 70:   PetscCall(BVNormColumn(X,0,NORM_2,&nrm));
 71:   zarray[k+5] = nrm;
 72:   PetscCall(VecRestoreArray(z,&zarray));

 74:   /* Use split operations */
 75:   PetscCall(VecCreateSeq(PETSC_COMM_SELF,k+6,&zsplit));
 76:   PetscCall(VecGetArray(zsplit,&zarray));
 77:   PetscCall(BVGetColumn(X,0,&w));
 78:   PetscCall(VecDotBegin(w,v,zarray));
 79:   PetscCall(BVDotVecBegin(X,v,zarray+1));
 80:   PetscCall(BVDotColumnBegin(X,2,zarray+1+k));
 81:   PetscCall(VecDotEnd(w,v,zarray));
 82:   PetscCall(BVRestoreColumn(X,0,&w));
 83:   PetscCall(BVDotVecEnd(X,v,zarray+1));
 84:   PetscCall(BVDotColumnEnd(X,2,zarray+1+k));

 86:   PetscCall(BVGetColumn(X,1,&y));
 87:   PetscCall(VecNormBegin(y,NORM_2,&nrm));
 88:   PetscCall(BVNormVecBegin(X,v,NORM_2,&nrm));
 89:   PetscCall(BVNormColumnBegin(X,0,NORM_2,&nrm));
 90:   PetscCall(VecNormEnd(y,NORM_2,&nrm));
 91:   PetscCall(BVRestoreColumn(X,1,&y));
 92:   zarray[k+3] = nrm;
 93:   PetscCall(BVNormVecEnd(X,v,NORM_2,&nrm));
 94:   zarray[k+4] = nrm;
 95:   PetscCall(BVNormColumnEnd(X,0,NORM_2,&nrm));
 96:   zarray[k+5] = nrm;
 97:   PetscCall(VecRestoreArray(zsplit,&zarray));

 99:   /* Show difference */
100:   PetscCall(VecAXPY(z,-1.0,zsplit));
101:   PetscCall(VecNorm(z,NORM_1,&nrm));
102:   PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"%g\n",(double)nrm));
103:   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,NULL));

105:   PetscCall(BVDestroy(&X));
106:   PetscCall(VecDestroy(&t));
107:   PetscCall(VecDestroy(&v));
108:   PetscCall(VecDestroy(&z));
109:   PetscCall(VecDestroy(&zsplit));
110:   PetscCall(SlepcFinalize());
111:   return 0;
112: }

114: /*TEST

116:    testset:
117:       nsize: 2
118:       output_file: output/test10_1.out
119:       test:
120:          suffix: 1
121:          args: -bv_type {{vecs contiguous svec mat}shared output}
122:       test:
123:          suffix: 1_cuda
124:          args: -bv_type {{svec mat}} -vec_type cuda
125:          requires: cuda

127: TEST*/