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