Actual source code: test7.c

slepc-3.19.2 2023-09-05
Report Typos and Errors
  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:   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(MatSetUp(B));
 57:     PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
 58:     for (i=Istart;i<Iend;i++) {
 59:       if (i>0 && i-1<n) PetscCall(MatSetValue(B,i,i-1,-1.0,INSERT_VALUES));
 60:       if (i+1<n) PetscCall(MatSetValue(B,i,i+1,-1.0,INSERT_VALUES));
 61:       if (i<n) PetscCall(MatSetValue(B,i,i,2.0,INSERT_VALUES));
 62:     }
 63:     PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
 64:     PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
 65:   }

 67:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Test BVMatMult (m=%" PetscInt_FMT ", n=%" PetscInt_FMT ", k=%" PetscInt_FMT ").\n",m,n,k));
 68:   PetscCall(MatCreateVecs(B,&t,&r));

 70:   /* Create BV object X */
 71:   PetscCall(BVCreate(PETSC_COMM_WORLD,&X));
 72:   PetscCall(PetscObjectSetName((PetscObject)X,"X"));
 73:   PetscCall(BVSetSizesFromVec(X,t,k));
 74:   PetscCall(BVSetMatMultMethod(X,BV_MATMULT_VECS));
 75:   PetscCall(BVSetFromOptions(X));
 76:   PetscCall(BVGetMatMultMethod(X,&vmm));
 77:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Using method: %s\n",BVMatMultTypes[vmm]));

 79:   /* Set up viewer */
 80:   PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view));
 81:   if (verbose) PetscCall(PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB));

 83:   /* Fill X entries */
 84:   for (j=0;j<k;j++) {
 85:     PetscCall(BVGetColumn(X,j,&v));
 86:     PetscCall(VecSet(v,0.0));
 87:     for (i=Istart;i<PetscMin(j+1,Iend);i++) PetscCall(VecSetValue(v,i,1.0,INSERT_VALUES));
 88:     PetscCall(VecAssemblyBegin(v));
 89:     PetscCall(VecAssemblyEnd(v));
 90:     PetscCall(BVRestoreColumn(X,j,&v));
 91:   }
 92:   if (verbose) {
 93:     PetscCall(MatView(B,view));
 94:     PetscCall(BVView(X,view));
 95:   }

 97:   /* Create BV object Y */
 98:   PetscCall(BVCreate(PETSC_COMM_WORLD,&Y));
 99:   PetscCall(PetscObjectSetName((PetscObject)Y,"Y"));
100:   PetscCall(BVSetSizesFromVec(Y,r,k+4));
101:   PetscCall(BVSetMatMultMethod(Y,BV_MATMULT_VECS));
102:   PetscCall(BVSetFromOptions(Y));
103:   PetscCall(BVSetActiveColumns(Y,2,k+2));

105:   /* Test BVMatMult */
106:   for (i=0;i<rep;i++) PetscCall(BVMatMult(X,B,Y));
107:   if (verbose) PetscCall(BVView(Y,view));

109:   if (fromfile) {
110:     /* Test BVMatMultTranspose */
111:     PetscCall(BVDuplicate(X,&Z));
112:     PetscCall(BVSetRandom(Z));
113:     for (i=0;i<rep;i++) PetscCall(BVMatMultTranspose(Z,B,Y));
114:     if (verbose) {
115:       PetscCall(BVView(Z,view));
116:       PetscCall(BVView(Y,view));
117:     }
118:     PetscCall(BVDestroy(&Z));
119:     PetscCall(BVMatMultTransposeColumn(Y,B,2));
120:     if (verbose) PetscCall(BVView(Y,view));
121:   }

123:   /* Test BVGetMat/RestoreMat */
124:   PetscCall(BVGetMat(Y,&Ymat));
125:   PetscCall(PetscObjectSetName((PetscObject)Ymat,"Ymat"));
126:   if (verbose) PetscCall(MatView(Ymat,view));
127:   PetscCall(BVRestoreMat(Y,&Ymat));

129:   if (!fromfile) {
130:     /* Create BV object Z */
131:     PetscCall(BVDuplicateResize(Y,k,&Z));
132:     PetscCall(PetscObjectSetName((PetscObject)Z,"Z"));

134:     /* Fill Z entries */
135:     for (j=0;j<k;j++) {
136:       PetscCall(BVGetColumn(Z,j,&v));
137:       PetscCall(VecSet(v,0.0));
138:       if (!Istart) PetscCall(VecSetValue(v,0,1.0,ADD_VALUES));
139:       if (j<n && j>=Istart && j<Iend) PetscCall(VecSetValue(v,j,1.0,ADD_VALUES));
140:       if (j+1<n && j>=Istart && j<Iend) PetscCall(VecSetValue(v,j+1,-1.0,ADD_VALUES));
141:       PetscCall(VecAssemblyBegin(v));
142:       PetscCall(VecAssemblyEnd(v));
143:       PetscCall(BVRestoreColumn(Z,j,&v));
144:     }
145:     if (verbose) PetscCall(BVView(Z,view));

147:     /* Save a copy of Z */
148:     PetscCall(BVDuplicate(Z,&Zcopy));
149:     PetscCall(BVCopy(Z,Zcopy));

151:     /* Test BVMult, check result of previous operations */
152:     PetscCall(BVMult(Z,-1.0,1.0,Y,NULL));
153:     PetscCall(BVNorm(Z,NORM_FROBENIUS,&norm));
154:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error: %g\n",(double)norm));
155:   }

157:   /* Test BVMatMultColumn, multiply Y(:,2), result in Y(:,3) */
158:   if (m==n) {
159:     PetscCall(BVMatMultColumn(Y,B,2));
160:     if (verbose) PetscCall(BVView(Y,view));

162:     if (!fromfile) {
163:       /* Test BVGetArray, modify Z to match Y */
164:       PetscCall(BVCopy(Zcopy,Z));
165:       PetscCall(BVGetArray(Z,&pZ));
166:       if (Istart==0) {
167:         PetscCheck(Iend>2,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"First process must have at least 3 rows");
168:         pZ[Iend]   = 5.0;   /* modify 3 first entries of second column */
169:         pZ[Iend+1] = -4.0;
170:         pZ[Iend+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:       args: -m 12
199:       test:
200:          suffix: 1
201:          args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult vecs
202:       test:
203:          suffix: 1_cuda
204:          args: -bv_type svec -mat_type aijcusparse -bv_matmult vecs
205:          requires: cuda
206:       test:
207:          suffix: 1_mat
208:          args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult mat

210:    testset:
211:       output_file: output/test7_2.out
212:       filter: grep -v "Using method"
213:       args: -m 40 -n 44 -k 9
214:       nsize: 2
215:       test:
216:          suffix: 2
217:          args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult vecs
218:       test:
219:          suffix: 2_cuda
220:          args: -bv_type svec -mat_type aijcusparse -bv_matmult vecs
221:          requires: cuda
222:       test:
223:          suffix: 2_mat
224:          args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult mat

226:    testset:
227:       output_file: output/test7_3.out
228:       filter: grep -v "Using method"
229:       args: -file ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -bv_reproducible_random
230:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
231:       nsize: 2
232:       test:
233:          suffix: 3
234:          args: -bv_type {{vecs contiguous svec mat}shared output} -bv_matmult {{vecs mat}}

236: TEST*/