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_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM, PetscInt*, PetscInt*, PetscInt*, PetscInt*, PetscBool);
10: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J)
11: {
12: PetscErrorCode ierr;
13: PetscInt innz = 0, ionz = 0, nlsiz;
14: DM_Moab *dmmoab = (DM_Moab*)dm->data;
15: PetscInt *nnz = 0, *onz = 0;
16: char *tmp = 0;
17: Mat A;
18: MatType mtype;
24: /* next, need to allocate the non-zero arrays to enable pre-allocation */
25: mtype = dm->mattype;
26: PetscStrstr(mtype, MATBAIJ, &tmp);
27: nlsiz = (tmp ? dmmoab->nloc : dmmoab->nloc * dmmoab->numFields);
29: /* allocate the nnz, onz arrays based on block size and local nodes */
30: PetscCalloc2(nlsiz, &nnz, nlsiz, &onz);
32: /* compute the nonzero pattern based on MOAB connectivity data for local elements */
33: DMMoab_Compute_NNZ_From_Connectivity(dm, &innz, nnz, &ionz, onz, (tmp ? PETSC_TRUE : PETSC_FALSE));
35: /* create the Matrix and set its type as specified by user */
36: MatCreate((((PetscObject)dm)->comm), &A);
37: MatSetSizes(A, dmmoab->nloc * dmmoab->numFields, dmmoab->nloc * dmmoab->numFields, PETSC_DETERMINE, PETSC_DETERMINE);
38: MatSetType(A, mtype);
39: MatSetBlockSize(A, dmmoab->bs);
40: MatSetDM(A, dm); /* set DM reference */
41: MatSetFromOptions(A);
43: if (!dmmoab->ltog_map) SETERRQ((((PetscObject)dm)->comm), PETSC_ERR_ORDER, "Cannot create a DMMoab Mat without calling DMSetUp first.");
44: MatSetLocalToGlobalMapping(A, dmmoab->ltog_map, dmmoab->ltog_map);
46: /* set preallocation based on different supported Mat types */
47: MatSeqAIJSetPreallocation(A, innz, nnz);
48: MatMPIAIJSetPreallocation(A, innz, nnz, ionz, onz);
49: MatSeqBAIJSetPreallocation(A, dmmoab->bs, innz, nnz);
50: MatMPIBAIJSetPreallocation(A, dmmoab->bs, innz, nnz, ionz, onz);
52: /* clean up temporary memory */
53: PetscFree2(nnz, onz);
55: /* set up internal matrix data-structures */
56: MatSetUp(A);
58: /* MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); */
60: *J = A;
61: return(0);
62: }
65: PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm, PetscInt* innz, PetscInt* nnz, PetscInt* ionz, PetscInt* onz, PetscBool isbaij)
66: {
67: PetscInt i, f, nloc, vpere, bs, n_nnz, n_onz, ivtx = 0;
68: PetscInt ibs, jbs, inbsize, iobsize, nfields, nlsiz;
69: DM_Moab *dmmoab = (DM_Moab*)dm->data;
70: moab::Range found;
71: std::vector<moab::EntityHandle> adjs, storage;
72: PetscBool isinterlaced;
73: moab::EntityHandle vtx;
74: moab::ErrorCode merr;
77: bs = dmmoab->bs;
78: nloc = dmmoab->nloc;
79: nfields = dmmoab->numFields;
80: isinterlaced = (isbaij || bs == nfields ? PETSC_TRUE : PETSC_FALSE);
81: nlsiz = (isinterlaced ? nloc : nloc * nfields);
83: /* loop over the locally owned vertices and figure out the NNZ pattern using connectivity information */
84: for (moab::Range::const_iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, ivtx++) {
86: vtx = *iter;
87: /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects
88: to the current vertex. We can then decipher if a vertex is ghosted or not and compute the
89: non-zero pattern accordingly. */
90: adjs.clear();
91: if (dmmoab->hlevel && (dmmoab->pcomm->size() == 1)) {
92: merr = dmmoab->hierarchy->get_adjacencies(vtx, dmmoab->dim, adjs); MBERRNM(merr);
93: }
94: else {
95: merr = dmmoab->mbiface->get_adjacencies(&vtx, 1, dmmoab->dim, true, adjs, moab::Interface::UNION); MBERRNM(merr);
96: }
98: /* reset counters */
99: n_nnz = n_onz = 0;
100: found.clear();
102: /* loop over vertices and update the number of connectivity */
103: for (unsigned jter = 0; jter < adjs.size(); ++jter) {
105: /* Get connectivity information in canonical ordering for the local element */
106: const moab::EntityHandle *connect;
107: std::vector<moab::EntityHandle> cconnect;
108: merr = dmmoab->mbiface->get_connectivity(adjs[jter], connect, vpere, false, &storage); MBERRNM(merr);
110: /* loop over each element connected to the adjacent vertex and update as needed */
111: for (i = 0; i < vpere; ++i) {
112: /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */
113: if (connect[i] == vtx || found.find(connect[i]) != found.end()) continue; /* make sure we don't double count shared vertices */
114: if (dmmoab->vghost->find(connect[i]) != dmmoab->vghost->end()) n_onz++; /* update out-of-proc onz */
115: else n_nnz++; /* else local vertex */
116: found.insert(connect[i]);
117: }
118: }
119: storage.clear();
121: if (isbaij) {
122: nnz[ivtx] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */
123: if (onz) {
124: onz[ivtx] = n_onz; /* add ghost non-owned nodes */
125: }
126: }
127: else { /* AIJ matrices */
128: if (!isinterlaced) {
129: for (f = 0; f < nfields; f++) {
130: nnz[f * nloc + ivtx] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */
131: if (onz)
132: onz[f * nloc + ivtx] = n_onz; /* add ghost non-owned nodes */
133: }
134: }
135: else {
136: for (f = 0; f < nfields; f++) {
137: nnz[nfields * ivtx + f] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */
138: if (onz)
139: onz[nfields * ivtx + f] = n_onz; /* add ghost non-owned nodes */
140: }
141: }
142: }
143: }
145: for (i = 0; i < nlsiz; i++)
146: nnz[i] += 1; /* self count the node */
148: for (ivtx = 0; ivtx < nloc; ivtx++) {
149: if (!isbaij) {
150: for (ibs = 0; ibs < nfields; ibs++) {
151: if (dmmoab->dfill) { /* first address the diagonal block */
152: /* just add up the ints -- easier/faster rather than branching based on "1" */
153: for (jbs = 0, inbsize = 0; jbs < nfields; jbs++)
154: inbsize += dmmoab->dfill[ibs * nfields + jbs];
155: }
156: else inbsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */
157: if (isinterlaced) nnz[ivtx * nfields + ibs] *= inbsize;
158: else nnz[ibs * nloc + ivtx] *= inbsize;
160: if (onz) {
161: if (dmmoab->ofill) { /* next address the off-diagonal block */
162: /* just add up the ints -- easier/faster rather than branching based on "1" */
163: for (jbs = 0, iobsize = 0; jbs < nfields; jbs++)
164: iobsize += dmmoab->dfill[ibs * nfields + jbs];
165: }
166: else iobsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */
167: if (isinterlaced) onz[ivtx * nfields + ibs] *= iobsize;
168: else onz[ibs * nloc + ivtx] *= iobsize;
169: }
170: }
171: }
172: else {
173: /* check if we got overzealous in our nnz and onz computations */
174: nnz[ivtx] = (nnz[ivtx] > dmmoab->nloc ? dmmoab->nloc : nnz[ivtx]);
175: if (onz) onz[ivtx] = (onz[ivtx] > dmmoab->nloc ? dmmoab->nloc : onz[ivtx]);
176: }
177: }
178: /* update innz and ionz based on local maxima */
179: if (innz || ionz) {
180: if (innz) *innz = 0;
181: if (ionz) *ionz = 0;
182: for (i = 0; i < nlsiz; i++) {
183: if (innz && (nnz[i] > *innz)) *innz = nnz[i];
184: if ((ionz && onz) && (onz[i] > *ionz)) *ionz = onz[i];
185: }
186: }
187: return(0);
188: }
191: static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w, const PetscInt *fill, PetscInt **rfill)
192: {
194: PetscInt i, j, *ifill;
197: if (!fill) return(0);
198: PetscMalloc1(w * w, &ifill);
199: for (i = 0; i < w; i++) {
200: for (j = 0; j < w; j++)
201: ifill[i * w + j] = fill[i * w + j];
202: }
204: *rfill = ifill;
205: return(0);
206: }
209: /*@C
210: DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem
211: of the matrix returned by DMCreateMatrix().
213: Logically Collective on da
215: Input Parameter:
216: + dm - the DMMoab object
217: . dfill - the fill pattern in the diagonal block (may be NULL, means use dense block)
218: - ofill - the fill pattern in the off-diagonal blocks
220: Level: developer
222: Notes:
223: This only makes sense when you are doing multicomponent problems but using the
224: MPIAIJ matrix format
226: The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries
227: representing coupling and 0 entries for missing coupling. For example
228: $ dfill[9] = {1, 0, 0,
229: $ 1, 1, 0,
230: $ 0, 1, 1}
231: means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
232: itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
233: diagonal block).
235: DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then
236: can be represented in the dfill, ofill format
238: Contributed by Glenn Hammond
240: .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly()
242: @*/
243: PetscErrorCode DMMoabSetBlockFills(DM dm, const PetscInt *dfill, const PetscInt *ofill)
244: {
245: DM_Moab *dmmoab = (DM_Moab*)dm->data;
250: DMMoabSetBlockFills_Private(dmmoab->numFields, dfill, &dmmoab->dfill);
251: DMMoabSetBlockFills_Private(dmmoab->numFields, ofill, &dmmoab->ofill);
252: return(0);
253: }