Actual source code: dmmbmat.cxx
1: #include <petsc/private/dmmbimpl.h>
2: #include <petsc/private/vecimpl.h>
4: #include <petscdmmoab.h>
5: #include <MBTagConventions.hpp>
6: #include <moab/NestedRefine.hpp>
8: PETSC_INTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM, PetscInt *, PetscInt *, PetscInt *, PetscInt *, PetscBool);
10: PETSC_INTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J)
11: {
12: PetscInt innz = 0, ionz = 0, nlsiz;
13: DM_Moab *dmmoab = (DM_Moab *)dm->data;
14: PetscInt *nnz = 0, *onz = 0;
15: char *tmp = 0;
16: Mat A;
17: MatType mtype;
19: PetscFunctionBegin;
21: PetscAssertPointer(J, 3);
23: /* next, need to allocate the non-zero arrays to enable pre-allocation */
24: mtype = dm->mattype;
25: PetscCall(PetscStrstr(mtype, MATBAIJ, &tmp));
26: nlsiz = (tmp ? dmmoab->nloc : dmmoab->nloc * dmmoab->numFields);
28: /* allocate the nnz, onz arrays based on block size and local nodes */
29: PetscCall(PetscCalloc2(nlsiz, &nnz, nlsiz, &onz));
31: /* compute the nonzero pattern based on MOAB connectivity data for local elements */
32: PetscCall(DMMoab_Compute_NNZ_From_Connectivity(dm, &innz, nnz, &ionz, onz, tmp ? PETSC_TRUE : PETSC_FALSE));
34: /* create the Matrix and set its type as specified by user */
35: PetscCall(MatCreate(((PetscObject)dm)->comm, &A));
36: PetscCall(MatSetSizes(A, dmmoab->nloc * dmmoab->numFields, dmmoab->nloc * dmmoab->numFields, PETSC_DETERMINE, PETSC_DETERMINE));
37: PetscCall(MatSetType(A, mtype));
38: PetscCall(MatSetBlockSize(A, dmmoab->bs));
39: PetscCall(MatSetDM(A, dm)); /* set DM reference */
40: PetscCall(MatSetFromOptions(A));
42: PetscCheck(dmmoab->ltog_map, ((PetscObject)dm)->comm, PETSC_ERR_ORDER, "Cannot create a DMMoab Mat without calling DMSetUp first.");
43: PetscCall(MatSetLocalToGlobalMapping(A, dmmoab->ltog_map, dmmoab->ltog_map));
45: /* set preallocation based on different supported Mat types */
46: PetscCall(MatSeqAIJSetPreallocation(A, innz, nnz));
47: PetscCall(MatMPIAIJSetPreallocation(A, innz, nnz, ionz, onz));
48: PetscCall(MatSeqBAIJSetPreallocation(A, dmmoab->bs, innz, nnz));
49: PetscCall(MatMPIBAIJSetPreallocation(A, dmmoab->bs, innz, nnz, ionz, onz));
51: /* clean up temporary memory */
52: PetscCall(PetscFree2(nnz, onz));
54: /* set up internal matrix data-structures */
55: PetscCall(MatSetUp(A));
57: /* MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); */
59: *J = A;
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: PETSC_INTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm, PetscInt *innz, PetscInt *nnz, PetscInt *ionz, PetscInt *onz, PetscBool isbaij)
64: {
65: PetscInt i, f, nloc, vpere, bs, n_nnz, n_onz, ivtx = 0;
66: PetscInt ibs, jbs, inbsize, iobsize, nfields, nlsiz;
67: DM_Moab *dmmoab = (DM_Moab *)dm->data;
68: moab::Range found;
69: std::vector<moab::EntityHandle> adjs, storage;
70: PetscBool isinterlaced;
71: moab::EntityHandle vtx;
73: PetscFunctionBegin;
74: bs = dmmoab->bs;
75: nloc = dmmoab->nloc;
76: nfields = dmmoab->numFields;
77: isinterlaced = (isbaij || bs == nfields ? PETSC_TRUE : PETSC_FALSE);
78: nlsiz = (isinterlaced ? nloc : nloc * nfields);
80: /* loop over the locally owned vertices and figure out the NNZ pattern using connectivity information */
81: for (moab::Range::const_iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, ivtx++) {
82: vtx = *iter;
83: /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects
84: to the current vertex. We can then decipher if a vertex is ghosted or not and compute the
85: non-zero pattern accordingly. */
86: adjs.clear();
87: if (dmmoab->hlevel && (dmmoab->pcomm->size() == 1)) PetscCallMOAB(dmmoab->hierarchy->get_adjacencies(vtx, dmmoab->dim, adjs));
88: else PetscCallMOAB(dmmoab->mbiface->get_adjacencies(&vtx, 1, dmmoab->dim, true, adjs, moab::Interface::UNION));
90: /* reset counters */
91: n_nnz = n_onz = 0;
92: found.clear();
94: /* loop over vertices and update the number of connectivity */
95: for (unsigned jter = 0; jter < adjs.size(); ++jter) {
96: /* Get connectivity information in canonical ordering for the local element */
97: const moab::EntityHandle *connect;
98: std::vector<moab::EntityHandle> cconnect;
99: PetscCallMOAB(dmmoab->mbiface->get_connectivity(adjs[jter], connect, vpere, false, &storage));
101: /* loop over each element connected to the adjacent vertex and update as needed */
102: for (i = 0; i < vpere; ++i) {
103: /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */
104: if (connect[i] == vtx || found.find(connect[i]) != found.end()) continue; /* make sure we don't double count shared vertices */
105: if (dmmoab->vghost->find(connect[i]) != dmmoab->vghost->end()) n_onz++; /* update out-of-proc onz */
106: else n_nnz++; /* else local vertex */
107: found.insert(connect[i]);
108: }
109: }
110: storage.clear();
112: if (isbaij) {
113: nnz[ivtx] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */
114: if (onz) onz[ivtx] = n_onz; /* add ghost non-owned nodes */
115: } else { /* AIJ matrices */
116: if (!isinterlaced) {
117: for (f = 0; f < nfields; f++) {
118: nnz[f * nloc + ivtx] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */
119: if (onz) onz[f * nloc + ivtx] = n_onz; /* add ghost non-owned nodes */
120: }
121: } else {
122: for (f = 0; f < nfields; f++) {
123: nnz[nfields * ivtx + f] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */
124: if (onz) onz[nfields * ivtx + f] = n_onz; /* add ghost non-owned nodes */
125: }
126: }
127: }
128: }
130: for (i = 0; i < nlsiz; i++) nnz[i] += 1; /* self count the node */
132: for (ivtx = 0; ivtx < nloc; ivtx++) {
133: if (!isbaij) {
134: for (ibs = 0; ibs < nfields; ibs++) {
135: if (dmmoab->dfill) { /* first address the diagonal block */
136: /* just add up the ints -- easier/faster rather than branching based on "1" */
137: for (jbs = 0, inbsize = 0; jbs < nfields; jbs++) inbsize += dmmoab->dfill[ibs * nfields + jbs];
138: } else inbsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */
139: if (isinterlaced) nnz[ivtx * nfields + ibs] *= inbsize;
140: else nnz[ibs * nloc + ivtx] *= inbsize;
142: if (onz) {
143: if (dmmoab->ofill) { /* next address the off-diagonal block */
144: /* just add up the ints -- easier/faster rather than branching based on "1" */
145: for (jbs = 0, iobsize = 0; jbs < nfields; jbs++) iobsize += dmmoab->dfill[ibs * nfields + jbs];
146: } else iobsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */
147: if (isinterlaced) onz[ivtx * nfields + ibs] *= iobsize;
148: else onz[ibs * nloc + ivtx] *= iobsize;
149: }
150: }
151: } else {
152: /* check if we got overzealous in our nnz and onz computations */
153: nnz[ivtx] = (nnz[ivtx] > dmmoab->nloc ? dmmoab->nloc : nnz[ivtx]);
154: if (onz) onz[ivtx] = (onz[ivtx] > dmmoab->nloc ? dmmoab->nloc : onz[ivtx]);
155: }
156: }
157: /* update innz and ionz based on local maxima */
158: if (innz || ionz) {
159: if (innz) *innz = 0;
160: if (ionz) *ionz = 0;
161: for (i = 0; i < nlsiz; i++) {
162: if (innz && nnz[i] > *innz) *innz = nnz[i];
163: if (ionz && onz && onz[i] > *ionz) *ionz = onz[i];
164: }
165: }
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w, const PetscInt *fill, PetscInt **rfill)
170: {
171: PetscInt i, j, *ifill;
173: PetscFunctionBegin;
174: if (!fill) PetscFunctionReturn(PETSC_SUCCESS);
175: PetscCall(PetscMalloc1(w * w, &ifill));
176: for (i = 0; i < w; i++) {
177: for (j = 0; j < w; j++) ifill[i * w + j] = fill[i * w + j];
178: }
180: *rfill = ifill;
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: /*@C
185: DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem
186: of the matrix returned by `DMCreateMatrix()`.
188: Logically Collective
190: Input Parameters:
191: + dm - the `DMMOAB` object
192: . dfill - the fill pattern in the diagonal block (may be `NULL`, means use dense block)
193: - ofill - the fill pattern in the off-diagonal blocks
195: Level: developer
197: Notes:
198: This only makes sense when you are doing multicomponent problems but using the
199: MPIAIJ matrix format
201: The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries
202: representing coupling and 0 entries for missing coupling. For example
203: .vb
204: dfill[9] = {1, 0, 0,
205: 1, 1, 0,
206: 0, 1, 1}
207: .ve
208: means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
209: itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
210: diagonal block).
212: `DMDASetGetMatrix()` allows you to provide general code for those more complicated nonzero patterns then
213: can be represented in the dfill, ofill format
215: Contributed by\:
216: Glenn Hammond
218: .seealso: `DMCreateMatrix()`, `DMDASetGetMatrix()`, `DMSetMatrixPreallocateOnly()`
219: @*/
220: PetscErrorCode DMMoabSetBlockFills(DM dm, const PetscInt *dfill, const PetscInt *ofill)
221: {
222: DM_Moab *dmmoab = (DM_Moab *)dm->data;
224: PetscFunctionBegin;
226: PetscCall(DMMoabSetBlockFills_Private(dmmoab->numFields, dfill, &dmmoab->dfill));
227: PetscCall(DMMoabSetBlockFills_Private(dmmoab->numFields, ofill, &dmmoab->ofill));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }