Actual source code: test7.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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: {
18: Vec t,v;
19: Mat B,Ymat;
20: BV X,Y,Z,Zcopy;
21: PetscInt i,j,n=10,k=5,rep=1,Istart,Iend;
22: PetscScalar *pZ;
23: PetscReal norm;
24: PetscViewer view;
25: PetscBool verbose,fromfile;
26: char filename[PETSC_MAX_PATH_LEN];
27: PetscViewer viewer;
28: BVMatMultType vmm;
30: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
31: PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL);
32: PetscOptionsGetInt(NULL,NULL,"-rep",&rep,NULL);
33: PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
34: PetscOptionsGetString(NULL,NULL,"-file",filename,PETSC_MAX_PATH_LEN,&fromfile);
35: MatCreate(PETSC_COMM_WORLD,&B);
36: PetscObjectSetName((PetscObject)B,"B");
37: if (fromfile) {
38: #if defined(PETSC_USE_COMPLEX)
39: PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n");
40: #else
41: PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n");
42: #endif
43: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
44: MatSetFromOptions(B);
45: MatLoad(B,viewer);
46: PetscViewerDestroy(&viewer);
47: MatGetSize(B,&n,NULL);
48: MatGetOwnershipRange(B,&Istart,&Iend);
49: } else {
50: /* Create 1-D Laplacian matrix */
51: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
52: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n);
53: MatSetFromOptions(B);
54: MatSetUp(B);
55: MatGetOwnershipRange(B,&Istart,&Iend);
56: for (i=Istart;i<Iend;i++) {
57: if (i>0) { MatSetValue(B,i,i-1,-1.0,INSERT_VALUES); }
58: if (i<n-1) { MatSetValue(B,i,i+1,-1.0,INSERT_VALUES); }
59: MatSetValue(B,i,i,2.0,INSERT_VALUES);
60: }
61: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
62: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
63: }
65: PetscPrintf(PETSC_COMM_WORLD,"Test BVMatMult (n=%D, k=%D).\n",n,k);
66: MatCreateVecs(B,&t,NULL);
68: /* Create BV object X */
69: BVCreate(PETSC_COMM_WORLD,&X);
70: PetscObjectSetName((PetscObject)X,"X");
71: BVSetSizesFromVec(X,t,k);
72: BVSetMatMultMethod(X,BV_MATMULT_VECS);
73: BVSetFromOptions(X);
74: BVGetMatMultMethod(X,&vmm);
75: PetscPrintf(PETSC_COMM_WORLD,"Using method: %s\n",BVMatMultTypes[vmm]);
77: /* Set up viewer */
78: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view);
79: if (verbose) {
80: PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB);
81: }
83: /* Fill X entries */
84: for (j=0;j<k;j++) {
85: BVGetColumn(X,j,&v);
86: VecSet(v,0.0);
87: for (i=Istart;i<PetscMin(j+1,Iend);i++) {
88: VecSetValue(v,i,1.0,INSERT_VALUES);
89: }
90: VecAssemblyBegin(v);
91: VecAssemblyEnd(v);
92: BVRestoreColumn(X,j,&v);
93: }
94: if (verbose) {
95: MatView(B,view);
96: BVView(X,view);
97: }
99: /* Create BV object Y */
100: BVDuplicateResize(X,k+4,&Y);
101: PetscObjectSetName((PetscObject)Y,"Y");
102: BVSetActiveColumns(Y,2,k+2);
104: /* Test BVMatMult */
105: for (i=0;i<rep;i++) {
106: BVMatMult(X,B,Y);
107: }
108: if (verbose) {
109: BVView(Y,view);
110: }
112: /* Test BVGetMat/RestoreMat */
113: BVGetMat(Y,&Ymat);
114: PetscObjectSetName((PetscObject)Ymat,"Ymat");
115: if (verbose) {
116: MatView(Ymat,view);
117: }
118: BVRestoreMat(Y,&Ymat);
120: if (!fromfile) {
121: /* Create BV object Z */
122: BVDuplicate(X,&Z);
123: PetscObjectSetName((PetscObject)Z,"Z");
125: /* Fill Z entries */
126: for (j=0;j<k;j++) {
127: BVGetColumn(Z,j,&v);
128: VecSet(v,0.0);
129: if (!Istart) { VecSetValue(v,0,1.0,ADD_VALUES); }
130: if (j<n && j>=Istart && j<Iend) { VecSetValue(v,j,1.0,ADD_VALUES); }
131: if (j+1<n && j>=Istart && j<Iend) { VecSetValue(v,j+1,-1.0,ADD_VALUES); }
132: VecAssemblyBegin(v);
133: VecAssemblyEnd(v);
134: BVRestoreColumn(Z,j,&v);
135: }
136: if (verbose) {
137: BVView(Z,view);
138: }
140: /* Save a copy of Z */
141: BVDuplicate(Z,&Zcopy);
142: BVCopy(Z,Zcopy);
144: /* Test BVMult, check result of previous operations */
145: BVMult(Z,-1.0,1.0,Y,NULL);
146: BVNorm(Z,NORM_FROBENIUS,&norm);
147: PetscPrintf(PETSC_COMM_WORLD,"Norm of error: %g\n",(double)norm);
148: }
150: /* Test BVMatMultColumn, multiply Y(:,2), result in Y(:,3) */
151: BVMatMultColumn(Y,B,2);
152: if (verbose) {
153: BVView(Y,view);
154: }
156: if (!fromfile) {
157: /* Test BVGetArray, modify Z to match Y */
158: BVCopy(Zcopy,Z);
159: BVGetArray(Z,&pZ);
160: if (Istart==0) {
161: if (Iend<3) SETERRQ(PETSC_COMM_WORLD,1,"First process must have at least 3 rows");
162: pZ[Iend] = 5.0; /* modify 3 first entries of second column */
163: pZ[Iend+1] = -4.0;
164: pZ[Iend+2] = 1.0;
165: }
166: BVRestoreArray(Z,&pZ);
167: if (verbose) {
168: BVView(Z,view);
169: }
171: /* Check result again with BVMult */
172: BVMult(Z,-1.0,1.0,Y,NULL);
173: BVNorm(Z,NORM_FROBENIUS,&norm);
174: PetscPrintf(PETSC_COMM_WORLD,"Norm of error: %g\n",(double)norm);
176: BVDestroy(&Z);
177: BVDestroy(&Zcopy);
178: }
180: BVDestroy(&X);
181: BVDestroy(&Y);
182: MatDestroy(&B);
183: VecDestroy(&t);
184: SlepcFinalize();
185: return ierr;
186: }
188: /*TEST
190: test:
191: suffix: 1
192: nsize: 1
193: args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult vecs
194: filter: grep -v "Using method"
195: output_file: output/test7_1.out
197: test:
198: suffix: 1_cuda
199: nsize: 1
200: args: -bv_type svec -mat_type aijcusparse -bv_matmult vecs
201: requires: cuda
202: filter: grep -v "Using method"
203: output_file: output/test7_1.out
205: test:
206: suffix: 2
207: nsize: 1
208: args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult mat
209: filter: grep -v "Using method"
210: output_file: output/test7_1.out
212: TEST*/