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