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,NULL,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 : test:
127 : suffix: 1_hip
128 : args: -bv_type {{svec mat}} -vec_type hip
129 : requires: hip
130 :
131 : TEST*/
|
