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;
17: PetscInitialize(&argc, &args, (char *)0, help);
18: MPI_Comm_size(PETSC_COMM_WORLD, &size);
19: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
21: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex184", "Mat");
22: M = 8;
23: PetscOptionsGetInt(NULL, NULL, "-mat_size", &M, NULL);
24: bs = 3;
25: PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL);
26: PetscOptionsEnd();
28: MatCreate(PETSC_COMM_WORLD, &A);
29: MatSetFromOptions(A);
30: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M * bs, M * bs);
31: MatSetBlockSize(A, bs);
32: MatSetUp(A);
33: MatGetLocalSize(A, &m, NULL);
34: PetscMalloc1(m / bs, &dnnz);
35: for (j = 0; j < m / bs; j++) dnnz[j] = 1;
36: MatXAIJSetPreallocation(A, bs, dnnz, NULL, NULL, NULL);
37: PetscFree(dnnz);
39: PetscMalloc1(bs * bs, &v);
40: 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: MatSetValuesBlocked(A, 1, &j, 1, &j, v, INSERT_VALUES);
52: }
53: PetscFree(v);
54: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
55: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
57: /* check that A = inv(inv(A)) */
58: MatCreate(PETSC_COMM_WORLD, &A_inv);
59: MatSetFromOptions(A_inv);
60: MatInvertBlockDiagonalMat(A, A_inv);
62: /* Test A_inv * A on a random vector */
63: MatCreateVecs(A, &X, &Y);
64: VecSetRandom(X, NULL);
65: MatMult(A, X, Y);
66: VecScale(X, -1);
67: MatMultAdd(A_inv, Y, X, X);
68: VecNorm(X, NORM_MAX, &norm);
69: if (norm > PETSC_SMALL) {
70: PetscPrintf(PETSC_COMM_WORLD, "Norm of error exceeds tolerance.\nInverse of block diagonal A\n");
71: MatView(A_inv, PETSC_VIEWER_STDOUT_WORLD);
72: }
74: MatDestroy(&A);
75: MatDestroy(&A_inv);
76: VecDestroy(&X);
77: VecDestroy(&Y);
79: 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*/