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;
14: PetscInitialize(&argc, &argv, NULL, help);
15: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
16: 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: ISCreateGeneral(PETSC_COMM_WORLD, 1, ix0a, PETSC_COPY_VALUES, &is0a);
27: ISCreateGeneral(PETSC_COMM_WORLD, 1, ix0b, PETSC_COPY_VALUES, &is0b);
28: ISCreateGeneral(PETSC_COMM_WORLD, 2, ix0, PETSC_COPY_VALUES, &is0);
29: ISCreateGeneral(PETSC_COMM_WORLD, 1, ix1, PETSC_COPY_VALUES, &is1);
30: }
31: {
32: ISCreateStride(PETSC_COMM_SELF, 6, 0, 1, &isl0);
33: ISCreateStride(PETSC_COMM_SELF, 3, 0, 1, &isl0a);
34: ISCreateStride(PETSC_COMM_SELF, 3, 3, 1, &isl0b);
35: ISCreateStride(PETSC_COMM_SELF, 3, 6, 1, &isl1);
36: }
38: usenest = PETSC_FALSE;
39: 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: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, 3, l2gind, PETSC_COPY_VALUES, &l2g);
49: for (i = 0; i < 9; i++) {
50: MatCreateAIJ(PETSC_COMM_WORLD, 1, 1, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, NULL, PETSC_DECIDE, NULL, &B[i]);
51: MatSetUp(B[i]);
52: 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: MatCreateNest(PETSC_COMM_WORLD, 2, isx, 2, isx, Bx00, &B00);
71: MatSetUp(B00);
72: MatCreateNest(PETSC_COMM_WORLD, 2, isx, 1, NULL, Bx01, &B01);
73: MatSetUp(B01);
74: MatCreateNest(PETSC_COMM_WORLD, 1, NULL, 2, isx, Bx10, &B10);
75: 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: MatCreateNest(PETSC_COMM_WORLD, 2, isy, 2, isy, By, &A);
88: MatSetUp(A);
89: }
90: MatDestroy(&B00);
91: MatDestroy(&B01);
92: MatDestroy(&B10);
93: }
94: for (i = 0; i < 9; i++) MatDestroy(&B[i]);
95: 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: ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, 9, l2gind, PETSC_COPY_VALUES, &l2g);
102: MatCreateAIJ(PETSC_COMM_WORLD, 3, 3, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, NULL, PETSC_DECIDE, NULL, &A);
103: MatSetLocalToGlobalMapping(A, l2g, l2g);
104: ISLocalToGlobalMappingDestroy(&l2g);
105: }
107: {
108: Mat A00, A11, A0a0a, A0a0b;
109: MatGetLocalSubMatrix(A, isl0, isl0, &A00);
110: MatGetLocalSubMatrix(A, isl1, isl1, &A11);
111: MatGetLocalSubMatrix(A00, isl0a, isl0a, &A0a0a);
112: MatGetLocalSubMatrix(A00, isl0a, isl0b, &A0a0b);
114: MatSetValueLocal(A0a0a, 0, 0, 100 * rank + 1, ADD_VALUES);
115: MatSetValueLocal(A0a0a, 0, 1, 100 * rank + 2, ADD_VALUES);
116: MatSetValueLocal(A0a0a, 2, 2, 100 * rank + 9, ADD_VALUES);
118: MatSetValueLocal(A0a0b, 1, 1, 100 * rank + 50 + 5, ADD_VALUES);
120: MatSetValueLocal(A11, 0, 0, 1000 * (rank + 1) + 1, ADD_VALUES);
121: MatSetValueLocal(A11, 1, 2, 1000 * (rank + 1) + 6, ADD_VALUES);
123: MatRestoreLocalSubMatrix(A00, isl0a, isl0a, &A0a0a);
124: MatRestoreLocalSubMatrix(A00, isl0a, isl0b, &A0a0b);
125: MatRestoreLocalSubMatrix(A, isl0, isl0, &A00);
126: MatRestoreLocalSubMatrix(A, isl1, isl1, &A11);
127: }
128: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
129: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
131: MatComputeOperator(A, MATAIJ, &Aexplicit);
132: MatView(Aexplicit, PETSC_VIEWER_STDOUT_WORLD);
134: MatDestroy(&A);
135: MatDestroy(&Aexplicit);
136: ISDestroy(&is0a);
137: ISDestroy(&is0b);
138: ISDestroy(&is0);
139: ISDestroy(&is1);
140: ISDestroy(&isl0a);
141: ISDestroy(&isl0b);
142: ISDestroy(&isl0);
143: ISDestroy(&isl1);
144: 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*/