Actual source code: ex317.c
1: static char help[] = "Tests MatGetLocalSubMatrix() on a MATIS with fields of different block sizes.\n\n";
3: #include <petscmat.h>
5: /* Deterministic entry k of the element matrix coupling node i of field fi with node j of field fj */
6: static PetscScalar ElementValue(PetscInt fi, PetscInt i, PetscInt fj, PetscInt j, PetscInt k)
7: {
8: return (PetscScalar)(1 + ((3 * fi + 5 * i + 7 * fj + 11 * j + 13 * k) % 17));
9: }
11: /*
12: The contract MatGetLocalSubMatrix() owes its caller, which MATIS and MATLOCALREF both honour: the
13: submatrix is sized by its index sets and carries their block size, on the matrix and on its local to
14: global maps alike. Which rows and columns those maps reach is left to the assembly below, where the same
15: operator is built twice through them, a node at a time and a degree of freedom at a time.
16: */
17: static PetscErrorCode CheckLocalSubMatrix(Mat sub, IS row, IS col)
18: {
19: ISLocalToGlobalMapping rl2g, cl2g;
20: PetscInt nrl, ncl, rbs, cbs, m, n, mbs, nbs;
22: PetscFunctionBeginUser;
23: PetscCall(ISGetLocalSize(row, &nrl));
24: PetscCall(ISGetLocalSize(col, &ncl));
25: PetscCall(ISGetBlockSize(row, &rbs));
26: PetscCall(ISGetBlockSize(col, &cbs));
27: PetscCall(MatGetLocalSize(sub, &m, &n));
28: PetscCheck(m == nrl && n == ncl, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Submatrix has local size %" PetscInt_FMT " x %" PetscInt_FMT ", its index sets have %" PetscInt_FMT " x %" PetscInt_FMT, m, n, nrl, ncl);
29: PetscCall(MatGetBlockSizes(sub, &mbs, &nbs));
30: PetscCheck(mbs == rbs && nbs == cbs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Submatrix has block size %" PetscInt_FMT " x %" PetscInt_FMT ", its index sets have %" PetscInt_FMT " x %" PetscInt_FMT, mbs, nbs, rbs, cbs);
31: PetscCall(MatGetLocalToGlobalMapping(sub, &rl2g, &cl2g));
32: PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g, &mbs));
33: PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g, &nbs));
34: PetscCheck(mbs == rbs && nbs == cbs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Maps of the submatrix have block size %" PetscInt_FMT " x %" PetscInt_FMT ", its index sets have %" PetscInt_FMT " x %" PetscInt_FMT, mbs, nbs, rbs, cbs);
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: /*
39: Assemble a two field operator through the local submatrices of a MATIS.
40: With blocked, each field is addressed one node at a time, through index sets that carry its block
41: size; otherwise every degree of freedom is addressed individually. Both paths insert the same values.
42: */
43: static PetscErrorCode AssembleMixed(Mat A, const PetscInt bs[], const PetscInt nn[], const PetscInt off[], PetscBool blocked)
44: {
45: IS is[2];
46: Mat sub[2][2];
47: PetscScalar *vals;
48: PetscInt *rows, *cols;
49: PetscInt fi, fj, i, j, k, mbs;
51: PetscFunctionBeginUser;
52: mbs = PetscMax(bs[0], bs[1]);
53: PetscCall(PetscMalloc3(mbs, &rows, mbs, &cols, mbs * mbs, &vals));
54: for (fi = 0; fi < 2; fi++) {
55: PetscCall(ISCreateStride(PETSC_COMM_SELF, nn[fi] * bs[fi], off[fi], 1, &is[fi]));
56: if (blocked) PetscCall(ISSetBlockSize(is[fi], bs[fi]));
57: }
58: for (fi = 0; fi < 2; fi++) {
59: for (fj = 0; fj < 2; fj++) {
60: PetscCall(MatGetLocalSubMatrix(A, is[fi], is[fj], &sub[fi][fj]));
61: PetscCall(CheckLocalSubMatrix(sub[fi][fj], is[fi], is[fj]));
62: }
63: }
65: /* every field couples its local nodes in a chain, and node i of a field to node i of the other one */
66: for (fi = 0; fi < 2; fi++) {
67: for (fj = 0; fj < 2; fj++) {
68: for (i = 0; i < nn[fi]; i++) {
69: for (j = PetscMax(i - 1, 0); j <= PetscMin(fi == fj ? i + 1 : i, nn[fj] - 1); j++) {
70: for (k = 0; k < bs[fi] * bs[fj]; k++) vals[k] = ElementValue(fi, i, fj, j, k);
71: /* the blocked path reaches the diagonal blocks by node and the off diagonal ones by degree of freedom */
72: if (blocked && fi == fj) PetscCall(MatSetValuesBlockedLocal(sub[fi][fj], 1, &i, 1, &j, vals, ADD_VALUES));
73: else {
74: for (k = 0; k < bs[fi]; k++) rows[k] = i * bs[fi] + k;
75: for (k = 0; k < bs[fj]; k++) cols[k] = j * bs[fj] + k;
76: PetscCall(MatSetValuesLocal(sub[fi][fj], bs[fi], rows, bs[fj], cols, vals, ADD_VALUES));
77: }
78: }
79: }
80: }
81: }
83: for (fi = 0; fi < 2; fi++) {
84: for (fj = 0; fj < 2; fj++) PetscCall(MatRestoreLocalSubMatrix(A, is[fi], is[fj], &sub[fi][fj]));
85: }
86: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
87: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
88: for (fi = 0; fi < 2; fi++) PetscCall(ISDestroy(&is[fi]));
89: PetscCall(PetscFree3(rows, cols, vals));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: int main(int argc, char **args)
94: {
95: Mat A, B, A2, B2;
96: ISLocalToGlobalMapping map;
97: PetscReal nrm, ref;
98: PetscInt *idxs, bs[2], nn[2], off[2], goff[2], gn, n, N, nl = 4, pbs, f, i, j;
99: PetscMPIInt rank, size;
101: PetscFunctionBeginUser;
102: PetscCall(PetscInitialize(&argc, &args, NULL, help));
103: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
104: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
105: bs[0] = 2;
106: bs[1] = 3;
107: PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs0", &bs[0], NULL));
108: PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs1", &bs[1], NULL));
109: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nl", &nl, NULL));
110: PetscCheck(nl > 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Need at least 3 local nodes per subdomain");
112: /*
113: Two fields laid out as in a mixed finite element space: neighboring subdomains share one node of
114: each field, and the second field starts wherever the first one ends. Neither that offset nor the
115: global size need be a multiple of the block size of a field, which is what its local submatrices
116: have to cope with.
117: */
118: nn[0] = nl;
119: nn[1] = nl - 1;
120: off[0] = 0;
121: off[1] = nn[0] * bs[0];
122: goff[0] = 0;
123: goff[1] = (size * (nn[0] - 1) + 1) * bs[0];
124: n = off[1] + nn[1] * bs[1];
125: N = goff[1] + (size * (nn[1] - 1) + 1) * bs[1];
126: PetscCall(PetscMalloc1(n, &idxs));
127: for (f = 0; f < 2; f++) {
128: for (i = 0; i < nn[f]; i++) {
129: gn = rank * (nn[f] - 1) + i;
130: for (j = 0; j < bs[f]; j++) idxs[off[f] + i * bs[f] + j] = goff[f] + gn * bs[f] + j;
131: }
132: }
133: /*
134: The whole space is blocked only when the two fields agree on a block size, and then the local
135: submatrices can be reached through the ordinary blocked maps; with fields of different block sizes
136: they cannot, and the two representations of a local submatrix are both exercised by the tests below.
137: */
138: pbs = bs[0] == bs[1] ? bs[0] : 1;
139: for (i = 0; i < n / pbs; i++) idxs[i] = idxs[i * pbs] / pbs;
140: PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, pbs, n / pbs, idxs, PETSC_OWN_POINTER, &map));
142: PetscCall(MatCreateIS(PETSC_COMM_WORLD, pbs, PETSC_DECIDE, PETSC_DECIDE, N, N, map, map, &A));
143: PetscCall(MatISSetPreallocation(A, 3 * (bs[0] + bs[1]), NULL, 3 * (bs[0] + bs[1]), NULL));
144: PetscCall(MatCreateIS(PETSC_COMM_WORLD, pbs, PETSC_DECIDE, PETSC_DECIDE, N, N, map, map, &B));
145: PetscCall(MatISSetPreallocation(B, 3 * (bs[0] + bs[1]), NULL, 3 * (bs[0] + bs[1]), NULL));
146: PetscCall(AssembleMixed(A, bs, nn, off, PETSC_TRUE));
147: PetscCall(AssembleMixed(B, bs, nn, off, PETSC_FALSE));
149: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &A2));
150: PetscCall(MatConvert(B, MATAIJ, MAT_INITIAL_MATRIX, &B2));
151: PetscCall(MatNorm(B2, NORM_INFINITY, &ref));
152: PetscCall(MatAXPY(A2, -1.0, B2, DIFFERENT_NONZERO_PATTERN));
153: PetscCall(MatNorm(A2, NORM_INFINITY, &nrm));
154: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Reference operator is nonzero: %s\n", ref > 0.0 ? "yes" : "no"));
155: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Blocked and scalar assembly agree: %s\n", nrm < PETSC_SMALL * ref ? "yes" : "no"));
157: PetscCall(MatDestroy(&A));
158: PetscCall(MatDestroy(&B));
159: PetscCall(MatDestroy(&A2));
160: PetscCall(MatDestroy(&B2));
161: PetscCall(ISLocalToGlobalMappingDestroy(&map));
162: PetscCall(PetscFinalize());
163: return 0;
164: }
166: /*TEST
168: test:
169: suffix: 1
170: nsize: {{1 3}}
171: args: -bs0 {{1 2 3}} -bs1 {{1 2 3}}
172: output_file: output/ex317_1.out
174: TEST*/