Actual source code: ex5.c


  2: static char help[] = "Tests MatMult(), MatMultAdd(), MatMultTranspose().\n\
  3: Also MatMultTransposeAdd(), MatScale(), MatGetDiagonal(), and MatDiagonalScale().\n\n";

  5: #include <petscmat.h>

  7: int main(int argc,char **args)
  8: {
  9:   Mat            C;
 10:   Vec            s,u,w,x,y,z;
 12:   PetscInt       i,j,m = 8,n,rstart,rend,vstart,vend;
 13:   PetscScalar    one = 1.0,negone = -1.0,v,alpha=0.1;
 14:   PetscReal      norm, tol = PETSC_SQRT_MACHINE_EPSILON;
 15:   PetscBool      flg;

 17:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 18:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_COMMON);
 19:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 20:   n    = m;
 21:   PetscOptionsHasName(NULL,NULL,"-rectA",&flg);
 22:   if (flg) n += 2;
 23:   PetscOptionsHasName(NULL,NULL,"-rectB",&flg);
 24:   if (flg) n -= 2;

 26:   /* ---------- Assemble matrix and vectors ----------- */

 28:   MatCreate(PETSC_COMM_WORLD,&C);
 29:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m,n);
 30:   MatSetFromOptions(C);
 31:   MatSetUp(C);
 32:   MatGetOwnershipRange(C,&rstart,&rend);
 33:   VecCreate(PETSC_COMM_WORLD,&x);
 34:   VecSetSizes(x,PETSC_DECIDE,m);
 35:   VecSetFromOptions(x);
 36:   VecDuplicate(x,&z);
 37:   VecDuplicate(x,&w);
 38:   VecCreate(PETSC_COMM_WORLD,&y);
 39:   VecSetSizes(y,PETSC_DECIDE,n);
 40:   VecSetFromOptions(y);
 41:   VecDuplicate(y,&u);
 42:   VecDuplicate(y,&s);
 43:   VecGetOwnershipRange(y,&vstart,&vend);

 45:   /* Assembly */
 46:   for (i=rstart; i<rend; i++) {
 47:     v    = 100*(i+1);
 48:     VecSetValues(z,1,&i,&v,INSERT_VALUES);
 49:     for (j=0; j<n; j++) {
 50:       v    = 10*(i+1)+j+1;
 51:       MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES);
 52:     }
 53:   }

 55:   /* Flush off proc Vec values and do more assembly */
 56:   VecAssemblyBegin(z);
 57:   for (i=vstart; i<vend; i++) {
 58:     v    = one*((PetscReal)i);
 59:     VecSetValues(y,1,&i,&v,INSERT_VALUES);
 60:     v    = 100.0*i;
 61:     VecSetValues(u,1,&i,&v,INSERT_VALUES);
 62:   }

 64:   /* Flush off proc Mat values and do more assembly */
 65:   MatAssemblyBegin(C,MAT_FLUSH_ASSEMBLY);
 66:   for (i=rstart; i<rend; i++) {
 67:     for (j=0; j<n; j++) {
 68:       v    = 10*(i+1)+j+1;
 69:       MatSetValues(C,1,&i,1,&j,&v,INSERT_VALUES);
 70:     }
 71:   }
 72:   /* Try overlap Coomunication with the next stage XXXSetValues */
 73:   VecAssemblyEnd(z);

 75:   MatAssemblyEnd(C,MAT_FLUSH_ASSEMBLY);
 76:   CHKMEMQ;
 77:   /* The Assembly for the second Stage */
 78:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 79:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
 80:   VecAssemblyBegin(y);
 81:   VecAssemblyEnd(y);
 82:   MatScale(C,alpha);
 83:   VecAssemblyBegin(u);
 84:   VecAssemblyEnd(u);
 85:   PetscPrintf(PETSC_COMM_WORLD,"testing MatMult()\n");
 86:   CHKMEMQ;
 87:   MatMult(C,y,x);
 88:   CHKMEMQ;
 89:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
 90:   PetscPrintf(PETSC_COMM_WORLD,"testing MatMultAdd()\n");
 91:   MatMultAdd(C,y,z,w);
 92:   VecAXPY(x,one,z);
 93:   VecAXPY(x,negone,w);
 94:   VecNorm(x,NORM_2,&norm);
 95:   if (norm > tol) {
 96:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm);
 97:   }

 99:   /* ------- Test MatMultTranspose(), MatMultTransposeAdd() ------- */

101:   for (i=rstart; i<rend; i++) {
102:     v    = one*((PetscReal)i);
103:     VecSetValues(x,1,&i,&v,INSERT_VALUES);
104:   }
105:   VecAssemblyBegin(x);
106:   VecAssemblyEnd(x);
107:   PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTranspose()\n");
108:   MatMultTranspose(C,x,y);
109:   VecView(y,PETSC_VIEWER_STDOUT_WORLD);

111:   PetscPrintf(PETSC_COMM_WORLD,"testing MatMultTransposeAdd()\n");
112:   MatMultTransposeAdd(C,x,u,s);
113:   VecAXPY(y,one,u);
114:   VecAXPY(y,negone,s);
115:   VecNorm(y,NORM_2,&norm);
116:   if (norm > tol) {
117:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error difference = %g\n",(double)norm);
118:   }

120:   /* -------------------- Test MatGetDiagonal() ------------------ */

122:   PetscPrintf(PETSC_COMM_WORLD,"testing MatGetDiagonal(), MatDiagonalScale()\n");
123:   MatView(C,PETSC_VIEWER_STDOUT_WORLD);
124:   VecSet(x,one);
125:   MatGetDiagonal(C,x);
126:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
127:   for (i=vstart; i<vend; i++) {
128:     v    = one*((PetscReal)(i+1));
129:     VecSetValues(y,1,&i,&v,INSERT_VALUES);
130:   }

132:   /* -------------------- Test () MatDiagonalScale ------------------ */
133:   PetscOptionsHasName(NULL,NULL,"-test_diagonalscale",&flg);
134:   if (flg) {
135:     MatDiagonalScale(C,x,y);
136:     MatView(C,PETSC_VIEWER_STDOUT_WORLD);
137:   }
138:   /* Free data structures */
139:   VecDestroy(&u); VecDestroy(&s);
140:   VecDestroy(&w); VecDestroy(&x);
141:   VecDestroy(&y); VecDestroy(&z);
142:   MatDestroy(&C);

144:   PetscFinalize();
145:   return ierr;
146: }

148: /*TEST

150:    test:
151:       suffix: 11_A
152:       args: -mat_type seqaij -rectA
153:       filter: grep -v type

155:    test:
156:       args: -mat_type seqdense -rectA
157:       suffix: 12_A

159:    test:
160:       args: -mat_type seqaij -rectB
161:       suffix: 11_B
162:       filter: grep -v type

164:    test:
165:       args: -mat_type seqdense -rectB
166:       suffix: 12_B

168:    test:
169:       suffix: 21
170:       args: -mat_type mpiaij
171:       filter: grep -v type

173:    test:
174:       suffix: 22
175:       args: -mat_type mpidense

177:    test:
178:       suffix: 23
179:       nsize: 3
180:       args: -mat_type mpiaij
181:       filter: grep -v type

183:    test:
184:       suffix: 24
185:       nsize: 3
186:       args: -mat_type mpidense

188:    test:
189:       suffix: 2_aijcusparse_1
190:       args: -mat_type mpiaijcusparse -vec_type cuda
191:       filter: grep -v type
192:       output_file: output/ex5_21.out
193:       requires: cuda

195:    test:
196:       nsize: 3
197:       suffix: 2_aijcusparse_2
198:       filter: grep -v type
199:       args: -mat_type mpiaijcusparse -vec_type cuda
200:       args: -sf_type {{basic neighbor}} -vecscatter_packongpu {{0 1}}
201:       output_file: output/ex5_23.out
202:       requires: cuda

204:    test:
205:       nsize: 3
206:       suffix: 2_aijcusparse_3
207:       filter: grep -v type
208:       args: -mat_type mpiaijcusparse -vec_type cuda
209:       args: -sf_type {{basic neighbor}}
210:       output_file: output/ex5_23.out
211:       requires: cuda defined(PETSC_HAVE_MPI_GPU_AWARE)

213:    test:
214:       suffix: 31
215:       args: -mat_type mpiaij -test_diagonalscale
216:       filter: grep -v type

218:    test:
219:       suffix: 32
220:       args: -mat_type mpibaij -test_diagonalscale
221:       filter: grep -v Mat_

223:    test:
224:       suffix: 33
225:       nsize: 3
226:       args: -mat_type mpiaij -test_diagonalscale
227:       filter: grep -v type

229:    test:
230:       suffix: 34
231:       nsize: 3
232:       args: -mat_type mpibaij -test_diagonalscale
233:       filter: grep -v Mat_

235:    test:
236:       suffix: 3_aijcusparse_1
237:       args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
238:       filter: grep -v type
239:       output_file: output/ex5_31.out
240:       requires: cuda

242:    test:
243:       suffix: 3_aijcusparse_2
244:       nsize: 3
245:       args: -mat_type mpiaijcusparse -vec_type cuda -test_diagonalscale
246:       filter: grep -v type
247:       output_file: output/ex5_33.out
248:       requires: cuda

250:    test:
251:       suffix: 3_kokkos
252:       nsize: 3
253:       args: -mat_type mpiaijkokkos -vec_type kokkos -test_diagonalscale
254:       filter: grep -v type
255:       output_file: output/ex5_33.out
256:       requires: kokkos_kernels

258:    test:
259:       suffix: aijcusparse_1
260:       args: -mat_type seqaijcusparse -vec_type cuda -rectA
261:       filter: grep -v type
262:       output_file: output/ex5_11_A.out
263:       requires: cuda

265:    test:
266:       suffix: aijcusparse_2
267:       args: -mat_type seqaijcusparse -vec_type cuda -rectB
268:       filter: grep -v type
269:       output_file: output/ex5_11_B.out
270:       requires: cuda

272:    test:
273:       suffix: sell_1
274:       args: -mat_type sell
275:       output_file: output/ex5_41.out

277:    test:
278:       suffix: sell_2
279:       nsize: 3
280:       args: -mat_type sell
281:       output_file: output/ex5_43.out

283:    test:
284:       suffix: sell_3
285:       args: -mat_type sell -test_diagonalscale
286:       output_file: output/ex5_51.out

288:    test:
289:       suffix: sell_4
290:       nsize: 3
291:       args: -mat_type sell -test_diagonalscale
292:       output_file: output/ex5_53.out

294: TEST*/