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:   PetscFunctionBeginUser;
 17:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 18:   PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON));
 19:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 20:   n = m;
 21:   PetscCall(PetscOptionsHasName(NULL, NULL, "-rectA", &flg));
 22:   if (flg) n += 2;
 23:   PetscCall(PetscOptionsHasName(NULL, NULL, "-rectB", &flg));
 24:   if (flg) n -= 2;

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

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

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

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

 64:   /* Flush off proc Mat values and do more assembly */
 65:   PetscCall(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:       PetscCall(MatSetValues(C, 1, &i, 1, &j, &v, INSERT_VALUES));
 70:     }
 71:   }
 72:   /* Try overlap Coomunication with the next stage XXXSetValues */
 73:   PetscCall(VecAssemblyEnd(z));

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

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

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

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

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

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

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

143:   PetscCall(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

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

227:    test:
228:       suffix: 34
229:       nsize: 3
230:       args: -mat_type mpibaij -test_diagonalscale

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

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

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

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

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

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

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

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

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

291: TEST*/