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