Actual source code: ex81.c
1: static char help[] = "Tests MatOption MAT_FORCE_DIAGONAL_ENTRIES.\n\n";
3: #include <petscmat.h>
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];
12: PetscFunctionBeginUser;
13: PetscCall(PetscInitialize(&argc, &args, NULL, help));
14: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
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));
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));
43: for (test = 0; test < 2; test++) {
44: PetscCall(MatProductCreate(A, A, NULL, &B));
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: }
56: /* Force allocate missing diagonal entries of B */
57: PetscCall(MatSetOption(B, MAT_FORCE_DIAGONAL_ENTRIES, PETSC_TRUE));
58: PetscCall(MatProductSetFromOptions(B));
60: PetscCall(MatProductSymbolic(B));
61: PetscCall(MatProductNumeric(B));
63: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
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));
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: }
78: PetscCall(MatDestroy(&A));
79: PetscCall(PetscFinalize());
80: return 0;
81: }
83: /*TEST
85: test:
86: output_file: output/ex81_1.out
88: test:
89: suffix: 2
90: args: -AtB_mat_product_algorithm at*b
91: output_file: output/ex81_1.out
93: test:
94: suffix: 3
95: args: -AtB_mat_product_algorithm outerproduct
96: output_file: output/ex81_1.out
98: test:
99: suffix: 4
100: nsize: 3
101: args: -AtB_mat_product_algorithm nonscalable
102: output_file: output/ex81_3.out
104: test:
105: suffix: 5
106: nsize: 3
107: args: -AtB_mat_product_algorithm scalable
108: output_file: output/ex81_3.out
110: test:
111: suffix: 6
112: nsize: 3
113: args: -AtB_mat_product_algorithm at*b
114: output_file: output/ex81_3.out
116: TEST*/