sell/mpi/mmsell.c

11    Thus certain data-structures must be rebuilt.
18   Mat_MPISELL *sell  = (Mat_MPISELL *)A->data;
19   Mat          B     = sell->B, Bnew;
20   Mat_SeqSELL *Bsell = (Mat_SeqSELL *)B->data;
21   PetscInt     i, j, totalslices, N = A->cmap->N, row;
25   /* free stuff related to matrix-vec multiply */
26   PetscCall(VecDestroy(&sell->lvec));
27   PetscCall(VecScatterDestroy(&sell->Mvctx));
28   if (sell->colmap) {
30     PetscCall(PetscHMapIDestroy(&sell->colmap));
32     PetscCall(PetscFree(sell->colmap));
36   /* make sure that B is assembled so we can access its values */
37   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
38   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
40   /* invent new B and copy stuff over */
42   PetscCall(MatSetSizes(Bnew, B->rmap->n, N, B->rmap->n, N));
44   PetscCall(MatSetType(Bnew, ((PetscObject)B)->type_name));
45   PetscCall(MatSeqSELLSetPreallocation(Bnew, 0, Bsell->rlen));
46   if (Bsell->nonew >= 0) { /* Inherit insertion error options (if positive). */
47     ((Mat_SeqSELL *)Bnew->data)->nonew = Bsell->nonew;
51    Ensure that B's nonzerostate is monotonically increasing.
52    Or should this follow the MatSetValues() loop to preserve B's nonzerstate across a MatDisAssemble() call?
54   Bnew->nonzerostate = B->nonzerostate;
56   totalslices = PetscCeilInt(B->rmap->n, Bsell->sliceheight);
58     for (j = Bsell->sliidx[i], row = 0; j < Bsell->sliidx[i + 1]; j++, row = (row + 1) % Bsell->sliceheight) {
59       isnonzero = (PetscBool)((j - Bsell->sliidx[i]) / Bsell->sliceheight < Bsell->rlen[Bsell->sliceheight * i + row]);
60       if (isnonzero) PetscCall(MatSetValue(Bnew, Bsell->sliceheight * i + row, sell->garray[Bsell->colidx[j]], Bsell->val[j], B->insertmode));
64   PetscCall(PetscFree(sell->garray));
65   PetscCall(MatDestroy(&B));
67   sell->B          = Bnew;
68   A->was_assembled = PETSC_FALSE;
69   A->assembled     = PETSC_FALSE;
75   Mat_MPISELL *sell = (Mat_MPISELL *)mat->data;
76   Mat_SeqSELL *B    = (Mat_SeqSELL *)sell->B->data;
77   PetscInt     i, j, *bcolidx = B->colidx, ec = 0, *garray, totalslices;
86   PetscInt N = mat->cmap->N, *indices;
90   totalslices = PetscCeilInt(sell->B->rmap->n, B->sliceheight);
95   PetscCall(PetscHMapICreateWithSize(sell->B->rmap->n, &gid1_lid1));
97     for (j = B->sliidx[i]; j < B->sliidx[i + 1]; j++) {
98       isnonzero = (PetscBool)((j - B->sliidx[i]) / B->sliceheight < B->rlen[i * B->sliceheight + j % B->sliceheight]);
116     gid--;
117     lid--;
124   /* compact out the extra columns in B */
126     for (j = B->sliidx[i]; j < B->sliidx[i + 1]; j++) {
127       isnonzero = (PetscBool)((j - B->sliidx[i]) / B->sliceheight < B->rlen[i * B->sliceheight + j % B->sliceheight]);
131         lid--;
136   PetscCall(PetscLayoutDestroy(&sell->B->cmap));
137   PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)sell->B), ec, ec, 1, &sell->B->cmap));
142   /* mark those columns that are in sell->B */
144     for (j = B->sliidx[i]; j < B->sliidx[i + 1]; j++) {
145       isnonzero = (PetscBool)((j - B->sliidx[i]) / B->sliceheight < B->rlen[i * B->sliceheight + j % B->sliceheight]);
163   /* compact out the extra columns in B */
165     for (j = B->sliidx[i]; j < B->sliidx[i + 1]; j++) {
166       isnonzero = (PetscBool)((j - B->sliidx[i]) / B->sliceheight < B->rlen[i * B->sliceheight + j % B->sliceheight]);
170   PetscCall(PetscLayoutDestroy(&sell->B->cmap));
171   PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)sell->B), ec, ec, 1, &sell->B->cmap));
175   PetscCall(VecCreateSeq(PETSC_COMM_SELF, ec, &sell->lvec));
182   PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat), 1, mat->cmap->n, mat->cmap->N, NULL, &gvec));
185   PetscCall(VecScatterCreate(gvec, from, sell->lvec, to, &sell->Mvctx));
186   PetscCall(VecScatterViewFromOptions(sell->Mvctx, (PetscObject)mat, "-matmult_vecscatter_view"));
188   sell->garray = garray;
197 static PetscInt *auglyrmapd = NULL, *auglyrmapo = NULL; /* mapping from the local ordering to the "diagonal" and "off-diagonal" parts of the local matrix */
202   Mat_MPISELL *ina = (Mat_MPISELL *)inA->data; /*access private part of matrix */
203   PetscInt     i, n, nt, cstart, cend, no, *garray = ina->garray, *lindices;
208   PetscCall(MatGetSize(ina->A, NULL, &n));
209   PetscCall(PetscCalloc1(inA->rmap->mapping->n, &r_rmapd));
211   for (i = 0; i < inA->rmap->mapping->n; i++) {
212     if (inA->rmap->mapping->indices[i] >= cstart && inA->rmap->mapping->indices[i] < cend) {
214       r_rmapd[i] = inA->rmap->mapping->indices[i] + 1;
219   for (i = 0; i < inA->rmap->mapping->n; i++) {
220     if (r_rmapd[i]) auglyrmapd[(r_rmapd[i] - 1) - cstart] = i;
224   PetscCall(PetscCalloc1(inA->cmap->N, &lindices));
225   for (i = 0; i < ina->B->cmap->n; i++) lindices[garray[i]] = i + 1;
226   no = inA->rmap->mapping->n - nt;
227   PetscCall(PetscCalloc1(inA->rmap->mapping->n, &r_rmapo));
229   for (i = 0; i < inA->rmap->mapping->n; i++) {
230     if (lindices[inA->rmap->mapping->indices[i]]) {
232       r_rmapo[i] = lindices[inA->rmap->mapping->indices[i]];
238   for (i = 0; i < inA->rmap->mapping->n; i++) {
239     if (r_rmapo[i]) auglyrmapo[(r_rmapo[i] - 1)] = i;
248   Mat_MPISELL       *a = (Mat_MPISELL *)A->data; /*access private part of matrix */
261   PetscCall(MatDiagonalScale(a->A, NULL, auglydd));
264   for (i = 0; i < n; i++) o[i] = s[auglyrmapo[i]]; /* copy "off-diagonal" portion of scale into oo vector */
267   /* column scale "off-diagonal" portion of local matrix */
268   PetscCall(MatDiagonalScale(a->B, NULL, auglyoo));