1 static char help[] = "Example of inverting a block diagonal matrix.\n" 2 "\n"; 3 4 #include <petscmat.h> 5 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, (char *)0, 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)); 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 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*/ 101