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, &ltogmap);
180:   MatSetLocalToGlobalMapping(*J, ltogmap, ltogmap);
181:   MatSetDM(*J, dm);
182:   return 0;
183: }