Actual source code: ex81.c
2: static char help[] = "Tests MatOption MAT_FORCE_DIAGONAL_ENTRIES.\n\n";
4: #include <petscmat.h>
6: int main(int argc, char **args)
7: {
8: Mat A, B;
9: Vec diag;
10: PetscInt i, n = 10, col[3], test;
11: PetscScalar v[3];
14: PetscInitialize(&argc, &args, (char *)0, help);
15: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
17: /* Create A which has empty 0-th row and column */
18: MatCreate(PETSC_COMM_WORLD, &A);
19: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
20: MatSetType(A, MATAIJ);
21: MatSetFromOptions(A);
22: MatSetUp(A);
24: v[0] = -1.;
25: v[1] = 2.;
26: v[2] = -1.;
27: for (i = 2; i < n - 1; i++) {
28: col[0] = i - 1;
29: col[1] = i;
30: col[2] = i + 1;
31: MatSetValues(A, 1, &i, 3, col, v, INSERT_VALUES);
32: }
33: i = 1;
34: col[0] = 1;
35: col[1] = 2;
36: MatSetValues(A, 1, &i, 2, col, v + 1, INSERT_VALUES);
37: i = n - 1;
38: col[0] = n - 2;
39: col[1] = n - 1;
40: MatSetValues(A, 1, &i, 2, col, v, INSERT_VALUES);
41: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
42: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
44: for (test = 0; test < 2; test++) {
45: MatProductCreate(A, A, NULL, &B);
47: if (test == 0) {
48: /* Compute B = A*A; B misses 0-th diagonal */
49: MatProductSetType(B, MATPRODUCT_AB);
50: MatSetOptionsPrefix(B, "AB_");
51: } else {
52: /* Compute B = A^t*A; B misses 0-th diagonal */
53: MatProductSetType(B, MATPRODUCT_AtB);
54: MatSetOptionsPrefix(B, "AtB_");
55: }
57: /* Force allocate missing diagonal entries of B */
58: MatSetOption(B, MAT_FORCE_DIAGONAL_ENTRIES, PETSC_TRUE);
59: MatProductSetFromOptions(B);
61: MatProductSymbolic(B);
62: MatProductNumeric(B);
64: MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
66: /* Insert entries to diagonal of B */
67: MatCreateVecs(B, NULL, &diag);
68: MatGetDiagonal(B, diag);
69: VecSetValue(diag, 0, 100.0, INSERT_VALUES);
70: VecAssemblyBegin(diag);
71: VecAssemblyEnd(diag);
73: MatDiagonalSet(B, diag, INSERT_VALUES);
74: if (test == 1) MatView(B, PETSC_VIEWER_STDOUT_WORLD);
75: MatDestroy(&B);
76: VecDestroy(&diag);
77: }
79: MatDestroy(&A);
80: PetscFinalize();
81: return 0;
82: }
84: /*TEST
86: test:
87: output_file: output/ex81_1.out
89: test:
90: suffix: 2
91: args: -AtB_mat_product_algorithm at*b
92: output_file: output/ex81_1.out
94: test:
95: suffix: 3
96: args: -AtB_mat_product_algorithm outerproduct
97: output_file: output/ex81_1.out
99: test:
100: suffix: 4
101: nsize: 3
102: args: -AtB_mat_product_algorithm nonscalable
103: output_file: output/ex81_3.out
105: test:
106: suffix: 5
107: nsize: 3
108: args: -AtB_mat_product_algorithm scalable
109: output_file: output/ex81_3.out
111: test:
112: suffix: 6
113: nsize: 3
114: args: -AtB_mat_product_algorithm at*b
115: output_file: output/ex81_3.out
117: TEST*/