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;
 11:   PetscInt       i,j,m = 8,n,rstart,rend,vstart,vend;
 12:   PetscScalar    one = 1.0,negone = -1.0,v,alpha=0.1;
 13:   PetscReal      norm, tol = PETSC_SQRT_MACHINE_EPSILON;
 14:   PetscBool      flg;

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

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

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

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

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

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

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

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

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

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

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

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

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

143:   PetscFinalize();
144:   return 0;
145: }

147: /*TEST

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

293: TEST*/