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 6d71ae5a4SJacob Faibussowitsch int main(int argc, char **args) 7d71ae5a4SJacob Faibussowitsch { 8c4762a1bSJed Brown Mat A, A_inv; 9c4762a1bSJed Brown PetscMPIInt rank, size; 10c4762a1bSJed Brown PetscInt M, m, bs, rstart, rend, j, x, y; 11c4762a1bSJed Brown PetscInt *dnnz; 12c4762a1bSJed Brown PetscScalar *v; 13c4762a1bSJed Brown Vec X, Y; 14c4762a1bSJed Brown PetscReal norm; 15c4762a1bSJed Brown 16327415f7SBarry Smith PetscFunctionBeginUser; 17c8025a54SPierre Jolivet PetscCall(PetscInitialize(&argc, &args, NULL, help)); 189566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 199566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 20c4762a1bSJed Brown 21d0609cedSBarry Smith PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex184", "Mat"); 22c4762a1bSJed Brown M = 8; 239566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &M, NULL)); 24c4762a1bSJed Brown bs = 3; 259566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL)); 26d0609cedSBarry Smith PetscOptionsEnd(); 27c4762a1bSJed Brown 289566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 299566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(A)); 309566063dSJacob Faibussowitsch PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M * bs, M * bs)); 319566063dSJacob Faibussowitsch PetscCall(MatSetBlockSize(A, bs)); 32ad79cf63SBarry Smith PetscCall(MatSetUp(A)); /* called so that MatGetLocalSize() will work */ 339566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(A, &m, NULL)); 349566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(m / bs, &dnnz)); 35ad540459SPierre Jolivet for (j = 0; j < m / bs; j++) dnnz[j] = 1; 369566063dSJacob Faibussowitsch PetscCall(MatXAIJSetPreallocation(A, bs, dnnz, NULL, NULL, NULL)); 379566063dSJacob Faibussowitsch PetscCall(PetscFree(dnnz)); 38c4762a1bSJed Brown 399566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(bs * bs, &v)); 409566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 41c4762a1bSJed Brown for (j = rstart / bs; j < rend / bs; j++) { 42c4762a1bSJed Brown for (x = 0; x < bs; x++) { 43c4762a1bSJed Brown for (y = 0; y < bs; y++) { 44c4762a1bSJed Brown if (x == y) { 45c4762a1bSJed Brown v[y + bs * x] = 2 * bs; 46c4762a1bSJed Brown } else { 47c4762a1bSJed Brown v[y + bs * x] = -1 * (x < y) - 2 * (x > y); 48c4762a1bSJed Brown } 49c4762a1bSJed Brown } 50c4762a1bSJed Brown } 519566063dSJacob Faibussowitsch PetscCall(MatSetValuesBlocked(A, 1, &j, 1, &j, v, INSERT_VALUES)); 52c4762a1bSJed Brown } 539566063dSJacob Faibussowitsch PetscCall(PetscFree(v)); 549566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 559566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 56c4762a1bSJed Brown 57c4762a1bSJed Brown /* check that A = inv(inv(A)) */ 589566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &A_inv)); 599566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(A_inv)); 609566063dSJacob Faibussowitsch PetscCall(MatInvertBlockDiagonalMat(A, A_inv)); 61c4762a1bSJed Brown 62c4762a1bSJed Brown /* Test A_inv * A on a random vector */ 639566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(A, &X, &Y)); 649566063dSJacob Faibussowitsch PetscCall(VecSetRandom(X, NULL)); 659566063dSJacob Faibussowitsch PetscCall(MatMult(A, X, Y)); 669566063dSJacob Faibussowitsch PetscCall(VecScale(X, -1)); 679566063dSJacob Faibussowitsch PetscCall(MatMultAdd(A_inv, Y, X, X)); 689566063dSJacob Faibussowitsch PetscCall(VecNorm(X, NORM_MAX, &norm)); 69c4762a1bSJed Brown if (norm > PETSC_SMALL) { 709566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error exceeds tolerance.\nInverse of block diagonal A\n")); 719566063dSJacob Faibussowitsch PetscCall(MatView(A_inv, PETSC_VIEWER_STDOUT_WORLD)); 72c4762a1bSJed Brown } 73c4762a1bSJed Brown 749566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 759566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A_inv)); 769566063dSJacob Faibussowitsch PetscCall(VecDestroy(&X)); 779566063dSJacob Faibussowitsch PetscCall(VecDestroy(&Y)); 78c4762a1bSJed Brown 799566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 80b122ec5aSJacob Faibussowitsch return 0; 81c4762a1bSJed Brown } 82c4762a1bSJed Brown 83c4762a1bSJed Brown /*TEST 84c4762a1bSJed Brown test: 85c4762a1bSJed Brown suffix: seqaij 86c4762a1bSJed Brown args: -mat_type seqaij -mat_size 12 -mat_block_size 3 87c4762a1bSJed Brown nsize: 1 88*3886731fSPierre Jolivet output_file: output/empty.out 89c4762a1bSJed Brown test: 90c4762a1bSJed Brown suffix: seqbaij 91c4762a1bSJed Brown args: -mat_type seqbaij -mat_size 12 -mat_block_size 3 92c4762a1bSJed Brown nsize: 1 93*3886731fSPierre Jolivet output_file: output/empty.out 94c4762a1bSJed Brown test: 95c4762a1bSJed Brown suffix: mpiaij 96c4762a1bSJed Brown args: -mat_type mpiaij -mat_size 12 -mat_block_size 3 97c4762a1bSJed Brown nsize: 2 98*3886731fSPierre Jolivet output_file: output/empty.out 99c4762a1bSJed Brown test: 100c4762a1bSJed Brown suffix: mpibaij 101c4762a1bSJed Brown args: -mat_type mpibaij -mat_size 12 -mat_block_size 3 102c4762a1bSJed Brown nsize: 2 103*3886731fSPierre Jolivet output_file: output/empty.out 104c4762a1bSJed Brown TEST*/ 105