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