Actual source code: ex159.c
1: static const char help[] = "Test MatGetLocalSubMatrix() with multiple levels of nesting.\n\n";
3: #include <petscmat.h>
5: int main(int argc, char *argv[])
6: {
7: IS is0a, is0b, is0, is1, isl0a, isl0b, isl0, isl1;
8: Mat A, Aexplicit;
9: PetscBool usenest;
10: PetscMPIInt rank, size;
11: PetscInt i, j;
13: PetscFunctionBeginUser;
14: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
15: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
16: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
18: {
19: PetscInt ix0a[1], ix0b[1], ix0[2], ix1[1];
21: ix0a[0] = rank * 2 + 0;
22: ix0b[0] = rank * 2 + 1;
23: ix0[0] = rank * 3 + 0;
24: ix0[1] = rank * 3 + 1;
25: ix1[0] = rank * 3 + 2;
26: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 1, ix0a, PETSC_COPY_VALUES, &is0a));
27: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 1, ix0b, PETSC_COPY_VALUES, &is0b));
28: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 2, ix0, PETSC_COPY_VALUES, &is0));
29: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 1, ix1, PETSC_COPY_VALUES, &is1));
30: }
31: {
32: PetscCall(ISCreateStride(PETSC_COMM_SELF, 6, 0, 1, &isl0));
33: PetscCall(ISCreateStride(PETSC_COMM_SELF, 3, 0, 1, &isl0a));
34: PetscCall(ISCreateStride(PETSC_COMM_SELF, 3, 3, 1, &isl0b));
35: PetscCall(ISCreateStride(PETSC_COMM_SELF, 3, 6, 1, &isl1));
36: }
38: usenest = PETSC_FALSE;
39: PetscCall(PetscOptionsGetBool(NULL, NULL, "-nest", &usenest, NULL));
40: if (usenest) {
41: ISLocalToGlobalMapping l2g;
42: PetscInt l2gind[3];
43: Mat B[9];
45: l2gind[0] = (rank - 1 + size) % size;
46: l2gind[1] = rank;
47: l2gind[2] = (rank + 1) % size;
48: PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, 3, l2gind, PETSC_COPY_VALUES, &l2g));
49: for (i = 0; i < 9; i++) {
50: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, 1, 1, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, NULL, PETSC_DECIDE, NULL, &B[i]));
51: PetscCall(MatSetUp(B[i]));
52: PetscCall(MatSetLocalToGlobalMapping(B[i], l2g, l2g));
53: }
54: {
55: IS isx[2];
56: Mat Bx00[4], Bx01[2], Bx10[2];
57: Mat B00, B01, B10;
59: isx[0] = is0a;
60: isx[1] = is0b;
61: Bx00[0] = B[0];
62: Bx00[1] = B[1];
63: Bx00[2] = B[3];
64: Bx00[3] = B[4];
65: Bx01[0] = B[2];
66: Bx01[1] = B[5];
67: Bx10[0] = B[6];
68: Bx10[1] = B[7];
70: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, isx, 2, isx, Bx00, &B00));
71: PetscCall(MatSetUp(B00));
72: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, isx, 1, NULL, Bx01, &B01));
73: PetscCall(MatSetUp(B01));
74: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 1, NULL, 2, isx, Bx10, &B10));
75: PetscCall(MatSetUp(B10));
76: {
77: Mat By[4];
78: IS isy[2];
80: By[0] = B00;
81: By[1] = B01;
82: By[2] = B10;
83: By[3] = B[8];
84: isy[0] = is0;
85: isy[1] = is1;
87: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, isy, 2, isy, By, &A));
88: PetscCall(MatSetUp(A));
89: }
90: PetscCall(MatDestroy(&B00));
91: PetscCall(MatDestroy(&B01));
92: PetscCall(MatDestroy(&B10));
93: }
94: for (i = 0; i < 9; i++) PetscCall(MatDestroy(&B[i]));
95: PetscCall(ISLocalToGlobalMappingDestroy(&l2g));
96: } else {
97: ISLocalToGlobalMapping l2g;
98: PetscInt l2gind[9];
99: for (i = 0; i < 3; i++)
100: for (j = 0; j < 3; j++) l2gind[3 * i + j] = ((rank - 1 + j + size) % size) * 3 + i;
101: PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, 9, l2gind, PETSC_COPY_VALUES, &l2g));
102: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, 3, 3, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, NULL, PETSC_DECIDE, NULL, &A));
103: PetscCall(MatSetLocalToGlobalMapping(A, l2g, l2g));
104: PetscCall(ISLocalToGlobalMappingDestroy(&l2g));
105: }
107: {
108: Mat A00, A11, A0a0a, A0a0b;
109: PetscCall(MatGetLocalSubMatrix(A, isl0, isl0, &A00));
110: PetscCall(MatGetLocalSubMatrix(A, isl1, isl1, &A11));
111: PetscCall(MatGetLocalSubMatrix(A00, isl0a, isl0a, &A0a0a));
112: PetscCall(MatGetLocalSubMatrix(A00, isl0a, isl0b, &A0a0b));
114: PetscCall(MatSetValueLocal(A0a0a, 0, 0, 100 * rank + 1, ADD_VALUES));
115: PetscCall(MatSetValueLocal(A0a0a, 0, 1, 100 * rank + 2, ADD_VALUES));
116: PetscCall(MatSetValueLocal(A0a0a, 2, 2, 100 * rank + 9, ADD_VALUES));
118: PetscCall(MatSetValueLocal(A0a0b, 1, 1, 100 * rank + 50 + 5, ADD_VALUES));
120: PetscCall(MatSetValueLocal(A11, 0, 0, 1000 * (rank + 1) + 1, ADD_VALUES));
121: PetscCall(MatSetValueLocal(A11, 1, 2, 1000 * (rank + 1) + 6, ADD_VALUES));
123: PetscCall(MatRestoreLocalSubMatrix(A00, isl0a, isl0a, &A0a0a));
124: PetscCall(MatRestoreLocalSubMatrix(A00, isl0a, isl0b, &A0a0b));
125: PetscCall(MatRestoreLocalSubMatrix(A, isl0, isl0, &A00));
126: PetscCall(MatRestoreLocalSubMatrix(A, isl1, isl1, &A11));
127: }
128: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
129: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
131: PetscCall(MatComputeOperator(A, MATAIJ, &Aexplicit));
132: PetscCall(MatView(Aexplicit, PETSC_VIEWER_STDOUT_WORLD));
134: PetscCall(MatDestroy(&A));
135: PetscCall(MatDestroy(&Aexplicit));
136: PetscCall(ISDestroy(&is0a));
137: PetscCall(ISDestroy(&is0b));
138: PetscCall(ISDestroy(&is0));
139: PetscCall(ISDestroy(&is1));
140: PetscCall(ISDestroy(&isl0a));
141: PetscCall(ISDestroy(&isl0b));
142: PetscCall(ISDestroy(&isl0));
143: PetscCall(ISDestroy(&isl1));
144: PetscCall(PetscFinalize());
145: return 0;
146: }
148: /*TEST
150: test:
151: nsize: 3
153: test:
154: suffix: nest
155: nsize: 3
156: args: -nest
158: TEST*/