Actual source code: test7.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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*/