Actual source code: ex212.c
1: static char help[] = "Test MatCreateSubMatrix with -mat_type nest and block sizes.\n";
3: #include <petscmat.h>
5: int main(int argc, char **argv)
6: {
7: Mat A, B, C, mats[6];
8: IS rows[2];
9: ISLocalToGlobalMapping cmap, rmap;
10: const PetscInt indices[3] = {0, 1, 2};
11: PetscInt i;
12: PetscMPIInt size;
14: PetscFunctionBeginUser;
15: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
16: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
17: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only coded for one process");
18: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
19: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
20: PetscCall(MatSetBlockSizes(A, 2, 1));
21: PetscCall(MatSetType(A, MATAIJ));
22: PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 2, 1, indices, PETSC_COPY_VALUES, &rmap));
23: PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, 1, indices, PETSC_COPY_VALUES, &cmap));
24: PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
25: PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
26: PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
27: PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, 2, indices, PETSC_COPY_VALUES, &rmap));
28: PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, 3, indices, PETSC_COPY_VALUES, &cmap));
29: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
30: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, 2, 3));
31: PetscCall(MatSetBlockSizes(B, 1, 1));
32: PetscCall(MatSetType(B, MATAIJ));
33: PetscCall(MatSetLocalToGlobalMapping(B, rmap, cmap));
34: PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
35: PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
36: PetscCall(MatSetUp(A));
37: PetscCall(MatSetUp(B));
38: mats[0] = A;
39: mats[1] = B;
40: mats[2] = A;
41: mats[3] = NULL;
42: mats[4] = B;
43: mats[5] = A;
44: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 3, NULL, mats, &C));
45: PetscCall(MatSetUp(C));
46: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
47: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
48: PetscCall(MatView(C, NULL));
49: PetscCall(MatNestGetISs(C, rows, NULL));
50: for (i = 0; i < 2; i++) {
51: Mat submat;
52: IS cols[3];
53: PetscInt j;
54: PetscCall(ISView(rows[i], NULL));
55: PetscCall(MatCreateSubMatrix(C, rows[i], NULL, MAT_INITIAL_MATRIX, &submat));
56: PetscCall(MatView(submat, NULL));
57: PetscCall(MatNestGetISs(submat, NULL, cols));
58: for (j = 0; j < 3; j++) PetscCall(ISView(cols[j], NULL));
59: PetscCall(MatDestroy(&submat));
60: }
61: PetscCall(MatDestroy(&C));
62: PetscCall(MatDestroy(&B));
63: PetscCall(MatDestroy(&A));
64: PetscCall(PetscFinalize());
65: return 0;
66: }
68: /*TEST
70: test:
72: TEST*/