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