Actual source code: dense.c
1: /*
2: Defines the basic matrix operations for sequential dense.
3: Portions of this code are under:
4: Copyright (c) 2022 Advanced Micro Devices, Inc. All rights reserved.
5: */
7: #include <../src/mat/impls/dense/seq/dense.h>
8: #include <../src/mat/impls/dense/mpi/mpidense.h>
9: #include <petscblaslapack.h>
10: #include <../src/mat/impls/aij/seq/aij.h>
11: #include <petsc/private/vecimpl.h>
13: PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian)
14: {
15: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
16: PetscInt j, k, n = A->rmap->n;
17: PetscScalar *v;
19: PetscFunctionBegin;
20: PetscCheck(A->rmap->n == A->cmap->n, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot symmetrize a rectangular matrix");
21: PetscCall(MatDenseGetArray(A, &v));
22: if (!hermitian) {
23: for (k = 0; k < n; k++) {
24: for (j = k; j < n; j++) v[j * mat->lda + k] = v[k * mat->lda + j];
25: }
26: } else {
27: for (k = 0; k < n; k++) {
28: for (j = k; j < n; j++) v[j * mat->lda + k] = PetscConj(v[k * mat->lda + j]);
29: }
30: }
31: PetscCall(MatDenseRestoreArray(A, &v));
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A)
36: {
37: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
38: PetscBLASInt info, n;
40: PetscFunctionBegin;
41: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
42: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
43: if (A->factortype == MAT_FACTOR_LU) {
44: PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present");
45: if (!mat->fwork) {
46: mat->lfwork = n;
47: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
48: }
49: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
50: PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
51: PetscCall(PetscFPTrapPop());
52: PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
53: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
54: if (A->spd == PETSC_BOOL3_TRUE) {
55: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
56: PetscCallBLAS("LAPACKpotri", LAPACKpotri_("L", &n, mat->v, &mat->lda, &info));
57: PetscCall(PetscFPTrapPop());
58: PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_TRUE));
59: #if defined(PETSC_USE_COMPLEX)
60: } else if (A->hermitian == PETSC_BOOL3_TRUE) {
61: PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present");
62: PetscCheck(mat->fwork, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Fwork not present");
63: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
64: PetscCallBLAS("LAPACKhetri", LAPACKhetri_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &info));
65: PetscCall(PetscFPTrapPop());
66: PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_TRUE));
67: #endif
68: } else { /* symmetric case */
69: PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present");
70: PetscCheck(mat->fwork, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Fwork not present");
71: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
72: PetscCallBLAS("LAPACKsytri", LAPACKsytri_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &info));
73: PetscCall(PetscFPTrapPop());
74: PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_FALSE));
75: }
76: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad Inversion: zero pivot in row %" PetscBLASInt_FMT, info - 1);
77: PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
78: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix must be factored to solve");
80: A->ops->solve = NULL;
81: A->ops->matsolve = NULL;
82: A->ops->solvetranspose = NULL;
83: A->ops->matsolvetranspose = NULL;
84: A->ops->solveadd = NULL;
85: A->ops->solvetransposeadd = NULL;
86: A->factortype = MAT_FACTOR_NONE;
87: PetscCall(PetscFree(A->solvertype));
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
91: static PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
92: {
93: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
94: PetscInt m = l->lda, n = A->cmap->n, r = A->rmap->n, i, j;
95: PetscScalar *slot, *bb, *v;
96: const PetscScalar *xx;
98: PetscFunctionBegin;
99: if (PetscDefined(USE_DEBUG)) {
100: for (i = 0; i < N; i++) {
101: PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
102: PetscCheck(rows[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " requested to be zeroed greater than or equal number of rows %" PetscInt_FMT, rows[i], A->rmap->n);
103: PetscCheck(rows[i] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Col %" PetscInt_FMT " requested to be zeroed greater than or equal number of cols %" PetscInt_FMT, rows[i], A->cmap->n);
104: }
105: }
106: if (!N) PetscFunctionReturn(PETSC_SUCCESS);
108: /* fix right-hand side if needed */
109: if (x && b) {
110: Vec xt;
112: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
113: PetscCall(VecDuplicate(x, &xt));
114: PetscCall(VecCopy(x, xt));
115: PetscCall(VecScale(xt, -1.0));
116: PetscCall(MatMultAdd(A, xt, b, b));
117: PetscCall(VecDestroy(&xt));
118: PetscCall(VecGetArrayRead(x, &xx));
119: PetscCall(VecGetArray(b, &bb));
120: for (i = 0; i < N; i++) bb[rows[i]] = diag * xx[rows[i]];
121: PetscCall(VecRestoreArrayRead(x, &xx));
122: PetscCall(VecRestoreArray(b, &bb));
123: }
125: PetscCall(MatDenseGetArray(A, &v));
126: for (i = 0; i < N; i++) {
127: slot = v + rows[i] * m;
128: PetscCall(PetscArrayzero(slot, r));
129: }
130: for (i = 0; i < N; i++) {
131: slot = v + rows[i];
132: for (j = 0; j < n; j++) {
133: *slot = 0.0;
134: slot += m;
135: }
136: }
137: if (diag != 0.0) {
138: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
139: for (i = 0; i < N; i++) {
140: slot = v + (m + 1) * rows[i];
141: *slot = diag;
142: }
143: }
144: PetscCall(MatDenseRestoreArray(A, &v));
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
149: {
150: Mat B = NULL;
151: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
152: Mat_SeqDense *b;
153: PetscInt *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->N, i;
154: const MatScalar *av;
155: PetscBool isseqdense;
157: PetscFunctionBegin;
158: if (reuse == MAT_REUSE_MATRIX) {
159: PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATSEQDENSE, &isseqdense));
160: PetscCheck(isseqdense, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_USER, "Cannot reuse matrix of type %s", ((PetscObject)*newmat)->type_name);
161: }
162: if (reuse != MAT_REUSE_MATRIX) {
163: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
164: PetscCall(MatSetSizes(B, m, n, m, n));
165: PetscCall(MatSetType(B, MATSEQDENSE));
166: PetscCall(MatSeqDenseSetPreallocation(B, NULL));
167: b = (Mat_SeqDense *)B->data;
168: } else {
169: b = (Mat_SeqDense *)((*newmat)->data);
170: for (i = 0; i < n; i++) PetscCall(PetscArrayzero(b->v + i * b->lda, m));
171: }
172: PetscCall(MatSeqAIJGetArrayRead(A, &av));
173: for (i = 0; i < m; i++) {
174: PetscInt j;
175: for (j = 0; j < ai[1] - ai[0]; j++) {
176: b->v[*aj * b->lda + i] = *av;
177: aj++;
178: av++;
179: }
180: ai++;
181: }
182: PetscCall(MatSeqAIJRestoreArrayRead(A, &av));
184: if (reuse == MAT_INPLACE_MATRIX) {
185: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
186: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
187: PetscCall(MatHeaderReplace(A, &B));
188: } else {
189: if (B) *newmat = B;
190: PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
191: PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
192: }
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
197: {
198: Mat B = NULL;
199: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
200: PetscInt i, j;
201: PetscInt *rows, *nnz;
202: MatScalar *aa = a->v, *vals;
204: PetscFunctionBegin;
205: PetscCall(PetscCalloc3(A->rmap->n, &rows, A->rmap->n, &nnz, A->rmap->n, &vals));
206: if (reuse != MAT_REUSE_MATRIX) {
207: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
208: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
209: PetscCall(MatSetType(B, MATSEQAIJ));
210: for (j = 0; j < A->cmap->n; j++) {
211: for (i = 0; i < A->rmap->n; i++)
212: if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) ++nnz[i];
213: aa += a->lda;
214: }
215: PetscCall(MatSeqAIJSetPreallocation(B, PETSC_DETERMINE, nnz));
216: } else B = *newmat;
217: aa = a->v;
218: for (j = 0; j < A->cmap->n; j++) {
219: PetscInt numRows = 0;
220: for (i = 0; i < A->rmap->n; i++)
221: if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) {
222: rows[numRows] = i;
223: vals[numRows++] = aa[i];
224: }
225: PetscCall(MatSetValues(B, numRows, rows, 1, &j, vals, INSERT_VALUES));
226: aa += a->lda;
227: }
228: PetscCall(PetscFree3(rows, nnz, vals));
229: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
230: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
232: if (reuse == MAT_INPLACE_MATRIX) {
233: PetscCall(MatHeaderReplace(A, &B));
234: } else if (reuse != MAT_REUSE_MATRIX) *newmat = B;
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: PetscErrorCode MatAXPY_SeqDense(Mat Y, PetscScalar alpha, Mat X, MatStructure str)
239: {
240: Mat_SeqDense *x = (Mat_SeqDense *)X->data, *y = (Mat_SeqDense *)Y->data;
241: const PetscScalar *xv;
242: PetscScalar *yv;
243: PetscBLASInt N, m, ldax = 0, lday = 0, one = 1;
245: PetscFunctionBegin;
246: PetscCall(MatDenseGetArrayRead(X, &xv));
247: PetscCall(MatDenseGetArray(Y, &yv));
248: PetscCall(PetscBLASIntCast(X->rmap->n * X->cmap->n, &N));
249: PetscCall(PetscBLASIntCast(X->rmap->n, &m));
250: PetscCall(PetscBLASIntCast(x->lda, &ldax));
251: PetscCall(PetscBLASIntCast(y->lda, &lday));
252: if (ldax > m || lday > m) {
253: for (PetscInt j = 0; j < X->cmap->n; j++) PetscCallBLAS("BLASaxpy", BLASaxpy_(&m, &alpha, PetscSafePointerPlusOffset(xv, j * ldax), &one, PetscSafePointerPlusOffset(yv, j * lday), &one));
254: } else {
255: PetscCallBLAS("BLASaxpy", BLASaxpy_(&N, &alpha, xv, &one, yv, &one));
256: }
257: PetscCall(MatDenseRestoreArrayRead(X, &xv));
258: PetscCall(MatDenseRestoreArray(Y, &yv));
259: PetscCall(PetscLogFlops(PetscMax(2.0 * N - 1, 0)));
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: static PetscErrorCode MatGetInfo_SeqDense(Mat A, MatInfoType flag, MatInfo *info)
264: {
265: PetscLogDouble N = A->rmap->n * A->cmap->n;
267: PetscFunctionBegin;
268: info->block_size = 1.0;
269: info->nz_allocated = N;
270: info->nz_used = N;
271: info->nz_unneeded = 0;
272: info->assemblies = A->num_ass;
273: info->mallocs = 0;
274: info->memory = 0; /* REVIEW ME */
275: info->fill_ratio_given = 0;
276: info->fill_ratio_needed = 0;
277: info->factor_mallocs = 0;
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: PetscErrorCode MatScale_SeqDense(Mat A, PetscScalar alpha)
282: {
283: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
284: PetscScalar *v;
285: PetscBLASInt one = 1, j, nz, lda = 0;
287: PetscFunctionBegin;
288: PetscCall(MatDenseGetArray(A, &v));
289: PetscCall(PetscBLASIntCast(a->lda, &lda));
290: if (lda > A->rmap->n) {
291: PetscCall(PetscBLASIntCast(A->rmap->n, &nz));
292: for (j = 0; j < A->cmap->n; j++) PetscCallBLAS("BLASscal", BLASscal_(&nz, &alpha, v + j * lda, &one));
293: } else {
294: PetscCall(PetscBLASIntCast(A->rmap->n * A->cmap->n, &nz));
295: PetscCallBLAS("BLASscal", BLASscal_(&nz, &alpha, v, &one));
296: }
297: PetscCall(PetscLogFlops(A->rmap->n * A->cmap->n));
298: PetscCall(MatDenseRestoreArray(A, &v));
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
302: PetscErrorCode MatShift_SeqDense(Mat A, PetscScalar alpha)
303: {
304: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
305: PetscScalar *v;
306: PetscInt j, k;
308: PetscFunctionBegin;
309: PetscCall(MatDenseGetArray(A, &v));
310: k = PetscMin(A->rmap->n, A->cmap->n);
311: for (j = 0; j < k; j++) v[j + j * a->lda] += alpha;
312: PetscCall(PetscLogFlops(k));
313: PetscCall(MatDenseRestoreArray(A, &v));
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: static PetscErrorCode MatIsHermitian_SeqDense(Mat A, PetscReal rtol, PetscBool *fl)
318: {
319: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
320: PetscInt i, j, m = A->rmap->n, N = a->lda;
321: const PetscScalar *v;
323: PetscFunctionBegin;
324: *fl = PETSC_FALSE;
325: if (A->rmap->n != A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
326: PetscCall(MatDenseGetArrayRead(A, &v));
327: for (i = 0; i < m; i++) {
328: for (j = i; j < m; j++) {
329: if (PetscAbsScalar(v[i + j * N] - PetscConj(v[j + i * N])) > rtol) goto restore;
330: }
331: }
332: *fl = PETSC_TRUE;
333: restore:
334: PetscCall(MatDenseRestoreArrayRead(A, &v));
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: static PetscErrorCode MatIsSymmetric_SeqDense(Mat A, PetscReal rtol, PetscBool *fl)
339: {
340: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
341: PetscInt i, j, m = A->rmap->n, N = a->lda;
342: const PetscScalar *v;
344: PetscFunctionBegin;
345: *fl = PETSC_FALSE;
346: if (A->rmap->n != A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
347: PetscCall(MatDenseGetArrayRead(A, &v));
348: for (i = 0; i < m; i++) {
349: for (j = i; j < m; j++) {
350: if (PetscAbsScalar(v[i + j * N] - v[j + i * N]) > rtol) goto restore;
351: }
352: }
353: *fl = PETSC_TRUE;
354: restore:
355: PetscCall(MatDenseRestoreArrayRead(A, &v));
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi, Mat A, MatDuplicateOption cpvalues)
360: {
361: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
362: PetscInt lda = mat->lda, j, m, nlda = lda;
363: PetscBool isdensecpu;
365: PetscFunctionBegin;
366: PetscCall(PetscLayoutReference(A->rmap, &newi->rmap));
367: PetscCall(PetscLayoutReference(A->cmap, &newi->cmap));
368: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */
369: PetscCall(MatDenseSetLDA(newi, lda));
370: }
371: PetscCall(PetscObjectTypeCompare((PetscObject)newi, MATSEQDENSE, &isdensecpu));
372: if (isdensecpu) PetscCall(MatSeqDenseSetPreallocation(newi, NULL));
373: if (cpvalues == MAT_COPY_VALUES) {
374: const PetscScalar *av;
375: PetscScalar *v;
377: PetscCall(MatDenseGetArrayRead(A, &av));
378: PetscCall(MatDenseGetArrayWrite(newi, &v));
379: PetscCall(MatDenseGetLDA(newi, &nlda));
380: m = A->rmap->n;
381: if (lda > m || nlda > m) {
382: for (j = 0; j < A->cmap->n; j++) PetscCall(PetscArraycpy(PetscSafePointerPlusOffset(v, j * nlda), PetscSafePointerPlusOffset(av, j * lda), m));
383: } else {
384: PetscCall(PetscArraycpy(v, av, A->rmap->n * A->cmap->n));
385: }
386: PetscCall(MatDenseRestoreArrayWrite(newi, &v));
387: PetscCall(MatDenseRestoreArrayRead(A, &av));
388: PetscCall(MatPropagateSymmetryOptions(A, newi));
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 %" PetscBLASInt_FMT, 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 %" PetscBLASInt_FMT, 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 %" PetscBLASInt_FMT, 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 %" PetscBLASInt_FMT, 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 %" PetscBLASInt_FMT, 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 %" PetscBLASInt_FMT, 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 %" PetscBLASInt_FMT, 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 %" PetscBLASInt_FMT, 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, PETSC_UNUSED 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 %" PetscBLASInt_FMT, info);
807: PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Bad LU factorization %" PetscBLASInt_FMT, 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)
823: {
824: PetscFunctionBegin;
825: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
826: PetscUseTypeMethod(fact, lufactor, NULL, NULL, info);
827: PetscFunctionReturn(PETSC_SUCCESS);
828: }
830: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, IS col, PETSC_UNUSED const MatFactorInfo *info)
831: {
832: PetscFunctionBegin;
833: fact->preallocated = PETSC_TRUE;
834: fact->assembled = PETSC_TRUE;
835: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
836: PetscFunctionReturn(PETSC_SUCCESS);
837: }
839: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
840: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A, IS perm, PETSC_UNUSED const MatFactorInfo *minfo)
841: {
842: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
843: PetscBLASInt info, n;
845: PetscFunctionBegin;
846: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
847: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
848: if (A->spd == PETSC_BOOL3_TRUE) {
849: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
850: PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("L", &n, mat->v, &mat->lda, &info));
851: PetscCall(PetscFPTrapPop());
852: #if defined(PETSC_USE_COMPLEX)
853: } else if (A->hermitian == PETSC_BOOL3_TRUE) {
854: if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); }
855: if (!mat->fwork) {
856: PetscScalar dummy;
858: mat->lfwork = -1;
859: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
860: PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info));
861: PetscCall(PetscFPTrapPop());
862: PetscCall(PetscBLASIntCast((PetscCount)(PetscRealPart(dummy)), &mat->lfwork));
863: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
864: }
865: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
866: PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
867: PetscCall(PetscFPTrapPop());
868: #endif
869: } else { /* symmetric case */
870: if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); }
871: if (!mat->fwork) {
872: PetscScalar dummy;
874: mat->lfwork = -1;
875: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
876: PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info));
877: PetscCall(PetscFPTrapPop());
878: PetscCall(PetscBLASIntCast((PetscCount)(PetscRealPart(dummy)), &mat->lfwork));
879: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
880: }
881: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
882: PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
883: PetscCall(PetscFPTrapPop());
884: }
885: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad factorization: zero pivot in row %" PetscBLASInt_FMT, info - 1);
887: A->ops->solve = MatSolve_SeqDense_Cholesky;
888: A->ops->matsolve = MatMatSolve_SeqDense_Cholesky;
889: A->ops->solvetranspose = MatSolveTranspose_SeqDense_Cholesky;
890: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky;
891: A->factortype = MAT_FACTOR_CHOLESKY;
893: PetscCall(PetscFree(A->solvertype));
894: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
896: PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info)
901: {
902: PetscFunctionBegin;
903: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
904: PetscUseTypeMethod(fact, choleskyfactor, NULL, info);
905: PetscFunctionReturn(PETSC_SUCCESS);
906: }
908: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info)
909: {
910: PetscFunctionBegin;
911: fact->assembled = PETSC_TRUE;
912: fact->preallocated = PETSC_TRUE;
913: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
914: PetscFunctionReturn(PETSC_SUCCESS);
915: }
917: PetscErrorCode MatQRFactor_SeqDense(Mat A, IS col, PETSC_UNUSED const MatFactorInfo *minfo)
918: {
919: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
920: PetscBLASInt n, m, info, min, max;
922: PetscFunctionBegin;
923: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
924: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
925: max = PetscMax(m, n);
926: min = PetscMin(m, n);
927: if (!mat->tau) { PetscCall(PetscMalloc1(min, &mat->tau)); }
928: if (!mat->pivots) { PetscCall(PetscMalloc1(n, &mat->pivots)); }
929: if (!mat->qrrhs) PetscCall(MatCreateVecs(A, NULL, &mat->qrrhs));
930: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
931: if (!mat->fwork) {
932: PetscScalar dummy;
934: mat->lfwork = -1;
935: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
936: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, &dummy, &mat->lfwork, &info));
937: PetscCall(PetscFPTrapPop());
938: PetscCall(PetscBLASIntCast((PetscCount)(PetscRealPart(dummy)), &mat->lfwork));
939: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
940: }
941: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
942: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, mat->fwork, &mat->lfwork, &info));
943: PetscCall(PetscFPTrapPop());
944: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to QR factorization %" PetscBLASInt_FMT, info);
945: // 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
946: mat->rank = min;
948: A->ops->solve = MatSolve_SeqDense_QR;
949: A->ops->matsolve = MatMatSolve_SeqDense_QR;
950: A->factortype = MAT_FACTOR_QR;
951: if (m == n) {
952: A->ops->solvetranspose = MatSolveTranspose_SeqDense_QR;
953: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
954: }
956: PetscCall(PetscFree(A->solvertype));
957: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
959: PetscCall(PetscLogFlops(2.0 * min * min * (max - min / 3.0)));
960: PetscFunctionReturn(PETSC_SUCCESS);
961: }
963: static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info)
964: {
965: PetscFunctionBegin;
966: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
967: PetscUseMethod(fact, "MatQRFactor_C", (Mat, IS, const MatFactorInfo *), (fact, NULL, info));
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info)
972: {
973: PetscFunctionBegin;
974: fact->assembled = PETSC_TRUE;
975: fact->preallocated = PETSC_TRUE;
976: PetscCall(PetscObjectComposeFunction((PetscObject)fact, "MatQRFactorNumeric_C", MatQRFactorNumeric_SeqDense));
977: PetscFunctionReturn(PETSC_SUCCESS);
978: }
980: /* uses LAPACK */
981: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A, MatFactorType ftype, Mat *fact)
982: {
983: PetscFunctionBegin;
984: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), fact));
985: PetscCall(MatSetSizes(*fact, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
986: PetscCall(MatSetType(*fact, MATDENSE));
987: (*fact)->trivialsymbolic = PETSC_TRUE;
988: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
989: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
990: (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
991: } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
992: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
993: } else if (ftype == MAT_FACTOR_QR) {
994: PetscCall(PetscObjectComposeFunction((PetscObject)*fact, "MatQRFactorSymbolic_C", MatQRFactorSymbolic_SeqDense));
995: }
996: (*fact)->factortype = ftype;
998: PetscCall(PetscFree((*fact)->solvertype));
999: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*fact)->solvertype));
1000: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_LU]));
1001: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ILU]));
1002: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_CHOLESKY]));
1003: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ICC]));
1004: PetscFunctionReturn(PETSC_SUCCESS);
1005: }
1007: static PetscErrorCode MatSOR_SeqDense(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal shift, PetscInt its, PetscInt lits, Vec xx)
1008: {
1009: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1010: PetscScalar *x, *v = mat->v, zero = 0.0, xt;
1011: const PetscScalar *b;
1012: PetscInt m = A->rmap->n, i;
1013: PetscBLASInt o = 1, bm = 0;
1015: PetscFunctionBegin;
1016: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1017: PetscCheck(A->offloadmask != PETSC_OFFLOAD_GPU, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented");
1018: #endif
1019: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
1020: PetscCall(PetscBLASIntCast(m, &bm));
1021: if (flag & SOR_ZERO_INITIAL_GUESS) {
1022: /* this is a hack fix, should have another version without the second BLASdotu */
1023: PetscCall(VecSet(xx, zero));
1024: }
1025: PetscCall(VecGetArray(xx, &x));
1026: PetscCall(VecGetArrayRead(bb, &b));
1027: its = its * lits;
1028: 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);
1029: while (its--) {
1030: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1031: for (i = 0; i < m; i++) {
1032: PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o));
1033: x[i] = (1. - omega) * x[i] + (xt + v[i + i * m] * x[i]) * omega / (v[i + i * m] + shift);
1034: }
1035: }
1036: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1037: for (i = m - 1; i >= 0; i--) {
1038: PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o));
1039: x[i] = (1. - omega) * x[i] + (xt + v[i + i * m] * x[i]) * omega / (v[i + i * m] + shift);
1040: }
1041: }
1042: }
1043: PetscCall(VecRestoreArrayRead(bb, &b));
1044: PetscCall(VecRestoreArray(xx, &x));
1045: PetscFunctionReturn(PETSC_SUCCESS);
1046: }
1048: static PetscErrorCode MatMultColumnRangeKernel_SeqDense(Mat A, Vec xx, Vec yy, PetscInt c_start, PetscInt c_end, PetscBool trans, PetscBool herm)
1049: {
1050: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1051: PetscScalar *y, _DOne = 1.0, _DZero = 0.0;
1052: PetscBLASInt m, n, _One = 1;
1053: const PetscScalar *v = mat->v, *x;
1055: PetscFunctionBegin;
1056: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1057: PetscCall(PetscBLASIntCast(c_end - c_start, &n));
1058: PetscCall(VecGetArrayRead(xx, &x));
1059: PetscCall(VecGetArrayWrite(yy, &y));
1060: if (!m || !n) {
1061: PetscBLASInt i;
1062: if (trans)
1063: for (i = 0; i < n; i++) y[i] = 0.0;
1064: else
1065: for (i = 0; i < m; i++) y[i] = 0.0;
1066: } else {
1067: if (trans) {
1068: if (herm) PetscCallBLAS("BLASgemv", BLASgemv_("C", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DZero, y + c_start, &_One));
1069: else PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DZero, y + c_start, &_One));
1070: } else {
1071: PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x + c_start, &_One, &_DZero, y, &_One));
1072: }
1073: PetscCall(PetscLogFlops(2.0 * m * n - n));
1074: }
1075: PetscCall(VecRestoreArrayRead(xx, &x));
1076: PetscCall(VecRestoreArrayWrite(yy, &y));
1077: PetscFunctionReturn(PETSC_SUCCESS);
1078: }
1080: PetscErrorCode MatMultHermitianTransposeColumnRange_SeqDense(Mat A, Vec xx, Vec yy, PetscInt c_start, PetscInt c_end)
1081: {
1082: PetscFunctionBegin;
1083: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, c_start, c_end, PETSC_TRUE, PETSC_TRUE));
1084: PetscFunctionReturn(PETSC_SUCCESS);
1085: }
1087: PetscErrorCode MatMult_SeqDense(Mat A, Vec xx, Vec yy)
1088: {
1089: PetscFunctionBegin;
1090: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_FALSE, PETSC_FALSE));
1091: PetscFunctionReturn(PETSC_SUCCESS);
1092: }
1094: PetscErrorCode MatMultTranspose_SeqDense(Mat A, Vec xx, Vec yy)
1095: {
1096: PetscFunctionBegin;
1097: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_FALSE));
1098: PetscFunctionReturn(PETSC_SUCCESS);
1099: }
1101: PetscErrorCode MatMultHermitianTranspose_SeqDense(Mat A, Vec xx, Vec yy)
1102: {
1103: PetscFunctionBegin;
1104: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_TRUE));
1105: PetscFunctionReturn(PETSC_SUCCESS);
1106: }
1108: static PetscErrorCode MatMultAddColumnRangeKernel_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end, PetscBool trans, PetscBool herm)
1109: {
1110: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1111: const PetscScalar *v = mat->v, *x;
1112: PetscScalar *y, _DOne = 1.0;
1113: PetscBLASInt m, n, _One = 1;
1115: PetscFunctionBegin;
1116: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1117: PetscCall(PetscBLASIntCast(c_end - c_start, &n));
1118: PetscCall(VecCopy(zz, yy));
1119: if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS);
1120: PetscCall(VecGetArray(yy, &y));
1121: PetscCall(VecGetArrayRead(xx, &x));
1122: if (trans) {
1123: if (herm) PetscCallBLAS("BLASgemv", BLASgemv_("C", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DOne, y + c_start, &_One));
1124: else PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DOne, y + c_start, &_One));
1125: } else {
1126: PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x + c_start, &_One, &_DOne, y, &_One));
1127: }
1128: PetscCall(VecRestoreArrayRead(xx, &x));
1129: PetscCall(VecRestoreArray(yy, &y));
1130: PetscCall(PetscLogFlops(2.0 * m * n));
1131: PetscFunctionReturn(PETSC_SUCCESS);
1132: }
1134: PetscErrorCode MatMultAddColumnRange_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end)
1135: {
1136: PetscFunctionBegin;
1137: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, c_start, c_end, PETSC_FALSE, PETSC_FALSE));
1138: PetscFunctionReturn(PETSC_SUCCESS);
1139: }
1141: PetscErrorCode MatMultHermitianTransposeAddColumnRange_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end)
1142: {
1143: PetscFunctionBegin;
1144: PetscMPIInt rank;
1145: PetscCallMPI(MPI_Comm_rank(MPI_COMM_WORLD, &rank));
1146: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, c_start, c_end, PETSC_TRUE, PETSC_TRUE));
1147: PetscFunctionReturn(PETSC_SUCCESS);
1148: }
1150: PetscErrorCode MatMultAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1151: {
1152: PetscFunctionBegin;
1153: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_FALSE, PETSC_FALSE));
1154: PetscFunctionReturn(PETSC_SUCCESS);
1155: }
1157: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1158: {
1159: PetscFunctionBegin;
1160: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_FALSE));
1161: PetscFunctionReturn(PETSC_SUCCESS);
1162: }
1164: PetscErrorCode MatMultHermitianTransposeAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1165: {
1166: PetscFunctionBegin;
1167: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_TRUE));
1168: PetscFunctionReturn(PETSC_SUCCESS);
1169: }
1171: static PetscErrorCode MatGetRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1172: {
1173: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1174: PetscInt i;
1176: PetscFunctionBegin;
1177: if (ncols) *ncols = A->cmap->n;
1178: if (cols) {
1179: PetscCall(PetscMalloc1(A->cmap->n, cols));
1180: for (i = 0; i < A->cmap->n; i++) (*cols)[i] = i;
1181: }
1182: if (vals) {
1183: const PetscScalar *v;
1185: PetscCall(MatDenseGetArrayRead(A, &v));
1186: PetscCall(PetscMalloc1(A->cmap->n, vals));
1187: v += row;
1188: for (i = 0; i < A->cmap->n; i++) {
1189: (*vals)[i] = *v;
1190: v += mat->lda;
1191: }
1192: PetscCall(MatDenseRestoreArrayRead(A, &v));
1193: }
1194: PetscFunctionReturn(PETSC_SUCCESS);
1195: }
1197: static PetscErrorCode MatRestoreRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1198: {
1199: PetscFunctionBegin;
1200: if (cols) PetscCall(PetscFree(*cols));
1201: if (vals) PetscCall(PetscFree(*vals));
1202: PetscFunctionReturn(PETSC_SUCCESS);
1203: }
1205: static PetscErrorCode MatSetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], const PetscScalar v[], InsertMode addv)
1206: {
1207: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1208: PetscScalar *av;
1209: PetscInt i, j, idx = 0;
1210: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1211: PetscOffloadMask oldf;
1212: #endif
1214: PetscFunctionBegin;
1215: PetscCall(MatDenseGetArray(A, &av));
1216: if (!mat->roworiented) {
1217: if (addv == INSERT_VALUES) {
1218: for (j = 0; j < n; j++) {
1219: if (indexn[j] < 0) {
1220: idx += m;
1221: continue;
1222: }
1223: 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);
1224: for (i = 0; i < m; i++) {
1225: if (indexm[i] < 0) {
1226: idx++;
1227: continue;
1228: }
1229: 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);
1230: av[indexn[j] * mat->lda + indexm[i]] = v ? v[idx++] : (idx++, 0.0);
1231: }
1232: }
1233: } else {
1234: for (j = 0; j < n; j++) {
1235: if (indexn[j] < 0) {
1236: idx += m;
1237: continue;
1238: }
1239: 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);
1240: for (i = 0; i < m; i++) {
1241: if (indexm[i] < 0) {
1242: idx++;
1243: continue;
1244: }
1245: 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);
1246: av[indexn[j] * mat->lda + indexm[i]] += v ? v[idx++] : (idx++, 0.0);
1247: }
1248: }
1249: }
1250: } else {
1251: if (addv == INSERT_VALUES) {
1252: for (i = 0; i < m; i++) {
1253: if (indexm[i] < 0) {
1254: idx += n;
1255: continue;
1256: }
1257: 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);
1258: for (j = 0; j < n; j++) {
1259: if (indexn[j] < 0) {
1260: idx++;
1261: continue;
1262: }
1263: 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);
1264: av[indexn[j] * mat->lda + indexm[i]] = v ? v[idx++] : (idx++, 0.0);
1265: }
1266: }
1267: } else {
1268: for (i = 0; i < m; i++) {
1269: if (indexm[i] < 0) {
1270: idx += n;
1271: continue;
1272: }
1273: 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);
1274: for (j = 0; j < n; j++) {
1275: if (indexn[j] < 0) {
1276: idx++;
1277: continue;
1278: }
1279: 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);
1280: av[indexn[j] * mat->lda + indexm[i]] += v ? v[idx++] : (idx++, 0.0);
1281: }
1282: }
1283: }
1284: }
1285: /* hack to prevent unneeded copy to the GPU while returning the array */
1286: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1287: oldf = A->offloadmask;
1288: A->offloadmask = PETSC_OFFLOAD_GPU;
1289: #endif
1290: PetscCall(MatDenseRestoreArray(A, &av));
1291: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1292: A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1293: #endif
1294: PetscFunctionReturn(PETSC_SUCCESS);
1295: }
1297: static PetscErrorCode MatGetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], PetscScalar v[])
1298: {
1299: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1300: const PetscScalar *vv;
1301: PetscInt i, j;
1303: PetscFunctionBegin;
1304: PetscCall(MatDenseGetArrayRead(A, &vv));
1305: /* row-oriented output */
1306: for (i = 0; i < m; i++) {
1307: if (indexm[i] < 0) {
1308: v += n;
1309: continue;
1310: }
1311: 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);
1312: for (j = 0; j < n; j++) {
1313: if (indexn[j] < 0) {
1314: v++;
1315: continue;
1316: }
1317: 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);
1318: *v++ = vv[indexn[j] * mat->lda + indexm[i]];
1319: }
1320: }
1321: PetscCall(MatDenseRestoreArrayRead(A, &vv));
1322: PetscFunctionReturn(PETSC_SUCCESS);
1323: }
1325: PetscErrorCode MatView_Dense_Binary(Mat mat, PetscViewer viewer)
1326: {
1327: PetscBool skipHeader;
1328: PetscViewerFormat format;
1329: PetscInt header[4], M, N, m, lda, i, j;
1330: PetscCount k;
1331: const PetscScalar *v;
1332: PetscScalar *vwork;
1334: PetscFunctionBegin;
1335: PetscCall(PetscViewerSetUp(viewer));
1336: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
1337: PetscCall(PetscViewerGetFormat(viewer, &format));
1338: if (skipHeader) format = PETSC_VIEWER_NATIVE;
1340: PetscCall(MatGetSize(mat, &M, &N));
1342: /* write matrix header */
1343: header[0] = MAT_FILE_CLASSID;
1344: header[1] = M;
1345: header[2] = N;
1346: header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M * N;
1347: if (!skipHeader) PetscCall(PetscViewerBinaryWrite(viewer, header, 4, PETSC_INT));
1349: PetscCall(MatGetLocalSize(mat, &m, NULL));
1350: if (format != PETSC_VIEWER_NATIVE) {
1351: PetscInt nnz = m * N, *iwork;
1352: /* store row lengths for each row */
1353: PetscCall(PetscMalloc1(nnz, &iwork));
1354: for (i = 0; i < m; i++) iwork[i] = N;
1355: PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1356: /* store column indices (zero start index) */
1357: for (k = 0, i = 0; i < m; i++)
1358: for (j = 0; j < N; j++, k++) iwork[k] = j;
1359: PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1360: PetscCall(PetscFree(iwork));
1361: }
1362: /* store matrix values as a dense matrix in row major order */
1363: PetscCall(PetscMalloc1(m * N, &vwork));
1364: PetscCall(MatDenseGetArrayRead(mat, &v));
1365: PetscCall(MatDenseGetLDA(mat, &lda));
1366: for (k = 0, i = 0; i < m; i++)
1367: for (j = 0; j < N; j++, k++) vwork[k] = v[i + (size_t)lda * j];
1368: PetscCall(MatDenseRestoreArrayRead(mat, &v));
1369: PetscCall(PetscViewerBinaryWriteAll(viewer, vwork, m * N, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1370: PetscCall(PetscFree(vwork));
1371: PetscFunctionReturn(PETSC_SUCCESS);
1372: }
1374: PetscErrorCode MatLoad_Dense_Binary(Mat mat, PetscViewer viewer)
1375: {
1376: PetscBool skipHeader;
1377: PetscInt header[4], M, N, m, nz, lda, i, j, k;
1378: PetscInt rows, cols;
1379: PetscScalar *v, *vwork;
1381: PetscFunctionBegin;
1382: PetscCall(PetscViewerSetUp(viewer));
1383: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
1385: if (!skipHeader) {
1386: PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT));
1387: PetscCheck(header[0] == MAT_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file");
1388: M = header[1];
1389: N = header[2];
1390: PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M);
1391: PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N);
1392: nz = header[3];
1393: PetscCheck(nz == MATRIX_BINARY_FORMAT_DENSE || nz >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Unknown matrix format %" PetscInt_FMT " in file", nz);
1394: } else {
1395: PetscCall(MatGetSize(mat, &M, &N));
1396: 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");
1397: nz = MATRIX_BINARY_FORMAT_DENSE;
1398: }
1400: /* setup global sizes if not set */
1401: if (mat->rmap->N < 0) mat->rmap->N = M;
1402: if (mat->cmap->N < 0) mat->cmap->N = N;
1403: PetscCall(MatSetUp(mat));
1404: /* check if global sizes are correct */
1405: PetscCall(MatGetSize(mat, &rows, &cols));
1406: 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);
1408: PetscCall(MatGetSize(mat, NULL, &N));
1409: PetscCall(MatGetLocalSize(mat, &m, NULL));
1410: PetscCall(MatDenseGetArray(mat, &v));
1411: PetscCall(MatDenseGetLDA(mat, &lda));
1412: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */
1413: PetscCount nnz = (size_t)m * N;
1414: /* read in matrix values */
1415: PetscCall(PetscMalloc1(nnz, &vwork));
1416: PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1417: /* store values in column major order */
1418: for (j = 0; j < N; j++)
1419: for (i = 0; i < m; i++) v[i + (size_t)lda * j] = vwork[(size_t)i * N + j];
1420: PetscCall(PetscFree(vwork));
1421: } else { /* matrix in file is sparse format */
1422: PetscInt nnz = 0, *rlens, *icols;
1423: /* read in row lengths */
1424: PetscCall(PetscMalloc1(m, &rlens));
1425: PetscCall(PetscViewerBinaryReadAll(viewer, rlens, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1426: for (i = 0; i < m; i++) nnz += rlens[i];
1427: /* read in column indices and values */
1428: PetscCall(PetscMalloc2(nnz, &icols, nnz, &vwork));
1429: PetscCall(PetscViewerBinaryReadAll(viewer, icols, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1430: PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1431: /* store values in column major order */
1432: for (k = 0, i = 0; i < m; i++)
1433: for (j = 0; j < rlens[i]; j++, k++) v[i + lda * icols[k]] = vwork[k];
1434: PetscCall(PetscFree(rlens));
1435: PetscCall(PetscFree2(icols, vwork));
1436: }
1437: PetscCall(MatDenseRestoreArray(mat, &v));
1438: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
1439: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
1440: PetscFunctionReturn(PETSC_SUCCESS);
1441: }
1443: static PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1444: {
1445: PetscBool isbinary, ishdf5;
1447: PetscFunctionBegin;
1450: /* force binary viewer to load .info file if it has not yet done so */
1451: PetscCall(PetscViewerSetUp(viewer));
1452: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1453: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1454: if (isbinary) {
1455: PetscCall(MatLoad_Dense_Binary(newMat, viewer));
1456: } else if (ishdf5) {
1457: #if defined(PETSC_HAVE_HDF5)
1458: PetscCall(MatLoad_Dense_HDF5(newMat, viewer));
1459: #else
1460: SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1461: #endif
1462: } else {
1463: 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);
1464: }
1465: PetscFunctionReturn(PETSC_SUCCESS);
1466: }
1468: static PetscErrorCode MatView_SeqDense_ASCII(Mat A, PetscViewer viewer)
1469: {
1470: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1471: PetscInt i, j;
1472: const char *name;
1473: PetscScalar *v, *av;
1474: PetscViewerFormat format;
1475: #if defined(PETSC_USE_COMPLEX)
1476: PetscBool allreal = PETSC_TRUE;
1477: #endif
1479: PetscFunctionBegin;
1480: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&av));
1481: PetscCall(PetscViewerGetFormat(viewer, &format));
1482: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1483: PetscFunctionReturn(PETSC_SUCCESS); /* do nothing for now */
1484: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1485: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1486: for (i = 0; i < A->rmap->n; i++) {
1487: v = av + i;
1488: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i));
1489: for (j = 0; j < A->cmap->n; j++) {
1490: #if defined(PETSC_USE_COMPLEX)
1491: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1492: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", j, (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1493: } else if (PetscRealPart(*v)) {
1494: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)PetscRealPart(*v)));
1495: }
1496: #else
1497: if (*v) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)*v));
1498: #endif
1499: v += a->lda;
1500: }
1501: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1502: }
1503: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1504: } else {
1505: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1506: #if defined(PETSC_USE_COMPLEX)
1507: /* determine if matrix has all real values */
1508: for (j = 0; j < A->cmap->n; j++) {
1509: v = av + j * a->lda;
1510: for (i = 0; i < A->rmap->n; i++) {
1511: if (PetscImaginaryPart(v[i])) {
1512: allreal = PETSC_FALSE;
1513: break;
1514: }
1515: }
1516: }
1517: #endif
1518: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1519: PetscCall(PetscObjectGetName((PetscObject)A, &name));
1520: PetscCall(PetscViewerASCIIPrintf(viewer, "%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n", A->rmap->n, A->cmap->n));
1521: PetscCall(PetscViewerASCIIPrintf(viewer, "%s = zeros(%" PetscInt_FMT ",%" PetscInt_FMT ");\n", name, A->rmap->n, A->cmap->n));
1522: PetscCall(PetscViewerASCIIPrintf(viewer, "%s = [\n", name));
1523: }
1525: for (i = 0; i < A->rmap->n; i++) {
1526: v = av + i;
1527: for (j = 0; j < A->cmap->n; j++) {
1528: #if defined(PETSC_USE_COMPLEX)
1529: if (allreal) {
1530: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)PetscRealPart(*v)));
1531: } else {
1532: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei ", (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1533: }
1534: #else
1535: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)*v));
1536: #endif
1537: v += a->lda;
1538: }
1539: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1540: }
1541: if (format == PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer, "];\n"));
1542: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1543: }
1544: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&av));
1545: PetscCall(PetscViewerFlush(viewer));
1546: PetscFunctionReturn(PETSC_SUCCESS);
1547: }
1549: #include <petscdraw.h>
1550: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw, void *Aa)
1551: {
1552: Mat A = (Mat)Aa;
1553: PetscInt m = A->rmap->n, n = A->cmap->n, i, j;
1554: int color = PETSC_DRAW_WHITE;
1555: const PetscScalar *v;
1556: PetscViewer viewer;
1557: PetscReal xl, yl, xr, yr, x_l, x_r, y_l, y_r;
1558: PetscViewerFormat format;
1560: PetscFunctionBegin;
1561: PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer));
1562: PetscCall(PetscViewerGetFormat(viewer, &format));
1563: PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));
1565: /* Loop over matrix elements drawing boxes */
1566: PetscCall(MatDenseGetArrayRead(A, &v));
1567: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1568: PetscDrawCollectiveBegin(draw);
1569: /* Blue for negative and Red for positive */
1570: for (j = 0; j < n; j++) {
1571: x_l = j;
1572: x_r = x_l + 1.0;
1573: for (i = 0; i < m; i++) {
1574: y_l = m - i - 1.0;
1575: y_r = y_l + 1.0;
1576: if (PetscRealPart(v[j * m + i]) > 0.) color = PETSC_DRAW_RED;
1577: else if (PetscRealPart(v[j * m + i]) < 0.) color = PETSC_DRAW_BLUE;
1578: else continue;
1579: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1580: }
1581: }
1582: PetscDrawCollectiveEnd(draw);
1583: } else {
1584: /* use contour shading to indicate magnitude of values */
1585: /* first determine max of all nonzero values */
1586: PetscReal minv = 0.0, maxv = 0.0;
1587: PetscDraw popup;
1589: for (i = 0; i < m * n; i++) {
1590: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1591: }
1592: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1593: PetscCall(PetscDrawGetPopup(draw, &popup));
1594: PetscCall(PetscDrawScalePopup(popup, minv, maxv));
1596: PetscDrawCollectiveBegin(draw);
1597: for (j = 0; j < n; j++) {
1598: x_l = j;
1599: x_r = x_l + 1.0;
1600: for (i = 0; i < m; i++) {
1601: y_l = m - i - 1.0;
1602: y_r = y_l + 1.0;
1603: color = PetscDrawRealToColor(PetscAbsScalar(v[j * m + i]), minv, maxv);
1604: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1605: }
1606: }
1607: PetscDrawCollectiveEnd(draw);
1608: }
1609: PetscCall(MatDenseRestoreArrayRead(A, &v));
1610: PetscFunctionReturn(PETSC_SUCCESS);
1611: }
1613: static PetscErrorCode MatView_SeqDense_Draw(Mat A, PetscViewer viewer)
1614: {
1615: PetscDraw draw;
1616: PetscBool isnull;
1617: PetscReal xr, yr, xl, yl, h, w;
1619: PetscFunctionBegin;
1620: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1621: PetscCall(PetscDrawIsNull(draw, &isnull));
1622: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
1624: xr = A->cmap->n;
1625: yr = A->rmap->n;
1626: h = yr / 10.0;
1627: w = xr / 10.0;
1628: xr += w;
1629: yr += h;
1630: xl = -w;
1631: yl = -h;
1632: PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr));
1633: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer));
1634: PetscCall(PetscDrawZoom(draw, MatView_SeqDense_Draw_Zoom, A));
1635: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL));
1636: PetscCall(PetscDrawSave(draw));
1637: PetscFunctionReturn(PETSC_SUCCESS);
1638: }
1640: PetscErrorCode MatView_SeqDense(Mat A, PetscViewer viewer)
1641: {
1642: PetscBool iascii, isbinary, isdraw;
1644: PetscFunctionBegin;
1645: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1646: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1647: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1648: if (iascii) PetscCall(MatView_SeqDense_ASCII(A, viewer));
1649: else if (isbinary) PetscCall(MatView_Dense_Binary(A, viewer));
1650: else if (isdraw) PetscCall(MatView_SeqDense_Draw(A, viewer));
1651: PetscFunctionReturn(PETSC_SUCCESS);
1652: }
1654: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A, const PetscScalar *array)
1655: {
1656: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1658: PetscFunctionBegin;
1659: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1660: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1661: PetscCheck(!a->unplacedarray, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseResetArray() first");
1662: a->unplacedarray = a->v;
1663: a->unplaced_user_alloc = a->user_alloc;
1664: a->v = (PetscScalar *)array;
1665: a->user_alloc = PETSC_TRUE;
1666: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1667: A->offloadmask = PETSC_OFFLOAD_CPU;
1668: #endif
1669: PetscFunctionReturn(PETSC_SUCCESS);
1670: }
1672: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1673: {
1674: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1676: PetscFunctionBegin;
1677: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1678: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1679: a->v = a->unplacedarray;
1680: a->user_alloc = a->unplaced_user_alloc;
1681: a->unplacedarray = NULL;
1682: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1683: A->offloadmask = PETSC_OFFLOAD_CPU;
1684: #endif
1685: PetscFunctionReturn(PETSC_SUCCESS);
1686: }
1688: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A, const PetscScalar *array)
1689: {
1690: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1692: PetscFunctionBegin;
1693: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1694: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1695: if (!a->user_alloc) PetscCall(PetscFree(a->v));
1696: a->v = (PetscScalar *)array;
1697: a->user_alloc = PETSC_FALSE;
1698: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1699: A->offloadmask = PETSC_OFFLOAD_CPU;
1700: #endif
1701: PetscFunctionReturn(PETSC_SUCCESS);
1702: }
1704: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1705: {
1706: Mat_SeqDense *l = (Mat_SeqDense *)mat->data;
1708: PetscFunctionBegin;
1709: PetscCall(PetscLogObjectState((PetscObject)mat, "Rows %" PetscInt_FMT " Cols %" PetscInt_FMT, mat->rmap->n, mat->cmap->n));
1710: PetscCall(VecDestroy(&l->qrrhs));
1711: PetscCall(PetscFree(l->tau));
1712: PetscCall(PetscFree(l->pivots));
1713: PetscCall(PetscFree(l->fwork));
1714: if (!l->user_alloc) PetscCall(PetscFree(l->v));
1715: if (!l->unplaced_user_alloc) PetscCall(PetscFree(l->unplacedarray));
1716: PetscCheck(!l->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1717: PetscCheck(!l->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1718: PetscCall(VecDestroy(&l->cvec));
1719: PetscCall(MatDestroy(&l->cmat));
1720: PetscCall(PetscFree(mat->data));
1722: PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL));
1723: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactor_C", NULL));
1724: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorSymbolic_C", NULL));
1725: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorNumeric_C", NULL));
1726: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetLDA_C", NULL));
1727: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseSetLDA_C", NULL));
1728: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArray_C", NULL));
1729: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArray_C", NULL));
1730: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDensePlaceArray_C", NULL));
1731: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseResetArray_C", NULL));
1732: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseReplaceArray_C", NULL));
1733: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayRead_C", NULL));
1734: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayRead_C", NULL));
1735: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayWrite_C", NULL));
1736: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayWrite_C", NULL));
1737: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqaij_C", NULL));
1738: #if defined(PETSC_HAVE_ELEMENTAL)
1739: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_elemental_C", NULL));
1740: #endif
1741: #if defined(PETSC_HAVE_SCALAPACK)
1742: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_scalapack_C", NULL));
1743: #endif
1744: #if defined(PETSC_HAVE_CUDA)
1745: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensecuda_C", NULL));
1746: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", NULL));
1747: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdense_C", NULL));
1748: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensecuda_C", NULL));
1749: #endif
1750: #if defined(PETSC_HAVE_HIP)
1751: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensehip_C", NULL));
1752: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", NULL));
1753: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdense_C", NULL));
1754: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensehip_C", NULL));
1755: #endif
1756: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatSeqDenseSetPreallocation_C", NULL));
1757: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqaij_seqdense_C", NULL));
1758: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdense_C", NULL));
1759: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqbaij_seqdense_C", NULL));
1760: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqsbaij_seqdense_C", NULL));
1762: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumn_C", NULL));
1763: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumn_C", NULL));
1764: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVec_C", NULL));
1765: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVec_C", NULL));
1766: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecRead_C", NULL));
1767: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecRead_C", NULL));
1768: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecWrite_C", NULL));
1769: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecWrite_C", NULL));
1770: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetSubMatrix_C", NULL));
1771: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreSubMatrix_C", NULL));
1772: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultAddColumnRange_C", NULL));
1773: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeColumnRange_C", NULL));
1774: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeAddColumnRange_C", NULL));
1775: PetscFunctionReturn(PETSC_SUCCESS);
1776: }
1778: static PetscErrorCode MatTranspose_SeqDense(Mat A, MatReuse reuse, Mat *matout)
1779: {
1780: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1781: PetscInt k, j, m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1782: PetscScalar *v, tmp;
1784: PetscFunctionBegin;
1785: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *matout));
1786: if (reuse == MAT_INPLACE_MATRIX) {
1787: if (m == n) { /* in place transpose */
1788: PetscCall(MatDenseGetArray(A, &v));
1789: for (j = 0; j < m; j++) {
1790: for (k = 0; k < j; k++) {
1791: tmp = v[j + k * M];
1792: v[j + k * M] = v[k + j * M];
1793: v[k + j * M] = tmp;
1794: }
1795: }
1796: PetscCall(MatDenseRestoreArray(A, &v));
1797: } else { /* reuse memory, temporary allocates new memory */
1798: PetscScalar *v2;
1799: PetscLayout tmplayout;
1801: PetscCall(PetscMalloc1((size_t)m * n, &v2));
1802: PetscCall(MatDenseGetArray(A, &v));
1803: for (j = 0; j < n; j++) {
1804: for (k = 0; k < m; k++) v2[j + (size_t)k * n] = v[k + (size_t)j * M];
1805: }
1806: PetscCall(PetscArraycpy(v, v2, (size_t)m * n));
1807: PetscCall(PetscFree(v2));
1808: PetscCall(MatDenseRestoreArray(A, &v));
1809: /* cleanup size dependent quantities */
1810: PetscCall(VecDestroy(&mat->cvec));
1811: PetscCall(MatDestroy(&mat->cmat));
1812: PetscCall(PetscFree(mat->pivots));
1813: PetscCall(PetscFree(mat->fwork));
1814: /* swap row/col layouts */
1815: PetscCall(PetscBLASIntCast(n, &mat->lda));
1816: tmplayout = A->rmap;
1817: A->rmap = A->cmap;
1818: A->cmap = tmplayout;
1819: }
1820: } else { /* out-of-place transpose */
1821: Mat tmat;
1822: Mat_SeqDense *tmatd;
1823: PetscScalar *v2;
1824: PetscInt M2;
1826: if (reuse == MAT_INITIAL_MATRIX) {
1827: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &tmat));
1828: PetscCall(MatSetSizes(tmat, A->cmap->n, A->rmap->n, A->cmap->n, A->rmap->n));
1829: PetscCall(MatSetType(tmat, ((PetscObject)A)->type_name));
1830: PetscCall(MatSeqDenseSetPreallocation(tmat, NULL));
1831: } else tmat = *matout;
1833: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&v));
1834: PetscCall(MatDenseGetArray(tmat, &v2));
1835: tmatd = (Mat_SeqDense *)tmat->data;
1836: M2 = tmatd->lda;
1837: for (j = 0; j < n; j++) {
1838: for (k = 0; k < m; k++) v2[j + k * M2] = v[k + j * M];
1839: }
1840: PetscCall(MatDenseRestoreArray(tmat, &v2));
1841: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&v));
1842: PetscCall(MatAssemblyBegin(tmat, MAT_FINAL_ASSEMBLY));
1843: PetscCall(MatAssemblyEnd(tmat, MAT_FINAL_ASSEMBLY));
1844: *matout = tmat;
1845: }
1846: PetscFunctionReturn(PETSC_SUCCESS);
1847: }
1849: static PetscErrorCode MatEqual_SeqDense(Mat A1, Mat A2, PetscBool *flg)
1850: {
1851: Mat_SeqDense *mat1 = (Mat_SeqDense *)A1->data;
1852: Mat_SeqDense *mat2 = (Mat_SeqDense *)A2->data;
1853: PetscInt i;
1854: const PetscScalar *v1, *v2;
1856: PetscFunctionBegin;
1857: if (A1->rmap->n != A2->rmap->n) {
1858: *flg = PETSC_FALSE;
1859: PetscFunctionReturn(PETSC_SUCCESS);
1860: }
1861: if (A1->cmap->n != A2->cmap->n) {
1862: *flg = PETSC_FALSE;
1863: PetscFunctionReturn(PETSC_SUCCESS);
1864: }
1865: PetscCall(MatDenseGetArrayRead(A1, &v1));
1866: PetscCall(MatDenseGetArrayRead(A2, &v2));
1867: for (i = 0; i < A1->cmap->n; i++) {
1868: PetscCall(PetscArraycmp(v1, v2, A1->rmap->n, flg));
1869: if (*flg == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS);
1870: v1 += mat1->lda;
1871: v2 += mat2->lda;
1872: }
1873: PetscCall(MatDenseRestoreArrayRead(A1, &v1));
1874: PetscCall(MatDenseRestoreArrayRead(A2, &v2));
1875: *flg = PETSC_TRUE;
1876: PetscFunctionReturn(PETSC_SUCCESS);
1877: }
1879: PetscErrorCode MatGetDiagonal_SeqDense(Mat A, Vec v)
1880: {
1881: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1882: PetscInt i, n, len;
1883: PetscScalar *x;
1884: const PetscScalar *vv;
1886: PetscFunctionBegin;
1887: PetscCall(VecGetSize(v, &n));
1888: PetscCall(VecGetArray(v, &x));
1889: len = PetscMin(A->rmap->n, A->cmap->n);
1890: PetscCall(MatDenseGetArrayRead(A, &vv));
1891: PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming mat and vec");
1892: for (i = 0; i < len; i++) x[i] = vv[i * mat->lda + i];
1893: PetscCall(MatDenseRestoreArrayRead(A, &vv));
1894: PetscCall(VecRestoreArray(v, &x));
1895: PetscFunctionReturn(PETSC_SUCCESS);
1896: }
1898: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A, Vec ll, Vec rr)
1899: {
1900: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1901: const PetscScalar *l, *r;
1902: PetscScalar x, *v, *vv;
1903: PetscInt i, j, m = A->rmap->n, n = A->cmap->n;
1905: PetscFunctionBegin;
1906: PetscCall(MatDenseGetArray(A, &vv));
1907: if (ll) {
1908: PetscCall(VecGetSize(ll, &m));
1909: PetscCall(VecGetArrayRead(ll, &l));
1910: PetscCheck(m == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vec wrong size");
1911: for (i = 0; i < m; i++) {
1912: x = l[i];
1913: v = vv + i;
1914: for (j = 0; j < n; j++) {
1915: (*v) *= x;
1916: v += mat->lda;
1917: }
1918: }
1919: PetscCall(VecRestoreArrayRead(ll, &l));
1920: PetscCall(PetscLogFlops(1.0 * n * m));
1921: }
1922: if (rr) {
1923: PetscCall(VecGetSize(rr, &n));
1924: PetscCall(VecGetArrayRead(rr, &r));
1925: PetscCheck(n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vec wrong size");
1926: for (i = 0; i < n; i++) {
1927: x = r[i];
1928: v = vv + i * mat->lda;
1929: for (j = 0; j < m; j++) (*v++) *= x;
1930: }
1931: PetscCall(VecRestoreArrayRead(rr, &r));
1932: PetscCall(PetscLogFlops(1.0 * n * m));
1933: }
1934: PetscCall(MatDenseRestoreArray(A, &vv));
1935: PetscFunctionReturn(PETSC_SUCCESS);
1936: }
1938: PetscErrorCode MatNorm_SeqDense(Mat A, NormType type, PetscReal *nrm)
1939: {
1940: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1941: PetscScalar *v, *vv;
1942: PetscReal sum = 0.0;
1943: PetscInt lda, m = A->rmap->n, i, j;
1945: PetscFunctionBegin;
1946: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&vv));
1947: PetscCall(MatDenseGetLDA(A, &lda));
1948: v = vv;
1949: if (type == NORM_FROBENIUS) {
1950: if (lda > m) {
1951: for (j = 0; j < A->cmap->n; j++) {
1952: v = vv + j * lda;
1953: for (i = 0; i < m; i++) {
1954: sum += PetscRealPart(PetscConj(*v) * (*v));
1955: v++;
1956: }
1957: }
1958: } else {
1959: #if defined(PETSC_USE_REAL___FP16)
1960: PetscBLASInt one = 1, cnt = A->cmap->n * A->rmap->n;
1961: PetscCallBLAS("BLASnrm2", *nrm = BLASnrm2_(&cnt, v, &one));
1962: }
1963: #else
1964: for (i = 0; i < A->cmap->n * A->rmap->n; i++) {
1965: sum += PetscRealPart(PetscConj(*v) * (*v));
1966: v++;
1967: }
1968: }
1969: *nrm = PetscSqrtReal(sum);
1970: #endif
1971: PetscCall(PetscLogFlops(2.0 * A->cmap->n * A->rmap->n));
1972: } else if (type == NORM_1) {
1973: *nrm = 0.0;
1974: for (j = 0; j < A->cmap->n; j++) {
1975: v = vv + j * mat->lda;
1976: sum = 0.0;
1977: for (i = 0; i < A->rmap->n; i++) {
1978: sum += PetscAbsScalar(*v);
1979: v++;
1980: }
1981: if (sum > *nrm) *nrm = sum;
1982: }
1983: PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
1984: } else if (type == NORM_INFINITY) {
1985: *nrm = 0.0;
1986: for (j = 0; j < A->rmap->n; j++) {
1987: v = vv + j;
1988: sum = 0.0;
1989: for (i = 0; i < A->cmap->n; i++) {
1990: sum += PetscAbsScalar(*v);
1991: v += mat->lda;
1992: }
1993: if (sum > *nrm) *nrm = sum;
1994: }
1995: PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
1996: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No two norm");
1997: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&vv));
1998: PetscFunctionReturn(PETSC_SUCCESS);
1999: }
2001: static PetscErrorCode MatSetOption_SeqDense(Mat A, MatOption op, PetscBool flg)
2002: {
2003: Mat_SeqDense *aij = (Mat_SeqDense *)A->data;
2005: PetscFunctionBegin;
2006: switch (op) {
2007: case MAT_ROW_ORIENTED:
2008: aij->roworiented = flg;
2009: break;
2010: default:
2011: break;
2012: }
2013: PetscFunctionReturn(PETSC_SUCCESS);
2014: }
2016: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
2017: {
2018: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
2019: PetscInt lda = l->lda, m = A->rmap->n, n = A->cmap->n, j;
2020: PetscScalar *v;
2022: PetscFunctionBegin;
2023: PetscCall(MatDenseGetArrayWrite(A, &v));
2024: if (lda > m) {
2025: for (j = 0; j < n; j++) PetscCall(PetscArrayzero(v + j * lda, m));
2026: } else {
2027: PetscCall(PetscArrayzero(v, PetscInt64Mult(m, n)));
2028: }
2029: PetscCall(MatDenseRestoreArrayWrite(A, &v));
2030: PetscFunctionReturn(PETSC_SUCCESS);
2031: }
2033: static PetscErrorCode MatZeroRows_SeqDense(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2034: {
2035: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
2036: PetscInt m = l->lda, n = A->cmap->n, i, j;
2037: PetscScalar *slot, *bb, *v;
2038: const PetscScalar *xx;
2040: PetscFunctionBegin;
2041: if (PetscDefined(USE_DEBUG)) {
2042: for (i = 0; i < N; i++) {
2043: PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
2044: 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);
2045: }
2046: }
2047: if (!N) PetscFunctionReturn(PETSC_SUCCESS);
2049: /* fix right-hand side if needed */
2050: if (x && b) {
2051: PetscCall(VecGetArrayRead(x, &xx));
2052: PetscCall(VecGetArray(b, &bb));
2053: for (i = 0; i < N; i++) bb[rows[i]] = diag * xx[rows[i]];
2054: PetscCall(VecRestoreArrayRead(x, &xx));
2055: PetscCall(VecRestoreArray(b, &bb));
2056: }
2058: PetscCall(MatDenseGetArray(A, &v));
2059: for (i = 0; i < N; i++) {
2060: slot = v + rows[i];
2061: for (j = 0; j < n; j++) {
2062: *slot = 0.0;
2063: slot += m;
2064: }
2065: }
2066: if (diag != 0.0) {
2067: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
2068: for (i = 0; i < N; i++) {
2069: slot = v + (m + 1) * rows[i];
2070: *slot = diag;
2071: }
2072: }
2073: PetscCall(MatDenseRestoreArray(A, &v));
2074: PetscFunctionReturn(PETSC_SUCCESS);
2075: }
2077: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A, PetscInt *lda)
2078: {
2079: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2081: PetscFunctionBegin;
2082: *lda = mat->lda;
2083: PetscFunctionReturn(PETSC_SUCCESS);
2084: }
2086: PetscErrorCode MatDenseGetArray_SeqDense(Mat A, PetscScalar **array)
2087: {
2088: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2090: PetscFunctionBegin;
2091: PetscCheck(!mat->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
2092: *array = mat->v;
2093: PetscFunctionReturn(PETSC_SUCCESS);
2094: }
2096: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A, PetscScalar **array)
2097: {
2098: PetscFunctionBegin;
2099: if (array) *array = NULL;
2100: PetscFunctionReturn(PETSC_SUCCESS);
2101: }
2103: /*@
2104: MatDenseGetLDA - gets the leading dimension of the array returned from `MatDenseGetArray()`
2106: Not Collective
2108: Input Parameter:
2109: . A - a `MATDENSE` or `MATDENSECUDA` matrix
2111: Output Parameter:
2112: . lda - the leading dimension
2114: Level: intermediate
2116: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseSetLDA()`
2117: @*/
2118: PetscErrorCode MatDenseGetLDA(Mat A, PetscInt *lda)
2119: {
2120: PetscFunctionBegin;
2122: PetscAssertPointer(lda, 2);
2123: MatCheckPreallocated(A, 1);
2124: PetscUseMethod(A, "MatDenseGetLDA_C", (Mat, PetscInt *), (A, lda));
2125: PetscFunctionReturn(PETSC_SUCCESS);
2126: }
2128: /*@
2129: MatDenseSetLDA - Sets the leading dimension of the array used by the `MATDENSE` matrix
2131: Collective if the matrix layouts have not yet been setup
2133: Input Parameters:
2134: + A - a `MATDENSE` or `MATDENSECUDA` matrix
2135: - lda - the leading dimension
2137: Level: intermediate
2139: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetLDA()`
2140: @*/
2141: PetscErrorCode MatDenseSetLDA(Mat A, PetscInt lda)
2142: {
2143: PetscFunctionBegin;
2145: PetscTryMethod(A, "MatDenseSetLDA_C", (Mat, PetscInt), (A, lda));
2146: PetscFunctionReturn(PETSC_SUCCESS);
2147: }
2149: /*@C
2150: MatDenseGetArray - gives read-write access to the array where the data for a `MATDENSE` matrix is stored
2152: Logically Collective
2154: Input Parameter:
2155: . A - a dense matrix
2157: Output Parameter:
2158: . array - pointer to the data
2160: Level: intermediate
2162: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2163: @*/
2164: PetscErrorCode MatDenseGetArray(Mat A, PetscScalar *array[]) PeNS
2165: {
2166: PetscFunctionBegin;
2168: PetscAssertPointer(array, 2);
2169: PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2170: PetscFunctionReturn(PETSC_SUCCESS);
2171: }
2173: /*@C
2174: MatDenseRestoreArray - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArray()`
2176: Logically Collective
2178: Input Parameters:
2179: + A - a dense matrix
2180: - array - pointer to the data (may be `NULL`)
2182: Level: intermediate
2184: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2185: @*/
2186: PetscErrorCode MatDenseRestoreArray(Mat A, PetscScalar *array[]) PeNS
2187: {
2188: PetscFunctionBegin;
2190: if (array) PetscAssertPointer(array, 2);
2191: PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2192: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2193: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2194: A->offloadmask = PETSC_OFFLOAD_CPU;
2195: #endif
2196: PetscFunctionReturn(PETSC_SUCCESS);
2197: }
2199: /*@C
2200: MatDenseGetArrayRead - gives read-only access to the array where the data for a `MATDENSE` matrix is stored
2202: Not Collective
2204: Input Parameter:
2205: . A - a dense matrix
2207: Output Parameter:
2208: . array - pointer to the data
2210: Level: intermediate
2212: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2213: @*/
2214: PetscErrorCode MatDenseGetArrayRead(Mat A, const PetscScalar *array[]) PeNS
2215: {
2216: PetscFunctionBegin;
2218: PetscAssertPointer(array, 2);
2219: PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2220: PetscFunctionReturn(PETSC_SUCCESS);
2221: }
2223: /*@C
2224: MatDenseRestoreArrayRead - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayRead()`
2226: Not Collective
2228: Input Parameters:
2229: + A - a dense matrix
2230: - array - pointer to the data (may be `NULL`)
2232: Level: intermediate
2234: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2235: @*/
2236: PetscErrorCode MatDenseRestoreArrayRead(Mat A, const PetscScalar *array[]) PeNS
2237: {
2238: PetscFunctionBegin;
2240: if (array) PetscAssertPointer(array, 2);
2241: PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2242: PetscFunctionReturn(PETSC_SUCCESS);
2243: }
2245: /*@C
2246: MatDenseGetArrayWrite - gives write-only access to the array where the data for a `MATDENSE` matrix is stored
2248: Not Collective
2250: Input Parameter:
2251: . A - a dense matrix
2253: Output Parameter:
2254: . array - pointer to the data
2256: Level: intermediate
2258: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2259: @*/
2260: PetscErrorCode MatDenseGetArrayWrite(Mat A, PetscScalar *array[]) PeNS
2261: {
2262: PetscFunctionBegin;
2264: PetscAssertPointer(array, 2);
2265: PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2266: PetscFunctionReturn(PETSC_SUCCESS);
2267: }
2269: /*@C
2270: MatDenseRestoreArrayWrite - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayWrite()`
2272: Not Collective
2274: Input Parameters:
2275: + A - a dense matrix
2276: - array - pointer to the data (may be `NULL`)
2278: Level: intermediate
2280: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2281: @*/
2282: PetscErrorCode MatDenseRestoreArrayWrite(Mat A, PetscScalar *array[]) PeNS
2283: {
2284: PetscFunctionBegin;
2286: if (array) PetscAssertPointer(array, 2);
2287: PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2288: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2289: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2290: A->offloadmask = PETSC_OFFLOAD_CPU;
2291: #endif
2292: PetscFunctionReturn(PETSC_SUCCESS);
2293: }
2295: /*@C
2296: MatDenseGetArrayAndMemType - gives read-write access to the array where the data for a `MATDENSE` matrix is stored
2298: Logically Collective
2300: Input Parameter:
2301: . A - a dense matrix
2303: Output Parameters:
2304: + array - pointer to the data
2305: - mtype - memory type of the returned pointer
2307: Level: intermediate
2309: Note:
2310: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2311: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2313: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArrayRead()`,
2314: `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2315: @*/
2316: PetscErrorCode MatDenseGetArrayAndMemType(Mat A, PetscScalar *array[], PetscMemType *mtype)
2317: {
2318: PetscBool isMPI;
2320: PetscFunctionBegin;
2322: PetscAssertPointer(array, 2);
2323: 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 */
2324: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2325: if (isMPI) {
2326: /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2327: PetscCall(MatDenseGetArrayAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2328: } else {
2329: PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);
2331: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayAndMemType_C", &fptr));
2332: if (fptr) {
2333: PetscCall((*fptr)(A, array, mtype));
2334: } else {
2335: PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2336: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2337: }
2338: }
2339: PetscFunctionReturn(PETSC_SUCCESS);
2340: }
2342: /*@C
2343: MatDenseRestoreArrayAndMemType - returns access to the array that is obtained by `MatDenseGetArrayAndMemType()`
2345: Logically Collective
2347: Input Parameters:
2348: + A - a dense matrix
2349: - array - pointer to the data
2351: Level: intermediate
2353: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2354: @*/
2355: PetscErrorCode MatDenseRestoreArrayAndMemType(Mat A, PetscScalar *array[])
2356: {
2357: PetscBool isMPI;
2359: PetscFunctionBegin;
2361: if (array) PetscAssertPointer(array, 2);
2362: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2363: if (isMPI) {
2364: PetscCall(MatDenseRestoreArrayAndMemType(((Mat_MPIDense *)A->data)->A, array));
2365: } else {
2366: PetscErrorCode (*fptr)(Mat, PetscScalar **);
2368: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayAndMemType_C", &fptr));
2369: if (fptr) {
2370: PetscCall((*fptr)(A, array));
2371: } else {
2372: PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2373: }
2374: if (array) *array = NULL;
2375: }
2376: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2377: PetscFunctionReturn(PETSC_SUCCESS);
2378: }
2380: /*@C
2381: MatDenseGetArrayReadAndMemType - gives read-only access to the array where the data for a `MATDENSE` matrix is stored
2383: Logically Collective
2385: Input Parameter:
2386: . A - a dense matrix
2388: Output Parameters:
2389: + array - pointer to the data
2390: - mtype - memory type of the returned pointer
2392: Level: intermediate
2394: Note:
2395: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2396: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2398: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`,
2399: `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2400: @*/
2401: PetscErrorCode MatDenseGetArrayReadAndMemType(Mat A, const PetscScalar *array[], PetscMemType *mtype)
2402: {
2403: PetscBool isMPI;
2405: PetscFunctionBegin;
2407: PetscAssertPointer(array, 2);
2408: 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 */
2409: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2410: if (isMPI) { /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2411: PetscCall(MatDenseGetArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2412: } else {
2413: PetscErrorCode (*fptr)(Mat, const PetscScalar **, PetscMemType *);
2415: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayReadAndMemType_C", &fptr));
2416: if (fptr) {
2417: PetscCall((*fptr)(A, array, mtype));
2418: } else {
2419: PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2420: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2421: }
2422: }
2423: PetscFunctionReturn(PETSC_SUCCESS);
2424: }
2426: /*@C
2427: MatDenseRestoreArrayReadAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()`
2429: Logically Collective
2431: Input Parameters:
2432: + A - a dense matrix
2433: - array - pointer to the data
2435: Level: intermediate
2437: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2438: @*/
2439: PetscErrorCode MatDenseRestoreArrayReadAndMemType(Mat A, const PetscScalar *array[])
2440: {
2441: PetscBool isMPI;
2443: PetscFunctionBegin;
2445: if (array) PetscAssertPointer(array, 2);
2446: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2447: if (isMPI) {
2448: PetscCall(MatDenseRestoreArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array));
2449: } else {
2450: PetscErrorCode (*fptr)(Mat, const PetscScalar **);
2452: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayReadAndMemType_C", &fptr));
2453: if (fptr) {
2454: PetscCall((*fptr)(A, array));
2455: } else {
2456: PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2457: }
2458: if (array) *array = NULL;
2459: }
2460: PetscFunctionReturn(PETSC_SUCCESS);
2461: }
2463: /*@C
2464: MatDenseGetArrayWriteAndMemType - gives write-only access to the array where the data for a `MATDENSE` matrix is stored
2466: Logically Collective
2468: Input Parameter:
2469: . A - a dense matrix
2471: Output Parameters:
2472: + array - pointer to the data
2473: - mtype - memory type of the returned pointer
2475: Level: intermediate
2477: Note:
2478: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2479: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2481: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWriteAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayRead()`,
2482: `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2483: @*/
2484: PetscErrorCode MatDenseGetArrayWriteAndMemType(Mat A, PetscScalar *array[], PetscMemType *mtype)
2485: {
2486: PetscBool isMPI;
2488: PetscFunctionBegin;
2490: PetscAssertPointer(array, 2);
2491: 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 */
2492: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2493: if (isMPI) {
2494: PetscCall(MatDenseGetArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2495: } else {
2496: PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);
2498: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayWriteAndMemType_C", &fptr));
2499: if (fptr) {
2500: PetscCall((*fptr)(A, array, mtype));
2501: } else {
2502: PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2503: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2504: }
2505: }
2506: PetscFunctionReturn(PETSC_SUCCESS);
2507: }
2509: /*@C
2510: MatDenseRestoreArrayWriteAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()`
2512: Logically Collective
2514: Input Parameters:
2515: + A - a dense matrix
2516: - array - pointer to the data
2518: Level: intermediate
2520: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2521: @*/
2522: PetscErrorCode MatDenseRestoreArrayWriteAndMemType(Mat A, PetscScalar *array[])
2523: {
2524: PetscBool isMPI;
2526: PetscFunctionBegin;
2528: if (array) PetscAssertPointer(array, 2);
2529: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2530: if (isMPI) {
2531: PetscCall(MatDenseRestoreArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array));
2532: } else {
2533: PetscErrorCode (*fptr)(Mat, PetscScalar **);
2535: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayWriteAndMemType_C", &fptr));
2536: if (fptr) {
2537: PetscCall((*fptr)(A, array));
2538: } else {
2539: PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2540: }
2541: if (array) *array = NULL;
2542: }
2543: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2544: PetscFunctionReturn(PETSC_SUCCESS);
2545: }
2547: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B)
2548: {
2549: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2550: PetscInt i, j, nrows, ncols, ldb;
2551: const PetscInt *irow, *icol;
2552: PetscScalar *av, *bv, *v = mat->v;
2553: Mat newmat;
2555: PetscFunctionBegin;
2556: PetscCall(ISGetIndices(isrow, &irow));
2557: PetscCall(ISGetIndices(iscol, &icol));
2558: PetscCall(ISGetLocalSize(isrow, &nrows));
2559: PetscCall(ISGetLocalSize(iscol, &ncols));
2561: /* Check submatrixcall */
2562: if (scall == MAT_REUSE_MATRIX) {
2563: PetscInt n_cols, n_rows;
2564: PetscCall(MatGetSize(*B, &n_rows, &n_cols));
2565: if (n_rows != nrows || n_cols != ncols) {
2566: /* resize the result matrix to match number of requested rows/columns */
2567: PetscCall(MatSetSizes(*B, nrows, ncols, nrows, ncols));
2568: }
2569: newmat = *B;
2570: } else {
2571: /* Create and fill new matrix */
2572: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &newmat));
2573: PetscCall(MatSetSizes(newmat, nrows, ncols, nrows, ncols));
2574: PetscCall(MatSetType(newmat, ((PetscObject)A)->type_name));
2575: PetscCall(MatSeqDenseSetPreallocation(newmat, NULL));
2576: }
2578: /* Now extract the data pointers and do the copy,column at a time */
2579: PetscCall(MatDenseGetArray(newmat, &bv));
2580: PetscCall(MatDenseGetLDA(newmat, &ldb));
2581: for (i = 0; i < ncols; i++) {
2582: av = v + mat->lda * icol[i];
2583: for (j = 0; j < nrows; j++) bv[j] = av[irow[j]];
2584: bv += ldb;
2585: }
2586: PetscCall(MatDenseRestoreArray(newmat, &bv));
2588: /* Assemble the matrices so that the correct flags are set */
2589: PetscCall(MatAssemblyBegin(newmat, MAT_FINAL_ASSEMBLY));
2590: PetscCall(MatAssemblyEnd(newmat, MAT_FINAL_ASSEMBLY));
2592: /* Free work space */
2593: PetscCall(ISRestoreIndices(isrow, &irow));
2594: PetscCall(ISRestoreIndices(iscol, &icol));
2595: *B = newmat;
2596: PetscFunctionReturn(PETSC_SUCCESS);
2597: }
2599: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[])
2600: {
2601: PetscInt i;
2603: PetscFunctionBegin;
2604: if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n, B));
2606: for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SeqDense(A, irow[i], icol[i], scall, &(*B)[i]));
2607: PetscFunctionReturn(PETSC_SUCCESS);
2608: }
2610: PetscErrorCode MatCopy_SeqDense(Mat A, Mat B, MatStructure str)
2611: {
2612: Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data;
2613: const PetscScalar *va;
2614: PetscScalar *vb;
2615: PetscInt lda1 = a->lda, lda2 = b->lda, m = A->rmap->n, n = A->cmap->n, j;
2617: PetscFunctionBegin;
2618: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2619: if (A->ops->copy != B->ops->copy) {
2620: PetscCall(MatCopy_Basic(A, B, str));
2621: PetscFunctionReturn(PETSC_SUCCESS);
2622: }
2623: PetscCheck(m == B->rmap->n && n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "size(B) != size(A)");
2624: PetscCall(MatDenseGetArrayRead(A, &va));
2625: PetscCall(MatDenseGetArray(B, &vb));
2626: if (lda1 > m || lda2 > m) {
2627: for (j = 0; j < n; j++) PetscCall(PetscArraycpy(vb + j * lda2, va + j * lda1, m));
2628: } else {
2629: PetscCall(PetscArraycpy(vb, va, A->rmap->n * A->cmap->n));
2630: }
2631: PetscCall(MatDenseRestoreArray(B, &vb));
2632: PetscCall(MatDenseRestoreArrayRead(A, &va));
2633: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2634: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2635: PetscFunctionReturn(PETSC_SUCCESS);
2636: }
2638: PetscErrorCode MatSetUp_SeqDense(Mat A)
2639: {
2640: PetscFunctionBegin;
2641: PetscCall(PetscLayoutSetUp(A->rmap));
2642: PetscCall(PetscLayoutSetUp(A->cmap));
2643: if (!A->preallocated) PetscCall(MatSeqDenseSetPreallocation(A, NULL));
2644: PetscFunctionReturn(PETSC_SUCCESS);
2645: }
2647: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2648: {
2649: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2650: PetscInt i, j;
2651: PetscInt min = PetscMin(A->rmap->n, A->cmap->n);
2652: PetscScalar *aa;
2654: PetscFunctionBegin;
2655: PetscCall(MatDenseGetArray(A, &aa));
2656: for (j = 0; j < A->cmap->n; j++) {
2657: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscConj(aa[i + j * mat->lda]);
2658: }
2659: PetscCall(MatDenseRestoreArray(A, &aa));
2660: if (mat->tau)
2661: for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
2662: PetscFunctionReturn(PETSC_SUCCESS);
2663: }
2665: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2666: {
2667: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2668: PetscInt i, j;
2669: PetscScalar *aa;
2671: PetscFunctionBegin;
2672: PetscCall(MatDenseGetArray(A, &aa));
2673: for (j = 0; j < A->cmap->n; j++) {
2674: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscRealPart(aa[i + j * mat->lda]);
2675: }
2676: PetscCall(MatDenseRestoreArray(A, &aa));
2677: PetscFunctionReturn(PETSC_SUCCESS);
2678: }
2680: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2681: {
2682: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2683: PetscInt i, j;
2684: PetscScalar *aa;
2686: PetscFunctionBegin;
2687: PetscCall(MatDenseGetArray(A, &aa));
2688: for (j = 0; j < A->cmap->n; j++) {
2689: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscImaginaryPart(aa[i + j * mat->lda]);
2690: }
2691: PetscCall(MatDenseRestoreArray(A, &aa));
2692: PetscFunctionReturn(PETSC_SUCCESS);
2693: }
2695: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2696: {
2697: PetscInt m = A->rmap->n, n = B->cmap->n;
2698: PetscBool cisdense = PETSC_FALSE;
2700: PetscFunctionBegin;
2701: PetscCall(MatSetSizes(C, m, n, m, n));
2702: #if defined(PETSC_HAVE_CUDA)
2703: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2704: #endif
2705: #if defined(PETSC_HAVE_HIP)
2706: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2707: #endif
2708: if (!cisdense) {
2709: PetscBool flg;
2711: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2712: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2713: }
2714: PetscCall(MatSetUp(C));
2715: PetscFunctionReturn(PETSC_SUCCESS);
2716: }
2718: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2719: {
2720: Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data, *c = (Mat_SeqDense *)C->data;
2721: PetscBLASInt m, n, k;
2722: const PetscScalar *av, *bv;
2723: PetscScalar *cv;
2724: PetscScalar _DOne = 1.0, _DZero = 0.0;
2726: PetscFunctionBegin;
2727: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2728: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2729: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2730: if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2731: PetscCall(MatDenseGetArrayRead(A, &av));
2732: PetscCall(MatDenseGetArrayRead(B, &bv));
2733: PetscCall(MatDenseGetArrayWrite(C, &cv));
2734: PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2735: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2736: PetscCall(MatDenseRestoreArrayRead(A, &av));
2737: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2738: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2739: PetscFunctionReturn(PETSC_SUCCESS);
2740: }
2742: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2743: {
2744: PetscInt m = A->rmap->n, n = B->rmap->n;
2745: PetscBool cisdense = PETSC_FALSE;
2747: PetscFunctionBegin;
2748: PetscCall(MatSetSizes(C, m, n, m, n));
2749: #if defined(PETSC_HAVE_CUDA)
2750: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2751: #endif
2752: #if defined(PETSC_HAVE_HIP)
2753: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2754: #endif
2755: if (!cisdense) {
2756: PetscBool flg;
2758: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2759: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2760: }
2761: PetscCall(MatSetUp(C));
2762: PetscFunctionReturn(PETSC_SUCCESS);
2763: }
2765: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2766: {
2767: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2768: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
2769: Mat_SeqDense *c = (Mat_SeqDense *)C->data;
2770: const PetscScalar *av, *bv;
2771: PetscScalar *cv;
2772: PetscBLASInt m, n, k;
2773: PetscScalar _DOne = 1.0, _DZero = 0.0;
2775: PetscFunctionBegin;
2776: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2777: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2778: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2779: if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2780: PetscCall(MatDenseGetArrayRead(A, &av));
2781: PetscCall(MatDenseGetArrayRead(B, &bv));
2782: PetscCall(MatDenseGetArrayWrite(C, &cv));
2783: PetscCallBLAS("BLASgemm", BLASgemm_("N", "T", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2784: PetscCall(MatDenseRestoreArrayRead(A, &av));
2785: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2786: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2787: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2788: PetscFunctionReturn(PETSC_SUCCESS);
2789: }
2791: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2792: {
2793: PetscInt m = A->cmap->n, n = B->cmap->n;
2794: PetscBool cisdense = PETSC_FALSE;
2796: PetscFunctionBegin;
2797: PetscCall(MatSetSizes(C, m, n, m, n));
2798: #if defined(PETSC_HAVE_CUDA)
2799: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2800: #endif
2801: #if defined(PETSC_HAVE_HIP)
2802: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2803: #endif
2804: if (!cisdense) {
2805: PetscBool flg;
2807: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2808: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2809: }
2810: PetscCall(MatSetUp(C));
2811: PetscFunctionReturn(PETSC_SUCCESS);
2812: }
2814: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2815: {
2816: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2817: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
2818: Mat_SeqDense *c = (Mat_SeqDense *)C->data;
2819: const PetscScalar *av, *bv;
2820: PetscScalar *cv;
2821: PetscBLASInt m, n, k;
2822: PetscScalar _DOne = 1.0, _DZero = 0.0;
2824: PetscFunctionBegin;
2825: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2826: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2827: PetscCall(PetscBLASIntCast(A->rmap->n, &k));
2828: if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2829: PetscCall(MatDenseGetArrayRead(A, &av));
2830: PetscCall(MatDenseGetArrayRead(B, &bv));
2831: PetscCall(MatDenseGetArrayWrite(C, &cv));
2832: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2833: PetscCall(MatDenseRestoreArrayRead(A, &av));
2834: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2835: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2836: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2837: PetscFunctionReturn(PETSC_SUCCESS);
2838: }
2840: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2841: {
2842: PetscFunctionBegin;
2843: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2844: C->ops->productsymbolic = MatProductSymbolic_AB;
2845: PetscFunctionReturn(PETSC_SUCCESS);
2846: }
2848: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2849: {
2850: PetscFunctionBegin;
2851: C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2852: C->ops->productsymbolic = MatProductSymbolic_AtB;
2853: PetscFunctionReturn(PETSC_SUCCESS);
2854: }
2856: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2857: {
2858: PetscFunctionBegin;
2859: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2860: C->ops->productsymbolic = MatProductSymbolic_ABt;
2861: PetscFunctionReturn(PETSC_SUCCESS);
2862: }
2864: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2865: {
2866: Mat_Product *product = C->product;
2868: PetscFunctionBegin;
2869: switch (product->type) {
2870: case MATPRODUCT_AB:
2871: PetscCall(MatProductSetFromOptions_SeqDense_AB(C));
2872: break;
2873: case MATPRODUCT_AtB:
2874: PetscCall(MatProductSetFromOptions_SeqDense_AtB(C));
2875: break;
2876: case MATPRODUCT_ABt:
2877: PetscCall(MatProductSetFromOptions_SeqDense_ABt(C));
2878: break;
2879: default:
2880: break;
2881: }
2882: PetscFunctionReturn(PETSC_SUCCESS);
2883: }
2885: static PetscErrorCode MatGetRowMax_SeqDense(Mat A, Vec v, PetscInt idx[])
2886: {
2887: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2888: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2889: PetscScalar *x;
2890: const PetscScalar *aa;
2892: PetscFunctionBegin;
2893: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2894: PetscCall(VecGetArray(v, &x));
2895: PetscCall(VecGetLocalSize(v, &p));
2896: PetscCall(MatDenseGetArrayRead(A, &aa));
2897: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2898: for (i = 0; i < m; i++) {
2899: x[i] = aa[i];
2900: if (idx) idx[i] = 0;
2901: for (j = 1; j < n; j++) {
2902: if (PetscRealPart(x[i]) < PetscRealPart(aa[i + a->lda * j])) {
2903: x[i] = aa[i + a->lda * j];
2904: if (idx) idx[i] = j;
2905: }
2906: }
2907: }
2908: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2909: PetscCall(VecRestoreArray(v, &x));
2910: PetscFunctionReturn(PETSC_SUCCESS);
2911: }
2913: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A, Vec v, PetscInt idx[])
2914: {
2915: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2916: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2917: PetscScalar *x;
2918: PetscReal atmp;
2919: const PetscScalar *aa;
2921: PetscFunctionBegin;
2922: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2923: PetscCall(VecGetArray(v, &x));
2924: PetscCall(VecGetLocalSize(v, &p));
2925: PetscCall(MatDenseGetArrayRead(A, &aa));
2926: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2927: for (i = 0; i < m; i++) {
2928: x[i] = PetscAbsScalar(aa[i]);
2929: for (j = 1; j < n; j++) {
2930: atmp = PetscAbsScalar(aa[i + a->lda * j]);
2931: if (PetscAbsScalar(x[i]) < atmp) {
2932: x[i] = atmp;
2933: if (idx) idx[i] = j;
2934: }
2935: }
2936: }
2937: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2938: PetscCall(VecRestoreArray(v, &x));
2939: PetscFunctionReturn(PETSC_SUCCESS);
2940: }
2942: static PetscErrorCode MatGetRowMin_SeqDense(Mat A, Vec v, PetscInt idx[])
2943: {
2944: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2945: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2946: PetscScalar *x;
2947: const PetscScalar *aa;
2949: PetscFunctionBegin;
2950: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2951: PetscCall(MatDenseGetArrayRead(A, &aa));
2952: PetscCall(VecGetArray(v, &x));
2953: PetscCall(VecGetLocalSize(v, &p));
2954: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2955: for (i = 0; i < m; i++) {
2956: x[i] = aa[i];
2957: if (idx) idx[i] = 0;
2958: for (j = 1; j < n; j++) {
2959: if (PetscRealPart(x[i]) > PetscRealPart(aa[i + a->lda * j])) {
2960: x[i] = aa[i + a->lda * j];
2961: if (idx) idx[i] = j;
2962: }
2963: }
2964: }
2965: PetscCall(VecRestoreArray(v, &x));
2966: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2967: PetscFunctionReturn(PETSC_SUCCESS);
2968: }
2970: PetscErrorCode MatGetColumnVector_SeqDense(Mat A, Vec v, PetscInt col)
2971: {
2972: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2973: PetscScalar *x;
2974: const PetscScalar *aa;
2976: PetscFunctionBegin;
2977: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2978: PetscCall(MatDenseGetArrayRead(A, &aa));
2979: PetscCall(VecGetArray(v, &x));
2980: PetscCall(PetscArraycpy(x, aa + col * a->lda, A->rmap->n));
2981: PetscCall(VecRestoreArray(v, &x));
2982: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2983: PetscFunctionReturn(PETSC_SUCCESS);
2984: }
2986: PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A, PetscInt type, PetscReal *reductions)
2987: {
2988: PetscInt i, j, m, n;
2989: const PetscScalar *a;
2991: PetscFunctionBegin;
2992: PetscCall(MatGetSize(A, &m, &n));
2993: PetscCall(PetscArrayzero(reductions, n));
2994: PetscCall(MatDenseGetArrayRead(A, &a));
2995: if (type == NORM_2) {
2996: for (i = 0; i < n; i++) {
2997: for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j] * a[j]);
2998: a = PetscSafePointerPlusOffset(a, m);
2999: }
3000: } else if (type == NORM_1) {
3001: for (i = 0; i < n; i++) {
3002: for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j]);
3003: a = PetscSafePointerPlusOffset(a, m);
3004: }
3005: } else if (type == NORM_INFINITY) {
3006: for (i = 0; i < n; i++) {
3007: for (j = 0; j < m; j++) reductions[i] = PetscMax(PetscAbsScalar(a[j]), reductions[i]);
3008: a = PetscSafePointerPlusOffset(a, m);
3009: }
3010: } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
3011: for (i = 0; i < n; i++) {
3012: for (j = 0; j < m; j++) reductions[i] += PetscRealPart(a[j]);
3013: a = PetscSafePointerPlusOffset(a, m);
3014: }
3015: } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3016: for (i = 0; i < n; i++) {
3017: for (j = 0; j < m; j++) reductions[i] += PetscImaginaryPart(a[j]);
3018: a = PetscSafePointerPlusOffset(a, m);
3019: }
3020: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown reduction type");
3021: PetscCall(MatDenseRestoreArrayRead(A, &a));
3022: if (type == NORM_2) {
3023: for (i = 0; i < n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
3024: } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3025: for (i = 0; i < n; i++) reductions[i] /= m;
3026: }
3027: PetscFunctionReturn(PETSC_SUCCESS);
3028: }
3030: PetscErrorCode MatSetRandom_SeqDense(Mat x, PetscRandom rctx)
3031: {
3032: PetscScalar *a;
3033: PetscInt lda, m, n, i, j;
3035: PetscFunctionBegin;
3036: PetscCall(MatGetSize(x, &m, &n));
3037: PetscCall(MatDenseGetLDA(x, &lda));
3038: PetscCall(MatDenseGetArrayWrite(x, &a));
3039: for (j = 0; j < n; j++) {
3040: for (i = 0; i < m; i++) PetscCall(PetscRandomGetValue(rctx, a + j * lda + i));
3041: }
3042: PetscCall(MatDenseRestoreArrayWrite(x, &a));
3043: PetscFunctionReturn(PETSC_SUCCESS);
3044: }
3046: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A, PetscBool *missing, PetscInt *d)
3047: {
3048: PetscFunctionBegin;
3049: *missing = PETSC_FALSE;
3050: PetscFunctionReturn(PETSC_SUCCESS);
3051: }
3053: /* vals is not const */
3054: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A, PetscInt col, PetscScalar **vals)
3055: {
3056: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3057: PetscScalar *v;
3059: PetscFunctionBegin;
3060: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3061: PetscCall(MatDenseGetArray(A, &v));
3062: *vals = v + col * a->lda;
3063: PetscCall(MatDenseRestoreArray(A, &v));
3064: PetscFunctionReturn(PETSC_SUCCESS);
3065: }
3067: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A, PetscScalar **vals)
3068: {
3069: PetscFunctionBegin;
3070: if (vals) *vals = NULL; /* user cannot accidentally use the array later */
3071: PetscFunctionReturn(PETSC_SUCCESS);
3072: }
3074: static struct _MatOps MatOps_Values = {MatSetValues_SeqDense,
3075: MatGetRow_SeqDense,
3076: MatRestoreRow_SeqDense,
3077: MatMult_SeqDense,
3078: /* 4*/ MatMultAdd_SeqDense,
3079: MatMultTranspose_SeqDense,
3080: MatMultTransposeAdd_SeqDense,
3081: NULL,
3082: NULL,
3083: NULL,
3084: /* 10*/ NULL,
3085: MatLUFactor_SeqDense,
3086: MatCholeskyFactor_SeqDense,
3087: MatSOR_SeqDense,
3088: MatTranspose_SeqDense,
3089: /* 15*/ MatGetInfo_SeqDense,
3090: MatEqual_SeqDense,
3091: MatGetDiagonal_SeqDense,
3092: MatDiagonalScale_SeqDense,
3093: MatNorm_SeqDense,
3094: /* 20*/ NULL,
3095: NULL,
3096: MatSetOption_SeqDense,
3097: MatZeroEntries_SeqDense,
3098: /* 24*/ MatZeroRows_SeqDense,
3099: NULL,
3100: NULL,
3101: NULL,
3102: NULL,
3103: /* 29*/ MatSetUp_SeqDense,
3104: NULL,
3105: NULL,
3106: NULL,
3107: NULL,
3108: /* 34*/ MatDuplicate_SeqDense,
3109: NULL,
3110: NULL,
3111: NULL,
3112: NULL,
3113: /* 39*/ MatAXPY_SeqDense,
3114: MatCreateSubMatrices_SeqDense,
3115: NULL,
3116: MatGetValues_SeqDense,
3117: MatCopy_SeqDense,
3118: /* 44*/ MatGetRowMax_SeqDense,
3119: MatScale_SeqDense,
3120: MatShift_SeqDense,
3121: NULL,
3122: MatZeroRowsColumns_SeqDense,
3123: /* 49*/ MatSetRandom_SeqDense,
3124: NULL,
3125: NULL,
3126: NULL,
3127: NULL,
3128: /* 54*/ NULL,
3129: NULL,
3130: NULL,
3131: NULL,
3132: NULL,
3133: /* 59*/ MatCreateSubMatrix_SeqDense,
3134: MatDestroy_SeqDense,
3135: MatView_SeqDense,
3136: NULL,
3137: NULL,
3138: /* 64*/ NULL,
3139: NULL,
3140: NULL,
3141: NULL,
3142: NULL,
3143: /* 69*/ MatGetRowMaxAbs_SeqDense,
3144: NULL,
3145: NULL,
3146: NULL,
3147: NULL,
3148: /* 74*/ NULL,
3149: NULL,
3150: NULL,
3151: NULL,
3152: NULL,
3153: /* 79*/ NULL,
3154: NULL,
3155: NULL,
3156: NULL,
3157: /* 83*/ MatLoad_SeqDense,
3158: MatIsSymmetric_SeqDense,
3159: MatIsHermitian_SeqDense,
3160: NULL,
3161: NULL,
3162: NULL,
3163: /* 89*/ NULL,
3164: NULL,
3165: MatMatMultNumeric_SeqDense_SeqDense,
3166: NULL,
3167: NULL,
3168: /* 94*/ NULL,
3169: NULL,
3170: NULL,
3171: MatMatTransposeMultNumeric_SeqDense_SeqDense,
3172: NULL,
3173: /* 99*/ MatProductSetFromOptions_SeqDense,
3174: NULL,
3175: NULL,
3176: MatConjugate_SeqDense,
3177: NULL,
3178: /*104*/ NULL,
3179: MatRealPart_SeqDense,
3180: MatImaginaryPart_SeqDense,
3181: NULL,
3182: NULL,
3183: /*109*/ NULL,
3184: NULL,
3185: MatGetRowMin_SeqDense,
3186: MatGetColumnVector_SeqDense,
3187: MatMissingDiagonal_SeqDense,
3188: /*114*/ NULL,
3189: NULL,
3190: NULL,
3191: NULL,
3192: NULL,
3193: /*119*/ NULL,
3194: NULL,
3195: MatMultHermitianTranspose_SeqDense,
3196: MatMultHermitianTransposeAdd_SeqDense,
3197: NULL,
3198: /*124*/ NULL,
3199: MatGetColumnReductions_SeqDense,
3200: NULL,
3201: NULL,
3202: NULL,
3203: /*129*/ NULL,
3204: NULL,
3205: NULL,
3206: MatTransposeMatMultNumeric_SeqDense_SeqDense,
3207: NULL,
3208: /*134*/ NULL,
3209: NULL,
3210: NULL,
3211: NULL,
3212: NULL,
3213: /*139*/ NULL,
3214: NULL,
3215: NULL,
3216: NULL,
3217: NULL,
3218: MatCreateMPIMatConcatenateSeqMat_SeqDense,
3219: /*145*/ NULL,
3220: NULL,
3221: NULL,
3222: NULL,
3223: NULL,
3224: /*150*/ NULL,
3225: NULL,
3226: NULL,
3227: NULL,
3228: NULL,
3229: /*155*/ NULL,
3230: NULL};
3232: /*@
3233: MatCreateSeqDense - Creates a `MATSEQDENSE` that
3234: is stored in column major order (the usual Fortran format).
3236: Collective
3238: Input Parameters:
3239: + comm - MPI communicator, set to `PETSC_COMM_SELF`
3240: . m - number of rows
3241: . n - number of columns
3242: - data - optional location of matrix data in column major order. Use `NULL` for PETSc
3243: to control all matrix memory allocation.
3245: Output Parameter:
3246: . A - the matrix
3248: Level: intermediate
3250: Note:
3251: The data input variable is intended primarily for Fortran programmers
3252: who wish to allocate their own matrix memory space. Most users should
3253: set `data` = `NULL`.
3255: Developer Note:
3256: Many of the matrix operations for this variant use the BLAS and LAPACK routines.
3258: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`
3259: @*/
3260: PetscErrorCode MatCreateSeqDense(MPI_Comm comm, PetscInt m, PetscInt n, PetscScalar data[], Mat *A)
3261: {
3262: PetscFunctionBegin;
3263: PetscCall(MatCreate(comm, A));
3264: PetscCall(MatSetSizes(*A, m, n, m, n));
3265: PetscCall(MatSetType(*A, MATSEQDENSE));
3266: PetscCall(MatSeqDenseSetPreallocation(*A, data));
3267: PetscFunctionReturn(PETSC_SUCCESS);
3268: }
3270: /*@
3271: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements of a `MATSEQDENSE` matrix
3273: Collective
3275: Input Parameters:
3276: + B - the matrix
3277: - data - the array (or `NULL`)
3279: Level: intermediate
3281: Note:
3282: The data input variable is intended primarily for Fortran programmers
3283: who wish to allocate their own matrix memory space. Most users should
3284: need not call this routine.
3286: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`, `MatDenseSetLDA()`
3287: @*/
3288: PetscErrorCode MatSeqDenseSetPreallocation(Mat B, PetscScalar data[])
3289: {
3290: PetscFunctionBegin;
3292: PetscTryMethod(B, "MatSeqDenseSetPreallocation_C", (Mat, PetscScalar[]), (B, data));
3293: PetscFunctionReturn(PETSC_SUCCESS);
3294: }
3296: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B, PetscScalar *data)
3297: {
3298: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
3300: PetscFunctionBegin;
3301: PetscCheck(!b->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3302: B->preallocated = PETSC_TRUE;
3304: PetscCall(PetscLayoutSetUp(B->rmap));
3305: PetscCall(PetscLayoutSetUp(B->cmap));
3307: if (b->lda <= 0) PetscCall(PetscBLASIntCast(B->rmap->n, &b->lda));
3309: if (!data) { /* petsc-allocated storage */
3310: if (!b->user_alloc) PetscCall(PetscFree(b->v));
3311: PetscCall(PetscCalloc1((size_t)b->lda * B->cmap->n, &b->v));
3313: b->user_alloc = PETSC_FALSE;
3314: } else { /* user-allocated storage */
3315: if (!b->user_alloc) PetscCall(PetscFree(b->v));
3316: b->v = data;
3317: b->user_alloc = PETSC_TRUE;
3318: }
3319: B->assembled = PETSC_TRUE;
3320: PetscFunctionReturn(PETSC_SUCCESS);
3321: }
3323: #if defined(PETSC_HAVE_ELEMENTAL)
3324: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
3325: {
3326: Mat mat_elemental;
3327: const PetscScalar *array;
3328: PetscScalar *v_colwise;
3329: PetscInt M = A->rmap->N, N = A->cmap->N, i, j, k, *rows, *cols;
3331: PetscFunctionBegin;
3332: PetscCall(PetscMalloc3(M * N, &v_colwise, M, &rows, N, &cols));
3333: PetscCall(MatDenseGetArrayRead(A, &array));
3334: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
3335: k = 0;
3336: for (j = 0; j < N; j++) {
3337: cols[j] = j;
3338: for (i = 0; i < M; i++) v_colwise[j * M + i] = array[k++];
3339: }
3340: for (i = 0; i < M; i++) rows[i] = i;
3341: PetscCall(MatDenseRestoreArrayRead(A, &array));
3343: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
3344: PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
3345: PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
3346: PetscCall(MatSetUp(mat_elemental));
3348: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
3349: PetscCall(MatSetValues(mat_elemental, M, rows, N, cols, v_colwise, ADD_VALUES));
3350: PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
3351: PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
3352: PetscCall(PetscFree3(v_colwise, rows, cols));
3354: if (reuse == MAT_INPLACE_MATRIX) {
3355: PetscCall(MatHeaderReplace(A, &mat_elemental));
3356: } else {
3357: *newmat = mat_elemental;
3358: }
3359: PetscFunctionReturn(PETSC_SUCCESS);
3360: }
3361: #endif
3363: PetscErrorCode MatDenseSetLDA_SeqDense(Mat B, PetscInt lda)
3364: {
3365: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
3366: PetscBool data;
3368: PetscFunctionBegin;
3369: data = (B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE;
3370: PetscCheck(b->user_alloc || !data || b->lda == lda, PETSC_COMM_SELF, PETSC_ERR_ORDER, "LDA cannot be changed after allocation of internal storage");
3371: 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);
3372: PetscCall(PetscBLASIntCast(lda, &b->lda));
3373: PetscFunctionReturn(PETSC_SUCCESS);
3374: }
3376: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
3377: {
3378: PetscFunctionBegin;
3379: PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIDense(comm, inmat, n, scall, outmat));
3380: PetscFunctionReturn(PETSC_SUCCESS);
3381: }
3383: PetscErrorCode MatDenseCreateColumnVec_Private(Mat A, Vec *v)
3384: {
3385: PetscBool isstd, iskok, iscuda, iship;
3386: PetscMPIInt size;
3387: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
3388: /* we pass the data of A, to prevent allocating needless GPU memory the first time VecCUPMPlaceArray is called. */
3389: const PetscScalar *a;
3390: #endif
3392: PetscFunctionBegin;
3393: *v = NULL;
3394: PetscCall(PetscStrcmpAny(A->defaultvectype, &isstd, VECSTANDARD, VECSEQ, VECMPI, ""));
3395: PetscCall(PetscStrcmpAny(A->defaultvectype, &iskok, VECKOKKOS, VECSEQKOKKOS, VECMPIKOKKOS, ""));
3396: PetscCall(PetscStrcmpAny(A->defaultvectype, &iscuda, VECCUDA, VECSEQCUDA, VECMPICUDA, ""));
3397: PetscCall(PetscStrcmpAny(A->defaultvectype, &iship, VECHIP, VECSEQHIP, VECMPIHIP, ""));
3398: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
3399: if (isstd) {
3400: if (size > 1) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, v));
3401: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, v));
3402: } else if (iskok) {
3403: PetscCheck(PetscDefined(HAVE_KOKKOS_KERNELS), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using KOKKOS kernels support");
3404: #if PetscDefined(HAVE_KOKKOS_KERNELS)
3405: if (size > 1) PetscCall(VecCreateMPIKokkosWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, v));
3406: else PetscCall(VecCreateSeqKokkosWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, v));
3407: #endif
3408: } else if (iscuda) {
3409: PetscCheck(PetscDefined(HAVE_CUDA), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using CUDA support");
3410: #if PetscDefined(HAVE_CUDA)
3411: PetscCall(MatDenseCUDAGetArrayRead(A, &a));
3412: if (size > 1) PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, a, v));
3413: else PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, a, v));
3414: #endif
3415: } else if (iship) {
3416: PetscCheck(PetscDefined(HAVE_HIP), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using HIP support");
3417: #if PetscDefined(HAVE_HIP)
3418: PetscCall(MatDenseHIPGetArrayRead(A, &a));
3419: if (size > 1) PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, a, v));
3420: else PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, a, v));
3421: #endif
3422: }
3423: PetscCheck(*v, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not coded for type %s", A->defaultvectype);
3424: PetscFunctionReturn(PETSC_SUCCESS);
3425: }
3427: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3428: {
3429: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3431: PetscFunctionBegin;
3432: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3433: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3434: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3435: a->vecinuse = col + 1;
3436: PetscCall(MatDenseGetArray(A, (PetscScalar **)&a->ptrinuse));
3437: PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda));
3438: *v = a->cvec;
3439: PetscFunctionReturn(PETSC_SUCCESS);
3440: }
3442: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3443: {
3444: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3446: PetscFunctionBegin;
3447: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3448: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3449: VecCheckAssembled(a->cvec);
3450: a->vecinuse = 0;
3451: PetscCall(MatDenseRestoreArray(A, (PetscScalar **)&a->ptrinuse));
3452: PetscCall(VecResetArray(a->cvec));
3453: if (v) *v = NULL;
3454: PetscFunctionReturn(PETSC_SUCCESS);
3455: }
3457: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3458: {
3459: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3461: PetscFunctionBegin;
3462: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3463: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3464: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3465: a->vecinuse = col + 1;
3466: PetscCall(MatDenseGetArrayRead(A, &a->ptrinuse));
3467: PetscCall(VecPlaceArray(a->cvec, PetscSafePointerPlusOffset(a->ptrinuse, (size_t)col * (size_t)a->lda)));
3468: PetscCall(VecLockReadPush(a->cvec));
3469: *v = a->cvec;
3470: PetscFunctionReturn(PETSC_SUCCESS);
3471: }
3473: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3474: {
3475: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3477: PetscFunctionBegin;
3478: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3479: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3480: VecCheckAssembled(a->cvec);
3481: a->vecinuse = 0;
3482: PetscCall(MatDenseRestoreArrayRead(A, &a->ptrinuse));
3483: PetscCall(VecLockReadPop(a->cvec));
3484: PetscCall(VecResetArray(a->cvec));
3485: if (v) *v = NULL;
3486: PetscFunctionReturn(PETSC_SUCCESS);
3487: }
3489: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3490: {
3491: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3493: PetscFunctionBegin;
3494: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3495: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3496: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3497: a->vecinuse = col + 1;
3498: PetscCall(MatDenseGetArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3499: PetscCall(VecPlaceArray(a->cvec, PetscSafePointerPlusOffset(a->ptrinuse, (size_t)col * (size_t)a->lda)));
3500: *v = a->cvec;
3501: PetscFunctionReturn(PETSC_SUCCESS);
3502: }
3504: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3505: {
3506: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3508: PetscFunctionBegin;
3509: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3510: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3511: VecCheckAssembled(a->cvec);
3512: a->vecinuse = 0;
3513: PetscCall(MatDenseRestoreArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3514: PetscCall(VecResetArray(a->cvec));
3515: if (v) *v = NULL;
3516: PetscFunctionReturn(PETSC_SUCCESS);
3517: }
3519: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3520: {
3521: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3523: PetscFunctionBegin;
3524: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3525: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3526: if (a->cmat && (cend - cbegin != a->cmat->cmap->N || rend - rbegin != a->cmat->rmap->N)) PetscCall(MatDestroy(&a->cmat));
3527: if (!a->cmat) {
3528: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), rend - rbegin, PETSC_DECIDE, rend - rbegin, cend - cbegin, PetscSafePointerPlusOffset(a->v, rbegin + (size_t)cbegin * a->lda), &a->cmat));
3529: } else {
3530: PetscCall(MatDensePlaceArray(a->cmat, PetscSafePointerPlusOffset(a->v, rbegin + (size_t)cbegin * a->lda)));
3531: }
3532: PetscCall(MatDenseSetLDA(a->cmat, a->lda));
3533: a->matinuse = cbegin + 1;
3534: *v = a->cmat;
3535: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
3536: A->offloadmask = PETSC_OFFLOAD_CPU;
3537: #endif
3538: PetscFunctionReturn(PETSC_SUCCESS);
3539: }
3541: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A, Mat *v)
3542: {
3543: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3545: PetscFunctionBegin;
3546: PetscCheck(a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetSubMatrix() first");
3547: PetscCheck(a->cmat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column matrix");
3548: PetscCheck(*v == a->cmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not the matrix obtained from MatDenseGetSubMatrix()");
3549: a->matinuse = 0;
3550: PetscCall(MatDenseResetArray(a->cmat));
3551: if (v) *v = NULL;
3552: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
3553: A->offloadmask = PETSC_OFFLOAD_CPU;
3554: #endif
3555: PetscFunctionReturn(PETSC_SUCCESS);
3556: }
3558: /*MC
3559: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
3561: Options Database Key:
3562: . -mat_type seqdense - sets the matrix type to `MATSEQDENSE` during a call to `MatSetFromOptions()`
3564: Level: beginner
3566: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreateSeqDense()`
3567: M*/
3568: PetscErrorCode MatCreate_SeqDense(Mat B)
3569: {
3570: Mat_SeqDense *b;
3571: PetscMPIInt size;
3573: PetscFunctionBegin;
3574: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
3575: PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1");
3577: PetscCall(PetscNew(&b));
3578: B->data = (void *)b;
3579: B->ops[0] = MatOps_Values;
3581: b->roworiented = PETSC_TRUE;
3583: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactor_C", MatQRFactor_SeqDense));
3584: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetLDA_C", MatDenseGetLDA_SeqDense));
3585: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseSetLDA_C", MatDenseSetLDA_SeqDense));
3586: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArray_C", MatDenseGetArray_SeqDense));
3587: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArray_C", MatDenseRestoreArray_SeqDense));
3588: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDensePlaceArray_C", MatDensePlaceArray_SeqDense));
3589: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseResetArray_C", MatDenseResetArray_SeqDense));
3590: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseReplaceArray_C", MatDenseReplaceArray_SeqDense));
3591: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayRead_C", MatDenseGetArray_SeqDense));
3592: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayRead_C", MatDenseRestoreArray_SeqDense));
3593: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayWrite_C", MatDenseGetArray_SeqDense));
3594: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayWrite_C", MatDenseRestoreArray_SeqDense));
3595: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqaij_C", MatConvert_SeqDense_SeqAIJ));
3596: #if defined(PETSC_HAVE_ELEMENTAL)
3597: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_elemental_C", MatConvert_SeqDense_Elemental));
3598: #endif
3599: #if defined(PETSC_HAVE_SCALAPACK)
3600: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_scalapack_C", MatConvert_Dense_ScaLAPACK));
3601: #endif
3602: #if defined(PETSC_HAVE_CUDA)
3603: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensecuda_C", MatConvert_SeqDense_SeqDenseCUDA));
3604: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3605: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdense_C", MatProductSetFromOptions_SeqDense));
3606: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3607: #endif
3608: #if defined(PETSC_HAVE_HIP)
3609: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensehip_C", MatConvert_SeqDense_SeqDenseHIP));
3610: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3611: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdense_C", MatProductSetFromOptions_SeqDense));
3612: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3613: #endif
3614: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqDenseSetPreallocation_C", MatSeqDenseSetPreallocation_SeqDense));
3615: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqdense_C", MatProductSetFromOptions_SeqAIJ_SeqDense));
3616: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdense_C", MatProductSetFromOptions_SeqDense));
3617: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));
3618: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqsbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));
3620: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumn_C", MatDenseGetColumn_SeqDense));
3621: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumn_C", MatDenseRestoreColumn_SeqDense));
3622: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVec_C", MatDenseGetColumnVec_SeqDense));
3623: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVec_C", MatDenseRestoreColumnVec_SeqDense));
3624: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecRead_C", MatDenseGetColumnVecRead_SeqDense));
3625: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecRead_C", MatDenseRestoreColumnVecRead_SeqDense));
3626: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecWrite_C", MatDenseGetColumnVecWrite_SeqDense));
3627: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecWrite_C", MatDenseRestoreColumnVecWrite_SeqDense));
3628: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetSubMatrix_C", MatDenseGetSubMatrix_SeqDense));
3629: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreSubMatrix_C", MatDenseRestoreSubMatrix_SeqDense));
3630: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultAddColumnRange_C", MatMultAddColumnRange_SeqDense));
3631: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultHermitianTransposeColumnRange_C", MatMultHermitianTransposeColumnRange_SeqDense));
3632: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultHermitianTransposeAddColumnRange_C", MatMultHermitianTransposeAddColumnRange_SeqDense));
3633: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQDENSE));
3634: PetscFunctionReturn(PETSC_SUCCESS);
3635: }
3637: /*@C
3638: 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.
3640: Not Collective
3642: Input Parameters:
3643: + A - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3644: - col - column index
3646: Output Parameter:
3647: . vals - pointer to the data
3649: Level: intermediate
3651: Note:
3652: Use `MatDenseGetColumnVec()` to get access to a column of a `MATDENSE` treated as a `Vec`
3654: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreColumn()`, `MatDenseGetColumnVec()`
3655: @*/
3656: PetscErrorCode MatDenseGetColumn(Mat A, PetscInt col, PetscScalar *vals[])
3657: {
3658: PetscFunctionBegin;
3661: PetscAssertPointer(vals, 3);
3662: PetscUseMethod(A, "MatDenseGetColumn_C", (Mat, PetscInt, PetscScalar **), (A, col, vals));
3663: PetscFunctionReturn(PETSC_SUCCESS);
3664: }
3666: /*@C
3667: MatDenseRestoreColumn - returns access to a column of a `MATDENSE` matrix which is returned by `MatDenseGetColumn()`.
3669: Not Collective
3671: Input Parameters:
3672: + A - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3673: - vals - pointer to the data (may be `NULL`)
3675: Level: intermediate
3677: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetColumn()`
3678: @*/
3679: PetscErrorCode MatDenseRestoreColumn(Mat A, PetscScalar *vals[])
3680: {
3681: PetscFunctionBegin;
3683: PetscAssertPointer(vals, 2);
3684: PetscUseMethod(A, "MatDenseRestoreColumn_C", (Mat, PetscScalar **), (A, vals));
3685: PetscFunctionReturn(PETSC_SUCCESS);
3686: }
3688: /*@
3689: MatDenseGetColumnVec - Gives read-write access to a column of a `MATDENSE` matrix, represented as a `Vec`.
3691: Collective
3693: Input Parameters:
3694: + A - the `Mat` object
3695: - col - the column index
3697: Output Parameter:
3698: . v - the vector
3700: Level: intermediate
3702: Notes:
3703: The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVec()` when the vector is no longer needed.
3705: Use `MatDenseGetColumnVecRead()` to obtain read-only access or `MatDenseGetColumnVecWrite()` for write-only access.
3707: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`, `MatDenseGetColumn()`
3708: @*/
3709: PetscErrorCode MatDenseGetColumnVec(Mat A, PetscInt col, Vec *v)
3710: {
3711: PetscFunctionBegin;
3715: PetscAssertPointer(v, 3);
3716: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3717: 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);
3718: PetscUseMethod(A, "MatDenseGetColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3719: PetscFunctionReturn(PETSC_SUCCESS);
3720: }
3722: /*@
3723: MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVec()`.
3725: Collective
3727: Input Parameters:
3728: + A - the `Mat` object
3729: . col - the column index
3730: - v - the `Vec` object (may be `NULL`)
3732: Level: intermediate
3734: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3735: @*/
3736: PetscErrorCode MatDenseRestoreColumnVec(Mat A, PetscInt col, Vec *v)
3737: {
3738: PetscFunctionBegin;
3742: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3743: 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);
3744: PetscUseMethod(A, "MatDenseRestoreColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3745: PetscFunctionReturn(PETSC_SUCCESS);
3746: }
3748: /*@
3749: MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a `Vec`.
3751: Collective
3753: Input Parameters:
3754: + A - the `Mat` object
3755: - col - the column index
3757: Output Parameter:
3758: . v - the vector
3760: Level: intermediate
3762: Notes:
3763: The vector is owned by PETSc and users cannot modify it.
3765: Users need to call `MatDenseRestoreColumnVecRead()` when the vector is no longer needed.
3767: Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecWrite()` for write-only access.
3769: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3770: @*/
3771: PetscErrorCode MatDenseGetColumnVecRead(Mat A, PetscInt col, Vec *v)
3772: {
3773: PetscFunctionBegin;
3777: PetscAssertPointer(v, 3);
3778: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3779: 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);
3780: PetscUseMethod(A, "MatDenseGetColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3781: PetscFunctionReturn(PETSC_SUCCESS);
3782: }
3784: /*@
3785: MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVecRead()`.
3787: Collective
3789: Input Parameters:
3790: + A - the `Mat` object
3791: . col - the column index
3792: - v - the `Vec` object (may be `NULL`)
3794: Level: intermediate
3796: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecWrite()`
3797: @*/
3798: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A, PetscInt col, Vec *v)
3799: {
3800: PetscFunctionBegin;
3804: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3805: 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);
3806: PetscUseMethod(A, "MatDenseRestoreColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3807: PetscFunctionReturn(PETSC_SUCCESS);
3808: }
3810: /*@
3811: MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a `Vec`.
3813: Collective
3815: Input Parameters:
3816: + A - the `Mat` object
3817: - col - the column index
3819: Output Parameter:
3820: . v - the vector
3822: Level: intermediate
3824: Notes:
3825: The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVecWrite()` when the vector is no longer needed.
3827: Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecRead()` for read-only access.
3829: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3830: @*/
3831: PetscErrorCode MatDenseGetColumnVecWrite(Mat A, PetscInt col, Vec *v)
3832: {
3833: PetscFunctionBegin;
3837: PetscAssertPointer(v, 3);
3838: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3839: 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);
3840: PetscUseMethod(A, "MatDenseGetColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3841: PetscFunctionReturn(PETSC_SUCCESS);
3842: }
3844: /*@
3845: MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVecWrite()`.
3847: Collective
3849: Input Parameters:
3850: + A - the `Mat` object
3851: . col - the column index
3852: - v - the `Vec` object (may be `NULL`)
3854: Level: intermediate
3856: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`
3857: @*/
3858: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A, PetscInt col, Vec *v)
3859: {
3860: PetscFunctionBegin;
3864: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3865: 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);
3866: PetscUseMethod(A, "MatDenseRestoreColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3867: PetscFunctionReturn(PETSC_SUCCESS);
3868: }
3870: /*@
3871: MatDenseGetSubMatrix - Gives access to a block of rows and columns of a dense matrix, represented as a `Mat`.
3873: Collective
3875: Input Parameters:
3876: + A - the `Mat` object
3877: . rbegin - the first global row index in the block (if `PETSC_DECIDE`, is 0)
3878: . rend - the global row index past the last one in the block (if `PETSC_DECIDE`, is `M`)
3879: . cbegin - the first global column index in the block (if `PETSC_DECIDE`, is 0)
3880: - cend - the global column index past the last one in the block (if `PETSC_DECIDE`, is `N`)
3882: Output Parameter:
3883: . v - the matrix
3885: Level: intermediate
3887: Notes:
3888: The matrix is owned by PETSc. Users need to call `MatDenseRestoreSubMatrix()` when the matrix is no longer needed.
3890: The output matrix is not redistributed by PETSc, so depending on the values of `rbegin` and `rend`, some processes may have no local rows.
3892: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreSubMatrix()`
3893: @*/
3894: PetscErrorCode MatDenseGetSubMatrix(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3895: {
3896: PetscFunctionBegin;
3903: PetscAssertPointer(v, 6);
3904: if (rbegin == PETSC_DECIDE) rbegin = 0;
3905: if (rend == PETSC_DECIDE) rend = A->rmap->N;
3906: if (cbegin == PETSC_DECIDE) cbegin = 0;
3907: if (cend == PETSC_DECIDE) cend = A->cmap->N;
3908: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3909: 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);
3910: 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);
3911: 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);
3912: 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);
3913: PetscUseMethod(A, "MatDenseGetSubMatrix_C", (Mat, PetscInt, PetscInt, PetscInt, PetscInt, Mat *), (A, rbegin, rend, cbegin, cend, v));
3914: PetscFunctionReturn(PETSC_SUCCESS);
3915: }
3917: /*@
3918: MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from `MatDenseGetSubMatrix()`.
3920: Collective
3922: Input Parameters:
3923: + A - the `Mat` object
3924: - v - the `Mat` object (may be `NULL`)
3926: Level: intermediate
3928: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseGetSubMatrix()`
3929: @*/
3930: PetscErrorCode MatDenseRestoreSubMatrix(Mat A, Mat *v)
3931: {
3932: PetscFunctionBegin;
3935: PetscAssertPointer(v, 2);
3936: PetscUseMethod(A, "MatDenseRestoreSubMatrix_C", (Mat, Mat *), (A, v));
3937: PetscFunctionReturn(PETSC_SUCCESS);
3938: }
3940: #include <petscblaslapack.h>
3941: #include <petsc/private/kernels/blockinvert.h>
3943: PetscErrorCode MatSeqDenseInvert(Mat A)
3944: {
3945: PetscInt m;
3946: const PetscReal shift = 0.0;
3947: PetscBool allowzeropivot, zeropivotdetected = PETSC_FALSE;
3948: PetscScalar *values;
3950: PetscFunctionBegin;
3952: PetscCall(MatDenseGetArray(A, &values));
3953: PetscCall(MatGetLocalSize(A, &m, NULL));
3954: allowzeropivot = PetscNot(A->erroriffailure);
3955: /* factor and invert each block */
3956: switch (m) {
3957: case 1:
3958: values[0] = (PetscScalar)1.0 / (values[0] + shift);
3959: break;
3960: case 2:
3961: PetscCall(PetscKernel_A_gets_inverse_A_2(values, shift, allowzeropivot, &zeropivotdetected));
3962: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3963: break;
3964: case 3:
3965: PetscCall(PetscKernel_A_gets_inverse_A_3(values, shift, allowzeropivot, &zeropivotdetected));
3966: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3967: break;
3968: case 4:
3969: PetscCall(PetscKernel_A_gets_inverse_A_4(values, shift, allowzeropivot, &zeropivotdetected));
3970: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3971: break;
3972: case 5: {
3973: PetscScalar work[25];
3974: PetscInt ipvt[5];
3976: PetscCall(PetscKernel_A_gets_inverse_A_5(values, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
3977: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3978: } break;
3979: case 6:
3980: PetscCall(PetscKernel_A_gets_inverse_A_6(values, shift, allowzeropivot, &zeropivotdetected));
3981: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3982: break;
3983: case 7:
3984: PetscCall(PetscKernel_A_gets_inverse_A_7(values, shift, allowzeropivot, &zeropivotdetected));
3985: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3986: break;
3987: default: {
3988: PetscInt *v_pivots, *IJ, j;
3989: PetscScalar *v_work;
3991: PetscCall(PetscMalloc3(m, &v_work, m, &v_pivots, m, &IJ));
3992: for (j = 0; j < m; j++) IJ[j] = j;
3993: PetscCall(PetscKernel_A_gets_inverse_A(m, values, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
3994: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3995: PetscCall(PetscFree3(v_work, v_pivots, IJ));
3996: }
3997: }
3998: PetscCall(MatDenseRestoreArray(A, &values));
3999: PetscFunctionReturn(PETSC_SUCCESS);
4000: }