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
main(int argc,char ** args)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