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