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>

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

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

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

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

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

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

 97:   PetscFunctionBegin;
 98:   if (PetscDefined(USE_DEBUG)) {
 99:     for (i = 0; i < N; i++) {
100:       PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
101:       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);
102:       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);
103:     }
104:   }
105:   if (!N) PetscFunctionReturn(PETSC_SUCCESS);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

415: static PetscErrorCode MatConjugate_SeqDense(Mat);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

938:     mat->lfwork = -1;
939:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
940:     PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, &dummy, &mat->lfwork, &info));
941:     PetscCall(PetscFPTrapPop());
942:     mat->lfwork = (PetscInt)PetscRealPart(dummy);
943:     PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
944:   }
945:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
946:   PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, mat->fwork, &mat->lfwork, &info));
947:   PetscCall(PetscFPTrapPop());
948:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to QR factorization %d", (int)info);
949:   // 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
950:   mat->rank = min;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1311:   PetscFunctionBegin;
1312:   PetscCall(MatDenseGetArrayRead(A, &vv));
1313:   /* row-oriented output */
1314:   for (i = 0; i < m; i++) {
1315:     if (indexm[i] < 0) {
1316:       v += n;
1317:       continue;
1318:     }
1319:     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);
1320:     for (j = 0; j < n; j++) {
1321:       if (indexn[j] < 0) {
1322:         v++;
1323:         continue;
1324:       }
1325:       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);
1326:       *v++ = vv[indexn[j] * mat->lda + indexm[i]];
1327:     }
1328:   }
1329:   PetscCall(MatDenseRestoreArrayRead(A, &vv));
1330:   PetscFunctionReturn(PETSC_SUCCESS);
1331: }

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

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

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

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

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

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

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

1392:   if (!skipHeader) {
1393:     PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT));
1394:     PetscCheck(header[0] == MAT_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file");
1395:     M = header[1];
1396:     N = header[2];
1397:     PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M);
1398:     PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N);
1399:     nz = header[3];
1400:     PetscCheck(nz == MATRIX_BINARY_FORMAT_DENSE || nz >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Unknown matrix format %" PetscInt_FMT " in file", nz);
1401:   } else {
1402:     PetscCall(MatGetSize(mat, &M, &N));
1403:     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");
1404:     nz = MATRIX_BINARY_FORMAT_DENSE;
1405:   }

1407:   /* setup global sizes if not set */
1408:   if (mat->rmap->N < 0) mat->rmap->N = M;
1409:   if (mat->cmap->N < 0) mat->cmap->N = N;
1410:   PetscCall(MatSetUp(mat));
1411:   /* check if global sizes are correct */
1412:   PetscCall(MatGetSize(mat, &rows, &cols));
1413:   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);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2067:   PetscFunctionBegin;
2068:   if (PetscDefined(USE_DEBUG)) {
2069:     for (i = 0; i < N; i++) {
2070:       PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
2071:       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);
2072:     }
2073:   }
2074:   if (!N) PetscFunctionReturn(PETSC_SUCCESS);

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

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

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

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

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

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

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

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

2133:   Not Collective

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

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

2141:   Level: intermediate

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

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

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

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

2164:   Level: intermediate

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

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

2179:   Logically Collective

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

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

2187:   Level: intermediate

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

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

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

2206:   Logically Collective

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

2212:   Level: intermediate

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

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

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

2235:   Not Collective

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

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

2243:   Level: intermediate

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

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

2259:   Not Collective

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

2265:   Level: intermediate

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

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

2281:   Not Collective

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

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

2289:   Level: intermediate

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

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

2305:   Not Collective

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

2311:   Level: intermediate

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

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

2331:   Logically Collective

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

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

2340:   Level: intermediate

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

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

2353:   PetscFunctionBegin;
2355:   PetscAssertPointer(array, 2);
2356:   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 */
2357:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2358:   if (isMPI) {
2359:     /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2360:     PetscCall(MatDenseGetArrayAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2361:   } else {
2362:     PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);

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

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

2378:   Logically Collective

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

2384:   Level: intermediate

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

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

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

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

2416:   Logically Collective

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

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

2425:   Level: intermediate

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

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

2438:   PetscFunctionBegin;
2440:   PetscAssertPointer(array, 2);
2441:   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 */
2442:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2443:   if (isMPI) { /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2444:     PetscCall(MatDenseGetArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2445:   } else {
2446:     PetscErrorCode (*fptr)(Mat, const PetscScalar **, PetscMemType *);

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

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

2462:   Logically Collective

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

2468:   Level: intermediate

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

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

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

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

2499:   Logically Collective

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

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

2508:   Level: intermediate

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

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

2521:   PetscFunctionBegin;
2523:   PetscAssertPointer(array, 2);
2524:   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 */
2525:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2526:   if (isMPI) {
2527:     PetscCall(MatDenseGetArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2528:   } else {
2529:     PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);

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

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

2545:   Logically Collective

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

2551:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3277:   Collective

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

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

3289:   Level: intermediate

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

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

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

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

3314:   Collective

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

3320:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3409:   PetscFunctionBegin;
3410:   data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
3411:   PetscCheck(b->user_alloc || !data || b->lda == lda, PETSC_COMM_SELF, PETSC_ERR_ORDER, "LDA cannot be changed after allocation of internal storage");
3412:   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);
3413:   b->lda = lda;
3414:   PetscFunctionReturn(PETSC_SUCCESS);
3415: }

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

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

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

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

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

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

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

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

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

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

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

3528: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3529: {
3530:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

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

3543: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3544: {
3545:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

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

3557: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3558: {
3559:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

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

3579: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A, Mat *v)
3580: {
3581:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

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

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

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

3602:   Level: beginner

3604: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreateSeqDense()`
3605: M*/
3606: PetscErrorCode MatCreate_SeqDense(Mat B)
3607: {
3608:   Mat_SeqDense *b;
3609:   PetscMPIInt   size;

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

3615:   PetscCall(PetscNew(&b));
3616:   B->data   = (void *)b;
3617:   B->ops[0] = MatOps_Values;

3619:   b->roworiented = PETSC_TRUE;

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

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

3675: /*@C
3676:   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.

3678:   Not Collective

3680:   Input Parameters:
3681: + A   - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3682: - col - column index

3684:   Output Parameter:
3685: . vals - pointer to the data

3687:   Level: intermediate

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

3692: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreColumn()`, `MatDenseGetColumnVec()`
3693: @*/
3694: PetscErrorCode MatDenseGetColumn(Mat A, PetscInt col, PetscScalar **vals)
3695: {
3696:   PetscFunctionBegin;
3699:   PetscAssertPointer(vals, 3);
3700:   PetscUseMethod(A, "MatDenseGetColumn_C", (Mat, PetscInt, PetscScalar **), (A, col, vals));
3701:   PetscFunctionReturn(PETSC_SUCCESS);
3702: }

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

3707:   Not Collective

3709:   Input Parameters:
3710: + A    - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3711: - vals - pointer to the data (may be `NULL`)

3713:   Level: intermediate

3715: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetColumn()`
3716: @*/
3717: PetscErrorCode MatDenseRestoreColumn(Mat A, PetscScalar **vals)
3718: {
3719:   PetscFunctionBegin;
3721:   PetscAssertPointer(vals, 2);
3722:   PetscUseMethod(A, "MatDenseRestoreColumn_C", (Mat, PetscScalar **), (A, vals));
3723:   PetscFunctionReturn(PETSC_SUCCESS);
3724: }

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

3729:   Collective

3731:   Input Parameters:
3732: + A   - the `Mat` object
3733: - col - the column index

3735:   Output Parameter:
3736: . v - the vector

3738:   Level: intermediate

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

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

3745: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`, `MatDenseGetColumn()`
3746: @*/
3747: PetscErrorCode MatDenseGetColumnVec(Mat A, PetscInt col, Vec *v)
3748: {
3749:   PetscFunctionBegin;
3753:   PetscAssertPointer(v, 3);
3754:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3755:   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);
3756:   PetscUseMethod(A, "MatDenseGetColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3757:   PetscFunctionReturn(PETSC_SUCCESS);
3758: }

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

3763:   Collective

3765:   Input Parameters:
3766: + A   - the `Mat` object
3767: . col - the column index
3768: - v   - the `Vec` object (may be `NULL`)

3770:   Level: intermediate

3772: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3773: @*/
3774: PetscErrorCode MatDenseRestoreColumnVec(Mat A, PetscInt col, Vec *v)
3775: {
3776:   PetscFunctionBegin;
3780:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3781:   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);
3782:   PetscUseMethod(A, "MatDenseRestoreColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3783:   PetscFunctionReturn(PETSC_SUCCESS);
3784: }

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

3789:   Collective

3791:   Input Parameters:
3792: + A   - the `Mat` object
3793: - col - the column index

3795:   Output Parameter:
3796: . v - the vector

3798:   Level: intermediate

3800:   Notes:
3801:   The vector is owned by PETSc and users cannot modify it.

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

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

3807: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3808: @*/
3809: PetscErrorCode MatDenseGetColumnVecRead(Mat A, PetscInt col, Vec *v)
3810: {
3811:   PetscFunctionBegin;
3815:   PetscAssertPointer(v, 3);
3816:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3817:   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);
3818:   PetscUseMethod(A, "MatDenseGetColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3819:   PetscFunctionReturn(PETSC_SUCCESS);
3820: }

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

3825:   Collective

3827:   Input Parameters:
3828: + A   - the `Mat` object
3829: . col - the column index
3830: - v   - the `Vec` object (may be `NULL`)

3832:   Level: intermediate

3834: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecWrite()`
3835: @*/
3836: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A, PetscInt col, Vec *v)
3837: {
3838:   PetscFunctionBegin;
3842:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3843:   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);
3844:   PetscUseMethod(A, "MatDenseRestoreColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3845:   PetscFunctionReturn(PETSC_SUCCESS);
3846: }

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

3851:   Collective

3853:   Input Parameters:
3854: + A   - the `Mat` object
3855: - col - the column index

3857:   Output Parameter:
3858: . v - the vector

3860:   Level: intermediate

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

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

3867: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3868: @*/
3869: PetscErrorCode MatDenseGetColumnVecWrite(Mat A, PetscInt col, Vec *v)
3870: {
3871:   PetscFunctionBegin;
3875:   PetscAssertPointer(v, 3);
3876:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3877:   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);
3878:   PetscUseMethod(A, "MatDenseGetColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3879:   PetscFunctionReturn(PETSC_SUCCESS);
3880: }

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

3885:   Collective

3887:   Input Parameters:
3888: + A   - the `Mat` object
3889: . col - the column index
3890: - v   - the `Vec` object (may be `NULL`)

3892:   Level: intermediate

3894: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`
3895: @*/
3896: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A, PetscInt col, Vec *v)
3897: {
3898:   PetscFunctionBegin;
3902:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3903:   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);
3904:   PetscUseMethod(A, "MatDenseRestoreColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3905:   PetscFunctionReturn(PETSC_SUCCESS);
3906: }

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

3911:   Collective

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

3920:   Output Parameter:
3921: . v - the matrix

3923:   Level: intermediate

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

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

3930: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreSubMatrix()`
3931: @*/
3932: PetscErrorCode MatDenseGetSubMatrix(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3933: {
3934:   PetscFunctionBegin;
3941:   PetscAssertPointer(v, 6);
3942:   if (rbegin == PETSC_DECIDE) rbegin = 0;
3943:   if (rend == PETSC_DECIDE) rend = A->rmap->N;
3944:   if (cbegin == PETSC_DECIDE) cbegin = 0;
3945:   if (cend == PETSC_DECIDE) cend = A->cmap->N;
3946:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3947:   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);
3948:   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);
3949:   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);
3950:   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);
3951:   PetscUseMethod(A, "MatDenseGetSubMatrix_C", (Mat, PetscInt, PetscInt, PetscInt, PetscInt, Mat *), (A, rbegin, rend, cbegin, cend, v));
3952:   PetscFunctionReturn(PETSC_SUCCESS);
3953: }

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

3958:   Collective

3960:   Input Parameters:
3961: + A - the `Mat` object
3962: - v - the `Mat` object (may be `NULL`)

3964:   Level: intermediate

3966: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseGetSubMatrix()`
3967: @*/
3968: PetscErrorCode MatDenseRestoreSubMatrix(Mat A, Mat *v)
3969: {
3970:   PetscFunctionBegin;
3973:   PetscAssertPointer(v, 2);
3974:   PetscUseMethod(A, "MatDenseRestoreSubMatrix_C", (Mat, Mat *), (A, v));
3975:   PetscFunctionReturn(PETSC_SUCCESS);
3976: }

3978: #include <petscblaslapack.h>
3979: #include <petsc/private/kernels/blockinvert.h>

3981: PetscErrorCode MatSeqDenseInvert(Mat A)
3982: {
3983:   PetscInt        m;
3984:   const PetscReal shift = 0.0;
3985:   PetscBool       allowzeropivot, zeropivotdetected = PETSC_FALSE;
3986:   PetscScalar    *values;

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

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

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