1 static char help[] = "Example of inverting a block diagonal matrix.\n"
2 "\n";
3
4 #include <petscmat.h>
5
main(int argc,char ** args)6 int main(int argc, char **args)
7 {
8 Mat A, A_inv;
9 PetscMPIInt rank, size;
10 PetscInt M, m, bs, rstart, rend, j, x, y;
11 PetscInt *dnnz;
12 PetscScalar *v;
13 Vec X, Y;
14 PetscReal norm;
15
16 PetscFunctionBeginUser;
17 PetscCall(PetscInitialize(&argc, &args, NULL, help));
18 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
19 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
20
21 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex184", "Mat");
22 M = 8;
23 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &M, NULL));
24 bs = 3;
25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL));
26 PetscOptionsEnd();
27
28 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
29 PetscCall(MatSetFromOptions(A));
30 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M * bs, M * bs));
31 PetscCall(MatSetBlockSize(A, bs));
32 PetscCall(MatSetUp(A)); /* called so that MatGetLocalSize() will work */
33 PetscCall(MatGetLocalSize(A, &m, NULL));
34 PetscCall(PetscMalloc1(m / bs, &dnnz));
35 for (j = 0; j < m / bs; j++) dnnz[j] = 1;
36 PetscCall(MatXAIJSetPreallocation(A, bs, dnnz, NULL, NULL, NULL));
37 PetscCall(PetscFree(dnnz));
38
39 PetscCall(PetscMalloc1(bs * bs, &v));
40 PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
41 for (j = rstart / bs; j < rend / bs; j++) {
42 for (x = 0; x < bs; x++) {
43 for (y = 0; y < bs; y++) {
44 if (x == y) {
45 v[y + bs * x] = 2 * bs;
46 } else {
47 v[y + bs * x] = -1 * (x < y) - 2 * (x > y);
48 }
49 }
50 }
51 PetscCall(MatSetValuesBlocked(A, 1, &j, 1, &j, v, INSERT_VALUES));
52 }
53 PetscCall(PetscFree(v));
54 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
55 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
56
57 /* check that A = inv(inv(A)) */
58 PetscCall(MatCreate(PETSC_COMM_WORLD, &A_inv));
59 PetscCall(MatSetFromOptions(A_inv));
60 PetscCall(MatInvertBlockDiagonalMat(A, A_inv));
61
62 /* Test A_inv * A on a random vector */
63 PetscCall(MatCreateVecs(A, &X, &Y));
64 PetscCall(VecSetRandom(X, NULL));
65 PetscCall(MatMult(A, X, Y));
66 PetscCall(VecScale(X, -1));
67 PetscCall(MatMultAdd(A_inv, Y, X, X));
68 PetscCall(VecNorm(X, NORM_MAX, &norm));
69 if (norm > PETSC_SMALL) {
70 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error exceeds tolerance.\nInverse of block diagonal A\n"));
71 PetscCall(MatView(A_inv, PETSC_VIEWER_STDOUT_WORLD));
72 }
73
74 PetscCall(MatDestroy(&A));
75 PetscCall(MatDestroy(&A_inv));
76 PetscCall(VecDestroy(&X));
77 PetscCall(VecDestroy(&Y));
78
79 PetscCall(PetscFinalize());
80 return 0;
81 }
82
83 /*TEST
84 test:
85 suffix: seqaij
86 args: -mat_type seqaij -mat_size 12 -mat_block_size 3
87 nsize: 1
88 output_file: output/empty.out
89 test:
90 suffix: seqbaij
91 args: -mat_type seqbaij -mat_size 12 -mat_block_size 3
92 nsize: 1
93 output_file: output/empty.out
94 test:
95 suffix: mpiaij
96 args: -mat_type mpiaij -mat_size 12 -mat_block_size 3
97 nsize: 2
98 output_file: output/empty.out
99 test:
100 suffix: mpibaij
101 args: -mat_type mpibaij -mat_size 12 -mat_block_size 3
102 nsize: 2
103 output_file: output/empty.out
104 TEST*/
105