Actual source code: mmsell.c

      1: /*
      2:    Support for the parallel SELL matrix vector multiply
      3: */
      4: #include <../src/mat/impls/sell/mpi/mpisell.h>
      5: #include <petsc/private/isimpl.h>

      7: /*
      8:    Takes the local part of an already assembled MPISELL matrix
      9:    and disassembles it. This is to allow new nonzeros into the matrix
     10:    that require more communication in the matrix vector multiply.
     11:    Thus certain data-structures must be rebuilt.

     13:    Kind of slow! But that's what application programmers get when
     14:    they are sloppy.
     15: */
     16: PetscErrorCode MatDisAssemble_MPISELL(Mat A)
     17: {
     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;
     22:   PetscBool    isnonzero;

     24:   PetscFunctionBegin;
     25:   /* free stuff related to matrix-vec multiply */
     26:   PetscCall(VecDestroy(&sell->lvec));
     27:   PetscCall(VecScatterDestroy(&sell->Mvctx));
     28:   if (sell->colmap) {
     29: #if defined(PETSC_USE_CTABLE)
     30:     PetscCall(PetscHMapIDestroy(&sell->colmap));
     31: #else
     32:     PetscCall(PetscFree(sell->colmap));
     33: #endif
     34:   }

     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 */
     41:   PetscCall(MatCreate(PETSC_COMM_SELF, &Bnew));
     42:   PetscCall(MatSetSizes(Bnew, B->rmap->n, N, B->rmap->n, N));
     43:   PetscCall(MatSetBlockSizesFromMats(Bnew, A, A));
     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;
     48:   }

     50:   /*
     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?
     53:    */
     54:   Bnew->nonzerostate = B->nonzerostate;

     56:   totalslices = PetscCeilInt(B->rmap->n, Bsell->sliceheight);
     57:   for (i = 0; i < totalslices; i++) { /* loop over slices */
     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));
     61:     }
     62:   }

     64:   PetscCall(PetscFree(sell->garray));
     65:   PetscCall(MatDestroy(&B));

     67:   sell->B          = Bnew;
     68:   A->was_assembled = PETSC_FALSE;
     69:   A->assembled     = PETSC_FALSE;
     70:   PetscFunctionReturn(PETSC_SUCCESS);
     71: }

     73: PetscErrorCode MatSetUpMultiply_MPISELL(Mat mat)
     74: {
     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;
     78:   IS           from, to;
     79:   Vec          gvec;
     80:   PetscBool    isnonzero;
     81: #if defined(PETSC_USE_CTABLE)
     82:   PetscHMapI    gid1_lid1 = NULL;
     83:   PetscHashIter tpos;
     84:   PetscInt      gid, lid;
     85: #else
     86:   PetscInt N = mat->cmap->N, *indices;
     87: #endif

     89:   PetscFunctionBegin;
     90:   totalslices = PetscCeilInt(sell->B->rmap->n, B->sliceheight);

     92:   /* ec counts the number of columns that contain nonzeros */
     93: #if defined(PETSC_USE_CTABLE)
     94:   /* use a table */
     95:   PetscCall(PetscHMapICreateWithSize(sell->B->rmap->n, &gid1_lid1));
     96:   for (i = 0; i < totalslices; i++) { /* loop over slices */
     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]);
     99:       if (isnonzero) { /* check the mask bit */
    100:         PetscInt data, gid1 = bcolidx[j] + 1;

    102:         PetscCall(PetscHMapIGetWithDefault(gid1_lid1, gid1, 0, &data));
    103:         /* one based table */
    104:         if (!data) PetscCall(PetscHMapISet(gid1_lid1, gid1, ++ec));
    105:       }
    106:     }
    107:   }

    109:   /* form array of columns we need */
    110:   PetscCall(PetscMalloc1(ec, &garray));
    111:   PetscHashIterBegin(gid1_lid1, tpos);
    112:   while (!PetscHashIterAtEnd(gid1_lid1, tpos)) {
    113:     PetscHashIterGetKey(gid1_lid1, tpos, gid);
    114:     PetscHashIterGetVal(gid1_lid1, tpos, lid);
    115:     PetscHashIterNext(gid1_lid1, tpos);
    116:     gid--;
    117:     lid--;
    118:     garray[lid] = gid;
    119:   }
    120:   PetscCall(PetscSortInt(ec, garray)); /* sort, and rebuild */
    121:   PetscCall(PetscHMapIClear(gid1_lid1));
    122:   for (i = 0; i < ec; i++) PetscCall(PetscHMapISet(gid1_lid1, garray[i] + 1, i + 1));

    124:   /* compact out the extra columns in B */
    125:   for (i = 0; i < totalslices; i++) { /* loop over slices */
    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]);
    128:       if (isnonzero) {
    129:         PetscInt gid1 = bcolidx[j] + 1;
    130:         PetscCall(PetscHMapIGetWithDefault(gid1_lid1, gid1, 0, &lid));
    131:         lid--;
    132:         bcolidx[j] = lid;
    133:       }
    134:     }
    135:   }
    136:   PetscCall(PetscLayoutDestroy(&sell->B->cmap));
    137:   PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)sell->B), ec, ec, 1, &sell->B->cmap));
    138:   PetscCall(PetscHMapIDestroy(&gid1_lid1));
    139: #else
    140:   /* Make an array as long as the number of columns */
    141:   PetscCall(PetscCalloc1(N, &indices));
    142:   /* mark those columns that are in sell->B */
    143:   for (i = 0; i < totalslices; i++) { /* loop over slices */
    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]);
    146:       if (isnonzero) {
    147:         if (!indices[bcolidx[j]]) ec++;
    148:         indices[bcolidx[j]] = 1;
    149:       }
    150:     }
    151:   }

    153:   /* form array of columns we need */
    154:   PetscCall(PetscMalloc1(ec, &garray));
    155:   ec = 0;
    156:   for (i = 0; i < N; i++) {
    157:     if (indices[i]) garray[ec++] = i;
    158:   }

    160:   /* make indices now point into garray */
    161:   for (i = 0; i < ec; i++) indices[garray[i]] = i;

    163:   /* compact out the extra columns in B */
    164:   for (i = 0; i < totalslices; i++) { /* loop over slices */
    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]);
    167:       if (isnonzero) bcolidx[j] = indices[bcolidx[j]];
    168:     }
    169:   }
    170:   PetscCall(PetscLayoutDestroy(&sell->B->cmap));
    171:   PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)sell->B), ec, ec, 1, &sell->B->cmap));
    172:   PetscCall(PetscFree(indices));
    173: #endif
    174:   /* create local vector that is used to scatter into */
    175:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, ec, &sell->lvec));
    176:   /* create two temporary Index sets for build scatter gather */
    177:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ec, garray, PETSC_COPY_VALUES, &from));
    178:   PetscCall(ISCreateStride(PETSC_COMM_SELF, ec, 0, 1, &to));

    180:   /* create temporary global vector to generate scatter context */
    181:   /* This does not allocate the array's memory so is efficient */
    182:   PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat), 1, mat->cmap->n, mat->cmap->N, NULL, &gvec));

    184:   /* generate the scatter context */
    185:   PetscCall(VecScatterCreate(gvec, from, sell->lvec, to, &sell->Mvctx));
    186:   PetscCall(VecScatterViewFromOptions(sell->Mvctx, (PetscObject)mat, "-matmult_vecscatter_view"));

    188:   sell->garray = garray;

    190:   PetscCall(ISDestroy(&from));
    191:   PetscCall(ISDestroy(&to));
    192:   PetscCall(VecDestroy(&gvec));
    193:   PetscFunctionReturn(PETSC_SUCCESS);
    194: }

    196: /*      ugly stuff added for Glenn someday we should fix this up */
    197: static PetscInt *auglyrmapd = NULL, *auglyrmapo = NULL; /* mapping from the local ordering to the "diagonal" and "off-diagonal" parts of the local matrix */
    198: static Vec       auglydd = NULL, auglyoo = NULL;        /* work vectors used to scale the two parts of the local matrix */

    200: static PetscErrorCode MatMPISELLDiagonalScaleLocalSetUp(Mat inA, Vec scale)
    201: {
    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;
    204:   PetscInt    *r_rmapd, *r_rmapo;

    206:   PetscFunctionBegin;
    207:   PetscCall(MatGetOwnershipRange(inA, &cstart, &cend));
    208:   PetscCall(MatGetSize(ina->A, NULL, &n));
    209:   PetscCall(PetscCalloc1(inA->rmap->mapping->n, &r_rmapd));
    210:   nt = 0;
    211:   for (i = 0; i < inA->rmap->mapping->n; i++) {
    212:     if (inA->rmap->mapping->indices[i] >= cstart && inA->rmap->mapping->indices[i] < cend) {
    213:       nt++;
    214:       r_rmapd[i] = inA->rmap->mapping->indices[i] + 1;
    215:     }
    216:   }
    217:   PetscCheck(nt == n, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Hmm nt %" PetscInt_FMT " n %" PetscInt_FMT, nt, n);
    218:   PetscCall(PetscMalloc1(n, &auglyrmapd));
    219:   for (i = 0; i < inA->rmap->mapping->n; i++) {
    220:     if (r_rmapd[i]) auglyrmapd[(r_rmapd[i] - 1) - cstart] = i;
    221:   }
    222:   PetscCall(PetscFree(r_rmapd));
    223:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &auglydd));
    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));
    228:   nt = 0;
    229:   for (i = 0; i < inA->rmap->mapping->n; i++) {
    230:     if (lindices[inA->rmap->mapping->indices[i]]) {
    231:       nt++;
    232:       r_rmapo[i] = lindices[inA->rmap->mapping->indices[i]];
    233:     }
    234:   }
    235:   PetscCheck(nt <= no, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Hmm nt %" PetscInt_FMT " no %" PetscInt_FMT, nt, n);
    236:   PetscCall(PetscFree(lindices));
    237:   PetscCall(PetscMalloc1(nt, &auglyrmapo));
    238:   for (i = 0; i < inA->rmap->mapping->n; i++) {
    239:     if (r_rmapo[i]) auglyrmapo[(r_rmapo[i] - 1)] = i;
    240:   }
    241:   PetscCall(PetscFree(r_rmapo));
    242:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, nt, &auglyoo));
    243:   PetscFunctionReturn(PETSC_SUCCESS);
    244: }

    246: PetscErrorCode MatDiagonalScaleLocal_MPISELL(Mat A, Vec scale)
    247: {
    248:   Mat_MPISELL       *a = (Mat_MPISELL *)A->data; /*access private part of matrix */
    249:   PetscInt           n, i;
    250:   PetscScalar       *d, *o;
    251:   const PetscScalar *s;

    253:   PetscFunctionBegin;
    254:   if (!auglyrmapd) PetscCall(MatMPISELLDiagonalScaleLocalSetUp(A, scale));
    255:   PetscCall(VecGetArrayRead(scale, &s));
    256:   PetscCall(VecGetLocalSize(auglydd, &n));
    257:   PetscCall(VecGetArray(auglydd, &d));
    258:   for (i = 0; i < n; i++) d[i] = s[auglyrmapd[i]]; /* copy "diagonal" (true local) portion of scale into dd vector */
    259:   PetscCall(VecRestoreArray(auglydd, &d));
    260:   /* column scale "diagonal" portion of local matrix */
    261:   PetscCall(MatDiagonalScale(a->A, NULL, auglydd));
    262:   PetscCall(VecGetLocalSize(auglyoo, &n));
    263:   PetscCall(VecGetArray(auglyoo, &o));
    264:   for (i = 0; i < n; i++) o[i] = s[auglyrmapo[i]]; /* copy "off-diagonal" portion of scale into oo vector */
    265:   PetscCall(VecRestoreArrayRead(scale, &s));
    266:   PetscCall(VecRestoreArray(auglyoo, &o));
    267:   /* column scale "off-diagonal" portion of local matrix */
    268:   PetscCall(MatDiagonalScale(a->B, NULL, auglyoo));
    269:   PetscFunctionReturn(PETSC_SUCCESS);
    270: }