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: 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:     PetscCall(PetscArrayzero(b->v, m * n));
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 * m + 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:     PetscInt j;

254:     for (j = 0; j < X->cmap->n; j++) PetscCallBLAS("BLASaxpy", BLASaxpy_(&m, &alpha, xv + j * ldax, &one, yv + j * lday, &one));
255:   } else {
256:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&N, &alpha, xv, &one, yv, &one));
257:   }
258:   PetscCall(MatDenseRestoreArrayRead(X, &xv));
259:   PetscCall(MatDenseRestoreArray(Y, &yv));
260:   PetscCall(PetscLogFlops(PetscMax(2.0 * N - 1, 0)));
261:   PetscFunctionReturn(PETSC_SUCCESS);
262: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

417: static PetscErrorCode MatConjugate_SeqDense(Mat);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1058: PetscErrorCode MatMultTranspose_SeqDense(Mat A, Vec xx, Vec yy)
1059: {
1060:   Mat_SeqDense      *mat = (Mat_SeqDense *)A->data;
1061:   const PetscScalar *v   = mat->v, *x;
1062:   PetscScalar       *y;
1063:   PetscBLASInt       m, n, _One = 1;
1064:   PetscScalar        _DOne = 1.0, _DZero = 0.0;

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

1083: PetscErrorCode MatMult_SeqDense(Mat A, Vec xx, Vec yy)
1084: {
1085:   Mat_SeqDense      *mat = (Mat_SeqDense *)A->data;
1086:   PetscScalar       *y, _DOne = 1.0, _DZero = 0.0;
1087:   PetscBLASInt       m, n, _One             = 1;
1088:   const PetscScalar *v = mat->v, *x;

1090:   PetscFunctionBegin;
1091:   PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1092:   PetscCall(PetscBLASIntCast(A->cmap->n, &n));
1093:   PetscCall(VecGetArrayRead(xx, &x));
1094:   PetscCall(VecGetArrayWrite(yy, &y));
1095:   if (!A->rmap->n || !A->cmap->n) {
1096:     PetscBLASInt i;
1097:     for (i = 0; i < m; i++) y[i] = 0.0;
1098:   } else {
1099:     PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v, &(mat->lda), x, &_One, &_DZero, y, &_One));
1100:     PetscCall(PetscLogFlops(2.0 * A->rmap->n * A->cmap->n - A->rmap->n));
1101:   }
1102:   PetscCall(VecRestoreArrayRead(xx, &x));
1103:   PetscCall(VecRestoreArrayWrite(yy, &y));
1104:   PetscFunctionReturn(PETSC_SUCCESS);
1105: }

1107: PetscErrorCode MatMultAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1108: {
1109:   Mat_SeqDense      *mat = (Mat_SeqDense *)A->data;
1110:   const PetscScalar *v   = mat->v, *x;
1111:   PetscScalar       *y, _DOne = 1.0;
1112:   PetscBLASInt       m, n, _One = 1;

1114:   PetscFunctionBegin;
1115:   PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1116:   PetscCall(PetscBLASIntCast(A->cmap->n, &n));
1117:   PetscCall(VecCopy(zz, yy));
1118:   if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
1119:   PetscCall(VecGetArrayRead(xx, &x));
1120:   PetscCall(VecGetArray(yy, &y));
1121:   PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v, &(mat->lda), x, &_One, &_DOne, y, &_One));
1122:   PetscCall(VecRestoreArrayRead(xx, &x));
1123:   PetscCall(VecRestoreArray(yy, &y));
1124:   PetscCall(PetscLogFlops(2.0 * A->rmap->n * A->cmap->n));
1125:   PetscFunctionReturn(PETSC_SUCCESS);
1126: }

1128: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1129: {
1130:   Mat_SeqDense      *mat = (Mat_SeqDense *)A->data;
1131:   const PetscScalar *v   = mat->v, *x;
1132:   PetscScalar       *y;
1133:   PetscBLASInt       m, n, _One = 1;
1134:   PetscScalar        _DOne = 1.0;

1136:   PetscFunctionBegin;
1137:   PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1138:   PetscCall(PetscBLASIntCast(A->cmap->n, &n));
1139:   PetscCall(VecCopy(zz, yy));
1140:   if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
1141:   PetscCall(VecGetArrayRead(xx, &x));
1142:   PetscCall(VecGetArray(yy, &y));
1143:   PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v, &(mat->lda), x, &_One, &_DOne, y, &_One));
1144:   PetscCall(VecRestoreArrayRead(xx, &x));
1145:   PetscCall(VecRestoreArray(yy, &y));
1146:   PetscCall(PetscLogFlops(2.0 * A->rmap->n * A->cmap->n));
1147:   PetscFunctionReturn(PETSC_SUCCESS);
1148: }

1150: static PetscErrorCode MatGetRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1151: {
1152:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1153:   PetscInt      i;

1155:   PetscFunctionBegin;
1156:   *ncols = A->cmap->n;
1157:   if (cols) {
1158:     PetscCall(PetscMalloc1(A->cmap->n, cols));
1159:     for (i = 0; i < A->cmap->n; i++) (*cols)[i] = i;
1160:   }
1161:   if (vals) {
1162:     const PetscScalar *v;

1164:     PetscCall(MatDenseGetArrayRead(A, &v));
1165:     PetscCall(PetscMalloc1(A->cmap->n, vals));
1166:     v += row;
1167:     for (i = 0; i < A->cmap->n; i++) {
1168:       (*vals)[i] = *v;
1169:       v += mat->lda;
1170:     }
1171:     PetscCall(MatDenseRestoreArrayRead(A, &v));
1172:   }
1173:   PetscFunctionReturn(PETSC_SUCCESS);
1174: }

1176: static PetscErrorCode MatRestoreRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1177: {
1178:   PetscFunctionBegin;
1179:   if (ncols) *ncols = 0;
1180:   if (cols) PetscCall(PetscFree(*cols));
1181:   if (vals) PetscCall(PetscFree(*vals));
1182:   PetscFunctionReturn(PETSC_SUCCESS);
1183: }

1185: static PetscErrorCode MatSetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], const PetscScalar v[], InsertMode addv)
1186: {
1187:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1188:   PetscScalar  *av;
1189:   PetscInt      i, j, idx = 0;
1190: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1191:   PetscOffloadMask oldf;
1192: #endif

1194:   PetscFunctionBegin;
1195:   PetscCall(MatDenseGetArray(A, &av));
1196:   if (!mat->roworiented) {
1197:     if (addv == INSERT_VALUES) {
1198:       for (j = 0; j < n; j++) {
1199:         if (indexn[j] < 0) {
1200:           idx += m;
1201:           continue;
1202:         }
1203:         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);
1204:         for (i = 0; i < m; i++) {
1205:           if (indexm[i] < 0) {
1206:             idx++;
1207:             continue;
1208:           }
1209:           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);
1210:           av[indexn[j] * mat->lda + indexm[i]] = v[idx++];
1211:         }
1212:       }
1213:     } else {
1214:       for (j = 0; j < n; j++) {
1215:         if (indexn[j] < 0) {
1216:           idx += m;
1217:           continue;
1218:         }
1219:         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);
1220:         for (i = 0; i < m; i++) {
1221:           if (indexm[i] < 0) {
1222:             idx++;
1223:             continue;
1224:           }
1225:           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);
1226:           av[indexn[j] * mat->lda + indexm[i]] += v[idx++];
1227:         }
1228:       }
1229:     }
1230:   } else {
1231:     if (addv == INSERT_VALUES) {
1232:       for (i = 0; i < m; i++) {
1233:         if (indexm[i] < 0) {
1234:           idx += n;
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:         for (j = 0; j < n; j++) {
1239:           if (indexn[j] < 0) {
1240:             idx++;
1241:             continue;
1242:           }
1243:           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);
1244:           av[indexn[j] * mat->lda + indexm[i]] = v[idx++];
1245:         }
1246:       }
1247:     } else {
1248:       for (i = 0; i < m; i++) {
1249:         if (indexm[i] < 0) {
1250:           idx += n;
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:         for (j = 0; j < n; j++) {
1255:           if (indexn[j] < 0) {
1256:             idx++;
1257:             continue;
1258:           }
1259:           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);
1260:           av[indexn[j] * mat->lda + indexm[i]] += v[idx++];
1261:         }
1262:       }
1263:     }
1264:   }
1265:   /* hack to prevent unneeded copy to the GPU while returning the array */
1266: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1267:   oldf           = A->offloadmask;
1268:   A->offloadmask = PETSC_OFFLOAD_GPU;
1269: #endif
1270:   PetscCall(MatDenseRestoreArray(A, &av));
1271: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1272:   A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1273: #endif
1274:   PetscFunctionReturn(PETSC_SUCCESS);
1275: }

1277: static PetscErrorCode MatGetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], PetscScalar v[])
1278: {
1279:   Mat_SeqDense      *mat = (Mat_SeqDense *)A->data;
1280:   const PetscScalar *vv;
1281:   PetscInt           i, j;

1283:   PetscFunctionBegin;
1284:   PetscCall(MatDenseGetArrayRead(A, &vv));
1285:   /* row-oriented output */
1286:   for (i = 0; i < m; i++) {
1287:     if (indexm[i] < 0) {
1288:       v += n;
1289:       continue;
1290:     }
1291:     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);
1292:     for (j = 0; j < n; j++) {
1293:       if (indexn[j] < 0) {
1294:         v++;
1295:         continue;
1296:       }
1297:       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);
1298:       *v++ = vv[indexn[j] * mat->lda + indexm[i]];
1299:     }
1300:   }
1301:   PetscCall(MatDenseRestoreArrayRead(A, &vv));
1302:   PetscFunctionReturn(PETSC_SUCCESS);
1303: }

1305: PetscErrorCode MatView_Dense_Binary(Mat mat, PetscViewer viewer)
1306: {
1307:   PetscBool          skipHeader;
1308:   PetscViewerFormat  format;
1309:   PetscInt           header[4], M, N, m, lda, i, j, k;
1310:   const PetscScalar *v;
1311:   PetscScalar       *vwork;

1313:   PetscFunctionBegin;
1314:   PetscCall(PetscViewerSetUp(viewer));
1315:   PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
1316:   PetscCall(PetscViewerGetFormat(viewer, &format));
1317:   if (skipHeader) format = PETSC_VIEWER_NATIVE;

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

1321:   /* write matrix header */
1322:   header[0] = MAT_FILE_CLASSID;
1323:   header[1] = M;
1324:   header[2] = N;
1325:   header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M * N;
1326:   if (!skipHeader) PetscCall(PetscViewerBinaryWrite(viewer, header, 4, PETSC_INT));

1328:   PetscCall(MatGetLocalSize(mat, &m, NULL));
1329:   if (format != PETSC_VIEWER_NATIVE) {
1330:     PetscInt nnz = m * N, *iwork;
1331:     /* store row lengths for each row */
1332:     PetscCall(PetscMalloc1(nnz, &iwork));
1333:     for (i = 0; i < m; i++) iwork[i] = N;
1334:     PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1335:     /* store column indices (zero start index) */
1336:     for (k = 0, i = 0; i < m; i++)
1337:       for (j = 0; j < N; j++, k++) iwork[k] = j;
1338:     PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1339:     PetscCall(PetscFree(iwork));
1340:   }
1341:   /* store matrix values as a dense matrix in row major order */
1342:   PetscCall(PetscMalloc1(m * N, &vwork));
1343:   PetscCall(MatDenseGetArrayRead(mat, &v));
1344:   PetscCall(MatDenseGetLDA(mat, &lda));
1345:   for (k = 0, i = 0; i < m; i++)
1346:     for (j = 0; j < N; j++, k++) vwork[k] = v[i + lda * j];
1347:   PetscCall(MatDenseRestoreArrayRead(mat, &v));
1348:   PetscCall(PetscViewerBinaryWriteAll(viewer, vwork, m * N, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1349:   PetscCall(PetscFree(vwork));
1350:   PetscFunctionReturn(PETSC_SUCCESS);
1351: }

1353: PetscErrorCode MatLoad_Dense_Binary(Mat mat, PetscViewer viewer)
1354: {
1355:   PetscBool    skipHeader;
1356:   PetscInt     header[4], M, N, m, nz, lda, i, j, k;
1357:   PetscInt     rows, cols;
1358:   PetscScalar *v, *vwork;

1360:   PetscFunctionBegin;
1361:   PetscCall(PetscViewerSetUp(viewer));
1362:   PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));

1364:   if (!skipHeader) {
1365:     PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT));
1366:     PetscCheck(header[0] == MAT_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file");
1367:     M = header[1];
1368:     N = header[2];
1369:     PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M);
1370:     PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N);
1371:     nz = header[3];
1372:     PetscCheck(nz == MATRIX_BINARY_FORMAT_DENSE || nz >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Unknown matrix format %" PetscInt_FMT " in file", nz);
1373:   } else {
1374:     PetscCall(MatGetSize(mat, &M, &N));
1375:     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");
1376:     nz = MATRIX_BINARY_FORMAT_DENSE;
1377:   }

1379:   /* setup global sizes if not set */
1380:   if (mat->rmap->N < 0) mat->rmap->N = M;
1381:   if (mat->cmap->N < 0) mat->cmap->N = N;
1382:   PetscCall(MatSetUp(mat));
1383:   /* check if global sizes are correct */
1384:   PetscCall(MatGetSize(mat, &rows, &cols));
1385:   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);

1387:   PetscCall(MatGetSize(mat, NULL, &N));
1388:   PetscCall(MatGetLocalSize(mat, &m, NULL));
1389:   PetscCall(MatDenseGetArray(mat, &v));
1390:   PetscCall(MatDenseGetLDA(mat, &lda));
1391:   if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */
1392:     PetscInt nnz = m * N;
1393:     /* read in matrix values */
1394:     PetscCall(PetscMalloc1(nnz, &vwork));
1395:     PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1396:     /* store values in column major order */
1397:     for (j = 0; j < N; j++)
1398:       for (i = 0; i < m; i++) v[i + lda * j] = vwork[i * N + j];
1399:     PetscCall(PetscFree(vwork));
1400:   } else { /* matrix in file is sparse format */
1401:     PetscInt nnz = 0, *rlens, *icols;
1402:     /* read in row lengths */
1403:     PetscCall(PetscMalloc1(m, &rlens));
1404:     PetscCall(PetscViewerBinaryReadAll(viewer, rlens, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1405:     for (i = 0; i < m; i++) nnz += rlens[i];
1406:     /* read in column indices and values */
1407:     PetscCall(PetscMalloc2(nnz, &icols, nnz, &vwork));
1408:     PetscCall(PetscViewerBinaryReadAll(viewer, icols, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1409:     PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1410:     /* store values in column major order */
1411:     for (k = 0, i = 0; i < m; i++)
1412:       for (j = 0; j < rlens[i]; j++, k++) v[i + lda * icols[k]] = vwork[k];
1413:     PetscCall(PetscFree(rlens));
1414:     PetscCall(PetscFree2(icols, vwork));
1415:   }
1416:   PetscCall(MatDenseRestoreArray(mat, &v));
1417:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
1418:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
1419:   PetscFunctionReturn(PETSC_SUCCESS);
1420: }

1422: PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1423: {
1424:   PetscBool isbinary, ishdf5;

1426:   PetscFunctionBegin;
1429:   /* force binary viewer to load .info file if it has not yet done so */
1430:   PetscCall(PetscViewerSetUp(viewer));
1431:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1432:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1433:   if (isbinary) {
1434:     PetscCall(MatLoad_Dense_Binary(newMat, viewer));
1435:   } else if (ishdf5) {
1436: #if defined(PETSC_HAVE_HDF5)
1437:     PetscCall(MatLoad_Dense_HDF5(newMat, viewer));
1438: #else
1439:     SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1440: #endif
1441:   } else {
1442:     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);
1443:   }
1444:   PetscFunctionReturn(PETSC_SUCCESS);
1445: }

1447: static PetscErrorCode MatView_SeqDense_ASCII(Mat A, PetscViewer viewer)
1448: {
1449:   Mat_SeqDense     *a = (Mat_SeqDense *)A->data;
1450:   PetscInt          i, j;
1451:   const char       *name;
1452:   PetscScalar      *v, *av;
1453:   PetscViewerFormat format;
1454: #if defined(PETSC_USE_COMPLEX)
1455:   PetscBool allreal = PETSC_TRUE;
1456: #endif

1458:   PetscFunctionBegin;
1459:   PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&av));
1460:   PetscCall(PetscViewerGetFormat(viewer, &format));
1461:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1462:     PetscFunctionReturn(PETSC_SUCCESS); /* do nothing for now */
1463:   } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1464:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1465:     for (i = 0; i < A->rmap->n; i++) {
1466:       v = av + i;
1467:       PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i));
1468:       for (j = 0; j < A->cmap->n; j++) {
1469: #if defined(PETSC_USE_COMPLEX)
1470:         if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1471:           PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", j, (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1472:         } else if (PetscRealPart(*v)) {
1473:           PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)PetscRealPart(*v)));
1474:         }
1475: #else
1476:         if (*v) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)*v));
1477: #endif
1478:         v += a->lda;
1479:       }
1480:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1481:     }
1482:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1483:   } else {
1484:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1485: #if defined(PETSC_USE_COMPLEX)
1486:     /* determine if matrix has all real values */
1487:     for (j = 0; j < A->cmap->n; j++) {
1488:       v = av + j * a->lda;
1489:       for (i = 0; i < A->rmap->n; i++) {
1490:         if (PetscImaginaryPart(v[i])) {
1491:           allreal = PETSC_FALSE;
1492:           break;
1493:         }
1494:       }
1495:     }
1496: #endif
1497:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
1498:       PetscCall(PetscObjectGetName((PetscObject)A, &name));
1499:       PetscCall(PetscViewerASCIIPrintf(viewer, "%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n", A->rmap->n, A->cmap->n));
1500:       PetscCall(PetscViewerASCIIPrintf(viewer, "%s = zeros(%" PetscInt_FMT ",%" PetscInt_FMT ");\n", name, A->rmap->n, A->cmap->n));
1501:       PetscCall(PetscViewerASCIIPrintf(viewer, "%s = [\n", name));
1502:     }

1504:     for (i = 0; i < A->rmap->n; i++) {
1505:       v = av + i;
1506:       for (j = 0; j < A->cmap->n; j++) {
1507: #if defined(PETSC_USE_COMPLEX)
1508:         if (allreal) {
1509:           PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)PetscRealPart(*v)));
1510:         } else {
1511:           PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei ", (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1512:         }
1513: #else
1514:         PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)*v));
1515: #endif
1516:         v += a->lda;
1517:       }
1518:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1519:     }
1520:     if (format == PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer, "];\n"));
1521:     PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1522:   }
1523:   PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&av));
1524:   PetscCall(PetscViewerFlush(viewer));
1525:   PetscFunctionReturn(PETSC_SUCCESS);
1526: }

1528: #include <petscdraw.h>
1529: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw, void *Aa)
1530: {
1531:   Mat                A = (Mat)Aa;
1532:   PetscInt           m = A->rmap->n, n = A->cmap->n, i, j;
1533:   int                color = PETSC_DRAW_WHITE;
1534:   const PetscScalar *v;
1535:   PetscViewer        viewer;
1536:   PetscReal          xl, yl, xr, yr, x_l, x_r, y_l, y_r;
1537:   PetscViewerFormat  format;

1539:   PetscFunctionBegin;
1540:   PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer));
1541:   PetscCall(PetscViewerGetFormat(viewer, &format));
1542:   PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));

1544:   /* Loop over matrix elements drawing boxes */
1545:   PetscCall(MatDenseGetArrayRead(A, &v));
1546:   if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1547:     PetscDrawCollectiveBegin(draw);
1548:     /* Blue for negative and Red for positive */
1549:     for (j = 0; j < n; j++) {
1550:       x_l = j;
1551:       x_r = x_l + 1.0;
1552:       for (i = 0; i < m; i++) {
1553:         y_l = m - i - 1.0;
1554:         y_r = y_l + 1.0;
1555:         if (PetscRealPart(v[j * m + i]) > 0.) color = PETSC_DRAW_RED;
1556:         else if (PetscRealPart(v[j * m + i]) < 0.) color = PETSC_DRAW_BLUE;
1557:         else continue;
1558:         PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1559:       }
1560:     }
1561:     PetscDrawCollectiveEnd(draw);
1562:   } else {
1563:     /* use contour shading to indicate magnitude of values */
1564:     /* first determine max of all nonzero values */
1565:     PetscReal minv = 0.0, maxv = 0.0;
1566:     PetscDraw popup;

1568:     for (i = 0; i < m * n; i++) {
1569:       if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1570:     }
1571:     if (minv >= maxv) maxv = minv + PETSC_SMALL;
1572:     PetscCall(PetscDrawGetPopup(draw, &popup));
1573:     PetscCall(PetscDrawScalePopup(popup, minv, maxv));

1575:     PetscDrawCollectiveBegin(draw);
1576:     for (j = 0; j < n; j++) {
1577:       x_l = j;
1578:       x_r = x_l + 1.0;
1579:       for (i = 0; i < m; i++) {
1580:         y_l   = m - i - 1.0;
1581:         y_r   = y_l + 1.0;
1582:         color = PetscDrawRealToColor(PetscAbsScalar(v[j * m + i]), minv, maxv);
1583:         PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1584:       }
1585:     }
1586:     PetscDrawCollectiveEnd(draw);
1587:   }
1588:   PetscCall(MatDenseRestoreArrayRead(A, &v));
1589:   PetscFunctionReturn(PETSC_SUCCESS);
1590: }

1592: static PetscErrorCode MatView_SeqDense_Draw(Mat A, PetscViewer viewer)
1593: {
1594:   PetscDraw draw;
1595:   PetscBool isnull;
1596:   PetscReal xr, yr, xl, yl, h, w;

1598:   PetscFunctionBegin;
1599:   PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1600:   PetscCall(PetscDrawIsNull(draw, &isnull));
1601:   if (isnull) PetscFunctionReturn(PETSC_SUCCESS);

1603:   xr = A->cmap->n;
1604:   yr = A->rmap->n;
1605:   h  = yr / 10.0;
1606:   w  = xr / 10.0;
1607:   xr += w;
1608:   yr += h;
1609:   xl = -w;
1610:   yl = -h;
1611:   PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr));
1612:   PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer));
1613:   PetscCall(PetscDrawZoom(draw, MatView_SeqDense_Draw_Zoom, A));
1614:   PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL));
1615:   PetscCall(PetscDrawSave(draw));
1616:   PetscFunctionReturn(PETSC_SUCCESS);
1617: }

1619: PetscErrorCode MatView_SeqDense(Mat A, PetscViewer viewer)
1620: {
1621:   PetscBool iascii, isbinary, isdraw;

1623:   PetscFunctionBegin;
1624:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1625:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1626:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1627:   if (iascii) PetscCall(MatView_SeqDense_ASCII(A, viewer));
1628:   else if (isbinary) PetscCall(MatView_Dense_Binary(A, viewer));
1629:   else if (isdraw) PetscCall(MatView_SeqDense_Draw(A, viewer));
1630:   PetscFunctionReturn(PETSC_SUCCESS);
1631: }

1633: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A, const PetscScalar *array)
1634: {
1635:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

1637:   PetscFunctionBegin;
1638:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1639:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1640:   PetscCheck(!a->unplacedarray, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreArray() first");
1641:   a->unplacedarray       = a->v;
1642:   a->unplaced_user_alloc = a->user_alloc;
1643:   a->v                   = (PetscScalar *)array;
1644:   a->user_alloc          = PETSC_TRUE;
1645: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1646:   A->offloadmask = PETSC_OFFLOAD_CPU;
1647: #endif
1648:   PetscFunctionReturn(PETSC_SUCCESS);
1649: }

1651: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1652: {
1653:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

1655:   PetscFunctionBegin;
1656:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1657:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1658:   a->v             = a->unplacedarray;
1659:   a->user_alloc    = a->unplaced_user_alloc;
1660:   a->unplacedarray = NULL;
1661: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1662:   A->offloadmask = PETSC_OFFLOAD_CPU;
1663: #endif
1664:   PetscFunctionReturn(PETSC_SUCCESS);
1665: }

1667: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A, const PetscScalar *array)
1668: {
1669:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

1671:   PetscFunctionBegin;
1672:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1673:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1674:   if (!a->user_alloc) PetscCall(PetscFree(a->v));
1675:   a->v          = (PetscScalar *)array;
1676:   a->user_alloc = PETSC_FALSE;
1677: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1678:   A->offloadmask = PETSC_OFFLOAD_CPU;
1679: #endif
1680:   PetscFunctionReturn(PETSC_SUCCESS);
1681: }

1683: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1684: {
1685:   Mat_SeqDense *l = (Mat_SeqDense *)mat->data;

1687:   PetscFunctionBegin;
1688: #if defined(PETSC_USE_LOG)
1689:   PetscCall(PetscLogObjectState((PetscObject)mat, "Rows %" PetscInt_FMT " Cols %" PetscInt_FMT, mat->rmap->n, mat->cmap->n));
1690: #endif
1691:   PetscCall(VecDestroy(&(l->qrrhs)));
1692:   PetscCall(PetscFree(l->tau));
1693:   PetscCall(PetscFree(l->pivots));
1694:   PetscCall(PetscFree(l->fwork));
1695:   if (!l->user_alloc) PetscCall(PetscFree(l->v));
1696:   if (!l->unplaced_user_alloc) PetscCall(PetscFree(l->unplacedarray));
1697:   PetscCheck(!l->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1698:   PetscCheck(!l->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1699:   PetscCall(VecDestroy(&l->cvec));
1700:   PetscCall(MatDestroy(&l->cmat));
1701:   PetscCall(PetscFree(mat->data));

1703:   PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL));
1704:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactor_C", NULL));
1705:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorSymbolic_C", NULL));
1706:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorNumeric_C", NULL));
1707:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetLDA_C", NULL));
1708:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseSetLDA_C", NULL));
1709:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArray_C", NULL));
1710:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArray_C", NULL));
1711:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDensePlaceArray_C", NULL));
1712:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseResetArray_C", NULL));
1713:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseReplaceArray_C", NULL));
1714:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayRead_C", NULL));
1715:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayRead_C", NULL));
1716:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayWrite_C", NULL));
1717:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayWrite_C", NULL));
1718:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqaij_C", NULL));
1719: #if defined(PETSC_HAVE_ELEMENTAL)
1720:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_elemental_C", NULL));
1721: #endif
1722: #if defined(PETSC_HAVE_SCALAPACK)
1723:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_scalapack_C", NULL));
1724: #endif
1725: #if defined(PETSC_HAVE_CUDA)
1726:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensecuda_C", NULL));
1727:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", NULL));
1728:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdense_C", NULL));
1729:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensecuda_C", NULL));
1730: #endif
1731: #if defined(PETSC_HAVE_HIP)
1732:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensehip_C", NULL));
1733:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", NULL));
1734:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdense_C", NULL));
1735:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensehip_C", NULL));
1736: #endif
1737:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatSeqDenseSetPreallocation_C", NULL));
1738:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqaij_seqdense_C", NULL));
1739:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdense_C", NULL));
1740:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqbaij_seqdense_C", NULL));
1741:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqsbaij_seqdense_C", NULL));

1743:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumn_C", NULL));
1744:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumn_C", NULL));
1745:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVec_C", NULL));
1746:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVec_C", NULL));
1747:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecRead_C", NULL));
1748:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecRead_C", NULL));
1749:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecWrite_C", NULL));
1750:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecWrite_C", NULL));
1751:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetSubMatrix_C", NULL));
1752:   PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreSubMatrix_C", NULL));
1753:   PetscFunctionReturn(PETSC_SUCCESS);
1754: }

1756: static PetscErrorCode MatTranspose_SeqDense(Mat A, MatReuse reuse, Mat *matout)
1757: {
1758:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1759:   PetscInt      k, j, m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1760:   PetscScalar  *v, tmp;

1762:   PetscFunctionBegin;
1763:   if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *matout));
1764:   if (reuse == MAT_INPLACE_MATRIX) {
1765:     if (m == n) { /* in place transpose */
1766:       PetscCall(MatDenseGetArray(A, &v));
1767:       for (j = 0; j < m; j++) {
1768:         for (k = 0; k < j; k++) {
1769:           tmp          = v[j + k * M];
1770:           v[j + k * M] = v[k + j * M];
1771:           v[k + j * M] = tmp;
1772:         }
1773:       }
1774:       PetscCall(MatDenseRestoreArray(A, &v));
1775:     } else { /* reuse memory, temporary allocates new memory */
1776:       PetscScalar *v2;
1777:       PetscLayout  tmplayout;

1779:       PetscCall(PetscMalloc1((size_t)m * n, &v2));
1780:       PetscCall(MatDenseGetArray(A, &v));
1781:       for (j = 0; j < n; j++) {
1782:         for (k = 0; k < m; k++) v2[j + (size_t)k * n] = v[k + (size_t)j * M];
1783:       }
1784:       PetscCall(PetscArraycpy(v, v2, (size_t)m * n));
1785:       PetscCall(PetscFree(v2));
1786:       PetscCall(MatDenseRestoreArray(A, &v));
1787:       /* cleanup size dependent quantities */
1788:       PetscCall(VecDestroy(&mat->cvec));
1789:       PetscCall(MatDestroy(&mat->cmat));
1790:       PetscCall(PetscFree(mat->pivots));
1791:       PetscCall(PetscFree(mat->fwork));
1792:       /* swap row/col layouts */
1793:       mat->lda  = n;
1794:       tmplayout = A->rmap;
1795:       A->rmap   = A->cmap;
1796:       A->cmap   = tmplayout;
1797:     }
1798:   } else { /* out-of-place transpose */
1799:     Mat           tmat;
1800:     Mat_SeqDense *tmatd;
1801:     PetscScalar  *v2;
1802:     PetscInt      M2;

1804:     if (reuse == MAT_INITIAL_MATRIX) {
1805:       PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &tmat));
1806:       PetscCall(MatSetSizes(tmat, A->cmap->n, A->rmap->n, A->cmap->n, A->rmap->n));
1807:       PetscCall(MatSetType(tmat, ((PetscObject)A)->type_name));
1808:       PetscCall(MatSeqDenseSetPreallocation(tmat, NULL));
1809:     } else tmat = *matout;

1811:     PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&v));
1812:     PetscCall(MatDenseGetArray(tmat, &v2));
1813:     tmatd = (Mat_SeqDense *)tmat->data;
1814:     M2    = tmatd->lda;
1815:     for (j = 0; j < n; j++) {
1816:       for (k = 0; k < m; k++) v2[j + k * M2] = v[k + j * M];
1817:     }
1818:     PetscCall(MatDenseRestoreArray(tmat, &v2));
1819:     PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&v));
1820:     PetscCall(MatAssemblyBegin(tmat, MAT_FINAL_ASSEMBLY));
1821:     PetscCall(MatAssemblyEnd(tmat, MAT_FINAL_ASSEMBLY));
1822:     *matout = tmat;
1823:   }
1824:   PetscFunctionReturn(PETSC_SUCCESS);
1825: }

1827: static PetscErrorCode MatEqual_SeqDense(Mat A1, Mat A2, PetscBool *flg)
1828: {
1829:   Mat_SeqDense      *mat1 = (Mat_SeqDense *)A1->data;
1830:   Mat_SeqDense      *mat2 = (Mat_SeqDense *)A2->data;
1831:   PetscInt           i;
1832:   const PetscScalar *v1, *v2;

1834:   PetscFunctionBegin;
1835:   if (A1->rmap->n != A2->rmap->n) {
1836:     *flg = PETSC_FALSE;
1837:     PetscFunctionReturn(PETSC_SUCCESS);
1838:   }
1839:   if (A1->cmap->n != A2->cmap->n) {
1840:     *flg = PETSC_FALSE;
1841:     PetscFunctionReturn(PETSC_SUCCESS);
1842:   }
1843:   PetscCall(MatDenseGetArrayRead(A1, &v1));
1844:   PetscCall(MatDenseGetArrayRead(A2, &v2));
1845:   for (i = 0; i < A1->cmap->n; i++) {
1846:     PetscCall(PetscArraycmp(v1, v2, A1->rmap->n, flg));
1847:     if (*flg == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS);
1848:     v1 += mat1->lda;
1849:     v2 += mat2->lda;
1850:   }
1851:   PetscCall(MatDenseRestoreArrayRead(A1, &v1));
1852:   PetscCall(MatDenseRestoreArrayRead(A2, &v2));
1853:   *flg = PETSC_TRUE;
1854:   PetscFunctionReturn(PETSC_SUCCESS);
1855: }

1857: static PetscErrorCode MatGetDiagonal_SeqDense(Mat A, Vec v)
1858: {
1859:   Mat_SeqDense      *mat = (Mat_SeqDense *)A->data;
1860:   PetscInt           i, n, len;
1861:   PetscScalar       *x;
1862:   const PetscScalar *vv;

1864:   PetscFunctionBegin;
1865:   PetscCall(VecGetSize(v, &n));
1866:   PetscCall(VecGetArray(v, &x));
1867:   len = PetscMin(A->rmap->n, A->cmap->n);
1868:   PetscCall(MatDenseGetArrayRead(A, &vv));
1869:   PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming mat and vec");
1870:   for (i = 0; i < len; i++) x[i] = vv[i * mat->lda + i];
1871:   PetscCall(MatDenseRestoreArrayRead(A, &vv));
1872:   PetscCall(VecRestoreArray(v, &x));
1873:   PetscFunctionReturn(PETSC_SUCCESS);
1874: }

1876: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A, Vec ll, Vec rr)
1877: {
1878:   Mat_SeqDense      *mat = (Mat_SeqDense *)A->data;
1879:   const PetscScalar *l, *r;
1880:   PetscScalar        x, *v, *vv;
1881:   PetscInt           i, j, m = A->rmap->n, n = A->cmap->n;

1883:   PetscFunctionBegin;
1884:   PetscCall(MatDenseGetArray(A, &vv));
1885:   if (ll) {
1886:     PetscCall(VecGetSize(ll, &m));
1887:     PetscCall(VecGetArrayRead(ll, &l));
1888:     PetscCheck(m == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vec wrong size");
1889:     for (i = 0; i < m; i++) {
1890:       x = l[i];
1891:       v = vv + i;
1892:       for (j = 0; j < n; j++) {
1893:         (*v) *= x;
1894:         v += mat->lda;
1895:       }
1896:     }
1897:     PetscCall(VecRestoreArrayRead(ll, &l));
1898:     PetscCall(PetscLogFlops(1.0 * n * m));
1899:   }
1900:   if (rr) {
1901:     PetscCall(VecGetSize(rr, &n));
1902:     PetscCall(VecGetArrayRead(rr, &r));
1903:     PetscCheck(n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vec wrong size");
1904:     for (i = 0; i < n; i++) {
1905:       x = r[i];
1906:       v = vv + i * mat->lda;
1907:       for (j = 0; j < m; j++) (*v++) *= x;
1908:     }
1909:     PetscCall(VecRestoreArrayRead(rr, &r));
1910:     PetscCall(PetscLogFlops(1.0 * n * m));
1911:   }
1912:   PetscCall(MatDenseRestoreArray(A, &vv));
1913:   PetscFunctionReturn(PETSC_SUCCESS);
1914: }

1916: PetscErrorCode MatNorm_SeqDense(Mat A, NormType type, PetscReal *nrm)
1917: {
1918:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1919:   PetscScalar  *v, *vv;
1920:   PetscReal     sum = 0.0;
1921:   PetscInt      lda, m = A->rmap->n, i, j;

1923:   PetscFunctionBegin;
1924:   PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&vv));
1925:   PetscCall(MatDenseGetLDA(A, &lda));
1926:   v = vv;
1927:   if (type == NORM_FROBENIUS) {
1928:     if (lda > m) {
1929:       for (j = 0; j < A->cmap->n; j++) {
1930:         v = vv + j * lda;
1931:         for (i = 0; i < m; i++) {
1932:           sum += PetscRealPart(PetscConj(*v) * (*v));
1933:           v++;
1934:         }
1935:       }
1936:     } else {
1937: #if defined(PETSC_USE_REAL___FP16)
1938:       PetscBLASInt one = 1, cnt = A->cmap->n * A->rmap->n;
1939:       PetscCallBLAS("BLASnrm2", *nrm = BLASnrm2_(&cnt, v, &one));
1940:     }
1941: #else
1942:       for (i = 0; i < A->cmap->n * A->rmap->n; i++) {
1943:         sum += PetscRealPart(PetscConj(*v) * (*v));
1944:         v++;
1945:       }
1946:     }
1947:     *nrm = PetscSqrtReal(sum);
1948: #endif
1949:     PetscCall(PetscLogFlops(2.0 * A->cmap->n * A->rmap->n));
1950:   } else if (type == NORM_1) {
1951:     *nrm = 0.0;
1952:     for (j = 0; j < A->cmap->n; j++) {
1953:       v   = vv + j * mat->lda;
1954:       sum = 0.0;
1955:       for (i = 0; i < A->rmap->n; i++) {
1956:         sum += PetscAbsScalar(*v);
1957:         v++;
1958:       }
1959:       if (sum > *nrm) *nrm = sum;
1960:     }
1961:     PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
1962:   } else if (type == NORM_INFINITY) {
1963:     *nrm = 0.0;
1964:     for (j = 0; j < A->rmap->n; j++) {
1965:       v   = vv + j;
1966:       sum = 0.0;
1967:       for (i = 0; i < A->cmap->n; i++) {
1968:         sum += PetscAbsScalar(*v);
1969:         v += mat->lda;
1970:       }
1971:       if (sum > *nrm) *nrm = sum;
1972:     }
1973:     PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
1974:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No two norm");
1975:   PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&vv));
1976:   PetscFunctionReturn(PETSC_SUCCESS);
1977: }

1979: static PetscErrorCode MatSetOption_SeqDense(Mat A, MatOption op, PetscBool flg)
1980: {
1981:   Mat_SeqDense *aij = (Mat_SeqDense *)A->data;

1983:   PetscFunctionBegin;
1984:   switch (op) {
1985:   case MAT_ROW_ORIENTED:
1986:     aij->roworiented = flg;
1987:     break;
1988:   case MAT_NEW_NONZERO_LOCATIONS:
1989:   case MAT_NEW_NONZERO_LOCATION_ERR:
1990:   case MAT_NEW_NONZERO_ALLOCATION_ERR:
1991:   case MAT_FORCE_DIAGONAL_ENTRIES:
1992:   case MAT_KEEP_NONZERO_PATTERN:
1993:   case MAT_IGNORE_OFF_PROC_ENTRIES:
1994:   case MAT_USE_HASH_TABLE:
1995:   case MAT_IGNORE_ZERO_ENTRIES:
1996:   case MAT_IGNORE_LOWER_TRIANGULAR:
1997:   case MAT_SORTED_FULL:
1998:     PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op]));
1999:     break;
2000:   case MAT_SPD:
2001:   case MAT_SYMMETRIC:
2002:   case MAT_STRUCTURALLY_SYMMETRIC:
2003:   case MAT_HERMITIAN:
2004:   case MAT_SYMMETRY_ETERNAL:
2005:   case MAT_STRUCTURAL_SYMMETRY_ETERNAL:
2006:   case MAT_SPD_ETERNAL:
2007:     break;
2008:   default:
2009:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unknown option %s", MatOptions[op]);
2010:   }
2011:   PetscFunctionReturn(PETSC_SUCCESS);
2012: }

2014: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
2015: {
2016:   Mat_SeqDense *l   = (Mat_SeqDense *)A->data;
2017:   PetscInt      lda = l->lda, m = A->rmap->n, n = A->cmap->n, j;
2018:   PetscScalar  *v;

2020:   PetscFunctionBegin;
2021:   PetscCall(MatDenseGetArrayWrite(A, &v));
2022:   if (lda > m) {
2023:     for (j = 0; j < n; j++) PetscCall(PetscArrayzero(v + j * lda, m));
2024:   } else {
2025:     PetscCall(PetscArrayzero(v, PetscInt64Mult(m, n)));
2026:   }
2027:   PetscCall(MatDenseRestoreArrayWrite(A, &v));
2028:   PetscFunctionReturn(PETSC_SUCCESS);
2029: }

2031: static PetscErrorCode MatZeroRows_SeqDense(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2032: {
2033:   Mat_SeqDense      *l = (Mat_SeqDense *)A->data;
2034:   PetscInt           m = l->lda, n = A->cmap->n, i, j;
2035:   PetscScalar       *slot, *bb, *v;
2036:   const PetscScalar *xx;

2038:   PetscFunctionBegin;
2039:   if (PetscDefined(USE_DEBUG)) {
2040:     for (i = 0; i < N; i++) {
2041:       PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
2042:       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);
2043:     }
2044:   }
2045:   if (!N) PetscFunctionReturn(PETSC_SUCCESS);

2047:   /* fix right hand side if needed */
2048:   if (x && b) {
2049:     PetscCall(VecGetArrayRead(x, &xx));
2050:     PetscCall(VecGetArray(b, &bb));
2051:     for (i = 0; i < N; i++) bb[rows[i]] = diag * xx[rows[i]];
2052:     PetscCall(VecRestoreArrayRead(x, &xx));
2053:     PetscCall(VecRestoreArray(b, &bb));
2054:   }

2056:   PetscCall(MatDenseGetArray(A, &v));
2057:   for (i = 0; i < N; i++) {
2058:     slot = v + rows[i];
2059:     for (j = 0; j < n; j++) {
2060:       *slot = 0.0;
2061:       slot += m;
2062:     }
2063:   }
2064:   if (diag != 0.0) {
2065:     PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
2066:     for (i = 0; i < N; i++) {
2067:       slot  = v + (m + 1) * rows[i];
2068:       *slot = diag;
2069:     }
2070:   }
2071:   PetscCall(MatDenseRestoreArray(A, &v));
2072:   PetscFunctionReturn(PETSC_SUCCESS);
2073: }

2075: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A, PetscInt *lda)
2076: {
2077:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;

2079:   PetscFunctionBegin;
2080:   *lda = mat->lda;
2081:   PetscFunctionReturn(PETSC_SUCCESS);
2082: }

2084: PetscErrorCode MatDenseGetArray_SeqDense(Mat A, PetscScalar **array)
2085: {
2086:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;

2088:   PetscFunctionBegin;
2089:   PetscCheck(!mat->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
2090:   *array = mat->v;
2091:   PetscFunctionReturn(PETSC_SUCCESS);
2092: }

2094: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A, PetscScalar **array)
2095: {
2096:   PetscFunctionBegin;
2097:   if (array) *array = NULL;
2098:   PetscFunctionReturn(PETSC_SUCCESS);
2099: }

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

2104:    Not Collective

2106:    Input Parameter:
2107: .  mat - a `MATDENSE` or `MATDENSECUDA` matrix

2109:    Output Parameter:
2110: .   lda - the leading dimension

2112:    Level: intermediate

2114: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseSetLDA()`
2115: @*/
2116: PetscErrorCode MatDenseGetLDA(Mat A, PetscInt *lda)
2117: {
2118:   PetscFunctionBegin;
2121:   MatCheckPreallocated(A, 1);
2122:   PetscUseMethod(A, "MatDenseGetLDA_C", (Mat, PetscInt *), (A, lda));
2123:   PetscFunctionReturn(PETSC_SUCCESS);
2124: }

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

2129:    Not Collective

2131:    Input Parameters:
2132: +  mat - a `MATDENSE` or `MATDENSECUDA` matrix
2133: -  lda - the leading dimension

2135:    Level: intermediate

2137: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetLDA()`
2138: @*/
2139: PetscErrorCode MatDenseSetLDA(Mat A, PetscInt lda)
2140: {
2141:   PetscFunctionBegin;
2143:   PetscTryMethod(A, "MatDenseSetLDA_C", (Mat, PetscInt), (A, lda));
2144:   PetscFunctionReturn(PETSC_SUCCESS);
2145: }

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

2150:    Logically Collective

2152:    Input Parameter:
2153: .  mat - a dense matrix

2155:    Output Parameter:
2156: .   array - pointer to the data

2158:    Level: intermediate

2160:    Fortran Note:
2161:    `MatDenseGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `MatDenseGetArrayF90()`

2163: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2164: @*/
2165: PetscErrorCode MatDenseGetArray(Mat A, PetscScalar **array)
2166: {
2167:   PetscFunctionBegin;
2170:   PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2171:   PetscFunctionReturn(PETSC_SUCCESS);
2172: }

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

2177:    Logically Collective

2179:    Input Parameters:
2180: +  mat - a dense matrix
2181: -  array - pointer to the data (may be `NULL`)

2183:    Level: intermediate

2185:    Fortran Note:
2186:    `MatDenseRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `MatDenseRestoreArrayF90()`

2188: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2189: @*/
2190: PetscErrorCode MatDenseRestoreArray(Mat A, PetscScalar **array)
2191: {
2192:   PetscFunctionBegin;
2195:   PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2196:   PetscCall(PetscObjectStateIncrease((PetscObject)A));
2197: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2198:   A->offloadmask = PETSC_OFFLOAD_CPU;
2199: #endif
2200:   PetscFunctionReturn(PETSC_SUCCESS);
2201: }

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

2206:    Not Collective; No Fortran Support

2208:    Input Parameter:
2209: .  mat - a dense matrix

2211:    Output Parameter:
2212: .   array - pointer to the data

2214:    Level: intermediate

2216: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2217: @*/
2218: PetscErrorCode MatDenseGetArrayRead(Mat A, const PetscScalar **array)
2219: {
2220:   PetscFunctionBegin;
2223:   PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, const PetscScalar **), (A, array));
2224:   PetscFunctionReturn(PETSC_SUCCESS);
2225: }

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

2230:    Not Collective; No Fortran Support

2232:    Input Parameters:
2233: +  mat - a dense matrix
2234: -  array - pointer to the data (may be `NULL`)

2236:    Level: intermediate

2238: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2239: @*/
2240: PetscErrorCode MatDenseRestoreArrayRead(Mat A, const PetscScalar **array)
2241: {
2242:   PetscFunctionBegin;
2245:   PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, const PetscScalar **), (A, array));
2246:   PetscFunctionReturn(PETSC_SUCCESS);
2247: }

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

2252:    Not Collective; No Fortran Support

2254:    Input Parameter:
2255: .  mat - a dense matrix

2257:    Output Parameter:
2258: .   array - pointer to the data

2260:    Level: intermediate

2262: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2263: @*/
2264: PetscErrorCode MatDenseGetArrayWrite(Mat A, PetscScalar **array)
2265: {
2266:   PetscFunctionBegin;
2269:   PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2270:   PetscFunctionReturn(PETSC_SUCCESS);
2271: }

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

2276:    Not Collective; No Fortran Support

2278:    Input Parameters:
2279: +  mat - a dense matrix
2280: -  array - pointer to the data (may be `NULL`)

2282:    Level: intermediate

2284: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2285: @*/
2286: PetscErrorCode MatDenseRestoreArrayWrite(Mat A, PetscScalar **array)
2287: {
2288:   PetscFunctionBegin;
2291:   PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2292:   PetscCall(PetscObjectStateIncrease((PetscObject)A));
2293: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2294:   A->offloadmask = PETSC_OFFLOAD_CPU;
2295: #endif
2296:   PetscFunctionReturn(PETSC_SUCCESS);
2297: }

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

2302:    Logically Collective

2304:    Input Parameter:
2305: .  mat - a dense matrix

2307:    Output Parameters:
2308: +  array - pointer to the data
2309: -  mtype - memory type of the returned pointer

2311:    Level: intermediate

2313:    Notes:
2314:    If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2315:    an array on device is always returned and is guaranteed to contain the matrix's latest data.

2317: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArrayRead()`,
2318:    `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2319: @*/
2320: PetscErrorCode MatDenseGetArrayAndMemType(Mat A, PetscScalar **array, PetscMemType *mtype)
2321: {
2322:   PetscBool isMPI;

2324:   PetscFunctionBegin;
2327:   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 */
2328:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2329:   if (isMPI) {
2330:     /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2331:     PetscCall(MatDenseGetArrayAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2332:   } else {
2333:     PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);

2335:     PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayAndMemType_C", &fptr));
2336:     if (fptr) {
2337:       PetscCall((*fptr)(A, array, mtype));
2338:     } else {
2339:       PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2340:       if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2341:     }
2342:   }
2343:   PetscFunctionReturn(PETSC_SUCCESS);
2344: }

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

2349:    Logically Collective

2351:    Input Parameters:
2352: +  mat - a dense matrix
2353: -  array - pointer to the data

2355:    Level: intermediate

2357: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2358: @*/
2359: PetscErrorCode MatDenseRestoreArrayAndMemType(Mat A, PetscScalar **array)
2360: {
2361:   PetscBool isMPI;

2363:   PetscFunctionBegin;
2366:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2367:   if (isMPI) {
2368:     PetscCall(MatDenseRestoreArrayAndMemType(((Mat_MPIDense *)A->data)->A, array));
2369:   } else {
2370:     PetscErrorCode (*fptr)(Mat, PetscScalar **);

2372:     PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayAndMemType_C", &fptr));
2373:     if (fptr) {
2374:       PetscCall((*fptr)(A, array));
2375:     } else {
2376:       PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2377:     }
2378:     *array = NULL;
2379:   }
2380:   PetscCall(PetscObjectStateIncrease((PetscObject)A));
2381:   PetscFunctionReturn(PETSC_SUCCESS);
2382: }

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

2387:    Logically Collective

2389:    Input Parameter:
2390: .  mat - a dense matrix

2392:    Output Parameters:
2393: +  array - pointer to the data
2394: -  mtype - memory type of the returned pointer

2396:    Level: intermediate

2398:    Notes:
2399:    If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2400:    an array on device is always returned and is guaranteed to contain the matrix's latest data.

2402: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`,
2403:    `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2404: @*/
2405: PetscErrorCode MatDenseGetArrayReadAndMemType(Mat A, const PetscScalar **array, PetscMemType *mtype)
2406: {
2407:   PetscBool isMPI;

2409:   PetscFunctionBegin;
2412:   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 */
2413:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2414:   if (isMPI) { /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2415:     PetscCall(MatDenseGetArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2416:   } else {
2417:     PetscErrorCode (*fptr)(Mat, const PetscScalar **, PetscMemType *);

2419:     PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayReadAndMemType_C", &fptr));
2420:     if (fptr) {
2421:       PetscCall((*fptr)(A, array, mtype));
2422:     } else {
2423:       PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, const PetscScalar **), (A, array));
2424:       if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2425:     }
2426:   }
2427:   PetscFunctionReturn(PETSC_SUCCESS);
2428: }

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

2433:    Logically Collective

2435:    Input Parameters:
2436: +  mat - a dense matrix
2437: -  array - pointer to the data

2439:    Level: intermediate

2441: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2442: @*/
2443: PetscErrorCode MatDenseRestoreArrayReadAndMemType(Mat A, const PetscScalar **array)
2444: {
2445:   PetscBool isMPI;

2447:   PetscFunctionBegin;
2450:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2451:   if (isMPI) {
2452:     PetscCall(MatDenseRestoreArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array));
2453:   } else {
2454:     PetscErrorCode (*fptr)(Mat, const PetscScalar **);

2456:     PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayReadAndMemType_C", &fptr));
2457:     if (fptr) {
2458:       PetscCall((*fptr)(A, array));
2459:     } else {
2460:       PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, const PetscScalar **), (A, array));
2461:     }
2462:     *array = NULL;
2463:   }
2464:   PetscFunctionReturn(PETSC_SUCCESS);
2465: }

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

2470:    Logically Collective

2472:    Input Parameter:
2473: .  mat - a dense matrix

2475:    Output Parameters:
2476: +  array - pointer to the data
2477: -  mtype - memory type of the returned pointer

2479:    Level: intermediate

2481:    Notes:
2482:    If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2483:    an array on device is always returned and is guaranteed to contain the matrix's latest data.

2485: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWriteAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayRead()`,
2486:   `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2487: @*/
2488: PetscErrorCode MatDenseGetArrayWriteAndMemType(Mat A, PetscScalar **array, PetscMemType *mtype)
2489: {
2490:   PetscBool isMPI;

2492:   PetscFunctionBegin;
2495:   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 */
2496:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2497:   if (isMPI) {
2498:     PetscCall(MatDenseGetArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2499:   } else {
2500:     PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);

2502:     PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayWriteAndMemType_C", &fptr));
2503:     if (fptr) {
2504:       PetscCall((*fptr)(A, array, mtype));
2505:     } else {
2506:       PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2507:       if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2508:     }
2509:   }
2510:   PetscFunctionReturn(PETSC_SUCCESS);
2511: }

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

2516:    Logically Collective

2518:    Input Parameters:
2519: +  mat - a dense matrix
2520: -  array - pointer to the data

2522:    Level: intermediate

2524: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2525: @*/
2526: PetscErrorCode MatDenseRestoreArrayWriteAndMemType(Mat A, PetscScalar **array)
2527: {
2528:   PetscBool isMPI;

2530:   PetscFunctionBegin;
2533:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2534:   if (isMPI) {
2535:     PetscCall(MatDenseRestoreArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array));
2536:   } else {
2537:     PetscErrorCode (*fptr)(Mat, PetscScalar **);

2539:     PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayWriteAndMemType_C", &fptr));
2540:     if (fptr) {
2541:       PetscCall((*fptr)(A, array));
2542:     } else {
2543:       PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2544:     }
2545:     *array = NULL;
2546:   }
2547:   PetscCall(PetscObjectStateIncrease((PetscObject)A));
2548:   PetscFunctionReturn(PETSC_SUCCESS);
2549: }

2551: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B)
2552: {
2553:   Mat_SeqDense   *mat = (Mat_SeqDense *)A->data;
2554:   PetscInt        i, j, nrows, ncols, ldb;
2555:   const PetscInt *irow, *icol;
2556:   PetscScalar    *av, *bv, *v = mat->v;
2557:   Mat             newmat;

2559:   PetscFunctionBegin;
2560:   PetscCall(ISGetIndices(isrow, &irow));
2561:   PetscCall(ISGetIndices(iscol, &icol));
2562:   PetscCall(ISGetLocalSize(isrow, &nrows));
2563:   PetscCall(ISGetLocalSize(iscol, &ncols));

2565:   /* Check submatrixcall */
2566:   if (scall == MAT_REUSE_MATRIX) {
2567:     PetscInt n_cols, n_rows;
2568:     PetscCall(MatGetSize(*B, &n_rows, &n_cols));
2569:     if (n_rows != nrows || n_cols != ncols) {
2570:       /* resize the result matrix to match number of requested rows/columns */
2571:       PetscCall(MatSetSizes(*B, nrows, ncols, nrows, ncols));
2572:     }
2573:     newmat = *B;
2574:   } else {
2575:     /* Create and fill new matrix */
2576:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &newmat));
2577:     PetscCall(MatSetSizes(newmat, nrows, ncols, nrows, ncols));
2578:     PetscCall(MatSetType(newmat, ((PetscObject)A)->type_name));
2579:     PetscCall(MatSeqDenseSetPreallocation(newmat, NULL));
2580:   }

2582:   /* Now extract the data pointers and do the copy,column at a time */
2583:   PetscCall(MatDenseGetArray(newmat, &bv));
2584:   PetscCall(MatDenseGetLDA(newmat, &ldb));
2585:   for (i = 0; i < ncols; i++) {
2586:     av = v + mat->lda * icol[i];
2587:     for (j = 0; j < nrows; j++) bv[j] = av[irow[j]];
2588:     bv += ldb;
2589:   }
2590:   PetscCall(MatDenseRestoreArray(newmat, &bv));

2592:   /* Assemble the matrices so that the correct flags are set */
2593:   PetscCall(MatAssemblyBegin(newmat, MAT_FINAL_ASSEMBLY));
2594:   PetscCall(MatAssemblyEnd(newmat, MAT_FINAL_ASSEMBLY));

2596:   /* Free work space */
2597:   PetscCall(ISRestoreIndices(isrow, &irow));
2598:   PetscCall(ISRestoreIndices(iscol, &icol));
2599:   *B = newmat;
2600:   PetscFunctionReturn(PETSC_SUCCESS);
2601: }

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

2607:   PetscFunctionBegin;
2608:   if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n, B));

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

2614: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat, MatAssemblyType mode)
2615: {
2616:   PetscFunctionBegin;
2617:   PetscFunctionReturn(PETSC_SUCCESS);
2618: }

2620: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat, MatAssemblyType mode)
2621: {
2622:   PetscFunctionBegin;
2623:   PetscFunctionReturn(PETSC_SUCCESS);
2624: }

2626: PetscErrorCode MatCopy_SeqDense(Mat A, Mat B, MatStructure str)
2627: {
2628:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data;
2629:   const PetscScalar *va;
2630:   PetscScalar       *vb;
2631:   PetscInt           lda1 = a->lda, lda2 = b->lda, m = A->rmap->n, n = A->cmap->n, j;

2633:   PetscFunctionBegin;
2634:   /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2635:   if (A->ops->copy != B->ops->copy) {
2636:     PetscCall(MatCopy_Basic(A, B, str));
2637:     PetscFunctionReturn(PETSC_SUCCESS);
2638:   }
2639:   PetscCheck(m == B->rmap->n && n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "size(B) != size(A)");
2640:   PetscCall(MatDenseGetArrayRead(A, &va));
2641:   PetscCall(MatDenseGetArray(B, &vb));
2642:   if (lda1 > m || lda2 > m) {
2643:     for (j = 0; j < n; j++) PetscCall(PetscArraycpy(vb + j * lda2, va + j * lda1, m));
2644:   } else {
2645:     PetscCall(PetscArraycpy(vb, va, A->rmap->n * A->cmap->n));
2646:   }
2647:   PetscCall(MatDenseRestoreArray(B, &vb));
2648:   PetscCall(MatDenseRestoreArrayRead(A, &va));
2649:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2650:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2651:   PetscFunctionReturn(PETSC_SUCCESS);
2652: }

2654: PetscErrorCode MatSetUp_SeqDense(Mat A)
2655: {
2656:   PetscFunctionBegin;
2657:   PetscCall(PetscLayoutSetUp(A->rmap));
2658:   PetscCall(PetscLayoutSetUp(A->cmap));
2659:   if (!A->preallocated) PetscCall(MatSeqDenseSetPreallocation(A, NULL));
2660:   PetscFunctionReturn(PETSC_SUCCESS);
2661: }

2663: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2664: {
2665:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2666:   PetscInt      i, j;
2667:   PetscInt      min = PetscMin(A->rmap->n, A->cmap->n);
2668:   PetscScalar  *aa;

2670:   PetscFunctionBegin;
2671:   PetscCall(MatDenseGetArray(A, &aa));
2672:   for (j = 0; j < A->cmap->n; j++) {
2673:     for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscConj(aa[i + j * mat->lda]);
2674:   }
2675:   PetscCall(MatDenseRestoreArray(A, &aa));
2676:   if (mat->tau)
2677:     for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
2678:   PetscFunctionReturn(PETSC_SUCCESS);
2679: }

2681: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2682: {
2683:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2684:   PetscInt      i, j;
2685:   PetscScalar  *aa;

2687:   PetscFunctionBegin;
2688:   PetscCall(MatDenseGetArray(A, &aa));
2689:   for (j = 0; j < A->cmap->n; j++) {
2690:     for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscRealPart(aa[i + j * mat->lda]);
2691:   }
2692:   PetscCall(MatDenseRestoreArray(A, &aa));
2693:   PetscFunctionReturn(PETSC_SUCCESS);
2694: }

2696: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2697: {
2698:   Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2699:   PetscInt      i, j;
2700:   PetscScalar  *aa;

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

2711: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2712: {
2713:   PetscInt  m = A->rmap->n, n = B->cmap->n;
2714:   PetscBool cisdense = PETSC_FALSE;

2716:   PetscFunctionBegin;
2717:   PetscCall(MatSetSizes(C, m, n, m, n));
2718: #if defined(PETSC_HAVE_CUDA)
2719:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2720: #endif
2721: #if defined(PETSC_HAVE_HIP)
2722:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2723: #endif
2724:   if (!cisdense) {
2725:     PetscBool flg;

2727:     PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2728:     PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2729:   }
2730:   PetscCall(MatSetUp(C));
2731:   PetscFunctionReturn(PETSC_SUCCESS);
2732: }

2734: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2735: {
2736:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data, *c = (Mat_SeqDense *)C->data;
2737:   PetscBLASInt       m, n, k;
2738:   const PetscScalar *av, *bv;
2739:   PetscScalar       *cv;
2740:   PetscScalar        _DOne = 1.0, _DZero = 0.0;

2742:   PetscFunctionBegin;
2743:   PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2744:   PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2745:   PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2746:   if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2747:   PetscCall(MatDenseGetArrayRead(A, &av));
2748:   PetscCall(MatDenseGetArrayRead(B, &bv));
2749:   PetscCall(MatDenseGetArrayWrite(C, &cv));
2750:   PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2751:   PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2752:   PetscCall(MatDenseRestoreArrayRead(A, &av));
2753:   PetscCall(MatDenseRestoreArrayRead(B, &bv));
2754:   PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2755:   PetscFunctionReturn(PETSC_SUCCESS);
2756: }

2758: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2759: {
2760:   PetscInt  m = A->rmap->n, n = B->rmap->n;
2761:   PetscBool cisdense = PETSC_FALSE;

2763:   PetscFunctionBegin;
2764:   PetscCall(MatSetSizes(C, m, n, m, n));
2765: #if defined(PETSC_HAVE_CUDA)
2766:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2767: #endif
2768: #if defined(PETSC_HAVE_HIP)
2769:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2770: #endif
2771:   if (!cisdense) {
2772:     PetscBool flg;

2774:     PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2775:     PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2776:   }
2777:   PetscCall(MatSetUp(C));
2778:   PetscFunctionReturn(PETSC_SUCCESS);
2779: }

2781: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2782: {
2783:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data;
2784:   Mat_SeqDense      *b = (Mat_SeqDense *)B->data;
2785:   Mat_SeqDense      *c = (Mat_SeqDense *)C->data;
2786:   const PetscScalar *av, *bv;
2787:   PetscScalar       *cv;
2788:   PetscBLASInt       m, n, k;
2789:   PetscScalar        _DOne = 1.0, _DZero = 0.0;

2791:   PetscFunctionBegin;
2792:   PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2793:   PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2794:   PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2795:   if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2796:   PetscCall(MatDenseGetArrayRead(A, &av));
2797:   PetscCall(MatDenseGetArrayRead(B, &bv));
2798:   PetscCall(MatDenseGetArrayWrite(C, &cv));
2799:   PetscCallBLAS("BLASgemm", BLASgemm_("N", "T", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2800:   PetscCall(MatDenseRestoreArrayRead(A, &av));
2801:   PetscCall(MatDenseRestoreArrayRead(B, &bv));
2802:   PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2803:   PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2804:   PetscFunctionReturn(PETSC_SUCCESS);
2805: }

2807: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2808: {
2809:   PetscInt  m = A->cmap->n, n = B->cmap->n;
2810:   PetscBool cisdense = PETSC_FALSE;

2812:   PetscFunctionBegin;
2813:   PetscCall(MatSetSizes(C, m, n, m, n));
2814: #if defined(PETSC_HAVE_CUDA)
2815:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2816: #endif
2817: #if defined(PETSC_HAVE_HIP)
2818:   PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2819: #endif
2820:   if (!cisdense) {
2821:     PetscBool flg;

2823:     PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2824:     PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2825:   }
2826:   PetscCall(MatSetUp(C));
2827:   PetscFunctionReturn(PETSC_SUCCESS);
2828: }

2830: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2831: {
2832:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data;
2833:   Mat_SeqDense      *b = (Mat_SeqDense *)B->data;
2834:   Mat_SeqDense      *c = (Mat_SeqDense *)C->data;
2835:   const PetscScalar *av, *bv;
2836:   PetscScalar       *cv;
2837:   PetscBLASInt       m, n, k;
2838:   PetscScalar        _DOne = 1.0, _DZero = 0.0;

2840:   PetscFunctionBegin;
2841:   PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2842:   PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2843:   PetscCall(PetscBLASIntCast(A->rmap->n, &k));
2844:   if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2845:   PetscCall(MatDenseGetArrayRead(A, &av));
2846:   PetscCall(MatDenseGetArrayRead(B, &bv));
2847:   PetscCall(MatDenseGetArrayWrite(C, &cv));
2848:   PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2849:   PetscCall(MatDenseRestoreArrayRead(A, &av));
2850:   PetscCall(MatDenseRestoreArrayRead(B, &bv));
2851:   PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2852:   PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2853:   PetscFunctionReturn(PETSC_SUCCESS);
2854: }

2856: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2857: {
2858:   PetscFunctionBegin;
2859:   C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2860:   C->ops->productsymbolic = MatProductSymbolic_AB;
2861:   PetscFunctionReturn(PETSC_SUCCESS);
2862: }

2864: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2865: {
2866:   PetscFunctionBegin;
2867:   C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2868:   C->ops->productsymbolic          = MatProductSymbolic_AtB;
2869:   PetscFunctionReturn(PETSC_SUCCESS);
2870: }

2872: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2873: {
2874:   PetscFunctionBegin;
2875:   C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2876:   C->ops->productsymbolic          = MatProductSymbolic_ABt;
2877:   PetscFunctionReturn(PETSC_SUCCESS);
2878: }

2880: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2881: {
2882:   Mat_Product *product = C->product;

2884:   PetscFunctionBegin;
2885:   switch (product->type) {
2886:   case MATPRODUCT_AB:
2887:     PetscCall(MatProductSetFromOptions_SeqDense_AB(C));
2888:     break;
2889:   case MATPRODUCT_AtB:
2890:     PetscCall(MatProductSetFromOptions_SeqDense_AtB(C));
2891:     break;
2892:   case MATPRODUCT_ABt:
2893:     PetscCall(MatProductSetFromOptions_SeqDense_ABt(C));
2894:     break;
2895:   default:
2896:     break;
2897:   }
2898:   PetscFunctionReturn(PETSC_SUCCESS);
2899: }

2901: static PetscErrorCode MatGetRowMax_SeqDense(Mat A, Vec v, PetscInt idx[])
2902: {
2903:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data;
2904:   PetscInt           i, j, m = A->rmap->n, n = A->cmap->n, p;
2905:   PetscScalar       *x;
2906:   const PetscScalar *aa;

2908:   PetscFunctionBegin;
2909:   PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2910:   PetscCall(VecGetArray(v, &x));
2911:   PetscCall(VecGetLocalSize(v, &p));
2912:   PetscCall(MatDenseGetArrayRead(A, &aa));
2913:   PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2914:   for (i = 0; i < m; i++) {
2915:     x[i] = aa[i];
2916:     if (idx) idx[i] = 0;
2917:     for (j = 1; j < n; j++) {
2918:       if (PetscRealPart(x[i]) < PetscRealPart(aa[i + a->lda * j])) {
2919:         x[i] = aa[i + a->lda * j];
2920:         if (idx) idx[i] = j;
2921:       }
2922:     }
2923:   }
2924:   PetscCall(MatDenseRestoreArrayRead(A, &aa));
2925:   PetscCall(VecRestoreArray(v, &x));
2926:   PetscFunctionReturn(PETSC_SUCCESS);
2927: }

2929: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A, Vec v, PetscInt idx[])
2930: {
2931:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data;
2932:   PetscInt           i, j, m = A->rmap->n, n = A->cmap->n, p;
2933:   PetscScalar       *x;
2934:   PetscReal          atmp;
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] = PetscAbsScalar(aa[i]);
2945:     for (j = 1; j < n; j++) {
2946:       atmp = PetscAbsScalar(aa[i + a->lda * j]);
2947:       if (PetscAbsScalar(x[i]) < atmp) {
2948:         x[i] = atmp;
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 MatGetRowMin_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:   const PetscScalar *aa;

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

2986: PetscErrorCode MatGetColumnVector_SeqDense(Mat A, Vec v, PetscInt col)
2987: {
2988:   Mat_SeqDense      *a = (Mat_SeqDense *)A->data;
2989:   PetscScalar       *x;
2990:   const PetscScalar *aa;

2992:   PetscFunctionBegin;
2993:   PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2994:   PetscCall(MatDenseGetArrayRead(A, &aa));
2995:   PetscCall(VecGetArray(v, &x));
2996:   PetscCall(PetscArraycpy(x, aa + col * a->lda, A->rmap->n));
2997:   PetscCall(VecRestoreArray(v, &x));
2998:   PetscCall(MatDenseRestoreArrayRead(A, &aa));
2999:   PetscFunctionReturn(PETSC_SUCCESS);
3000: }

3002: PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A, PetscInt type, PetscReal *reductions)
3003: {
3004:   PetscInt           i, j, m, n;
3005:   const PetscScalar *a;

3007:   PetscFunctionBegin;
3008:   PetscCall(MatGetSize(A, &m, &n));
3009:   PetscCall(PetscArrayzero(reductions, n));
3010:   PetscCall(MatDenseGetArrayRead(A, &a));
3011:   if (type == NORM_2) {
3012:     for (i = 0; i < n; i++) {
3013:       for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j] * a[j]);
3014:       a += m;
3015:     }
3016:   } else if (type == NORM_1) {
3017:     for (i = 0; i < n; i++) {
3018:       for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j]);
3019:       a += m;
3020:     }
3021:   } else if (type == NORM_INFINITY) {
3022:     for (i = 0; i < n; i++) {
3023:       for (j = 0; j < m; j++) reductions[i] = PetscMax(PetscAbsScalar(a[j]), reductions[i]);
3024:       a += m;
3025:     }
3026:   } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
3027:     for (i = 0; i < n; i++) {
3028:       for (j = 0; j < m; j++) reductions[i] += PetscRealPart(a[j]);
3029:       a += m;
3030:     }
3031:   } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3032:     for (i = 0; i < n; i++) {
3033:       for (j = 0; j < m; j++) reductions[i] += PetscImaginaryPart(a[j]);
3034:       a += m;
3035:     }
3036:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown reduction type");
3037:   PetscCall(MatDenseRestoreArrayRead(A, &a));
3038:   if (type == NORM_2) {
3039:     for (i = 0; i < n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
3040:   } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3041:     for (i = 0; i < n; i++) reductions[i] /= m;
3042:   }
3043:   PetscFunctionReturn(PETSC_SUCCESS);
3044: }

3046: PetscErrorCode MatSetRandom_SeqDense(Mat x, PetscRandom rctx)
3047: {
3048:   PetscScalar *a;
3049:   PetscInt     lda, m, n, i, j;

3051:   PetscFunctionBegin;
3052:   PetscCall(MatGetSize(x, &m, &n));
3053:   PetscCall(MatDenseGetLDA(x, &lda));
3054:   PetscCall(MatDenseGetArrayWrite(x, &a));
3055:   for (j = 0; j < n; j++) {
3056:     for (i = 0; i < m; i++) PetscCall(PetscRandomGetValue(rctx, a + j * lda + i));
3057:   }
3058:   PetscCall(MatDenseRestoreArrayWrite(x, &a));
3059:   PetscFunctionReturn(PETSC_SUCCESS);
3060: }

3062: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A, PetscBool *missing, PetscInt *d)
3063: {
3064:   PetscFunctionBegin;
3065:   *missing = PETSC_FALSE;
3066:   PetscFunctionReturn(PETSC_SUCCESS);
3067: }

3069: /* vals is not const */
3070: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A, PetscInt col, PetscScalar **vals)
3071: {
3072:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3073:   PetscScalar  *v;

3075:   PetscFunctionBegin;
3076:   PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3077:   PetscCall(MatDenseGetArray(A, &v));
3078:   *vals = v + col * a->lda;
3079:   PetscCall(MatDenseRestoreArray(A, &v));
3080:   PetscFunctionReturn(PETSC_SUCCESS);
3081: }

3083: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A, PetscScalar **vals)
3084: {
3085:   PetscFunctionBegin;
3086:   if (vals) *vals = NULL; /* user cannot accidentally use the array later */
3087:   PetscFunctionReturn(PETSC_SUCCESS);
3088: }

3090: static struct _MatOps MatOps_Values = {MatSetValues_SeqDense,
3091:                                        MatGetRow_SeqDense,
3092:                                        MatRestoreRow_SeqDense,
3093:                                        MatMult_SeqDense,
3094:                                        /*  4*/ MatMultAdd_SeqDense,
3095:                                        MatMultTranspose_SeqDense,
3096:                                        MatMultTransposeAdd_SeqDense,
3097:                                        NULL,
3098:                                        NULL,
3099:                                        NULL,
3100:                                        /* 10*/ NULL,
3101:                                        MatLUFactor_SeqDense,
3102:                                        MatCholeskyFactor_SeqDense,
3103:                                        MatSOR_SeqDense,
3104:                                        MatTranspose_SeqDense,
3105:                                        /* 15*/ MatGetInfo_SeqDense,
3106:                                        MatEqual_SeqDense,
3107:                                        MatGetDiagonal_SeqDense,
3108:                                        MatDiagonalScale_SeqDense,
3109:                                        MatNorm_SeqDense,
3110:                                        /* 20*/ MatAssemblyBegin_SeqDense,
3111:                                        MatAssemblyEnd_SeqDense,
3112:                                        MatSetOption_SeqDense,
3113:                                        MatZeroEntries_SeqDense,
3114:                                        /* 24*/ MatZeroRows_SeqDense,
3115:                                        NULL,
3116:                                        NULL,
3117:                                        NULL,
3118:                                        NULL,
3119:                                        /* 29*/ MatSetUp_SeqDense,
3120:                                        NULL,
3121:                                        NULL,
3122:                                        NULL,
3123:                                        NULL,
3124:                                        /* 34*/ MatDuplicate_SeqDense,
3125:                                        NULL,
3126:                                        NULL,
3127:                                        NULL,
3128:                                        NULL,
3129:                                        /* 39*/ MatAXPY_SeqDense,
3130:                                        MatCreateSubMatrices_SeqDense,
3131:                                        NULL,
3132:                                        MatGetValues_SeqDense,
3133:                                        MatCopy_SeqDense,
3134:                                        /* 44*/ MatGetRowMax_SeqDense,
3135:                                        MatScale_SeqDense,
3136:                                        MatShift_SeqDense,
3137:                                        NULL,
3138:                                        MatZeroRowsColumns_SeqDense,
3139:                                        /* 49*/ MatSetRandom_SeqDense,
3140:                                        NULL,
3141:                                        NULL,
3142:                                        NULL,
3143:                                        NULL,
3144:                                        /* 54*/ NULL,
3145:                                        NULL,
3146:                                        NULL,
3147:                                        NULL,
3148:                                        NULL,
3149:                                        /* 59*/ MatCreateSubMatrix_SeqDense,
3150:                                        MatDestroy_SeqDense,
3151:                                        MatView_SeqDense,
3152:                                        NULL,
3153:                                        NULL,
3154:                                        /* 64*/ NULL,
3155:                                        NULL,
3156:                                        NULL,
3157:                                        NULL,
3158:                                        NULL,
3159:                                        /* 69*/ MatGetRowMaxAbs_SeqDense,
3160:                                        NULL,
3161:                                        NULL,
3162:                                        NULL,
3163:                                        NULL,
3164:                                        /* 74*/ NULL,
3165:                                        NULL,
3166:                                        NULL,
3167:                                        NULL,
3168:                                        NULL,
3169:                                        /* 79*/ NULL,
3170:                                        NULL,
3171:                                        NULL,
3172:                                        NULL,
3173:                                        /* 83*/ MatLoad_SeqDense,
3174:                                        MatIsSymmetric_SeqDense,
3175:                                        MatIsHermitian_SeqDense,
3176:                                        NULL,
3177:                                        NULL,
3178:                                        NULL,
3179:                                        /* 89*/ NULL,
3180:                                        NULL,
3181:                                        MatMatMultNumeric_SeqDense_SeqDense,
3182:                                        NULL,
3183:                                        NULL,
3184:                                        /* 94*/ NULL,
3185:                                        NULL,
3186:                                        NULL,
3187:                                        MatMatTransposeMultNumeric_SeqDense_SeqDense,
3188:                                        NULL,
3189:                                        /* 99*/ MatProductSetFromOptions_SeqDense,
3190:                                        NULL,
3191:                                        NULL,
3192:                                        MatConjugate_SeqDense,
3193:                                        NULL,
3194:                                        /*104*/ NULL,
3195:                                        MatRealPart_SeqDense,
3196:                                        MatImaginaryPart_SeqDense,
3197:                                        NULL,
3198:                                        NULL,
3199:                                        /*109*/ NULL,
3200:                                        NULL,
3201:                                        MatGetRowMin_SeqDense,
3202:                                        MatGetColumnVector_SeqDense,
3203:                                        MatMissingDiagonal_SeqDense,
3204:                                        /*114*/ NULL,
3205:                                        NULL,
3206:                                        NULL,
3207:                                        NULL,
3208:                                        NULL,
3209:                                        /*119*/ NULL,
3210:                                        NULL,
3211:                                        NULL,
3212:                                        NULL,
3213:                                        NULL,
3214:                                        /*124*/ NULL,
3215:                                        MatGetColumnReductions_SeqDense,
3216:                                        NULL,
3217:                                        NULL,
3218:                                        NULL,
3219:                                        /*129*/ NULL,
3220:                                        NULL,
3221:                                        NULL,
3222:                                        MatTransposeMatMultNumeric_SeqDense_SeqDense,
3223:                                        NULL,
3224:                                        /*134*/ NULL,
3225:                                        NULL,
3226:                                        NULL,
3227:                                        NULL,
3228:                                        NULL,
3229:                                        /*139*/ NULL,
3230:                                        NULL,
3231:                                        NULL,
3232:                                        NULL,
3233:                                        NULL,
3234:                                        MatCreateMPIMatConcatenateSeqMat_SeqDense,
3235:                                        /*145*/ NULL,
3236:                                        NULL,
3237:                                        NULL,
3238:                                        NULL,
3239:                                        NULL,
3240:                                        /*150*/ NULL,
3241:                                        NULL};

3243: /*@C
3244:    MatCreateSeqDense - Creates a `MATSEQDENSE` that
3245:    is stored in column major order (the usual Fortran manner). Many
3246:    of the matrix operations use the BLAS and LAPACK routines.

3248:    Collective

3250:    Input Parameters:
3251: +  comm - MPI communicator, set to `PETSC_COMM_SELF`
3252: .  m - number of rows
3253: .  n - number of columns
3254: -  data - optional location of matrix data in column major order.  Use `NULL` for PETSc
3255:    to control all matrix memory allocation.

3257:    Output Parameter:
3258: .  A - the matrix

3260:    Level: intermediate

3262:    Note:
3263:    The data input variable is intended primarily for Fortran programmers
3264:    who wish to allocate their own matrix memory space.  Most users should
3265:    set `data` = `NULL`.

3267: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`
3268: @*/
3269: PetscErrorCode MatCreateSeqDense(MPI_Comm comm, PetscInt m, PetscInt n, PetscScalar *data, Mat *A)
3270: {
3271:   PetscFunctionBegin;
3272:   PetscCall(MatCreate(comm, A));
3273:   PetscCall(MatSetSizes(*A, m, n, m, n));
3274:   PetscCall(MatSetType(*A, MATSEQDENSE));
3275:   PetscCall(MatSeqDenseSetPreallocation(*A, data));
3276:   PetscFunctionReturn(PETSC_SUCCESS);
3277: }

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

3282:    Collective

3284:    Input Parameters:
3285: +  B - the matrix
3286: -  data - the array (or `NULL`)

3288:    Level: intermediate

3290:    Note:
3291:    The data input variable is intended primarily for Fortran programmers
3292:    who wish to allocate their own matrix memory space.  Most users should
3293:    need not call this routine.

3295: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`, `MatDenseSetLDA()`
3296: @*/
3297: PetscErrorCode MatSeqDenseSetPreallocation(Mat B, PetscScalar data[])
3298: {
3299:   PetscFunctionBegin;
3301:   PetscTryMethod(B, "MatSeqDenseSetPreallocation_C", (Mat, PetscScalar[]), (B, data));
3302:   PetscFunctionReturn(PETSC_SUCCESS);
3303: }

3305: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B, PetscScalar *data)
3306: {
3307:   Mat_SeqDense *b = (Mat_SeqDense *)B->data;

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

3313:   PetscCall(PetscLayoutSetUp(B->rmap));
3314:   PetscCall(PetscLayoutSetUp(B->cmap));

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

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

3322:     b->user_alloc = PETSC_FALSE;
3323:   } else { /* user-allocated storage */
3324:     if (!b->user_alloc) PetscCall(PetscFree(b->v));
3325:     b->v          = data;
3326:     b->user_alloc = PETSC_TRUE;
3327:   }
3328:   B->assembled = PETSC_TRUE;
3329:   PetscFunctionReturn(PETSC_SUCCESS);
3330: }

3332: #if defined(PETSC_HAVE_ELEMENTAL)
3333: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
3334: {
3335:   Mat                mat_elemental;
3336:   const PetscScalar *array;
3337:   PetscScalar       *v_colwise;
3338:   PetscInt           M = A->rmap->N, N = A->cmap->N, i, j, k, *rows, *cols;

3340:   PetscFunctionBegin;
3341:   PetscCall(PetscMalloc3(M * N, &v_colwise, M, &rows, N, &cols));
3342:   PetscCall(MatDenseGetArrayRead(A, &array));
3343:   /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
3344:   k = 0;
3345:   for (j = 0; j < N; j++) {
3346:     cols[j] = j;
3347:     for (i = 0; i < M; i++) v_colwise[j * M + i] = array[k++];
3348:   }
3349:   for (i = 0; i < M; i++) rows[i] = i;
3350:   PetscCall(MatDenseRestoreArrayRead(A, &array));

3352:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
3353:   PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
3354:   PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
3355:   PetscCall(MatSetUp(mat_elemental));

3357:   /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
3358:   PetscCall(MatSetValues(mat_elemental, M, rows, N, cols, v_colwise, ADD_VALUES));
3359:   PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
3360:   PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
3361:   PetscCall(PetscFree3(v_colwise, rows, cols));

3363:   if (reuse == MAT_INPLACE_MATRIX) {
3364:     PetscCall(MatHeaderReplace(A, &mat_elemental));
3365:   } else {
3366:     *newmat = mat_elemental;
3367:   }
3368:   PetscFunctionReturn(PETSC_SUCCESS);
3369: }
3370: #endif

3372: PetscErrorCode MatDenseSetLDA_SeqDense(Mat B, PetscInt lda)
3373: {
3374:   Mat_SeqDense *b = (Mat_SeqDense *)B->data;
3375:   PetscBool     data;

3377:   PetscFunctionBegin;
3378:   data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
3379:   PetscCheck(b->user_alloc || !data || b->lda == lda, PETSC_COMM_SELF, PETSC_ERR_ORDER, "LDA cannot be changed after allocation of internal storage");
3380:   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);
3381:   b->lda = lda;
3382:   PetscFunctionReturn(PETSC_SUCCESS);
3383: }

3385: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
3386: {
3387:   PetscFunctionBegin;
3388:   PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIDense(comm, inmat, n, scall, outmat));
3389:   PetscFunctionReturn(PETSC_SUCCESS);
3390: }

3392: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3393: {
3394:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3396:   PetscFunctionBegin;
3397:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3398:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3399:   if (!a->cvec) { PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, &a->cvec)); }
3400:   a->vecinuse = col + 1;
3401:   PetscCall(MatDenseGetArray(A, (PetscScalar **)&a->ptrinuse));
3402:   PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda));
3403:   *v = a->cvec;
3404:   PetscFunctionReturn(PETSC_SUCCESS);
3405: }

3407: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3408: {
3409:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3411:   PetscFunctionBegin;
3412:   PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3413:   PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3414:   a->vecinuse = 0;
3415:   PetscCall(MatDenseRestoreArray(A, (PetscScalar **)&a->ptrinuse));
3416:   PetscCall(VecResetArray(a->cvec));
3417:   if (v) *v = NULL;
3418:   PetscFunctionReturn(PETSC_SUCCESS);
3419: }

3421: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3422: {
3423:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3425:   PetscFunctionBegin;
3426:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3427:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3428:   if (!a->cvec) { PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, &a->cvec)); }
3429:   a->vecinuse = col + 1;
3430:   PetscCall(MatDenseGetArrayRead(A, &a->ptrinuse));
3431:   PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda));
3432:   PetscCall(VecLockReadPush(a->cvec));
3433:   *v = a->cvec;
3434:   PetscFunctionReturn(PETSC_SUCCESS);
3435: }

3437: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3438: {
3439:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3441:   PetscFunctionBegin;
3442:   PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3443:   PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3444:   a->vecinuse = 0;
3445:   PetscCall(MatDenseRestoreArrayRead(A, &a->ptrinuse));
3446:   PetscCall(VecLockReadPop(a->cvec));
3447:   PetscCall(VecResetArray(a->cvec));
3448:   if (v) *v = NULL;
3449:   PetscFunctionReturn(PETSC_SUCCESS);
3450: }

3452: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3453: {
3454:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3456:   PetscFunctionBegin;
3457:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3458:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3459:   if (!a->cvec) PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, &a->cvec));
3460:   a->vecinuse = col + 1;
3461:   PetscCall(MatDenseGetArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3462:   PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda));
3463:   *v = a->cvec;
3464:   PetscFunctionReturn(PETSC_SUCCESS);
3465: }

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

3471:   PetscFunctionBegin;
3472:   PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3473:   PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3474:   a->vecinuse = 0;
3475:   PetscCall(MatDenseRestoreArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3476:   PetscCall(VecResetArray(a->cvec));
3477:   if (v) *v = NULL;
3478:   PetscFunctionReturn(PETSC_SUCCESS);
3479: }

3481: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3482: {
3483:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3485:   PetscFunctionBegin;
3486:   PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3487:   PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3488:   if (a->cmat && (cend - cbegin != a->cmat->cmap->N || rend - rbegin != a->cmat->rmap->N)) PetscCall(MatDestroy(&a->cmat));
3489:   if (!a->cmat) {
3490:     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), rend - rbegin, PETSC_DECIDE, rend - rbegin, cend - cbegin, a->v + rbegin + (size_t)cbegin * a->lda, &a->cmat));
3491:   } else {
3492:     PetscCall(MatDensePlaceArray(a->cmat, a->v + rbegin + (size_t)cbegin * a->lda));
3493:   }
3494:   PetscCall(MatDenseSetLDA(a->cmat, a->lda));
3495:   a->matinuse = cbegin + 1;
3496:   *v          = a->cmat;
3497: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
3498:   A->offloadmask = PETSC_OFFLOAD_CPU;
3499: #endif
3500:   PetscFunctionReturn(PETSC_SUCCESS);
3501: }

3503: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A, Mat *v)
3504: {
3505:   Mat_SeqDense *a = (Mat_SeqDense *)A->data;

3507:   PetscFunctionBegin;
3508:   PetscCheck(a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetSubMatrix() first");
3509:   PetscCheck(a->cmat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column matrix");
3510:   PetscCheck(*v == a->cmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not the matrix obtained from MatDenseGetSubMatrix()");
3511:   a->matinuse = 0;
3512:   PetscCall(MatDenseResetArray(a->cmat));
3513:   if (v) *v = NULL;
3514: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
3515:   A->offloadmask = PETSC_OFFLOAD_CPU;
3516: #endif
3517:   PetscFunctionReturn(PETSC_SUCCESS);
3518: }

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

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

3526:   Level: beginner

3528: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreateSeqDense()`
3529: M*/
3530: PetscErrorCode MatCreate_SeqDense(Mat B)
3531: {
3532:   Mat_SeqDense *b;
3533:   PetscMPIInt   size;

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

3539:   PetscCall(PetscNew(&b));
3540:   PetscCall(PetscMemcpy(B->ops, &MatOps_Values, sizeof(struct _MatOps)));
3541:   B->data = (void *)b;

3543:   b->roworiented = PETSC_TRUE;

3545:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactor_C", MatQRFactor_SeqDense));
3546:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetLDA_C", MatDenseGetLDA_SeqDense));
3547:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseSetLDA_C", MatDenseSetLDA_SeqDense));
3548:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArray_C", MatDenseGetArray_SeqDense));
3549:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArray_C", MatDenseRestoreArray_SeqDense));
3550:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDensePlaceArray_C", MatDensePlaceArray_SeqDense));
3551:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseResetArray_C", MatDenseResetArray_SeqDense));
3552:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseReplaceArray_C", MatDenseReplaceArray_SeqDense));
3553:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayRead_C", MatDenseGetArray_SeqDense));
3554:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayRead_C", MatDenseRestoreArray_SeqDense));
3555:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayWrite_C", MatDenseGetArray_SeqDense));
3556:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayWrite_C", MatDenseRestoreArray_SeqDense));
3557:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqaij_C", MatConvert_SeqDense_SeqAIJ));
3558: #if defined(PETSC_HAVE_ELEMENTAL)
3559:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_elemental_C", MatConvert_SeqDense_Elemental));
3560: #endif
3561: #if defined(PETSC_HAVE_SCALAPACK)
3562:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_scalapack_C", MatConvert_Dense_ScaLAPACK));
3563: #endif
3564: #if defined(PETSC_HAVE_CUDA)
3565:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensecuda_C", MatConvert_SeqDense_SeqDenseCUDA));
3566:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3567:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdense_C", MatProductSetFromOptions_SeqDense));
3568:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3569: #endif
3570: #if defined(PETSC_HAVE_HIP)
3571:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensehip_C", MatConvert_SeqDense_SeqDenseHIP));
3572:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3573:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdense_C", MatProductSetFromOptions_SeqDense));
3574:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3575: #endif
3576:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqDenseSetPreallocation_C", MatSeqDenseSetPreallocation_SeqDense));
3577:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqdense_C", MatProductSetFromOptions_SeqAIJ_SeqDense));
3578:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdense_C", MatProductSetFromOptions_SeqDense));
3579:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));
3580:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqsbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));

3582:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumn_C", MatDenseGetColumn_SeqDense));
3583:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumn_C", MatDenseRestoreColumn_SeqDense));
3584:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVec_C", MatDenseGetColumnVec_SeqDense));
3585:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVec_C", MatDenseRestoreColumnVec_SeqDense));
3586:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecRead_C", MatDenseGetColumnVecRead_SeqDense));
3587:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecRead_C", MatDenseRestoreColumnVecRead_SeqDense));
3588:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecWrite_C", MatDenseGetColumnVecWrite_SeqDense));
3589:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecWrite_C", MatDenseRestoreColumnVecWrite_SeqDense));
3590:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetSubMatrix_C", MatDenseGetSubMatrix_SeqDense));
3591:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreSubMatrix_C", MatDenseRestoreSubMatrix_SeqDense));
3592:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQDENSE));
3593:   PetscFunctionReturn(PETSC_SUCCESS);
3594: }

3596: /*@C
3597:    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.

3599:    Not Collective

3601:    Input Parameters:
3602: +  mat - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3603: -  col - column index

3605:    Output Parameter:
3606: .  vals - pointer to the data

3608:    Level: intermediate

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

3613: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreColumn()`, `MatDenseGetColumnVec()`
3614: @*/
3615: PetscErrorCode MatDenseGetColumn(Mat A, PetscInt col, PetscScalar **vals)
3616: {
3617:   PetscFunctionBegin;
3621:   PetscUseMethod(A, "MatDenseGetColumn_C", (Mat, PetscInt, PetscScalar **), (A, col, vals));
3622:   PetscFunctionReturn(PETSC_SUCCESS);
3623: }

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

3628:    Not Collective

3630:    Input Parameters:
3631: +  mat - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3632: -  vals - pointer to the data (may be `NULL`)

3634:    Level: intermediate

3636: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetColumn()`
3637: @*/
3638: PetscErrorCode MatDenseRestoreColumn(Mat A, PetscScalar **vals)
3639: {
3640:   PetscFunctionBegin;
3643:   PetscUseMethod(A, "MatDenseRestoreColumn_C", (Mat, PetscScalar **), (A, vals));
3644:   PetscFunctionReturn(PETSC_SUCCESS);
3645: }

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

3650:    Collective

3652:    Input Parameters:
3653: +  mat - the `Mat` object
3654: -  col - the column index

3656:    Output Parameter:
3657: .  v - the vector

3659:    Level: intermediate

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

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

3666: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`, `MatDenseGetColumn()`
3667: @*/
3668: PetscErrorCode MatDenseGetColumnVec(Mat A, PetscInt col, Vec *v)
3669: {
3670:   PetscFunctionBegin;
3675:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3676:   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);
3677:   PetscUseMethod(A, "MatDenseGetColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3678:   PetscFunctionReturn(PETSC_SUCCESS);
3679: }

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

3684:    Collective

3686:    Input Parameters:
3687: +  mat - the Mat object
3688: .  col - the column index
3689: -  v - the Vec object (may be `NULL`)

3691:    Level: intermediate

3693: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3694: @*/
3695: PetscErrorCode MatDenseRestoreColumnVec(Mat A, PetscInt col, Vec *v)
3696: {
3697:   PetscFunctionBegin;
3701:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3702:   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);
3703:   PetscUseMethod(A, "MatDenseRestoreColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3704:   PetscFunctionReturn(PETSC_SUCCESS);
3705: }

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

3710:    Collective

3712:    Input Parameters:
3713: +  mat - the `Mat` object
3714: -  col - the column index

3716:    Output Parameter:
3717: .  v - the vector

3719:    Level: intermediate

3721:    Notes:
3722:      The vector is owned by PETSc and users cannot modify it.

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

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

3728: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3729: @*/
3730: PetscErrorCode MatDenseGetColumnVecRead(Mat A, PetscInt col, Vec *v)
3731: {
3732:   PetscFunctionBegin;
3737:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3738:   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);
3739:   PetscUseMethod(A, "MatDenseGetColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3740:   PetscFunctionReturn(PETSC_SUCCESS);
3741: }

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

3746:    Collective

3748:    Input Parameters:
3749: +  mat - the `Mat` object
3750: .  col - the column index
3751: -  v - the Vec object (may be `NULL`)

3753:    Level: intermediate

3755: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecWrite()`
3756: @*/
3757: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A, PetscInt col, Vec *v)
3758: {
3759:   PetscFunctionBegin;
3763:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3764:   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);
3765:   PetscUseMethod(A, "MatDenseRestoreColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3766:   PetscFunctionReturn(PETSC_SUCCESS);
3767: }

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

3772:    Collective

3774:    Input Parameters:
3775: +  mat - the `Mat` object
3776: -  col - the column index

3778:    Output Parameter:
3779: .  v - the vector

3781:    Level: intermediate

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

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

3788: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3789: @*/
3790: PetscErrorCode MatDenseGetColumnVecWrite(Mat A, PetscInt col, Vec *v)
3791: {
3792:   PetscFunctionBegin;
3797:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3798:   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);
3799:   PetscUseMethod(A, "MatDenseGetColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3800:   PetscFunctionReturn(PETSC_SUCCESS);
3801: }

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

3806:    Collective

3808:    Input Parameters:
3809: +  mat - the `Mat` object
3810: .  col - the column index
3811: -  v - the `Vec` object (may be `NULL`)

3813:    Level: intermediate

3815: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`
3816: @*/
3817: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A, PetscInt col, Vec *v)
3818: {
3819:   PetscFunctionBegin;
3823:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3824:   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);
3825:   PetscUseMethod(A, "MatDenseRestoreColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3826:   PetscFunctionReturn(PETSC_SUCCESS);
3827: }

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

3832:    Collective

3834:    Input Parameters:
3835: +  mat - the Mat object
3836: .  rbegin - the first global row index in the block (if `PETSC_DECIDE`, is 0)
3837: .  rend - the global row index past the last one in the block (if `PETSC_DECIDE`, is `M`)
3838: .  cbegin - the first global column index in the block (if `PETSC_DECIDE`, is 0)
3839: -  cend - the global column index past the last one in the block (if `PETSC_DECIDE`, is `N`)

3841:    Output Parameter:
3842: .  v - the matrix

3844:    Level: intermediate

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

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

3851: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreSubMatrix()`
3852: @*/
3853: PetscErrorCode MatDenseGetSubMatrix(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3854: {
3855:   PetscFunctionBegin;
3863:   if (rbegin == PETSC_DECIDE) rbegin = 0;
3864:   if (rend == PETSC_DECIDE) rend = A->rmap->N;
3865:   if (cbegin == PETSC_DECIDE) cbegin = 0;
3866:   if (cend == PETSC_DECIDE) cend = A->cmap->N;
3867:   PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3868:   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);
3869:   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);
3870:   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);
3871:   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);
3872:   PetscUseMethod(A, "MatDenseGetSubMatrix_C", (Mat, PetscInt, PetscInt, PetscInt, PetscInt, Mat *), (A, rbegin, rend, cbegin, cend, v));
3873:   PetscFunctionReturn(PETSC_SUCCESS);
3874: }

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

3879:    Collective

3881:    Input Parameters:
3882: +  mat - the `Mat` object
3883: -  v - the `Mat` object (may be `NULL`)

3885:    Level: intermediate

3887: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseGetSubMatrix()`
3888: @*/
3889: PetscErrorCode MatDenseRestoreSubMatrix(Mat A, Mat *v)
3890: {
3891:   PetscFunctionBegin;
3895:   PetscUseMethod(A, "MatDenseRestoreSubMatrix_C", (Mat, Mat *), (A, v));
3896:   PetscFunctionReturn(PETSC_SUCCESS);
3897: }

3899: #include <petscblaslapack.h>
3900: #include <petsc/private/kernels/blockinvert.h>

3902: PetscErrorCode MatSeqDenseInvert(Mat A)
3903: {
3904:   PetscInt        m;
3905:   const PetscReal shift = 0.0;
3906:   PetscBool       allowzeropivot, zeropivotdetected = PETSC_FALSE;
3907:   PetscScalar    *values;

3909:   PetscFunctionBegin;
3911:   PetscCall(MatDenseGetArray(A, &values));
3912:   PetscCall(MatGetLocalSize(A, &m, NULL));
3913:   allowzeropivot = PetscNot(A->erroriffailure);
3914:   /* factor and invert each block */
3915:   switch (m) {
3916:   case 1:
3917:     values[0] = (PetscScalar)1.0 / (values[0] + shift);
3918:     break;
3919:   case 2:
3920:     PetscCall(PetscKernel_A_gets_inverse_A_2(values, shift, allowzeropivot, &zeropivotdetected));
3921:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3922:     break;
3923:   case 3:
3924:     PetscCall(PetscKernel_A_gets_inverse_A_3(values, shift, allowzeropivot, &zeropivotdetected));
3925:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3926:     break;
3927:   case 4:
3928:     PetscCall(PetscKernel_A_gets_inverse_A_4(values, shift, allowzeropivot, &zeropivotdetected));
3929:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3930:     break;
3931:   case 5: {
3932:     PetscScalar work[25];
3933:     PetscInt    ipvt[5];

3935:     PetscCall(PetscKernel_A_gets_inverse_A_5(values, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
3936:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3937:   } break;
3938:   case 6:
3939:     PetscCall(PetscKernel_A_gets_inverse_A_6(values, shift, allowzeropivot, &zeropivotdetected));
3940:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3941:     break;
3942:   case 7:
3943:     PetscCall(PetscKernel_A_gets_inverse_A_7(values, shift, allowzeropivot, &zeropivotdetected));
3944:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3945:     break;
3946:   default: {
3947:     PetscInt    *v_pivots, *IJ, j;
3948:     PetscScalar *v_work;

3950:     PetscCall(PetscMalloc3(m, &v_work, m, &v_pivots, m, &IJ));
3951:     for (j = 0; j < m; j++) IJ[j] = j;
3952:     PetscCall(PetscKernel_A_gets_inverse_A(m, values, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
3953:     if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3954:     PetscCall(PetscFree3(v_work, v_pivots, IJ));
3955:   }
3956:   }
3957:   PetscCall(MatDenseRestoreArray(A, &values));
3958:   PetscFunctionReturn(PETSC_SUCCESS);
3959: }