Actual source code: mpidense.c

  1: /*
  2:    Basic functions for basic parallel dense matrices.
  3:    Portions of this code are under:
  4:    Copyright (c) 2022 Advanced Micro Devices, Inc. All rights reserved.
  5: */

  7: #include <../src/mat/impls/dense/mpi/mpidense.h>
  8: #include <../src/mat/impls/aij/mpi/mpiaij.h>
  9: #include <petscblaslapack.h>

 11: /*@
 12:   MatDenseGetLocalMatrix - For a `MATMPIDENSE` or `MATSEQDENSE` matrix returns the sequential
 13:   matrix that represents the operator. For sequential matrices it returns itself.

 15:   Input Parameter:
 16: . A - the sequential or MPI `MATDENSE` matrix

 18:   Output Parameter:
 19: . B - the inner matrix

 21:   Level: intermediate

 23: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATMPIDENSE`, `MATSEQDENSE`
 24: @*/
 25: PetscErrorCode MatDenseGetLocalMatrix(Mat A, Mat *B)
 26: {
 27:   Mat_MPIDense *mat = (Mat_MPIDense *)A->data;
 28:   PetscBool     flg;

 30:   PetscFunctionBegin;
 32:   PetscAssertPointer(B, 2);
 33:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &flg));
 34:   if (flg) *B = mat->A;
 35:   else {
 36:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQDENSE, &flg));
 37:     PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for matrix type %s", ((PetscObject)A)->type_name);
 38:     *B = A;
 39:   }
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

 43: static PetscErrorCode MatCopy_MPIDense(Mat A, Mat B, MatStructure s)
 44: {
 45:   Mat_MPIDense *Amat = (Mat_MPIDense *)A->data;
 46:   Mat_MPIDense *Bmat = (Mat_MPIDense *)B->data;

 48:   PetscFunctionBegin;
 49:   PetscCall(MatCopy(Amat->A, Bmat->A, s));
 50:   PetscFunctionReturn(PETSC_SUCCESS);
 51: }

 53: PetscErrorCode MatShift_MPIDense(Mat A, PetscScalar alpha)
 54: {
 55:   Mat_MPIDense *mat = (Mat_MPIDense *)A->data;
 56:   PetscInt      j, lda, rstart = A->rmap->rstart, rend = A->rmap->rend, rend2;
 57:   PetscScalar  *v;

 59:   PetscFunctionBegin;
 60:   PetscCall(MatDenseGetArray(mat->A, &v));
 61:   PetscCall(MatDenseGetLDA(mat->A, &lda));
 62:   rend2 = PetscMin(rend, A->cmap->N);
 63:   if (rend2 > rstart) {
 64:     for (j = rstart; j < rend2; j++) v[j - rstart + j * lda] += alpha;
 65:     PetscCall(PetscLogFlops(rend2 - rstart));
 66:   }
 67:   PetscCall(MatDenseRestoreArray(mat->A, &v));
 68:   PetscFunctionReturn(PETSC_SUCCESS);
 69: }

 71: static PetscErrorCode MatGetRow_MPIDense(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
 72: {
 73:   Mat_MPIDense *mat = (Mat_MPIDense *)A->data;
 74:   PetscInt      lrow, rstart = A->rmap->rstart, rend = A->rmap->rend;

 76:   PetscFunctionBegin;
 77:   PetscCheck(row >= rstart && row < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "only local rows");
 78:   lrow = row - rstart;
 79:   PetscCall(MatGetRow(mat->A, lrow, nz, (const PetscInt **)idx, (const PetscScalar **)v));
 80:   PetscFunctionReturn(PETSC_SUCCESS);
 81: }

 83: static PetscErrorCode MatRestoreRow_MPIDense(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
 84: {
 85:   Mat_MPIDense *mat = (Mat_MPIDense *)A->data;
 86:   PetscInt      lrow, rstart = A->rmap->rstart, rend = A->rmap->rend;

 88:   PetscFunctionBegin;
 89:   PetscCheck(row >= rstart && row < rend, PETSC_COMM_SELF, PETSC_ERR_SUP, "only local rows");
 90:   lrow = row - rstart;
 91:   PetscCall(MatRestoreRow(mat->A, lrow, nz, (const PetscInt **)idx, (const PetscScalar **)v));
 92:   PetscFunctionReturn(PETSC_SUCCESS);
 93: }

 95: static PetscErrorCode MatGetDiagonalBlock_MPIDense(Mat A, Mat *a)
 96: {
 97:   Mat_MPIDense *mdn = (Mat_MPIDense *)A->data;
 98:   PetscInt      m = A->rmap->n, rstart = A->rmap->rstart;
 99:   PetscScalar  *array;
100:   MPI_Comm      comm;
101:   PetscBool     flg;
102:   Mat           B;

104:   PetscFunctionBegin;
105:   PetscCall(MatHasCongruentLayouts(A, &flg));
106:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only square matrices supported.");
107:   PetscCall(PetscObjectQuery((PetscObject)A, "DiagonalBlock", (PetscObject *)&B));
108:   if (!B) { /* This should use MatDenseGetSubMatrix (not create), but we would need a call like MatRestoreDiagonalBlock */
109: #if PetscDefined(HAVE_CUDA)
110:     PetscCall(PetscObjectTypeCompare((PetscObject)mdn->A, MATSEQDENSECUDA, &flg));
111:     PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not coded for %s. Send an email to petsc-dev@mcs.anl.gov to request this feature", MATSEQDENSECUDA);
112: #elif PetscDefined(HAVE_HIP)
113:     PetscCall(PetscObjectTypeCompare((PetscObject)mdn->A, MATSEQDENSEHIP, &flg));
114:     PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not coded for %s. Send an email to petsc-dev@mcs.anl.gov to request this feature", MATSEQDENSEHIP);
115: #endif
116:     PetscCall(PetscObjectGetComm((PetscObject)mdn->A, &comm));
117:     PetscCall(MatCreate(comm, &B));
118:     PetscCall(MatSetSizes(B, m, m, m, m));
119:     PetscCall(MatSetType(B, ((PetscObject)mdn->A)->type_name));
120:     PetscCall(MatDenseGetArrayRead(mdn->A, (const PetscScalar **)&array));
121:     PetscCall(MatSeqDenseSetPreallocation(B, array + m * rstart));
122:     PetscCall(MatDenseRestoreArrayRead(mdn->A, (const PetscScalar **)&array));
123:     PetscCall(PetscObjectCompose((PetscObject)A, "DiagonalBlock", (PetscObject)B));
124:     *a = B;
125:     PetscCall(MatDestroy(&B));
126:   } else *a = B;
127:   PetscFunctionReturn(PETSC_SUCCESS);
128: }

130: static PetscErrorCode MatSetValues_MPIDense(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], const PetscScalar v[], InsertMode addv)
131: {
132:   Mat_MPIDense *A = (Mat_MPIDense *)mat->data;
133:   PetscInt      i, j, rstart = mat->rmap->rstart, rend = mat->rmap->rend, row;
134:   PetscBool     roworiented = A->roworiented;

136:   PetscFunctionBegin;
137:   for (i = 0; i < m; i++) {
138:     if (idxm[i] < 0) continue;
139:     PetscCheck(idxm[i] < mat->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large");
140:     if (idxm[i] >= rstart && idxm[i] < rend) {
141:       row = idxm[i] - rstart;
142:       if (roworiented) {
143:         PetscCall(MatSetValues(A->A, 1, &row, n, idxn, v ? v + i * n : NULL, addv));
144:       } else {
145:         for (j = 0; j < n; j++) {
146:           if (idxn[j] < 0) continue;
147:           PetscCheck(idxn[j] < mat->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large");
148:           PetscCall(MatSetValues(A->A, 1, &row, 1, &idxn[j], v ? v + i + j * m : NULL, addv));
149:         }
150:       }
151:     } else if (!A->donotstash) {
152:       mat->assembled = PETSC_FALSE;
153:       if (roworiented) {
154:         PetscCall(MatStashValuesRow_Private(&mat->stash, idxm[i], n, idxn, v ? v + i * n : NULL, PETSC_FALSE));
155:       } else {
156:         PetscCall(MatStashValuesCol_Private(&mat->stash, idxm[i], n, idxn, v ? v + i : NULL, m, PETSC_FALSE));
157:       }
158:     }
159:   }
160:   PetscFunctionReturn(PETSC_SUCCESS);
161: }

163: static PetscErrorCode MatGetValues_MPIDense(Mat mat, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], PetscScalar v[])
164: {
165:   Mat_MPIDense *mdn = (Mat_MPIDense *)mat->data;
166:   PetscInt      i, j, rstart = mat->rmap->rstart, rend = mat->rmap->rend, row;

168:   PetscFunctionBegin;
169:   for (i = 0; i < m; i++) {
170:     if (idxm[i] < 0) continue; /* negative row */
171:     PetscCheck(idxm[i] < mat->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large");
172:     if (idxm[i] >= rstart && idxm[i] < rend) {
173:       row = idxm[i] - rstart;
174:       for (j = 0; j < n; j++) {
175:         if (idxn[j] < 0) continue; /* negative column */
176:         PetscCheck(idxn[j] < mat->cmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large");
177:         PetscCall(MatGetValues(mdn->A, 1, &row, 1, &idxn[j], v + i * n + j));
178:       }
179:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only local values currently supported");
180:   }
181:   PetscFunctionReturn(PETSC_SUCCESS);
182: }

184: static PetscErrorCode MatDenseGetLDA_MPIDense(Mat A, PetscInt *lda)
185: {
186:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

188:   PetscFunctionBegin;
189:   PetscCall(MatDenseGetLDA(a->A, lda));
190:   PetscFunctionReturn(PETSC_SUCCESS);
191: }

193: static PetscErrorCode MatDenseSetLDA_MPIDense(Mat A, PetscInt lda)
194: {
195:   Mat_MPIDense *a     = (Mat_MPIDense *)A->data;
196:   MatType       mtype = MATSEQDENSE;

198:   PetscFunctionBegin;
199:   if (!a->A) {
200:     PetscCheck(!a->matinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
201:     PetscCall(PetscLayoutSetUp(A->rmap));
202:     PetscCall(PetscLayoutSetUp(A->cmap));
203:     PetscCall(MatCreate(PETSC_COMM_SELF, &a->A));
204:     PetscCall(MatSetSizes(a->A, A->rmap->n, A->cmap->N, A->rmap->n, A->cmap->N));
205: #if PetscDefined(HAVE_CUDA)
206:     PetscBool iscuda;
207:     PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIDENSECUDA, &iscuda));
208:     if (iscuda) mtype = MATSEQDENSECUDA;
209: #elif PetscDefined(HAVE_HIP)
210:     PetscBool iship;
211:     PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIDENSEHIP, &iship));
212:     if (iship) mtype = MATSEQDENSEHIP;
213: #endif
214:     PetscCall(MatSetType(a->A, mtype));
215:   }
216:   PetscCall(MatDenseSetLDA(a->A, lda));
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: static PetscErrorCode MatDenseGetArray_MPIDense(Mat A, PetscScalar **array)
221: {
222:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

224:   PetscFunctionBegin;
225:   PetscCheck(!a->matinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
226:   PetscCall(MatDenseGetArray(a->A, array));
227:   PetscFunctionReturn(PETSC_SUCCESS);
228: }

230: static PetscErrorCode MatDenseGetArrayRead_MPIDense(Mat A, PetscScalar **array)
231: {
232:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

234:   PetscFunctionBegin;
235:   PetscCheck(!a->matinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
236:   PetscCall(MatDenseGetArrayRead(a->A, (const PetscScalar **)array));
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }

240: static PetscErrorCode MatDenseGetArrayWrite_MPIDense(Mat A, PetscScalar **array)
241: {
242:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

244:   PetscFunctionBegin;
245:   PetscCheck(!a->matinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
246:   PetscCall(MatDenseGetArrayWrite(a->A, array));
247:   PetscFunctionReturn(PETSC_SUCCESS);
248: }

250: static PetscErrorCode MatDensePlaceArray_MPIDense(Mat A, const PetscScalar *array)
251: {
252:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

254:   PetscFunctionBegin;
255:   PetscCheck(!a->vecinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
256:   PetscCheck(!a->matinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
257:   PetscCall(MatDensePlaceArray(a->A, array));
258:   PetscFunctionReturn(PETSC_SUCCESS);
259: }

261: static PetscErrorCode MatDenseResetArray_MPIDense(Mat A)
262: {
263:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

265:   PetscFunctionBegin;
266:   PetscCheck(!a->vecinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
267:   PetscCheck(!a->matinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
268:   PetscCall(MatDenseResetArray(a->A));
269:   PetscFunctionReturn(PETSC_SUCCESS);
270: }

272: static PetscErrorCode MatDenseReplaceArray_MPIDense(Mat A, const PetscScalar *array)
273: {
274:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

276:   PetscFunctionBegin;
277:   PetscCheck(!a->vecinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
278:   PetscCheck(!a->matinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
279:   PetscCall(MatDenseReplaceArray(a->A, array));
280:   PetscFunctionReturn(PETSC_SUCCESS);
281: }

283: static PetscErrorCode MatCreateSubMatrix_MPIDense(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B)
284: {
285:   Mat_MPIDense      *mat = (Mat_MPIDense *)A->data, *newmatd;
286:   PetscInt           lda, i, j, rstart, rend, nrows, ncols, Ncols, nlrows, nlcols;
287:   const PetscInt    *irow, *icol;
288:   const PetscScalar *v;
289:   PetscScalar       *bv;
290:   Mat                newmat;
291:   IS                 iscol_local;
292:   MPI_Comm           comm_is, comm_mat;

294:   PetscFunctionBegin;
295:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm_mat));
296:   PetscCall(PetscObjectGetComm((PetscObject)iscol, &comm_is));
297:   PetscCheck(comm_mat == comm_is, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMECOMM, "IS communicator must match matrix communicator");

299:   PetscCall(ISAllGather(iscol, &iscol_local));
300:   PetscCall(ISGetIndices(isrow, &irow));
301:   PetscCall(ISGetIndices(iscol_local, &icol));
302:   PetscCall(ISGetLocalSize(isrow, &nrows));
303:   PetscCall(ISGetLocalSize(iscol, &ncols));
304:   PetscCall(ISGetSize(iscol, &Ncols)); /* global number of columns, size of iscol_local */

306:   /* No parallel redistribution currently supported! Should really check each index set
307:      to confirm that it is OK.  ... Currently supports only submatrix same partitioning as
308:      original matrix! */

310:   PetscCall(MatGetLocalSize(A, &nlrows, &nlcols));
311:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));

313:   /* Check submatrix call */
314:   if (scall == MAT_REUSE_MATRIX) {
315:     /* SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); */
316:     /* Really need to test rows and column sizes! */
317:     newmat = *B;
318:   } else {
319:     /* Create and fill new matrix */
320:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &newmat));
321:     PetscCall(MatSetSizes(newmat, nrows, ncols, PETSC_DECIDE, Ncols));
322:     PetscCall(MatSetType(newmat, ((PetscObject)A)->type_name));
323:     PetscCall(MatMPIDenseSetPreallocation(newmat, NULL));
324:   }

326:   /* Now extract the data pointers and do the copy, column at a time */
327:   newmatd = (Mat_MPIDense *)newmat->data;
328:   PetscCall(MatDenseGetArray(newmatd->A, &bv));
329:   PetscCall(MatDenseGetArrayRead(mat->A, &v));
330:   PetscCall(MatDenseGetLDA(mat->A, &lda));
331:   for (i = 0; i < Ncols; i++) {
332:     const PetscScalar *av = v + lda * icol[i];
333:     for (j = 0; j < nrows; j++) *bv++ = av[irow[j] - rstart];
334:   }
335:   PetscCall(MatDenseRestoreArrayRead(mat->A, &v));
336:   PetscCall(MatDenseRestoreArray(newmatd->A, &bv));

338:   /* Assemble the matrices so that the correct flags are set */
339:   PetscCall(MatAssemblyBegin(newmat, MAT_FINAL_ASSEMBLY));
340:   PetscCall(MatAssemblyEnd(newmat, MAT_FINAL_ASSEMBLY));

342:   /* Free work space */
343:   PetscCall(ISRestoreIndices(isrow, &irow));
344:   PetscCall(ISRestoreIndices(iscol_local, &icol));
345:   PetscCall(ISDestroy(&iscol_local));
346:   *B = newmat;
347:   PetscFunctionReturn(PETSC_SUCCESS);
348: }

350: static PetscErrorCode MatDenseRestoreArray_MPIDense(Mat A, PetscScalar **array)
351: {
352:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

354:   PetscFunctionBegin;
355:   PetscCall(MatDenseRestoreArray(a->A, array));
356:   PetscFunctionReturn(PETSC_SUCCESS);
357: }

359: static PetscErrorCode MatDenseRestoreArrayRead_MPIDense(Mat A, PetscScalar **array)
360: {
361:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

363:   PetscFunctionBegin;
364:   PetscCall(MatDenseRestoreArrayRead(a->A, (const PetscScalar **)array));
365:   PetscFunctionReturn(PETSC_SUCCESS);
366: }

368: static PetscErrorCode MatDenseRestoreArrayWrite_MPIDense(Mat A, PetscScalar **array)
369: {
370:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

372:   PetscFunctionBegin;
373:   PetscCall(MatDenseRestoreArrayWrite(a->A, array));
374:   PetscFunctionReturn(PETSC_SUCCESS);
375: }

377: static PetscErrorCode MatAssemblyBegin_MPIDense(Mat mat, MatAssemblyType mode)
378: {
379:   Mat_MPIDense *mdn = (Mat_MPIDense *)mat->data;
380:   PetscInt      nstash, reallocs;

382:   PetscFunctionBegin;
383:   if (mdn->donotstash || mat->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS);

385:   PetscCall(MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range));
386:   PetscCall(MatStashGetInfo_Private(&mat->stash, &nstash, &reallocs));
387:   PetscCall(PetscInfo(mdn->A, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
388:   PetscFunctionReturn(PETSC_SUCCESS);
389: }

391: static PetscErrorCode MatAssemblyEnd_MPIDense(Mat mat, MatAssemblyType mode)
392: {
393:   Mat_MPIDense *mdn = (Mat_MPIDense *)mat->data;
394:   PetscInt      i, *row, *col, flg, j, rstart, ncols;
395:   PetscMPIInt   n;
396:   PetscScalar  *val;

398:   PetscFunctionBegin;
399:   if (!mdn->donotstash && !mat->nooffprocentries) {
400:     /*  wait on receives */
401:     while (1) {
402:       PetscCall(MatStashScatterGetMesg_Private(&mat->stash, &n, &row, &col, &val, &flg));
403:       if (!flg) break;

405:       for (i = 0; i < n;) {
406:         /* Now identify the consecutive vals belonging to the same row */
407:         for (j = i, rstart = row[j]; j < n; j++) {
408:           if (row[j] != rstart) break;
409:         }
410:         if (j < n) ncols = j - i;
411:         else ncols = n - i;
412:         /* Now assemble all these values with a single function call */
413:         PetscCall(MatSetValues_MPIDense(mat, 1, row + i, ncols, col + i, val + i, mat->insertmode));
414:         i = j;
415:       }
416:     }
417:     PetscCall(MatStashScatterEnd_Private(&mat->stash));
418:   }

420:   PetscCall(MatAssemblyBegin(mdn->A, mode));
421:   PetscCall(MatAssemblyEnd(mdn->A, mode));
422:   PetscFunctionReturn(PETSC_SUCCESS);
423: }

425: static PetscErrorCode MatZeroEntries_MPIDense(Mat A)
426: {
427:   Mat_MPIDense *l = (Mat_MPIDense *)A->data;

429:   PetscFunctionBegin;
430:   PetscCall(MatZeroEntries(l->A));
431:   PetscFunctionReturn(PETSC_SUCCESS);
432: }

434: static PetscErrorCode MatZeroRows_MPIDense(Mat A, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
435: {
436:   Mat_MPIDense *l = (Mat_MPIDense *)A->data;
437:   PetscInt      i, len, *lrows;

439:   PetscFunctionBegin;
440:   /* get locally owned rows */
441:   PetscCall(PetscLayoutMapLocal(A->rmap, n, rows, &len, &lrows, NULL));
442:   /* fix right-hand side if needed */
443:   if (x && b) {
444:     const PetscScalar *xx;
445:     PetscScalar       *bb;

447:     PetscCall(VecGetArrayRead(x, &xx));
448:     PetscCall(VecGetArrayWrite(b, &bb));
449:     for (i = 0; i < len; ++i) bb[lrows[i]] = diag * xx[lrows[i]];
450:     PetscCall(VecRestoreArrayRead(x, &xx));
451:     PetscCall(VecRestoreArrayWrite(b, &bb));
452:   }
453:   PetscCall(MatZeroRows(l->A, len, lrows, 0.0, NULL, NULL));
454:   if (diag != 0.0) {
455:     Vec d;

457:     PetscCall(MatCreateVecs(A, NULL, &d));
458:     PetscCall(VecSet(d, diag));
459:     PetscCall(MatDiagonalSet(A, d, INSERT_VALUES));
460:     PetscCall(VecDestroy(&d));
461:   }
462:   PetscCall(PetscFree(lrows));
463:   PetscFunctionReturn(PETSC_SUCCESS);
464: }

466: PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat, Vec, Vec);
467: PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat, Vec, Vec, Vec);
468: PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat, Vec, Vec);
469: PETSC_INTERN PetscErrorCode MatMultTransposeAdd_SeqDense(Mat, Vec, Vec, Vec);

471: static PetscErrorCode MatMult_MPIDense(Mat mat, Vec xx, Vec yy)
472: {
473:   Mat_MPIDense      *mdn = (Mat_MPIDense *)mat->data;
474:   const PetscScalar *ax;
475:   PetscScalar       *ay;
476:   PetscMemType       axmtype, aymtype;

478:   PetscFunctionBegin;
479:   if (!mdn->Mvctx) PetscCall(MatSetUpMultiply_MPIDense(mat));
480:   PetscCall(VecGetArrayReadAndMemType(xx, &ax, &axmtype));
481:   PetscCall(VecGetArrayAndMemType(mdn->lvec, &ay, &aymtype));
482:   PetscCall(PetscSFBcastWithMemTypeBegin(mdn->Mvctx, MPIU_SCALAR, axmtype, ax, aymtype, ay, MPI_REPLACE));
483:   PetscCall(PetscSFBcastEnd(mdn->Mvctx, MPIU_SCALAR, ax, ay, MPI_REPLACE));
484:   PetscCall(VecRestoreArrayAndMemType(mdn->lvec, &ay));
485:   PetscCall(VecRestoreArrayReadAndMemType(xx, &ax));
486:   PetscCall((*mdn->A->ops->mult)(mdn->A, mdn->lvec, yy));
487:   PetscFunctionReturn(PETSC_SUCCESS);
488: }

490: static PetscErrorCode MatMultAddColumnRange_MPIDense(Mat mat, Vec xx, Vec yy, Vec zz, PetscInt c_start, PetscInt c_end)
491: {
492:   Mat_MPIDense      *mdn = (Mat_MPIDense *)mat->data;
493:   const PetscScalar *ax;
494:   PetscScalar       *ay;
495:   PetscMemType       axmtype, aymtype;

497:   PetscFunctionBegin;
498:   if (!mdn->Mvctx) PetscCall(MatSetUpMultiply_MPIDense(mat));
499:   PetscCall(VecGetArrayReadAndMemType(xx, &ax, &axmtype));
500:   PetscCall(VecGetArrayAndMemType(mdn->lvec, &ay, &aymtype));
501:   PetscCall(PetscSFBcastWithMemTypeBegin(mdn->Mvctx, MPIU_SCALAR, axmtype, ax, aymtype, ay, MPI_REPLACE));
502:   PetscCall(PetscSFBcastEnd(mdn->Mvctx, MPIU_SCALAR, ax, ay, MPI_REPLACE));
503:   PetscCall(VecRestoreArrayAndMemType(mdn->lvec, &ay));
504:   PetscCall(VecRestoreArrayReadAndMemType(xx, &ax));
505:   PetscUseMethod(mdn->A, "MatMultAddColumnRange_C", (Mat, Vec, Vec, Vec, PetscInt, PetscInt), (mdn->A, mdn->lvec, yy, zz, c_start, c_end));
506:   PetscFunctionReturn(PETSC_SUCCESS);
507: }

509: static PetscErrorCode MatMultAdd_MPIDense(Mat mat, Vec xx, Vec yy, Vec zz)
510: {
511:   Mat_MPIDense      *mdn = (Mat_MPIDense *)mat->data;
512:   const PetscScalar *ax;
513:   PetscScalar       *ay;
514:   PetscMemType       axmtype, aymtype;

516:   PetscFunctionBegin;
517:   if (!mdn->Mvctx) PetscCall(MatSetUpMultiply_MPIDense(mat));
518:   PetscCall(VecGetArrayReadAndMemType(xx, &ax, &axmtype));
519:   PetscCall(VecGetArrayAndMemType(mdn->lvec, &ay, &aymtype));
520:   PetscCall(PetscSFBcastWithMemTypeBegin(mdn->Mvctx, MPIU_SCALAR, axmtype, ax, aymtype, ay, MPI_REPLACE));
521:   PetscCall(PetscSFBcastEnd(mdn->Mvctx, MPIU_SCALAR, ax, ay, MPI_REPLACE));
522:   PetscCall(VecRestoreArrayAndMemType(mdn->lvec, &ay));
523:   PetscCall(VecRestoreArrayReadAndMemType(xx, &ax));
524:   PetscCall((*mdn->A->ops->multadd)(mdn->A, mdn->lvec, yy, zz));
525:   PetscFunctionReturn(PETSC_SUCCESS);
526: }

528: static PetscErrorCode MatMultHermitianTransposeColumnRange_MPIDense(Mat A, Vec xx, Vec yy, PetscInt c_start, PetscInt c_end)
529: {
530:   Mat_MPIDense      *a = (Mat_MPIDense *)A->data;
531:   const PetscScalar *ax;
532:   PetscScalar       *ay;
533:   PetscMemType       axmtype, aymtype;

535:   PetscFunctionBegin;
536:   if (!a->Mvctx) PetscCall(MatSetUpMultiply_MPIDense(A));
537:   PetscCall(VecSet(yy, 0.0));
538:   PetscUseMethod(a->A, "MatMultHermitianTransposeColumnRange_C", (Mat, Vec, Vec, PetscInt, PetscInt), (a->A, xx, a->lvec, c_start, c_end));
539:   PetscCall(VecGetArrayReadAndMemType(a->lvec, &ax, &axmtype));
540:   PetscCall(VecGetArrayAndMemType(yy, &ay, &aymtype));
541:   PetscCall(PetscSFReduceWithMemTypeBegin(a->Mvctx, MPIU_SCALAR, axmtype, ax, aymtype, ay, MPIU_SUM));
542:   PetscCall(PetscSFReduceEnd(a->Mvctx, MPIU_SCALAR, ax, ay, MPIU_SUM));
543:   PetscCall(VecRestoreArrayReadAndMemType(a->lvec, &ax));
544:   PetscCall(VecRestoreArrayAndMemType(yy, &ay));
545:   PetscFunctionReturn(PETSC_SUCCESS);
546: }

548: static PetscErrorCode MatMultTransposeKernel_MPIDense(Mat A, Vec xx, Vec yy, PetscBool herm)
549: {
550:   Mat_MPIDense      *a = (Mat_MPIDense *)A->data;
551:   const PetscScalar *ax;
552:   PetscScalar       *ay;
553:   PetscMemType       axmtype, aymtype;

555:   PetscFunctionBegin;
556:   if (!a->Mvctx) PetscCall(MatSetUpMultiply_MPIDense(A));
557:   PetscCall(VecSet(yy, 0.0));
558:   if (herm) PetscCall((*a->A->ops->multhermitiantranspose)(a->A, xx, a->lvec));
559:   else PetscCall((*a->A->ops->multtranspose)(a->A, xx, a->lvec));
560:   PetscCall(VecGetArrayReadAndMemType(a->lvec, &ax, &axmtype));
561:   PetscCall(VecGetArrayAndMemType(yy, &ay, &aymtype));
562:   PetscCall(PetscSFReduceWithMemTypeBegin(a->Mvctx, MPIU_SCALAR, axmtype, ax, aymtype, ay, MPIU_SUM));
563:   PetscCall(PetscSFReduceEnd(a->Mvctx, MPIU_SCALAR, ax, ay, MPIU_SUM));
564:   PetscCall(VecRestoreArrayReadAndMemType(a->lvec, &ax));
565:   PetscCall(VecRestoreArrayAndMemType(yy, &ay));
566:   PetscFunctionReturn(PETSC_SUCCESS);
567: }

569: static PetscErrorCode MatMultHermitianTransposeAddColumnRange_MPIDense(Mat A, Vec xx, Vec yy, Vec zz, PetscInt c_start, PetscInt c_end)
570: {
571:   Mat_MPIDense      *a = (Mat_MPIDense *)A->data;
572:   const PetscScalar *ax;
573:   PetscScalar       *ay;
574:   PetscMemType       axmtype, aymtype;

576:   PetscFunctionBegin;
577:   if (!a->Mvctx) PetscCall(MatSetUpMultiply_MPIDense(A));
578:   PetscCall(VecCopy(yy, zz));
579:   PetscMPIInt rank;
580:   PetscCallMPI(MPI_Comm_rank(MPI_COMM_WORLD, &rank));
581:   PetscUseMethod(a->A, "MatMultHermitianTransposeColumnRange_C", (Mat, Vec, Vec, PetscInt, PetscInt), (a->A, xx, a->lvec, c_start, c_end));
582:   PetscCall(VecGetArrayReadAndMemType(a->lvec, &ax, &axmtype));
583:   PetscCall(VecGetArrayAndMemType(zz, &ay, &aymtype));
584:   PetscCall(PetscSFReduceWithMemTypeBegin(a->Mvctx, MPIU_SCALAR, axmtype, ax, aymtype, ay, MPIU_SUM));
585:   PetscCall(PetscSFReduceEnd(a->Mvctx, MPIU_SCALAR, ax, ay, MPIU_SUM));
586:   PetscCall(VecRestoreArrayReadAndMemType(a->lvec, &ax));
587:   PetscCall(VecRestoreArrayAndMemType(zz, &ay));
588:   PetscFunctionReturn(PETSC_SUCCESS);
589: }

591: static PetscErrorCode MatMultTransposeAddKernel_MPIDense(Mat A, Vec xx, Vec yy, Vec zz, PetscBool herm)
592: {
593:   Mat_MPIDense      *a = (Mat_MPIDense *)A->data;
594:   const PetscScalar *ax;
595:   PetscScalar       *ay;
596:   PetscMemType       axmtype, aymtype;

598:   PetscFunctionBegin;
599:   if (!a->Mvctx) PetscCall(MatSetUpMultiply_MPIDense(A));
600:   PetscCall(VecCopy(yy, zz));
601:   if (herm) PetscCall((*a->A->ops->multhermitiantranspose)(a->A, xx, a->lvec));
602:   else PetscCall((*a->A->ops->multtranspose)(a->A, xx, a->lvec));
603:   PetscCall(VecGetArrayReadAndMemType(a->lvec, &ax, &axmtype));
604:   PetscCall(VecGetArrayAndMemType(zz, &ay, &aymtype));
605:   PetscCall(PetscSFReduceWithMemTypeBegin(a->Mvctx, MPIU_SCALAR, axmtype, ax, aymtype, ay, MPIU_SUM));
606:   PetscCall(PetscSFReduceEnd(a->Mvctx, MPIU_SCALAR, ax, ay, MPIU_SUM));
607:   PetscCall(VecRestoreArrayReadAndMemType(a->lvec, &ax));
608:   PetscCall(VecRestoreArrayAndMemType(zz, &ay));
609:   PetscFunctionReturn(PETSC_SUCCESS);
610: }

612: static PetscErrorCode MatMultTranspose_MPIDense(Mat A, Vec xx, Vec yy)
613: {
614:   PetscFunctionBegin;
615:   PetscCall(MatMultTransposeKernel_MPIDense(A, xx, yy, PETSC_FALSE));
616:   PetscFunctionReturn(PETSC_SUCCESS);
617: }

619: static PetscErrorCode MatMultTransposeAdd_MPIDense(Mat A, Vec xx, Vec yy, Vec zz)
620: {
621:   PetscFunctionBegin;
622:   PetscCall(MatMultTransposeAddKernel_MPIDense(A, xx, yy, zz, PETSC_FALSE));
623:   PetscFunctionReturn(PETSC_SUCCESS);
624: }

626: static PetscErrorCode MatMultHermitianTranspose_MPIDense(Mat A, Vec xx, Vec yy)
627: {
628:   PetscFunctionBegin;
629:   PetscCall(MatMultTransposeKernel_MPIDense(A, xx, yy, PETSC_TRUE));
630:   PetscFunctionReturn(PETSC_SUCCESS);
631: }

633: static PetscErrorCode MatMultHermitianTransposeAdd_MPIDense(Mat A, Vec xx, Vec yy, Vec zz)
634: {
635:   PetscFunctionBegin;
636:   PetscCall(MatMultTransposeAddKernel_MPIDense(A, xx, yy, zz, PETSC_TRUE));
637:   PetscFunctionReturn(PETSC_SUCCESS);
638: }

640: PetscErrorCode MatGetDiagonal_MPIDense(Mat A, Vec v)
641: {
642:   Mat_MPIDense      *a = (Mat_MPIDense *)A->data;
643:   PetscInt           lda, len, i, nl, ng, m = A->rmap->n, radd;
644:   PetscScalar       *x;
645:   const PetscScalar *av;

647:   PetscFunctionBegin;
648:   PetscCall(VecGetArray(v, &x));
649:   PetscCall(VecGetSize(v, &ng));
650:   PetscCheck(ng == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming mat and vec");
651:   PetscCall(VecGetLocalSize(v, &nl));
652:   len  = PetscMin(a->A->rmap->n, a->A->cmap->n);
653:   radd = A->rmap->rstart * m;
654:   PetscCall(MatDenseGetArrayRead(a->A, &av));
655:   PetscCall(MatDenseGetLDA(a->A, &lda));
656:   for (i = 0; i < len; i++) x[i] = av[radd + i * lda + i];
657:   PetscCall(MatDenseRestoreArrayRead(a->A, &av));
658:   if (nl - i > 0) PetscCall(PetscArrayzero(x + i, nl - i));
659:   PetscCall(VecRestoreArray(v, &x));
660:   PetscFunctionReturn(PETSC_SUCCESS);
661: }

663: static PetscErrorCode MatDestroy_MPIDense(Mat mat)
664: {
665:   Mat_MPIDense *mdn = (Mat_MPIDense *)mat->data;

667:   PetscFunctionBegin;
668:   PetscCall(PetscLogObjectState((PetscObject)mat, "Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT, mat->rmap->N, mat->cmap->N));
669:   PetscCall(MatStashDestroy_Private(&mat->stash));
670:   PetscCheck(!mdn->vecinuse, PetscObjectComm((PetscObject)mat), PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
671:   PetscCheck(!mdn->matinuse, PetscObjectComm((PetscObject)mat), PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
672:   PetscCall(MatDestroy(&mdn->A));
673:   PetscCall(VecDestroy(&mdn->lvec));
674:   PetscCall(PetscSFDestroy(&mdn->Mvctx));
675:   PetscCall(VecDestroy(&mdn->cvec));
676:   PetscCall(MatDestroy(&mdn->cmat));

678:   PetscCall(PetscFree(mat->data));
679:   PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL));

681:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetLDA_C", NULL));
682:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseSetLDA_C", NULL));
683:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArray_C", NULL));
684:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArray_C", NULL));
685:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayRead_C", NULL));
686:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayRead_C", NULL));
687:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayWrite_C", NULL));
688:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayWrite_C", NULL));
689:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDensePlaceArray_C", NULL));
690:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseResetArray_C", NULL));
691:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseReplaceArray_C", NULL));
692:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpiaij_mpidense_C", NULL));
693:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpidense_mpiaij_C", NULL));
694: #if defined(PETSC_HAVE_ELEMENTAL)
695:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpidense_elemental_C", NULL));
696: #endif
697: #if defined(PETSC_HAVE_SCALAPACK)
698:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpidense_scalapack_C", NULL));
699: #endif
700:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPIDenseSetPreallocation_C", NULL));
701:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpiaij_mpidense_C", NULL));
702:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpidense_mpiaij_C", NULL));
703: #if defined(PETSC_HAVE_CUDA)
704:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpiaijcusparse_mpidense_C", NULL));
705:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpidense_mpiaijcusparse_C", NULL));
706:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpidense_mpidensecuda_C", NULL));
707:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpidensecuda_mpidense_C", NULL));
708:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpiaij_mpidensecuda_C", NULL));
709:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpiaijcusparse_mpidensecuda_C", NULL));
710:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpidensecuda_mpiaij_C", NULL));
711:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpidensecuda_mpiaijcusparse_C", NULL));
712:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseCUDAGetArray_C", NULL));
713:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseCUDAGetArrayRead_C", NULL));
714:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseCUDAGetArrayWrite_C", NULL));
715:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseCUDARestoreArray_C", NULL));
716:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseCUDARestoreArrayRead_C", NULL));
717:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseCUDARestoreArrayWrite_C", NULL));
718:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseCUDAPlaceArray_C", NULL));
719:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseCUDAResetArray_C", NULL));
720:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseCUDAReplaceArray_C", NULL));
721:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseCUDASetPreallocation_C", NULL));
722: #endif
723: #if defined(PETSC_HAVE_HIP)
724:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpiaijhipsparse_mpidense_C", NULL));
725:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpidense_mpiaijhipsparse_C", NULL));
726:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpidense_mpidensehip_C", NULL));
727:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpidensehip_mpidense_C", NULL));
728:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpiaij_mpidensehip_C", NULL));
729:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpiaijhipsparse_mpidensehip_C", NULL));
730:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpidensehip_mpiaij_C", NULL));
731:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpidensehip_mpiaijhipsparse_C", NULL));
732:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseHIPGetArray_C", NULL));
733:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseHIPGetArrayRead_C", NULL));
734:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseHIPGetArrayWrite_C", NULL));
735:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseHIPRestoreArray_C", NULL));
736:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseHIPRestoreArrayRead_C", NULL));
737:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseHIPRestoreArrayWrite_C", NULL));
738:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseHIPPlaceArray_C", NULL));
739:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseHIPResetArray_C", NULL));
740:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseHIPReplaceArray_C", NULL));
741:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseHIPSetPreallocation_C", NULL));
742: #endif
743:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumn_C", NULL));
744:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumn_C", NULL));
745:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVec_C", NULL));
746:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVec_C", NULL));
747:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecRead_C", NULL));
748:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecRead_C", NULL));
749:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecWrite_C", NULL));
750:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecWrite_C", NULL));
751:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetSubMatrix_C", NULL));
752:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreSubMatrix_C", NULL));
753:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultAddColumnRange_C", NULL));
754:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeColumnRange_C", NULL));
755:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeAddColumnRange_C", NULL));

757:   PetscCall(PetscObjectCompose((PetscObject)mat, "DiagonalBlock", NULL));
758:   PetscFunctionReturn(PETSC_SUCCESS);
759: }

761: #include <petscdraw.h>
762: static PetscErrorCode MatView_MPIDense_ASCIIorDraworSocket(Mat mat, PetscViewer viewer)
763: {
764:   Mat_MPIDense     *mdn = (Mat_MPIDense *)mat->data;
765:   PetscMPIInt       rank;
766:   PetscViewerType   vtype;
767:   PetscBool         iascii, isdraw;
768:   PetscViewer       sviewer;
769:   PetscViewerFormat format;

771:   PetscFunctionBegin;
772:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mat), &rank));
773:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
774:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
775:   if (iascii) {
776:     PetscCall(PetscViewerGetType(viewer, &vtype));
777:     PetscCall(PetscViewerGetFormat(viewer, &format));
778:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
779:       MatInfo info;
780:       PetscCall(MatGetInfo(mat, MAT_LOCAL, &info));
781:       PetscCall(PetscViewerASCIIPushSynchronized(viewer));
782:       PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "  [%d] local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT " \n", rank, mat->rmap->n, (PetscInt)info.nz_used, (PetscInt)info.nz_allocated,
783:                                                    (PetscInt)info.memory));
784:       PetscCall(PetscViewerFlush(viewer));
785:       PetscCall(PetscViewerASCIIPopSynchronized(viewer));
786:       if (mdn->Mvctx) PetscCall(PetscSFView(mdn->Mvctx, viewer));
787:       PetscFunctionReturn(PETSC_SUCCESS);
788:     } else if (format == PETSC_VIEWER_ASCII_INFO) {
789:       PetscFunctionReturn(PETSC_SUCCESS);
790:     }
791:   } else if (isdraw) {
792:     PetscDraw draw;
793:     PetscBool isnull;

795:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
796:     PetscCall(PetscDrawIsNull(draw, &isnull));
797:     if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
798:   }

800:   {
801:     /* assemble the entire matrix onto first processor. */
802:     Mat          A;
803:     PetscInt     M = mat->rmap->N, N = mat->cmap->N, m, row, i, nz;
804:     PetscInt    *cols;
805:     PetscScalar *vals;

807:     PetscCall(MatCreate(PetscObjectComm((PetscObject)mat), &A));
808:     if (rank == 0) {
809:       PetscCall(MatSetSizes(A, M, N, M, N));
810:     } else {
811:       PetscCall(MatSetSizes(A, 0, 0, M, N));
812:     }
813:     /* Since this is a temporary matrix, MATMPIDENSE instead of ((PetscObject)A)->type_name here is probably acceptable. */
814:     PetscCall(MatSetType(A, MATMPIDENSE));
815:     PetscCall(MatMPIDenseSetPreallocation(A, NULL));

817:     /* Copy the matrix ... This isn't the most efficient means,
818:        but it's quick for now */
819:     A->insertmode = INSERT_VALUES;

821:     row = mat->rmap->rstart;
822:     m   = mdn->A->rmap->n;
823:     for (i = 0; i < m; i++) {
824:       PetscCall(MatGetRow_MPIDense(mat, row, &nz, &cols, &vals));
825:       PetscCall(MatSetValues_MPIDense(A, 1, &row, nz, cols, vals, INSERT_VALUES));
826:       PetscCall(MatRestoreRow_MPIDense(mat, row, &nz, &cols, &vals));
827:       row++;
828:     }

830:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
831:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
832:     PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
833:     if (rank == 0) {
834:       PetscCall(PetscObjectSetName((PetscObject)((Mat_MPIDense *)A->data)->A, ((PetscObject)mat)->name));
835:       PetscCall(MatView_SeqDense(((Mat_MPIDense *)A->data)->A, sviewer));
836:     }
837:     PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
838:     PetscCall(MatDestroy(&A));
839:   }
840:   PetscFunctionReturn(PETSC_SUCCESS);
841: }

843: static PetscErrorCode MatView_MPIDense(Mat mat, PetscViewer viewer)
844: {
845:   PetscBool iascii, isbinary, isdraw, issocket;

847:   PetscFunctionBegin;
848:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
849:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
850:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSOCKET, &issocket));
851:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));

853:   if (iascii || issocket || isdraw) {
854:     PetscCall(MatView_MPIDense_ASCIIorDraworSocket(mat, viewer));
855:   } else if (isbinary) PetscCall(MatView_Dense_Binary(mat, viewer));
856:   PetscFunctionReturn(PETSC_SUCCESS);
857: }

859: static PetscErrorCode MatGetInfo_MPIDense(Mat A, MatInfoType flag, MatInfo *info)
860: {
861:   Mat_MPIDense  *mat = (Mat_MPIDense *)A->data;
862:   Mat            mdn = mat->A;
863:   PetscLogDouble isend[5], irecv[5];

865:   PetscFunctionBegin;
866:   info->block_size = 1.0;

868:   PetscCall(MatGetInfo(mdn, MAT_LOCAL, info));

870:   isend[0] = info->nz_used;
871:   isend[1] = info->nz_allocated;
872:   isend[2] = info->nz_unneeded;
873:   isend[3] = info->memory;
874:   isend[4] = info->mallocs;
875:   if (flag == MAT_LOCAL) {
876:     info->nz_used      = isend[0];
877:     info->nz_allocated = isend[1];
878:     info->nz_unneeded  = isend[2];
879:     info->memory       = isend[3];
880:     info->mallocs      = isend[4];
881:   } else if (flag == MAT_GLOBAL_MAX) {
882:     PetscCall(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)A)));

884:     info->nz_used      = irecv[0];
885:     info->nz_allocated = irecv[1];
886:     info->nz_unneeded  = irecv[2];
887:     info->memory       = irecv[3];
888:     info->mallocs      = irecv[4];
889:   } else if (flag == MAT_GLOBAL_SUM) {
890:     PetscCall(MPIU_Allreduce(isend, irecv, 5, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)A)));

892:     info->nz_used      = irecv[0];
893:     info->nz_allocated = irecv[1];
894:     info->nz_unneeded  = irecv[2];
895:     info->memory       = irecv[3];
896:     info->mallocs      = irecv[4];
897:   }
898:   info->fill_ratio_given  = 0; /* no parallel LU/ILU/Cholesky */
899:   info->fill_ratio_needed = 0;
900:   info->factor_mallocs    = 0;
901:   PetscFunctionReturn(PETSC_SUCCESS);
902: }

904: static PetscErrorCode MatSetOption_MPIDense(Mat A, MatOption op, PetscBool flg)
905: {
906:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

908:   PetscFunctionBegin;
909:   switch (op) {
910:   case MAT_NEW_NONZERO_LOCATIONS:
911:   case MAT_NEW_NONZERO_LOCATION_ERR:
912:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
913:     MatCheckPreallocated(A, 1);
914:     PetscCall(MatSetOption(a->A, op, flg));
915:     break;
916:   case MAT_ROW_ORIENTED:
917:     MatCheckPreallocated(A, 1);
918:     a->roworiented = flg;
919:     PetscCall(MatSetOption(a->A, op, flg));
920:     break;
921:   case MAT_FORCE_DIAGONAL_ENTRIES:
922:   case MAT_KEEP_NONZERO_PATTERN:
923:   case MAT_USE_HASH_TABLE:
924:   case MAT_SORTED_FULL:
925:     PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op]));
926:     break;
927:   case MAT_IGNORE_OFF_PROC_ENTRIES:
928:     a->donotstash = flg;
929:     break;
930:   case MAT_SYMMETRIC:
931:   case MAT_STRUCTURALLY_SYMMETRIC:
932:   case MAT_HERMITIAN:
933:   case MAT_SYMMETRY_ETERNAL:
934:   case MAT_STRUCTURAL_SYMMETRY_ETERNAL:
935:   case MAT_SPD:
936:   case MAT_IGNORE_LOWER_TRIANGULAR:
937:   case MAT_IGNORE_ZERO_ENTRIES:
938:   case MAT_SPD_ETERNAL:
939:     /* if the diagonal matrix is square it inherits some of the properties above */
940:     PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op]));
941:     break;
942:   default:
943:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unknown option %s", MatOptions[op]);
944:   }
945:   PetscFunctionReturn(PETSC_SUCCESS);
946: }

948: static PetscErrorCode MatDiagonalScale_MPIDense(Mat A, Vec ll, Vec rr)
949: {
950:   Mat_MPIDense      *mdn = (Mat_MPIDense *)A->data;
951:   const PetscScalar *l;
952:   PetscScalar        x, *v, *vv, *r;
953:   PetscInt           i, j, s2a, s3a, s2, s3, m = mdn->A->rmap->n, n = mdn->A->cmap->n, lda;

955:   PetscFunctionBegin;
956:   PetscCall(MatDenseGetArray(mdn->A, &vv));
957:   PetscCall(MatDenseGetLDA(mdn->A, &lda));
958:   PetscCall(MatGetLocalSize(A, &s2, &s3));
959:   if (ll) {
960:     PetscCall(VecGetLocalSize(ll, &s2a));
961:     PetscCheck(s2a == s2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vector non-conforming local size, %" PetscInt_FMT " != %" PetscInt_FMT, s2a, s2);
962:     PetscCall(VecGetArrayRead(ll, &l));
963:     for (i = 0; i < m; i++) {
964:       x = l[i];
965:       v = vv + i;
966:       for (j = 0; j < n; j++) {
967:         (*v) *= x;
968:         v += lda;
969:       }
970:     }
971:     PetscCall(VecRestoreArrayRead(ll, &l));
972:     PetscCall(PetscLogFlops(1.0 * n * m));
973:   }
974:   if (rr) {
975:     const PetscScalar *ar;

977:     PetscCall(VecGetLocalSize(rr, &s3a));
978:     PetscCheck(s3a == s3, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vec non-conforming local size, %" PetscInt_FMT " != %" PetscInt_FMT ".", s3a, s3);
979:     PetscCall(VecGetArrayRead(rr, &ar));
980:     if (!mdn->Mvctx) PetscCall(MatSetUpMultiply_MPIDense(A));
981:     PetscCall(VecGetArray(mdn->lvec, &r));
982:     PetscCall(PetscSFBcastBegin(mdn->Mvctx, MPIU_SCALAR, ar, r, MPI_REPLACE));
983:     PetscCall(PetscSFBcastEnd(mdn->Mvctx, MPIU_SCALAR, ar, r, MPI_REPLACE));
984:     PetscCall(VecRestoreArrayRead(rr, &ar));
985:     for (i = 0; i < n; i++) {
986:       x = r[i];
987:       v = vv + i * lda;
988:       for (j = 0; j < m; j++) (*v++) *= x;
989:     }
990:     PetscCall(VecRestoreArray(mdn->lvec, &r));
991:     PetscCall(PetscLogFlops(1.0 * n * m));
992:   }
993:   PetscCall(MatDenseRestoreArray(mdn->A, &vv));
994:   PetscFunctionReturn(PETSC_SUCCESS);
995: }

997: static PetscErrorCode MatNorm_MPIDense(Mat A, NormType type, PetscReal *nrm)
998: {
999:   Mat_MPIDense      *mdn = (Mat_MPIDense *)A->data;
1000:   PetscInt           i, j;
1001:   PetscMPIInt        size;
1002:   PetscReal          sum = 0.0;
1003:   const PetscScalar *av, *v;

1005:   PetscFunctionBegin;
1006:   PetscCall(MatDenseGetArrayRead(mdn->A, &av));
1007:   v = av;
1008:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1009:   if (size == 1) {
1010:     PetscCall(MatNorm(mdn->A, type, nrm));
1011:   } else {
1012:     if (type == NORM_FROBENIUS) {
1013:       for (i = 0; i < mdn->A->cmap->n * mdn->A->rmap->n; i++) {
1014:         sum += PetscRealPart(PetscConj(*v) * (*v));
1015:         v++;
1016:       }
1017:       PetscCall(MPIU_Allreduce(&sum, nrm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)A)));
1018:       *nrm = PetscSqrtReal(*nrm);
1019:       PetscCall(PetscLogFlops(2.0 * mdn->A->cmap->n * mdn->A->rmap->n));
1020:     } else if (type == NORM_1) {
1021:       PetscReal *tmp, *tmp2;
1022:       PetscCall(PetscCalloc2(A->cmap->N, &tmp, A->cmap->N, &tmp2));
1023:       *nrm = 0.0;
1024:       v    = av;
1025:       for (j = 0; j < mdn->A->cmap->n; j++) {
1026:         for (i = 0; i < mdn->A->rmap->n; i++) {
1027:           tmp[j] += PetscAbsScalar(*v);
1028:           v++;
1029:         }
1030:       }
1031:       PetscCall(MPIU_Allreduce(tmp, tmp2, A->cmap->N, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)A)));
1032:       for (j = 0; j < A->cmap->N; j++) {
1033:         if (tmp2[j] > *nrm) *nrm = tmp2[j];
1034:       }
1035:       PetscCall(PetscFree2(tmp, tmp2));
1036:       PetscCall(PetscLogFlops(A->cmap->n * A->rmap->n));
1037:     } else if (type == NORM_INFINITY) { /* max row norm */
1038:       PetscReal ntemp;
1039:       PetscCall(MatNorm(mdn->A, type, &ntemp));
1040:       PetscCall(MPIU_Allreduce(&ntemp, nrm, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)A)));
1041:     } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for two norm");
1042:   }
1043:   PetscCall(MatDenseRestoreArrayRead(mdn->A, &av));
1044:   PetscFunctionReturn(PETSC_SUCCESS);
1045: }

1047: static PetscErrorCode MatTranspose_MPIDense(Mat A, MatReuse reuse, Mat *matout)
1048: {
1049:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;
1050:   Mat           B;
1051:   PetscInt      M = A->rmap->N, N = A->cmap->N, m, n, *rwork, rstart = A->rmap->rstart;
1052:   PetscInt      j, i, lda;
1053:   PetscScalar  *v;

1055:   PetscFunctionBegin;
1056:   if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *matout));
1057:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1058:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
1059:     PetscCall(MatSetSizes(B, A->cmap->n, A->rmap->n, N, M));
1060:     PetscCall(MatSetType(B, ((PetscObject)A)->type_name));
1061:     PetscCall(MatMPIDenseSetPreallocation(B, NULL));
1062:   } else B = *matout;

1064:   m = a->A->rmap->n;
1065:   n = a->A->cmap->n;
1066:   PetscCall(MatDenseGetArrayRead(a->A, (const PetscScalar **)&v));
1067:   PetscCall(MatDenseGetLDA(a->A, &lda));
1068:   PetscCall(PetscMalloc1(m, &rwork));
1069:   for (i = 0; i < m; i++) rwork[i] = rstart + i;
1070:   for (j = 0; j < n; j++) {
1071:     PetscCall(MatSetValues(B, 1, &j, m, rwork, v, INSERT_VALUES));
1072:     v = PetscSafePointerPlusOffset(v, lda);
1073:   }
1074:   PetscCall(MatDenseRestoreArrayRead(a->A, (const PetscScalar **)&v));
1075:   PetscCall(PetscFree(rwork));
1076:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1077:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1078:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1079:     *matout = B;
1080:   } else {
1081:     PetscCall(MatHeaderMerge(A, &B));
1082:   }
1083:   PetscFunctionReturn(PETSC_SUCCESS);
1084: }

1086: static PetscErrorCode       MatDuplicate_MPIDense(Mat, MatDuplicateOption, Mat *);
1087: PETSC_INTERN PetscErrorCode MatScale_MPIDense(Mat, PetscScalar);

1089: static PetscErrorCode MatSetUp_MPIDense(Mat A)
1090: {
1091:   PetscFunctionBegin;
1092:   PetscCall(PetscLayoutSetUp(A->rmap));
1093:   PetscCall(PetscLayoutSetUp(A->cmap));
1094:   if (!A->preallocated) PetscCall(MatMPIDenseSetPreallocation(A, NULL));
1095:   PetscFunctionReturn(PETSC_SUCCESS);
1096: }

1098: static PetscErrorCode MatAXPY_MPIDense(Mat Y, PetscScalar alpha, Mat X, MatStructure str)
1099: {
1100:   Mat_MPIDense *A = (Mat_MPIDense *)Y->data, *B = (Mat_MPIDense *)X->data;

1102:   PetscFunctionBegin;
1103:   PetscCall(MatAXPY(A->A, alpha, B->A, str));
1104:   PetscFunctionReturn(PETSC_SUCCESS);
1105: }

1107: static PetscErrorCode MatConjugate_MPIDense(Mat mat)
1108: {
1109:   Mat_MPIDense *a = (Mat_MPIDense *)mat->data;

1111:   PetscFunctionBegin;
1112:   PetscCall(MatConjugate(a->A));
1113:   PetscFunctionReturn(PETSC_SUCCESS);
1114: }

1116: static PetscErrorCode MatRealPart_MPIDense(Mat A)
1117: {
1118:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

1120:   PetscFunctionBegin;
1121:   PetscCall(MatRealPart(a->A));
1122:   PetscFunctionReturn(PETSC_SUCCESS);
1123: }

1125: static PetscErrorCode MatImaginaryPart_MPIDense(Mat A)
1126: {
1127:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

1129:   PetscFunctionBegin;
1130:   PetscCall(MatImaginaryPart(a->A));
1131:   PetscFunctionReturn(PETSC_SUCCESS);
1132: }

1134: static PetscErrorCode MatGetColumnVector_MPIDense(Mat A, Vec v, PetscInt col)
1135: {
1136:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

1138:   PetscFunctionBegin;
1139:   PetscCheck(a->A, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Missing local matrix");
1140:   PetscCheck(a->A->ops->getcolumnvector, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Missing get column operation");
1141:   PetscCall((*a->A->ops->getcolumnvector)(a->A, v, col));
1142:   PetscFunctionReturn(PETSC_SUCCESS);
1143: }

1145: PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat, PetscInt, PetscReal *);

1147: static PetscErrorCode MatGetColumnReductions_MPIDense(Mat A, PetscInt type, PetscReal *reductions)
1148: {
1149:   PetscInt      i, m, n;
1150:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;
1151:   PetscReal    *work;

1153:   PetscFunctionBegin;
1154:   PetscCall(MatGetSize(A, &m, &n));
1155:   PetscCall(PetscMalloc1(n, &work));
1156:   if (type == REDUCTION_MEAN_REALPART) {
1157:     PetscCall(MatGetColumnReductions_SeqDense(a->A, (PetscInt)REDUCTION_SUM_REALPART, work));
1158:   } else if (type == REDUCTION_MEAN_IMAGINARYPART) {
1159:     PetscCall(MatGetColumnReductions_SeqDense(a->A, (PetscInt)REDUCTION_SUM_IMAGINARYPART, work));
1160:   } else {
1161:     PetscCall(MatGetColumnReductions_SeqDense(a->A, type, work));
1162:   }
1163:   if (type == NORM_2) {
1164:     for (i = 0; i < n; i++) work[i] *= work[i];
1165:   }
1166:   if (type == NORM_INFINITY) {
1167:     PetscCall(MPIU_Allreduce(work, reductions, n, MPIU_REAL, MPIU_MAX, A->hdr.comm));
1168:   } else {
1169:     PetscCall(MPIU_Allreduce(work, reductions, n, MPIU_REAL, MPIU_SUM, A->hdr.comm));
1170:   }
1171:   PetscCall(PetscFree(work));
1172:   if (type == NORM_2) {
1173:     for (i = 0; i < n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
1174:   } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
1175:     for (i = 0; i < n; i++) reductions[i] /= m;
1176:   }
1177:   PetscFunctionReturn(PETSC_SUCCESS);
1178: }

1180: static PetscErrorCode MatSetRandom_MPIDense(Mat x, PetscRandom rctx)
1181: {
1182:   Mat_MPIDense *d = (Mat_MPIDense *)x->data;

1184:   PetscFunctionBegin;
1185:   PetscCall(MatSetRandom(d->A, rctx));
1186: #if defined(PETSC_HAVE_DEVICE)
1187:   x->offloadmask = d->A->offloadmask;
1188: #endif
1189:   PetscFunctionReturn(PETSC_SUCCESS);
1190: }

1192: static PetscErrorCode MatMissingDiagonal_MPIDense(Mat A, PetscBool *missing, PetscInt *d)
1193: {
1194:   PetscFunctionBegin;
1195:   *missing = PETSC_FALSE;
1196:   PetscFunctionReturn(PETSC_SUCCESS);
1197: }

1199: static PetscErrorCode MatMatTransposeMultSymbolic_MPIDense_MPIDense(Mat, Mat, PetscReal, Mat);
1200: static PetscErrorCode MatMatTransposeMultNumeric_MPIDense_MPIDense(Mat, Mat, Mat);
1201: static PetscErrorCode MatTransposeMatMultSymbolic_MPIDense_MPIDense(Mat, Mat, PetscReal, Mat);
1202: static PetscErrorCode MatTransposeMatMultNumeric_MPIDense_MPIDense(Mat, Mat, Mat);
1203: static PetscErrorCode MatEqual_MPIDense(Mat, Mat, PetscBool *);
1204: static PetscErrorCode MatLoad_MPIDense(Mat, PetscViewer);
1205: static PetscErrorCode MatProductSetFromOptions_MPIDense(Mat);

1207: static struct _MatOps MatOps_Values = {MatSetValues_MPIDense,
1208:                                        MatGetRow_MPIDense,
1209:                                        MatRestoreRow_MPIDense,
1210:                                        MatMult_MPIDense,
1211:                                        /*  4*/ MatMultAdd_MPIDense,
1212:                                        MatMultTranspose_MPIDense,
1213:                                        MatMultTransposeAdd_MPIDense,
1214:                                        NULL,
1215:                                        NULL,
1216:                                        NULL,
1217:                                        /* 10*/ NULL,
1218:                                        NULL,
1219:                                        NULL,
1220:                                        NULL,
1221:                                        MatTranspose_MPIDense,
1222:                                        /* 15*/ MatGetInfo_MPIDense,
1223:                                        MatEqual_MPIDense,
1224:                                        MatGetDiagonal_MPIDense,
1225:                                        MatDiagonalScale_MPIDense,
1226:                                        MatNorm_MPIDense,
1227:                                        /* 20*/ MatAssemblyBegin_MPIDense,
1228:                                        MatAssemblyEnd_MPIDense,
1229:                                        MatSetOption_MPIDense,
1230:                                        MatZeroEntries_MPIDense,
1231:                                        /* 24*/ MatZeroRows_MPIDense,
1232:                                        NULL,
1233:                                        NULL,
1234:                                        NULL,
1235:                                        NULL,
1236:                                        /* 29*/ MatSetUp_MPIDense,
1237:                                        NULL,
1238:                                        NULL,
1239:                                        MatGetDiagonalBlock_MPIDense,
1240:                                        NULL,
1241:                                        /* 34*/ MatDuplicate_MPIDense,
1242:                                        NULL,
1243:                                        NULL,
1244:                                        NULL,
1245:                                        NULL,
1246:                                        /* 39*/ MatAXPY_MPIDense,
1247:                                        MatCreateSubMatrices_MPIDense,
1248:                                        NULL,
1249:                                        MatGetValues_MPIDense,
1250:                                        MatCopy_MPIDense,
1251:                                        /* 44*/ NULL,
1252:                                        MatScale_MPIDense,
1253:                                        MatShift_MPIDense,
1254:                                        NULL,
1255:                                        NULL,
1256:                                        /* 49*/ MatSetRandom_MPIDense,
1257:                                        NULL,
1258:                                        NULL,
1259:                                        NULL,
1260:                                        NULL,
1261:                                        /* 54*/ NULL,
1262:                                        NULL,
1263:                                        NULL,
1264:                                        NULL,
1265:                                        NULL,
1266:                                        /* 59*/ MatCreateSubMatrix_MPIDense,
1267:                                        MatDestroy_MPIDense,
1268:                                        MatView_MPIDense,
1269:                                        NULL,
1270:                                        NULL,
1271:                                        /* 64*/ NULL,
1272:                                        NULL,
1273:                                        NULL,
1274:                                        NULL,
1275:                                        NULL,
1276:                                        /* 69*/ NULL,
1277:                                        NULL,
1278:                                        NULL,
1279:                                        NULL,
1280:                                        NULL,
1281:                                        /* 74*/ NULL,
1282:                                        NULL,
1283:                                        NULL,
1284:                                        NULL,
1285:                                        NULL,
1286:                                        /* 79*/ NULL,
1287:                                        NULL,
1288:                                        NULL,
1289:                                        NULL,
1290:                                        /* 83*/ MatLoad_MPIDense,
1291:                                        NULL,
1292:                                        NULL,
1293:                                        NULL,
1294:                                        NULL,
1295:                                        NULL,
1296:                                        /* 89*/ NULL,
1297:                                        NULL,
1298:                                        NULL,
1299:                                        NULL,
1300:                                        NULL,
1301:                                        /* 94*/ NULL,
1302:                                        NULL,
1303:                                        MatMatTransposeMultSymbolic_MPIDense_MPIDense,
1304:                                        MatMatTransposeMultNumeric_MPIDense_MPIDense,
1305:                                        NULL,
1306:                                        /* 99*/ MatProductSetFromOptions_MPIDense,
1307:                                        NULL,
1308:                                        NULL,
1309:                                        MatConjugate_MPIDense,
1310:                                        NULL,
1311:                                        /*104*/ NULL,
1312:                                        MatRealPart_MPIDense,
1313:                                        MatImaginaryPart_MPIDense,
1314:                                        NULL,
1315:                                        NULL,
1316:                                        /*109*/ NULL,
1317:                                        NULL,
1318:                                        NULL,
1319:                                        MatGetColumnVector_MPIDense,
1320:                                        MatMissingDiagonal_MPIDense,
1321:                                        /*114*/ NULL,
1322:                                        NULL,
1323:                                        NULL,
1324:                                        NULL,
1325:                                        NULL,
1326:                                        /*119*/ NULL,
1327:                                        NULL,
1328:                                        MatMultHermitianTranspose_MPIDense,
1329:                                        MatMultHermitianTransposeAdd_MPIDense,
1330:                                        NULL,
1331:                                        /*124*/ NULL,
1332:                                        MatGetColumnReductions_MPIDense,
1333:                                        NULL,
1334:                                        NULL,
1335:                                        NULL,
1336:                                        /*129*/ NULL,
1337:                                        NULL,
1338:                                        MatTransposeMatMultSymbolic_MPIDense_MPIDense,
1339:                                        MatTransposeMatMultNumeric_MPIDense_MPIDense,
1340:                                        NULL,
1341:                                        /*134*/ NULL,
1342:                                        NULL,
1343:                                        NULL,
1344:                                        NULL,
1345:                                        NULL,
1346:                                        /*139*/ NULL,
1347:                                        NULL,
1348:                                        NULL,
1349:                                        NULL,
1350:                                        NULL,
1351:                                        MatCreateMPIMatConcatenateSeqMat_MPIDense,
1352:                                        /*145*/ NULL,
1353:                                        NULL,
1354:                                        NULL,
1355:                                        NULL,
1356:                                        NULL,
1357:                                        /*150*/ NULL,
1358:                                        NULL,
1359:                                        NULL};

1361: static PetscErrorCode MatMPIDenseSetPreallocation_MPIDense(Mat mat, PetscScalar *data)
1362: {
1363:   Mat_MPIDense *a     = (Mat_MPIDense *)mat->data;
1364:   MatType       mtype = MATSEQDENSE;

1366:   PetscFunctionBegin;
1367:   PetscCheck(!a->matinuse, PetscObjectComm((PetscObject)mat), PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1368:   PetscCall(PetscLayoutSetUp(mat->rmap));
1369:   PetscCall(PetscLayoutSetUp(mat->cmap));
1370:   if (!a->A) {
1371:     PetscCall(MatCreate(PETSC_COMM_SELF, &a->A));
1372:     PetscCall(MatSetSizes(a->A, mat->rmap->n, mat->cmap->N, mat->rmap->n, mat->cmap->N));
1373:   }
1374: #if defined(PETSC_HAVE_CUDA)
1375:   PetscBool iscuda;
1376:   PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPIDENSECUDA, &iscuda));
1377:   if (iscuda) mtype = MATSEQDENSECUDA;
1378: #endif
1379: #if defined(PETSC_HAVE_HIP)
1380:   PetscBool iship;
1381:   PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATMPIDENSEHIP, &iship));
1382:   if (iship) mtype = MATSEQDENSEHIP;
1383: #endif
1384:   PetscCall(MatSetType(a->A, mtype));
1385:   PetscCall(MatSeqDenseSetPreallocation(a->A, data));
1386: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1387:   mat->offloadmask = a->A->offloadmask;
1388: #endif
1389:   mat->preallocated = PETSC_TRUE;
1390:   mat->assembled    = PETSC_TRUE;
1391:   PetscFunctionReturn(PETSC_SUCCESS);
1392: }

1394: PETSC_INTERN PetscErrorCode MatConvert_MPIAIJ_MPIDense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1395: {
1396:   Mat B, C;

1398:   PetscFunctionBegin;
1399:   PetscCall(MatMPIAIJGetLocalMat(A, MAT_INITIAL_MATRIX, &C));
1400:   PetscCall(MatConvert_SeqAIJ_SeqDense(C, MATSEQDENSE, MAT_INITIAL_MATRIX, &B));
1401:   PetscCall(MatDestroy(&C));
1402:   if (reuse == MAT_REUSE_MATRIX) {
1403:     C = *newmat;
1404:   } else C = NULL;
1405:   PetscCall(MatCreateMPIMatConcatenateSeqMat(PetscObjectComm((PetscObject)A), B, A->cmap->n, !C ? MAT_INITIAL_MATRIX : MAT_REUSE_MATRIX, &C));
1406:   PetscCall(MatDestroy(&B));
1407:   if (reuse == MAT_INPLACE_MATRIX) {
1408:     PetscCall(MatHeaderReplace(A, &C));
1409:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = C;
1410:   PetscFunctionReturn(PETSC_SUCCESS);
1411: }

1413: static PetscErrorCode MatConvert_MPIDense_MPIAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1414: {
1415:   Mat B, C;

1417:   PetscFunctionBegin;
1418:   PetscCall(MatDenseGetLocalMatrix(A, &C));
1419:   PetscCall(MatConvert_SeqDense_SeqAIJ(C, MATSEQAIJ, MAT_INITIAL_MATRIX, &B));
1420:   if (reuse == MAT_REUSE_MATRIX) {
1421:     C = *newmat;
1422:   } else C = NULL;
1423:   PetscCall(MatCreateMPIMatConcatenateSeqMat(PetscObjectComm((PetscObject)A), B, A->cmap->n, !C ? MAT_INITIAL_MATRIX : MAT_REUSE_MATRIX, &C));
1424:   PetscCall(MatDestroy(&B));
1425:   if (reuse == MAT_INPLACE_MATRIX) {
1426:     PetscCall(MatHeaderReplace(A, &C));
1427:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = C;
1428:   PetscFunctionReturn(PETSC_SUCCESS);
1429: }

1431: #if defined(PETSC_HAVE_ELEMENTAL)
1432: PETSC_INTERN PetscErrorCode MatConvert_MPIDense_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1433: {
1434:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;
1435:   Mat           mat_elemental;
1436:   PetscScalar  *v;
1437:   PetscInt      m = A->rmap->n, N = A->cmap->N, rstart = A->rmap->rstart, i, *rows, *cols, lda;

1439:   PetscFunctionBegin;
1440:   if (reuse == MAT_REUSE_MATRIX) {
1441:     mat_elemental = *newmat;
1442:     PetscCall(MatZeroEntries(*newmat));
1443:   } else {
1444:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
1445:     PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, A->rmap->N, A->cmap->N));
1446:     PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
1447:     PetscCall(MatSetUp(mat_elemental));
1448:     PetscCall(MatSetOption(mat_elemental, MAT_ROW_ORIENTED, PETSC_FALSE));
1449:   }

1451:   PetscCall(PetscMalloc2(m, &rows, N, &cols));
1452:   for (i = 0; i < N; i++) cols[i] = i;
1453:   for (i = 0; i < m; i++) rows[i] = rstart + i;

1455:   /* PETSc-Elemental interface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
1456:   PetscCall(MatDenseGetArray(A, &v));
1457:   PetscCall(MatDenseGetLDA(a->A, &lda));
1458:   if (lda == m) PetscCall(MatSetValues(mat_elemental, m, rows, N, cols, v, ADD_VALUES));
1459:   else {
1460:     for (i = 0; i < N; i++) PetscCall(MatSetValues(mat_elemental, m, rows, 1, &i, v + lda * i, ADD_VALUES));
1461:   }
1462:   PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
1463:   PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
1464:   PetscCall(MatDenseRestoreArray(A, &v));
1465:   PetscCall(PetscFree2(rows, cols));

1467:   if (reuse == MAT_INPLACE_MATRIX) {
1468:     PetscCall(MatHeaderReplace(A, &mat_elemental));
1469:   } else {
1470:     *newmat = mat_elemental;
1471:   }
1472:   PetscFunctionReturn(PETSC_SUCCESS);
1473: }
1474: #endif

1476: static PetscErrorCode MatDenseGetColumn_MPIDense(Mat A, PetscInt col, PetscScalar **vals)
1477: {
1478:   Mat_MPIDense *mat = (Mat_MPIDense *)A->data;

1480:   PetscFunctionBegin;
1481:   PetscCall(MatDenseGetColumn(mat->A, col, vals));
1482:   PetscFunctionReturn(PETSC_SUCCESS);
1483: }

1485: static PetscErrorCode MatDenseRestoreColumn_MPIDense(Mat A, PetscScalar **vals)
1486: {
1487:   Mat_MPIDense *mat = (Mat_MPIDense *)A->data;

1489:   PetscFunctionBegin;
1490:   PetscCall(MatDenseRestoreColumn(mat->A, vals));
1491:   PetscFunctionReturn(PETSC_SUCCESS);
1492: }

1494: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_MPIDense(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
1495: {
1496:   Mat_MPIDense *mat;
1497:   PetscInt      m, nloc, N;

1499:   PetscFunctionBegin;
1500:   PetscCall(MatGetSize(inmat, &m, &N));
1501:   PetscCall(MatGetLocalSize(inmat, NULL, &nloc));
1502:   if (scall == MAT_INITIAL_MATRIX) { /* symbolic phase */
1503:     PetscInt sum;

1505:     if (n == PETSC_DECIDE) PetscCall(PetscSplitOwnership(comm, &n, &N));
1506:     /* Check sum(n) = N */
1507:     PetscCall(MPIU_Allreduce(&n, &sum, 1, MPIU_INT, MPI_SUM, comm));
1508:     PetscCheck(sum == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Sum of local columns %" PetscInt_FMT " != global columns %" PetscInt_FMT, sum, N);

1510:     PetscCall(MatCreateDense(comm, m, n, PETSC_DETERMINE, N, NULL, outmat));
1511:     PetscCall(MatSetOption(*outmat, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
1512:   }

1514:   /* numeric phase */
1515:   mat = (Mat_MPIDense *)(*outmat)->data;
1516:   PetscCall(MatCopy(inmat, mat->A, SAME_NONZERO_PATTERN));
1517:   PetscFunctionReturn(PETSC_SUCCESS);
1518: }

1520: PetscErrorCode MatDenseGetColumnVec_MPIDense(Mat A, PetscInt col, Vec *v)
1521: {
1522:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;
1523:   PetscInt      lda;

1525:   PetscFunctionBegin;
1526:   PetscCheck(!a->vecinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1527:   PetscCheck(!a->matinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1528:   if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
1529:   a->vecinuse = col + 1;
1530:   PetscCall(MatDenseGetLDA(a->A, &lda));
1531:   PetscCall(MatDenseGetArray(a->A, (PetscScalar **)&a->ptrinuse));
1532:   PetscCall(VecPlaceArray(a->cvec, PetscSafePointerPlusOffset(a->ptrinuse, (size_t)col * (size_t)lda)));
1533:   *v = a->cvec;
1534:   PetscFunctionReturn(PETSC_SUCCESS);
1535: }

1537: PetscErrorCode MatDenseRestoreColumnVec_MPIDense(Mat A, PetscInt col, Vec *v)
1538: {
1539:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

1541:   PetscFunctionBegin;
1542:   PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
1543:   PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
1544:   a->vecinuse = 0;
1545:   PetscCall(MatDenseRestoreArray(a->A, (PetscScalar **)&a->ptrinuse));
1546:   PetscCall(VecResetArray(a->cvec));
1547:   if (v) *v = NULL;
1548:   PetscFunctionReturn(PETSC_SUCCESS);
1549: }

1551: PetscErrorCode MatDenseGetColumnVecRead_MPIDense(Mat A, PetscInt col, Vec *v)
1552: {
1553:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;
1554:   PetscInt      lda;

1556:   PetscFunctionBegin;
1557:   PetscCheck(!a->vecinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1558:   PetscCheck(!a->matinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1559:   if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
1560:   a->vecinuse = col + 1;
1561:   PetscCall(MatDenseGetLDA(a->A, &lda));
1562:   PetscCall(MatDenseGetArrayRead(a->A, &a->ptrinuse));
1563:   PetscCall(VecPlaceArray(a->cvec, PetscSafePointerPlusOffset(a->ptrinuse, (size_t)col * (size_t)lda)));
1564:   PetscCall(VecLockReadPush(a->cvec));
1565:   *v = a->cvec;
1566:   PetscFunctionReturn(PETSC_SUCCESS);
1567: }

1569: PetscErrorCode MatDenseRestoreColumnVecRead_MPIDense(Mat A, PetscInt col, Vec *v)
1570: {
1571:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

1573:   PetscFunctionBegin;
1574:   PetscCheck(a->vecinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
1575:   PetscCheck(a->cvec, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing internal column vector");
1576:   a->vecinuse = 0;
1577:   PetscCall(MatDenseRestoreArrayRead(a->A, &a->ptrinuse));
1578:   PetscCall(VecLockReadPop(a->cvec));
1579:   PetscCall(VecResetArray(a->cvec));
1580:   if (v) *v = NULL;
1581:   PetscFunctionReturn(PETSC_SUCCESS);
1582: }

1584: PetscErrorCode MatDenseGetColumnVecWrite_MPIDense(Mat A, PetscInt col, Vec *v)
1585: {
1586:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;
1587:   PetscInt      lda;

1589:   PetscFunctionBegin;
1590:   PetscCheck(!a->vecinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1591:   PetscCheck(!a->matinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1592:   if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
1593:   a->vecinuse = col + 1;
1594:   PetscCall(MatDenseGetLDA(a->A, &lda));
1595:   PetscCall(MatDenseGetArrayWrite(a->A, (PetscScalar **)&a->ptrinuse));
1596:   PetscCall(VecPlaceArray(a->cvec, PetscSafePointerPlusOffset(a->ptrinuse, (size_t)col * (size_t)lda)));
1597:   *v = a->cvec;
1598:   PetscFunctionReturn(PETSC_SUCCESS);
1599: }

1601: PetscErrorCode MatDenseRestoreColumnVecWrite_MPIDense(Mat A, PetscInt col, Vec *v)
1602: {
1603:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;

1605:   PetscFunctionBegin;
1606:   PetscCheck(a->vecinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
1607:   PetscCheck(a->cvec, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing internal column vector");
1608:   a->vecinuse = 0;
1609:   PetscCall(MatDenseRestoreArrayWrite(a->A, (PetscScalar **)&a->ptrinuse));
1610:   PetscCall(VecResetArray(a->cvec));
1611:   if (v) *v = NULL;
1612:   PetscFunctionReturn(PETSC_SUCCESS);
1613: }

1615: static PetscErrorCode MatDenseGetSubMatrix_MPIDense(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
1616: {
1617:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;
1618:   Mat_MPIDense *c;
1619:   MPI_Comm      comm;
1620:   PetscInt      pbegin, pend;

1622:   PetscFunctionBegin;
1623:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
1624:   PetscCheck(!a->vecinuse, comm, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1625:   PetscCheck(!a->matinuse, comm, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1626:   pbegin = PetscMax(0, PetscMin(A->rmap->rend, rbegin) - A->rmap->rstart);
1627:   pend   = PetscMin(A->rmap->n, PetscMax(0, rend - A->rmap->rstart));
1628:   if (!a->cmat) {
1629:     PetscCall(MatCreate(comm, &a->cmat));
1630:     PetscCall(MatSetType(a->cmat, ((PetscObject)A)->type_name));
1631:     if (rend - rbegin == A->rmap->N) PetscCall(PetscLayoutReference(A->rmap, &a->cmat->rmap));
1632:     else {
1633:       PetscCall(PetscLayoutSetLocalSize(a->cmat->rmap, pend - pbegin));
1634:       PetscCall(PetscLayoutSetSize(a->cmat->rmap, rend - rbegin));
1635:       PetscCall(PetscLayoutSetUp(a->cmat->rmap));
1636:     }
1637:     PetscCall(PetscLayoutSetSize(a->cmat->cmap, cend - cbegin));
1638:     PetscCall(PetscLayoutSetUp(a->cmat->cmap));
1639:   } else {
1640:     PetscBool same = (PetscBool)(rend - rbegin == a->cmat->rmap->N);
1641:     if (same && a->cmat->rmap->N != A->rmap->N) {
1642:       same = (PetscBool)(pend - pbegin == a->cmat->rmap->n);
1643:       PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &same, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
1644:     }
1645:     if (!same) {
1646:       PetscCall(PetscLayoutDestroy(&a->cmat->rmap));
1647:       PetscCall(PetscLayoutCreate(comm, &a->cmat->rmap));
1648:       PetscCall(PetscLayoutSetLocalSize(a->cmat->rmap, pend - pbegin));
1649:       PetscCall(PetscLayoutSetSize(a->cmat->rmap, rend - rbegin));
1650:       PetscCall(PetscLayoutSetUp(a->cmat->rmap));
1651:     }
1652:     if (cend - cbegin != a->cmat->cmap->N) {
1653:       PetscCall(PetscLayoutDestroy(&a->cmat->cmap));
1654:       PetscCall(PetscLayoutCreate(comm, &a->cmat->cmap));
1655:       PetscCall(PetscLayoutSetSize(a->cmat->cmap, cend - cbegin));
1656:       PetscCall(PetscLayoutSetUp(a->cmat->cmap));
1657:     }
1658:   }
1659:   c = (Mat_MPIDense *)a->cmat->data;
1660:   PetscCheck(!c->A, comm, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1661:   PetscCall(MatDenseGetSubMatrix(a->A, pbegin, pend, cbegin, cend, &c->A));

1663:   a->cmat->preallocated = PETSC_TRUE;
1664:   a->cmat->assembled    = PETSC_TRUE;
1665: #if defined(PETSC_HAVE_DEVICE)
1666:   a->cmat->offloadmask = c->A->offloadmask;
1667: #endif
1668:   a->matinuse = cbegin + 1;
1669:   *v          = a->cmat;
1670:   PetscFunctionReturn(PETSC_SUCCESS);
1671: }

1673: static PetscErrorCode MatDenseRestoreSubMatrix_MPIDense(Mat A, Mat *v)
1674: {
1675:   Mat_MPIDense *a = (Mat_MPIDense *)A->data;
1676:   Mat_MPIDense *c;

1678:   PetscFunctionBegin;
1679:   PetscCheck(a->matinuse, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Need to call MatDenseGetSubMatrix() first");
1680:   PetscCheck(a->cmat, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing internal matrix");
1681:   PetscCheck(*v == a->cmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not the matrix obtained from MatDenseGetSubMatrix()");
1682:   a->matinuse = 0;
1683:   c           = (Mat_MPIDense *)a->cmat->data;
1684:   PetscCall(MatDenseRestoreSubMatrix(a->A, &c->A));
1685:   if (v) *v = NULL;
1686: #if defined(PETSC_HAVE_DEVICE)
1687:   A->offloadmask = a->A->offloadmask;
1688: #endif
1689:   PetscFunctionReturn(PETSC_SUCCESS);
1690: }

1692: /*MC
1693:    MATMPIDENSE - MATMPIDENSE = "mpidense" - A matrix type to be used for distributed dense matrices.

1695:    Options Database Key:
1696: . -mat_type mpidense - sets the matrix type to `MATMPIDENSE` during a call to `MatSetFromOptions()`

1698:   Level: beginner

1700: .seealso: [](ch_matrices), `Mat`, `MatCreateDense()`, `MATSEQDENSE`, `MATDENSE`
1701: M*/
1702: PetscErrorCode MatCreate_MPIDense(Mat mat)
1703: {
1704:   Mat_MPIDense *a;

1706:   PetscFunctionBegin;
1707:   PetscCall(PetscNew(&a));
1708:   mat->data   = (void *)a;
1709:   mat->ops[0] = MatOps_Values;

1711:   mat->insertmode = NOT_SET_VALUES;

1713:   /* build cache for off array entries formed */
1714:   a->donotstash = PETSC_FALSE;

1716:   PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)mat), 1, &mat->stash));

1718:   /* stuff used for matrix vector multiply */
1719:   a->lvec        = NULL;
1720:   a->Mvctx       = NULL;
1721:   a->roworiented = PETSC_TRUE;

1723:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetLDA_C", MatDenseGetLDA_MPIDense));
1724:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseSetLDA_C", MatDenseSetLDA_MPIDense));
1725:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArray_C", MatDenseGetArray_MPIDense));
1726:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArray_C", MatDenseRestoreArray_MPIDense));
1727:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayRead_C", MatDenseGetArrayRead_MPIDense));
1728:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayRead_C", MatDenseRestoreArrayRead_MPIDense));
1729:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayWrite_C", MatDenseGetArrayWrite_MPIDense));
1730:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayWrite_C", MatDenseRestoreArrayWrite_MPIDense));
1731:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDensePlaceArray_C", MatDensePlaceArray_MPIDense));
1732:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseResetArray_C", MatDenseResetArray_MPIDense));
1733:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseReplaceArray_C", MatDenseReplaceArray_MPIDense));
1734:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVec_C", MatDenseGetColumnVec_MPIDense));
1735:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVec_C", MatDenseRestoreColumnVec_MPIDense));
1736:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecRead_C", MatDenseGetColumnVecRead_MPIDense));
1737:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecRead_C", MatDenseRestoreColumnVecRead_MPIDense));
1738:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecWrite_C", MatDenseGetColumnVecWrite_MPIDense));
1739:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecWrite_C", MatDenseRestoreColumnVecWrite_MPIDense));
1740:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetSubMatrix_C", MatDenseGetSubMatrix_MPIDense));
1741:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreSubMatrix_C", MatDenseRestoreSubMatrix_MPIDense));
1742:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpiaij_mpidense_C", MatConvert_MPIAIJ_MPIDense));
1743:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpidense_mpiaij_C", MatConvert_MPIDense_MPIAIJ));
1744: #if defined(PETSC_HAVE_ELEMENTAL)
1745:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpidense_elemental_C", MatConvert_MPIDense_Elemental));
1746: #endif
1747: #if defined(PETSC_HAVE_SCALAPACK)
1748:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpidense_scalapack_C", MatConvert_Dense_ScaLAPACK));
1749: #endif
1750: #if defined(PETSC_HAVE_CUDA)
1751:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpidense_mpidensecuda_C", MatConvert_MPIDense_MPIDenseCUDA));
1752: #endif
1753:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMPIDenseSetPreallocation_C", MatMPIDenseSetPreallocation_MPIDense));
1754:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpiaij_mpidense_C", MatProductSetFromOptions_MPIAIJ_MPIDense));
1755:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpidense_mpiaij_C", MatProductSetFromOptions_MPIDense_MPIAIJ));
1756: #if defined(PETSC_HAVE_CUDA)
1757:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpiaijcusparse_mpidense_C", MatProductSetFromOptions_MPIAIJ_MPIDense));
1758:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpidense_mpiaijcusparse_C", MatProductSetFromOptions_MPIDense_MPIAIJ));
1759: #endif
1760: #if defined(PETSC_HAVE_HIP)
1761:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_mpidense_mpidensehip_C", MatConvert_MPIDense_MPIDenseHIP));
1762:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpiaijhipsparse_mpidense_C", MatProductSetFromOptions_MPIAIJ_MPIDense));
1763:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_mpidense_mpiaijhipsparse_C", MatProductSetFromOptions_MPIDense_MPIAIJ));
1764: #endif
1765:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumn_C", MatDenseGetColumn_MPIDense));
1766:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumn_C", MatDenseRestoreColumn_MPIDense));
1767:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultAddColumnRange_C", MatMultAddColumnRange_MPIDense));
1768:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeColumnRange_C", MatMultHermitianTransposeColumnRange_MPIDense));
1769:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeAddColumnRange_C", MatMultHermitianTransposeAddColumnRange_MPIDense));
1770:   PetscCall(PetscObjectChangeTypeName((PetscObject)mat, MATMPIDENSE));
1771:   PetscFunctionReturn(PETSC_SUCCESS);
1772: }

1774: /*MC
1775:    MATDENSE - MATDENSE = "dense" - A matrix type to be used for dense matrices.

1777:    This matrix type is identical to `MATSEQDENSE` when constructed with a single process communicator,
1778:    and `MATMPIDENSE` otherwise.

1780:    Options Database Key:
1781: . -mat_type dense - sets the matrix type to `MATDENSE` during a call to `MatSetFromOptions()`

1783:   Level: beginner

1785: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MATMPIDENSE`, `MATDENSECUDA`, `MATDENSEHIP`
1786: M*/

1788: /*@C
1789:   MatMPIDenseSetPreallocation - Sets the array used to store the matrix entries

1791:   Collective

1793:   Input Parameters:
1794: + B    - the matrix
1795: - data - optional location of matrix data.  Set to `NULL` for PETSc
1796:          to control all matrix memory allocation.

1798:   Level: intermediate

1800:   Notes:
1801:   The dense format is fully compatible with standard Fortran
1802:   storage by columns.

1804:   The data input variable is intended primarily for Fortran programmers
1805:   who wish to allocate their own matrix memory space.  Most users should
1806:   set `data` to `NULL`.

1808: .seealso: [](ch_matrices), `Mat`, `MATMPIDENSE`, `MatCreate()`, `MatCreateSeqDense()`, `MatSetValues()`
1809: @*/
1810: PetscErrorCode MatMPIDenseSetPreallocation(Mat B, PetscScalar *data)
1811: {
1812:   PetscFunctionBegin;
1814:   PetscTryMethod(B, "MatMPIDenseSetPreallocation_C", (Mat, PetscScalar *), (B, data));
1815:   PetscFunctionReturn(PETSC_SUCCESS);
1816: }

1818: /*@
1819:   MatDensePlaceArray - Allows one to replace the array in a `MATDENSE` matrix with an
1820:   array provided by the user. This is useful to avoid copying an array
1821:   into a matrix

1823:   Not Collective

1825:   Input Parameters:
1826: + mat   - the matrix
1827: - array - the array in column major order

1829:   Level: developer

1831:   Note:
1832:   You can return to the original array with a call to `MatDenseResetArray()`. The user is responsible for freeing this array; it will not be
1833:   freed when the matrix is destroyed.

1835: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArray()`, `MatDenseResetArray()`, `VecPlaceArray()`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`,
1836:           `MatDenseReplaceArray()`
1837: @*/
1838: PetscErrorCode MatDensePlaceArray(Mat mat, const PetscScalar *array)
1839: {
1840:   PetscFunctionBegin;
1842:   PetscUseMethod(mat, "MatDensePlaceArray_C", (Mat, const PetscScalar *), (mat, array));
1843:   PetscCall(PetscObjectStateIncrease((PetscObject)mat));
1844: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1845:   mat->offloadmask = PETSC_OFFLOAD_CPU;
1846: #endif
1847:   PetscFunctionReturn(PETSC_SUCCESS);
1848: }

1850: /*@
1851:   MatDenseResetArray - Resets the matrix array to that it previously had before the call to `MatDensePlaceArray()`

1853:   Not Collective

1855:   Input Parameter:
1856: . mat - the matrix

1858:   Level: developer

1860:   Note:
1861:   You can only call this after a call to `MatDensePlaceArray()`

1863: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArray()`, `MatDensePlaceArray()`, `VecPlaceArray()`, `VecGetArray()`, `VecRestoreArray()`, `VecReplaceArray()`, `VecResetArray()`
1864: @*/
1865: PetscErrorCode MatDenseResetArray(Mat mat)
1866: {
1867:   PetscFunctionBegin;
1869:   PetscUseMethod(mat, "MatDenseResetArray_C", (Mat), (mat));
1870:   PetscCall(PetscObjectStateIncrease((PetscObject)mat));
1871:   PetscFunctionReturn(PETSC_SUCCESS);
1872: }

1874: /*@
1875:   MatDenseReplaceArray - Allows one to replace the array in a dense matrix with an
1876:   array provided by the user. This is useful to avoid copying an array
1877:   into a matrix

1879:   Not Collective

1881:   Input Parameters:
1882: + mat   - the matrix
1883: - array - the array in column major order

1885:   Level: developer

1887:   Note:
1888:   The memory passed in MUST be obtained with `PetscMalloc()` and CANNOT be
1889:   freed by the user. It will be freed when the matrix is destroyed.

1891: .seealso: [](ch_matrices), `Mat`, `MatDensePlaceArray()`, `MatDenseGetArray()`, `VecReplaceArray()`
1892: @*/
1893: PetscErrorCode MatDenseReplaceArray(Mat mat, const PetscScalar *array)
1894: {
1895:   PetscFunctionBegin;
1897:   PetscUseMethod(mat, "MatDenseReplaceArray_C", (Mat, const PetscScalar *), (mat, array));
1898:   PetscCall(PetscObjectStateIncrease((PetscObject)mat));
1899: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1900:   mat->offloadmask = PETSC_OFFLOAD_CPU;
1901: #endif
1902:   PetscFunctionReturn(PETSC_SUCCESS);
1903: }

1905: /*@C
1906:   MatCreateDense - Creates a matrix in `MATDENSE` format.

1908:   Collective

1910:   Input Parameters:
1911: + comm - MPI communicator
1912: . m    - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
1913: . n    - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
1914: . M    - number of global rows (or `PETSC_DECIDE` to have calculated if `m` is given)
1915: . N    - number of global columns (or `PETSC_DECIDE` to have calculated if `n` is given)
1916: - data - optional location of matrix data.  Set data to `NULL` (`PETSC_NULL_SCALAR` for Fortran users) for PETSc
1917:    to control all matrix memory allocation.

1919:   Output Parameter:
1920: . A - the matrix

1922:   Level: intermediate

1924:   Notes:
1925:   The dense format is fully compatible with standard Fortran
1926:   storage by columns.

1928:   Although local portions of the matrix are stored in column-major
1929:   order, the matrix is partitioned across MPI ranks by row.

1931:   The data input variable is intended primarily for Fortran programmers
1932:   who wish to allocate their own matrix memory space.  Most users should
1933:   set `data` to `NULL` (`PETSC_NULL_SCALAR` for Fortran users).

1935:   The user MUST specify either the local or global matrix dimensions
1936:   (possibly both).

1938: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatCreate()`, `MatCreateSeqDense()`, `MatSetValues()`
1939: @*/
1940: PetscErrorCode MatCreateDense(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscScalar *data, Mat *A)
1941: {
1942:   PetscFunctionBegin;
1943:   PetscCall(MatCreate(comm, A));
1944:   PetscCall(MatSetSizes(*A, m, n, M, N));
1945:   PetscCall(MatSetType(*A, MATDENSE));
1946:   PetscCall(MatSeqDenseSetPreallocation(*A, data));
1947:   PetscCall(MatMPIDenseSetPreallocation(*A, data));
1948:   PetscFunctionReturn(PETSC_SUCCESS);
1949: }

1951: static PetscErrorCode MatDuplicate_MPIDense(Mat A, MatDuplicateOption cpvalues, Mat *newmat)
1952: {
1953:   Mat           mat;
1954:   Mat_MPIDense *a, *oldmat = (Mat_MPIDense *)A->data;

1956:   PetscFunctionBegin;
1957:   *newmat = NULL;
1958:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat));
1959:   PetscCall(MatSetSizes(mat, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
1960:   PetscCall(MatSetType(mat, ((PetscObject)A)->type_name));
1961:   a = (Mat_MPIDense *)mat->data;

1963:   mat->factortype   = A->factortype;
1964:   mat->assembled    = PETSC_TRUE;
1965:   mat->preallocated = PETSC_TRUE;

1967:   mat->insertmode = NOT_SET_VALUES;
1968:   a->donotstash   = oldmat->donotstash;

1970:   PetscCall(PetscLayoutReference(A->rmap, &mat->rmap));
1971:   PetscCall(PetscLayoutReference(A->cmap, &mat->cmap));

1973:   PetscCall(MatDuplicate(oldmat->A, cpvalues, &a->A));

1975:   *newmat = mat;
1976:   PetscFunctionReturn(PETSC_SUCCESS);
1977: }

1979: static PetscErrorCode MatLoad_MPIDense(Mat newMat, PetscViewer viewer)
1980: {
1981:   PetscBool isbinary;
1982: #if defined(PETSC_HAVE_HDF5)
1983:   PetscBool ishdf5;
1984: #endif

1986:   PetscFunctionBegin;
1989:   /* force binary viewer to load .info file if it has not yet done so */
1990:   PetscCall(PetscViewerSetUp(viewer));
1991:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1992: #if defined(PETSC_HAVE_HDF5)
1993:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1994: #endif
1995:   if (isbinary) {
1996:     PetscCall(MatLoad_Dense_Binary(newMat, viewer));
1997: #if defined(PETSC_HAVE_HDF5)
1998:   } else if (ishdf5) {
1999:     PetscCall(MatLoad_Dense_HDF5(newMat, viewer));
2000: #endif
2001:   } else SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "Viewer type %s not yet supported for reading %s matrices", ((PetscObject)viewer)->type_name, ((PetscObject)newMat)->type_name);
2002:   PetscFunctionReturn(PETSC_SUCCESS);
2003: }

2005: static PetscErrorCode MatEqual_MPIDense(Mat A, Mat B, PetscBool *flag)
2006: {
2007:   Mat_MPIDense *matB = (Mat_MPIDense *)B->data, *matA = (Mat_MPIDense *)A->data;
2008:   Mat           a, b;

2010:   PetscFunctionBegin;
2011:   a = matA->A;
2012:   b = matB->A;
2013:   PetscCall(MatEqual(a, b, flag));
2014:   PetscCall(MPIU_Allreduce(MPI_IN_PLACE, flag, 1, MPIU_BOOL, MPI_LAND, PetscObjectComm((PetscObject)A)));
2015:   PetscFunctionReturn(PETSC_SUCCESS);
2016: }

2018: static PetscErrorCode MatDestroy_MatTransMatMult_MPIDense_MPIDense(void *data)
2019: {
2020:   Mat_TransMatMultDense *atb = (Mat_TransMatMultDense *)data;

2022:   PetscFunctionBegin;
2023:   PetscCall(PetscFree2(atb->sendbuf, atb->recvcounts));
2024:   PetscCall(MatDestroy(&atb->atb));
2025:   PetscCall(PetscFree(atb));
2026:   PetscFunctionReturn(PETSC_SUCCESS);
2027: }

2029: static PetscErrorCode MatDestroy_MatMatTransMult_MPIDense_MPIDense(void *data)
2030: {
2031:   Mat_MatTransMultDense *abt = (Mat_MatTransMultDense *)data;

2033:   PetscFunctionBegin;
2034:   PetscCall(PetscFree2(abt->buf[0], abt->buf[1]));
2035:   PetscCall(PetscFree2(abt->recvcounts, abt->recvdispls));
2036:   PetscCall(PetscFree(abt));
2037:   PetscFunctionReturn(PETSC_SUCCESS);
2038: }

2040: static PetscErrorCode MatTransposeMatMultNumeric_MPIDense_MPIDense(Mat A, Mat B, Mat C)
2041: {
2042:   Mat_MPIDense          *a = (Mat_MPIDense *)A->data, *b = (Mat_MPIDense *)B->data, *c = (Mat_MPIDense *)C->data;
2043:   Mat_TransMatMultDense *atb;
2044:   MPI_Comm               comm;
2045:   PetscMPIInt            size, *recvcounts;
2046:   PetscScalar           *carray, *sendbuf;
2047:   const PetscScalar     *atbarray;
2048:   PetscInt               i, cN = C->cmap->N, proc, k, j, lda;
2049:   const PetscInt        *ranges;

2051:   PetscFunctionBegin;
2052:   MatCheckProduct(C, 3);
2053:   PetscCheck(C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty");
2054:   atb        = (Mat_TransMatMultDense *)C->product->data;
2055:   recvcounts = atb->recvcounts;
2056:   sendbuf    = atb->sendbuf;

2058:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
2059:   PetscCallMPI(MPI_Comm_size(comm, &size));

2061:   /* compute atbarray = aseq^T * bseq */
2062:   PetscCall(MatTransposeMatMult(a->A, b->A, atb->atb ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, PETSC_DEFAULT, &atb->atb));

2064:   PetscCall(MatGetOwnershipRanges(C, &ranges));

2066:   /* arrange atbarray into sendbuf */
2067:   PetscCall(MatDenseGetArrayRead(atb->atb, &atbarray));
2068:   PetscCall(MatDenseGetLDA(atb->atb, &lda));
2069:   for (proc = 0, k = 0; proc < size; proc++) {
2070:     for (j = 0; j < cN; j++) {
2071:       for (i = ranges[proc]; i < ranges[proc + 1]; i++) sendbuf[k++] = atbarray[i + j * lda];
2072:     }
2073:   }
2074:   PetscCall(MatDenseRestoreArrayRead(atb->atb, &atbarray));

2076:   /* sum all atbarray to local values of C */
2077:   PetscCall(MatDenseGetArrayWrite(c->A, &carray));
2078:   PetscCallMPI(MPI_Reduce_scatter(sendbuf, carray, recvcounts, MPIU_SCALAR, MPIU_SUM, comm));
2079:   PetscCall(MatDenseRestoreArrayWrite(c->A, &carray));
2080:   PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
2081:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
2082:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
2083:   PetscFunctionReturn(PETSC_SUCCESS);
2084: }

2086: static PetscErrorCode MatTransposeMatMultSymbolic_MPIDense_MPIDense(Mat A, Mat B, PetscReal fill, Mat C)
2087: {
2088:   MPI_Comm               comm;
2089:   PetscMPIInt            size;
2090:   PetscInt               cm = A->cmap->n, cM, cN = B->cmap->N;
2091:   Mat_TransMatMultDense *atb;
2092:   PetscBool              cisdense = PETSC_FALSE;
2093:   PetscInt               i;
2094:   const PetscInt        *ranges;

2096:   PetscFunctionBegin;
2097:   MatCheckProduct(C, 4);
2098:   PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
2099:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
2100:   if (A->rmap->rstart != B->rmap->rstart || A->rmap->rend != B->rmap->rend) {
2101:     SETERRQ(comm, PETSC_ERR_ARG_SIZ, "Matrix local dimensions are incompatible, A (%" PetscInt_FMT ", %" PetscInt_FMT ") != B (%" PetscInt_FMT ",%" PetscInt_FMT ")", A->rmap->rstart, A->rmap->rend, B->rmap->rstart, B->rmap->rend);
2102:   }

2104:   /* create matrix product C */
2105:   PetscCall(MatSetSizes(C, cm, B->cmap->n, A->cmap->N, B->cmap->N));
2106: #if defined(PETSC_HAVE_CUDA)
2107:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATMPIDENSE, MATMPIDENSECUDA, ""));
2108: #endif
2109: #if defined(PETSC_HAVE_HIP)
2110:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATMPIDENSE, MATMPIDENSEHIP, ""));
2111: #endif
2112:   if (!cisdense) PetscCall(MatSetType(C, ((PetscObject)A)->type_name));
2113:   PetscCall(MatSetUp(C));

2115:   /* create data structure for reuse C */
2116:   PetscCallMPI(MPI_Comm_size(comm, &size));
2117:   PetscCall(PetscNew(&atb));
2118:   cM = C->rmap->N;
2119:   PetscCall(PetscMalloc2(cM * cN, &atb->sendbuf, size, &atb->recvcounts));
2120:   PetscCall(MatGetOwnershipRanges(C, &ranges));
2121:   for (i = 0; i < size; i++) atb->recvcounts[i] = (ranges[i + 1] - ranges[i]) * cN;

2123:   C->product->data    = atb;
2124:   C->product->destroy = MatDestroy_MatTransMatMult_MPIDense_MPIDense;
2125:   PetscFunctionReturn(PETSC_SUCCESS);
2126: }

2128: static PetscErrorCode MatMatTransposeMultSymbolic_MPIDense_MPIDense(Mat A, Mat B, PetscReal fill, Mat C)
2129: {
2130:   MPI_Comm               comm;
2131:   PetscMPIInt            i, size;
2132:   PetscInt               maxRows, bufsiz;
2133:   PetscMPIInt            tag;
2134:   PetscInt               alg;
2135:   Mat_MatTransMultDense *abt;
2136:   Mat_Product           *product = C->product;
2137:   PetscBool              flg;

2139:   PetscFunctionBegin;
2140:   MatCheckProduct(C, 4);
2141:   PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
2142:   /* check local size of A and B */
2143:   PetscCheck(A->cmap->n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Matrix local column dimensions are incompatible, A (%" PetscInt_FMT ") != B (%" PetscInt_FMT ")", A->cmap->n, B->cmap->n);

2145:   PetscCall(PetscStrcmp(product->alg, "allgatherv", &flg));
2146:   alg = flg ? 0 : 1;

2148:   /* setup matrix product C */
2149:   PetscCall(MatSetSizes(C, A->rmap->n, B->rmap->n, A->rmap->N, B->rmap->N));
2150:   PetscCall(MatSetType(C, MATMPIDENSE));
2151:   PetscCall(MatSetUp(C));
2152:   PetscCall(PetscObjectGetNewTag((PetscObject)C, &tag));

2154:   /* create data structure for reuse C */
2155:   PetscCall(PetscObjectGetComm((PetscObject)C, &comm));
2156:   PetscCallMPI(MPI_Comm_size(comm, &size));
2157:   PetscCall(PetscNew(&abt));
2158:   abt->tag = tag;
2159:   abt->alg = alg;
2160:   switch (alg) {
2161:   case 1: /* alg: "cyclic" */
2162:     for (maxRows = 0, i = 0; i < size; i++) maxRows = PetscMax(maxRows, (B->rmap->range[i + 1] - B->rmap->range[i]));
2163:     bufsiz = A->cmap->N * maxRows;
2164:     PetscCall(PetscMalloc2(bufsiz, &abt->buf[0], bufsiz, &abt->buf[1]));
2165:     break;
2166:   default: /* alg: "allgatherv" */
2167:     PetscCall(PetscMalloc2(B->rmap->n * B->cmap->N, &abt->buf[0], B->rmap->N * B->cmap->N, &abt->buf[1]));
2168:     PetscCall(PetscMalloc2(size, &abt->recvcounts, size + 1, &abt->recvdispls));
2169:     for (i = 0; i <= size; i++) abt->recvdispls[i] = B->rmap->range[i] * A->cmap->N;
2170:     for (i = 0; i < size; i++) abt->recvcounts[i] = abt->recvdispls[i + 1] - abt->recvdispls[i];
2171:     break;
2172:   }

2174:   C->product->data                = abt;
2175:   C->product->destroy             = MatDestroy_MatMatTransMult_MPIDense_MPIDense;
2176:   C->ops->mattransposemultnumeric = MatMatTransposeMultNumeric_MPIDense_MPIDense;
2177:   PetscFunctionReturn(PETSC_SUCCESS);
2178: }

2180: static PetscErrorCode MatMatTransposeMultNumeric_MPIDense_MPIDense_Cyclic(Mat A, Mat B, Mat C)
2181: {
2182:   Mat_MPIDense          *a = (Mat_MPIDense *)A->data, *b = (Mat_MPIDense *)B->data, *c = (Mat_MPIDense *)C->data;
2183:   Mat_MatTransMultDense *abt;
2184:   MPI_Comm               comm;
2185:   PetscMPIInt            rank, size, sendsiz, recvsiz, sendto, recvfrom, recvisfrom;
2186:   PetscScalar           *sendbuf, *recvbuf = NULL, *cv;
2187:   PetscInt               i, cK             = A->cmap->N, k, j, bn;
2188:   PetscScalar            _DOne = 1.0, _DZero = 0.0;
2189:   const PetscScalar     *av, *bv;
2190:   PetscBLASInt           cm, cn, ck, alda, blda = 0, clda;
2191:   MPI_Request            reqs[2];
2192:   const PetscInt        *ranges;

2194:   PetscFunctionBegin;
2195:   MatCheckProduct(C, 3);
2196:   PetscCheck(C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty");
2197:   abt = (Mat_MatTransMultDense *)C->product->data;
2198:   PetscCall(PetscObjectGetComm((PetscObject)C, &comm));
2199:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
2200:   PetscCallMPI(MPI_Comm_size(comm, &size));
2201:   PetscCall(MatDenseGetArrayRead(a->A, &av));
2202:   PetscCall(MatDenseGetArrayRead(b->A, &bv));
2203:   PetscCall(MatDenseGetArrayWrite(c->A, &cv));
2204:   PetscCall(MatDenseGetLDA(a->A, &i));
2205:   PetscCall(PetscBLASIntCast(i, &alda));
2206:   PetscCall(MatDenseGetLDA(b->A, &i));
2207:   PetscCall(PetscBLASIntCast(i, &blda));
2208:   PetscCall(MatDenseGetLDA(c->A, &i));
2209:   PetscCall(PetscBLASIntCast(i, &clda));
2210:   PetscCall(MatGetOwnershipRanges(B, &ranges));
2211:   bn = B->rmap->n;
2212:   if (blda == bn) {
2213:     sendbuf = (PetscScalar *)bv;
2214:   } else {
2215:     sendbuf = abt->buf[0];
2216:     for (k = 0, i = 0; i < cK; i++) {
2217:       for (j = 0; j < bn; j++, k++) sendbuf[k] = bv[i * blda + j];
2218:     }
2219:   }
2220:   if (size > 1) {
2221:     sendto   = (rank + size - 1) % size;
2222:     recvfrom = (rank + size + 1) % size;
2223:   } else {
2224:     sendto = recvfrom = 0;
2225:   }
2226:   PetscCall(PetscBLASIntCast(cK, &ck));
2227:   PetscCall(PetscBLASIntCast(c->A->rmap->n, &cm));
2228:   recvisfrom = rank;
2229:   for (i = 0; i < size; i++) {
2230:     /* we have finished receiving in sending, bufs can be read/modified */
2231:     PetscInt nextrecvisfrom = (recvisfrom + 1) % size; /* which process the next recvbuf will originate on */
2232:     PetscInt nextbn         = ranges[nextrecvisfrom + 1] - ranges[nextrecvisfrom];

2234:     if (nextrecvisfrom != rank) {
2235:       /* start the cyclic sends from sendbuf, to recvbuf (which will switch to sendbuf) */
2236:       sendsiz = cK * bn;
2237:       recvsiz = cK * nextbn;
2238:       recvbuf = (i & 1) ? abt->buf[0] : abt->buf[1];
2239:       PetscCallMPI(MPI_Isend(sendbuf, sendsiz, MPIU_SCALAR, sendto, abt->tag, comm, &reqs[0]));
2240:       PetscCallMPI(MPI_Irecv(recvbuf, recvsiz, MPIU_SCALAR, recvfrom, abt->tag, comm, &reqs[1]));
2241:     }

2243:     /* local aseq * sendbuf^T */
2244:     PetscCall(PetscBLASIntCast(ranges[recvisfrom + 1] - ranges[recvisfrom], &cn));
2245:     if (cm && cn && ck) PetscCallBLAS("BLASgemm", BLASgemm_("N", "T", &cm, &cn, &ck, &_DOne, av, &alda, sendbuf, &cn, &_DZero, cv + clda * ranges[recvisfrom], &clda));

2247:     if (nextrecvisfrom != rank) {
2248:       /* wait for the sends and receives to complete, swap sendbuf and recvbuf */
2249:       PetscCallMPI(MPI_Waitall(2, reqs, MPI_STATUSES_IGNORE));
2250:     }
2251:     bn         = nextbn;
2252:     recvisfrom = nextrecvisfrom;
2253:     sendbuf    = recvbuf;
2254:   }
2255:   PetscCall(MatDenseRestoreArrayRead(a->A, &av));
2256:   PetscCall(MatDenseRestoreArrayRead(b->A, &bv));
2257:   PetscCall(MatDenseRestoreArrayWrite(c->A, &cv));
2258:   PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
2259:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
2260:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
2261:   PetscFunctionReturn(PETSC_SUCCESS);
2262: }

2264: static PetscErrorCode MatMatTransposeMultNumeric_MPIDense_MPIDense_Allgatherv(Mat A, Mat B, Mat C)
2265: {
2266:   Mat_MPIDense          *a = (Mat_MPIDense *)A->data, *b = (Mat_MPIDense *)B->data, *c = (Mat_MPIDense *)C->data;
2267:   Mat_MatTransMultDense *abt;
2268:   MPI_Comm               comm;
2269:   PetscMPIInt            size;
2270:   PetscScalar           *cv, *sendbuf, *recvbuf;
2271:   const PetscScalar     *av, *bv;
2272:   PetscInt               blda, i, cK = A->cmap->N, k, j, bn;
2273:   PetscScalar            _DOne = 1.0, _DZero = 0.0;
2274:   PetscBLASInt           cm, cn, ck, alda, clda;

2276:   PetscFunctionBegin;
2277:   MatCheckProduct(C, 3);
2278:   PetscCheck(C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty");
2279:   abt = (Mat_MatTransMultDense *)C->product->data;
2280:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
2281:   PetscCallMPI(MPI_Comm_size(comm, &size));
2282:   PetscCall(MatDenseGetArrayRead(a->A, &av));
2283:   PetscCall(MatDenseGetArrayRead(b->A, &bv));
2284:   PetscCall(MatDenseGetArrayWrite(c->A, &cv));
2285:   PetscCall(MatDenseGetLDA(a->A, &i));
2286:   PetscCall(PetscBLASIntCast(i, &alda));
2287:   PetscCall(MatDenseGetLDA(b->A, &blda));
2288:   PetscCall(MatDenseGetLDA(c->A, &i));
2289:   PetscCall(PetscBLASIntCast(i, &clda));
2290:   /* copy transpose of B into buf[0] */
2291:   bn      = B->rmap->n;
2292:   sendbuf = abt->buf[0];
2293:   recvbuf = abt->buf[1];
2294:   for (k = 0, j = 0; j < bn; j++) {
2295:     for (i = 0; i < cK; i++, k++) sendbuf[k] = bv[i * blda + j];
2296:   }
2297:   PetscCall(MatDenseRestoreArrayRead(b->A, &bv));
2298:   PetscCallMPI(MPI_Allgatherv(sendbuf, bn * cK, MPIU_SCALAR, recvbuf, abt->recvcounts, abt->recvdispls, MPIU_SCALAR, comm));
2299:   PetscCall(PetscBLASIntCast(cK, &ck));
2300:   PetscCall(PetscBLASIntCast(c->A->rmap->n, &cm));
2301:   PetscCall(PetscBLASIntCast(c->A->cmap->n, &cn));
2302:   if (cm && cn && ck) PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &cm, &cn, &ck, &_DOne, av, &alda, recvbuf, &ck, &_DZero, cv, &clda));
2303:   PetscCall(MatDenseRestoreArrayRead(a->A, &av));
2304:   PetscCall(MatDenseRestoreArrayRead(b->A, &bv));
2305:   PetscCall(MatDenseRestoreArrayWrite(c->A, &cv));
2306:   PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
2307:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
2308:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
2309:   PetscFunctionReturn(PETSC_SUCCESS);
2310: }

2312: static PetscErrorCode MatMatTransposeMultNumeric_MPIDense_MPIDense(Mat A, Mat B, Mat C)
2313: {
2314:   Mat_MatTransMultDense *abt;

2316:   PetscFunctionBegin;
2317:   MatCheckProduct(C, 3);
2318:   PetscCheck(C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty");
2319:   abt = (Mat_MatTransMultDense *)C->product->data;
2320:   switch (abt->alg) {
2321:   case 1:
2322:     PetscCall(MatMatTransposeMultNumeric_MPIDense_MPIDense_Cyclic(A, B, C));
2323:     break;
2324:   default:
2325:     PetscCall(MatMatTransposeMultNumeric_MPIDense_MPIDense_Allgatherv(A, B, C));
2326:     break;
2327:   }
2328:   PetscFunctionReturn(PETSC_SUCCESS);
2329: }

2331: static PetscErrorCode MatDestroy_MatMatMult_MPIDense_MPIDense(void *data)
2332: {
2333:   Mat_MatMultDense *ab = (Mat_MatMultDense *)data;

2335:   PetscFunctionBegin;
2336:   PetscCall(MatDestroy(&ab->Ce));
2337:   PetscCall(MatDestroy(&ab->Ae));
2338:   PetscCall(MatDestroy(&ab->Be));
2339:   PetscCall(PetscFree(ab));
2340:   PetscFunctionReturn(PETSC_SUCCESS);
2341: }

2343: static PetscErrorCode MatMatMultNumeric_MPIDense_MPIDense(Mat A, Mat B, Mat C)
2344: {
2345:   Mat_MatMultDense *ab;
2346:   Mat_MPIDense     *mdn = (Mat_MPIDense *)A->data;

2348:   PetscFunctionBegin;
2349:   MatCheckProduct(C, 3);
2350:   PetscCheck(C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Missing product data");
2351:   ab = (Mat_MatMultDense *)C->product->data;
2352:   if (ab->Ae && ab->Ce) {
2353: #if PetscDefined(HAVE_ELEMENTAL)
2354:     PetscCall(MatConvert_MPIDense_Elemental(A, MATELEMENTAL, MAT_REUSE_MATRIX, &ab->Ae));
2355:     PetscCall(MatConvert_MPIDense_Elemental(B, MATELEMENTAL, MAT_REUSE_MATRIX, &ab->Be));
2356:     PetscCall(MatMatMultNumeric_Elemental(ab->Ae, ab->Be, ab->Ce));
2357:     PetscCall(MatConvert(ab->Ce, MATMPIDENSE, MAT_REUSE_MATRIX, &C));
2358: #else
2359:     SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "PETSC_HAVE_ELEMENTAL not defined");
2360: #endif
2361:   } else {
2362:     const PetscScalar *read;
2363:     PetscScalar       *write;
2364:     PetscInt           lda;

2366:     PetscCall(MatDenseGetLDA(B, &lda));
2367:     PetscCall(MatDenseGetArrayRead(B, &read));
2368:     PetscCall(MatDenseGetArrayWrite(ab->Be, &write));
2369:     if (!mdn->Mvctx) PetscCall(MatSetUpMultiply_MPIDense(A)); /* cannot be done during the symbolic phase because of possible calls to MatProductReplaceMats() */
2370:     for (PetscInt i = 0; i < C->cmap->N; ++i) {
2371:       PetscCall(PetscSFBcastBegin(mdn->Mvctx, MPIU_SCALAR, read + i * lda, write + i * ab->Be->rmap->n, MPI_REPLACE));
2372:       PetscCall(PetscSFBcastEnd(mdn->Mvctx, MPIU_SCALAR, read + i * lda, write + i * ab->Be->rmap->n, MPI_REPLACE));
2373:     }
2374:     PetscCall(MatDenseRestoreArrayWrite(ab->Be, &write));
2375:     PetscCall(MatDenseRestoreArrayRead(B, &read));
2376:     PetscCall(MatMatMultNumeric_SeqDense_SeqDense(((Mat_MPIDense *)A->data)->A, ab->Be, ((Mat_MPIDense *)C->data)->A));
2377:   }
2378:   PetscFunctionReturn(PETSC_SUCCESS);
2379: }

2381: static PetscErrorCode MatMatMultSymbolic_MPIDense_MPIDense(Mat A, Mat B, PetscReal fill, Mat C)
2382: {
2383:   Mat_Product      *product = C->product;
2384:   PetscInt          alg;
2385:   Mat_MatMultDense *ab;
2386:   PetscBool         flg;

2388:   PetscFunctionBegin;
2389:   MatCheckProduct(C, 4);
2390:   PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
2391:   /* check local size of A and B */
2392:   PetscCheck(A->cmap->rstart == B->rmap->rstart && A->cmap->rend == B->rmap->rend, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_SIZ, "Matrix local dimensions are incompatible, A (%" PetscInt_FMT ", %" PetscInt_FMT ") != B (%" PetscInt_FMT ", %" PetscInt_FMT ")",
2393:              A->rmap->rstart, A->rmap->rend, B->rmap->rstart, B->rmap->rend);

2395:   PetscCall(PetscStrcmp(product->alg, "petsc", &flg));
2396:   alg = flg ? 0 : 1;

2398:   /* setup C */
2399:   PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
2400:   PetscCall(MatSetType(C, MATMPIDENSE));
2401:   PetscCall(MatSetUp(C));

2403:   /* create data structure for reuse Cdense */
2404:   PetscCall(PetscNew(&ab));

2406:   switch (alg) {
2407:   case 1: /* alg: "elemental" */
2408: #if PetscDefined(HAVE_ELEMENTAL)
2409:     /* create elemental matrices Ae and Be */
2410:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &ab->Ae));
2411:     PetscCall(MatSetSizes(ab->Ae, PETSC_DECIDE, PETSC_DECIDE, A->rmap->N, A->cmap->N));
2412:     PetscCall(MatSetType(ab->Ae, MATELEMENTAL));
2413:     PetscCall(MatSetUp(ab->Ae));
2414:     PetscCall(MatSetOption(ab->Ae, MAT_ROW_ORIENTED, PETSC_FALSE));

2416:     PetscCall(MatCreate(PetscObjectComm((PetscObject)B), &ab->Be));
2417:     PetscCall(MatSetSizes(ab->Be, PETSC_DECIDE, PETSC_DECIDE, B->rmap->N, B->cmap->N));
2418:     PetscCall(MatSetType(ab->Be, MATELEMENTAL));
2419:     PetscCall(MatSetUp(ab->Be));
2420:     PetscCall(MatSetOption(ab->Be, MAT_ROW_ORIENTED, PETSC_FALSE));

2422:     /* compute symbolic Ce = Ae*Be */
2423:     PetscCall(MatCreate(PetscObjectComm((PetscObject)C), &ab->Ce));
2424:     PetscCall(MatMatMultSymbolic_Elemental(ab->Ae, ab->Be, fill, ab->Ce));
2425: #else
2426:     SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "PETSC_HAVE_ELEMENTAL not defined");
2427: #endif
2428:     break;
2429:   default: /* alg: "petsc" */
2430:     ab->Ae = NULL;
2431:     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, A->cmap->N, B->cmap->N, NULL, &ab->Be));
2432:     ab->Ce = NULL;
2433:     break;
2434:   }

2436:   C->product->data       = ab;
2437:   C->product->destroy    = MatDestroy_MatMatMult_MPIDense_MPIDense;
2438:   C->ops->matmultnumeric = MatMatMultNumeric_MPIDense_MPIDense;
2439:   PetscFunctionReturn(PETSC_SUCCESS);
2440: }

2442: static PetscErrorCode MatProductSetFromOptions_MPIDense_AB(Mat C)
2443: {
2444:   Mat_Product *product     = C->product;
2445:   const char  *algTypes[2] = {"petsc", "elemental"};
2446:   PetscInt     alg, nalg = PetscDefined(HAVE_ELEMENTAL) ? 2 : 1;
2447:   PetscBool    flg = PETSC_FALSE;

2449:   PetscFunctionBegin;
2450:   /* Set default algorithm */
2451:   alg = 0; /* default is petsc */
2452:   PetscCall(PetscStrcmp(product->alg, "default", &flg));
2453:   if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));

2455:   /* Get runtime option */
2456:   PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_AB", "Mat");
2457:   PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_AB", algTypes, nalg, algTypes[alg], &alg, &flg));
2458:   PetscOptionsEnd();
2459:   if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));

2461:   C->ops->matmultsymbolic = MatMatMultSymbolic_MPIDense_MPIDense;
2462:   C->ops->productsymbolic = MatProductSymbolic_AB;
2463:   PetscFunctionReturn(PETSC_SUCCESS);
2464: }

2466: static PetscErrorCode MatProductSetFromOptions_MPIDense_AtB(Mat C)
2467: {
2468:   Mat_Product *product = C->product;
2469:   Mat          A = product->A, B = product->B;

2471:   PetscFunctionBegin;
2472:   PetscCheck(A->rmap->rstart == B->rmap->rstart && A->rmap->rend == B->rmap->rend, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Matrix local dimensions are incompatible, (%" PetscInt_FMT ", %" PetscInt_FMT ") != (%" PetscInt_FMT ",%" PetscInt_FMT ")",
2473:              A->rmap->rstart, A->rmap->rend, B->rmap->rstart, B->rmap->rend);
2474:   C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_MPIDense_MPIDense;
2475:   C->ops->productsymbolic          = MatProductSymbolic_AtB;
2476:   PetscFunctionReturn(PETSC_SUCCESS);
2477: }

2479: static PetscErrorCode MatProductSetFromOptions_MPIDense_ABt(Mat C)
2480: {
2481:   Mat_Product *product     = C->product;
2482:   const char  *algTypes[2] = {"allgatherv", "cyclic"};
2483:   PetscInt     alg, nalg = 2;
2484:   PetscBool    flg = PETSC_FALSE;

2486:   PetscFunctionBegin;
2487:   /* Set default algorithm */
2488:   alg = 0; /* default is allgatherv */
2489:   PetscCall(PetscStrcmp(product->alg, "default", &flg));
2490:   if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));

2492:   /* Get runtime option */
2493:   if (product->api_user) {
2494:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatMatTransposeMult", "Mat");
2495:     PetscCall(PetscOptionsEList("-matmattransmult_mpidense_mpidense_via", "Algorithmic approach", "MatMatTransposeMult", algTypes, nalg, algTypes[alg], &alg, &flg));
2496:     PetscOptionsEnd();
2497:   } else {
2498:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_ABt", "Mat");
2499:     PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_ABt", algTypes, nalg, algTypes[alg], &alg, &flg));
2500:     PetscOptionsEnd();
2501:   }
2502:   if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));

2504:   C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_MPIDense_MPIDense;
2505:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
2506:   PetscFunctionReturn(PETSC_SUCCESS);
2507: }

2509: static PetscErrorCode MatProductSetFromOptions_MPIDense(Mat C)
2510: {
2511:   Mat_Product *product = C->product;

2513:   PetscFunctionBegin;
2514:   switch (product->type) {
2515:   case MATPRODUCT_AB:
2516:     PetscCall(MatProductSetFromOptions_MPIDense_AB(C));
2517:     break;
2518:   case MATPRODUCT_AtB:
2519:     PetscCall(MatProductSetFromOptions_MPIDense_AtB(C));
2520:     break;
2521:   case MATPRODUCT_ABt:
2522:     PetscCall(MatProductSetFromOptions_MPIDense_ABt(C));
2523:     break;
2524:   default:
2525:     break;
2526:   }
2527:   PetscFunctionReturn(PETSC_SUCCESS);
2528: }