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 PetscErrorCode ierr; 13 PetscScalar *v; 14 Vec X, Y; 15 PetscReal norm; 16 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 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex184","Mat");PetscCall(ierr); 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 ierr = PetscOptionsEnd();PetscCall(ierr); 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++) { 36 dnnz[j] = 1; 37 } 38 PetscCall(MatXAIJSetPreallocation(A,bs,dnnz,NULL,NULL,NULL)); 39 PetscCall(PetscFree(dnnz)); 40 41 PetscCall(PetscMalloc1(bs*bs,&v)); 42 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 43 for (j = rstart/bs; j < rend/bs; j++) { 44 for (x = 0; x < bs; x++) { 45 for (y = 0; y < bs; y++) { 46 if (x == y) { 47 v[y+bs*x] = 2*bs; 48 } else { 49 v[y+bs*x] = -1 * (x < y) - 2 * (x > y); 50 } 51 } 52 } 53 PetscCall(MatSetValuesBlocked(A,1,&j,1,&j,v,INSERT_VALUES)); 54 } 55 PetscCall(PetscFree(v)); 56 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 57 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 58 59 /* check that A = inv(inv(A)) */ 60 PetscCall(MatCreate(PETSC_COMM_WORLD,&A_inv)); 61 PetscCall(MatSetFromOptions(A_inv)); 62 PetscCall(MatInvertBlockDiagonalMat(A,A_inv)); 63 64 /* Test A_inv * A on a random vector */ 65 PetscCall(MatCreateVecs(A, &X, &Y)); 66 PetscCall(VecSetRandom(X, NULL)); 67 PetscCall(MatMult(A, X, Y)); 68 PetscCall(VecScale(X, -1)); 69 PetscCall(MatMultAdd(A_inv, Y, X, X)); 70 PetscCall(VecNorm(X, NORM_MAX, &norm)); 71 if (norm > PETSC_SMALL) { 72 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error exceeds tolerance.\nInverse of block diagonal A\n")); 73 PetscCall(MatView(A_inv,PETSC_VIEWER_STDOUT_WORLD)); 74 } 75 76 PetscCall(MatDestroy(&A)); 77 PetscCall(MatDestroy(&A_inv)); 78 PetscCall(VecDestroy(&X)); 79 PetscCall(VecDestroy(&Y)); 80 81 PetscCall(PetscFinalize()); 82 return 0; 83 } 84 85 /*TEST 86 test: 87 suffix: seqaij 88 args: -mat_type seqaij -mat_size 12 -mat_block_size 3 89 nsize: 1 90 test: 91 suffix: seqbaij 92 args: -mat_type seqbaij -mat_size 12 -mat_block_size 3 93 nsize: 1 94 test: 95 suffix: mpiaij 96 args: -mat_type mpiaij -mat_size 12 -mat_block_size 3 97 nsize: 2 98 test: 99 suffix: mpibaij 100 args: -mat_type mpibaij -mat_size 12 -mat_block_size 3 101 nsize: 2 102 TEST*/ 103