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