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