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];
13: PetscFunctionBeginUser;
14: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
15: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
17: /* Create A which has empty 0-th row and column */
18: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
19: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
20: PetscCall(MatSetType(A, MATAIJ));
21: PetscCall(MatSetFromOptions(A));
22: PetscCall(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: PetscCall(MatSetValues(A, 1, &i, 3, col, v, INSERT_VALUES));
32: }
33: i = 1;
34: col[0] = 1;
35: col[1] = 2;
36: PetscCall(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: PetscCall(MatSetValues(A, 1, &i, 2, col, v, INSERT_VALUES));
41: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
42: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
44: for (test = 0; test < 2; test++) {
45: PetscCall(MatProductCreate(A, A, NULL, &B));
47: if (test == 0) {
48: /* Compute B = A*A; B misses 0-th diagonal */
49: PetscCall(MatProductSetType(B, MATPRODUCT_AB));
50: PetscCall(MatSetOptionsPrefix(B, "AB_"));
51: } else {
52: /* Compute B = A^t*A; B misses 0-th diagonal */
53: PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
54: PetscCall(MatSetOptionsPrefix(B, "AtB_"));
55: }
57: /* Force allocate missing diagonal entries of B */
58: PetscCall(MatSetOption(B, MAT_FORCE_DIAGONAL_ENTRIES, PETSC_TRUE));
59: PetscCall(MatProductSetFromOptions(B));
61: PetscCall(MatProductSymbolic(B));
62: PetscCall(MatProductNumeric(B));
64: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
66: /* Insert entries to diagonal of B */
67: PetscCall(MatCreateVecs(B, NULL, &diag));
68: PetscCall(MatGetDiagonal(B, diag));
69: PetscCall(VecSetValue(diag, 0, 100.0, INSERT_VALUES));
70: PetscCall(VecAssemblyBegin(diag));
71: PetscCall(VecAssemblyEnd(diag));
73: PetscCall(MatDiagonalSet(B, diag, INSERT_VALUES));
74: if (test == 1) PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
75: PetscCall(MatDestroy(&B));
76: PetscCall(VecDestroy(&diag));
77: }
79: PetscCall(MatDestroy(&A));
80: PetscCall(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*/