1c4762a1bSJed Brown static char help[] = "Example of inverting a block diagonal matrix.\n" 2c4762a1bSJed Brown "\n"; 3c4762a1bSJed Brown 4c4762a1bSJed Brown #include <petscmat.h> 5c4762a1bSJed Brown 6c4762a1bSJed Brown /*T 7c4762a1bSJed Brown Concepts: Mat 8c4762a1bSJed Brown T*/ 9c4762a1bSJed Brown 10c4762a1bSJed Brown int main(int argc, char **args) 11c4762a1bSJed Brown { 12c4762a1bSJed Brown Mat A,A_inv; 13c4762a1bSJed Brown PetscMPIInt rank,size; 14c4762a1bSJed Brown PetscInt M,m,bs,rstart,rend,j,x,y; 15c4762a1bSJed Brown PetscInt* dnnz; 16c4762a1bSJed Brown PetscErrorCode ierr; 17c4762a1bSJed Brown PetscScalar *v; 18c4762a1bSJed Brown Vec X, Y; 19c4762a1bSJed Brown PetscReal norm; 20c4762a1bSJed Brown 21*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscInitialize(&argc,&args,(char*)0,help)); 225f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 235f80ce2aSJacob Faibussowitsch CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 24c4762a1bSJed Brown 25c4762a1bSJed Brown ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex184","Mat");CHKERRQ(ierr); 26c4762a1bSJed Brown M=8; 275f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-mat_size",&M,NULL)); 28c4762a1bSJed Brown bs=3; 295f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL)); 30c4762a1bSJed Brown ierr = PetscOptionsEnd();CHKERRQ(ierr); 31c4762a1bSJed Brown 325f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreate(PETSC_COMM_WORLD, &A)); 335f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetFromOptions(A)); 345f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M*bs,M*bs)); 355f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetBlockSize(A,bs)); 365f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetUp(A)); 375f80ce2aSJacob Faibussowitsch CHKERRQ(MatGetLocalSize(A,&m,NULL)); 385f80ce2aSJacob Faibussowitsch CHKERRQ(PetscMalloc1(m/bs,&dnnz)); 39c4762a1bSJed Brown for (j = 0; j < m/bs; j++) { 40c4762a1bSJed Brown dnnz[j] = 1; 41c4762a1bSJed Brown } 425f80ce2aSJacob Faibussowitsch CHKERRQ(MatXAIJSetPreallocation(A,bs,dnnz,NULL,NULL,NULL)); 435f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFree(dnnz)); 44c4762a1bSJed Brown 455f80ce2aSJacob Faibussowitsch CHKERRQ(PetscMalloc1(bs*bs,&v)); 465f80ce2aSJacob Faibussowitsch CHKERRQ(MatGetOwnershipRange(A,&rstart,&rend)); 47c4762a1bSJed Brown for (j = rstart/bs; j < rend/bs; j++) { 48c4762a1bSJed Brown for (x = 0; x < bs; x++) { 49c4762a1bSJed Brown for (y = 0; y < bs; y++) { 50c4762a1bSJed Brown if (x == y) { 51c4762a1bSJed Brown v[y+bs*x] = 2*bs; 52c4762a1bSJed Brown } else { 53c4762a1bSJed Brown v[y+bs*x] = -1 * (x < y) - 2 * (x > y); 54c4762a1bSJed Brown } 55c4762a1bSJed Brown } 56c4762a1bSJed Brown } 575f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValuesBlocked(A,1,&j,1,&j,v,INSERT_VALUES)); 58c4762a1bSJed Brown } 595f80ce2aSJacob Faibussowitsch CHKERRQ(PetscFree(v)); 605f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 615f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 62c4762a1bSJed Brown 63c4762a1bSJed Brown /* check that A = inv(inv(A)) */ 645f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A_inv)); 655f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetFromOptions(A_inv)); 665f80ce2aSJacob Faibussowitsch CHKERRQ(MatInvertBlockDiagonalMat(A,A_inv)); 67c4762a1bSJed Brown 68c4762a1bSJed Brown /* Test A_inv * A on a random vector */ 695f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreateVecs(A, &X, &Y)); 705f80ce2aSJacob Faibussowitsch CHKERRQ(VecSetRandom(X, NULL)); 715f80ce2aSJacob Faibussowitsch CHKERRQ(MatMult(A, X, Y)); 725f80ce2aSJacob Faibussowitsch CHKERRQ(VecScale(X, -1)); 735f80ce2aSJacob Faibussowitsch CHKERRQ(MatMultAdd(A_inv, Y, X, X)); 745f80ce2aSJacob Faibussowitsch CHKERRQ(VecNorm(X, NORM_MAX, &norm)); 75c4762a1bSJed Brown if (norm > PETSC_SMALL) { 765f80ce2aSJacob Faibussowitsch CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Norm of error exceeds tolerance.\nInverse of block diagonal A\n")); 775f80ce2aSJacob Faibussowitsch CHKERRQ(MatView(A_inv,PETSC_VIEWER_STDOUT_WORLD)); 78c4762a1bSJed Brown } 79c4762a1bSJed Brown 805f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&A)); 815f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&A_inv)); 825f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&X)); 835f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&Y)); 84c4762a1bSJed Brown 85*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscFinalize()); 86*b122ec5aSJacob Faibussowitsch return 0; 87c4762a1bSJed Brown } 88c4762a1bSJed Brown 89c4762a1bSJed Brown /*TEST 90c4762a1bSJed Brown test: 91c4762a1bSJed Brown suffix: seqaij 92c4762a1bSJed Brown args: -mat_type seqaij -mat_size 12 -mat_block_size 3 93c4762a1bSJed Brown nsize: 1 94c4762a1bSJed Brown test: 95c4762a1bSJed Brown suffix: seqbaij 96c4762a1bSJed Brown args: -mat_type seqbaij -mat_size 12 -mat_block_size 3 97c4762a1bSJed Brown nsize: 1 98c4762a1bSJed Brown test: 99c4762a1bSJed Brown suffix: mpiaij 100c4762a1bSJed Brown args: -mat_type mpiaij -mat_size 12 -mat_block_size 3 101c4762a1bSJed Brown nsize: 2 102c4762a1bSJed Brown test: 103c4762a1bSJed Brown suffix: mpibaij 104c4762a1bSJed Brown args: -mat_type mpibaij -mat_size 12 -mat_block_size 3 105c4762a1bSJed Brown nsize: 2 106c4762a1bSJed Brown TEST*/ 107