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 multiplication of a Mat times a BV.\n\n";
12 :
13 : #include <slepcbv.h>
14 :
15 40 : int main(int argc,char **argv)
16 : {
17 40 : Vec t,r,v;
18 40 : Mat B,Ymat;
19 40 : BV X,Y,Z=NULL,Zcopy=NULL;
20 40 : PetscInt i,j,m=10,n,k=5,rep=1,Istart,Iend,ld;
21 40 : PetscScalar *pZ;
22 40 : PetscReal norm;
23 40 : PetscViewer view;
24 40 : PetscBool flg,verbose,fromfile;
25 40 : char filename[PETSC_MAX_PATH_LEN];
26 40 : PetscViewer viewer;
27 40 : BVMatMultType vmm;
28 :
29 40 : PetscFunctionBeginUser;
30 40 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
31 40 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL));
32 40 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-rep",&rep,NULL));
33 40 : PetscCall(PetscOptionsHasName(NULL,NULL,"-verbose",&verbose));
34 40 : PetscCall(PetscOptionsGetString(NULL,NULL,"-file",filename,sizeof(filename),&fromfile));
35 40 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
36 40 : PetscCall(PetscObjectSetName((PetscObject)B,"B"));
37 40 : if (fromfile) {
38 : #if defined(PETSC_USE_COMPLEX)
39 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n"));
40 : #else
41 16 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n"));
42 : #endif
43 16 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
44 16 : PetscCall(MatSetFromOptions(B));
45 16 : PetscCall(MatLoad(B,viewer));
46 16 : PetscCall(PetscViewerDestroy(&viewer));
47 16 : PetscCall(MatGetSize(B,&m,&n));
48 16 : PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
49 : } else {
50 : /* Create 1-D Laplacian matrix */
51 24 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
52 24 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,&flg));
53 24 : if (!flg) n = m;
54 24 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n));
55 24 : PetscCall(MatSetFromOptions(B));
56 24 : PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
57 376 : for (i=Istart;i<Iend;i++) {
58 352 : if (i>0 && i-1<n) PetscCall(MatSetValue(B,i,i-1,-1.0,INSERT_VALUES));
59 352 : if (i+1<n) PetscCall(MatSetValue(B,i,i+1,-1.0,INSERT_VALUES));
60 352 : if (i<n) PetscCall(MatSetValue(B,i,i,2.0,INSERT_VALUES));
61 : }
62 24 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
63 24 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
64 : }
65 :
66 40 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BVMatMult (m=%" PetscInt_FMT ", n=%" PetscInt_FMT ", k=%" PetscInt_FMT ").\n",m,n,k));
67 40 : PetscCall(MatCreateVecs(B,&t,&r));
68 :
69 : /* Create BV object X */
70 40 : PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
71 40 : PetscCall(PetscObjectSetName((PetscObject)X,"X"));
72 40 : PetscCall(BVSetSizesFromVec(X,t,k));
73 40 : PetscCall(BVSetMatMultMethod(X,BV_MATMULT_VECS));
74 40 : PetscCall(BVSetFromOptions(X));
75 40 : PetscCall(BVGetMatMultMethod(X,&vmm));
76 40 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Using method: %s\n",BVMatMultTypes[vmm]));
77 :
78 : /* Set up viewer */
79 40 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
80 40 : if (verbose) PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));
81 :
82 : /* Fill X entries */
83 304 : for (j=0;j<k;j++) {
84 264 : PetscCall(BVGetColumn(X,j,&v));
85 264 : PetscCall(VecSet(v,0.0));
86 864 : for (i=Istart;i<PetscMin(j+1,Iend);i++) PetscCall(VecSetValue(v,i,1.0,INSERT_VALUES));
87 264 : PetscCall(VecAssemblyBegin(v));
88 264 : PetscCall(VecAssemblyEnd(v));
89 264 : PetscCall(BVRestoreColumn(X,j,&v));
90 : }
91 40 : if (verbose) {
92 0 : PetscCall(MatView(B,view));
93 0 : PetscCall(BVView(X,view));
94 : }
95 :
96 : /* Create BV object Y */
97 40 : PetscCall(BVCreate(PETSC_COMM_WORLD,&Y));
98 40 : PetscCall(PetscObjectSetName((PetscObject)Y,"Y"));
99 40 : PetscCall(BVSetSizesFromVec(Y,r,k+4));
100 40 : PetscCall(BVSetMatMultMethod(Y,BV_MATMULT_VECS));
101 40 : PetscCall(BVSetFromOptions(Y));
102 40 : PetscCall(BVSetActiveColumns(Y,2,k+2));
103 :
104 : /* Test BVMatMult */
105 80 : for (i=0;i<rep;i++) PetscCall(BVMatMult(X,B,Y));
106 40 : if (verbose) PetscCall(BVView(Y,view));
107 :
108 40 : if (fromfile) {
109 : /* Test BVMatMultTranspose */
110 16 : PetscCall(BVDuplicate(X,&Z));
111 16 : PetscCall(BVSetRandom(Z));
112 32 : for (i=0;i<rep;i++) PetscCall(BVMatMultTranspose(Z,B,Y));
113 16 : if (verbose) {
114 0 : PetscCall(BVView(Z,view));
115 0 : PetscCall(BVView(Y,view));
116 : }
117 16 : PetscCall(BVDestroy(&Z));
118 16 : PetscCall(BVMatMultTransposeColumn(Y,B,2));
119 16 : if (verbose) PetscCall(BVView(Y,view));
120 : }
121 :
122 : /* Test BVGetMat/RestoreMat */
123 40 : PetscCall(BVGetMat(Y,&Ymat));
124 40 : PetscCall(PetscObjectSetName((PetscObject)Ymat,"Ymat"));
125 40 : if (verbose) PetscCall(MatView(Ymat,view));
126 40 : PetscCall(BVRestoreMat(Y,&Ymat));
127 :
128 40 : if (!fromfile) {
129 : /* Create BV object Z */
130 24 : PetscCall(BVDuplicateResize(Y,k,&Z));
131 24 : PetscCall(PetscObjectSetName((PetscObject)Z,"Z"));
132 :
133 : /* Fill Z entries */
134 208 : for (j=0;j<k;j++) {
135 184 : PetscCall(BVGetColumn(Z,j,&v));
136 184 : PetscCall(VecSet(v,0.0));
137 184 : if (!Istart) PetscCall(VecSetValue(v,0,1.0,ADD_VALUES));
138 184 : if (j<n && j>=Istart && j<Iend) PetscCall(VecSetValue(v,j,1.0,ADD_VALUES));
139 184 : if (j+1<n && j>=Istart && j<Iend) PetscCall(VecSetValue(v,j+1,-1.0,ADD_VALUES));
140 184 : PetscCall(VecAssemblyBegin(v));
141 184 : PetscCall(VecAssemblyEnd(v));
142 184 : PetscCall(BVRestoreColumn(Z,j,&v));
143 : }
144 24 : if (verbose) PetscCall(BVView(Z,view));
145 :
146 : /* Save a copy of Z */
147 24 : PetscCall(BVDuplicate(Z,&Zcopy));
148 24 : PetscCall(BVCopy(Z,Zcopy));
149 :
150 : /* Test BVMult, check result of previous operations */
151 24 : PetscCall(BVMult(Z,-1.0,1.0,Y,NULL));
152 24 : PetscCall(BVNorm(Z,NORM_FROBENIUS,&norm));
153 24 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error: %g\n",(double)norm));
154 : }
155 :
156 : /* Test BVMatMultColumn, multiply Y(:,2), result in Y(:,3) */
157 40 : if (m==n) {
158 24 : PetscCall(BVMatMultColumn(Y,B,2));
159 24 : if (verbose) PetscCall(BVView(Y,view));
160 :
161 24 : if (!fromfile) {
162 : /* Test BVGetArray, modify Z to match Y */
163 8 : PetscCall(BVCopy(Zcopy,Z));
164 8 : PetscCall(BVGetArray(Z,&pZ));
165 8 : if (Istart==0) {
166 8 : PetscCheck(Iend>2,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"First process must have at least 3 rows");
167 8 : PetscCall(BVGetLeadingDimension(Z,&ld));
168 8 : pZ[ld] = 5.0; /* modify 3 first entries of second column */
169 8 : pZ[ld+1] = -4.0;
170 8 : pZ[ld+2] = 1.0;
171 : }
172 8 : PetscCall(BVRestoreArray(Z,&pZ));
173 8 : if (verbose) PetscCall(BVView(Z,view));
174 :
175 : /* Check result again with BVMult */
176 8 : PetscCall(BVMult(Z,-1.0,1.0,Y,NULL));
177 8 : PetscCall(BVNorm(Z,NORM_FROBENIUS,&norm));
178 8 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error: %g\n",(double)norm));
179 : }
180 : }
181 :
182 40 : PetscCall(BVDestroy(&Z));
183 40 : PetscCall(BVDestroy(&Zcopy));
184 40 : PetscCall(BVDestroy(&X));
185 40 : PetscCall(BVDestroy(&Y));
186 40 : PetscCall(MatDestroy(&B));
187 40 : PetscCall(VecDestroy(&t));
188 40 : PetscCall(VecDestroy(&r));
189 40 : PetscCall(SlepcFinalize());
190 : return 0;
191 : }
192 :
193 : /*TEST
194 :
195 : testset:
196 : output_file: output/test7_1.out
197 : filter: grep -v "Using method"
198 : test:
199 : suffix: 1
200 : args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult vecs
201 : test:
202 : suffix: 1_cuda
203 : args: -bv_type {{svec mat}} -mat_type aijcusparse -bv_matmult vecs
204 : requires: cuda
205 : test:
206 : suffix: 1_hip
207 : args: -bv_type {{svec mat}} -mat_type aijhipsparse -bv_matmult vecs
208 : requires: hip
209 : test:
210 : suffix: 1_mat
211 : args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult mat
212 :
213 : testset:
214 : output_file: output/test7_2.out
215 : filter: grep -v "Using method" | sed -e "s/error: -0\./error: 0./"
216 : args: -m 34 -n 38 -k 9
217 : nsize: 2
218 : test:
219 : suffix: 2
220 : args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult vecs
221 : test:
222 : suffix: 2_cuda
223 : args: -bv_type {{svec mat}} -mat_type aijcusparse -bv_matmult vecs
224 : requires: cuda
225 : test:
226 : suffix: 2_hip
227 : args: -bv_type {{svec mat}} -mat_type aijhipsparse -bv_matmult vecs
228 : requires: hip
229 : test:
230 : suffix: 2_mat
231 : args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult mat
232 :
233 : testset:
234 : output_file: output/test7_3.out
235 : filter: grep -v "Using method"
236 : args: -file ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -bv_reproducible_random
237 : requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
238 : nsize: 2
239 : test:
240 : suffix: 3
241 : args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult {{vecs mat}}
242 :
243 : TEST*/
|
