Actual source code: ex184.c
1: static char help[] = "Example of inverting a block diagonal matrix.\n"
2: "\n";
4: #include <petscmat.h>
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;
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));
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();
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));
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));
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));
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: }
74: PetscCall(MatDestroy(&A));
75: PetscCall(MatDestroy(&A_inv));
76: PetscCall(VecDestroy(&X));
77: PetscCall(VecDestroy(&Y));
79: PetscCall(PetscFinalize());
80: return 0;
81: }
83: /*TEST
84: test:
85: suffix: seqaij
86: args: -mat_type seqaij -mat_size 12 -mat_block_size 3
87: nsize: 1
88: test:
89: suffix: seqbaij
90: args: -mat_type seqbaij -mat_size 12 -mat_block_size 3
91: nsize: 1
92: test:
93: suffix: mpiaij
94: args: -mat_type mpiaij -mat_size 12 -mat_block_size 3
95: nsize: 2
96: test:
97: suffix: mpibaij
98: args: -mat_type mpibaij -mat_size 12 -mat_block_size 3
99: nsize: 2
100: TEST*/