Actual source code: packm.c
2: #include <../src/dm/impls/composite/packimpl.h>
4: static PetscErrorCode DMCreateMatrix_Composite_Nest(DM dm, Mat *J)
5: {
6: const DM_Composite *com = (DM_Composite *)dm->data;
7: const struct DMCompositeLink *rlink, *clink;
8: IS *isg;
9: Mat *submats;
10: PetscInt i, j, n;
12: n = com->nDM; /* Total number of entries */
14: /* Explicit index sets are not required for MatCreateNest, but getting them here allows MatNest to do consistency
15: * checking and allows ISEqual to compare by identity instead of by contents. */
16: DMCompositeGetGlobalISs(dm, &isg);
18: /* Get submatrices */
19: PetscMalloc1(n * n, &submats);
20: for (i = 0, rlink = com->next; rlink; i++, rlink = rlink->next) {
21: for (j = 0, clink = com->next; clink; j++, clink = clink->next) {
22: Mat sub = NULL;
23: if (i == j) {
24: DMCreateMatrix(rlink->dm, &sub);
26: submats[i * n + j] = sub;
27: }
28: }
30: MatCreateNest(PetscObjectComm((PetscObject)dm), n, isg, n, isg, submats, J);
32: /* Disown references */
33: for (i = 0; i < n; i++) ISDestroy(&isg[i]);
34: PetscFree(isg);
36: for (i = 0; i < n * n; i++) {
37: if (submats[i]) MatDestroy(&submats[i]);
38: }
39: PetscFree(submats);
40: return 0;
41: }
43: static PetscErrorCode DMCreateMatrix_Composite_AIJ(DM dm, Mat *J)
44: {
45: DM_Composite *com = (DM_Composite *)dm->data;
46: struct DMCompositeLink *next;
47: PetscInt m, *dnz, *onz, i, j, mA;
48: Mat Atmp;
49: PetscMPIInt rank;
50: PetscBool dense = PETSC_FALSE;
52: /* use global vector to determine layout needed for matrix */
53: m = com->n;
55: MatCreate(PetscObjectComm((PetscObject)dm), J);
56: MatSetSizes(*J, m, m, PETSC_DETERMINE, PETSC_DETERMINE);
57: MatSetType(*J, dm->mattype);
59: /*
60: Extremely inefficient but will compute entire Jacobian for testing
61: */
62: PetscOptionsGetBool(((PetscObject)dm)->options, ((PetscObject)dm)->prefix, "-dmcomposite_dense_jacobian", &dense, NULL);
63: if (dense) {
64: PetscInt rstart, rend, *indices;
65: PetscScalar *values;
67: mA = com->N;
68: MatMPIAIJSetPreallocation(*J, mA, NULL, mA - m, NULL);
69: MatSeqAIJSetPreallocation(*J, mA, NULL);
71: MatGetOwnershipRange(*J, &rstart, &rend);
72: PetscMalloc2(mA, &values, mA, &indices);
73: PetscArrayzero(values, mA);
74: for (i = 0; i < mA; i++) indices[i] = i;
75: for (i = rstart; i < rend; i++) MatSetValues(*J, 1, &i, mA, indices, values, INSERT_VALUES);
76: PetscFree2(values, indices);
77: MatAssemblyBegin(*J, MAT_FINAL_ASSEMBLY);
78: MatAssemblyEnd(*J, MAT_FINAL_ASSEMBLY);
79: return 0;
80: }
82: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
83: MatPreallocateBegin(PetscObjectComm((PetscObject)dm), m, m, dnz, onz);
84: /* loop over packed objects, handling one at at time */
85: next = com->next;
86: while (next) {
87: PetscInt nc, rstart, *ccols, maxnc;
88: const PetscInt *cols, *rstarts;
89: PetscMPIInt proc;
91: DMCreateMatrix(next->dm, &Atmp);
92: MatGetOwnershipRange(Atmp, &rstart, NULL);
93: MatGetOwnershipRanges(Atmp, &rstarts);
94: MatGetLocalSize(Atmp, &mA, NULL);
96: maxnc = 0;
97: for (i = 0; i < mA; i++) {
98: MatGetRow(Atmp, rstart + i, &nc, NULL, NULL);
99: maxnc = PetscMax(nc, maxnc);
100: MatRestoreRow(Atmp, rstart + i, &nc, NULL, NULL);
101: }
102: PetscMalloc1(maxnc, &ccols);
103: for (i = 0; i < mA; i++) {
104: MatGetRow(Atmp, rstart + i, &nc, &cols, NULL);
105: /* remap the columns taking into how much they are shifted on each process */
106: for (j = 0; j < nc; j++) {
107: proc = 0;
108: while (cols[j] >= rstarts[proc + 1]) proc++;
109: ccols[j] = cols[j] + next->grstarts[proc] - rstarts[proc];
110: }
111: MatPreallocateSet(com->rstart + next->rstart + i, nc, ccols, dnz, onz);
112: MatRestoreRow(Atmp, rstart + i, &nc, &cols, NULL);
113: }
114: PetscFree(ccols);
115: MatDestroy(&Atmp);
116: next = next->next;
117: }
118: if (com->FormCoupleLocations) (*com->FormCoupleLocations)(dm, NULL, dnz, onz, __rstart, __nrows, __start, __end);
119: MatMPIAIJSetPreallocation(*J, 0, dnz, 0, onz);
120: MatSeqAIJSetPreallocation(*J, 0, dnz);
121: MatPreallocateEnd(dnz, onz);
123: if (dm->prealloc_only) return 0;
125: next = com->next;
126: while (next) {
127: PetscInt nc, rstart, row, maxnc, *ccols;
128: const PetscInt *cols, *rstarts;
129: const PetscScalar *values;
130: PetscMPIInt proc;
132: DMCreateMatrix(next->dm, &Atmp);
133: MatGetOwnershipRange(Atmp, &rstart, NULL);
134: MatGetOwnershipRanges(Atmp, &rstarts);
135: MatGetLocalSize(Atmp, &mA, NULL);
136: maxnc = 0;
137: for (i = 0; i < mA; i++) {
138: MatGetRow(Atmp, rstart + i, &nc, NULL, NULL);
139: maxnc = PetscMax(nc, maxnc);
140: MatRestoreRow(Atmp, rstart + i, &nc, NULL, NULL);
141: }
142: PetscMalloc1(maxnc, &ccols);
143: for (i = 0; i < mA; i++) {
144: MatGetRow(Atmp, rstart + i, &nc, (const PetscInt **)&cols, &values);
145: for (j = 0; j < nc; j++) {
146: proc = 0;
147: while (cols[j] >= rstarts[proc + 1]) proc++;
148: ccols[j] = cols[j] + next->grstarts[proc] - rstarts[proc];
149: }
150: row = com->rstart + next->rstart + i;
151: MatSetValues(*J, 1, &row, nc, ccols, values, INSERT_VALUES);
152: MatRestoreRow(Atmp, rstart + i, &nc, (const PetscInt **)&cols, &values);
153: }
154: PetscFree(ccols);
155: MatDestroy(&Atmp);
156: next = next->next;
157: }
158: if (com->FormCoupleLocations) {
159: PetscInt __rstart;
160: MatGetOwnershipRange(*J, &__rstart, NULL);
161: (*com->FormCoupleLocations)(dm, *J, NULL, NULL, __rstart, 0, 0, 0);
162: }
163: MatAssemblyBegin(*J, MAT_FINAL_ASSEMBLY);
164: MatAssemblyEnd(*J, MAT_FINAL_ASSEMBLY);
165: return 0;
166: }
168: PetscErrorCode DMCreateMatrix_Composite(DM dm, Mat *J)
169: {
170: PetscBool usenest;
171: ISLocalToGlobalMapping ltogmap;
173: DMSetFromOptions(dm);
174: DMSetUp(dm);
175: PetscStrcmp(dm->mattype, MATNEST, &usenest);
176: if (usenest) DMCreateMatrix_Composite_Nest(dm, J);
177: else DMCreateMatrix_Composite_AIJ(dm, J);
179: DMGetLocalToGlobalMapping(dm, <ogmap);
180: MatSetLocalToGlobalMapping(*J, ltogmap, ltogmap);
181: MatSetDM(*J, dm);
182: return 0;
183: }