Actual source code: dense.c

  1: /*
  2:      Defines the basic matrix operations for sequential dense.
  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/seq/dense.h>
  8: #include <../src/mat/impls/dense/mpi/mpidense.h>
  9: #include <petscblaslapack.h>
 10: #include <../src/mat/impls/aij/seq/aij.h>
 11: #include <petsc/private/vecimpl.h>

 13: PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian)
 14: {
 15:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
 16:   PetscInt      j, k, n = A->rmap->n;
 17:   PetscScalar  *v;

 19:   PetscFunctionBegin;
 20:   PetscCheck(A->rmap->n == A->cmap->n, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot symmetrize a rectangular matrix");
 21:   PetscCall(MatDenseGetArray(A, &v));
 22:   if (!hermitian) {
 23:     for (k = 0; k < n; k++) {
 24:       for (j = k; j < n; j++) v[j * mat->lda + k] = v[k * mat->lda + j];
 25:     }
 26:   } else {
 27:     for (k = 0; k < n; k++) {
 28:       for (j = k; j < n; j++) v[j * mat->lda + k] = PetscConj(v[k * mat->lda + j]);
 29:     }
 30:   }
 31:   PetscCall(MatDenseRestoreArray(A, &v));
 32:   PetscFunctionReturn(PETSC_SUCCESS);
 33: }

 35: PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A)
 36: {
 37:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
 38:   PetscBLASInt  info, n;

 40:   PetscFunctionBegin;
 41:   if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
 42:   PetscCall(PetscBLASIntCast(A->cmap->n, &n));
 43:   if (A->factortype == MAT_FACTOR_LU) {
 44:     PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present");
 45:     if (!mat->fwork) {
 46:       mat->lfwork = n;
 47:       PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
 48:     }
 49:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
 50:     PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
 51:     PetscCall(PetscFPTrapPop());
 52:     PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
 53:   } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
 54:     if (A->spd == PETSC_BOOL3_TRUE) {
 55:       PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
 56:       PetscCallBLAS("LAPACKpotri", LAPACKpotri_("L", &n, mat->v, &mat->lda, &info));
 57:       PetscCall(PetscFPTrapPop());
 58:       PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_TRUE));
 59: #if defined(PETSC_USE_COMPLEX)
 60:     } else if (A->hermitian == PETSC_BOOL3_TRUE) {
 61:       PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present");
 62:       PetscCheck(mat->fwork, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Fwork not present");
 63:       PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
 64:       PetscCallBLAS("LAPACKhetri", LAPACKhetri_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &info));
 65:       PetscCall(PetscFPTrapPop());
 66:       PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_TRUE));
 67: #endif
 68:     } else { /* symmetric case */
 69:       PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present");
 70:       PetscCheck(mat->fwork, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Fwork not present");
 71:       PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
 72:       PetscCallBLAS("LAPACKsytri", LAPACKsytri_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &info));
 73:       PetscCall(PetscFPTrapPop());
 74:       PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_FALSE));
 75:     }
 76:     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad Inversion: zero pivot in row %" PetscInt_FMT, (PetscInt)info - 1);
 77:     PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
 78:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix must be factored to solve");

 80:   A->ops->solve             = NULL;
 81:   A->ops->matsolve          = NULL;
 82:   A->ops->solvetranspose    = NULL;
 83:   A->ops->matsolvetranspose = NULL;
 84:   A->ops->solveadd          = NULL;
 85:   A->ops->solvetransposeadd = NULL;
 86:   A->factortype             = MAT_FACTOR_NONE;
 87:   PetscCall(PetscFree(A->solvertype));
 88:   PetscFunctionReturn(PETSC_SUCCESS);
 89: }

 91: static PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
 92: {
 93:   Mat_SeqDense      *l = (Mat_SeqDense *)A->data;
 94:   PetscInt           m = l->lda, n = A->cmap->n, r = A->rmap->n, i, j;
 95:   PetscScalar       *slot, *bb, *v;
 96:   const PetscScalar *xx;

 98:   PetscFunctionBegin;
 99:   if (PetscDefined(USE_DEBUG)) {
100:     for (i = 0; i < N; i++) {
101:       PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
102:       PetscCheck(rows[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " requested to be zeroed greater than or equal number of rows %" PetscInt_FMT, rows[i], A->rmap->n);
103:       PetscCheck(rows[i] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Col %" PetscInt_FMT " requested to be zeroed greater than or equal number of cols %" PetscInt_FMT, rows[i], A->cmap->n);
104:     }
105:   }
106:   if (!N) PetscFunctionReturn(PETSC_SUCCESS);

108:   /* fix right-hand side if needed */
109:   if (x && b) {
110:     Vec xt;

112:     PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
113:     PetscCall(VecDuplicate(x, &xt));
114:     PetscCall(VecCopy(x, xt));
115:     PetscCall(VecScale(xt, -1.0));
116:     PetscCall(MatMultAdd(A, xt, b, b));
117:     PetscCall(VecDestroy(&xt));
118:     PetscCall(VecGetArrayRead(x, &xx));
119:     PetscCall(VecGetArray(b, &bb));
120:     for (i = 0; i < N; i++) bb[rows[i]] = diag * xx[rows[i]];
121:     PetscCall(VecRestoreArrayRead(x, &xx));
122:     PetscCall(VecRestoreArray(b, &bb));
123:   }

125:   PetscCall(MatDenseGetArray(A, &v));
126:   for (i = 0; i < N; i++) {
127:     slot = v + rows[i] * m;
128:     PetscCall(PetscArrayzero(slot, r));
129:   }
130:   for (i = 0; i < N; i++) {
131:     slot = v + rows[i];
132:     for (j = 0; j < n; j++) {
133:       *slot = 0.0;
134:       slot += m;
135:     }
136:   }
137:   if (diag != 0.0) {
138:     PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
139:     for (i = 0; i < N; i++) {
140:       slot  = v + (m + 1) * rows[i];
141:       *slot = diag;
142:     }
143:   }
144:   PetscCall(MatDenseRestoreArray(A, &v));
145:   PetscFunctionReturn(PETSC_SUCCESS);
146: }

148: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
149: {
150:   Mat              B = NULL;
151:   Mat_SeqAIJ      *a = (Mat_SeqAIJ *)A->data;
152:   Mat_SeqDense    *b;
153:   PetscInt        *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->N, i;
154:   const MatScalar *av;
155:   PetscBool        isseqdense;

157:   PetscFunctionBegin;
158:   if (reuse == MAT_REUSE_MATRIX) {
159:     PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATSEQDENSE, &isseqdense));
160:     PetscCheck(isseqdense, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_USER, "Cannot reuse matrix of type %s", ((PetscObject)*newmat)->type_name);
161:   }
162:   if (reuse != MAT_REUSE_MATRIX) {
163:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
164:     PetscCall(MatSetSizes(B, m, n, m, n));
165:     PetscCall(MatSetType(B, MATSEQDENSE));
166:     PetscCall(MatSeqDenseSetPreallocation(B, NULL));
167:     b = (Mat_SeqDense *)B->data;
168:   } else {
169:     b = (Mat_SeqDense *)((*newmat)->data);
170:     for (i = 0; i < n; i++) PetscCall(PetscArrayzero(b->v + i * b->lda, m));
171:   }
172:   PetscCall(MatSeqAIJGetArrayRead(A, &av));
173:   for (i = 0; i < m; i++) {
174:     PetscInt j;
175:     for (j = 0; j < ai[1] - ai[0]; j++) {
176:       b->v[*aj * b->lda + i] = *av;
177:       aj++;
178:       av++;
179:     }
180:     ai++;
181:   }
182:   PetscCall(MatSeqAIJRestoreArrayRead(A, &av));

184:   if (reuse == MAT_INPLACE_MATRIX) {
185:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
186:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
187:     PetscCall(MatHeaderReplace(A, &B));
188:   } else {
189:     if (B) *newmat = B;
190:     PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
191:     PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
192:   }
193:   PetscFunctionReturn(PETSC_SUCCESS);
194: }

196: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
197: {
198:   Mat           B = NULL;
199:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;
200:   PetscInt      i, j;
201:   PetscInt     *rows, *nnz;
202:   MatScalar    *aa = a->v, *vals;

204:   PetscFunctionBegin;
205:   PetscCall(PetscCalloc3(A->rmap->n, &rows, A->rmap->n, &nnz, A->rmap->n, &vals));
206:   if (reuse != MAT_REUSE_MATRIX) {
207:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
208:     PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
209:     PetscCall(MatSetType(B, MATSEQAIJ));
210:     for (j = 0; j < A->cmap->n; j++) {
211:       for (i = 0; i < A->rmap->n; i++)
212:         if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) ++nnz[i];
213:       aa += a->lda;
214:     }
215:     PetscCall(MatSeqAIJSetPreallocation(B, PETSC_DETERMINE, nnz));
216:   } else B = *newmat;
217:   aa = a->v;
218:   for (j = 0; j < A->cmap->n; j++) {
219:     PetscInt numRows = 0;
220:     for (i = 0; i < A->rmap->n; i++)
221:       if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) {
222:         rows[numRows]   = i;
223:         vals[numRows++] = aa[i];
224:       }
225:     PetscCall(MatSetValues(B, numRows, rows, 1, &j, vals, INSERT_VALUES));
226:     aa += a->lda;
227:   }
228:   PetscCall(PetscFree3(rows, nnz, vals));
229:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
230:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

232:   if (reuse == MAT_INPLACE_MATRIX) {
233:     PetscCall(MatHeaderReplace(A, &B));
234:   } else if (reuse != MAT_REUSE_MATRIX) *newmat = B;
235:   PetscFunctionReturn(PETSC_SUCCESS);
236: }

238: PetscErrorCode MatAXPY_SeqDense(Mat Y, PetscScalar alpha, Mat X, MatStructure str)
239: {
240:   Mat_SeqDense      *x = (Mat_SeqDense *)X->data, *y = (Mat_SeqDense *)Y->data;
241:   const PetscScalar *xv;
242:   PetscScalar       *yv;
243:   PetscBLASInt       N, m, ldax = 0, lday = 0, one = 1;

245:   PetscFunctionBegin;
246:   PetscCall(MatDenseGetArrayRead(X, &xv));
247:   PetscCall(MatDenseGetArray(Y, &yv));
248:   PetscCall(PetscBLASIntCast(X->rmap->n * X->cmap->n, &N));
249:   PetscCall(PetscBLASIntCast(X->rmap->n, &m));
250:   PetscCall(PetscBLASIntCast(x->lda, &ldax));
251:   PetscCall(PetscBLASIntCast(y->lda, &lday));
252:   if (ldax > m || lday > m) {
253:     for (PetscInt j = 0; j < X->cmap->n; j++) PetscCallBLAS("BLASaxpy", BLASaxpy_(&m, &alpha, PetscSafePointerPlusOffset(xv, j * ldax), &one, PetscSafePointerPlusOffset(yv, j * lday), &one));
254:   } else {
255:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&N, &alpha, xv, &one, yv, &one));
256:   }
257:   PetscCall(MatDenseRestoreArrayRead(X, &xv));
258:   PetscCall(MatDenseRestoreArray(Y, &yv));
259:   PetscCall(PetscLogFlops(PetscMax(2.0 * N - 1, 0)));
260:   PetscFunctionReturn(PETSC_SUCCESS);
261: }

263: static PetscErrorCode MatGetInfo_SeqDense(Mat A, MatInfoType flag, MatInfo *info)
264: {
265:   PetscLogDouble N = A->rmap->n * A->cmap->n;

267:   PetscFunctionBegin;
268:   info->block_size        = 1.0;
269:   info->nz_allocated      = N;
270:   info->nz_used           = N;
271:   info->nz_unneeded       = 0;
272:   info->assemblies        = A->num_ass;
273:   info->mallocs           = 0;
274:   info->memory            = 0; /* REVIEW ME */
275:   info->fill_ratio_given  = 0;
276:   info->fill_ratio_needed = 0;
277:   info->factor_mallocs    = 0;
278:   PetscFunctionReturn(PETSC_SUCCESS);
279: }

281: PetscErrorCode MatScale_SeqDense(Mat A, PetscScalar alpha)
282: {
283:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;
284:   PetscScalar  *v;
285:   PetscBLASInt  one = 1, j, nz, lda = 0;

287:   PetscFunctionBegin;
288:   PetscCall(MatDenseGetArray(A, &v));
289:   PetscCall(PetscBLASIntCast(a->lda, &lda));
290:   if (lda > A->rmap->n) {
291:     PetscCall(PetscBLASIntCast(A->rmap->n, &nz));
292:     for (j = 0; j < A->cmap->n; j++) PetscCallBLAS("BLASscal", BLASscal_(&nz, &alpha, v + j * lda, &one));
293:   } else {
294:     PetscCall(PetscBLASIntCast(A->rmap->n * A->cmap->n, &nz));
295:     PetscCallBLAS("BLASscal", BLASscal_(&nz, &alpha, v, &one));
296:   }
297:   PetscCall(PetscLogFlops(A->rmap->n * A->cmap->n));
298:   PetscCall(MatDenseRestoreArray(A, &v));
299:   PetscFunctionReturn(PETSC_SUCCESS);
300: }

302: PetscErrorCode MatShift_SeqDense(Mat A, PetscScalar alpha)
303: {
304:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;
305:   PetscScalar  *v;
306:   PetscInt      j, k;

308:   PetscFunctionBegin;
309:   PetscCall(MatDenseGetArray(A, &v));
310:   k = PetscMin(A->rmap->n, A->cmap->n);
311:   for (j = 0; j < k; j++) v[j + j * a->lda] += alpha;
312:   PetscCall(PetscLogFlops(k));
313:   PetscCall(MatDenseRestoreArray(A, &v));
314:   PetscFunctionReturn(PETSC_SUCCESS);
315: }

317: static PetscErrorCode MatIsHermitian_SeqDense(Mat A, PetscReal rtol, PetscBool *fl)
318: {
319:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data;
320:   PetscInt           i, j, m = A->rmap->n, N = a->lda;
321:   const PetscScalar *v;

323:   PetscFunctionBegin;
324:   *fl = PETSC_FALSE;
325:   if (A->rmap->n != A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
326:   PetscCall(MatDenseGetArrayRead(A, &v));
327:   for (i = 0; i < m; i++) {
328:     for (j = i; j < m; j++) {
329:       if (PetscAbsScalar(v[i + j * N] - PetscConj(v[j + i * N])) > rtol) goto restore;
330:     }
331:   }
332:   *fl = PETSC_TRUE;
333: restore:
334:   PetscCall(MatDenseRestoreArrayRead(A, &v));
335:   PetscFunctionReturn(PETSC_SUCCESS);
336: }

338: static PetscErrorCode MatIsSymmetric_SeqDense(Mat A, PetscReal rtol, PetscBool *fl)
339: {
340:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data;
341:   PetscInt           i, j, m = A->rmap->n, N = a->lda;
342:   const PetscScalar *v;

344:   PetscFunctionBegin;
345:   *fl = PETSC_FALSE;
346:   if (A->rmap->n != A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
347:   PetscCall(MatDenseGetArrayRead(A, &v));
348:   for (i = 0; i < m; i++) {
349:     for (j = i; j < m; j++) {
350:       if (PetscAbsScalar(v[i + j * N] - v[j + i * N]) > rtol) goto restore;
351:     }
352:   }
353:   *fl = PETSC_TRUE;
354: restore:
355:   PetscCall(MatDenseRestoreArrayRead(A, &v));
356:   PetscFunctionReturn(PETSC_SUCCESS);
357: }

359: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi, Mat A, MatDuplicateOption cpvalues)
360: {
361:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
362:   PetscInt      lda = (PetscInt)mat->lda, j, m, nlda = lda;
363:   PetscBool     isdensecpu;

365:   PetscFunctionBegin;
366:   PetscCall(PetscLayoutReference(A->rmap, &newi->rmap));
367:   PetscCall(PetscLayoutReference(A->cmap, &newi->cmap));
368:   if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */
369:     PetscCall(MatDenseSetLDA(newi, lda));
370:   }
371:   PetscCall(PetscObjectTypeCompare((PetscObject)newi, MATSEQDENSE, &isdensecpu));
372:   if (isdensecpu) PetscCall(MatSeqDenseSetPreallocation(newi, NULL));
373:   if (cpvalues == MAT_COPY_VALUES) {
374:     const PetscScalar *av;
375:     PetscScalar       *v;

377:     PetscCall(MatDenseGetArrayRead(A, &av));
378:     PetscCall(MatDenseGetArrayWrite(newi, &v));
379:     PetscCall(MatDenseGetLDA(newi, &nlda));
380:     m = A->rmap->n;
381:     if (lda > m || nlda > m) {
382:       for (j = 0; j < A->cmap->n; j++) PetscCall(PetscArraycpy(PetscSafePointerPlusOffset(v, j * nlda), PetscSafePointerPlusOffset(av, j * lda), m));
383:     } else {
384:       PetscCall(PetscArraycpy(v, av, A->rmap->n * A->cmap->n));
385:     }
386:     PetscCall(MatDenseRestoreArrayWrite(newi, &v));
387:     PetscCall(MatDenseRestoreArrayRead(A, &av));
388:   }
389:   PetscFunctionReturn(PETSC_SUCCESS);
390: }

392: PetscErrorCode MatDuplicate_SeqDense(Mat A, MatDuplicateOption cpvalues, Mat *newmat)
393: {
394:   PetscFunctionBegin;
395:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), newmat));
396:   PetscCall(MatSetSizes(*newmat, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
397:   PetscCall(MatSetType(*newmat, ((PetscObject)A)->type_name));
398:   PetscCall(MatDuplicateNoCreate_SeqDense(*newmat, A, cpvalues));
399:   PetscFunctionReturn(PETSC_SUCCESS);
400: }

402: static PetscErrorCode MatSolve_SeqDense_Internal_LU(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
403: {
404:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
405:   PetscBLASInt  info;

407:   PetscFunctionBegin;
408:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
409:   PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_(T ? "T" : "N", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
410:   PetscCall(PetscFPTrapPop());
411:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "GETRS - Bad solve %d", (int)info);
412:   PetscCall(PetscLogFlops(nrhs * (2.0 * m * m - m)));
413:   PetscFunctionReturn(PETSC_SUCCESS);
414: }

416: static PetscErrorCode MatConjugate_SeqDense(Mat);

418: static PetscErrorCode MatSolve_SeqDense_Internal_Cholesky(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
419: {
420:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
421:   PetscBLASInt  info;

423:   PetscFunctionBegin;
424:   if (A->spd == PETSC_BOOL3_TRUE) {
425:     if (PetscDefined(USE_COMPLEX) && T) PetscCall(MatConjugate_SeqDense(A));
426:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
427:     PetscCallBLAS("LAPACKpotrs", LAPACKpotrs_("L", &m, &nrhs, mat->v, &mat->lda, x, &m, &info));
428:     PetscCall(PetscFPTrapPop());
429:     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "POTRS Bad solve %d", (int)info);
430:     if (PetscDefined(USE_COMPLEX) && T) PetscCall(MatConjugate_SeqDense(A));
431: #if defined(PETSC_USE_COMPLEX)
432:   } else if (A->hermitian == PETSC_BOOL3_TRUE) {
433:     if (T) PetscCall(MatConjugate_SeqDense(A));
434:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
435:     PetscCallBLAS("LAPACKhetrs", LAPACKhetrs_("L", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
436:     PetscCall(PetscFPTrapPop());
437:     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "HETRS Bad solve %d", (int)info);
438:     if (T) PetscCall(MatConjugate_SeqDense(A));
439: #endif
440:   } else { /* symmetric case */
441:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
442:     PetscCallBLAS("LAPACKsytrs", LAPACKsytrs_("L", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
443:     PetscCall(PetscFPTrapPop());
444:     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "SYTRS Bad solve %d", (int)info);
445:   }
446:   PetscCall(PetscLogFlops(nrhs * (2.0 * m * m - m)));
447:   PetscFunctionReturn(PETSC_SUCCESS);
448: }

450: static PetscErrorCode MatSolve_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
451: {
452:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
453:   PetscBLASInt  info;
454:   char          trans;

456:   PetscFunctionBegin;
457:   if (PetscDefined(USE_COMPLEX)) {
458:     trans = 'C';
459:   } else {
460:     trans = 'T';
461:   }
462:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
463:   { /* lwork depends on the number of right-hand sides */
464:     PetscBLASInt nlfwork, lfwork = -1;
465:     PetscScalar  fwork;

467:     PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", &trans, &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, &fwork, &lfwork, &info));
468:     nlfwork = (PetscBLASInt)PetscRealPart(fwork);
469:     if (nlfwork > mat->lfwork) {
470:       mat->lfwork = nlfwork;
471:       PetscCall(PetscFree(mat->fwork));
472:       PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
473:     }
474:   }
475:   PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", &trans, &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, mat->fwork, &mat->lfwork, &info));
476:   PetscCall(PetscFPTrapPop());
477:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "ORMQR - Bad orthogonal transform %d", (int)info);
478:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
479:   PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "N", "N", &mat->rank, &nrhs, mat->v, &mat->lda, x, &ldx, &info));
480:   PetscCall(PetscFPTrapPop());
481:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "TRTRS - Bad triangular solve %d", (int)info);
482:   for (PetscInt j = 0; j < nrhs; j++) {
483:     for (PetscInt i = mat->rank; i < k; i++) x[j * ldx + i] = 0.;
484:   }
485:   PetscCall(PetscLogFlops(nrhs * (4.0 * m * mat->rank - PetscSqr(mat->rank))));
486:   PetscFunctionReturn(PETSC_SUCCESS);
487: }

489: static PetscErrorCode MatSolveTranspose_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
490: {
491:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
492:   PetscBLASInt  info;

494:   PetscFunctionBegin;
495:   if (A->rmap->n == A->cmap->n && mat->rank == A->rmap->n) {
496:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
497:     PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "T", "N", &m, &nrhs, mat->v, &mat->lda, x, &ldx, &info));
498:     PetscCall(PetscFPTrapPop());
499:     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "TRTRS - Bad triangular solve %d", (int)info);
500:     if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate_SeqDense(A));
501:     { /* lwork depends on the number of right-hand sides */
502:       PetscBLASInt nlfwork, lfwork = -1;
503:       PetscScalar  fwork;

505:       PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", "N", &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, &fwork, &lfwork, &info));
506:       nlfwork = (PetscBLASInt)PetscRealPart(fwork);
507:       if (nlfwork > mat->lfwork) {
508:         mat->lfwork = nlfwork;
509:         PetscCall(PetscFree(mat->fwork));
510:         PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
511:       }
512:     }
513:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
514:     PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", "N", &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, mat->fwork, &mat->lfwork, &info));
515:     PetscCall(PetscFPTrapPop());
516:     PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "ORMQR - Bad orthogonal transform %d", (int)info);
517:     if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate_SeqDense(A));
518:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "QR factored matrix cannot be used for transpose solve");
519:   PetscCall(PetscLogFlops(nrhs * (4.0 * m * mat->rank - PetscSqr(mat->rank))));
520:   PetscFunctionReturn(PETSC_SUCCESS);
521: }

523: static PetscErrorCode MatSolve_SeqDense_SetUp(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
524: {
525:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
526:   PetscScalar  *y;
527:   PetscBLASInt  m = 0, k = 0;

529:   PetscFunctionBegin;
530:   PetscCall(PetscBLASIntCast(A->rmap->n, &m));
531:   PetscCall(PetscBLASIntCast(A->cmap->n, &k));
532:   if (k < m) {
533:     PetscCall(VecCopy(xx, mat->qrrhs));
534:     PetscCall(VecGetArray(mat->qrrhs, &y));
535:   } else {
536:     PetscCall(VecCopy(xx, yy));
537:     PetscCall(VecGetArray(yy, &y));
538:   }
539:   *_y = y;
540:   *_k = k;
541:   *_m = m;
542:   PetscFunctionReturn(PETSC_SUCCESS);
543: }

545: static PetscErrorCode MatSolve_SeqDense_TearDown(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
546: {
547:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
548:   PetscScalar  *y   = NULL;
549:   PetscBLASInt  m, k;

551:   PetscFunctionBegin;
552:   y   = *_y;
553:   *_y = NULL;
554:   k   = *_k;
555:   m   = *_m;
556:   if (k < m) {
557:     PetscScalar *yv;
558:     PetscCall(VecGetArray(yy, &yv));
559:     PetscCall(PetscArraycpy(yv, y, k));
560:     PetscCall(VecRestoreArray(yy, &yv));
561:     PetscCall(VecRestoreArray(mat->qrrhs, &y));
562:   } else {
563:     PetscCall(VecRestoreArray(yy, &y));
564:   }
565:   PetscFunctionReturn(PETSC_SUCCESS);
566: }

568: static PetscErrorCode MatSolve_SeqDense_LU(Mat A, Vec xx, Vec yy)
569: {
570:   PetscScalar *y = NULL;
571:   PetscBLASInt m = 0, k = 0;

573:   PetscFunctionBegin;
574:   PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
575:   PetscCall(MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_FALSE));
576:   PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
577:   PetscFunctionReturn(PETSC_SUCCESS);
578: }

580: static PetscErrorCode MatSolveTranspose_SeqDense_LU(Mat A, Vec xx, Vec yy)
581: {
582:   PetscScalar *y = NULL;
583:   PetscBLASInt m = 0, k = 0;

585:   PetscFunctionBegin;
586:   PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
587:   PetscCall(MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_TRUE));
588:   PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
589:   PetscFunctionReturn(PETSC_SUCCESS);
590: }

592: static PetscErrorCode MatSolve_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
593: {
594:   PetscScalar *y = NULL;
595:   PetscBLASInt m = 0, k = 0;

597:   PetscFunctionBegin;
598:   PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
599:   PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_FALSE));
600:   PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
601:   PetscFunctionReturn(PETSC_SUCCESS);
602: }

604: static PetscErrorCode MatSolveTranspose_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
605: {
606:   PetscScalar *y = NULL;
607:   PetscBLASInt m = 0, k = 0;

609:   PetscFunctionBegin;
610:   PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
611:   PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_TRUE));
612:   PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
613:   PetscFunctionReturn(PETSC_SUCCESS);
614: }

616: static PetscErrorCode MatSolve_SeqDense_QR(Mat A, Vec xx, Vec yy)
617: {
618:   PetscScalar *y = NULL;
619:   PetscBLASInt m = 0, k = 0;

621:   PetscFunctionBegin;
622:   PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
623:   PetscCall(MatSolve_SeqDense_Internal_QR(A, y, PetscMax(m, k), m, 1, k));
624:   PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
625:   PetscFunctionReturn(PETSC_SUCCESS);
626: }

628: static PetscErrorCode MatSolveTranspose_SeqDense_QR(Mat A, Vec xx, Vec yy)
629: {
630:   PetscScalar *y = NULL;
631:   PetscBLASInt m = 0, k = 0;

633:   PetscFunctionBegin;
634:   PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
635:   PetscCall(MatSolveTranspose_SeqDense_Internal_QR(A, y, PetscMax(m, k), m, 1, k));
636:   PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
637:   PetscFunctionReturn(PETSC_SUCCESS);
638: }

640: static PetscErrorCode MatMatSolve_SeqDense_SetUp(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
641: {
642:   const PetscScalar *b;
643:   PetscScalar       *y;
644:   PetscInt           n, _ldb, _ldx;
645:   PetscBLASInt       nrhs = 0, m = 0, k = 0, ldb = 0, ldx = 0, ldy = 0;

647:   PetscFunctionBegin;
648:   *_ldy  = 0;
649:   *_m    = 0;
650:   *_nrhs = 0;
651:   *_k    = 0;
652:   *_y    = NULL;
653:   PetscCall(PetscBLASIntCast(A->rmap->n, &m));
654:   PetscCall(PetscBLASIntCast(A->cmap->n, &k));
655:   PetscCall(MatGetSize(B, NULL, &n));
656:   PetscCall(PetscBLASIntCast(n, &nrhs));
657:   PetscCall(MatDenseGetLDA(B, &_ldb));
658:   PetscCall(PetscBLASIntCast(_ldb, &ldb));
659:   PetscCall(MatDenseGetLDA(X, &_ldx));
660:   PetscCall(PetscBLASIntCast(_ldx, &ldx));
661:   if (ldx < m) {
662:     PetscCall(MatDenseGetArrayRead(B, &b));
663:     PetscCall(PetscMalloc1(nrhs * m, &y));
664:     if (ldb == m) {
665:       PetscCall(PetscArraycpy(y, b, ldb * nrhs));
666:     } else {
667:       for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&y[j * m], &b[j * ldb], m));
668:     }
669:     ldy = m;
670:     PetscCall(MatDenseRestoreArrayRead(B, &b));
671:   } else {
672:     if (ldb == ldx) {
673:       PetscCall(MatCopy(B, X, SAME_NONZERO_PATTERN));
674:       PetscCall(MatDenseGetArray(X, &y));
675:     } else {
676:       PetscCall(MatDenseGetArray(X, &y));
677:       PetscCall(MatDenseGetArrayRead(B, &b));
678:       for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&y[j * ldx], &b[j * ldb], m));
679:       PetscCall(MatDenseRestoreArrayRead(B, &b));
680:     }
681:     ldy = ldx;
682:   }
683:   *_y    = y;
684:   *_ldy  = ldy;
685:   *_k    = k;
686:   *_m    = m;
687:   *_nrhs = nrhs;
688:   PetscFunctionReturn(PETSC_SUCCESS);
689: }

691: static PetscErrorCode MatMatSolve_SeqDense_TearDown(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
692: {
693:   PetscScalar *y;
694:   PetscInt     _ldx;
695:   PetscBLASInt k, ldy, nrhs, ldx = 0;

697:   PetscFunctionBegin;
698:   y    = *_y;
699:   *_y  = NULL;
700:   k    = *_k;
701:   ldy  = *_ldy;
702:   nrhs = *_nrhs;
703:   PetscCall(MatDenseGetLDA(X, &_ldx));
704:   PetscCall(PetscBLASIntCast(_ldx, &ldx));
705:   if (ldx != ldy) {
706:     PetscScalar *xv;
707:     PetscCall(MatDenseGetArray(X, &xv));
708:     for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&xv[j * ldx], &y[j * ldy], k));
709:     PetscCall(MatDenseRestoreArray(X, &xv));
710:     PetscCall(PetscFree(y));
711:   } else {
712:     PetscCall(MatDenseRestoreArray(X, &y));
713:   }
714:   PetscFunctionReturn(PETSC_SUCCESS);
715: }

717: static PetscErrorCode MatMatSolve_SeqDense_LU(Mat A, Mat B, Mat X)
718: {
719:   PetscScalar *y;
720:   PetscBLASInt m, k, ldy, nrhs;

722:   PetscFunctionBegin;
723:   PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
724:   PetscCall(MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_FALSE));
725:   PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
726:   PetscFunctionReturn(PETSC_SUCCESS);
727: }

729: static PetscErrorCode MatMatSolveTranspose_SeqDense_LU(Mat A, Mat B, Mat X)
730: {
731:   PetscScalar *y;
732:   PetscBLASInt m, k, ldy, nrhs;

734:   PetscFunctionBegin;
735:   PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
736:   PetscCall(MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_TRUE));
737:   PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
738:   PetscFunctionReturn(PETSC_SUCCESS);
739: }

741: static PetscErrorCode MatMatSolve_SeqDense_Cholesky(Mat A, Mat B, Mat X)
742: {
743:   PetscScalar *y;
744:   PetscBLASInt m, k, ldy, nrhs;

746:   PetscFunctionBegin;
747:   PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
748:   PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_FALSE));
749:   PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
750:   PetscFunctionReturn(PETSC_SUCCESS);
751: }

753: static PetscErrorCode MatMatSolveTranspose_SeqDense_Cholesky(Mat A, Mat B, Mat X)
754: {
755:   PetscScalar *y;
756:   PetscBLASInt m, k, ldy, nrhs;

758:   PetscFunctionBegin;
759:   PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
760:   PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_TRUE));
761:   PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
762:   PetscFunctionReturn(PETSC_SUCCESS);
763: }

765: static PetscErrorCode MatMatSolve_SeqDense_QR(Mat A, Mat B, Mat X)
766: {
767:   PetscScalar *y;
768:   PetscBLASInt m, k, ldy, nrhs;

770:   PetscFunctionBegin;
771:   PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
772:   PetscCall(MatSolve_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k));
773:   PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
774:   PetscFunctionReturn(PETSC_SUCCESS);
775: }

777: static PetscErrorCode MatMatSolveTranspose_SeqDense_QR(Mat A, Mat B, Mat X)
778: {
779:   PetscScalar *y;
780:   PetscBLASInt m, k, ldy, nrhs;

782:   PetscFunctionBegin;
783:   PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
784:   PetscCall(MatSolveTranspose_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k));
785:   PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
786:   PetscFunctionReturn(PETSC_SUCCESS);
787: }

789: /* COMMENT: I have chosen to hide row permutation in the pivots,
790:    rather than put it in the Mat->row slot.*/
791: PetscErrorCode MatLUFactor_SeqDense(Mat A, IS row, IS col, const MatFactorInfo *minfo)
792: {
793:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
794:   PetscBLASInt  n, m, info;

796:   PetscFunctionBegin;
797:   PetscCall(PetscBLASIntCast(A->cmap->n, &n));
798:   PetscCall(PetscBLASIntCast(A->rmap->n, &m));
799:   if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); }
800:   if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
801:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
802:   PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&m, &n, mat->v, &mat->lda, mat->pivots, &info));
803:   PetscCall(PetscFPTrapPop());

805:   PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to LU factorization %d", (int)info);
806:   PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Bad LU factorization %d", (int)info);

808:   A->ops->solve             = MatSolve_SeqDense_LU;
809:   A->ops->matsolve          = MatMatSolve_SeqDense_LU;
810:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense_LU;
811:   A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_LU;
812:   A->factortype             = MAT_FACTOR_LU;

814:   PetscCall(PetscFree(A->solvertype));
815:   PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));

817:   PetscCall(PetscLogFlops((2.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3));
818:   PetscFunctionReturn(PETSC_SUCCESS);
819: }

821: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info_dummy)
822: {
823:   MatFactorInfo info;

825:   PetscFunctionBegin;
826:   PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
827:   PetscUseTypeMethod(fact, lufactor, NULL, NULL, &info);
828:   PetscFunctionReturn(PETSC_SUCCESS);
829: }

831: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, IS col, const MatFactorInfo *info)
832: {
833:   PetscFunctionBegin;
834:   fact->preallocated         = PETSC_TRUE;
835:   fact->assembled            = PETSC_TRUE;
836:   fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
837:   PetscFunctionReturn(PETSC_SUCCESS);
838: }

840: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
841: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A, IS perm, const MatFactorInfo *factinfo)
842: {
843:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
844:   PetscBLASInt  info, n;

846:   PetscFunctionBegin;
847:   PetscCall(PetscBLASIntCast(A->cmap->n, &n));
848:   if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
849:   if (A->spd == PETSC_BOOL3_TRUE) {
850:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
851:     PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("L", &n, mat->v, &mat->lda, &info));
852:     PetscCall(PetscFPTrapPop());
853: #if defined(PETSC_USE_COMPLEX)
854:   } else if (A->hermitian == PETSC_BOOL3_TRUE) {
855:     if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); }
856:     if (!mat->fwork) {
857:       PetscScalar dummy;

859:       mat->lfwork = -1;
860:       PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
861:       PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info));
862:       PetscCall(PetscFPTrapPop());
863:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
864:       PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
865:     }
866:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
867:     PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
868:     PetscCall(PetscFPTrapPop());
869: #endif
870:   } else { /* symmetric case */
871:     if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); }
872:     if (!mat->fwork) {
873:       PetscScalar dummy;

875:       mat->lfwork = -1;
876:       PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
877:       PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info));
878:       PetscCall(PetscFPTrapPop());
879:       mat->lfwork = (PetscInt)PetscRealPart(dummy);
880:       PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
881:     }
882:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
883:     PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
884:     PetscCall(PetscFPTrapPop());
885:   }
886:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad factorization: zero pivot in row %" PetscInt_FMT, (PetscInt)info - 1);

888:   A->ops->solve             = MatSolve_SeqDense_Cholesky;
889:   A->ops->matsolve          = MatMatSolve_SeqDense_Cholesky;
890:   A->ops->solvetranspose    = MatSolveTranspose_SeqDense_Cholesky;
891:   A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky;
892:   A->factortype             = MAT_FACTOR_CHOLESKY;

894:   PetscCall(PetscFree(A->solvertype));
895:   PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));

897:   PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
898:   PetscFunctionReturn(PETSC_SUCCESS);
899: }

901: static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info_dummy)
902: {
903:   MatFactorInfo info;

905:   PetscFunctionBegin;
906:   info.fill = 1.0;

908:   PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
909:   PetscUseTypeMethod(fact, choleskyfactor, NULL, &info);
910:   PetscFunctionReturn(PETSC_SUCCESS);
911: }

913: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info)
914: {
915:   PetscFunctionBegin;
916:   fact->assembled                  = PETSC_TRUE;
917:   fact->preallocated               = PETSC_TRUE;
918:   fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
919:   PetscFunctionReturn(PETSC_SUCCESS);
920: }

922: PetscErrorCode MatQRFactor_SeqDense(Mat A, IS col, const MatFactorInfo *minfo)
923: {
924:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
925:   PetscBLASInt  n, m, info, min, max;

927:   PetscFunctionBegin;
928:   PetscCall(PetscBLASIntCast(A->cmap->n, &n));
929:   PetscCall(PetscBLASIntCast(A->rmap->n, &m));
930:   max = PetscMax(m, n);
931:   min = PetscMin(m, n);
932:   if (!mat->tau) { PetscCall(PetscMalloc1(min, &mat->tau)); }
933:   if (!mat->pivots) { PetscCall(PetscMalloc1(n, &mat->pivots)); }
934:   if (!mat->qrrhs) PetscCall(MatCreateVecs(A, NULL, &mat->qrrhs));
935:   if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
936:   if (!mat->fwork) {
937:     PetscScalar dummy;

939:     mat->lfwork = -1;
940:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
941:     PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, &dummy, &mat->lfwork, &info));
942:     PetscCall(PetscFPTrapPop());
943:     mat->lfwork = (PetscInt)PetscRealPart(dummy);
944:     PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
945:   }
946:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
947:   PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, mat->fwork, &mat->lfwork, &info));
948:   PetscCall(PetscFPTrapPop());
949:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to QR factorization %d", (int)info);
950:   // TODO: try to estimate rank or test for and use geqp3 for rank revealing QR.  For now just say rank is min of m and n
951:   mat->rank = min;

953:   A->ops->solve    = MatSolve_SeqDense_QR;
954:   A->ops->matsolve = MatMatSolve_SeqDense_QR;
955:   A->factortype    = MAT_FACTOR_QR;
956:   if (m == n) {
957:     A->ops->solvetranspose    = MatSolveTranspose_SeqDense_QR;
958:     A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
959:   }

961:   PetscCall(PetscFree(A->solvertype));
962:   PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));

964:   PetscCall(PetscLogFlops(2.0 * min * min * (max - min / 3.0)));
965:   PetscFunctionReturn(PETSC_SUCCESS);
966: }

968: static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info_dummy)
969: {
970:   MatFactorInfo info;

972:   PetscFunctionBegin;
973:   info.fill = 1.0;

975:   PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
976:   PetscUseMethod(fact, "MatQRFactor_C", (Mat, IS, const MatFactorInfo *), (fact, NULL, &info));
977:   PetscFunctionReturn(PETSC_SUCCESS);
978: }

980: PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info)
981: {
982:   PetscFunctionBegin;
983:   fact->assembled    = PETSC_TRUE;
984:   fact->preallocated = PETSC_TRUE;
985:   PetscCall(PetscObjectComposeFunction((PetscObject)fact, "MatQRFactorNumeric_C", MatQRFactorNumeric_SeqDense));
986:   PetscFunctionReturn(PETSC_SUCCESS);
987: }

989: /* uses LAPACK */
990: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A, MatFactorType ftype, Mat *fact)
991: {
992:   PetscFunctionBegin;
993:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), fact));
994:   PetscCall(MatSetSizes(*fact, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
995:   PetscCall(MatSetType(*fact, MATDENSE));
996:   (*fact)->trivialsymbolic = PETSC_TRUE;
997:   if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
998:     (*fact)->ops->lufactorsymbolic  = MatLUFactorSymbolic_SeqDense;
999:     (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
1000:   } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
1001:     (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
1002:   } else if (ftype == MAT_FACTOR_QR) {
1003:     PetscCall(PetscObjectComposeFunction((PetscObject)*fact, "MatQRFactorSymbolic_C", MatQRFactorSymbolic_SeqDense));
1004:   }
1005:   (*fact)->factortype = ftype;

1007:   PetscCall(PetscFree((*fact)->solvertype));
1008:   PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*fact)->solvertype));
1009:   PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_LU]));
1010:   PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ILU]));
1011:   PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_CHOLESKY]));
1012:   PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ICC]));
1013:   PetscFunctionReturn(PETSC_SUCCESS);
1014: }

1016: static PetscErrorCode MatSOR_SeqDense(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal shift, PetscInt its, PetscInt lits, Vec xx)
1017: {
1018:   Mat_SeqDense      *mat = (Mat_SeqDense *)A->data;
1019:   PetscScalar       *x, *v = mat->v, zero = 0.0, xt;
1020:   const PetscScalar *b;
1021:   PetscInt           m = A->rmap->n, i;
1022:   PetscBLASInt       o = 1, bm = 0;

1024:   PetscFunctionBegin;
1025: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1026:   PetscCheck(A->offloadmask != PETSC_OFFLOAD_GPU, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented");
1027: #endif
1028:   if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
1029:   PetscCall(PetscBLASIntCast(m, &bm));
1030:   if (flag & SOR_ZERO_INITIAL_GUESS) {
1031:     /* this is a hack fix, should have another version without the second BLASdotu */
1032:     PetscCall(VecSet(xx, zero));
1033:   }
1034:   PetscCall(VecGetArray(xx, &x));
1035:   PetscCall(VecGetArrayRead(bb, &b));
1036:   its = its * lits;
1037:   PetscCheck(its > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Relaxation requires global its %" PetscInt_FMT " and local its %" PetscInt_FMT " both positive", its, lits);
1038:   while (its--) {
1039:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1040:       for (i = 0; i < m; i++) {
1041:         PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o));
1042:         x[i] = (1. - omega) * x[i] + omega * (xt + v[i + i * m] * x[i]) / (v[i + i * m] + shift);
1043:       }
1044:     }
1045:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1046:       for (i = m - 1; i >= 0; i--) {
1047:         PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o));
1048:         x[i] = (1. - omega) * x[i] + omega * (xt + v[i + i * m] * x[i]) / (v[i + i * m] + shift);
1049:       }
1050:     }
1051:   }
1052:   PetscCall(VecRestoreArrayRead(bb, &b));
1053:   PetscCall(VecRestoreArray(xx, &x));
1054:   PetscFunctionReturn(PETSC_SUCCESS);
1055: }

1057: static PetscErrorCode MatMultColumnRangeKernel_SeqDense(Mat A, Vec xx, Vec yy, PetscInt c_start, PetscInt c_end, PetscBool trans, PetscBool herm)
1058: {
1059:   Mat_SeqDense      *mat = (Mat_SeqDense *)A->data;
1060:   PetscScalar       *y, _DOne = 1.0, _DZero = 0.0;
1061:   PetscBLASInt       m, n, _One             = 1;
1062:   const PetscScalar *v = mat->v, *x;

1064:   PetscFunctionBegin;
1065:   PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1066:   PetscCall(PetscBLASIntCast(c_end - c_start, &n));
1067:   PetscCall(VecGetArrayRead(xx, &x));
1068:   PetscCall(VecGetArrayWrite(yy, &y));
1069:   if (!m || !n) {
1070:     PetscBLASInt i;
1071:     if (trans)
1072:       for (i = 0; i < n; i++) y[i] = 0.0;
1073:     else
1074:       for (i = 0; i < m; i++) y[i] = 0.0;
1075:   } else {
1076:     if (trans) {
1077:       if (herm) PetscCallBLAS("BLASgemv", BLASgemv_("C", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DZero, y + c_start, &_One));
1078:       else PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DZero, y + c_start, &_One));
1079:     } else {
1080:       PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x + c_start, &_One, &_DZero, y, &_One));
1081:     }
1082:     PetscCall(PetscLogFlops(2.0 * m * n - n));
1083:   }
1084:   PetscCall(VecRestoreArrayRead(xx, &x));
1085:   PetscCall(VecRestoreArrayWrite(yy, &y));
1086:   PetscFunctionReturn(PETSC_SUCCESS);
1087: }

1089: PetscErrorCode MatMultHermitianTransposeColumnRange_SeqDense(Mat A, Vec xx, Vec yy, PetscInt c_start, PetscInt c_end)
1090: {
1091:   PetscFunctionBegin;
1092:   PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, c_start, c_end, PETSC_TRUE, PETSC_TRUE));
1093:   PetscFunctionReturn(PETSC_SUCCESS);
1094: }

1096: PetscErrorCode MatMult_SeqDense(Mat A, Vec xx, Vec yy)
1097: {
1098:   PetscFunctionBegin;
1099:   PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_FALSE, PETSC_FALSE));
1100:   PetscFunctionReturn(PETSC_SUCCESS);
1101: }

1103: PetscErrorCode MatMultTranspose_SeqDense(Mat A, Vec xx, Vec yy)
1104: {
1105:   PetscFunctionBegin;
1106:   PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_FALSE));
1107:   PetscFunctionReturn(PETSC_SUCCESS);
1108: }

1110: PetscErrorCode MatMultHermitianTranspose_SeqDense(Mat A, Vec xx, Vec yy)
1111: {
1112:   PetscFunctionBegin;
1113:   PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_TRUE));
1114:   PetscFunctionReturn(PETSC_SUCCESS);
1115: }

1117: static PetscErrorCode MatMultAddColumnRangeKernel_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end, PetscBool trans, PetscBool herm)
1118: {
1119:   Mat_SeqDense      *mat = (Mat_SeqDense *)A->data;
1120:   const PetscScalar *v   = mat->v, *x;
1121:   PetscScalar       *y, _DOne = 1.0;
1122:   PetscBLASInt       m, n, _One = 1;

1124:   PetscFunctionBegin;
1125:   PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1126:   PetscCall(PetscBLASIntCast(c_end - c_start, &n));
1127:   PetscCall(VecCopy(zz, yy));
1128:   if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS);
1129:   PetscCall(VecGetArray(yy, &y));
1130:   PetscCall(VecGetArrayRead(xx, &x));
1131:   if (trans) {
1132:     if (herm) PetscCallBLAS("BLASgemv", BLASgemv_("C", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DOne, y + c_start, &_One));
1133:     else PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DOne, y + c_start, &_One));
1134:   } else {
1135:     PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x + c_start, &_One, &_DOne, y, &_One));
1136:   }
1137:   PetscCall(VecRestoreArrayRead(xx, &x));
1138:   PetscCall(VecRestoreArray(yy, &y));
1139:   PetscCall(PetscLogFlops(2.0 * m * n));
1140:   PetscFunctionReturn(PETSC_SUCCESS);
1141: }

1143: PetscErrorCode MatMultAddColumnRange_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end)
1144: {
1145:   PetscFunctionBegin;
1146:   PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, c_start, c_end, PETSC_FALSE, PETSC_FALSE));
1147:   PetscFunctionReturn(PETSC_SUCCESS);
1148: }

1150: PetscErrorCode MatMultHermitianTransposeAddColumnRange_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end)
1151: {
1152:   PetscFunctionBegin;
1153:   PetscMPIInt rank;
1154:   PetscCallMPI(MPI_Comm_rank(MPI_COMM_WORLD, &rank));
1155:   PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, c_start, c_end, PETSC_TRUE, PETSC_TRUE));
1156:   PetscFunctionReturn(PETSC_SUCCESS);
1157: }

1159: PetscErrorCode MatMultAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1160: {
1161:   PetscFunctionBegin;
1162:   PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_FALSE, PETSC_FALSE));
1163:   PetscFunctionReturn(PETSC_SUCCESS);
1164: }

1166: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1167: {
1168:   PetscFunctionBegin;
1169:   PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_FALSE));
1170:   PetscFunctionReturn(PETSC_SUCCESS);
1171: }

1173: PetscErrorCode MatMultHermitianTransposeAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1174: {
1175:   PetscFunctionBegin;
1176:   PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_TRUE));
1177:   PetscFunctionReturn(PETSC_SUCCESS);
1178: }

1180: static PetscErrorCode MatGetRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1181: {
1182:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1183:   PetscInt      i;

1185:   PetscFunctionBegin;
1186:   if (ncols) *ncols = A->cmap->n;
1187:   if (cols) {
1188:     PetscCall(PetscMalloc1(A->cmap->n, cols));
1189:     for (i = 0; i < A->cmap->n; i++) (*cols)[i] = i;
1190:   }
1191:   if (vals) {
1192:     const PetscScalar *v;

1194:     PetscCall(MatDenseGetArrayRead(A, &v));
1195:     PetscCall(PetscMalloc1(A->cmap->n, vals));
1196:     v += row;
1197:     for (i = 0; i < A->cmap->n; i++) {
1198:       (*vals)[i] = *v;
1199:       v += mat->lda;
1200:     }
1201:     PetscCall(MatDenseRestoreArrayRead(A, &v));
1202:   }
1203:   PetscFunctionReturn(PETSC_SUCCESS);
1204: }

1206: static PetscErrorCode MatRestoreRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1207: {
1208:   PetscFunctionBegin;
1209:   if (cols) PetscCall(PetscFree(*cols));
1210:   if (vals) PetscCall(PetscFree(*vals));
1211:   PetscFunctionReturn(PETSC_SUCCESS);
1212: }

1214: static PetscErrorCode MatSetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], const PetscScalar v[], InsertMode addv)
1215: {
1216:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1217:   PetscScalar  *av;
1218:   PetscInt      i, j, idx = 0;
1219: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1220:   PetscOffloadMask oldf;
1221: #endif

1223:   PetscFunctionBegin;
1224:   PetscCall(MatDenseGetArray(A, &av));
1225:   if (!mat->roworiented) {
1226:     if (addv == INSERT_VALUES) {
1227:       for (j = 0; j < n; j++) {
1228:         if (indexn[j] < 0) {
1229:           idx += m;
1230:           continue;
1231:         }
1232:         PetscCheck(indexn[j] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, indexn[j], A->cmap->n - 1);
1233:         for (i = 0; i < m; i++) {
1234:           if (indexm[i] < 0) {
1235:             idx++;
1236:             continue;
1237:           }
1238:           PetscCheck(indexm[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, indexm[i], A->rmap->n - 1);
1239:           av[indexn[j] * mat->lda + indexm[i]] = v[idx++];
1240:         }
1241:       }
1242:     } else {
1243:       for (j = 0; j < n; j++) {
1244:         if (indexn[j] < 0) {
1245:           idx += m;
1246:           continue;
1247:         }
1248:         PetscCheck(indexn[j] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, indexn[j], A->cmap->n - 1);
1249:         for (i = 0; i < m; i++) {
1250:           if (indexm[i] < 0) {
1251:             idx++;
1252:             continue;
1253:           }
1254:           PetscCheck(indexm[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, indexm[i], A->rmap->n - 1);
1255:           av[indexn[j] * mat->lda + indexm[i]] += v[idx++];
1256:         }
1257:       }
1258:     }
1259:   } else {
1260:     if (addv == INSERT_VALUES) {
1261:       for (i = 0; i < m; i++) {
1262:         if (indexm[i] < 0) {
1263:           idx += n;
1264:           continue;
1265:         }
1266:         PetscCheck(indexm[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, indexm[i], A->rmap->n - 1);
1267:         for (j = 0; j < n; j++) {
1268:           if (indexn[j] < 0) {
1269:             idx++;
1270:             continue;
1271:           }
1272:           PetscCheck(indexn[j] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, indexn[j], A->cmap->n - 1);
1273:           av[indexn[j] * mat->lda + indexm[i]] = v[idx++];
1274:         }
1275:       }
1276:     } else {
1277:       for (i = 0; i < m; i++) {
1278:         if (indexm[i] < 0) {
1279:           idx += n;
1280:           continue;
1281:         }
1282:         PetscCheck(indexm[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, indexm[i], A->rmap->n - 1);
1283:         for (j = 0; j < n; j++) {
1284:           if (indexn[j] < 0) {
1285:             idx++;
1286:             continue;
1287:           }
1288:           PetscCheck(indexn[j] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, indexn[j], A->cmap->n - 1);
1289:           av[indexn[j] * mat->lda + indexm[i]] += v[idx++];
1290:         }
1291:       }
1292:     }
1293:   }
1294:   /* hack to prevent unneeded copy to the GPU while returning the array */
1295: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1296:   oldf           = A->offloadmask;
1297:   A->offloadmask = PETSC_OFFLOAD_GPU;
1298: #endif
1299:   PetscCall(MatDenseRestoreArray(A, &av));
1300: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1301:   A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1302: #endif
1303:   PetscFunctionReturn(PETSC_SUCCESS);
1304: }

1306: static PetscErrorCode MatGetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], PetscScalar v[])
1307: {
1308:   Mat_SeqDense      *mat = (Mat_SeqDense *)A->data;
1309:   const PetscScalar *vv;
1310:   PetscInt           i, j;

1312:   PetscFunctionBegin;
1313:   PetscCall(MatDenseGetArrayRead(A, &vv));
1314:   /* row-oriented output */
1315:   for (i = 0; i < m; i++) {
1316:     if (indexm[i] < 0) {
1317:       v += n;
1318:       continue;
1319:     }
1320:     PetscCheck(indexm[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " requested larger than number rows %" PetscInt_FMT, indexm[i], A->rmap->n);
1321:     for (j = 0; j < n; j++) {
1322:       if (indexn[j] < 0) {
1323:         v++;
1324:         continue;
1325:       }
1326:       PetscCheck(indexn[j] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column %" PetscInt_FMT " requested larger than number columns %" PetscInt_FMT, indexn[j], A->cmap->n);
1327:       *v++ = vv[indexn[j] * mat->lda + indexm[i]];
1328:     }
1329:   }
1330:   PetscCall(MatDenseRestoreArrayRead(A, &vv));
1331:   PetscFunctionReturn(PETSC_SUCCESS);
1332: }

1334: PetscErrorCode MatView_Dense_Binary(Mat mat, PetscViewer viewer)
1335: {
1336:   PetscBool          skipHeader;
1337:   PetscViewerFormat  format;
1338:   PetscInt           header[4], M, N, m, lda, i, j, k;
1339:   const PetscScalar *v;
1340:   PetscScalar       *vwork;

1342:   PetscFunctionBegin;
1343:   PetscCall(PetscViewerSetUp(viewer));
1344:   PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
1345:   PetscCall(PetscViewerGetFormat(viewer, &format));
1346:   if (skipHeader) format = PETSC_VIEWER_NATIVE;

1348:   PetscCall(MatGetSize(mat, &M, &N));

1350:   /* write matrix header */
1351:   header[0] = MAT_FILE_CLASSID;
1352:   header[1] = M;
1353:   header[2] = N;
1354:   header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M * N;
1355:   if (!skipHeader) PetscCall(PetscViewerBinaryWrite(viewer, header, 4, PETSC_INT));

1357:   PetscCall(MatGetLocalSize(mat, &m, NULL));
1358:   if (format != PETSC_VIEWER_NATIVE) {
1359:     PetscInt nnz = m * N, *iwork;
1360:     /* store row lengths for each row */
1361:     PetscCall(PetscMalloc1(nnz, &iwork));
1362:     for (i = 0; i < m; i++) iwork[i] = N;
1363:     PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1364:     /* store column indices (zero start index) */
1365:     for (k = 0, i = 0; i < m; i++)
1366:       for (j = 0; j < N; j++, k++) iwork[k] = j;
1367:     PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1368:     PetscCall(PetscFree(iwork));
1369:   }
1370:   /* store matrix values as a dense matrix in row major order */
1371:   PetscCall(PetscMalloc1(m * N, &vwork));
1372:   PetscCall(MatDenseGetArrayRead(mat, &v));
1373:   PetscCall(MatDenseGetLDA(mat, &lda));
1374:   for (k = 0, i = 0; i < m; i++)
1375:     for (j = 0; j < N; j++, k++) vwork[k] = v[i + lda * j];
1376:   PetscCall(MatDenseRestoreArrayRead(mat, &v));
1377:   PetscCall(PetscViewerBinaryWriteAll(viewer, vwork, m * N, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1378:   PetscCall(PetscFree(vwork));
1379:   PetscFunctionReturn(PETSC_SUCCESS);
1380: }

1382: PetscErrorCode MatLoad_Dense_Binary(Mat mat, PetscViewer viewer)
1383: {
1384:   PetscBool    skipHeader;
1385:   PetscInt     header[4], M, N, m, nz, lda, i, j, k;
1386:   PetscInt     rows, cols;
1387:   PetscScalar *v, *vwork;

1389:   PetscFunctionBegin;
1390:   PetscCall(PetscViewerSetUp(viewer));
1391:   PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));

1393:   if (!skipHeader) {
1394:     PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT));
1395:     PetscCheck(header[0] == MAT_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file");
1396:     M = header[1];
1397:     N = header[2];
1398:     PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M);
1399:     PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N);
1400:     nz = header[3];
1401:     PetscCheck(nz == MATRIX_BINARY_FORMAT_DENSE || nz >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Unknown matrix format %" PetscInt_FMT " in file", nz);
1402:   } else {
1403:     PetscCall(MatGetSize(mat, &M, &N));
1404:     PetscCheck(M >= 0 && N >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "Matrix binary file header was skipped, thus the user must specify the global sizes of input matrix");
1405:     nz = MATRIX_BINARY_FORMAT_DENSE;
1406:   }

1408:   /* setup global sizes if not set */
1409:   if (mat->rmap->N < 0) mat->rmap->N = M;
1410:   if (mat->cmap->N < 0) mat->cmap->N = N;
1411:   PetscCall(MatSetUp(mat));
1412:   /* check if global sizes are correct */
1413:   PetscCall(MatGetSize(mat, &rows, &cols));
1414:   PetscCheck(M == rows && N == cols, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")", M, N, rows, cols);

1416:   PetscCall(MatGetSize(mat, NULL, &N));
1417:   PetscCall(MatGetLocalSize(mat, &m, NULL));
1418:   PetscCall(MatDenseGetArray(mat, &v));
1419:   PetscCall(MatDenseGetLDA(mat, &lda));
1420:   if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */
1421:     PetscInt nnz = m * N;
1422:     /* read in matrix values */
1423:     PetscCall(PetscMalloc1(nnz, &vwork));
1424:     PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1425:     /* store values in column major order */
1426:     for (j = 0; j < N; j++)
1427:       for (i = 0; i < m; i++) v[i + lda * j] = vwork[i * N + j];
1428:     PetscCall(PetscFree(vwork));
1429:   } else { /* matrix in file is sparse format */
1430:     PetscInt nnz = 0, *rlens, *icols;
1431:     /* read in row lengths */
1432:     PetscCall(PetscMalloc1(m, &rlens));
1433:     PetscCall(PetscViewerBinaryReadAll(viewer, rlens, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1434:     for (i = 0; i < m; i++) nnz += rlens[i];
1435:     /* read in column indices and values */
1436:     PetscCall(PetscMalloc2(nnz, &icols, nnz, &vwork));
1437:     PetscCall(PetscViewerBinaryReadAll(viewer, icols, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1438:     PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1439:     /* store values in column major order */
1440:     for (k = 0, i = 0; i < m; i++)
1441:       for (j = 0; j < rlens[i]; j++, k++) v[i + lda * icols[k]] = vwork[k];
1442:     PetscCall(PetscFree(rlens));
1443:     PetscCall(PetscFree2(icols, vwork));
1444:   }
1445:   PetscCall(MatDenseRestoreArray(mat, &v));
1446:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
1447:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
1448:   PetscFunctionReturn(PETSC_SUCCESS);
1449: }

1451: static PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1452: {
1453:   PetscBool isbinary, ishdf5;

1455:   PetscFunctionBegin;
1458:   /* force binary viewer to load .info file if it has not yet done so */
1459:   PetscCall(PetscViewerSetUp(viewer));
1460:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1461:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1462:   if (isbinary) {
1463:     PetscCall(MatLoad_Dense_Binary(newMat, viewer));
1464:   } else if (ishdf5) {
1465: #if defined(PETSC_HAVE_HDF5)
1466:     PetscCall(MatLoad_Dense_HDF5(newMat, viewer));
1467: #else
1468:     SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1469: #endif
1470:   } else {
1471:     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);
1472:   }
1473:   PetscFunctionReturn(PETSC_SUCCESS);
1474: }

1476: static PetscErrorCode MatView_SeqDense_ASCII(Mat A, PetscViewer viewer)
1477: {
1478:   Mat_SeqDense     *a = (Mat_SeqDense *)A->data;
1479:   PetscInt          i, j;
1480:   const char       *name;
1481:   PetscScalar      *v, *av;
1482:   PetscViewerFormat format;
1483: #if defined(PETSC_USE_COMPLEX)
1484:   PetscBool allreal = PETSC_TRUE;
1485: #endif

1487:   PetscFunctionBegin;
1488:   PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&av));
1489:   PetscCall(PetscViewerGetFormat(viewer, &format));
1490:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1491:     PetscFunctionReturn(PETSC_SUCCESS); /* do nothing for now */
1492:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1493:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1494:     for (i = 0; i < A->rmap->n; i++) {
1495:       v = av + i;
1496:       PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i));
1497:       for (j = 0; j < A->cmap->n; j++) {
1498: #if defined(PETSC_USE_COMPLEX)
1499:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1500:           PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", j, (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1501:         } else if (PetscRealPart(*v)) {
1502:           PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)PetscRealPart(*v)));
1503:         }
1504: #else
1505:         if (*v) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)*v));
1506: #endif
1507:         v += a->lda;
1508:       }
1509:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1510:     }
1511:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1512:   } else {
1513:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1514: #if defined(PETSC_USE_COMPLEX)
1515:     /* determine if matrix has all real values */
1516:     for (j = 0; j < A->cmap->n; j++) {
1517:       v = av + j * a->lda;
1518:       for (i = 0; i < A->rmap->n; i++) {
1519:         if (PetscImaginaryPart(v[i])) {
1520:           allreal = PETSC_FALSE;
1521:           break;
1522:         }
1523:       }
1524:     }
1525: #endif
1526:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1527:       PetscCall(PetscObjectGetName((PetscObject)A, &name));
1528:       PetscCall(PetscViewerASCIIPrintf(viewer, "%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n", A->rmap->n, A->cmap->n));
1529:       PetscCall(PetscViewerASCIIPrintf(viewer, "%s = zeros(%" PetscInt_FMT ",%" PetscInt_FMT ");\n", name, A->rmap->n, A->cmap->n));
1530:       PetscCall(PetscViewerASCIIPrintf(viewer, "%s = [\n", name));
1531:     }

1533:     for (i = 0; i < A->rmap->n; i++) {
1534:       v = av + i;
1535:       for (j = 0; j < A->cmap->n; j++) {
1536: #if defined(PETSC_USE_COMPLEX)
1537:         if (allreal) {
1538:           PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)PetscRealPart(*v)));
1539:         } else {
1540:           PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei ", (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1541:         }
1542: #else
1543:         PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)*v));
1544: #endif
1545:         v += a->lda;
1546:       }
1547:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1548:     }
1549:     if (format == PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer, "];\n"));
1550:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1551:   }
1552:   PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&av));
1553:   PetscCall(PetscViewerFlush(viewer));
1554:   PetscFunctionReturn(PETSC_SUCCESS);
1555: }

1557: #include <petscdraw.h>
1558: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw, void *Aa)
1559: {
1560:   Mat                A = (Mat)Aa;
1561:   PetscInt           m = A->rmap->n, n = A->cmap->n, i, j;
1562:   int                color = PETSC_DRAW_WHITE;
1563:   const PetscScalar *v;
1564:   PetscViewer        viewer;
1565:   PetscReal          xl, yl, xr, yr, x_l, x_r, y_l, y_r;
1566:   PetscViewerFormat  format;

1568:   PetscFunctionBegin;
1569:   PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer));
1570:   PetscCall(PetscViewerGetFormat(viewer, &format));
1571:   PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));

1573:   /* Loop over matrix elements drawing boxes */
1574:   PetscCall(MatDenseGetArrayRead(A, &v));
1575:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1576:     PetscDrawCollectiveBegin(draw);
1577:     /* Blue for negative and Red for positive */
1578:     for (j = 0; j < n; j++) {
1579:       x_l = j;
1580:       x_r = x_l + 1.0;
1581:       for (i = 0; i < m; i++) {
1582:         y_l = m - i - 1.0;
1583:         y_r = y_l + 1.0;
1584:         if (PetscRealPart(v[j * m + i]) > 0.) color = PETSC_DRAW_RED;
1585:         else if (PetscRealPart(v[j * m + i]) < 0.) color = PETSC_DRAW_BLUE;
1586:         else continue;
1587:         PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1588:       }
1589:     }
1590:     PetscDrawCollectiveEnd(draw);
1591:   } else {
1592:     /* use contour shading to indicate magnitude of values */
1593:     /* first determine max of all nonzero values */
1594:     PetscReal minv = 0.0, maxv = 0.0;
1595:     PetscDraw popup;

1597:     for (i = 0; i < m * n; i++) {
1598:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1599:     }
1600:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1601:     PetscCall(PetscDrawGetPopup(draw, &popup));
1602:     PetscCall(PetscDrawScalePopup(popup, minv, maxv));

1604:     PetscDrawCollectiveBegin(draw);
1605:     for (j = 0; j < n; j++) {
1606:       x_l = j;
1607:       x_r = x_l + 1.0;
1608:       for (i = 0; i < m; i++) {
1609:         y_l   = m - i - 1.0;
1610:         y_r   = y_l + 1.0;
1611:         color = PetscDrawRealToColor(PetscAbsScalar(v[j * m + i]), minv, maxv);
1612:         PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1613:       }
1614:     }
1615:     PetscDrawCollectiveEnd(draw);
1616:   }
1617:   PetscCall(MatDenseRestoreArrayRead(A, &v));
1618:   PetscFunctionReturn(PETSC_SUCCESS);
1619: }

1621: static PetscErrorCode MatView_SeqDense_Draw(Mat A, PetscViewer viewer)
1622: {
1623:   PetscDraw draw;
1624:   PetscBool isnull;
1625:   PetscReal xr, yr, xl, yl, h, w;

1627:   PetscFunctionBegin;
1628:   PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1629:   PetscCall(PetscDrawIsNull(draw, &isnull));
1630:   if (isnull) PetscFunctionReturn(PETSC_SUCCESS);

1632:   xr = A->cmap->n;
1633:   yr = A->rmap->n;
1634:   h  = yr / 10.0;
1635:   w  = xr / 10.0;
1636:   xr += w;
1637:   yr += h;
1638:   xl = -w;
1639:   yl = -h;
1640:   PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr));
1641:   PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer));
1642:   PetscCall(PetscDrawZoom(draw, MatView_SeqDense_Draw_Zoom, A));
1643:   PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL));
1644:   PetscCall(PetscDrawSave(draw));
1645:   PetscFunctionReturn(PETSC_SUCCESS);
1646: }

1648: PetscErrorCode MatView_SeqDense(Mat A, PetscViewer viewer)
1649: {
1650:   PetscBool iascii, isbinary, isdraw;

1652:   PetscFunctionBegin;
1653:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1654:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1655:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1656:   if (iascii) PetscCall(MatView_SeqDense_ASCII(A, viewer));
1657:   else if (isbinary) PetscCall(MatView_Dense_Binary(A, viewer));
1658:   else if (isdraw) PetscCall(MatView_SeqDense_Draw(A, viewer));
1659:   PetscFunctionReturn(PETSC_SUCCESS);
1660: }

1662: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A, const PetscScalar *array)
1663: {
1664:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

1666:   PetscFunctionBegin;
1667:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1668:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1669:   PetscCheck(!a->unplacedarray, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseResetArray() first");
1670:   a->unplacedarray       = a->v;
1671:   a->unplaced_user_alloc = a->user_alloc;
1672:   a->v                   = (PetscScalar *)array;
1673:   a->user_alloc          = PETSC_TRUE;
1674: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1675:   A->offloadmask = PETSC_OFFLOAD_CPU;
1676: #endif
1677:   PetscFunctionReturn(PETSC_SUCCESS);
1678: }

1680: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1681: {
1682:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

1684:   PetscFunctionBegin;
1685:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1686:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1687:   a->v             = a->unplacedarray;
1688:   a->user_alloc    = a->unplaced_user_alloc;
1689:   a->unplacedarray = NULL;
1690: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1691:   A->offloadmask = PETSC_OFFLOAD_CPU;
1692: #endif
1693:   PetscFunctionReturn(PETSC_SUCCESS);
1694: }

1696: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A, const PetscScalar *array)
1697: {
1698:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

1700:   PetscFunctionBegin;
1701:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1702:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1703:   if (!a->user_alloc) PetscCall(PetscFree(a->v));
1704:   a->v          = (PetscScalar *)array;
1705:   a->user_alloc = PETSC_FALSE;
1706: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1707:   A->offloadmask = PETSC_OFFLOAD_CPU;
1708: #endif
1709:   PetscFunctionReturn(PETSC_SUCCESS);
1710: }

1712: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1713: {
1714:   Mat_SeqDense *l = (Mat_SeqDense *)mat->data;

1716:   PetscFunctionBegin;
1717:   PetscCall(PetscLogObjectState((PetscObject)mat, "Rows %" PetscInt_FMT " Cols %" PetscInt_FMT, mat->rmap->n, mat->cmap->n));
1718:   PetscCall(VecDestroy(&l->qrrhs));
1719:   PetscCall(PetscFree(l->tau));
1720:   PetscCall(PetscFree(l->pivots));
1721:   PetscCall(PetscFree(l->fwork));
1722:   if (!l->user_alloc) PetscCall(PetscFree(l->v));
1723:   if (!l->unplaced_user_alloc) PetscCall(PetscFree(l->unplacedarray));
1724:   PetscCheck(!l->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1725:   PetscCheck(!l->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1726:   PetscCall(VecDestroy(&l->cvec));
1727:   PetscCall(MatDestroy(&l->cmat));
1728:   PetscCall(PetscFree(mat->data));

1730:   PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL));
1731:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactor_C", NULL));
1732:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorSymbolic_C", NULL));
1733:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorNumeric_C", NULL));
1734:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetLDA_C", NULL));
1735:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseSetLDA_C", NULL));
1736:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArray_C", NULL));
1737:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArray_C", NULL));
1738:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDensePlaceArray_C", NULL));
1739:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseResetArray_C", NULL));
1740:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseReplaceArray_C", NULL));
1741:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayRead_C", NULL));
1742:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayRead_C", NULL));
1743:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayWrite_C", NULL));
1744:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayWrite_C", NULL));
1745:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqaij_C", NULL));
1746: #if defined(PETSC_HAVE_ELEMENTAL)
1747:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_elemental_C", NULL));
1748: #endif
1749: #if defined(PETSC_HAVE_SCALAPACK)
1750:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_scalapack_C", NULL));
1751: #endif
1752: #if defined(PETSC_HAVE_CUDA)
1753:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensecuda_C", NULL));
1754:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", NULL));
1755:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdense_C", NULL));
1756:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensecuda_C", NULL));
1757: #endif
1758: #if defined(PETSC_HAVE_HIP)
1759:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensehip_C", NULL));
1760:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", NULL));
1761:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdense_C", NULL));
1762:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensehip_C", NULL));
1763: #endif
1764:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatSeqDenseSetPreallocation_C", NULL));
1765:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqaij_seqdense_C", NULL));
1766:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdense_C", NULL));
1767:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqbaij_seqdense_C", NULL));
1768:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqsbaij_seqdense_C", NULL));

1770:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumn_C", NULL));
1771:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumn_C", NULL));
1772:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVec_C", NULL));
1773:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVec_C", NULL));
1774:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecRead_C", NULL));
1775:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecRead_C", NULL));
1776:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecWrite_C", NULL));
1777:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecWrite_C", NULL));
1778:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetSubMatrix_C", NULL));
1779:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreSubMatrix_C", NULL));
1780:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultAddColumnRange_C", NULL));
1781:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeColumnRange_C", NULL));
1782:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeAddColumnRange_C", NULL));
1783:   PetscFunctionReturn(PETSC_SUCCESS);
1784: }

1786: static PetscErrorCode MatTranspose_SeqDense(Mat A, MatReuse reuse, Mat *matout)
1787: {
1788:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1789:   PetscInt      k, j, m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1790:   PetscScalar  *v, tmp;

1792:   PetscFunctionBegin;
1793:   if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *matout));
1794:   if (reuse == MAT_INPLACE_MATRIX) {
1795:     if (m == n) { /* in place transpose */
1796:       PetscCall(MatDenseGetArray(A, &v));
1797:       for (j = 0; j < m; j++) {
1798:         for (k = 0; k < j; k++) {
1799:           tmp          = v[j + k * M];
1800:           v[j + k * M] = v[k + j * M];
1801:           v[k + j * M] = tmp;
1802:         }
1803:       }
1804:       PetscCall(MatDenseRestoreArray(A, &v));
1805:     } else { /* reuse memory, temporary allocates new memory */
1806:       PetscScalar *v2;
1807:       PetscLayout  tmplayout;

1809:       PetscCall(PetscMalloc1((size_t)m * n, &v2));
1810:       PetscCall(MatDenseGetArray(A, &v));
1811:       for (j = 0; j < n; j++) {
1812:         for (k = 0; k < m; k++) v2[j + (size_t)k * n] = v[k + (size_t)j * M];
1813:       }
1814:       PetscCall(PetscArraycpy(v, v2, (size_t)m * n));
1815:       PetscCall(PetscFree(v2));
1816:       PetscCall(MatDenseRestoreArray(A, &v));
1817:       /* cleanup size dependent quantities */
1818:       PetscCall(VecDestroy(&mat->cvec));
1819:       PetscCall(MatDestroy(&mat->cmat));
1820:       PetscCall(PetscFree(mat->pivots));
1821:       PetscCall(PetscFree(mat->fwork));
1822:       /* swap row/col layouts */
1823:       mat->lda  = n;
1824:       tmplayout = A->rmap;
1825:       A->rmap   = A->cmap;
1826:       A->cmap   = tmplayout;
1827:     }
1828:   } else { /* out-of-place transpose */
1829:     Mat           tmat;
1830:     Mat_SeqDense *tmatd;
1831:     PetscScalar  *v2;
1832:     PetscInt      M2;

1834:     if (reuse == MAT_INITIAL_MATRIX) {
1835:       PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &tmat));
1836:       PetscCall(MatSetSizes(tmat, A->cmap->n, A->rmap->n, A->cmap->n, A->rmap->n));
1837:       PetscCall(MatSetType(tmat, ((PetscObject)A)->type_name));
1838:       PetscCall(MatSeqDenseSetPreallocation(tmat, NULL));
1839:     } else tmat = *matout;

1841:     PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&v));
1842:     PetscCall(MatDenseGetArray(tmat, &v2));
1843:     tmatd = (Mat_SeqDense *)tmat->data;
1844:     M2    = tmatd->lda;
1845:     for (j = 0; j < n; j++) {
1846:       for (k = 0; k < m; k++) v2[j + k * M2] = v[k + j * M];
1847:     }
1848:     PetscCall(MatDenseRestoreArray(tmat, &v2));
1849:     PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&v));
1850:     PetscCall(MatAssemblyBegin(tmat, MAT_FINAL_ASSEMBLY));
1851:     PetscCall(MatAssemblyEnd(tmat, MAT_FINAL_ASSEMBLY));
1852:     *matout = tmat;
1853:   }
1854:   PetscFunctionReturn(PETSC_SUCCESS);
1855: }

1857: static PetscErrorCode MatEqual_SeqDense(Mat A1, Mat A2, PetscBool *flg)
1858: {
1859:   Mat_SeqDense      *mat1 = (Mat_SeqDense *)A1->data;
1860:   Mat_SeqDense      *mat2 = (Mat_SeqDense *)A2->data;
1861:   PetscInt           i;
1862:   const PetscScalar *v1, *v2;

1864:   PetscFunctionBegin;
1865:   if (A1->rmap->n != A2->rmap->n) {
1866:     *flg = PETSC_FALSE;
1867:     PetscFunctionReturn(PETSC_SUCCESS);
1868:   }
1869:   if (A1->cmap->n != A2->cmap->n) {
1870:     *flg = PETSC_FALSE;
1871:     PetscFunctionReturn(PETSC_SUCCESS);
1872:   }
1873:   PetscCall(MatDenseGetArrayRead(A1, &v1));
1874:   PetscCall(MatDenseGetArrayRead(A2, &v2));
1875:   for (i = 0; i < A1->cmap->n; i++) {
1876:     PetscCall(PetscArraycmp(v1, v2, A1->rmap->n, flg));
1877:     if (*flg == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS);
1878:     v1 += mat1->lda;
1879:     v2 += mat2->lda;
1880:   }
1881:   PetscCall(MatDenseRestoreArrayRead(A1, &v1));
1882:   PetscCall(MatDenseRestoreArrayRead(A2, &v2));
1883:   *flg = PETSC_TRUE;
1884:   PetscFunctionReturn(PETSC_SUCCESS);
1885: }

1887: PetscErrorCode MatGetDiagonal_SeqDense(Mat A, Vec v)
1888: {
1889:   Mat_SeqDense      *mat = (Mat_SeqDense *)A->data;
1890:   PetscInt           i, n, len;
1891:   PetscScalar       *x;
1892:   const PetscScalar *vv;

1894:   PetscFunctionBegin;
1895:   PetscCall(VecGetSize(v, &n));
1896:   PetscCall(VecGetArray(v, &x));
1897:   len = PetscMin(A->rmap->n, A->cmap->n);
1898:   PetscCall(MatDenseGetArrayRead(A, &vv));
1899:   PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming mat and vec");
1900:   for (i = 0; i < len; i++) x[i] = vv[i * mat->lda + i];
1901:   PetscCall(MatDenseRestoreArrayRead(A, &vv));
1902:   PetscCall(VecRestoreArray(v, &x));
1903:   PetscFunctionReturn(PETSC_SUCCESS);
1904: }

1906: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A, Vec ll, Vec rr)
1907: {
1908:   Mat_SeqDense      *mat = (Mat_SeqDense *)A->data;
1909:   const PetscScalar *l, *r;
1910:   PetscScalar        x, *v, *vv;
1911:   PetscInt           i, j, m = A->rmap->n, n = A->cmap->n;

1913:   PetscFunctionBegin;
1914:   PetscCall(MatDenseGetArray(A, &vv));
1915:   if (ll) {
1916:     PetscCall(VecGetSize(ll, &m));
1917:     PetscCall(VecGetArrayRead(ll, &l));
1918:     PetscCheck(m == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vec wrong size");
1919:     for (i = 0; i < m; i++) {
1920:       x = l[i];
1921:       v = vv + i;
1922:       for (j = 0; j < n; j++) {
1923:         (*v) *= x;
1924:         v += mat->lda;
1925:       }
1926:     }
1927:     PetscCall(VecRestoreArrayRead(ll, &l));
1928:     PetscCall(PetscLogFlops(1.0 * n * m));
1929:   }
1930:   if (rr) {
1931:     PetscCall(VecGetSize(rr, &n));
1932:     PetscCall(VecGetArrayRead(rr, &r));
1933:     PetscCheck(n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vec wrong size");
1934:     for (i = 0; i < n; i++) {
1935:       x = r[i];
1936:       v = vv + i * mat->lda;
1937:       for (j = 0; j < m; j++) (*v++) *= x;
1938:     }
1939:     PetscCall(VecRestoreArrayRead(rr, &r));
1940:     PetscCall(PetscLogFlops(1.0 * n * m));
1941:   }
1942:   PetscCall(MatDenseRestoreArray(A, &vv));
1943:   PetscFunctionReturn(PETSC_SUCCESS);
1944: }

1946: PetscErrorCode MatNorm_SeqDense(Mat A, NormType type, PetscReal *nrm)
1947: {
1948:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1949:   PetscScalar  *v, *vv;
1950:   PetscReal     sum = 0.0;
1951:   PetscInt      lda, m = A->rmap->n, i, j;

1953:   PetscFunctionBegin;
1954:   PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&vv));
1955:   PetscCall(MatDenseGetLDA(A, &lda));
1956:   v = vv;
1957:   if (type == NORM_FROBENIUS) {
1958:     if (lda > m) {
1959:       for (j = 0; j < A->cmap->n; j++) {
1960:         v = vv + j * lda;
1961:         for (i = 0; i < m; i++) {
1962:           sum += PetscRealPart(PetscConj(*v) * (*v));
1963:           v++;
1964:         }
1965:       }
1966:     } else {
1967: #if defined(PETSC_USE_REAL___FP16)
1968:       PetscBLASInt one = 1, cnt = A->cmap->n * A->rmap->n;
1969:       PetscCallBLAS("BLASnrm2", *nrm = BLASnrm2_(&cnt, v, &one));
1970:     }
1971: #else
1972:       for (i = 0; i < A->cmap->n * A->rmap->n; i++) {
1973:         sum += PetscRealPart(PetscConj(*v) * (*v));
1974:         v++;
1975:       }
1976:     }
1977:     *nrm = PetscSqrtReal(sum);
1978: #endif
1979:     PetscCall(PetscLogFlops(2.0 * A->cmap->n * A->rmap->n));
1980:   } else if (type == NORM_1) {
1981:     *nrm = 0.0;
1982:     for (j = 0; j < A->cmap->n; j++) {
1983:       v   = vv + j * mat->lda;
1984:       sum = 0.0;
1985:       for (i = 0; i < A->rmap->n; i++) {
1986:         sum += PetscAbsScalar(*v);
1987:         v++;
1988:       }
1989:       if (sum > *nrm) *nrm = sum;
1990:     }
1991:     PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
1992:   } else if (type == NORM_INFINITY) {
1993:     *nrm = 0.0;
1994:     for (j = 0; j < A->rmap->n; j++) {
1995:       v   = vv + j;
1996:       sum = 0.0;
1997:       for (i = 0; i < A->cmap->n; i++) {
1998:         sum += PetscAbsScalar(*v);
1999:         v += mat->lda;
2000:       }
2001:       if (sum > *nrm) *nrm = sum;
2002:     }
2003:     PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
2004:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No two norm");
2005:   PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&vv));
2006:   PetscFunctionReturn(PETSC_SUCCESS);
2007: }

2009: static PetscErrorCode MatSetOption_SeqDense(Mat A, MatOption op, PetscBool flg)
2010: {
2011:   Mat_SeqDense *aij = (Mat_SeqDense *)A->data;

2013:   PetscFunctionBegin;
2014:   switch (op) {
2015:   case MAT_ROW_ORIENTED:
2016:     aij->roworiented = flg;
2017:     break;
2018:   case MAT_NEW_NONZERO_LOCATIONS:
2019:   case MAT_NEW_NONZERO_LOCATION_ERR:
2020:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
2021:   case MAT_FORCE_DIAGONAL_ENTRIES:
2022:   case MAT_KEEP_NONZERO_PATTERN:
2023:   case MAT_IGNORE_OFF_PROC_ENTRIES:
2024:   case MAT_USE_HASH_TABLE:
2025:   case MAT_IGNORE_ZERO_ENTRIES:
2026:   case MAT_IGNORE_LOWER_TRIANGULAR:
2027:   case MAT_SORTED_FULL:
2028:     PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op]));
2029:     break;
2030:   case MAT_SPD:
2031:   case MAT_SYMMETRIC:
2032:   case MAT_STRUCTURALLY_SYMMETRIC:
2033:   case MAT_HERMITIAN:
2034:   case MAT_SYMMETRY_ETERNAL:
2035:   case MAT_STRUCTURAL_SYMMETRY_ETERNAL:
2036:   case MAT_SPD_ETERNAL:
2037:     break;
2038:   default:
2039:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unknown option %s", MatOptions[op]);
2040:   }
2041:   PetscFunctionReturn(PETSC_SUCCESS);
2042: }

2044: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
2045: {
2046:   Mat_SeqDense *l   = (Mat_SeqDense *)A->data;
2047:   PetscInt      lda = l->lda, m = A->rmap->n, n = A->cmap->n, j;
2048:   PetscScalar  *v;

2050:   PetscFunctionBegin;
2051:   PetscCall(MatDenseGetArrayWrite(A, &v));
2052:   if (lda > m) {
2053:     for (j = 0; j < n; j++) PetscCall(PetscArrayzero(v + j * lda, m));
2054:   } else {
2055:     PetscCall(PetscArrayzero(v, PetscInt64Mult(m, n)));
2056:   }
2057:   PetscCall(MatDenseRestoreArrayWrite(A, &v));
2058:   PetscFunctionReturn(PETSC_SUCCESS);
2059: }

2061: static PetscErrorCode MatZeroRows_SeqDense(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2062: {
2063:   Mat_SeqDense      *l = (Mat_SeqDense *)A->data;
2064:   PetscInt           m = l->lda, n = A->cmap->n, i, j;
2065:   PetscScalar       *slot, *bb, *v;
2066:   const PetscScalar *xx;

2068:   PetscFunctionBegin;
2069:   if (PetscDefined(USE_DEBUG)) {
2070:     for (i = 0; i < N; i++) {
2071:       PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
2072:       PetscCheck(rows[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " requested to be zeroed greater than or equal number of rows %" PetscInt_FMT, rows[i], A->rmap->n);
2073:     }
2074:   }
2075:   if (!N) PetscFunctionReturn(PETSC_SUCCESS);

2077:   /* fix right-hand side if needed */
2078:   if (x && b) {
2079:     PetscCall(VecGetArrayRead(x, &xx));
2080:     PetscCall(VecGetArray(b, &bb));
2081:     for (i = 0; i < N; i++) bb[rows[i]] = diag * xx[rows[i]];
2082:     PetscCall(VecRestoreArrayRead(x, &xx));
2083:     PetscCall(VecRestoreArray(b, &bb));
2084:   }

2086:   PetscCall(MatDenseGetArray(A, &v));
2087:   for (i = 0; i < N; i++) {
2088:     slot = v + rows[i];
2089:     for (j = 0; j < n; j++) {
2090:       *slot = 0.0;
2091:       slot += m;
2092:     }
2093:   }
2094:   if (diag != 0.0) {
2095:     PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
2096:     for (i = 0; i < N; i++) {
2097:       slot  = v + (m + 1) * rows[i];
2098:       *slot = diag;
2099:     }
2100:   }
2101:   PetscCall(MatDenseRestoreArray(A, &v));
2102:   PetscFunctionReturn(PETSC_SUCCESS);
2103: }

2105: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A, PetscInt *lda)
2106: {
2107:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;

2109:   PetscFunctionBegin;
2110:   *lda = mat->lda;
2111:   PetscFunctionReturn(PETSC_SUCCESS);
2112: }

2114: PetscErrorCode MatDenseGetArray_SeqDense(Mat A, PetscScalar **array)
2115: {
2116:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;

2118:   PetscFunctionBegin;
2119:   PetscCheck(!mat->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
2120:   *array = mat->v;
2121:   PetscFunctionReturn(PETSC_SUCCESS);
2122: }

2124: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A, PetscScalar **array)
2125: {
2126:   PetscFunctionBegin;
2127:   if (array) *array = NULL;
2128:   PetscFunctionReturn(PETSC_SUCCESS);
2129: }

2131: /*@
2132:   MatDenseGetLDA - gets the leading dimension of the array returned from `MatDenseGetArray()`

2134:   Not Collective

2136:   Input Parameter:
2137: . A - a `MATDENSE` or `MATDENSECUDA` matrix

2139:   Output Parameter:
2140: . lda - the leading dimension

2142:   Level: intermediate

2144: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseSetLDA()`
2145: @*/
2146: PetscErrorCode MatDenseGetLDA(Mat A, PetscInt *lda)
2147: {
2148:   PetscFunctionBegin;
2150:   PetscAssertPointer(lda, 2);
2151:   MatCheckPreallocated(A, 1);
2152:   PetscUseMethod(A, "MatDenseGetLDA_C", (Mat, PetscInt *), (A, lda));
2153:   PetscFunctionReturn(PETSC_SUCCESS);
2154: }

2156: /*@
2157:   MatDenseSetLDA - Sets the leading dimension of the array used by the `MATDENSE` matrix

2159:   Collective if the matrix layouts have not yet been setup

2161:   Input Parameters:
2162: + A   - a `MATDENSE` or `MATDENSECUDA` matrix
2163: - lda - the leading dimension

2165:   Level: intermediate

2167: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetLDA()`
2168: @*/
2169: PetscErrorCode MatDenseSetLDA(Mat A, PetscInt lda)
2170: {
2171:   PetscFunctionBegin;
2173:   PetscTryMethod(A, "MatDenseSetLDA_C", (Mat, PetscInt), (A, lda));
2174:   PetscFunctionReturn(PETSC_SUCCESS);
2175: }

2177: /*@C
2178:   MatDenseGetArray - gives read-write access to the array where the data for a `MATDENSE` matrix is stored

2180:   Logically Collective

2182:   Input Parameter:
2183: . A - a dense matrix

2185:   Output Parameter:
2186: . array - pointer to the data

2188:   Level: intermediate

2190:   Fortran Notes:
2191:   `MatDenseGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `MatDenseGetArrayF90()`

2193: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2194: @*/
2195: PetscErrorCode MatDenseGetArray(Mat A, PetscScalar **array)
2196: {
2197:   PetscFunctionBegin;
2199:   PetscAssertPointer(array, 2);
2200:   PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2201:   PetscFunctionReturn(PETSC_SUCCESS);
2202: }

2204: /*@C
2205:   MatDenseRestoreArray - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArray()`

2207:   Logically Collective

2209:   Input Parameters:
2210: + A     - a dense matrix
2211: - array - pointer to the data (may be `NULL`)

2213:   Level: intermediate

2215:   Fortran Notes:
2216:   `MatDenseRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `MatDenseRestoreArrayF90()`

2218: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2219: @*/
2220: PetscErrorCode MatDenseRestoreArray(Mat A, PetscScalar **array)
2221: {
2222:   PetscFunctionBegin;
2224:   if (array) PetscAssertPointer(array, 2);
2225:   PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2226:   PetscCall(PetscObjectStateIncrease((PetscObject)A));
2227: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2228:   A->offloadmask = PETSC_OFFLOAD_CPU;
2229: #endif
2230:   PetscFunctionReturn(PETSC_SUCCESS);
2231: }

2233: /*@C
2234:   MatDenseGetArrayRead - gives read-only access to the array where the data for a `MATDENSE` matrix is stored

2236:   Not Collective

2238:   Input Parameter:
2239: . A - a dense matrix

2241:   Output Parameter:
2242: . array - pointer to the data

2244:   Level: intermediate

2246: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2247: @*/
2248: PetscErrorCode MatDenseGetArrayRead(Mat A, const PetscScalar **array)
2249: {
2250:   PetscFunctionBegin;
2252:   PetscAssertPointer(array, 2);
2253:   PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2254:   PetscFunctionReturn(PETSC_SUCCESS);
2255: }

2257: /*@C
2258:   MatDenseRestoreArrayRead - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayRead()`

2260:   Not Collective

2262:   Input Parameters:
2263: + A     - a dense matrix
2264: - array - pointer to the data (may be `NULL`)

2266:   Level: intermediate

2268: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2269: @*/
2270: PetscErrorCode MatDenseRestoreArrayRead(Mat A, const PetscScalar **array)
2271: {
2272:   PetscFunctionBegin;
2274:   if (array) PetscAssertPointer(array, 2);
2275:   PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2276:   PetscFunctionReturn(PETSC_SUCCESS);
2277: }

2279: /*@C
2280:   MatDenseGetArrayWrite - gives write-only access to the array where the data for a `MATDENSE` matrix is stored

2282:   Not Collective

2284:   Input Parameter:
2285: . A - a dense matrix

2287:   Output Parameter:
2288: . array - pointer to the data

2290:   Level: intermediate

2292: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2293: @*/
2294: PetscErrorCode MatDenseGetArrayWrite(Mat A, PetscScalar **array)
2295: {
2296:   PetscFunctionBegin;
2298:   PetscAssertPointer(array, 2);
2299:   PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2300:   PetscFunctionReturn(PETSC_SUCCESS);
2301: }

2303: /*@C
2304:   MatDenseRestoreArrayWrite - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayWrite()`

2306:   Not Collective

2308:   Input Parameters:
2309: + A     - a dense matrix
2310: - array - pointer to the data (may be `NULL`)

2312:   Level: intermediate

2314: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2315: @*/
2316: PetscErrorCode MatDenseRestoreArrayWrite(Mat A, PetscScalar **array)
2317: {
2318:   PetscFunctionBegin;
2320:   if (array) PetscAssertPointer(array, 2);
2321:   PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2322:   PetscCall(PetscObjectStateIncrease((PetscObject)A));
2323: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2324:   A->offloadmask = PETSC_OFFLOAD_CPU;
2325: #endif
2326:   PetscFunctionReturn(PETSC_SUCCESS);
2327: }

2329: /*@C
2330:   MatDenseGetArrayAndMemType - gives read-write access to the array where the data for a `MATDENSE` matrix is stored

2332:   Logically Collective

2334:   Input Parameter:
2335: . A - a dense matrix

2337:   Output Parameters:
2338: + array - pointer to the data
2339: - mtype - memory type of the returned pointer

2341:   Level: intermediate

2343:   Note:
2344:   If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2345:   an array on device is always returned and is guaranteed to contain the matrix's latest data.

2347: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArrayRead()`,
2348:    `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2349: @*/
2350: PetscErrorCode MatDenseGetArrayAndMemType(Mat A, PetscScalar **array, PetscMemType *mtype)
2351: {
2352:   PetscBool isMPI;

2354:   PetscFunctionBegin;
2356:   PetscAssertPointer(array, 2);
2357:   PetscCall(MatBindToCPU(A, PETSC_FALSE)); /* We want device matrices to always return device arrays, so we unbind the matrix if it is bound to CPU */
2358:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2359:   if (isMPI) {
2360:     /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2361:     PetscCall(MatDenseGetArrayAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2362:   } else {
2363:     PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);

2365:     PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayAndMemType_C", &fptr));
2366:     if (fptr) {
2367:       PetscCall((*fptr)(A, array, mtype));
2368:     } else {
2369:       PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2370:       if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2371:     }
2372:   }
2373:   PetscFunctionReturn(PETSC_SUCCESS);
2374: }

2376: /*@C
2377:   MatDenseRestoreArrayAndMemType - returns access to the array that is obtained by `MatDenseGetArrayAndMemType()`

2379:   Logically Collective

2381:   Input Parameters:
2382: + A     - a dense matrix
2383: - array - pointer to the data

2385:   Level: intermediate

2387: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2388: @*/
2389: PetscErrorCode MatDenseRestoreArrayAndMemType(Mat A, PetscScalar **array)
2390: {
2391:   PetscBool isMPI;

2393:   PetscFunctionBegin;
2395:   PetscAssertPointer(array, 2);
2396:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2397:   if (isMPI) {
2398:     PetscCall(MatDenseRestoreArrayAndMemType(((Mat_MPIDense *)A->data)->A, array));
2399:   } else {
2400:     PetscErrorCode (*fptr)(Mat, PetscScalar **);

2402:     PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayAndMemType_C", &fptr));
2403:     if (fptr) {
2404:       PetscCall((*fptr)(A, array));
2405:     } else {
2406:       PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2407:     }
2408:     *array = NULL;
2409:   }
2410:   PetscCall(PetscObjectStateIncrease((PetscObject)A));
2411:   PetscFunctionReturn(PETSC_SUCCESS);
2412: }

2414: /*@C
2415:   MatDenseGetArrayReadAndMemType - gives read-only access to the array where the data for a `MATDENSE` matrix is stored

2417:   Logically Collective

2419:   Input Parameter:
2420: . A - a dense matrix

2422:   Output Parameters:
2423: + array - pointer to the data
2424: - mtype - memory type of the returned pointer

2426:   Level: intermediate

2428:   Note:
2429:   If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2430:   an array on device is always returned and is guaranteed to contain the matrix's latest data.

2432: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`,
2433:    `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2434: @*/
2435: PetscErrorCode MatDenseGetArrayReadAndMemType(Mat A, const PetscScalar **array, PetscMemType *mtype)
2436: {
2437:   PetscBool isMPI;

2439:   PetscFunctionBegin;
2441:   PetscAssertPointer(array, 2);
2442:   PetscCall(MatBindToCPU(A, PETSC_FALSE)); /* We want device matrices to always return device arrays, so we unbind the matrix if it is bound to CPU */
2443:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2444:   if (isMPI) { /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2445:     PetscCall(MatDenseGetArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2446:   } else {
2447:     PetscErrorCode (*fptr)(Mat, const PetscScalar **, PetscMemType *);

2449:     PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayReadAndMemType_C", &fptr));
2450:     if (fptr) {
2451:       PetscCall((*fptr)(A, array, mtype));
2452:     } else {
2453:       PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2454:       if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2455:     }
2456:   }
2457:   PetscFunctionReturn(PETSC_SUCCESS);
2458: }

2460: /*@C
2461:   MatDenseRestoreArrayReadAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()`

2463:   Logically Collective

2465:   Input Parameters:
2466: + A     - a dense matrix
2467: - array - pointer to the data

2469:   Level: intermediate

2471: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2472: @*/
2473: PetscErrorCode MatDenseRestoreArrayReadAndMemType(Mat A, const PetscScalar **array)
2474: {
2475:   PetscBool isMPI;

2477:   PetscFunctionBegin;
2479:   PetscAssertPointer(array, 2);
2480:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2481:   if (isMPI) {
2482:     PetscCall(MatDenseRestoreArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array));
2483:   } else {
2484:     PetscErrorCode (*fptr)(Mat, const PetscScalar **);

2486:     PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayReadAndMemType_C", &fptr));
2487:     if (fptr) {
2488:       PetscCall((*fptr)(A, array));
2489:     } else {
2490:       PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2491:     }
2492:     *array = NULL;
2493:   }
2494:   PetscFunctionReturn(PETSC_SUCCESS);
2495: }

2497: /*@C
2498:   MatDenseGetArrayWriteAndMemType - gives write-only access to the array where the data for a `MATDENSE` matrix is stored

2500:   Logically Collective

2502:   Input Parameter:
2503: . A - a dense matrix

2505:   Output Parameters:
2506: + array - pointer to the data
2507: - mtype - memory type of the returned pointer

2509:   Level: intermediate

2511:   Note:
2512:   If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2513:   an array on device is always returned and is guaranteed to contain the matrix's latest data.

2515: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWriteAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayRead()`,
2516:   `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2517: @*/
2518: PetscErrorCode MatDenseGetArrayWriteAndMemType(Mat A, PetscScalar **array, PetscMemType *mtype)
2519: {
2520:   PetscBool isMPI;

2522:   PetscFunctionBegin;
2524:   PetscAssertPointer(array, 2);
2525:   PetscCall(MatBindToCPU(A, PETSC_FALSE)); /* We want device matrices to always return device arrays, so we unbind the matrix if it is bound to CPU */
2526:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2527:   if (isMPI) {
2528:     PetscCall(MatDenseGetArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2529:   } else {
2530:     PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);

2532:     PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayWriteAndMemType_C", &fptr));
2533:     if (fptr) {
2534:       PetscCall((*fptr)(A, array, mtype));
2535:     } else {
2536:       PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2537:       if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2538:     }
2539:   }
2540:   PetscFunctionReturn(PETSC_SUCCESS);
2541: }

2543: /*@C
2544:   MatDenseRestoreArrayWriteAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()`

2546:   Logically Collective

2548:   Input Parameters:
2549: + A     - a dense matrix
2550: - array - pointer to the data

2552:   Level: intermediate

2554: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2555: @*/
2556: PetscErrorCode MatDenseRestoreArrayWriteAndMemType(Mat A, PetscScalar **array)
2557: {
2558:   PetscBool isMPI;

2560:   PetscFunctionBegin;
2562:   PetscAssertPointer(array, 2);
2563:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2564:   if (isMPI) {
2565:     PetscCall(MatDenseRestoreArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array));
2566:   } else {
2567:     PetscErrorCode (*fptr)(Mat, PetscScalar **);

2569:     PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayWriteAndMemType_C", &fptr));
2570:     if (fptr) {
2571:       PetscCall((*fptr)(A, array));
2572:     } else {
2573:       PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2574:     }
2575:     *array = NULL;
2576:   }
2577:   PetscCall(PetscObjectStateIncrease((PetscObject)A));
2578:   PetscFunctionReturn(PETSC_SUCCESS);
2579: }

2581: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B)
2582: {
2583:   Mat_SeqDense   *mat = (Mat_SeqDense *)A->data;
2584:   PetscInt        i, j, nrows, ncols, ldb;
2585:   const PetscInt *irow, *icol;
2586:   PetscScalar    *av, *bv, *v = mat->v;
2587:   Mat             newmat;

2589:   PetscFunctionBegin;
2590:   PetscCall(ISGetIndices(isrow, &irow));
2591:   PetscCall(ISGetIndices(iscol, &icol));
2592:   PetscCall(ISGetLocalSize(isrow, &nrows));
2593:   PetscCall(ISGetLocalSize(iscol, &ncols));

2595:   /* Check submatrixcall */
2596:   if (scall == MAT_REUSE_MATRIX) {
2597:     PetscInt n_cols, n_rows;
2598:     PetscCall(MatGetSize(*B, &n_rows, &n_cols));
2599:     if (n_rows != nrows || n_cols != ncols) {
2600:       /* resize the result matrix to match number of requested rows/columns */
2601:       PetscCall(MatSetSizes(*B, nrows, ncols, nrows, ncols));
2602:     }
2603:     newmat = *B;
2604:   } else {
2605:     /* Create and fill new matrix */
2606:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &newmat));
2607:     PetscCall(MatSetSizes(newmat, nrows, ncols, nrows, ncols));
2608:     PetscCall(MatSetType(newmat, ((PetscObject)A)->type_name));
2609:     PetscCall(MatSeqDenseSetPreallocation(newmat, NULL));
2610:   }

2612:   /* Now extract the data pointers and do the copy,column at a time */
2613:   PetscCall(MatDenseGetArray(newmat, &bv));
2614:   PetscCall(MatDenseGetLDA(newmat, &ldb));
2615:   for (i = 0; i < ncols; i++) {
2616:     av = v + mat->lda * icol[i];
2617:     for (j = 0; j < nrows; j++) bv[j] = av[irow[j]];
2618:     bv += ldb;
2619:   }
2620:   PetscCall(MatDenseRestoreArray(newmat, &bv));

2622:   /* Assemble the matrices so that the correct flags are set */
2623:   PetscCall(MatAssemblyBegin(newmat, MAT_FINAL_ASSEMBLY));
2624:   PetscCall(MatAssemblyEnd(newmat, MAT_FINAL_ASSEMBLY));

2626:   /* Free work space */
2627:   PetscCall(ISRestoreIndices(isrow, &irow));
2628:   PetscCall(ISRestoreIndices(iscol, &icol));
2629:   *B = newmat;
2630:   PetscFunctionReturn(PETSC_SUCCESS);
2631: }

2633: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[])
2634: {
2635:   PetscInt i;

2637:   PetscFunctionBegin;
2638:   if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n, B));

2640:   for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SeqDense(A, irow[i], icol[i], scall, &(*B)[i]));
2641:   PetscFunctionReturn(PETSC_SUCCESS);
2642: }

2644: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat, MatAssemblyType mode)
2645: {
2646:   PetscFunctionBegin;
2647:   PetscFunctionReturn(PETSC_SUCCESS);
2648: }

2650: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat, MatAssemblyType mode)
2651: {
2652:   PetscFunctionBegin;
2653:   PetscFunctionReturn(PETSC_SUCCESS);
2654: }

2656: PetscErrorCode MatCopy_SeqDense(Mat A, Mat B, MatStructure str)
2657: {
2658:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data;
2659:   const PetscScalar *va;
2660:   PetscScalar       *vb;
2661:   PetscInt           lda1 = a->lda, lda2 = b->lda, m = A->rmap->n, n = A->cmap->n, j;

2663:   PetscFunctionBegin;
2664:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2665:   if (A->ops->copy != B->ops->copy) {
2666:     PetscCall(MatCopy_Basic(A, B, str));
2667:     PetscFunctionReturn(PETSC_SUCCESS);
2668:   }
2669:   PetscCheck(m == B->rmap->n && n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "size(B) != size(A)");
2670:   PetscCall(MatDenseGetArrayRead(A, &va));
2671:   PetscCall(MatDenseGetArray(B, &vb));
2672:   if (lda1 > m || lda2 > m) {
2673:     for (j = 0; j < n; j++) PetscCall(PetscArraycpy(vb + j * lda2, va + j * lda1, m));
2674:   } else {
2675:     PetscCall(PetscArraycpy(vb, va, A->rmap->n * A->cmap->n));
2676:   }
2677:   PetscCall(MatDenseRestoreArray(B, &vb));
2678:   PetscCall(MatDenseRestoreArrayRead(A, &va));
2679:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2680:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2681:   PetscFunctionReturn(PETSC_SUCCESS);
2682: }

2684: PetscErrorCode MatSetUp_SeqDense(Mat A)
2685: {
2686:   PetscFunctionBegin;
2687:   PetscCall(PetscLayoutSetUp(A->rmap));
2688:   PetscCall(PetscLayoutSetUp(A->cmap));
2689:   if (!A->preallocated) PetscCall(MatSeqDenseSetPreallocation(A, NULL));
2690:   PetscFunctionReturn(PETSC_SUCCESS);
2691: }

2693: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2694: {
2695:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2696:   PetscInt      i, j;
2697:   PetscInt      min = PetscMin(A->rmap->n, A->cmap->n);
2698:   PetscScalar  *aa;

2700:   PetscFunctionBegin;
2701:   PetscCall(MatDenseGetArray(A, &aa));
2702:   for (j = 0; j < A->cmap->n; j++) {
2703:     for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscConj(aa[i + j * mat->lda]);
2704:   }
2705:   PetscCall(MatDenseRestoreArray(A, &aa));
2706:   if (mat->tau)
2707:     for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
2708:   PetscFunctionReturn(PETSC_SUCCESS);
2709: }

2711: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2712: {
2713:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2714:   PetscInt      i, j;
2715:   PetscScalar  *aa;

2717:   PetscFunctionBegin;
2718:   PetscCall(MatDenseGetArray(A, &aa));
2719:   for (j = 0; j < A->cmap->n; j++) {
2720:     for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscRealPart(aa[i + j * mat->lda]);
2721:   }
2722:   PetscCall(MatDenseRestoreArray(A, &aa));
2723:   PetscFunctionReturn(PETSC_SUCCESS);
2724: }

2726: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2727: {
2728:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2729:   PetscInt      i, j;
2730:   PetscScalar  *aa;

2732:   PetscFunctionBegin;
2733:   PetscCall(MatDenseGetArray(A, &aa));
2734:   for (j = 0; j < A->cmap->n; j++) {
2735:     for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscImaginaryPart(aa[i + j * mat->lda]);
2736:   }
2737:   PetscCall(MatDenseRestoreArray(A, &aa));
2738:   PetscFunctionReturn(PETSC_SUCCESS);
2739: }

2741: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2742: {
2743:   PetscInt  m = A->rmap->n, n = B->cmap->n;
2744:   PetscBool cisdense = PETSC_FALSE;

2746:   PetscFunctionBegin;
2747:   PetscCall(MatSetSizes(C, m, n, m, n));
2748: #if defined(PETSC_HAVE_CUDA)
2749:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2750: #endif
2751: #if defined(PETSC_HAVE_HIP)
2752:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2753: #endif
2754:   if (!cisdense) {
2755:     PetscBool flg;

2757:     PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2758:     PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2759:   }
2760:   PetscCall(MatSetUp(C));
2761:   PetscFunctionReturn(PETSC_SUCCESS);
2762: }

2764: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2765: {
2766:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data, *c = (Mat_SeqDense *)C->data;
2767:   PetscBLASInt       m, n, k;
2768:   const PetscScalar *av, *bv;
2769:   PetscScalar       *cv;
2770:   PetscScalar        _DOne = 1.0, _DZero = 0.0;

2772:   PetscFunctionBegin;
2773:   PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2774:   PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2775:   PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2776:   if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2777:   PetscCall(MatDenseGetArrayRead(A, &av));
2778:   PetscCall(MatDenseGetArrayRead(B, &bv));
2779:   PetscCall(MatDenseGetArrayWrite(C, &cv));
2780:   PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2781:   PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2782:   PetscCall(MatDenseRestoreArrayRead(A, &av));
2783:   PetscCall(MatDenseRestoreArrayRead(B, &bv));
2784:   PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2785:   PetscFunctionReturn(PETSC_SUCCESS);
2786: }

2788: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2789: {
2790:   PetscInt  m = A->rmap->n, n = B->rmap->n;
2791:   PetscBool cisdense = PETSC_FALSE;

2793:   PetscFunctionBegin;
2794:   PetscCall(MatSetSizes(C, m, n, m, n));
2795: #if defined(PETSC_HAVE_CUDA)
2796:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2797: #endif
2798: #if defined(PETSC_HAVE_HIP)
2799:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2800: #endif
2801:   if (!cisdense) {
2802:     PetscBool flg;

2804:     PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2805:     PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2806:   }
2807:   PetscCall(MatSetUp(C));
2808:   PetscFunctionReturn(PETSC_SUCCESS);
2809: }

2811: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2812: {
2813:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data;
2814:   Mat_SeqDense      *b = (Mat_SeqDense *)B->data;
2815:   Mat_SeqDense      *c = (Mat_SeqDense *)C->data;
2816:   const PetscScalar *av, *bv;
2817:   PetscScalar       *cv;
2818:   PetscBLASInt       m, n, k;
2819:   PetscScalar        _DOne = 1.0, _DZero = 0.0;

2821:   PetscFunctionBegin;
2822:   PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2823:   PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2824:   PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2825:   if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2826:   PetscCall(MatDenseGetArrayRead(A, &av));
2827:   PetscCall(MatDenseGetArrayRead(B, &bv));
2828:   PetscCall(MatDenseGetArrayWrite(C, &cv));
2829:   PetscCallBLAS("BLASgemm", BLASgemm_("N", "T", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2830:   PetscCall(MatDenseRestoreArrayRead(A, &av));
2831:   PetscCall(MatDenseRestoreArrayRead(B, &bv));
2832:   PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2833:   PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2834:   PetscFunctionReturn(PETSC_SUCCESS);
2835: }

2837: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2838: {
2839:   PetscInt  m = A->cmap->n, n = B->cmap->n;
2840:   PetscBool cisdense = PETSC_FALSE;

2842:   PetscFunctionBegin;
2843:   PetscCall(MatSetSizes(C, m, n, m, n));
2844: #if defined(PETSC_HAVE_CUDA)
2845:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2846: #endif
2847: #if defined(PETSC_HAVE_HIP)
2848:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2849: #endif
2850:   if (!cisdense) {
2851:     PetscBool flg;

2853:     PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2854:     PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2855:   }
2856:   PetscCall(MatSetUp(C));
2857:   PetscFunctionReturn(PETSC_SUCCESS);
2858: }

2860: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2861: {
2862:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data;
2863:   Mat_SeqDense      *b = (Mat_SeqDense *)B->data;
2864:   Mat_SeqDense      *c = (Mat_SeqDense *)C->data;
2865:   const PetscScalar *av, *bv;
2866:   PetscScalar       *cv;
2867:   PetscBLASInt       m, n, k;
2868:   PetscScalar        _DOne = 1.0, _DZero = 0.0;

2870:   PetscFunctionBegin;
2871:   PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2872:   PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2873:   PetscCall(PetscBLASIntCast(A->rmap->n, &k));
2874:   if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2875:   PetscCall(MatDenseGetArrayRead(A, &av));
2876:   PetscCall(MatDenseGetArrayRead(B, &bv));
2877:   PetscCall(MatDenseGetArrayWrite(C, &cv));
2878:   PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2879:   PetscCall(MatDenseRestoreArrayRead(A, &av));
2880:   PetscCall(MatDenseRestoreArrayRead(B, &bv));
2881:   PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2882:   PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2883:   PetscFunctionReturn(PETSC_SUCCESS);
2884: }

2886: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2887: {
2888:   PetscFunctionBegin;
2889:   C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2890:   C->ops->productsymbolic = MatProductSymbolic_AB;
2891:   PetscFunctionReturn(PETSC_SUCCESS);
2892: }

2894: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2895: {
2896:   PetscFunctionBegin;
2897:   C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2898:   C->ops->productsymbolic          = MatProductSymbolic_AtB;
2899:   PetscFunctionReturn(PETSC_SUCCESS);
2900: }

2902: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2903: {
2904:   PetscFunctionBegin;
2905:   C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2906:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
2907:   PetscFunctionReturn(PETSC_SUCCESS);
2908: }

2910: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2911: {
2912:   Mat_Product *product = C->product;

2914:   PetscFunctionBegin;
2915:   switch (product->type) {
2916:   case MATPRODUCT_AB:
2917:     PetscCall(MatProductSetFromOptions_SeqDense_AB(C));
2918:     break;
2919:   case MATPRODUCT_AtB:
2920:     PetscCall(MatProductSetFromOptions_SeqDense_AtB(C));
2921:     break;
2922:   case MATPRODUCT_ABt:
2923:     PetscCall(MatProductSetFromOptions_SeqDense_ABt(C));
2924:     break;
2925:   default:
2926:     break;
2927:   }
2928:   PetscFunctionReturn(PETSC_SUCCESS);
2929: }

2931: static PetscErrorCode MatGetRowMax_SeqDense(Mat A, Vec v, PetscInt idx[])
2932: {
2933:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data;
2934:   PetscInt           i, j, m = A->rmap->n, n = A->cmap->n, p;
2935:   PetscScalar       *x;
2936:   const PetscScalar *aa;

2938:   PetscFunctionBegin;
2939:   PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2940:   PetscCall(VecGetArray(v, &x));
2941:   PetscCall(VecGetLocalSize(v, &p));
2942:   PetscCall(MatDenseGetArrayRead(A, &aa));
2943:   PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2944:   for (i = 0; i < m; i++) {
2945:     x[i] = aa[i];
2946:     if (idx) idx[i] = 0;
2947:     for (j = 1; j < n; j++) {
2948:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i + a->lda * j])) {
2949:         x[i] = aa[i + a->lda * j];
2950:         if (idx) idx[i] = j;
2951:       }
2952:     }
2953:   }
2954:   PetscCall(MatDenseRestoreArrayRead(A, &aa));
2955:   PetscCall(VecRestoreArray(v, &x));
2956:   PetscFunctionReturn(PETSC_SUCCESS);
2957: }

2959: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A, Vec v, PetscInt idx[])
2960: {
2961:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data;
2962:   PetscInt           i, j, m = A->rmap->n, n = A->cmap->n, p;
2963:   PetscScalar       *x;
2964:   PetscReal          atmp;
2965:   const PetscScalar *aa;

2967:   PetscFunctionBegin;
2968:   PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2969:   PetscCall(VecGetArray(v, &x));
2970:   PetscCall(VecGetLocalSize(v, &p));
2971:   PetscCall(MatDenseGetArrayRead(A, &aa));
2972:   PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2973:   for (i = 0; i < m; i++) {
2974:     x[i] = PetscAbsScalar(aa[i]);
2975:     for (j = 1; j < n; j++) {
2976:       atmp = PetscAbsScalar(aa[i + a->lda * j]);
2977:       if (PetscAbsScalar(x[i]) < atmp) {
2978:         x[i] = atmp;
2979:         if (idx) idx[i] = j;
2980:       }
2981:     }
2982:   }
2983:   PetscCall(MatDenseRestoreArrayRead(A, &aa));
2984:   PetscCall(VecRestoreArray(v, &x));
2985:   PetscFunctionReturn(PETSC_SUCCESS);
2986: }

2988: static PetscErrorCode MatGetRowMin_SeqDense(Mat A, Vec v, PetscInt idx[])
2989: {
2990:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data;
2991:   PetscInt           i, j, m = A->rmap->n, n = A->cmap->n, p;
2992:   PetscScalar       *x;
2993:   const PetscScalar *aa;

2995:   PetscFunctionBegin;
2996:   PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2997:   PetscCall(MatDenseGetArrayRead(A, &aa));
2998:   PetscCall(VecGetArray(v, &x));
2999:   PetscCall(VecGetLocalSize(v, &p));
3000:   PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
3001:   for (i = 0; i < m; i++) {
3002:     x[i] = aa[i];
3003:     if (idx) idx[i] = 0;
3004:     for (j = 1; j < n; j++) {
3005:       if (PetscRealPart(x[i]) > PetscRealPart(aa[i + a->lda * j])) {
3006:         x[i] = aa[i + a->lda * j];
3007:         if (idx) idx[i] = j;
3008:       }
3009:     }
3010:   }
3011:   PetscCall(VecRestoreArray(v, &x));
3012:   PetscCall(MatDenseRestoreArrayRead(A, &aa));
3013:   PetscFunctionReturn(PETSC_SUCCESS);
3014: }

3016: PetscErrorCode MatGetColumnVector_SeqDense(Mat A, Vec v, PetscInt col)
3017: {
3018:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data;
3019:   PetscScalar       *x;
3020:   const PetscScalar *aa;

3022:   PetscFunctionBegin;
3023:   PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3024:   PetscCall(MatDenseGetArrayRead(A, &aa));
3025:   PetscCall(VecGetArray(v, &x));
3026:   PetscCall(PetscArraycpy(x, aa + col * a->lda, A->rmap->n));
3027:   PetscCall(VecRestoreArray(v, &x));
3028:   PetscCall(MatDenseRestoreArrayRead(A, &aa));
3029:   PetscFunctionReturn(PETSC_SUCCESS);
3030: }

3032: PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A, PetscInt type, PetscReal *reductions)
3033: {
3034:   PetscInt           i, j, m, n;
3035:   const PetscScalar *a;

3037:   PetscFunctionBegin;
3038:   PetscCall(MatGetSize(A, &m, &n));
3039:   PetscCall(PetscArrayzero(reductions, n));
3040:   PetscCall(MatDenseGetArrayRead(A, &a));
3041:   if (type == NORM_2) {
3042:     for (i = 0; i < n; i++) {
3043:       for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j] * a[j]);
3044:       a = PetscSafePointerPlusOffset(a, m);
3045:     }
3046:   } else if (type == NORM_1) {
3047:     for (i = 0; i < n; i++) {
3048:       for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j]);
3049:       a = PetscSafePointerPlusOffset(a, m);
3050:     }
3051:   } else if (type == NORM_INFINITY) {
3052:     for (i = 0; i < n; i++) {
3053:       for (j = 0; j < m; j++) reductions[i] = PetscMax(PetscAbsScalar(a[j]), reductions[i]);
3054:       a = PetscSafePointerPlusOffset(a, m);
3055:     }
3056:   } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
3057:     for (i = 0; i < n; i++) {
3058:       for (j = 0; j < m; j++) reductions[i] += PetscRealPart(a[j]);
3059:       a = PetscSafePointerPlusOffset(a, m);
3060:     }
3061:   } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3062:     for (i = 0; i < n; i++) {
3063:       for (j = 0; j < m; j++) reductions[i] += PetscImaginaryPart(a[j]);
3064:       a = PetscSafePointerPlusOffset(a, m);
3065:     }
3066:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown reduction type");
3067:   PetscCall(MatDenseRestoreArrayRead(A, &a));
3068:   if (type == NORM_2) {
3069:     for (i = 0; i < n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
3070:   } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3071:     for (i = 0; i < n; i++) reductions[i] /= m;
3072:   }
3073:   PetscFunctionReturn(PETSC_SUCCESS);
3074: }

3076: PetscErrorCode MatSetRandom_SeqDense(Mat x, PetscRandom rctx)
3077: {
3078:   PetscScalar *a;
3079:   PetscInt     lda, m, n, i, j;

3081:   PetscFunctionBegin;
3082:   PetscCall(MatGetSize(x, &m, &n));
3083:   PetscCall(MatDenseGetLDA(x, &lda));
3084:   PetscCall(MatDenseGetArrayWrite(x, &a));
3085:   for (j = 0; j < n; j++) {
3086:     for (i = 0; i < m; i++) PetscCall(PetscRandomGetValue(rctx, a + j * lda + i));
3087:   }
3088:   PetscCall(MatDenseRestoreArrayWrite(x, &a));
3089:   PetscFunctionReturn(PETSC_SUCCESS);
3090: }

3092: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A, PetscBool *missing, PetscInt *d)
3093: {
3094:   PetscFunctionBegin;
3095:   *missing = PETSC_FALSE;
3096:   PetscFunctionReturn(PETSC_SUCCESS);
3097: }

3099: /* vals is not const */
3100: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A, PetscInt col, PetscScalar **vals)
3101: {
3102:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3103:   PetscScalar  *v;

3105:   PetscFunctionBegin;
3106:   PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3107:   PetscCall(MatDenseGetArray(A, &v));
3108:   *vals = v + col * a->lda;
3109:   PetscCall(MatDenseRestoreArray(A, &v));
3110:   PetscFunctionReturn(PETSC_SUCCESS);
3111: }

3113: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A, PetscScalar **vals)
3114: {
3115:   PetscFunctionBegin;
3116:   if (vals) *vals = NULL; /* user cannot accidentally use the array later */
3117:   PetscFunctionReturn(PETSC_SUCCESS);
3118: }

3120: static struct _MatOps MatOps_Values = {MatSetValues_SeqDense,
3121:                                        MatGetRow_SeqDense,
3122:                                        MatRestoreRow_SeqDense,
3123:                                        MatMult_SeqDense,
3124:                                        /*  4*/ MatMultAdd_SeqDense,
3125:                                        MatMultTranspose_SeqDense,
3126:                                        MatMultTransposeAdd_SeqDense,
3127:                                        NULL,
3128:                                        NULL,
3129:                                        NULL,
3130:                                        /* 10*/ NULL,
3131:                                        MatLUFactor_SeqDense,
3132:                                        MatCholeskyFactor_SeqDense,
3133:                                        MatSOR_SeqDense,
3134:                                        MatTranspose_SeqDense,
3135:                                        /* 15*/ MatGetInfo_SeqDense,
3136:                                        MatEqual_SeqDense,
3137:                                        MatGetDiagonal_SeqDense,
3138:                                        MatDiagonalScale_SeqDense,
3139:                                        MatNorm_SeqDense,
3140:                                        /* 20*/ MatAssemblyBegin_SeqDense,
3141:                                        MatAssemblyEnd_SeqDense,
3142:                                        MatSetOption_SeqDense,
3143:                                        MatZeroEntries_SeqDense,
3144:                                        /* 24*/ MatZeroRows_SeqDense,
3145:                                        NULL,
3146:                                        NULL,
3147:                                        NULL,
3148:                                        NULL,
3149:                                        /* 29*/ MatSetUp_SeqDense,
3150:                                        NULL,
3151:                                        NULL,
3152:                                        NULL,
3153:                                        NULL,
3154:                                        /* 34*/ MatDuplicate_SeqDense,
3155:                                        NULL,
3156:                                        NULL,
3157:                                        NULL,
3158:                                        NULL,
3159:                                        /* 39*/ MatAXPY_SeqDense,
3160:                                        MatCreateSubMatrices_SeqDense,
3161:                                        NULL,
3162:                                        MatGetValues_SeqDense,
3163:                                        MatCopy_SeqDense,
3164:                                        /* 44*/ MatGetRowMax_SeqDense,
3165:                                        MatScale_SeqDense,
3166:                                        MatShift_SeqDense,
3167:                                        NULL,
3168:                                        MatZeroRowsColumns_SeqDense,
3169:                                        /* 49*/ MatSetRandom_SeqDense,
3170:                                        NULL,
3171:                                        NULL,
3172:                                        NULL,
3173:                                        NULL,
3174:                                        /* 54*/ NULL,
3175:                                        NULL,
3176:                                        NULL,
3177:                                        NULL,
3178:                                        NULL,
3179:                                        /* 59*/ MatCreateSubMatrix_SeqDense,
3180:                                        MatDestroy_SeqDense,
3181:                                        MatView_SeqDense,
3182:                                        NULL,
3183:                                        NULL,
3184:                                        /* 64*/ NULL,
3185:                                        NULL,
3186:                                        NULL,
3187:                                        NULL,
3188:                                        NULL,
3189:                                        /* 69*/ MatGetRowMaxAbs_SeqDense,
3190:                                        NULL,
3191:                                        NULL,
3192:                                        NULL,
3193:                                        NULL,
3194:                                        /* 74*/ NULL,
3195:                                        NULL,
3196:                                        NULL,
3197:                                        NULL,
3198:                                        NULL,
3199:                                        /* 79*/ NULL,
3200:                                        NULL,
3201:                                        NULL,
3202:                                        NULL,
3203:                                        /* 83*/ MatLoad_SeqDense,
3204:                                        MatIsSymmetric_SeqDense,
3205:                                        MatIsHermitian_SeqDense,
3206:                                        NULL,
3207:                                        NULL,
3208:                                        NULL,
3209:                                        /* 89*/ NULL,
3210:                                        NULL,
3211:                                        MatMatMultNumeric_SeqDense_SeqDense,
3212:                                        NULL,
3213:                                        NULL,
3214:                                        /* 94*/ NULL,
3215:                                        NULL,
3216:                                        NULL,
3217:                                        MatMatTransposeMultNumeric_SeqDense_SeqDense,
3218:                                        NULL,
3219:                                        /* 99*/ MatProductSetFromOptions_SeqDense,
3220:                                        NULL,
3221:                                        NULL,
3222:                                        MatConjugate_SeqDense,
3223:                                        NULL,
3224:                                        /*104*/ NULL,
3225:                                        MatRealPart_SeqDense,
3226:                                        MatImaginaryPart_SeqDense,
3227:                                        NULL,
3228:                                        NULL,
3229:                                        /*109*/ NULL,
3230:                                        NULL,
3231:                                        MatGetRowMin_SeqDense,
3232:                                        MatGetColumnVector_SeqDense,
3233:                                        MatMissingDiagonal_SeqDense,
3234:                                        /*114*/ NULL,
3235:                                        NULL,
3236:                                        NULL,
3237:                                        NULL,
3238:                                        NULL,
3239:                                        /*119*/ NULL,
3240:                                        NULL,
3241:                                        MatMultHermitianTranspose_SeqDense,
3242:                                        MatMultHermitianTransposeAdd_SeqDense,
3243:                                        NULL,
3244:                                        /*124*/ NULL,
3245:                                        MatGetColumnReductions_SeqDense,
3246:                                        NULL,
3247:                                        NULL,
3248:                                        NULL,
3249:                                        /*129*/ NULL,
3250:                                        NULL,
3251:                                        NULL,
3252:                                        MatTransposeMatMultNumeric_SeqDense_SeqDense,
3253:                                        NULL,
3254:                                        /*134*/ NULL,
3255:                                        NULL,
3256:                                        NULL,
3257:                                        NULL,
3258:                                        NULL,
3259:                                        /*139*/ NULL,
3260:                                        NULL,
3261:                                        NULL,
3262:                                        NULL,
3263:                                        NULL,
3264:                                        MatCreateMPIMatConcatenateSeqMat_SeqDense,
3265:                                        /*145*/ NULL,
3266:                                        NULL,
3267:                                        NULL,
3268:                                        NULL,
3269:                                        NULL,
3270:                                        /*150*/ NULL,
3271:                                        NULL,
3272:                                        NULL};

3274: /*@C
3275:   MatCreateSeqDense - Creates a `MATSEQDENSE` that
3276:   is stored in column major order (the usual Fortran format).

3278:   Collective

3280:   Input Parameters:
3281: + comm - MPI communicator, set to `PETSC_COMM_SELF`
3282: . m    - number of rows
3283: . n    - number of columns
3284: - data - optional location of matrix data in column major order.  Use `NULL` for PETSc
3285:    to control all matrix memory allocation.

3287:   Output Parameter:
3288: . A - the matrix

3290:   Level: intermediate

3292:   Note:
3293:   The data input variable is intended primarily for Fortran programmers
3294:   who wish to allocate their own matrix memory space.  Most users should
3295:   set `data` = `NULL`.

3297:   Developer Note:
3298:   Many of the matrix operations for this variant use the BLAS and LAPACK routines.

3300: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`
3301: @*/
3302: PetscErrorCode MatCreateSeqDense(MPI_Comm comm, PetscInt m, PetscInt n, PetscScalar *data, Mat *A)
3303: {
3304:   PetscFunctionBegin;
3305:   PetscCall(MatCreate(comm, A));
3306:   PetscCall(MatSetSizes(*A, m, n, m, n));
3307:   PetscCall(MatSetType(*A, MATSEQDENSE));
3308:   PetscCall(MatSeqDenseSetPreallocation(*A, data));
3309:   PetscFunctionReturn(PETSC_SUCCESS);
3310: }

3312: /*@C
3313:   MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements of a `MATSEQDENSE` matrix

3315:   Collective

3317:   Input Parameters:
3318: + B    - the matrix
3319: - data - the array (or `NULL`)

3321:   Level: intermediate

3323:   Note:
3324:   The data input variable is intended primarily for Fortran programmers
3325:   who wish to allocate their own matrix memory space.  Most users should
3326:   need not call this routine.

3328: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`, `MatDenseSetLDA()`
3329: @*/
3330: PetscErrorCode MatSeqDenseSetPreallocation(Mat B, PetscScalar data[])
3331: {
3332:   PetscFunctionBegin;
3334:   PetscTryMethod(B, "MatSeqDenseSetPreallocation_C", (Mat, PetscScalar[]), (B, data));
3335:   PetscFunctionReturn(PETSC_SUCCESS);
3336: }

3338: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B, PetscScalar *data)
3339: {
3340:   Mat_SeqDense *b = (Mat_SeqDense *)B->data;

3342:   PetscFunctionBegin;
3343:   PetscCheck(!b->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3344:   B->preallocated = PETSC_TRUE;

3346:   PetscCall(PetscLayoutSetUp(B->rmap));
3347:   PetscCall(PetscLayoutSetUp(B->cmap));

3349:   if (b->lda <= 0) b->lda = B->rmap->n;

3351:   if (!data) { /* petsc-allocated storage */
3352:     if (!b->user_alloc) PetscCall(PetscFree(b->v));
3353:     PetscCall(PetscCalloc1((size_t)b->lda * B->cmap->n, &b->v));

3355:     b->user_alloc = PETSC_FALSE;
3356:   } else { /* user-allocated storage */
3357:     if (!b->user_alloc) PetscCall(PetscFree(b->v));
3358:     b->v          = data;
3359:     b->user_alloc = PETSC_TRUE;
3360:   }
3361:   B->assembled = PETSC_TRUE;
3362:   PetscFunctionReturn(PETSC_SUCCESS);
3363: }

3365: #if defined(PETSC_HAVE_ELEMENTAL)
3366: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
3367: {
3368:   Mat                mat_elemental;
3369:   const PetscScalar *array;
3370:   PetscScalar       *v_colwise;
3371:   PetscInt           M = A->rmap->N, N = A->cmap->N, i, j, k, *rows, *cols;

3373:   PetscFunctionBegin;
3374:   PetscCall(PetscMalloc3(M * N, &v_colwise, M, &rows, N, &cols));
3375:   PetscCall(MatDenseGetArrayRead(A, &array));
3376:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
3377:   k = 0;
3378:   for (j = 0; j < N; j++) {
3379:     cols[j] = j;
3380:     for (i = 0; i < M; i++) v_colwise[j * M + i] = array[k++];
3381:   }
3382:   for (i = 0; i < M; i++) rows[i] = i;
3383:   PetscCall(MatDenseRestoreArrayRead(A, &array));

3385:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
3386:   PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
3387:   PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
3388:   PetscCall(MatSetUp(mat_elemental));

3390:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
3391:   PetscCall(MatSetValues(mat_elemental, M, rows, N, cols, v_colwise, ADD_VALUES));
3392:   PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
3393:   PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
3394:   PetscCall(PetscFree3(v_colwise, rows, cols));

3396:   if (reuse == MAT_INPLACE_MATRIX) {
3397:     PetscCall(MatHeaderReplace(A, &mat_elemental));
3398:   } else {
3399:     *newmat = mat_elemental;
3400:   }
3401:   PetscFunctionReturn(PETSC_SUCCESS);
3402: }
3403: #endif

3405: PetscErrorCode MatDenseSetLDA_SeqDense(Mat B, PetscInt lda)
3406: {
3407:   Mat_SeqDense *b = (Mat_SeqDense *)B->data;
3408:   PetscBool     data;

3410:   PetscFunctionBegin;
3411:   data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
3412:   PetscCheck(b->user_alloc || !data || b->lda == lda, PETSC_COMM_SELF, PETSC_ERR_ORDER, "LDA cannot be changed after allocation of internal storage");
3413:   PetscCheck(lda >= B->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "LDA %" PetscInt_FMT " must be at least matrix dimension %" PetscInt_FMT, lda, B->rmap->n);
3414:   b->lda = lda;
3415:   PetscFunctionReturn(PETSC_SUCCESS);
3416: }

3418: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
3419: {
3420:   PetscFunctionBegin;
3421:   PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIDense(comm, inmat, n, scall, outmat));
3422:   PetscFunctionReturn(PETSC_SUCCESS);
3423: }

3425: PetscErrorCode MatDenseCreateColumnVec_Private(Mat A, Vec *v)
3426: {
3427:   PetscBool   isstd, iskok, iscuda, iship;
3428:   PetscMPIInt size;
3429: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
3430:   /* we pass the data of A, to prevent allocating needless GPU memory the first time VecCUPMPlaceArray is called. */
3431:   const PetscScalar *a;
3432: #endif

3434:   PetscFunctionBegin;
3435:   *v = NULL;
3436:   PetscCall(PetscStrcmpAny(A->defaultvectype, &isstd, VECSTANDARD, VECSEQ, VECMPI, ""));
3437:   PetscCall(PetscStrcmpAny(A->defaultvectype, &iskok, VECKOKKOS, VECSEQKOKKOS, VECMPIKOKKOS, ""));
3438:   PetscCall(PetscStrcmpAny(A->defaultvectype, &iscuda, VECCUDA, VECSEQCUDA, VECMPICUDA, ""));
3439:   PetscCall(PetscStrcmpAny(A->defaultvectype, &iship, VECHIP, VECSEQHIP, VECMPIHIP, ""));
3440:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
3441:   if (isstd) {
3442:     if (size > 1) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, v));
3443:     else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, v));
3444:   } else if (iskok) {
3445:     PetscCheck(PetscDefined(HAVE_KOKKOS_KERNELS), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using KOKKOS kernels support");
3446: #if PetscDefined(HAVE_KOKKOS_KERNELS)
3447:     if (size > 1) PetscCall(VecCreateMPIKokkosWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, v));
3448:     else PetscCall(VecCreateSeqKokkosWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, v));
3449: #endif
3450:   } else if (iscuda) {
3451:     PetscCheck(PetscDefined(HAVE_CUDA), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using CUDA support");
3452: #if PetscDefined(HAVE_CUDA)
3453:     PetscCall(MatDenseCUDAGetArrayRead(A, &a));
3454:     if (size > 1) PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, a, v));
3455:     else PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, a, v));
3456: #endif
3457:   } else if (iship) {
3458:     PetscCheck(PetscDefined(HAVE_HIP), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using HIP support");
3459: #if PetscDefined(HAVE_HIP)
3460:     PetscCall(MatDenseHIPGetArrayRead(A, &a));
3461:     if (size > 1) PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, a, v));
3462:     else PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, a, v));
3463: #endif
3464:   }
3465:   PetscCheck(*v, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not coded for type %s", A->defaultvectype);
3466:   PetscFunctionReturn(PETSC_SUCCESS);
3467: }

3469: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3470: {
3471:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3473:   PetscFunctionBegin;
3474:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3475:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3476:   if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3477:   a->vecinuse = col + 1;
3478:   PetscCall(MatDenseGetArray(A, (PetscScalar **)&a->ptrinuse));
3479:   PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda));
3480:   *v = a->cvec;
3481:   PetscFunctionReturn(PETSC_SUCCESS);
3482: }

3484: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3485: {
3486:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3488:   PetscFunctionBegin;
3489:   PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3490:   PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3491:   VecCheckAssembled(a->cvec);
3492:   a->vecinuse = 0;
3493:   PetscCall(MatDenseRestoreArray(A, (PetscScalar **)&a->ptrinuse));
3494:   PetscCall(VecResetArray(a->cvec));
3495:   if (v) *v = NULL;
3496:   PetscFunctionReturn(PETSC_SUCCESS);
3497: }

3499: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3500: {
3501:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3503:   PetscFunctionBegin;
3504:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3505:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3506:   if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3507:   a->vecinuse = col + 1;
3508:   PetscCall(MatDenseGetArrayRead(A, &a->ptrinuse));
3509:   PetscCall(VecPlaceArray(a->cvec, PetscSafePointerPlusOffset(a->ptrinuse, (size_t)col * (size_t)a->lda)));
3510:   PetscCall(VecLockReadPush(a->cvec));
3511:   *v = a->cvec;
3512:   PetscFunctionReturn(PETSC_SUCCESS);
3513: }

3515: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3516: {
3517:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3519:   PetscFunctionBegin;
3520:   PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3521:   PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3522:   VecCheckAssembled(a->cvec);
3523:   a->vecinuse = 0;
3524:   PetscCall(MatDenseRestoreArrayRead(A, &a->ptrinuse));
3525:   PetscCall(VecLockReadPop(a->cvec));
3526:   PetscCall(VecResetArray(a->cvec));
3527:   if (v) *v = NULL;
3528:   PetscFunctionReturn(PETSC_SUCCESS);
3529: }

3531: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3532: {
3533:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3535:   PetscFunctionBegin;
3536:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3537:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3538:   if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3539:   a->vecinuse = col + 1;
3540:   PetscCall(MatDenseGetArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3541:   PetscCall(VecPlaceArray(a->cvec, PetscSafePointerPlusOffset(a->ptrinuse, (size_t)col * (size_t)a->lda)));
3542:   *v = a->cvec;
3543:   PetscFunctionReturn(PETSC_SUCCESS);
3544: }

3546: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3547: {
3548:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3550:   PetscFunctionBegin;
3551:   PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3552:   PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3553:   VecCheckAssembled(a->cvec);
3554:   a->vecinuse = 0;
3555:   PetscCall(MatDenseRestoreArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3556:   PetscCall(VecResetArray(a->cvec));
3557:   if (v) *v = NULL;
3558:   PetscFunctionReturn(PETSC_SUCCESS);
3559: }

3561: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3562: {
3563:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3565:   PetscFunctionBegin;
3566:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3567:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3568:   if (a->cmat && (cend - cbegin != a->cmat->cmap->N || rend - rbegin != a->cmat->rmap->N)) PetscCall(MatDestroy(&a->cmat));
3569:   if (!a->cmat) {
3570:     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), rend - rbegin, PETSC_DECIDE, rend - rbegin, cend - cbegin, PetscSafePointerPlusOffset(a->v, rbegin + (size_t)cbegin * a->lda), &a->cmat));
3571:   } else {
3572:     PetscCall(MatDensePlaceArray(a->cmat, PetscSafePointerPlusOffset(a->v, rbegin + (size_t)cbegin * a->lda)));
3573:   }
3574:   PetscCall(MatDenseSetLDA(a->cmat, a->lda));
3575:   a->matinuse = cbegin + 1;
3576:   *v          = a->cmat;
3577: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
3578:   A->offloadmask = PETSC_OFFLOAD_CPU;
3579: #endif
3580:   PetscFunctionReturn(PETSC_SUCCESS);
3581: }

3583: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A, Mat *v)
3584: {
3585:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3587:   PetscFunctionBegin;
3588:   PetscCheck(a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetSubMatrix() first");
3589:   PetscCheck(a->cmat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column matrix");
3590:   PetscCheck(*v == a->cmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not the matrix obtained from MatDenseGetSubMatrix()");
3591:   a->matinuse = 0;
3592:   PetscCall(MatDenseResetArray(a->cmat));
3593:   if (v) *v = NULL;
3594: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
3595:   A->offloadmask = PETSC_OFFLOAD_CPU;
3596: #endif
3597:   PetscFunctionReturn(PETSC_SUCCESS);
3598: }

3600: /*MC
3601:    MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.

3603:    Options Database Key:
3604: . -mat_type seqdense - sets the matrix type to `MATSEQDENSE` during a call to `MatSetFromOptions()`

3606:   Level: beginner

3608: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreateSeqDense()`
3609: M*/
3610: PetscErrorCode MatCreate_SeqDense(Mat B)
3611: {
3612:   Mat_SeqDense *b;
3613:   PetscMPIInt   size;

3615:   PetscFunctionBegin;
3616:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
3617:   PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1");

3619:   PetscCall(PetscNew(&b));
3620:   B->data   = (void *)b;
3621:   B->ops[0] = MatOps_Values;

3623:   b->roworiented = PETSC_TRUE;

3625:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactor_C", MatQRFactor_SeqDense));
3626:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetLDA_C", MatDenseGetLDA_SeqDense));
3627:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseSetLDA_C", MatDenseSetLDA_SeqDense));
3628:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArray_C", MatDenseGetArray_SeqDense));
3629:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArray_C", MatDenseRestoreArray_SeqDense));
3630:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDensePlaceArray_C", MatDensePlaceArray_SeqDense));
3631:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseResetArray_C", MatDenseResetArray_SeqDense));
3632:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseReplaceArray_C", MatDenseReplaceArray_SeqDense));
3633:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayRead_C", MatDenseGetArray_SeqDense));
3634:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayRead_C", MatDenseRestoreArray_SeqDense));
3635:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayWrite_C", MatDenseGetArray_SeqDense));
3636:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayWrite_C", MatDenseRestoreArray_SeqDense));
3637:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqaij_C", MatConvert_SeqDense_SeqAIJ));
3638: #if defined(PETSC_HAVE_ELEMENTAL)
3639:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_elemental_C", MatConvert_SeqDense_Elemental));
3640: #endif
3641: #if defined(PETSC_HAVE_SCALAPACK)
3642:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_scalapack_C", MatConvert_Dense_ScaLAPACK));
3643: #endif
3644: #if defined(PETSC_HAVE_CUDA)
3645:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensecuda_C", MatConvert_SeqDense_SeqDenseCUDA));
3646:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3647:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdense_C", MatProductSetFromOptions_SeqDense));
3648:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3649: #endif
3650: #if defined(PETSC_HAVE_HIP)
3651:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensehip_C", MatConvert_SeqDense_SeqDenseHIP));
3652:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3653:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdense_C", MatProductSetFromOptions_SeqDense));
3654:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3655: #endif
3656:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqDenseSetPreallocation_C", MatSeqDenseSetPreallocation_SeqDense));
3657:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqdense_C", MatProductSetFromOptions_SeqAIJ_SeqDense));
3658:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdense_C", MatProductSetFromOptions_SeqDense));
3659:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));
3660:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqsbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));

3662:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumn_C", MatDenseGetColumn_SeqDense));
3663:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumn_C", MatDenseRestoreColumn_SeqDense));
3664:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVec_C", MatDenseGetColumnVec_SeqDense));
3665:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVec_C", MatDenseRestoreColumnVec_SeqDense));
3666:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecRead_C", MatDenseGetColumnVecRead_SeqDense));
3667:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecRead_C", MatDenseRestoreColumnVecRead_SeqDense));
3668:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecWrite_C", MatDenseGetColumnVecWrite_SeqDense));
3669:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecWrite_C", MatDenseRestoreColumnVecWrite_SeqDense));
3670:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetSubMatrix_C", MatDenseGetSubMatrix_SeqDense));
3671:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreSubMatrix_C", MatDenseRestoreSubMatrix_SeqDense));
3672:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultAddColumnRange_C", MatMultAddColumnRange_SeqDense));
3673:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultHermitianTransposeColumnRange_C", MatMultHermitianTransposeColumnRange_SeqDense));
3674:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultHermitianTransposeAddColumnRange_C", MatMultHermitianTransposeAddColumnRange_SeqDense));
3675:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQDENSE));
3676:   PetscFunctionReturn(PETSC_SUCCESS);
3677: }

3679: /*@C
3680:   MatDenseGetColumn - gives access to a column of a dense matrix. This is only the local part of the column. You MUST call `MatDenseRestoreColumn()` to avoid memory bleeding.

3682:   Not Collective

3684:   Input Parameters:
3685: + A   - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3686: - col - column index

3688:   Output Parameter:
3689: . vals - pointer to the data

3691:   Level: intermediate

3693:   Note:
3694:   Use `MatDenseGetColumnVec()` to get access to a column of a `MATDENSE` treated as a `Vec`

3696: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreColumn()`, `MatDenseGetColumnVec()`
3697: @*/
3698: PetscErrorCode MatDenseGetColumn(Mat A, PetscInt col, PetscScalar **vals)
3699: {
3700:   PetscFunctionBegin;
3703:   PetscAssertPointer(vals, 3);
3704:   PetscUseMethod(A, "MatDenseGetColumn_C", (Mat, PetscInt, PetscScalar **), (A, col, vals));
3705:   PetscFunctionReturn(PETSC_SUCCESS);
3706: }

3708: /*@C
3709:   MatDenseRestoreColumn - returns access to a column of a `MATDENSE` matrix which is returned by `MatDenseGetColumn()`.

3711:   Not Collective

3713:   Input Parameters:
3714: + A    - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3715: - vals - pointer to the data (may be `NULL`)

3717:   Level: intermediate

3719: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetColumn()`
3720: @*/
3721: PetscErrorCode MatDenseRestoreColumn(Mat A, PetscScalar **vals)
3722: {
3723:   PetscFunctionBegin;
3725:   PetscAssertPointer(vals, 2);
3726:   PetscUseMethod(A, "MatDenseRestoreColumn_C", (Mat, PetscScalar **), (A, vals));
3727:   PetscFunctionReturn(PETSC_SUCCESS);
3728: }

3730: /*@
3731:   MatDenseGetColumnVec - Gives read-write access to a column of a `MATDENSE` matrix, represented as a `Vec`.

3733:   Collective

3735:   Input Parameters:
3736: + A   - the `Mat` object
3737: - col - the column index

3739:   Output Parameter:
3740: . v - the vector

3742:   Level: intermediate

3744:   Notes:
3745:   The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVec()` when the vector is no longer needed.

3747:   Use `MatDenseGetColumnVecRead()` to obtain read-only access or `MatDenseGetColumnVecWrite()` for write-only access.

3749: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`, `MatDenseGetColumn()`
3750: @*/
3751: PetscErrorCode MatDenseGetColumnVec(Mat A, PetscInt col, Vec *v)
3752: {
3753:   PetscFunctionBegin;
3757:   PetscAssertPointer(v, 3);
3758:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3759:   PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N);
3760:   PetscUseMethod(A, "MatDenseGetColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3761:   PetscFunctionReturn(PETSC_SUCCESS);
3762: }

3764: /*@
3765:   MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVec()`.

3767:   Collective

3769:   Input Parameters:
3770: + A   - the `Mat` object
3771: . col - the column index
3772: - v   - the `Vec` object (may be `NULL`)

3774:   Level: intermediate

3776: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3777: @*/
3778: PetscErrorCode MatDenseRestoreColumnVec(Mat A, PetscInt col, Vec *v)
3779: {
3780:   PetscFunctionBegin;
3784:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3785:   PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N);
3786:   PetscUseMethod(A, "MatDenseRestoreColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3787:   PetscFunctionReturn(PETSC_SUCCESS);
3788: }

3790: /*@
3791:   MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a `Vec`.

3793:   Collective

3795:   Input Parameters:
3796: + A   - the `Mat` object
3797: - col - the column index

3799:   Output Parameter:
3800: . v - the vector

3802:   Level: intermediate

3804:   Notes:
3805:   The vector is owned by PETSc and users cannot modify it.

3807:   Users need to call `MatDenseRestoreColumnVecRead()` when the vector is no longer needed.

3809:   Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecWrite()` for write-only access.

3811: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3812: @*/
3813: PetscErrorCode MatDenseGetColumnVecRead(Mat A, PetscInt col, Vec *v)
3814: {
3815:   PetscFunctionBegin;
3819:   PetscAssertPointer(v, 3);
3820:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3821:   PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N);
3822:   PetscUseMethod(A, "MatDenseGetColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3823:   PetscFunctionReturn(PETSC_SUCCESS);
3824: }

3826: /*@
3827:   MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVecRead()`.

3829:   Collective

3831:   Input Parameters:
3832: + A   - the `Mat` object
3833: . col - the column index
3834: - v   - the `Vec` object (may be `NULL`)

3836:   Level: intermediate

3838: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecWrite()`
3839: @*/
3840: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A, PetscInt col, Vec *v)
3841: {
3842:   PetscFunctionBegin;
3846:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3847:   PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N);
3848:   PetscUseMethod(A, "MatDenseRestoreColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3849:   PetscFunctionReturn(PETSC_SUCCESS);
3850: }

3852: /*@
3853:   MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a `Vec`.

3855:   Collective

3857:   Input Parameters:
3858: + A   - the `Mat` object
3859: - col - the column index

3861:   Output Parameter:
3862: . v - the vector

3864:   Level: intermediate

3866:   Notes:
3867:   The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVecWrite()` when the vector is no longer needed.

3869:   Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecRead()` for read-only access.

3871: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3872: @*/
3873: PetscErrorCode MatDenseGetColumnVecWrite(Mat A, PetscInt col, Vec *v)
3874: {
3875:   PetscFunctionBegin;
3879:   PetscAssertPointer(v, 3);
3880:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3881:   PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N);
3882:   PetscUseMethod(A, "MatDenseGetColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3883:   PetscFunctionReturn(PETSC_SUCCESS);
3884: }

3886: /*@
3887:   MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVecWrite()`.

3889:   Collective

3891:   Input Parameters:
3892: + A   - the `Mat` object
3893: . col - the column index
3894: - v   - the `Vec` object (may be `NULL`)

3896:   Level: intermediate

3898: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`
3899: @*/
3900: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A, PetscInt col, Vec *v)
3901: {
3902:   PetscFunctionBegin;
3906:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3907:   PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N);
3908:   PetscUseMethod(A, "MatDenseRestoreColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3909:   PetscFunctionReturn(PETSC_SUCCESS);
3910: }

3912: /*@
3913:   MatDenseGetSubMatrix - Gives access to a block of rows and columns of a dense matrix, represented as a `Mat`.

3915:   Collective

3917:   Input Parameters:
3918: + A      - the `Mat` object
3919: . rbegin - the first global row index in the block (if `PETSC_DECIDE`, is 0)
3920: . rend   - the global row index past the last one in the block (if `PETSC_DECIDE`, is `M`)
3921: . cbegin - the first global column index in the block (if `PETSC_DECIDE`, is 0)
3922: - cend   - the global column index past the last one in the block (if `PETSC_DECIDE`, is `N`)

3924:   Output Parameter:
3925: . v - the matrix

3927:   Level: intermediate

3929:   Notes:
3930:   The matrix is owned by PETSc. Users need to call `MatDenseRestoreSubMatrix()` when the matrix is no longer needed.

3932:   The output matrix is not redistributed by PETSc, so depending on the values of `rbegin` and `rend`, some processes may have no local rows.

3934: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreSubMatrix()`
3935: @*/
3936: PetscErrorCode MatDenseGetSubMatrix(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3937: {
3938:   PetscFunctionBegin;
3945:   PetscAssertPointer(v, 6);
3946:   if (rbegin == PETSC_DECIDE) rbegin = 0;
3947:   if (rend == PETSC_DECIDE) rend = A->rmap->N;
3948:   if (cbegin == PETSC_DECIDE) cbegin = 0;
3949:   if (cend == PETSC_DECIDE) cend = A->cmap->N;
3950:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3951:   PetscCheck(rbegin >= 0 && rbegin <= A->rmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid rbegin %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT "]", rbegin, A->rmap->N);
3952:   PetscCheck(rend >= rbegin && rend <= A->rmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid rend %" PetscInt_FMT ", should be in [%" PetscInt_FMT ",%" PetscInt_FMT "]", rend, rbegin, A->rmap->N);
3953:   PetscCheck(cbegin >= 0 && cbegin <= A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid cbegin %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT "]", cbegin, A->cmap->N);
3954:   PetscCheck(cend >= cbegin && cend <= A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid cend %" PetscInt_FMT ", should be in [%" PetscInt_FMT ",%" PetscInt_FMT "]", cend, cbegin, A->cmap->N);
3955:   PetscUseMethod(A, "MatDenseGetSubMatrix_C", (Mat, PetscInt, PetscInt, PetscInt, PetscInt, Mat *), (A, rbegin, rend, cbegin, cend, v));
3956:   PetscFunctionReturn(PETSC_SUCCESS);
3957: }

3959: /*@
3960:   MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from `MatDenseGetSubMatrix()`.

3962:   Collective

3964:   Input Parameters:
3965: + A - the `Mat` object
3966: - v - the `Mat` object (may be `NULL`)

3968:   Level: intermediate

3970: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseGetSubMatrix()`
3971: @*/
3972: PetscErrorCode MatDenseRestoreSubMatrix(Mat A, Mat *v)
3973: {
3974:   PetscFunctionBegin;
3977:   PetscAssertPointer(v, 2);
3978:   PetscUseMethod(A, "MatDenseRestoreSubMatrix_C", (Mat, Mat *), (A, v));
3979:   PetscFunctionReturn(PETSC_SUCCESS);
3980: }

3982: #include <petscblaslapack.h>
3983: #include <petsc/private/kernels/blockinvert.h>

3985: PetscErrorCode MatSeqDenseInvert(Mat A)
3986: {
3987:   PetscInt        m;
3988:   const PetscReal shift = 0.0;
3989:   PetscBool       allowzeropivot, zeropivotdetected = PETSC_FALSE;
3990:   PetscScalar    *values;

3992:   PetscFunctionBegin;
3994:   PetscCall(MatDenseGetArray(A, &values));
3995:   PetscCall(MatGetLocalSize(A, &m, NULL));
3996:   allowzeropivot = PetscNot(A->erroriffailure);
3997:   /* factor and invert each block */
3998:   switch (m) {
3999:   case 1:
4000:     values[0] = (PetscScalar)1.0 / (values[0] + shift);
4001:     break;
4002:   case 2:
4003:     PetscCall(PetscKernel_A_gets_inverse_A_2(values, shift, allowzeropivot, &zeropivotdetected));
4004:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4005:     break;
4006:   case 3:
4007:     PetscCall(PetscKernel_A_gets_inverse_A_3(values, shift, allowzeropivot, &zeropivotdetected));
4008:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4009:     break;
4010:   case 4:
4011:     PetscCall(PetscKernel_A_gets_inverse_A_4(values, shift, allowzeropivot, &zeropivotdetected));
4012:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4013:     break;
4014:   case 5: {
4015:     PetscScalar work[25];
4016:     PetscInt    ipvt[5];

4018:     PetscCall(PetscKernel_A_gets_inverse_A_5(values, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
4019:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4020:   } break;
4021:   case 6:
4022:     PetscCall(PetscKernel_A_gets_inverse_A_6(values, shift, allowzeropivot, &zeropivotdetected));
4023:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4024:     break;
4025:   case 7:
4026:     PetscCall(PetscKernel_A_gets_inverse_A_7(values, shift, allowzeropivot, &zeropivotdetected));
4027:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4028:     break;
4029:   default: {
4030:     PetscInt    *v_pivots, *IJ, j;
4031:     PetscScalar *v_work;

4033:     PetscCall(PetscMalloc3(m, &v_work, m, &v_pivots, m, &IJ));
4034:     for (j = 0; j < m; j++) IJ[j] = j;
4035:     PetscCall(PetscKernel_A_gets_inverse_A(m, values, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
4036:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4037:     PetscCall(PetscFree3(v_work, v_pivots, IJ));
4038:   }
4039:   }
4040:   PetscCall(MatDenseRestoreArray(A, &values));
4041:   PetscFunctionReturn(PETSC_SUCCESS);
4042: }