Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : 
      11             : static char help[] = "Test split reductions in BV.\n\n";
      12             : 
      13             : #include <slepcbv.h>
      14             : 
      15           8 : int main(int argc,char **argv)
      16             : {
      17           8 :   Vec            t,v,w,y,z,zsplit;
      18           8 :   BV             X;
      19           8 :   PetscInt       i,j,n=10,k=5;
      20           8 :   PetscScalar    *zarray;
      21           8 :   PetscReal      nrm;
      22             : 
      23           8 :   PetscFunctionBeginUser;
      24           8 :   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
      25           8 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      26           8 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
      27           8 :   PetscCheck(k>2,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Should specify at least k=3 columns");
      28           8 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"BV split ops (%" PetscInt_FMT " columns of dimension %" PetscInt_FMT ").\n",k,n));
      29             : 
      30             :   /* Create template vector */
      31           8 :   PetscCall(VecCreate(PETSC_COMM_WORLD,&t));
      32           8 :   PetscCall(VecSetSizes(t,PETSC_DECIDE,n));
      33           8 :   PetscCall(VecSetFromOptions(t));
      34           8 :   PetscCall(VecDuplicate(t,&v));
      35           8 :   PetscCall(VecSet(v,1.0));
      36             : 
      37             :   /* Create BV object X */
      38           8 :   PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
      39           8 :   PetscCall(PetscObjectSetName((PetscObject)X,"X"));
      40           8 :   PetscCall(BVSetSizesFromVec(X,t,k));
      41           8 :   PetscCall(BVSetFromOptions(X));
      42             : 
      43             :   /* Fill X entries */
      44          48 :   for (j=0;j<k;j++) {
      45          40 :     PetscCall(BVGetColumn(X,j,&w));
      46          40 :     PetscCall(VecSet(w,0.0));
      47         200 :     for (i=0;i<4;i++) {
      48         160 :       if (i+j<n) PetscCall(VecSetValue(w,i+j,(PetscScalar)(3*i+j-2),INSERT_VALUES));
      49             :     }
      50          40 :     PetscCall(VecAssemblyBegin(w));
      51          40 :     PetscCall(VecAssemblyEnd(w));
      52          40 :     PetscCall(BVRestoreColumn(X,j,&w));
      53             :   }
      54             : 
      55             :   /* Use regular operations */
      56           8 :   PetscCall(VecCreateSeq(PETSC_COMM_SELF,k+6,&z));
      57           8 :   PetscCall(VecGetArray(z,&zarray));
      58           8 :   PetscCall(BVGetColumn(X,0,&w));
      59           8 :   PetscCall(VecDot(w,v,zarray));
      60           8 :   PetscCall(BVRestoreColumn(X,0,&w));
      61           8 :   PetscCall(BVDotVec(X,v,zarray+1));
      62           8 :   PetscCall(BVDotColumn(X,2,zarray+1+k));
      63             : 
      64           8 :   PetscCall(BVGetColumn(X,1,&y));
      65           8 :   PetscCall(VecNorm(y,NORM_2,&nrm));
      66           8 :   PetscCall(BVRestoreColumn(X,1,&y));
      67           8 :   zarray[k+3] = nrm;
      68           8 :   PetscCall(BVNormVec(X,v,NORM_2,&nrm));
      69           8 :   zarray[k+4] = nrm;
      70           8 :   PetscCall(BVNormColumn(X,0,NORM_2,&nrm));
      71           8 :   zarray[k+5] = nrm;
      72           8 :   PetscCall(VecRestoreArray(z,&zarray));
      73             : 
      74             :   /* Use split operations */
      75           8 :   PetscCall(VecCreateSeq(PETSC_COMM_SELF,k+6,&zsplit));
      76           8 :   PetscCall(VecGetArray(zsplit,&zarray));
      77           8 :   PetscCall(BVGetColumn(X,0,&w));
      78           8 :   PetscCall(VecDotBegin(w,v,zarray));
      79           8 :   PetscCall(BVDotVecBegin(X,v,zarray+1));
      80           8 :   PetscCall(BVDotColumnBegin(X,2,zarray+1+k));
      81           8 :   PetscCall(VecDotEnd(w,v,zarray));
      82           8 :   PetscCall(BVRestoreColumn(X,0,&w));
      83           8 :   PetscCall(BVDotVecEnd(X,v,zarray+1));
      84           8 :   PetscCall(BVDotColumnEnd(X,2,zarray+1+k));
      85             : 
      86           8 :   PetscCall(BVGetColumn(X,1,&y));
      87           8 :   PetscCall(VecNormBegin(y,NORM_2,&nrm));
      88           8 :   PetscCall(BVNormVecBegin(X,v,NORM_2,&nrm));
      89           8 :   PetscCall(BVNormColumnBegin(X,0,NORM_2,&nrm));
      90           8 :   PetscCall(VecNormEnd(y,NORM_2,&nrm));
      91           8 :   PetscCall(BVRestoreColumn(X,1,&y));
      92           8 :   zarray[k+3] = nrm;
      93           8 :   PetscCall(BVNormVecEnd(X,v,NORM_2,&nrm));
      94           8 :   zarray[k+4] = nrm;
      95           8 :   PetscCall(BVNormColumnEnd(X,0,NORM_2,&nrm));
      96           8 :   zarray[k+5] = nrm;
      97           8 :   PetscCall(VecRestoreArray(zsplit,&zarray));
      98             : 
      99             :   /* Show difference */
     100           8 :   PetscCall(VecAXPY(z,-1.0,zsplit));
     101           8 :   PetscCall(VecNorm(z,NORM_1,&nrm));
     102           8 :   PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"%g\n",(double)nrm));
     103           8 :   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,NULL));
     104             : 
     105           8 :   PetscCall(BVDestroy(&X));
     106           8 :   PetscCall(VecDestroy(&t));
     107           8 :   PetscCall(VecDestroy(&v));
     108           8 :   PetscCall(VecDestroy(&z));
     109           8 :   PetscCall(VecDestroy(&zsplit));
     110           8 :   PetscCall(SlepcFinalize());
     111             :   return 0;
     112             : }
     113             : 
     114             : /*TEST
     115             : 
     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
     126             : 
     127             : TEST*/