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 PetscDefined(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 >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in LAPACK argument %" PetscBLASInt_FMT, -info);
77: PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad Inversion: zero pivot in row %" PetscBLASInt_FMT, info - 1);
78: PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
79: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix must be factored to solve");
81: A->ops->solve = NULL;
82: A->ops->matsolve = NULL;
83: A->ops->solvetranspose = NULL;
84: A->ops->matsolvetranspose = NULL;
85: A->ops->solveadd = NULL;
86: A->ops->solvetransposeadd = NULL;
87: A->factortype = MAT_FACTOR_NONE;
88: PetscCall(PetscFree(A->solvertype));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: static PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
93: {
94: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
95: PetscInt m = l->lda, n = A->cmap->n, r = A->rmap->n, i, j;
96: PetscScalar *slot, *bb, *v;
97: const PetscScalar *xx;
99: PetscFunctionBegin;
100: if (PetscDefined(USE_DEBUG)) {
101: for (i = 0; i < N; i++) {
102: PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
103: 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);
104: 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);
105: }
106: }
107: if (!N) PetscFunctionReturn(PETSC_SUCCESS);
109: /* fix right-hand side if needed */
110: if (x && b) {
111: Vec xt;
113: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
114: PetscCall(VecDuplicate(x, &xt));
115: PetscCall(VecCopy(x, xt));
116: PetscCall(VecScale(xt, -1.0));
117: PetscCall(MatMultAdd(A, xt, b, b));
118: PetscCall(VecDestroy(&xt));
119: PetscCall(VecGetArrayRead(x, &xx));
120: PetscCall(VecGetArray(b, &bb));
121: for (i = 0; i < N; i++) bb[rows[i]] = diag * xx[rows[i]];
122: PetscCall(VecRestoreArrayRead(x, &xx));
123: PetscCall(VecRestoreArray(b, &bb));
124: }
126: PetscCall(MatDenseGetArray(A, &v));
127: for (i = 0; i < N; i++) {
128: slot = v + rows[i] * m;
129: PetscCall(PetscArrayzero(slot, r));
130: }
131: for (i = 0; i < N; i++) {
132: slot = v + rows[i];
133: for (j = 0; j < n; j++) {
134: *slot = 0.0;
135: slot += m;
136: }
137: }
138: if (diag != 0.0) {
139: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
140: for (i = 0; i < N; i++) {
141: slot = v + (m + 1) * rows[i];
142: *slot = diag;
143: }
144: }
145: PetscCall(MatDenseRestoreArray(A, &v));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
150: {
151: Mat B = NULL;
152: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
153: Mat_SeqDense *b;
154: PetscInt *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->N, i;
155: const MatScalar *av;
156: PetscBool isseqdense;
158: PetscFunctionBegin;
159: if (reuse == MAT_REUSE_MATRIX) {
160: PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATSEQDENSE, &isseqdense));
161: PetscCheck(isseqdense, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_USER, "Cannot reuse matrix of type %s", ((PetscObject)*newmat)->type_name);
162: }
163: if (reuse != MAT_REUSE_MATRIX) {
164: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
165: PetscCall(MatSetSizes(B, m, n, m, n));
166: PetscCall(MatSetType(B, MATSEQDENSE));
167: PetscCall(MatSeqDenseSetPreallocation(B, NULL));
168: b = (Mat_SeqDense *)B->data;
169: } else {
170: b = (Mat_SeqDense *)(*newmat)->data;
171: for (i = 0; i < n; i++) PetscCall(PetscArrayzero(b->v + i * b->lda, m));
172: }
173: PetscCall(MatSeqAIJGetArrayRead(A, &av));
174: for (i = 0; i < m; i++) {
175: PetscInt j;
176: for (j = 0; j < ai[1] - ai[0]; j++) {
177: b->v[*aj * b->lda + i] = *av;
178: aj++;
179: av++;
180: }
181: ai++;
182: }
183: PetscCall(MatSeqAIJRestoreArrayRead(A, &av));
185: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
186: else if (B) *newmat = B;
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
191: {
192: Mat B = NULL;
193: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
194: PetscInt i, j;
195: PetscInt *rows, *nnz;
196: MatScalar *aa = a->v, *vals;
198: PetscFunctionBegin;
199: PetscCall(PetscCalloc3(A->rmap->n, &rows, A->rmap->n, &nnz, A->rmap->n, &vals));
200: if (reuse != MAT_REUSE_MATRIX) {
201: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
202: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
203: PetscCall(MatSetType(B, MATSEQAIJ));
204: for (j = 0; j < A->cmap->n; j++) {
205: for (i = 0; i < A->rmap->n; i++)
206: if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) ++nnz[i];
207: aa += a->lda;
208: }
209: PetscCall(MatSeqAIJSetPreallocation(B, PETSC_DETERMINE, nnz));
210: } else B = *newmat;
211: aa = a->v;
212: for (j = 0; j < A->cmap->n; j++) {
213: PetscInt numRows = 0;
214: for (i = 0; i < A->rmap->n; i++)
215: if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) {
216: rows[numRows] = i;
217: vals[numRows++] = aa[i];
218: }
219: PetscCall(MatSetValues(B, numRows, rows, 1, &j, vals, INSERT_VALUES));
220: aa += a->lda;
221: }
222: PetscCall(PetscFree3(rows, nnz, vals));
223: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
224: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
226: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
227: else if (reuse != MAT_REUSE_MATRIX) *newmat = B;
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: PetscErrorCode MatAXPY_SeqDense(Mat Y, PetscScalar alpha, Mat X, MatStructure str)
232: {
233: Mat_SeqDense *x = (Mat_SeqDense *)X->data, *y = (Mat_SeqDense *)Y->data;
234: const PetscScalar *xv;
235: PetscScalar *yv;
236: PetscBLASInt N, m, ldax = 0, lday = 0, one = 1;
238: PetscFunctionBegin;
239: PetscCall(MatDenseGetArrayRead(X, &xv));
240: PetscCall(MatDenseGetArray(Y, &yv));
241: PetscCall(PetscBLASIntCast(X->rmap->n * X->cmap->n, &N));
242: PetscCall(PetscBLASIntCast(X->rmap->n, &m));
243: PetscCall(PetscBLASIntCast(x->lda, &ldax));
244: PetscCall(PetscBLASIntCast(y->lda, &lday));
245: if (ldax > m || lday > m) {
246: for (PetscInt j = 0; j < X->cmap->n; j++) PetscCallBLAS("BLASaxpy", BLASaxpy_(&m, &alpha, PetscSafePointerPlusOffset(xv, j * ldax), &one, PetscSafePointerPlusOffset(yv, j * lday), &one));
247: } else {
248: PetscCallBLAS("BLASaxpy", BLASaxpy_(&N, &alpha, xv, &one, yv, &one));
249: }
250: PetscCall(MatDenseRestoreArrayRead(X, &xv));
251: PetscCall(MatDenseRestoreArray(Y, &yv));
252: PetscCall(PetscLogFlops(PetscMax(2.0 * N - 1, 0)));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: static PetscErrorCode MatGetInfo_SeqDense(Mat A, MatInfoType flag, MatInfo *info)
257: {
258: PetscLogDouble N = A->rmap->n * A->cmap->n;
260: PetscFunctionBegin;
261: info->block_size = 1.0;
262: info->nz_allocated = N;
263: info->nz_used = N;
264: info->nz_unneeded = 0;
265: info->assemblies = A->num_ass;
266: info->mallocs = 0;
267: info->memory = 0; /* REVIEW ME */
268: info->fill_ratio_given = 0;
269: info->fill_ratio_needed = 0;
270: info->factor_mallocs = 0;
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: PetscErrorCode MatScale_SeqDense(Mat A, PetscScalar alpha)
275: {
276: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
277: PetscScalar *v;
278: PetscBLASInt one = 1, j, nz, lda = 0;
280: PetscFunctionBegin;
281: PetscCall(MatDenseGetArray(A, &v));
282: PetscCall(PetscBLASIntCast(a->lda, &lda));
283: if (lda > A->rmap->n) {
284: PetscCall(PetscBLASIntCast(A->rmap->n, &nz));
285: for (j = 0; j < A->cmap->n; j++) PetscCallBLAS("BLASscal", BLASscal_(&nz, &alpha, v + j * lda, &one));
286: } else {
287: PetscCall(PetscBLASIntCast(A->rmap->n * A->cmap->n, &nz));
288: PetscCallBLAS("BLASscal", BLASscal_(&nz, &alpha, v, &one));
289: }
290: PetscCall(PetscLogFlops(A->rmap->n * A->cmap->n));
291: PetscCall(MatDenseRestoreArray(A, &v));
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: PetscErrorCode MatShift_SeqDense(Mat A, PetscScalar alpha)
296: {
297: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
298: PetscScalar *v;
299: PetscInt j, k;
301: PetscFunctionBegin;
302: PetscCall(MatDenseGetArray(A, &v));
303: k = PetscMin(A->rmap->n, A->cmap->n);
304: for (j = 0; j < k; j++) v[j + j * a->lda] += alpha;
305: PetscCall(PetscLogFlops(k));
306: PetscCall(MatDenseRestoreArray(A, &v));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: static PetscErrorCode MatIsHermitian_SeqDense(Mat A, PetscReal rtol, PetscBool *fl)
311: {
312: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
313: PetscInt i, j, m = A->rmap->n, N = a->lda;
314: const PetscScalar *v;
316: PetscFunctionBegin;
317: *fl = PETSC_FALSE;
318: if (A->rmap->n != A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
319: PetscCall(MatDenseGetArrayRead(A, &v));
320: for (i = 0; i < m; i++) {
321: for (j = i; j < m; j++) {
322: if (PetscAbsScalar(v[i + j * N] - PetscConj(v[j + i * N])) > rtol) goto restore;
323: }
324: }
325: *fl = PETSC_TRUE;
326: restore:
327: PetscCall(MatDenseRestoreArrayRead(A, &v));
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: static PetscErrorCode MatIsSymmetric_SeqDense(Mat A, PetscReal rtol, PetscBool *fl)
332: {
333: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
334: PetscInt i, j, m = A->rmap->n, N = a->lda;
335: const PetscScalar *v;
337: PetscFunctionBegin;
338: *fl = PETSC_FALSE;
339: if (A->rmap->n != A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
340: PetscCall(MatDenseGetArrayRead(A, &v));
341: for (i = 0; i < m; i++) {
342: for (j = i; j < m; j++) {
343: if (PetscAbsScalar(v[i + j * N] - v[j + i * N]) > rtol) goto restore;
344: }
345: }
346: *fl = PETSC_TRUE;
347: restore:
348: PetscCall(MatDenseRestoreArrayRead(A, &v));
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi, Mat A, MatDuplicateOption cpvalues)
353: {
354: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
355: PetscInt lda = mat->lda, j, m, nlda = lda;
356: PetscBool isdensecpu;
358: PetscFunctionBegin;
359: PetscCall(PetscLayoutReference(A->rmap, &newi->rmap));
360: PetscCall(PetscLayoutReference(A->cmap, &newi->cmap));
361: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */
362: PetscCall(MatDenseSetLDA(newi, lda));
363: }
364: PetscCall(PetscObjectTypeCompare((PetscObject)newi, MATSEQDENSE, &isdensecpu));
365: if (isdensecpu) PetscCall(MatSeqDenseSetPreallocation(newi, NULL));
366: if (cpvalues == MAT_COPY_VALUES) {
367: const PetscScalar *av;
368: PetscScalar *v;
370: PetscCall(MatDenseGetArrayRead(A, &av));
371: PetscCall(MatDenseGetArrayWrite(newi, &v));
372: PetscCall(MatDenseGetLDA(newi, &nlda));
373: m = A->rmap->n;
374: if (lda > m || nlda > m) {
375: for (j = 0; j < A->cmap->n; j++) PetscCall(PetscArraycpy(PetscSafePointerPlusOffset(v, j * nlda), PetscSafePointerPlusOffset(av, j * lda), m));
376: } else {
377: PetscCall(PetscArraycpy(v, av, A->rmap->n * A->cmap->n));
378: }
379: PetscCall(MatDenseRestoreArrayWrite(newi, &v));
380: PetscCall(MatDenseRestoreArrayRead(A, &av));
381: PetscCall(MatPropagateSymmetryOptions(A, newi));
382: }
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: PetscErrorCode MatDuplicate_SeqDense(Mat A, MatDuplicateOption cpvalues, Mat *newmat)
387: {
388: PetscFunctionBegin;
389: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), newmat));
390: PetscCall(MatSetSizes(*newmat, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
391: PetscCall(MatSetType(*newmat, ((PetscObject)A)->type_name));
392: PetscCall(MatDuplicateNoCreate_SeqDense(*newmat, A, cpvalues));
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: static PetscErrorCode MatSolve_SeqDense_Internal_LU(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
397: {
398: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
400: PetscFunctionBegin;
401: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
402: PetscCallLAPACKInfo("LAPACKgetrs", LAPACKgetrs_(T ? "T" : "N", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
403: PetscCall(PetscFPTrapPop());
404: PetscCall(PetscLogFlops(nrhs * (2.0 * m * m - m)));
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: static PetscErrorCode MatSolve_SeqDense_Internal_Cholesky(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
409: {
410: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
412: PetscFunctionBegin;
413: if (A->spd == PETSC_BOOL3_TRUE) {
414: if (PetscDefined(USE_COMPLEX) && T) PetscCall(MatConjugate_SeqDense(A));
415: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
416: PetscCallLAPACKInfo("LAPACKpotrs", LAPACKpotrs_("L", &m, &nrhs, mat->v, &mat->lda, x, &m, &info));
417: PetscCall(PetscFPTrapPop());
418: if (PetscDefined(USE_COMPLEX) && T) PetscCall(MatConjugate_SeqDense(A));
419: #if PetscDefined(USE_COMPLEX)
420: } else if (A->hermitian == PETSC_BOOL3_TRUE) {
421: if (T) PetscCall(MatConjugate_SeqDense(A));
422: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
423: PetscCallLAPACKInfo("LAPACKhetrs", LAPACKhetrs_("L", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
424: PetscCall(PetscFPTrapPop());
425: if (T) PetscCall(MatConjugate_SeqDense(A));
426: #endif
427: } else { /* symmetric case */
428: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
429: PetscCallLAPACKInfo("LAPACKsytrs", LAPACKsytrs_("L", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
430: PetscCall(PetscFPTrapPop());
431: }
432: PetscCall(PetscLogFlops(nrhs * (2.0 * m * m - m)));
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: static PetscErrorCode MatSolve_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
437: {
438: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
439: char trans;
441: PetscFunctionBegin;
442: if (PetscDefined(USE_COMPLEX)) {
443: trans = 'C';
444: } else {
445: trans = 'T';
446: }
447: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
448: { /* lwork depends on the number of right-hand sides */
449: PetscBLASInt nlfwork, lfwork = -1;
450: PetscScalar fwork;
452: PetscCallLAPACKInfo("LAPACKormqr", LAPACKormqr_("L", &trans, &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, &fwork, &lfwork, &info));
453: nlfwork = (PetscBLASInt)PetscRealPart(fwork);
454: if (nlfwork > mat->lfwork) {
455: mat->lfwork = nlfwork;
456: PetscCall(PetscFree(mat->fwork));
457: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
458: }
459: }
460: PetscCallLAPACKInfo("LAPACKormqr", LAPACKormqr_("L", &trans, &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, mat->fwork, &mat->lfwork, &info));
461: PetscCall(PetscFPTrapPop());
462: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
463: PetscCallLAPACKInfo("LAPACKtrtrs", LAPACKtrtrs_("U", "N", "N", &mat->rank, &nrhs, mat->v, &mat->lda, x, &ldx, &info));
464: PetscCall(PetscFPTrapPop());
465: for (PetscInt j = 0; j < nrhs; j++) {
466: for (PetscInt i = mat->rank; i < k; i++) x[j * ldx + i] = 0.;
467: }
468: PetscCall(PetscLogFlops(nrhs * (4.0 * m * mat->rank - PetscSqr(mat->rank))));
469: PetscFunctionReturn(PETSC_SUCCESS);
470: }
472: static PetscErrorCode MatSolveTranspose_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
473: {
474: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
476: PetscFunctionBegin;
477: if (A->rmap->n == A->cmap->n && mat->rank == A->rmap->n) {
478: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
479: PetscCallLAPACKInfo("LAPACKtrtrs", LAPACKtrtrs_("U", "T", "N", &m, &nrhs, mat->v, &mat->lda, x, &ldx, &info));
480: PetscCall(PetscFPTrapPop());
481: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate_SeqDense(A));
482: { /* lwork depends on the number of right-hand sides */
483: PetscBLASInt nlfwork, lfwork = -1;
484: PetscScalar fwork;
486: PetscCallLAPACKInfo("LAPACKormqr", LAPACKormqr_("L", "N", &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, &fwork, &lfwork, &info));
487: nlfwork = (PetscBLASInt)PetscRealPart(fwork);
488: if (nlfwork > mat->lfwork) {
489: mat->lfwork = nlfwork;
490: PetscCall(PetscFree(mat->fwork));
491: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
492: }
493: }
494: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
495: PetscCallLAPACKInfo("LAPACKormqr", LAPACKormqr_("L", "N", &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, mat->fwork, &mat->lfwork, &info));
496: PetscCall(PetscFPTrapPop());
497: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate_SeqDense(A));
498: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "QR factored matrix cannot be used for transpose solve");
499: PetscCall(PetscLogFlops(nrhs * (4.0 * m * mat->rank - PetscSqr(mat->rank))));
500: PetscFunctionReturn(PETSC_SUCCESS);
501: }
503: static PetscErrorCode MatSolve_SeqDense_SetUp(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
504: {
505: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
506: PetscScalar *y;
507: PetscBLASInt m = 0, k = 0;
509: PetscFunctionBegin;
510: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
511: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
512: if (k < m) {
513: PetscCall(VecCopy(xx, mat->qrrhs));
514: PetscCall(VecGetArray(mat->qrrhs, &y));
515: } else {
516: PetscCall(VecCopy(xx, yy));
517: PetscCall(VecGetArray(yy, &y));
518: }
519: *_y = y;
520: *_k = k;
521: *_m = m;
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: static PetscErrorCode MatSolve_SeqDense_TearDown(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
526: {
527: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
528: PetscScalar *y = NULL;
529: PetscBLASInt m, k;
531: PetscFunctionBegin;
532: y = *_y;
533: *_y = NULL;
534: k = *_k;
535: m = *_m;
536: if (k < m) {
537: PetscScalar *yv;
538: PetscCall(VecGetArray(yy, &yv));
539: PetscCall(PetscArraycpy(yv, y, k));
540: PetscCall(VecRestoreArray(yy, &yv));
541: PetscCall(VecRestoreArray(mat->qrrhs, &y));
542: } else {
543: PetscCall(VecRestoreArray(yy, &y));
544: }
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: static PetscErrorCode MatSolve_SeqDense_LU(Mat A, Vec xx, Vec yy)
549: {
550: PetscScalar *y = NULL;
551: PetscBLASInt m = 0, k = 0;
553: PetscFunctionBegin;
554: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
555: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_FALSE));
556: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: static PetscErrorCode MatSolveTranspose_SeqDense_LU(Mat A, Vec xx, Vec yy)
561: {
562: PetscScalar *y = NULL;
563: PetscBLASInt m = 0, k = 0;
565: PetscFunctionBegin;
566: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
567: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_TRUE));
568: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: static PetscErrorCode MatSolve_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
573: {
574: PetscScalar *y = NULL;
575: PetscBLASInt m = 0, k = 0;
577: PetscFunctionBegin;
578: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
579: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_FALSE));
580: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: static PetscErrorCode MatSolveTranspose_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
585: {
586: PetscScalar *y = NULL;
587: PetscBLASInt m = 0, k = 0;
589: PetscFunctionBegin;
590: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
591: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_TRUE));
592: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: static PetscErrorCode MatSolve_SeqDense_QR(Mat A, Vec xx, Vec yy)
597: {
598: PetscScalar *y = NULL;
599: PetscBLASInt m = 0, k = 0;
601: PetscFunctionBegin;
602: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
603: PetscCall(MatSolve_SeqDense_Internal_QR(A, y, PetscMax(m, k), m, 1, k));
604: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: static PetscErrorCode MatSolveTranspose_SeqDense_QR(Mat A, Vec xx, Vec yy)
609: {
610: PetscScalar *y = NULL;
611: PetscBLASInt m = 0, k = 0;
613: PetscFunctionBegin;
614: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
615: PetscCall(MatSolveTranspose_SeqDense_Internal_QR(A, y, PetscMax(m, k), m, 1, k));
616: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
620: static PetscErrorCode MatMatSolve_SeqDense_SetUp(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
621: {
622: const PetscScalar *b;
623: PetscScalar *y;
624: PetscInt n, _ldb, _ldx;
625: PetscBLASInt nrhs = 0, m = 0, k = 0, ldb = 0, ldx = 0, ldy = 0;
627: PetscFunctionBegin;
628: *_ldy = 0;
629: *_m = 0;
630: *_nrhs = 0;
631: *_k = 0;
632: *_y = NULL;
633: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
634: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
635: PetscCall(MatGetSize(B, NULL, &n));
636: PetscCall(PetscBLASIntCast(n, &nrhs));
637: PetscCall(MatDenseGetLDA(B, &_ldb));
638: PetscCall(PetscBLASIntCast(_ldb, &ldb));
639: PetscCall(MatDenseGetLDA(X, &_ldx));
640: PetscCall(PetscBLASIntCast(_ldx, &ldx));
641: if (ldx < m) {
642: PetscCall(MatDenseGetArrayRead(B, &b));
643: PetscCall(PetscMalloc1(nrhs * m, &y));
644: if (ldb == m) {
645: PetscCall(PetscArraycpy(y, b, ldb * nrhs));
646: } else {
647: for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&y[j * m], &b[j * ldb], m));
648: }
649: ldy = m;
650: PetscCall(MatDenseRestoreArrayRead(B, &b));
651: } else {
652: if (ldb == ldx) {
653: PetscCall(MatCopy(B, X, SAME_NONZERO_PATTERN));
654: PetscCall(MatDenseGetArray(X, &y));
655: } else {
656: PetscCall(MatDenseGetArray(X, &y));
657: PetscCall(MatDenseGetArrayRead(B, &b));
658: for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&y[j * ldx], &b[j * ldb], m));
659: PetscCall(MatDenseRestoreArrayRead(B, &b));
660: }
661: ldy = ldx;
662: }
663: *_y = y;
664: *_ldy = ldy;
665: *_k = k;
666: *_m = m;
667: *_nrhs = nrhs;
668: PetscFunctionReturn(PETSC_SUCCESS);
669: }
671: static PetscErrorCode MatMatSolve_SeqDense_TearDown(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
672: {
673: PetscScalar *y;
674: PetscInt _ldx;
675: PetscBLASInt k, ldy, nrhs, ldx = 0;
677: PetscFunctionBegin;
678: y = *_y;
679: *_y = NULL;
680: k = *_k;
681: ldy = *_ldy;
682: nrhs = *_nrhs;
683: PetscCall(MatDenseGetLDA(X, &_ldx));
684: PetscCall(PetscBLASIntCast(_ldx, &ldx));
685: if (ldx != ldy) {
686: PetscScalar *xv;
687: PetscCall(MatDenseGetArray(X, &xv));
688: for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&xv[j * ldx], &y[j * ldy], k));
689: PetscCall(MatDenseRestoreArray(X, &xv));
690: PetscCall(PetscFree(y));
691: } else {
692: PetscCall(MatDenseRestoreArray(X, &y));
693: }
694: PetscFunctionReturn(PETSC_SUCCESS);
695: }
697: static PetscErrorCode MatMatSolve_SeqDense_LU(Mat A, Mat B, Mat X)
698: {
699: PetscScalar *y;
700: PetscBLASInt m, k, ldy, nrhs;
702: PetscFunctionBegin;
703: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
704: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_FALSE));
705: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: static PetscErrorCode MatMatSolveTranspose_SeqDense_LU(Mat A, Mat B, Mat X)
710: {
711: PetscScalar *y;
712: PetscBLASInt m, k, ldy, nrhs;
714: PetscFunctionBegin;
715: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
716: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_TRUE));
717: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: static PetscErrorCode MatMatSolve_SeqDense_Cholesky(Mat A, Mat B, Mat X)
722: {
723: PetscScalar *y;
724: PetscBLASInt m, k, ldy, nrhs;
726: PetscFunctionBegin;
727: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
728: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_FALSE));
729: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
730: PetscFunctionReturn(PETSC_SUCCESS);
731: }
733: static PetscErrorCode MatMatSolveTranspose_SeqDense_Cholesky(Mat A, Mat B, Mat X)
734: {
735: PetscScalar *y;
736: PetscBLASInt m, k, ldy, nrhs;
738: PetscFunctionBegin;
739: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
740: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_TRUE));
741: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }
745: static PetscErrorCode MatMatSolve_SeqDense_QR(Mat A, Mat B, Mat X)
746: {
747: PetscScalar *y;
748: PetscBLASInt m, k, ldy, nrhs;
750: PetscFunctionBegin;
751: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
752: PetscCall(MatSolve_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k));
753: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
754: PetscFunctionReturn(PETSC_SUCCESS);
755: }
757: static PetscErrorCode MatMatSolveTranspose_SeqDense_QR(Mat A, Mat B, Mat X)
758: {
759: PetscScalar *y;
760: PetscBLASInt m, k, ldy, nrhs;
762: PetscFunctionBegin;
763: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
764: PetscCall(MatSolveTranspose_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k));
765: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
766: PetscFunctionReturn(PETSC_SUCCESS);
767: }
769: /* COMMENT: I have chosen to hide row permutation in the pivots,
770: rather than put it in the Mat->row slot.*/
771: PetscErrorCode MatLUFactor_SeqDense(Mat A, IS row, IS col, PETSC_UNUSED const MatFactorInfo *minfo)
772: {
773: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
774: PetscBLASInt n, m, info;
776: PetscFunctionBegin;
777: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
778: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
779: if (!mat->pivots) PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots));
780: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
781: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
782: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&m, &n, mat->v, &mat->lda, mat->pivots, &info));
783: PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in LAPACK argument %" PetscBLASInt_FMT, -info);
784: PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Bad factorization: zero pivot in row %" PetscBLASInt_FMT, info - 1);
785: PetscCall(PetscFPTrapPop());
787: A->ops->solve = MatSolve_SeqDense_LU;
788: A->ops->matsolve = MatMatSolve_SeqDense_LU;
789: A->ops->solvetranspose = MatSolveTranspose_SeqDense_LU;
790: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_LU;
791: A->factortype = MAT_FACTOR_LU;
793: PetscCall(PetscFree(A->solvertype));
794: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
796: PetscCall(PetscLogFlops((2.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3));
797: PetscFunctionReturn(PETSC_SUCCESS);
798: }
800: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info)
801: {
802: PetscFunctionBegin;
803: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
804: PetscUseTypeMethod(fact, lufactor, NULL, NULL, info);
805: PetscFunctionReturn(PETSC_SUCCESS);
806: }
808: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, IS col, PETSC_UNUSED const MatFactorInfo *info)
809: {
810: PetscFunctionBegin;
811: fact->preallocated = PETSC_TRUE;
812: fact->assembled = PETSC_TRUE;
813: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
814: PetscFunctionReturn(PETSC_SUCCESS);
815: }
817: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
818: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A, IS perm, PETSC_UNUSED const MatFactorInfo *minfo)
819: {
820: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
821: PetscBLASInt info, n;
823: PetscFunctionBegin;
824: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
825: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
826: if (A->spd == PETSC_BOOL3_TRUE) {
827: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
828: PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("L", &n, mat->v, &mat->lda, &info));
829: PetscCall(PetscFPTrapPop());
830: #if PetscDefined(USE_COMPLEX)
831: } else if (A->hermitian == PETSC_BOOL3_TRUE) {
832: if (!mat->pivots) PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots));
833: if (!mat->fwork) {
834: PetscScalar dummy;
836: mat->lfwork = -1;
837: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
838: PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info));
839: PetscCall(PetscFPTrapPop());
840: PetscCall(PetscBLASIntCast((PetscCount)(PetscRealPart(dummy)), &mat->lfwork));
841: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
842: }
843: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
844: PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
845: PetscCall(PetscFPTrapPop());
846: #endif
847: } else { /* symmetric case */
848: if (!mat->pivots) PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots));
849: if (!mat->fwork) {
850: PetscScalar dummy;
852: mat->lfwork = -1;
853: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
854: PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info));
855: PetscCall(PetscFPTrapPop());
856: PetscCall(PetscBLASIntCast((PetscCount)(PetscRealPart(dummy)), &mat->lfwork));
857: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
858: }
859: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
860: PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
861: PetscCall(PetscFPTrapPop());
862: }
863: PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in LAPACK argument %" PetscBLASInt_FMT, -info);
864: PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad factorization: zero pivot in row %" PetscBLASInt_FMT, info - 1);
866: A->ops->solve = MatSolve_SeqDense_Cholesky;
867: A->ops->matsolve = MatMatSolve_SeqDense_Cholesky;
868: A->ops->solvetranspose = MatSolveTranspose_SeqDense_Cholesky;
869: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky;
870: A->factortype = MAT_FACTOR_CHOLESKY;
872: PetscCall(PetscFree(A->solvertype));
873: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
875: PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
876: PetscFunctionReturn(PETSC_SUCCESS);
877: }
879: static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info)
880: {
881: PetscFunctionBegin;
882: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
883: PetscUseTypeMethod(fact, choleskyfactor, NULL, info);
884: PetscFunctionReturn(PETSC_SUCCESS);
885: }
887: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info)
888: {
889: PetscFunctionBegin;
890: fact->assembled = PETSC_TRUE;
891: fact->preallocated = PETSC_TRUE;
892: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
893: PetscFunctionReturn(PETSC_SUCCESS);
894: }
896: PetscErrorCode MatQRFactor_SeqDense(Mat A, IS col, PETSC_UNUSED const MatFactorInfo *minfo)
897: {
898: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
899: PetscBLASInt n, m, min, max;
901: PetscFunctionBegin;
902: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
903: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
904: max = PetscMax(m, n);
905: min = PetscMin(m, n);
906: if (!mat->tau) PetscCall(PetscMalloc1(min, &mat->tau));
907: if (!mat->pivots) PetscCall(PetscMalloc1(n, &mat->pivots));
908: if (!mat->qrrhs) PetscCall(MatCreateVecs(A, NULL, &mat->qrrhs));
909: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
910: if (!mat->fwork) {
911: PetscScalar dummy;
913: mat->lfwork = -1;
914: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
915: PetscCallLAPACKInfo("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, &dummy, &mat->lfwork, &info));
916: PetscCall(PetscFPTrapPop());
917: PetscCall(PetscBLASIntCast((PetscCount)(PetscRealPart(dummy)), &mat->lfwork));
918: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
919: }
920: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
921: PetscCallLAPACKInfo("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, mat->fwork, &mat->lfwork, &info));
922: PetscCall(PetscFPTrapPop());
923: // 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
924: mat->rank = min;
926: A->ops->solve = MatSolve_SeqDense_QR;
927: A->ops->matsolve = MatMatSolve_SeqDense_QR;
928: A->factortype = MAT_FACTOR_QR;
929: if (m == n) {
930: A->ops->solvetranspose = MatSolveTranspose_SeqDense_QR;
931: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
932: }
934: PetscCall(PetscFree(A->solvertype));
935: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
937: PetscCall(PetscLogFlops(2.0 * min * min * (max - min / 3.0)));
938: PetscFunctionReturn(PETSC_SUCCESS);
939: }
941: static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info)
942: {
943: PetscFunctionBegin;
944: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
945: PetscUseMethod(fact, "MatQRFactor_C", (Mat, IS, const MatFactorInfo *), (fact, NULL, info));
946: PetscFunctionReturn(PETSC_SUCCESS);
947: }
949: PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info)
950: {
951: PetscFunctionBegin;
952: fact->assembled = PETSC_TRUE;
953: fact->preallocated = PETSC_TRUE;
954: PetscCall(PetscObjectComposeFunction((PetscObject)fact, "MatQRFactorNumeric_C", MatQRFactorNumeric_SeqDense));
955: PetscFunctionReturn(PETSC_SUCCESS);
956: }
958: /* uses LAPACK */
959: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A, MatFactorType ftype, Mat *fact)
960: {
961: PetscFunctionBegin;
962: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), fact));
963: PetscCall(MatSetSizes(*fact, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
964: PetscCall(MatSetType(*fact, MATDENSE));
965: (*fact)->trivialsymbolic = PETSC_TRUE;
966: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
967: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
968: (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
969: } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
970: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
971: } else if (ftype == MAT_FACTOR_QR) {
972: PetscCall(PetscObjectComposeFunction((PetscObject)*fact, "MatQRFactorSymbolic_C", MatQRFactorSymbolic_SeqDense));
973: }
974: (*fact)->factortype = ftype;
976: PetscCall(PetscFree((*fact)->solvertype));
977: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*fact)->solvertype));
978: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_LU]));
979: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ILU]));
980: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_CHOLESKY]));
981: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ICC]));
982: PetscFunctionReturn(PETSC_SUCCESS);
983: }
985: static PetscErrorCode MatSOR_SeqDense(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal shift, PetscInt its, PetscInt lits, Vec xx)
986: {
987: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
988: PetscScalar *x, *v = mat->v, zero = 0.0, xt;
989: const PetscScalar *b;
990: PetscInt m = A->rmap->n, i;
991: PetscBLASInt o = 1, bm = 0;
993: PetscFunctionBegin;
994: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
995: PetscCheck(A->offloadmask != PETSC_OFFLOAD_GPU, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented");
996: #endif
997: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
998: PetscCall(PetscBLASIntCast(m, &bm));
999: if (flag & SOR_ZERO_INITIAL_GUESS) {
1000: /* this is a hack fix, should have another version without the second BLASdotu */
1001: PetscCall(VecSet(xx, zero));
1002: }
1003: PetscCall(VecGetArray(xx, &x));
1004: PetscCall(VecGetArrayRead(bb, &b));
1005: its = its * lits;
1006: 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);
1007: while (its--) {
1008: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1009: for (i = 0; i < m; i++) {
1010: PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o));
1011: x[i] = (1. - omega) * x[i] + (xt + v[i + i * m] * x[i]) * omega / (v[i + i * m] + shift);
1012: }
1013: }
1014: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1015: for (i = m - 1; i >= 0; i--) {
1016: PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o));
1017: x[i] = (1. - omega) * x[i] + (xt + v[i + i * m] * x[i]) * omega / (v[i + i * m] + shift);
1018: }
1019: }
1020: }
1021: PetscCall(VecRestoreArrayRead(bb, &b));
1022: PetscCall(VecRestoreArray(xx, &x));
1023: PetscFunctionReturn(PETSC_SUCCESS);
1024: }
1026: PETSC_INTERN PetscErrorCode MatMultColumnRangeKernel_SeqDense(Mat A, Vec xx, Vec yy, PetscInt c_start, PetscInt c_end, PetscBool trans, PetscBool herm)
1027: {
1028: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1029: PetscScalar *y, _DOne = 1.0, _DZero = 0.0;
1030: PetscBLASInt m, n, _One = 1;
1031: const PetscScalar *v = mat->v, *x;
1033: PetscFunctionBegin;
1034: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1035: PetscCall(PetscBLASIntCast(c_end - c_start, &n));
1036: PetscCall(VecGetArrayRead(xx, &x));
1037: PetscCall(VecGetArrayWrite(yy, &y));
1038: if (!m || !n) {
1039: PetscBLASInt i;
1040: if (trans)
1041: for (i = 0; i < n; i++) y[i] = 0.0;
1042: else
1043: for (i = 0; i < m; i++) y[i] = 0.0;
1044: } else {
1045: if (trans) {
1046: if (herm) PetscCallBLAS("BLASgemv", BLASgemv_("C", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DZero, y + c_start, &_One));
1047: else PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DZero, y + c_start, &_One));
1048: } else {
1049: PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x + c_start, &_One, &_DZero, y, &_One));
1050: }
1051: PetscCall(PetscLogFlops(2.0 * m * n - n));
1052: }
1053: PetscCall(VecRestoreArrayRead(xx, &x));
1054: PetscCall(VecRestoreArrayWrite(yy, &y));
1055: PetscFunctionReturn(PETSC_SUCCESS);
1056: }
1058: PetscErrorCode MatMultHermitianTransposeColumnRange_SeqDense(Mat A, Vec xx, Vec yy, PetscInt c_start, PetscInt c_end)
1059: {
1060: PetscFunctionBegin;
1061: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, c_start, c_end, PETSC_TRUE, PETSC_TRUE));
1062: PetscFunctionReturn(PETSC_SUCCESS);
1063: }
1065: PetscErrorCode MatMult_SeqDense(Mat A, Vec xx, Vec yy)
1066: {
1067: PetscFunctionBegin;
1068: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_FALSE, PETSC_FALSE));
1069: PetscFunctionReturn(PETSC_SUCCESS);
1070: }
1072: PetscErrorCode MatMultTranspose_SeqDense(Mat A, Vec xx, Vec yy)
1073: {
1074: PetscFunctionBegin;
1075: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_FALSE));
1076: PetscFunctionReturn(PETSC_SUCCESS);
1077: }
1079: PetscErrorCode MatMultHermitianTranspose_SeqDense(Mat A, Vec xx, Vec yy)
1080: {
1081: PetscFunctionBegin;
1082: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_TRUE));
1083: PetscFunctionReturn(PETSC_SUCCESS);
1084: }
1086: PETSC_INTERN PetscErrorCode MatMultAddColumnRangeKernel_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end, PetscBool trans, PetscBool herm)
1087: {
1088: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1089: const PetscScalar *v = mat->v, *x;
1090: PetscScalar *y, _DOne = 1.0;
1091: PetscBLASInt m, n, _One = 1;
1093: PetscFunctionBegin;
1094: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1095: PetscCall(PetscBLASIntCast(c_end - c_start, &n));
1096: PetscCall(VecCopy(zz, yy));
1097: if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS);
1098: PetscCall(VecGetArray(yy, &y));
1099: PetscCall(VecGetArrayRead(xx, &x));
1100: if (trans) {
1101: if (herm) PetscCallBLAS("BLASgemv", BLASgemv_("C", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DOne, y + c_start, &_One));
1102: else PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DOne, y + c_start, &_One));
1103: } else {
1104: PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x + c_start, &_One, &_DOne, y, &_One));
1105: }
1106: PetscCall(VecRestoreArrayRead(xx, &x));
1107: PetscCall(VecRestoreArray(yy, &y));
1108: PetscCall(PetscLogFlops(2.0 * m * n));
1109: PetscFunctionReturn(PETSC_SUCCESS);
1110: }
1112: PetscErrorCode MatMultColumnRange_SeqDense(Mat A, Vec xx, Vec yy, PetscInt c_start, PetscInt c_end)
1113: {
1114: PetscFunctionBegin;
1115: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, c_start, c_end, PETSC_FALSE, PETSC_FALSE));
1116: PetscFunctionReturn(PETSC_SUCCESS);
1117: }
1119: PetscErrorCode MatMultAddColumnRange_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end)
1120: {
1121: PetscFunctionBegin;
1122: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, c_start, c_end, PETSC_FALSE, PETSC_FALSE));
1123: PetscFunctionReturn(PETSC_SUCCESS);
1124: }
1126: PetscErrorCode MatMultHermitianTransposeAddColumnRange_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end)
1127: {
1128: PetscFunctionBegin;
1129: PetscMPIInt rank;
1130: PetscCallMPI(MPI_Comm_rank(MPI_COMM_WORLD, &rank));
1131: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, c_start, c_end, PETSC_TRUE, PETSC_TRUE));
1132: PetscFunctionReturn(PETSC_SUCCESS);
1133: }
1135: PetscErrorCode MatMultAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1136: {
1137: PetscFunctionBegin;
1138: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_FALSE, PETSC_FALSE));
1139: PetscFunctionReturn(PETSC_SUCCESS);
1140: }
1142: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1143: {
1144: PetscFunctionBegin;
1145: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_FALSE));
1146: PetscFunctionReturn(PETSC_SUCCESS);
1147: }
1149: PetscErrorCode MatMultHermitianTransposeAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1150: {
1151: PetscFunctionBegin;
1152: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_TRUE));
1153: PetscFunctionReturn(PETSC_SUCCESS);
1154: }
1156: static PetscErrorCode MatGetRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1157: {
1158: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1159: PetscInt i;
1161: PetscFunctionBegin;
1162: if (ncols) *ncols = A->cmap->n;
1163: if (cols) {
1164: PetscCall(PetscMalloc1(A->cmap->n, cols));
1165: for (i = 0; i < A->cmap->n; i++) (*cols)[i] = i;
1166: }
1167: if (vals) {
1168: const PetscScalar *v;
1170: PetscCall(MatDenseGetArrayRead(A, &v));
1171: PetscCall(PetscMalloc1(A->cmap->n, vals));
1172: v += row;
1173: for (i = 0; i < A->cmap->n; i++) {
1174: (*vals)[i] = *v;
1175: v += mat->lda;
1176: }
1177: PetscCall(MatDenseRestoreArrayRead(A, &v));
1178: }
1179: PetscFunctionReturn(PETSC_SUCCESS);
1180: }
1182: static PetscErrorCode MatRestoreRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1183: {
1184: PetscFunctionBegin;
1185: if (cols) PetscCall(PetscFree(*cols));
1186: if (vals) PetscCall(PetscFree(*vals));
1187: PetscFunctionReturn(PETSC_SUCCESS);
1188: }
1190: static PetscErrorCode MatSetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], const PetscScalar v[], InsertMode addv)
1191: {
1192: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1193: PetscScalar *av;
1194: PetscInt i, j, idx = 0;
1195: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
1196: PetscOffloadMask oldf;
1197: #endif
1199: PetscFunctionBegin;
1200: PetscCall(MatDenseGetArray(A, &av));
1201: if (!mat->roworiented) {
1202: if (addv == INSERT_VALUES) {
1203: for (j = 0; j < n; j++) {
1204: if (indexn[j] < 0) {
1205: idx += m;
1206: continue;
1207: }
1208: 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);
1209: for (i = 0; i < m; i++) {
1210: if (indexm[i] < 0) {
1211: idx++;
1212: continue;
1213: }
1214: 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);
1215: av[indexn[j] * mat->lda + indexm[i]] = v ? v[idx++] : (idx++, 0.0);
1216: }
1217: }
1218: } else {
1219: for (j = 0; j < n; j++) {
1220: if (indexn[j] < 0) {
1221: idx += m;
1222: continue;
1223: }
1224: 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);
1225: for (i = 0; i < m; i++) {
1226: if (indexm[i] < 0) {
1227: idx++;
1228: continue;
1229: }
1230: 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);
1231: av[indexn[j] * mat->lda + indexm[i]] += v ? v[idx++] : (idx++, 0.0);
1232: }
1233: }
1234: }
1235: } else {
1236: if (addv == INSERT_VALUES) {
1237: for (i = 0; i < m; i++) {
1238: if (indexm[i] < 0) {
1239: idx += n;
1240: continue;
1241: }
1242: 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);
1243: for (j = 0; j < n; j++) {
1244: if (indexn[j] < 0) {
1245: idx++;
1246: continue;
1247: }
1248: PetscCheck(indexn[j] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, indexn[j], A->cmap->n - 1);
1249: av[indexn[j] * mat->lda + indexm[i]] = v ? v[idx++] : (idx++, 0.0);
1250: }
1251: }
1252: } else {
1253: for (i = 0; i < m; i++) {
1254: if (indexm[i] < 0) {
1255: idx += n;
1256: continue;
1257: }
1258: 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);
1259: for (j = 0; j < n; j++) {
1260: if (indexn[j] < 0) {
1261: idx++;
1262: continue;
1263: }
1264: 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);
1265: av[indexn[j] * mat->lda + indexm[i]] += v ? v[idx++] : (idx++, 0.0);
1266: }
1267: }
1268: }
1269: }
1270: /* hack to prevent unneeded copy to the GPU while returning the array */
1271: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
1272: oldf = A->offloadmask;
1273: A->offloadmask = PETSC_OFFLOAD_GPU;
1274: #endif
1275: PetscCall(MatDenseRestoreArray(A, &av));
1276: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
1277: A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1278: #endif
1279: PetscFunctionReturn(PETSC_SUCCESS);
1280: }
1282: static PetscErrorCode MatGetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], PetscScalar v[])
1283: {
1284: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1285: const PetscScalar *vv;
1286: PetscInt i, j;
1287: PetscBool roworiented = mat->roworiented;
1288: PetscScalar *value;
1290: PetscFunctionBegin;
1291: PetscCall(MatDenseGetArrayRead(A, &vv));
1292: for (i = 0; i < m; i++) {
1293: if (indexm[i] < 0) continue;
1294: 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);
1295: for (j = 0; j < n; j++) {
1296: if (indexn[j] < 0) continue;
1297: PetscCheck(indexn[j] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column %" PetscInt_FMT " requested larger than number columns %" PetscInt_FMT, indexn[j], A->cmap->n);
1298: value = roworiented ? &v[j + i * n] : &v[i + j * m];
1299: *value = vv[indexn[j] * mat->lda + indexm[i]];
1300: }
1301: }
1302: PetscCall(MatDenseRestoreArrayRead(A, &vv));
1303: PetscFunctionReturn(PETSC_SUCCESS);
1304: }
1306: PetscErrorCode MatView_Dense_Binary(Mat mat, PetscViewer viewer)
1307: {
1308: PetscBool skipHeader;
1309: PetscViewerFormat format;
1310: PetscInt header[4], M, N, m, lda, i, j;
1311: PetscCount k;
1312: const PetscScalar *v;
1313: PetscScalar *vwork;
1315: PetscFunctionBegin;
1316: PetscCall(PetscViewerSetUp(viewer));
1317: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
1318: PetscCall(PetscViewerGetFormat(viewer, &format));
1319: if (skipHeader) format = PETSC_VIEWER_NATIVE;
1321: PetscCall(MatGetSize(mat, &M, &N));
1323: /* write matrix header */
1324: header[0] = MAT_FILE_CLASSID;
1325: header[1] = M;
1326: header[2] = N;
1327: header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M * N;
1328: if (!skipHeader) PetscCall(PetscViewerBinaryWrite(viewer, header, 4, PETSC_INT));
1330: PetscCall(MatGetLocalSize(mat, &m, NULL));
1331: if (format != PETSC_VIEWER_NATIVE) {
1332: PetscInt nnz = m * N, *iwork;
1333: /* store row lengths for each row */
1334: PetscCall(PetscMalloc1(nnz, &iwork));
1335: for (i = 0; i < m; i++) iwork[i] = N;
1336: PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1337: /* store column indices (zero start index) */
1338: for (k = 0, i = 0; i < m; i++)
1339: for (j = 0; j < N; j++, k++) iwork[k] = j;
1340: PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1341: PetscCall(PetscFree(iwork));
1342: }
1343: /* store matrix values as a dense matrix in row major order */
1344: PetscCall(PetscMalloc1(m * N, &vwork));
1345: PetscCall(MatDenseGetArrayRead(mat, &v));
1346: PetscCall(MatDenseGetLDA(mat, &lda));
1347: for (k = 0, i = 0; i < m; i++)
1348: for (j = 0; j < N; j++, k++) vwork[k] = v[i + (size_t)lda * j];
1349: PetscCall(MatDenseRestoreArrayRead(mat, &v));
1350: PetscCall(PetscViewerBinaryWriteAll(viewer, vwork, m * N, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1351: PetscCall(PetscFree(vwork));
1352: PetscFunctionReturn(PETSC_SUCCESS);
1353: }
1355: PetscErrorCode MatLoad_Dense_Binary(Mat mat, PetscViewer viewer)
1356: {
1357: PetscBool skipHeader;
1358: PetscInt header[4], M, N, m, nz, lda, i, j, k;
1359: PetscInt rows, cols;
1360: PetscScalar *v, *vwork;
1362: PetscFunctionBegin;
1363: PetscCall(PetscViewerSetUp(viewer));
1364: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
1366: if (!skipHeader) {
1367: PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT));
1368: PetscCheck(header[0] == MAT_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file");
1369: M = header[1];
1370: N = header[2];
1371: PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M);
1372: PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N);
1373: nz = header[3];
1374: PetscCheck(nz == MATRIX_BINARY_FORMAT_DENSE || nz >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Unknown matrix format %" PetscInt_FMT " in file", nz);
1375: } else {
1376: PetscCall(MatGetSize(mat, &M, &N));
1377: 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");
1378: nz = MATRIX_BINARY_FORMAT_DENSE;
1379: }
1381: /* setup global sizes if not set */
1382: if (mat->rmap->N < 0) mat->rmap->N = M;
1383: if (mat->cmap->N < 0) mat->cmap->N = N;
1384: PetscCall(MatSetUp(mat));
1385: /* check if global sizes are correct */
1386: PetscCall(MatGetSize(mat, &rows, &cols));
1387: 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);
1389: PetscCall(MatGetSize(mat, NULL, &N));
1390: PetscCall(MatGetLocalSize(mat, &m, NULL));
1391: PetscCall(MatDenseGetArray(mat, &v));
1392: PetscCall(MatDenseGetLDA(mat, &lda));
1393: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */
1394: PetscCount nnz = (size_t)m * N;
1395: /* read in matrix values */
1396: PetscCall(PetscMalloc1(nnz, &vwork));
1397: PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1398: /* store values in column major order */
1399: for (j = 0; j < N; j++)
1400: for (i = 0; i < m; i++) v[i + (size_t)lda * j] = vwork[(size_t)i * N + j];
1401: PetscCall(PetscFree(vwork));
1402: } else { /* matrix in file is sparse format */
1403: PetscInt nnz = 0, *rlens, *icols;
1404: /* read in row lengths */
1405: PetscCall(PetscMalloc1(m, &rlens));
1406: PetscCall(PetscViewerBinaryReadAll(viewer, rlens, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1407: for (i = 0; i < m; i++) nnz += rlens[i];
1408: /* read in column indices and values */
1409: PetscCall(PetscMalloc2(nnz, &icols, nnz, &vwork));
1410: PetscCall(PetscViewerBinaryReadAll(viewer, icols, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1411: PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1412: /* store values in column major order */
1413: for (k = 0, i = 0; i < m; i++)
1414: for (j = 0; j < rlens[i]; j++, k++) v[i + lda * icols[k]] = vwork[k];
1415: PetscCall(PetscFree(rlens));
1416: PetscCall(PetscFree2(icols, vwork));
1417: }
1418: PetscCall(MatDenseRestoreArray(mat, &v));
1419: PetscFunctionReturn(PETSC_SUCCESS);
1420: }
1422: static PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1423: {
1424: PetscBool isbinary, ishdf5;
1426: PetscFunctionBegin;
1429: /* force binary viewer to load .info file if it has not yet done so */
1430: PetscCall(PetscViewerSetUp(viewer));
1431: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1432: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1433: if (isbinary) {
1434: PetscCall(MatLoad_Dense_Binary(newMat, viewer));
1435: } else if (ishdf5) {
1436: #if PetscDefined(HAVE_HDF5)
1437: PetscCall(MatLoad_Dense_HDF5(newMat, viewer));
1438: #else
1439: SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1440: #endif
1441: } else {
1442: SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "Viewer type %s not yet supported for reading %s matrices", ((PetscObject)viewer)->type_name, ((PetscObject)newMat)->type_name);
1443: }
1444: PetscFunctionReturn(PETSC_SUCCESS);
1445: }
1447: static PetscErrorCode MatView_SeqDense_ASCII(Mat A, PetscViewer viewer)
1448: {
1449: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1450: PetscInt i, j;
1451: const char *name;
1452: PetscScalar *v, *av;
1453: PetscViewerFormat format;
1454: PetscBool allreal = PETSC_TRUE;
1456: PetscFunctionBegin;
1457: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&av));
1458: PetscCall(PetscViewerGetFormat(viewer, &format));
1459: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1460: PetscFunctionReturn(PETSC_SUCCESS); /* do nothing for now */
1461: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1462: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1463: for (i = 0; i < A->rmap->n; i++) {
1464: v = av + i;
1465: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i));
1466: for (j = 0; j < A->cmap->n; j++) {
1467: if (PetscDefined(USE_COMPLEX) && PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1468: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", j, (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1469: } else if (PetscRealPart(*v)) {
1470: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)PetscRealPart(*v)));
1471: }
1472: v += a->lda;
1473: }
1474: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1475: }
1476: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1477: } else {
1478: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1479: if (PetscDefined(USE_COMPLEX)) {
1480: /* determine if matrix has all real values */
1481: for (j = 0; j < A->cmap->n; j++) {
1482: v = av + j * a->lda;
1483: for (i = 0; i < A->rmap->n; i++) {
1484: if (PetscImaginaryPart(v[i])) {
1485: allreal = PETSC_FALSE;
1486: break;
1487: }
1488: }
1489: }
1490: }
1491: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1492: PetscCall(PetscObjectGetName((PetscObject)A, &name));
1493: PetscCall(PetscViewerASCIIPrintf(viewer, "%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n", A->rmap->n, A->cmap->n));
1494: PetscCall(PetscViewerASCIIPrintf(viewer, "%s = zeros(%" PetscInt_FMT ",%" PetscInt_FMT ");\n", name, A->rmap->n, A->cmap->n));
1495: PetscCall(PetscViewerASCIIPrintf(viewer, "%s = [\n", name));
1496: }
1498: for (i = 0; i < A->rmap->n; i++) {
1499: v = av + i;
1500: for (j = 0; j < A->cmap->n; j++) {
1501: if (allreal) {
1502: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)PetscRealPart(*v)));
1503: } else {
1504: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei ", (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1505: }
1506: v += a->lda;
1507: }
1508: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1509: }
1510: if (format == PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer, "];\n"));
1511: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1512: }
1513: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&av));
1514: PetscCall(PetscViewerFlush(viewer));
1515: PetscFunctionReturn(PETSC_SUCCESS);
1516: }
1518: #include <petscdraw.h>
1519: #if defined(__GNUC__) && !defined(__clang__)
1520: #pragma GCC diagnostic push
1521: #pragma GCC diagnostic ignored "-Wclobbered"
1522: #endif
1523: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw, void *Aa)
1524: {
1525: Mat A = (Mat)Aa;
1526: PetscInt m = A->rmap->n, n = A->cmap->n, i, j;
1527: int color = PETSC_DRAW_WHITE;
1528: const PetscScalar *v;
1529: PetscViewer viewer;
1530: PetscReal xl, yl, xr, yr, x_l, x_r, y_l, y_r;
1531: PetscViewerFormat format;
1533: PetscFunctionBegin;
1534: PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer));
1535: PetscCall(PetscViewerGetFormat(viewer, &format));
1536: PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));
1538: /* Loop over matrix elements drawing boxes */
1539: PetscCall(MatDenseGetArrayRead(A, &v));
1540: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1541: PetscDrawCollectiveBegin(draw);
1542: /* Blue for negative and Red for positive */
1543: for (j = 0; j < n; j++) {
1544: x_l = j;
1545: x_r = x_l + 1.0;
1546: for (i = 0; i < m; i++) {
1547: y_l = m - i - 1.0;
1548: y_r = y_l + 1.0;
1549: if (PetscRealPart(v[j * m + i]) > 0.) color = PETSC_DRAW_RED;
1550: else if (PetscRealPart(v[j * m + i]) < 0.) color = PETSC_DRAW_BLUE;
1551: else continue;
1552: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1553: }
1554: }
1555: PetscDrawCollectiveEnd(draw);
1556: } else {
1557: /* use contour shading to indicate magnitude of values */
1558: /* first determine max of all nonzero values */
1559: PetscReal minv = 0.0, maxv = 0.0;
1560: PetscDraw popup;
1562: for (i = 0; i < m * n; i++) {
1563: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1564: }
1565: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1566: PetscCall(PetscDrawGetPopup(draw, &popup));
1567: PetscCall(PetscDrawScalePopup(popup, minv, maxv));
1569: PetscDrawCollectiveBegin(draw);
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: color = PetscDrawRealToColor(PetscAbsScalar(v[j * m + i]), minv, maxv);
1577: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1578: }
1579: }
1580: PetscDrawCollectiveEnd(draw);
1581: }
1582: PetscCall(MatDenseRestoreArrayRead(A, &v));
1583: PetscFunctionReturn(PETSC_SUCCESS);
1584: }
1585: #if defined(__GNUC__) && !defined(__clang__)
1586: #pragma GCC diagnostic pop
1587: #endif
1589: static PetscErrorCode MatView_SeqDense_Draw(Mat A, PetscViewer viewer)
1590: {
1591: PetscDraw draw;
1592: PetscBool isnull;
1593: PetscReal xr, yr, xl, yl, h, w;
1595: PetscFunctionBegin;
1596: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1597: PetscCall(PetscDrawIsNull(draw, &isnull));
1598: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
1600: xr = A->cmap->n;
1601: yr = A->rmap->n;
1602: h = yr / 10.0;
1603: w = xr / 10.0;
1604: xr += w;
1605: yr += h;
1606: xl = -w;
1607: yl = -h;
1608: PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr));
1609: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer));
1610: PetscCall(PetscDrawZoom(draw, MatView_SeqDense_Draw_Zoom, A));
1611: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL));
1612: PetscCall(PetscDrawSave(draw));
1613: PetscFunctionReturn(PETSC_SUCCESS);
1614: }
1616: PetscErrorCode MatView_SeqDense(Mat A, PetscViewer viewer)
1617: {
1618: PetscBool isascii, isbinary, isdraw;
1620: PetscFunctionBegin;
1621: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1622: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1623: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1624: if (isascii) PetscCall(MatView_SeqDense_ASCII(A, viewer));
1625: else if (isbinary) PetscCall(MatView_Dense_Binary(A, viewer));
1626: else if (isdraw) PetscCall(MatView_SeqDense_Draw(A, viewer));
1627: PetscFunctionReturn(PETSC_SUCCESS);
1628: }
1630: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A, const PetscScalar *array)
1631: {
1632: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1634: PetscFunctionBegin;
1635: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1636: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1637: PetscCheck(!a->unplacedarray, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseResetArray() first");
1638: a->unplacedarray = a->v;
1639: a->unplaced_user_alloc = a->user_alloc;
1640: a->v = (PetscScalar *)array;
1641: a->user_alloc = PETSC_TRUE;
1642: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
1643: A->offloadmask = PETSC_OFFLOAD_CPU;
1644: #endif
1645: PetscFunctionReturn(PETSC_SUCCESS);
1646: }
1648: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1649: {
1650: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1652: PetscFunctionBegin;
1653: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1654: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1655: a->v = a->unplacedarray;
1656: a->user_alloc = a->unplaced_user_alloc;
1657: a->unplacedarray = NULL;
1658: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
1659: A->offloadmask = PETSC_OFFLOAD_CPU;
1660: #endif
1661: PetscFunctionReturn(PETSC_SUCCESS);
1662: }
1664: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A, const PetscScalar *array)
1665: {
1666: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1668: PetscFunctionBegin;
1669: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1670: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1671: if (!a->user_alloc) PetscCall(PetscFree(a->v));
1672: a->v = (PetscScalar *)array;
1673: a->user_alloc = PETSC_FALSE;
1674: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
1675: A->offloadmask = PETSC_OFFLOAD_CPU;
1676: #endif
1677: PetscFunctionReturn(PETSC_SUCCESS);
1678: }
1680: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1681: {
1682: Mat_SeqDense *l = (Mat_SeqDense *)mat->data;
1684: PetscFunctionBegin;
1685: PetscCall(PetscLogObjectState((PetscObject)mat, "Rows %" PetscInt_FMT " Cols %" PetscInt_FMT, mat->rmap->n, mat->cmap->n));
1686: PetscCall(VecDestroy(&l->qrrhs));
1687: PetscCall(PetscFree(l->tau));
1688: PetscCall(PetscFree(l->pivots));
1689: PetscCall(PetscFree(l->fwork));
1690: if (!l->user_alloc) PetscCall(PetscFree(l->v));
1691: if (!l->unplaced_user_alloc) PetscCall(PetscFree(l->unplacedarray));
1692: PetscCheck(!l->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1693: PetscCheck(!l->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1694: PetscCall(VecDestroy(&l->cvec));
1695: PetscCall(MatDestroy(&l->cmat));
1696: PetscCall(PetscFree(mat->data));
1698: PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL));
1699: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactor_C", NULL));
1700: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorSymbolic_C", NULL));
1701: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorNumeric_C", NULL));
1702: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetLDA_C", NULL));
1703: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseSetLDA_C", NULL));
1704: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArray_C", NULL));
1705: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArray_C", NULL));
1706: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDensePlaceArray_C", NULL));
1707: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseResetArray_C", NULL));
1708: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseReplaceArray_C", NULL));
1709: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayRead_C", NULL));
1710: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayRead_C", NULL));
1711: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayWrite_C", NULL));
1712: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayWrite_C", NULL));
1713: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqaij_C", NULL));
1714: #if PetscDefined(HAVE_ELEMENTAL)
1715: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_elemental_C", NULL));
1716: #endif
1717: #if PetscDefined(HAVE_SCALAPACK) && (PetscDefined(USE_REAL_SINGLE) || PetscDefined(USE_REAL_DOUBLE))
1718: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_scalapack_C", NULL));
1719: #endif
1720: #if PetscDefined(HAVE_CUDA)
1721: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensecuda_C", NULL));
1722: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", NULL));
1723: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdense_C", NULL));
1724: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensecuda_C", NULL));
1725: #endif
1726: #if PetscDefined(HAVE_HIP)
1727: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensehip_C", NULL));
1728: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", NULL));
1729: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdense_C", NULL));
1730: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensehip_C", NULL));
1731: #endif
1732: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatSeqDenseSetPreallocation_C", NULL));
1733: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqaij_seqdense_C", NULL));
1734: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdense_C", NULL));
1735: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqbaij_seqdense_C", NULL));
1736: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqsbaij_seqdense_C", NULL));
1738: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumn_C", NULL));
1739: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumn_C", NULL));
1740: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVec_C", NULL));
1741: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVec_C", NULL));
1742: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecRead_C", NULL));
1743: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecRead_C", NULL));
1744: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecWrite_C", NULL));
1745: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecWrite_C", NULL));
1746: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetSubMatrix_C", NULL));
1747: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreSubMatrix_C", NULL));
1748: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultColumnRange_C", NULL));
1749: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultAddColumnRange_C", NULL));
1750: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeColumnRange_C", NULL));
1751: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeAddColumnRange_C", NULL));
1752: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseUpdateColumnLayout_C", NULL));
1753: PetscFunctionReturn(PETSC_SUCCESS);
1754: }
1756: static PetscErrorCode MatTranspose_SeqDense(Mat A, MatReuse reuse, Mat *matout)
1757: {
1758: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1759: PetscInt k, j, m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1760: PetscScalar *v, tmp;
1762: PetscFunctionBegin;
1763: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *matout));
1764: if (reuse == MAT_INPLACE_MATRIX) {
1765: if (m == n) { /* in place transpose */
1766: PetscCall(MatDenseGetArray(A, &v));
1767: for (j = 0; j < m; j++) {
1768: for (k = 0; k < j; k++) {
1769: tmp = v[j + k * M];
1770: v[j + k * M] = v[k + j * M];
1771: v[k + j * M] = tmp;
1772: }
1773: }
1774: PetscCall(MatDenseRestoreArray(A, &v));
1775: } else { /* reuse memory, temporary allocates new memory */
1776: PetscScalar *v2;
1777: PetscLayout tmplayout;
1779: PetscCall(PetscMalloc1((size_t)m * n, &v2));
1780: PetscCall(MatDenseGetArray(A, &v));
1781: for (j = 0; j < n; j++) {
1782: for (k = 0; k < m; k++) v2[j + (size_t)k * n] = v[k + (size_t)j * M];
1783: }
1784: PetscCall(PetscArraycpy(v, v2, (size_t)m * n));
1785: PetscCall(PetscFree(v2));
1786: PetscCall(MatDenseRestoreArray(A, &v));
1787: /* cleanup size dependent quantities */
1788: PetscCall(VecDestroy(&mat->cvec));
1789: PetscCall(MatDestroy(&mat->cmat));
1790: PetscCall(PetscFree(mat->pivots));
1791: PetscCall(PetscFree(mat->fwork));
1792: /* swap row/col layouts */
1793: PetscCall(PetscBLASIntCast(n, &mat->lda));
1794: tmplayout = A->rmap;
1795: A->rmap = A->cmap;
1796: A->cmap = tmplayout;
1797: }
1798: } else { /* out-of-place transpose */
1799: Mat tmat;
1800: Mat_SeqDense *tmatd;
1801: PetscScalar *v2;
1802: PetscInt M2;
1804: if (reuse == MAT_INITIAL_MATRIX) {
1805: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &tmat));
1806: PetscCall(MatSetSizes(tmat, A->cmap->n, A->rmap->n, A->cmap->n, A->rmap->n));
1807: PetscCall(MatSetType(tmat, ((PetscObject)A)->type_name));
1808: PetscCall(MatSeqDenseSetPreallocation(tmat, NULL));
1809: } else tmat = *matout;
1811: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&v));
1812: PetscCall(MatDenseGetArray(tmat, &v2));
1813: tmatd = (Mat_SeqDense *)tmat->data;
1814: M2 = tmatd->lda;
1815: for (j = 0; j < n; j++) {
1816: for (k = 0; k < m; k++) v2[j + k * M2] = v[k + j * M];
1817: }
1818: PetscCall(MatDenseRestoreArray(tmat, &v2));
1819: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&v));
1820: *matout = tmat;
1821: }
1822: PetscFunctionReturn(PETSC_SUCCESS);
1823: }
1825: static PetscErrorCode MatEqual_SeqDense(Mat A1, Mat A2, PetscBool *flg)
1826: {
1827: Mat_SeqDense *mat1 = (Mat_SeqDense *)A1->data;
1828: Mat_SeqDense *mat2 = (Mat_SeqDense *)A2->data;
1829: PetscInt i;
1830: const PetscScalar *v1, *v2;
1832: PetscFunctionBegin;
1833: if (A1->rmap->n != A2->rmap->n) {
1834: *flg = PETSC_FALSE;
1835: PetscFunctionReturn(PETSC_SUCCESS);
1836: }
1837: if (A1->cmap->n != A2->cmap->n) {
1838: *flg = PETSC_FALSE;
1839: PetscFunctionReturn(PETSC_SUCCESS);
1840: }
1841: PetscCall(MatDenseGetArrayRead(A1, &v1));
1842: PetscCall(MatDenseGetArrayRead(A2, &v2));
1843: for (i = 0; i < A1->cmap->n; i++) {
1844: PetscCall(PetscArraycmp(v1, v2, A1->rmap->n, flg));
1845: if (*flg == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS);
1846: v1 += mat1->lda;
1847: v2 += mat2->lda;
1848: }
1849: PetscCall(MatDenseRestoreArrayRead(A1, &v1));
1850: PetscCall(MatDenseRestoreArrayRead(A2, &v2));
1851: *flg = PETSC_TRUE;
1852: PetscFunctionReturn(PETSC_SUCCESS);
1853: }
1855: PetscErrorCode MatGetDiagonal_SeqDense(Mat A, Vec v)
1856: {
1857: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1858: PetscInt i, n, len;
1859: PetscScalar *x;
1860: const PetscScalar *vv;
1862: PetscFunctionBegin;
1863: PetscCall(VecGetSize(v, &n));
1864: PetscCall(VecGetArray(v, &x));
1865: len = PetscMin(A->rmap->n, A->cmap->n);
1866: PetscCall(MatDenseGetArrayRead(A, &vv));
1867: PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming mat and vec");
1868: for (i = 0; i < len; i++) x[i] = vv[i * mat->lda + i];
1869: PetscCall(MatDenseRestoreArrayRead(A, &vv));
1870: PetscCall(VecRestoreArray(v, &x));
1871: PetscFunctionReturn(PETSC_SUCCESS);
1872: }
1874: PetscErrorCode MatDiagonalScale_SeqDense(Mat A, Vec ll, Vec rr)
1875: {
1876: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1877: const PetscScalar *l, *r;
1878: PetscScalar x, *v, *vv;
1879: PetscInt i, j, m = A->rmap->n, n = A->cmap->n;
1881: PetscFunctionBegin;
1882: PetscCall(MatDenseGetArray(A, &vv));
1883: if (ll) {
1884: PetscCall(VecGetLocalSize(ll, &m));
1885: PetscCall(VecGetArrayRead(ll, &l));
1886: PetscCheck(m == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vec wrong size");
1887: for (i = 0; i < m; i++) {
1888: x = l[i];
1889: v = vv + i;
1890: for (j = 0; j < n; j++) {
1891: (*v) *= x;
1892: v += mat->lda;
1893: }
1894: }
1895: PetscCall(VecRestoreArrayRead(ll, &l));
1896: PetscCall(PetscLogFlops(1.0 * n * m));
1897: }
1898: if (rr) {
1899: PetscCall(VecGetLocalSize(rr, &n));
1900: PetscCall(VecGetArrayRead(rr, &r));
1901: PetscCheck(n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vec wrong size");
1902: for (i = 0; i < n; i++) {
1903: x = r[i];
1904: v = vv + i * mat->lda;
1905: for (j = 0; j < m; j++) (*v++) *= x;
1906: }
1907: PetscCall(VecRestoreArrayRead(rr, &r));
1908: PetscCall(PetscLogFlops(1.0 * n * m));
1909: }
1910: PetscCall(MatDenseRestoreArray(A, &vv));
1911: PetscFunctionReturn(PETSC_SUCCESS);
1912: }
1914: PetscErrorCode MatNorm_SeqDense(Mat A, NormType type, PetscReal *nrm)
1915: {
1916: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1917: PetscScalar *v, *vv, *work, *av = NULL;
1918: PetscReal sum = 0.0;
1919: PetscInt lda, m = A->rmap->n, i, j;
1921: PetscFunctionBegin;
1922: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&vv));
1923: PetscCall(MatDenseGetLDA(A, &lda));
1924: v = vv;
1925: if (type == NORM_FROBENIUS) {
1926: if (lda > m) {
1927: for (j = 0; j < A->cmap->n; j++) {
1928: v = vv + j * lda;
1929: for (i = 0; i < m; i++) {
1930: sum += PetscRealPart(PetscConj(*v) * (*v));
1931: v++;
1932: }
1933: }
1934: } else {
1935: #if PetscDefined(USE_REAL___FP16)
1936: PetscBLASInt one = 1, cnt = A->cmap->n * A->rmap->n;
1937: PetscCallBLAS("BLASnrm2", *nrm = BLASnrm2_(&cnt, v, &one));
1938: }
1939: #else
1940: for (i = 0; i < A->cmap->n * A->rmap->n; i++) {
1941: sum += PetscRealPart(PetscConj(*v) * (*v));
1942: v++;
1943: }
1944: }
1945: *nrm = PetscSqrtReal(sum);
1946: #endif
1947: PetscCall(PetscLogFlops(2.0 * A->cmap->n * A->rmap->n));
1948: } else if (type == NORM_1) {
1949: *nrm = 0.0;
1950: for (j = 0; j < A->cmap->n; j++) {
1951: v = vv + j * mat->lda;
1952: sum = 0.0;
1953: for (i = 0; i < A->rmap->n; i++) {
1954: sum += PetscAbsScalar(*v);
1955: v++;
1956: }
1957: if (sum > *nrm) *nrm = sum;
1958: }
1959: PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
1960: } else if (type == NORM_INFINITY) {
1961: *nrm = 0.0;
1962: for (j = 0; j < A->rmap->n; j++) {
1963: v = vv + j;
1964: sum = 0.0;
1965: for (i = 0; i < A->cmap->n; i++) {
1966: sum += PetscAbsScalar(*v);
1967: v += mat->lda;
1968: }
1969: if (sum > *nrm) *nrm = sum;
1970: }
1971: PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
1972: } else if (type == NORM_2) {
1973: PetscReal *s;
1974: PetscBLASInt bm, bn, blda, min, lwork;
1976: PetscCall(PetscBLASIntCast(A->rmap->n, &bm));
1977: PetscCall(PetscBLASIntCast(A->cmap->n, &bn));
1978: PetscCall(PetscBLASIntCast(PetscMax(A->rmap->n, 1), &blda));
1979: min = PetscMin(bm, bn);
1980: if (!min) {
1981: *nrm = 0.0;
1982: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&vv));
1983: PetscFunctionReturn(PETSC_SUCCESS);
1984: }
1985: PetscCall(PetscMalloc2(A->rmap->n * A->cmap->n, &av, min, &s));
1986: for (j = 0; j < A->cmap->n; j++) PetscCall(PetscArraycpy(av + j * A->rmap->n, vv + j * lda, A->rmap->n));
1988: lwork = -1;
1989: {
1990: PetscScalar workquery;
1991: #if PetscDefined(USE_COMPLEX)
1992: PetscReal *rwork;
1994: PetscCall(PetscMalloc1(5 * min, &rwork));
1995: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
1996: PetscCallLAPACKInfo("LAPACKgesvd", LAPACKgesvd_("N", "N", &bm, &bn, av, &blda, s, NULL, &bm, NULL, &min, &workquery, &lwork, rwork, &info));
1997: lwork = (PetscBLASInt)PetscRealPart(workquery);
1998: PetscCall(PetscMalloc1(lwork, &work));
1999: PetscCallLAPACKInfo("LAPACKgesvd", LAPACKgesvd_("N", "N", &bm, &bn, av, &blda, s, NULL, &bm, NULL, &min, work, &lwork, rwork, &info));
2000: PetscCall(PetscFPTrapPop());
2001: PetscCall(PetscFree(rwork));
2002: #else
2003: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2004: PetscCallLAPACKInfo("LAPACKgesvd", LAPACKgesvd_("N", "N", &bm, &bn, av, &blda, s, NULL, &bm, NULL, &min, &workquery, &lwork, &info));
2005: lwork = (PetscBLASInt)PetscRealPart(workquery);
2006: PetscCall(PetscMalloc1(lwork, &work));
2007: PetscCallLAPACKInfo("LAPACKgesvd", LAPACKgesvd_("N", "N", &bm, &bn, av, &blda, s, NULL, &bm, NULL, &min, work, &lwork, &info));
2008: PetscCall(PetscFPTrapPop());
2009: #endif
2010: }
2011: *nrm = s[0];
2012: PetscCall(PetscFree(work));
2013: PetscCall(PetscFree2(av, s));
2014: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported norm type %s", NormTypes[type]);
2015: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&vv));
2016: PetscFunctionReturn(PETSC_SUCCESS);
2017: }
2019: static PetscErrorCode MatSetOption_SeqDense(Mat A, MatOption op, PetscBool flg)
2020: {
2021: Mat_SeqDense *aij = (Mat_SeqDense *)A->data;
2023: PetscFunctionBegin;
2024: switch (op) {
2025: case MAT_ROW_ORIENTED:
2026: aij->roworiented = flg;
2027: break;
2028: default:
2029: break;
2030: }
2031: PetscFunctionReturn(PETSC_SUCCESS);
2032: }
2034: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
2035: {
2036: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
2037: PetscInt lda = l->lda, m = A->rmap->n, n = A->cmap->n, j;
2038: PetscScalar *v;
2040: PetscFunctionBegin;
2041: PetscCall(MatDenseGetArrayWrite(A, &v));
2042: if (lda > m) {
2043: for (j = 0; j < n; j++) PetscCall(PetscArrayzero(v + PetscInt64Mult(j, lda), m));
2044: } else {
2045: PetscCall(PetscArrayzero(v, PetscInt64Mult(m, n)));
2046: }
2047: PetscCall(MatDenseRestoreArrayWrite(A, &v));
2048: PetscFunctionReturn(PETSC_SUCCESS);
2049: }
2051: static PetscErrorCode MatSetInf_SeqDense(Mat A)
2052: {
2053: // MSVC gives "divide by zero" error at compile time - so declare as volatile to skip this check.
2054: volatile PetscReal one = 1.0, zero = 0.0;
2055: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
2056: PetscInt lda = l->lda, m = A->rmap->n, n = A->cmap->n, i, j;
2057: PetscScalar *v;
2058: PetscScalar inf;
2060: PetscFunctionBegin;
2061: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
2062: inf = one / zero;
2063: PetscCall(PetscFPTrapPop());
2064: PetscCall(MatDenseGetArrayWrite(A, &v));
2065: if (lda > m) {
2066: for (j = 0; j < n; j++) {
2067: PetscScalar *vj = v + PetscInt64Mult(j, lda);
2069: for (i = 0; i < m; i++) vj[i] = inf;
2070: }
2071: } else {
2072: const PetscInt64 mn = PetscInt64Mult(m, n);
2074: for (PetscInt64 k = 0; k < mn; k++) v[k] = inf;
2075: }
2076: PetscCall(MatDenseRestoreArrayWrite(A, &v));
2077: PetscFunctionReturn(PETSC_SUCCESS);
2078: }
2080: static PetscErrorCode MatZeroRows_SeqDense(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2081: {
2082: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
2083: PetscInt m = l->lda, n = A->cmap->n, i, j;
2084: PetscScalar *slot, *bb, *v;
2085: const PetscScalar *xx;
2087: PetscFunctionBegin;
2088: if (PetscDefined(USE_DEBUG)) {
2089: for (i = 0; i < N; i++) {
2090: PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
2091: 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);
2092: }
2093: }
2094: if (!N) PetscFunctionReturn(PETSC_SUCCESS);
2096: /* fix right-hand side if needed */
2097: if (x && b) {
2098: PetscCall(VecGetArrayRead(x, &xx));
2099: PetscCall(VecGetArray(b, &bb));
2100: for (i = 0; i < N; i++) bb[rows[i]] = diag * xx[rows[i]];
2101: PetscCall(VecRestoreArrayRead(x, &xx));
2102: PetscCall(VecRestoreArray(b, &bb));
2103: }
2105: PetscCall(MatDenseGetArray(A, &v));
2106: for (i = 0; i < N; i++) {
2107: slot = v + rows[i];
2108: for (j = 0; j < n; j++) {
2109: *slot = 0.0;
2110: slot += m;
2111: }
2112: }
2113: if (diag != 0.0) {
2114: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
2115: for (i = 0; i < N; i++) {
2116: slot = v + (m + 1) * rows[i];
2117: *slot = diag;
2118: }
2119: }
2120: PetscCall(MatDenseRestoreArray(A, &v));
2121: PetscFunctionReturn(PETSC_SUCCESS);
2122: }
2124: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A, PetscInt *lda)
2125: {
2126: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2128: PetscFunctionBegin;
2129: *lda = mat->lda;
2130: PetscFunctionReturn(PETSC_SUCCESS);
2131: }
2133: PetscErrorCode MatDenseGetArray_SeqDense(Mat A, PetscScalar **array)
2134: {
2135: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2137: PetscFunctionBegin;
2138: PetscCheck(!mat->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
2139: *array = mat->v;
2140: PetscFunctionReturn(PETSC_SUCCESS);
2141: }
2143: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A, PetscScalar **array)
2144: {
2145: PetscFunctionBegin;
2146: if (array) *array = NULL;
2147: PetscFunctionReturn(PETSC_SUCCESS);
2148: }
2150: /*@
2151: MatDenseGetLDA - gets the leading dimension of the array returned from `MatDenseGetArray()`
2153: Not Collective
2155: Input Parameter:
2156: . A - a `MATDENSE` or `MATDENSECUDA` matrix
2158: Output Parameter:
2159: . lda - the leading dimension
2161: Level: intermediate
2163: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseSetLDA()`
2164: @*/
2165: PetscErrorCode MatDenseGetLDA(Mat A, PetscInt *lda)
2166: {
2167: PetscFunctionBegin;
2169: PetscAssertPointer(lda, 2);
2170: MatCheckPreallocated(A, 1);
2171: PetscUseMethod(A, "MatDenseGetLDA_C", (Mat, PetscInt *), (A, lda));
2172: PetscFunctionReturn(PETSC_SUCCESS);
2173: }
2175: /*@
2176: MatDenseSetLDA - Sets the leading dimension of the array used by the `MATDENSE` matrix
2178: Collective if the matrix layouts have not yet been setup
2180: Input Parameters:
2181: + A - a `MATDENSE` or `MATDENSECUDA` matrix
2182: - lda - the leading dimension
2184: Level: intermediate
2186: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetLDA()`
2187: @*/
2188: PetscErrorCode MatDenseSetLDA(Mat A, PetscInt lda)
2189: {
2190: PetscFunctionBegin;
2192: PetscTryMethod(A, "MatDenseSetLDA_C", (Mat, PetscInt), (A, lda));
2193: PetscFunctionReturn(PETSC_SUCCESS);
2194: }
2196: /*@
2197: MatDenseGetArray - gives read-write access to the array where the data for a `MATDENSE` matrix is stored
2199: Logically Collective
2201: Input Parameter:
2202: . A - a dense matrix
2204: Output Parameter:
2205: . array - pointer to the data
2207: Level: intermediate
2209: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2210: @*/
2211: PetscErrorCode MatDenseGetArray(Mat A, PetscScalar *array[]) PeNS
2212: {
2213: PetscFunctionBegin;
2215: PetscAssertPointer(array, 2);
2216: PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2217: PetscFunctionReturn(PETSC_SUCCESS);
2218: }
2220: /*@
2221: MatDenseRestoreArray - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArray()`
2223: Logically Collective
2225: Input Parameters:
2226: + A - a dense matrix
2227: - array - pointer to the data (may be `NULL`)
2229: Level: intermediate
2231: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2232: @*/
2233: PetscErrorCode MatDenseRestoreArray(Mat A, PetscScalar *array[]) PeNS
2234: {
2235: PetscFunctionBegin;
2237: if (array) PetscAssertPointer(array, 2);
2238: PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2239: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2240: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
2241: A->offloadmask = PETSC_OFFLOAD_CPU;
2242: #endif
2243: PetscFunctionReturn(PETSC_SUCCESS);
2244: }
2246: /*@
2247: MatDenseGetArrayRead - gives read-only access to the array where the data for a `MATDENSE` matrix is stored
2249: Not Collective
2251: Input Parameter:
2252: . A - a dense matrix
2254: Output Parameter:
2255: . array - pointer to the data
2257: Level: intermediate
2259: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2260: @*/
2261: PetscErrorCode MatDenseGetArrayRead(Mat A, const PetscScalar *array[]) PeNS
2262: {
2263: PetscFunctionBegin;
2265: PetscAssertPointer(array, 2);
2266: PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2267: PetscFunctionReturn(PETSC_SUCCESS);
2268: }
2270: /*@
2271: MatDenseRestoreArrayRead - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayRead()`
2273: Not Collective
2275: Input Parameters:
2276: + A - a dense matrix
2277: - array - pointer to the data (may be `NULL`)
2279: Level: intermediate
2281: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2282: @*/
2283: PetscErrorCode MatDenseRestoreArrayRead(Mat A, const PetscScalar *array[]) PeNS
2284: {
2285: PetscFunctionBegin;
2287: if (array) PetscAssertPointer(array, 2);
2288: PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2289: PetscFunctionReturn(PETSC_SUCCESS);
2290: }
2292: /*@
2293: MatDenseGetArrayWrite - gives write-only access to the array where the data for a `MATDENSE` matrix is stored
2295: Not Collective
2297: Input Parameter:
2298: . A - a dense matrix
2300: Output Parameter:
2301: . array - pointer to the data
2303: Level: intermediate
2305: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2306: @*/
2307: PetscErrorCode MatDenseGetArrayWrite(Mat A, PetscScalar *array[]) PeNS
2308: {
2309: PetscFunctionBegin;
2311: PetscAssertPointer(array, 2);
2312: PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2313: PetscFunctionReturn(PETSC_SUCCESS);
2314: }
2316: /*@
2317: MatDenseRestoreArrayWrite - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayWrite()`
2319: Not Collective
2321: Input Parameters:
2322: + A - a dense matrix
2323: - array - pointer to the data (may be `NULL`)
2325: Level: intermediate
2327: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2328: @*/
2329: PetscErrorCode MatDenseRestoreArrayWrite(Mat A, PetscScalar *array[]) PeNS
2330: {
2331: PetscFunctionBegin;
2333: if (array) PetscAssertPointer(array, 2);
2334: PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2335: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2336: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
2337: A->offloadmask = PETSC_OFFLOAD_CPU;
2338: #endif
2339: PetscFunctionReturn(PETSC_SUCCESS);
2340: }
2342: /*@
2343: MatDenseGetArrayAndMemType - gives read-write access to the array where the data for a `MATDENSE` matrix is stored
2345: Logically Collective
2347: Input Parameter:
2348: . A - a dense matrix
2350: Output Parameters:
2351: + array - pointer to the data
2352: - mtype - memory type of the returned pointer
2354: Level: intermediate
2356: Note:
2357: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2358: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2360: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArrayRead()`,
2361: `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2362: @*/
2363: PetscErrorCode MatDenseGetArrayAndMemType(Mat A, PetscScalar *array[], PetscMemType *mtype)
2364: {
2365: PetscBool isMPI;
2367: PetscFunctionBegin;
2369: PetscAssertPointer(array, 2);
2370: 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 */
2371: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2372: if (isMPI) {
2373: /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2374: PetscCall(MatDenseGetArrayAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2375: } else {
2376: PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);
2378: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayAndMemType_C", &fptr));
2379: if (fptr) {
2380: PetscCall((*fptr)(A, array, mtype));
2381: } else {
2382: PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2383: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2384: }
2385: }
2386: PetscFunctionReturn(PETSC_SUCCESS);
2387: }
2389: /*@
2390: MatDenseRestoreArrayAndMemType - returns access to the array that is obtained by `MatDenseGetArrayAndMemType()`
2392: Logically Collective
2394: Input Parameters:
2395: + A - a dense matrix
2396: - array - pointer to the data
2398: Level: intermediate
2400: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2401: @*/
2402: PetscErrorCode MatDenseRestoreArrayAndMemType(Mat A, PetscScalar *array[])
2403: {
2404: PetscBool isMPI;
2406: PetscFunctionBegin;
2408: if (array) PetscAssertPointer(array, 2);
2409: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2410: if (isMPI) {
2411: PetscCall(MatDenseRestoreArrayAndMemType(((Mat_MPIDense *)A->data)->A, array));
2412: } else {
2413: PetscErrorCode (*fptr)(Mat, PetscScalar **);
2415: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayAndMemType_C", &fptr));
2416: if (fptr) {
2417: PetscCall((*fptr)(A, array));
2418: } else {
2419: PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2420: }
2421: if (array) *array = NULL;
2422: }
2423: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2424: PetscFunctionReturn(PETSC_SUCCESS);
2425: }
2427: /*@
2428: MatDenseGetArrayReadAndMemType - gives read-only access to the array where the data for a `MATDENSE` matrix is stored
2430: Logically Collective
2432: Input Parameter:
2433: . A - a dense matrix
2435: Output Parameters:
2436: + array - pointer to the data
2437: - mtype - memory type of the returned pointer
2439: Level: intermediate
2441: Note:
2442: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2443: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2445: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`,
2446: `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2447: @*/
2448: PetscErrorCode MatDenseGetArrayReadAndMemType(Mat A, const PetscScalar *array[], PetscMemType *mtype)
2449: {
2450: PetscBool isMPI;
2452: PetscFunctionBegin;
2454: PetscAssertPointer(array, 2);
2455: 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 */
2456: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2457: if (isMPI) { /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2458: PetscCall(MatDenseGetArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2459: } else {
2460: PetscErrorCode (*fptr)(Mat, const PetscScalar **, PetscMemType *);
2462: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayReadAndMemType_C", &fptr));
2463: if (fptr) {
2464: PetscCall((*fptr)(A, array, mtype));
2465: } else {
2466: PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2467: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2468: }
2469: }
2470: PetscFunctionReturn(PETSC_SUCCESS);
2471: }
2473: /*@
2474: MatDenseRestoreArrayReadAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()`
2476: Logically Collective
2478: Input Parameters:
2479: + A - a dense matrix
2480: - array - pointer to the data
2482: Level: intermediate
2484: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2485: @*/
2486: PetscErrorCode MatDenseRestoreArrayReadAndMemType(Mat A, const PetscScalar *array[])
2487: {
2488: PetscBool isMPI;
2490: PetscFunctionBegin;
2492: if (array) PetscAssertPointer(array, 2);
2493: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2494: if (isMPI) {
2495: PetscCall(MatDenseRestoreArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array));
2496: } else {
2497: PetscErrorCode (*fptr)(Mat, const PetscScalar **);
2499: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayReadAndMemType_C", &fptr));
2500: if (fptr) {
2501: PetscCall((*fptr)(A, array));
2502: } else {
2503: PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2504: }
2505: if (array) *array = NULL;
2506: }
2507: PetscFunctionReturn(PETSC_SUCCESS);
2508: }
2510: /*@
2511: MatDenseGetArrayWriteAndMemType - gives write-only access to the array where the data for a `MATDENSE` matrix is stored
2513: Logically Collective
2515: Input Parameter:
2516: . A - a dense matrix
2518: Output Parameters:
2519: + array - pointer to the data
2520: - mtype - memory type of the returned pointer
2522: Level: intermediate
2524: Note:
2525: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2526: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2528: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWriteAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayRead()`,
2529: `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2530: @*/
2531: PetscErrorCode MatDenseGetArrayWriteAndMemType(Mat A, PetscScalar *array[], PetscMemType *mtype)
2532: {
2533: PetscBool isMPI;
2535: PetscFunctionBegin;
2537: PetscAssertPointer(array, 2);
2538: 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 */
2539: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2540: if (isMPI) {
2541: PetscCall(MatDenseGetArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2542: } else {
2543: PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);
2545: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayWriteAndMemType_C", &fptr));
2546: if (fptr) {
2547: PetscCall((*fptr)(A, array, mtype));
2548: } else {
2549: PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2550: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2551: }
2552: }
2553: PetscFunctionReturn(PETSC_SUCCESS);
2554: }
2556: /*@
2557: MatDenseRestoreArrayWriteAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()`
2559: Logically Collective
2561: Input Parameters:
2562: + A - a dense matrix
2563: - array - pointer to the data
2565: Level: intermediate
2567: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2568: @*/
2569: PetscErrorCode MatDenseRestoreArrayWriteAndMemType(Mat A, PetscScalar *array[])
2570: {
2571: PetscBool isMPI;
2573: PetscFunctionBegin;
2575: if (array) PetscAssertPointer(array, 2);
2576: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2577: if (isMPI) {
2578: PetscCall(MatDenseRestoreArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array));
2579: } else {
2580: PetscErrorCode (*fptr)(Mat, PetscScalar **);
2582: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayWriteAndMemType_C", &fptr));
2583: if (fptr) {
2584: PetscCall((*fptr)(A, array));
2585: } else {
2586: PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2587: }
2588: if (array) *array = NULL;
2589: }
2590: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2591: PetscFunctionReturn(PETSC_SUCCESS);
2592: }
2594: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B)
2595: {
2596: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2597: PetscInt i, j, nrows, ncols, ldb;
2598: const PetscInt *irow, *icol;
2599: PetscScalar *av, *bv, *v = mat->v;
2600: Mat newmat;
2602: PetscFunctionBegin;
2603: PetscCall(ISGetIndices(isrow, &irow));
2604: PetscCall(ISGetIndices(iscol, &icol));
2605: PetscCall(ISGetLocalSize(isrow, &nrows));
2606: PetscCall(ISGetLocalSize(iscol, &ncols));
2608: /* Check submatrixcall */
2609: if (scall == MAT_REUSE_MATRIX) {
2610: PetscInt n_cols, n_rows;
2611: PetscCall(MatGetSize(*B, &n_rows, &n_cols));
2612: if (n_rows != nrows || n_cols != ncols) {
2613: /* resize the result matrix to match number of requested rows/columns */
2614: PetscCall(MatSetSizes(*B, nrows, ncols, nrows, ncols));
2615: }
2616: newmat = *B;
2617: } else {
2618: /* Create and fill new matrix */
2619: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &newmat));
2620: PetscCall(MatSetSizes(newmat, nrows, ncols, nrows, ncols));
2621: PetscCall(MatSetType(newmat, ((PetscObject)A)->type_name));
2622: PetscCall(MatSeqDenseSetPreallocation(newmat, NULL));
2623: }
2625: /* Now extract the data pointers and do the copy,column at a time */
2626: PetscCall(MatDenseGetArray(newmat, &bv));
2627: PetscCall(MatDenseGetLDA(newmat, &ldb));
2628: for (i = 0; i < ncols; i++) {
2629: av = v + mat->lda * icol[i];
2630: for (j = 0; j < nrows; j++) bv[j] = av[irow[j]];
2631: bv += ldb;
2632: }
2633: PetscCall(MatDenseRestoreArray(newmat, &bv));
2635: /* Free work space */
2636: PetscCall(ISRestoreIndices(isrow, &irow));
2637: PetscCall(ISRestoreIndices(iscol, &icol));
2638: *B = newmat;
2639: PetscFunctionReturn(PETSC_SUCCESS);
2640: }
2642: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[])
2643: {
2644: PetscInt i;
2646: PetscFunctionBegin;
2647: if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n, B));
2649: for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SeqDense(A, irow[i], icol[i], scall, &(*B)[i]));
2650: PetscFunctionReturn(PETSC_SUCCESS);
2651: }
2653: PetscErrorCode MatCopy_SeqDense(Mat A, Mat B, MatStructure str)
2654: {
2655: Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data;
2656: const PetscScalar *va;
2657: PetscScalar *vb;
2658: PetscInt lda1 = a->lda, lda2 = b->lda, m = A->rmap->n, n = A->cmap->n, j;
2660: PetscFunctionBegin;
2661: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2662: if (A->ops->copy != B->ops->copy) {
2663: PetscCall(MatCopy_Basic(A, B, str));
2664: PetscFunctionReturn(PETSC_SUCCESS);
2665: }
2666: PetscCheck(m == B->rmap->n && n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "size(B) != size(A)");
2667: PetscCall(MatDenseGetArrayRead(A, &va));
2668: PetscCall(MatDenseGetArray(B, &vb));
2669: if (lda1 > m || lda2 > m) {
2670: for (j = 0; j < n; j++) PetscCall(PetscArraycpy(vb + j * lda2, va + j * lda1, m));
2671: } else {
2672: PetscCall(PetscArraycpy(vb, va, A->rmap->n * A->cmap->n));
2673: }
2674: PetscCall(MatDenseRestoreArray(B, &vb));
2675: PetscCall(MatDenseRestoreArrayRead(A, &va));
2676: PetscFunctionReturn(PETSC_SUCCESS);
2677: }
2679: PetscErrorCode MatSetUp_SeqDense(Mat A)
2680: {
2681: PetscFunctionBegin;
2682: PetscCall(PetscLayoutSetUp(A->rmap));
2683: PetscCall(PetscLayoutSetUp(A->cmap));
2684: if (!A->preallocated) PetscCall(MatSeqDenseSetPreallocation(A, NULL));
2685: PetscFunctionReturn(PETSC_SUCCESS);
2686: }
2688: PetscErrorCode MatConjugate_SeqDense(Mat A)
2689: {
2690: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2691: PetscInt i, j;
2692: PetscInt min = PetscMin(A->rmap->n, A->cmap->n);
2693: PetscScalar *aa;
2695: PetscFunctionBegin;
2696: PetscCall(MatDenseGetArray(A, &aa));
2697: for (j = 0; j < A->cmap->n; j++)
2698: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscConj(aa[i + j * mat->lda]);
2699: PetscCall(MatDenseRestoreArray(A, &aa));
2700: if (mat->tau)
2701: for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
2702: PetscFunctionReturn(PETSC_SUCCESS);
2703: }
2705: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2706: {
2707: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2708: PetscInt i, j;
2709: PetscScalar *aa;
2711: PetscFunctionBegin;
2712: PetscCall(MatDenseGetArray(A, &aa));
2713: for (j = 0; j < A->cmap->n; j++) {
2714: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscRealPart(aa[i + j * mat->lda]);
2715: }
2716: PetscCall(MatDenseRestoreArray(A, &aa));
2717: PetscFunctionReturn(PETSC_SUCCESS);
2718: }
2720: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2721: {
2722: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2723: PetscInt i, j;
2724: PetscScalar *aa;
2726: PetscFunctionBegin;
2727: PetscCall(MatDenseGetArray(A, &aa));
2728: for (j = 0; j < A->cmap->n; j++) {
2729: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscImaginaryPart(aa[i + j * mat->lda]);
2730: }
2731: PetscCall(MatDenseRestoreArray(A, &aa));
2732: PetscFunctionReturn(PETSC_SUCCESS);
2733: }
2735: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2736: {
2737: PetscInt m = A->rmap->n, n = B->cmap->n;
2738: PetscBool cisdense = PETSC_FALSE;
2740: PetscFunctionBegin;
2741: PetscCall(MatSetSizes(C, m, n, m, n));
2742: #if PetscDefined(HAVE_CUDA)
2743: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2744: #endif
2745: #if PetscDefined(HAVE_HIP)
2746: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2747: #endif
2748: if (!cisdense) {
2749: PetscBool flg;
2751: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2752: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2753: /* keep the VecType of A, or else of B since A may be sparse, e.g. VECKOKKOS from MatCreateDenseFromVecType(), but only from one with the same MatType as C, since e.g. the VECCUDA of a MATSEQDENSECUDA is not valid for a MATSEQDENSE C */
2754: PetscCall(PetscObjectTypeCompare((PetscObject)A, ((PetscObject)C)->type_name, &flg));
2755: if (flg) PetscCall(MatSetVecType(C, A->defaultvectype));
2756: else {
2757: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)C)->type_name, &flg));
2758: if (flg) PetscCall(MatSetVecType(C, B->defaultvectype));
2759: }
2760: }
2761: PetscCall(MatSetUp(C));
2762: PetscFunctionReturn(PETSC_SUCCESS);
2763: }
2765: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2766: {
2767: Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data, *c = (Mat_SeqDense *)C->data;
2768: const PetscScalar *av, *bv;
2769: PetscScalar *cv;
2770: PetscBLASInt m, n, k;
2771: PetscScalar _DOne = 1.0, _DZero = 0.0;
2773: PetscFunctionBegin;
2774: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2775: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2776: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2777: if (!m || !n || !k) {
2778: PetscCall(MatZeroEntries(C));
2779: PetscFunctionReturn(PETSC_SUCCESS);
2780: }
2781: PetscCall(MatDenseGetArrayRead(A, &av));
2782: PetscCall(MatDenseGetArrayRead(B, &bv));
2783: PetscCall(MatDenseGetArrayWrite(C, &cv));
2784: PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2785: PetscCall(MatDenseRestoreArrayRead(A, &av));
2786: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2787: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2788: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2789: PetscFunctionReturn(PETSC_SUCCESS);
2790: }
2792: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2793: {
2794: PetscInt m = A->rmap->n, n = B->rmap->n;
2795: PetscBool cisdense = PETSC_FALSE;
2797: PetscFunctionBegin;
2798: PetscCall(MatSetSizes(C, m, n, m, n));
2799: #if PetscDefined(HAVE_CUDA)
2800: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2801: #endif
2802: #if PetscDefined(HAVE_HIP)
2803: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2804: #endif
2805: if (!cisdense) {
2806: PetscBool flg;
2808: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2809: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2810: /* keep the VecType of A, or else of B since A may be sparse, e.g. VECKOKKOS from MatCreateDenseFromVecType(), but only from one with the same MatType as C, since e.g. the VECCUDA of a MATSEQDENSECUDA is not valid for a MATSEQDENSE C */
2811: PetscCall(PetscObjectTypeCompare((PetscObject)A, ((PetscObject)C)->type_name, &flg));
2812: if (flg) PetscCall(MatSetVecType(C, A->defaultvectype));
2813: else {
2814: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)C)->type_name, &flg));
2815: if (flg) PetscCall(MatSetVecType(C, B->defaultvectype));
2816: }
2817: }
2818: PetscCall(MatSetUp(C));
2819: PetscFunctionReturn(PETSC_SUCCESS);
2820: }
2822: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2823: {
2824: Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data, *c = (Mat_SeqDense *)C->data;
2825: const PetscScalar *av, *bv;
2826: PetscScalar *cv;
2827: PetscBLASInt m, n, k;
2828: PetscScalar _DOne = 1.0, _DZero = 0.0;
2830: PetscFunctionBegin;
2831: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2832: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2833: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2834: if (!m || !n || !k) {
2835: PetscCall(MatZeroEntries(C));
2836: PetscFunctionReturn(PETSC_SUCCESS);
2837: }
2838: PetscCall(MatDenseGetArrayRead(A, &av));
2839: PetscCall(MatDenseGetArrayRead(B, &bv));
2840: PetscCall(MatDenseGetArrayWrite(C, &cv));
2841: PetscCallBLAS("BLASgemm", BLASgemm_("N", "T", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2842: PetscCall(MatDenseRestoreArrayRead(A, &av));
2843: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2844: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2845: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2846: PetscFunctionReturn(PETSC_SUCCESS);
2847: }
2849: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2850: {
2851: PetscInt m = A->cmap->n, n = B->cmap->n;
2852: PetscBool cisdense = PETSC_FALSE;
2854: PetscFunctionBegin;
2855: PetscCall(MatSetSizes(C, m, n, m, n));
2856: #if PetscDefined(HAVE_CUDA)
2857: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2858: #endif
2859: #if PetscDefined(HAVE_HIP)
2860: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2861: #endif
2862: if (!cisdense) {
2863: PetscBool flg;
2865: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2866: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2867: /* keep the VecType of A, or else of B since A may be sparse, e.g. VECKOKKOS from MatCreateDenseFromVecType(), but only from one with the same MatType as C, since e.g. the VECCUDA of a MATSEQDENSECUDA is not valid for a MATSEQDENSE C */
2868: PetscCall(PetscObjectTypeCompare((PetscObject)A, ((PetscObject)C)->type_name, &flg));
2869: if (flg) PetscCall(MatSetVecType(C, A->defaultvectype));
2870: else {
2871: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)C)->type_name, &flg));
2872: if (flg) PetscCall(MatSetVecType(C, B->defaultvectype));
2873: }
2874: }
2875: PetscCall(MatSetUp(C));
2876: PetscFunctionReturn(PETSC_SUCCESS);
2877: }
2879: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2880: {
2881: Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data, *c = (Mat_SeqDense *)C->data;
2882: const PetscScalar *av, *bv;
2883: PetscScalar *cv;
2884: PetscBLASInt m, n, k;
2885: PetscScalar _DOne = 1.0, _DZero = 0.0;
2887: PetscFunctionBegin;
2888: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2889: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2890: PetscCall(PetscBLASIntCast(A->rmap->n, &k));
2891: if (!m || !n || !k) {
2892: PetscCall(MatZeroEntries(C));
2893: PetscFunctionReturn(PETSC_SUCCESS);
2894: }
2895: PetscCall(MatDenseGetArrayRead(A, &av));
2896: PetscCall(MatDenseGetArrayRead(B, &bv));
2897: PetscCall(MatDenseGetArrayWrite(C, &cv));
2898: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2899: PetscCall(MatDenseRestoreArrayRead(A, &av));
2900: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2901: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2902: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2903: PetscFunctionReturn(PETSC_SUCCESS);
2904: }
2906: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2907: {
2908: PetscFunctionBegin;
2909: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2910: C->ops->productsymbolic = MatProductSymbolic_AB;
2911: PetscFunctionReturn(PETSC_SUCCESS);
2912: }
2914: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2915: {
2916: PetscFunctionBegin;
2917: C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2918: C->ops->productsymbolic = MatProductSymbolic_AtB;
2919: PetscFunctionReturn(PETSC_SUCCESS);
2920: }
2922: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2923: {
2924: PetscFunctionBegin;
2925: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2926: C->ops->productsymbolic = MatProductSymbolic_ABt;
2927: PetscFunctionReturn(PETSC_SUCCESS);
2928: }
2930: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2931: {
2932: Mat_Product *product = C->product;
2934: PetscFunctionBegin;
2935: switch (product->type) {
2936: case MATPRODUCT_AB:
2937: PetscCall(MatProductSetFromOptions_SeqDense_AB(C));
2938: break;
2939: case MATPRODUCT_AtB:
2940: PetscCall(MatProductSetFromOptions_SeqDense_AtB(C));
2941: break;
2942: case MATPRODUCT_ABt:
2943: PetscCall(MatProductSetFromOptions_SeqDense_ABt(C));
2944: break;
2945: default:
2946: break;
2947: }
2948: PetscFunctionReturn(PETSC_SUCCESS);
2949: }
2951: static PetscErrorCode MatGetRowMax_SeqDense(Mat A, Vec v, PetscInt idx[])
2952: {
2953: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2954: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2955: PetscScalar *x;
2956: const PetscScalar *aa;
2958: PetscFunctionBegin;
2959: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2960: PetscCall(VecGetArray(v, &x));
2961: PetscCall(VecGetLocalSize(v, &p));
2962: PetscCall(MatDenseGetArrayRead(A, &aa));
2963: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2964: for (i = 0; i < m; i++) {
2965: x[i] = aa[i];
2966: if (idx) idx[i] = 0;
2967: for (j = 1; j < n; j++) {
2968: if (PetscRealPart(x[i]) < PetscRealPart(aa[i + a->lda * j])) {
2969: x[i] = aa[i + a->lda * j];
2970: if (idx) idx[i] = j;
2971: }
2972: }
2973: }
2974: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2975: PetscCall(VecRestoreArray(v, &x));
2976: PetscFunctionReturn(PETSC_SUCCESS);
2977: }
2979: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A, Vec v, PetscInt idx[])
2980: {
2981: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2982: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2983: PetscScalar *x;
2984: PetscReal atmp;
2985: const PetscScalar *aa;
2987: PetscFunctionBegin;
2988: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2989: PetscCall(VecGetArray(v, &x));
2990: PetscCall(VecGetLocalSize(v, &p));
2991: PetscCall(MatDenseGetArrayRead(A, &aa));
2992: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2993: for (i = 0; i < m; i++) {
2994: x[i] = PetscAbsScalar(aa[i]);
2995: for (j = 1; j < n; j++) {
2996: atmp = PetscAbsScalar(aa[i + a->lda * j]);
2997: if (PetscAbsScalar(x[i]) < atmp) {
2998: x[i] = atmp;
2999: if (idx) idx[i] = j;
3000: }
3001: }
3002: }
3003: PetscCall(MatDenseRestoreArrayRead(A, &aa));
3004: PetscCall(VecRestoreArray(v, &x));
3005: PetscFunctionReturn(PETSC_SUCCESS);
3006: }
3008: static PetscErrorCode MatGetRowMin_SeqDense(Mat A, Vec v, PetscInt idx[])
3009: {
3010: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3011: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
3012: PetscScalar *x;
3013: const PetscScalar *aa;
3015: PetscFunctionBegin;
3016: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3017: PetscCall(MatDenseGetArrayRead(A, &aa));
3018: PetscCall(VecGetArray(v, &x));
3019: PetscCall(VecGetLocalSize(v, &p));
3020: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
3021: for (i = 0; i < m; i++) {
3022: x[i] = aa[i];
3023: if (idx) idx[i] = 0;
3024: for (j = 1; j < n; j++) {
3025: if (PetscRealPart(x[i]) > PetscRealPart(aa[i + a->lda * j])) {
3026: x[i] = aa[i + a->lda * j];
3027: if (idx) idx[i] = j;
3028: }
3029: }
3030: }
3031: PetscCall(VecRestoreArray(v, &x));
3032: PetscCall(MatDenseRestoreArrayRead(A, &aa));
3033: PetscFunctionReturn(PETSC_SUCCESS);
3034: }
3036: PetscErrorCode MatGetColumnVector_SeqDense(Mat A, Vec v, PetscInt col)
3037: {
3038: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3039: PetscScalar *x;
3040: const PetscScalar *aa;
3042: PetscFunctionBegin;
3043: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3044: PetscCall(MatDenseGetArrayRead(A, &aa));
3045: PetscCall(VecGetArray(v, &x));
3046: PetscCall(PetscArraycpy(x, aa + col * a->lda, A->rmap->n));
3047: PetscCall(VecRestoreArray(v, &x));
3048: PetscCall(MatDenseRestoreArrayRead(A, &aa));
3049: PetscFunctionReturn(PETSC_SUCCESS);
3050: }
3052: PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A, PetscInt type, PetscReal *reductions)
3053: {
3054: PetscInt i, j, m, n;
3055: const PetscScalar *a;
3057: PetscFunctionBegin;
3058: PetscCall(MatGetSize(A, &m, &n));
3059: PetscCall(PetscArrayzero(reductions, n));
3060: PetscCall(MatDenseGetArrayRead(A, &a));
3061: if (type == NORM_2) {
3062: for (i = 0; i < n; i++) {
3063: for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j] * a[j]);
3064: a = PetscSafePointerPlusOffset(a, m);
3065: }
3066: } else if (type == NORM_1) {
3067: for (i = 0; i < n; i++) {
3068: for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j]);
3069: a = PetscSafePointerPlusOffset(a, m);
3070: }
3071: } else if (type == NORM_INFINITY) {
3072: for (i = 0; i < n; i++) {
3073: for (j = 0; j < m; j++) reductions[i] = PetscMax(PetscAbsScalar(a[j]), reductions[i]);
3074: a = PetscSafePointerPlusOffset(a, m);
3075: }
3076: } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
3077: for (i = 0; i < n; i++) {
3078: for (j = 0; j < m; j++) reductions[i] += PetscRealPart(a[j]);
3079: a = PetscSafePointerPlusOffset(a, m);
3080: }
3081: } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3082: for (i = 0; i < n; i++) {
3083: for (j = 0; j < m; j++) reductions[i] += PetscImaginaryPart(a[j]);
3084: a = PetscSafePointerPlusOffset(a, m);
3085: }
3086: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown reduction type");
3087: PetscCall(MatDenseRestoreArrayRead(A, &a));
3088: if (type == NORM_2) {
3089: for (i = 0; i < n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
3090: } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3091: for (i = 0; i < n; i++) reductions[i] /= m;
3092: }
3093: PetscFunctionReturn(PETSC_SUCCESS);
3094: }
3096: PetscErrorCode MatSetRandom_SeqDense(Mat x, PetscRandom rctx)
3097: {
3098: PetscScalar *a;
3099: PetscInt lda, m, n, i, j;
3101: PetscFunctionBegin;
3102: PetscCall(MatGetSize(x, &m, &n));
3103: PetscCall(MatDenseGetLDA(x, &lda));
3104: PetscCall(MatDenseGetArrayWrite(x, &a));
3105: for (j = 0; j < n; j++) {
3106: for (i = 0; i < m; i++) PetscCall(PetscRandomGetValue(rctx, a + j * lda + i));
3107: }
3108: PetscCall(MatDenseRestoreArrayWrite(x, &a));
3109: PetscFunctionReturn(PETSC_SUCCESS);
3110: }
3112: /* vals is not const */
3113: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A, PetscInt col, PetscScalar **vals)
3114: {
3115: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3116: PetscScalar *v;
3118: PetscFunctionBegin;
3119: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3120: PetscCall(MatDenseGetArray(A, &v));
3121: *vals = v + col * a->lda;
3122: PetscCall(MatDenseRestoreArray(A, &v));
3123: PetscFunctionReturn(PETSC_SUCCESS);
3124: }
3126: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A, PetscScalar **vals)
3127: {
3128: PetscFunctionBegin;
3129: if (vals) *vals = NULL; /* user cannot accidentally use the array later */
3130: PetscFunctionReturn(PETSC_SUCCESS);
3131: }
3133: static struct _MatOps MatOps_Values = {MatSetValues_SeqDense,
3134: MatGetRow_SeqDense,
3135: MatRestoreRow_SeqDense,
3136: MatMult_SeqDense,
3137: /* 4*/ MatMultAdd_SeqDense,
3138: MatMultTranspose_SeqDense,
3139: MatMultTransposeAdd_SeqDense,
3140: NULL,
3141: NULL,
3142: NULL,
3143: /* 10*/ NULL,
3144: MatLUFactor_SeqDense,
3145: MatCholeskyFactor_SeqDense,
3146: MatSOR_SeqDense,
3147: MatTranspose_SeqDense,
3148: /* 15*/ MatGetInfo_SeqDense,
3149: MatEqual_SeqDense,
3150: MatGetDiagonal_SeqDense,
3151: MatDiagonalScale_SeqDense,
3152: MatNorm_SeqDense,
3153: /* 20*/ NULL,
3154: NULL,
3155: MatSetOption_SeqDense,
3156: MatZeroEntries_SeqDense,
3157: /* 24*/ MatZeroRows_SeqDense,
3158: NULL,
3159: NULL,
3160: NULL,
3161: NULL,
3162: /* 29*/ MatSetUp_SeqDense,
3163: NULL,
3164: NULL,
3165: NULL,
3166: MatSetInf_SeqDense,
3167: /* 34*/ MatDuplicate_SeqDense,
3168: NULL,
3169: NULL,
3170: NULL,
3171: NULL,
3172: /* 39*/ MatAXPY_SeqDense,
3173: MatCreateSubMatrices_SeqDense,
3174: NULL,
3175: MatGetValues_SeqDense,
3176: MatCopy_SeqDense,
3177: /* 44*/ MatGetRowMax_SeqDense,
3178: MatScale_SeqDense,
3179: MatShift_SeqDense,
3180: NULL,
3181: MatZeroRowsColumns_SeqDense,
3182: /* 49*/ MatSetRandom_SeqDense,
3183: NULL,
3184: NULL,
3185: NULL,
3186: NULL,
3187: /* 54*/ NULL,
3188: NULL,
3189: NULL,
3190: NULL,
3191: NULL,
3192: /* 59*/ MatCreateSubMatrix_SeqDense,
3193: MatDestroy_SeqDense,
3194: MatView_SeqDense,
3195: NULL,
3196: NULL,
3197: /* 64*/ NULL,
3198: NULL,
3199: NULL,
3200: NULL,
3201: MatGetRowMaxAbs_SeqDense,
3202: /* 69*/ NULL,
3203: NULL,
3204: NULL,
3205: NULL,
3206: NULL,
3207: /* 74*/ NULL,
3208: NULL,
3209: NULL,
3210: NULL,
3211: MatLoad_SeqDense,
3212: /* 79*/ MatIsSymmetric_SeqDense,
3213: MatIsHermitian_SeqDense,
3214: NULL,
3215: NULL,
3216: NULL,
3217: /* 84*/ NULL,
3218: MatMatMultNumeric_SeqDense_SeqDense,
3219: NULL,
3220: NULL,
3221: MatMatTransposeMultNumeric_SeqDense_SeqDense,
3222: /* 89*/ NULL,
3223: MatProductSetFromOptions_SeqDense,
3224: NULL,
3225: NULL,
3226: MatConjugate_SeqDense,
3227: /* 94*/ NULL,
3228: NULL,
3229: MatRealPart_SeqDense,
3230: MatImaginaryPart_SeqDense,
3231: NULL,
3232: /* 99*/ NULL,
3233: NULL,
3234: NULL,
3235: MatGetRowMin_SeqDense,
3236: MatGetColumnVector_SeqDense,
3237: /*104*/ NULL,
3238: NULL,
3239: NULL,
3240: NULL,
3241: NULL,
3242: /*109*/ NULL,
3243: NULL,
3244: MatMultHermitianTranspose_SeqDense,
3245: MatMultHermitianTransposeAdd_SeqDense,
3246: NULL,
3247: /*114*/ NULL,
3248: MatGetColumnReductions_SeqDense,
3249: NULL,
3250: NULL,
3251: NULL,
3252: /*119*/ NULL,
3253: MatTransposeMatMultNumeric_SeqDense_SeqDense,
3254: NULL,
3255: NULL,
3256: NULL,
3257: /*124*/ NULL,
3258: NULL,
3259: NULL,
3260: NULL,
3261: NULL,
3262: /*129*/ MatCreateMPIMatConcatenateSeqMat_SeqDense,
3263: NULL,
3264: NULL,
3265: NULL,
3266: NULL,
3267: /*134*/ NULL,
3268: NULL,
3269: NULL,
3270: NULL,
3271: NULL,
3272: /*139*/ NULL,
3273: NULL,
3274: NULL,
3275: NULL,
3276: NULL,
3277: /*144*/ NULL,
3278: NULL,
3279: NULL,
3280: NULL};
3282: /*@
3283: MatCreateSeqDense - Creates a `MATSEQDENSE` that
3284: is stored in column major order (the usual Fortran format).
3286: Collective
3288: Input Parameters:
3289: + comm - MPI communicator, set to `PETSC_COMM_SELF`
3290: . m - number of rows
3291: . n - number of columns
3292: - data - optional location of matrix data in column major order. Use `NULL` for PETSc
3293: to control all matrix memory allocation.
3295: Output Parameter:
3296: . A - the matrix
3298: Level: intermediate
3300: Note:
3301: The data input variable is intended primarily for Fortran programmers
3302: who wish to allocate their own matrix memory space. Most users should
3303: set `data` = `NULL`.
3305: Developer Note:
3306: Many of the matrix operations for this variant use the BLAS and LAPACK routines.
3308: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`
3309: @*/
3310: PetscErrorCode MatCreateSeqDense(MPI_Comm comm, PetscInt m, PetscInt n, PetscScalar data[], Mat *A)
3311: {
3312: PetscFunctionBegin;
3313: PetscCall(MatCreate(comm, A));
3314: PetscCall(MatSetSizes(*A, m, n, m, n));
3315: PetscCall(MatSetType(*A, MATSEQDENSE));
3316: PetscCall(MatSeqDenseSetPreallocation(*A, data));
3317: PetscFunctionReturn(PETSC_SUCCESS);
3318: }
3320: /*@
3321: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements of a `MATSEQDENSE` matrix
3323: Collective
3325: Input Parameters:
3326: + B - the matrix
3327: - data - the array (or `NULL`)
3329: Level: intermediate
3331: Note:
3332: The data input variable is intended primarily for Fortran programmers
3333: who wish to allocate their own matrix memory space. Most users should
3334: need not call this routine.
3336: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`, `MatDenseSetLDA()`
3337: @*/
3338: PetscErrorCode MatSeqDenseSetPreallocation(Mat B, PetscScalar data[])
3339: {
3340: PetscFunctionBegin;
3342: PetscTryMethod(B, "MatSeqDenseSetPreallocation_C", (Mat, PetscScalar[]), (B, data));
3343: PetscFunctionReturn(PETSC_SUCCESS);
3344: }
3346: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B, PetscScalar *data)
3347: {
3348: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
3350: PetscFunctionBegin;
3351: PetscCheck(!b->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3352: B->preallocated = PETSC_TRUE;
3354: PetscCall(PetscLayoutSetUp(B->rmap));
3355: PetscCall(PetscLayoutSetUp(B->cmap));
3357: if (b->lda <= 0) PetscCall(PetscBLASIntCast(B->rmap->n, &b->lda));
3359: if (!data) { /* petsc-allocated storage */
3360: if (!b->user_alloc) PetscCall(PetscFree(b->v));
3361: PetscCall(PetscCalloc1((size_t)b->lda * B->cmap->n, &b->v));
3363: b->user_alloc = PETSC_FALSE;
3364: } else { /* user-allocated storage */
3365: if (!b->user_alloc) PetscCall(PetscFree(b->v));
3366: b->v = data;
3367: b->user_alloc = PETSC_TRUE;
3368: }
3369: B->assembled = PETSC_TRUE;
3370: PetscFunctionReturn(PETSC_SUCCESS);
3371: }
3373: #if PetscDefined(HAVE_ELEMENTAL)
3374: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
3375: {
3376: Mat mat_elemental;
3377: const PetscScalar *array;
3378: PetscScalar *v_colwise;
3379: PetscInt M = A->rmap->N, N = A->cmap->N, i, j, k, *rows, *cols;
3381: PetscFunctionBegin;
3382: PetscCall(PetscMalloc3(M * N, &v_colwise, M, &rows, N, &cols));
3383: PetscCall(MatDenseGetArrayRead(A, &array));
3384: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
3385: k = 0;
3386: for (j = 0; j < N; j++) {
3387: cols[j] = j;
3388: for (i = 0; i < M; i++) v_colwise[j * M + i] = array[k++];
3389: }
3390: for (i = 0; i < M; i++) rows[i] = i;
3391: PetscCall(MatDenseRestoreArrayRead(A, &array));
3393: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
3394: PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
3395: PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
3396: PetscCall(MatSetUp(mat_elemental));
3398: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
3399: PetscCall(MatSetValues(mat_elemental, M, rows, N, cols, v_colwise, ADD_VALUES));
3400: PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
3401: PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
3402: PetscCall(PetscFree3(v_colwise, rows, cols));
3404: if (reuse == MAT_INPLACE_MATRIX) {
3405: PetscCall(MatHeaderReplace(A, &mat_elemental));
3406: } else {
3407: *newmat = mat_elemental;
3408: }
3409: PetscFunctionReturn(PETSC_SUCCESS);
3410: }
3411: #endif
3413: PetscErrorCode MatDenseSetLDA_SeqDense(Mat B, PetscInt lda)
3414: {
3415: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
3416: PetscBool data;
3418: PetscFunctionBegin;
3419: data = (B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE;
3420: PetscCheck(b->user_alloc || !data || b->lda == lda, PETSC_COMM_SELF, PETSC_ERR_ORDER, "LDA cannot be changed after allocation of internal storage");
3421: 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);
3422: PetscCall(PetscBLASIntCast(lda, &b->lda));
3423: PetscFunctionReturn(PETSC_SUCCESS);
3424: }
3426: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
3427: {
3428: PetscFunctionBegin;
3429: PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIDense(comm, inmat, n, scall, outmat));
3430: PetscFunctionReturn(PETSC_SUCCESS);
3431: }
3433: PetscErrorCode MatDenseCreateColumnVec_Private(Mat A, Vec *v)
3434: {
3435: PetscBool isstd, iskok, isdevice;
3436: PetscMPIInt size;
3438: PetscFunctionBegin;
3439: *v = NULL;
3440: PetscCall(PetscStrcmpAny(A->defaultvectype, &isstd, VECSTANDARD, VECSEQ, VECMPI, ""));
3441: PetscCall(PetscStrcmpAny(A->defaultvectype, &iskok, VECKOKKOS, VECSEQKOKKOS, VECMPIKOKKOS, ""));
3442: PetscCall(PetscStrcmpAny(A->defaultvectype, &isdevice, VECCUDA, VECSEQCUDA, VECMPICUDA, VECHIP, VECSEQHIP, VECMPIHIP, ""));
3443: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
3444: PetscCheck(isstd || iskok || isdevice, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not coded for type %s", A->defaultvectype);
3445: if (iskok) {
3446: PetscCheck(PetscDefined(HAVE_KOKKOS_KERNELS), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using KOKKOS kernels support");
3447: #if PetscDefined(HAVE_KOKKOS_KERNELS)
3448: if (size > 1) PetscCall(VecCreateMPIKokkosWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, v));
3449: else PetscCall(VecCreateSeqKokkosWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, v));
3450: #endif
3451: } else {
3452: PetscMemType mtype = PETSC_MEMTYPE_HOST;
3453: const PetscScalar *a = NULL;
3454: PetscBool boundtocpu = PETSC_FALSE;
3456: if (isdevice) {
3457: PetscCall(MatBoundToCPU(A, &boundtocpu));
3458: if (!boundtocpu) {
3459: /* Pass A's device data to avoid an allocation when VecCUPMPlaceArray() is first called. */
3460: PetscCall(MatDenseGetArrayReadAndMemType(A, &a, &mtype));
3461: }
3462: }
3463: if (size > 1) PetscCall(VecCreateMPIWithArrayAndMemType(PetscObjectComm((PetscObject)A), mtype, A->rmap->bs, A->rmap->n, A->rmap->N, a, v));
3464: else PetscCall(VecCreateSeqWithArrayAndMemType(PetscObjectComm((PetscObject)A), mtype, A->rmap->bs, A->rmap->n, a, v));
3465: if (isdevice && !boundtocpu) PetscCall(MatDenseRestoreArrayReadAndMemType(A, &a));
3466: }
3467: PetscFunctionReturn(PETSC_SUCCESS);
3468: }
3470: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3471: {
3472: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3474: PetscFunctionBegin;
3475: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3476: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3477: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3478: a->vecinuse = col + 1;
3479: PetscCall(MatDenseGetArray(A, (PetscScalar **)&a->ptrinuse));
3480: PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda));
3481: *v = a->cvec;
3482: PetscFunctionReturn(PETSC_SUCCESS);
3483: }
3485: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3486: {
3487: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3489: PetscFunctionBegin;
3490: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3491: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3492: VecCheckAssembled(a->cvec);
3493: a->vecinuse = 0;
3494: PetscCall(MatDenseRestoreArray(A, (PetscScalar **)&a->ptrinuse));
3495: PetscCall(VecResetArray(a->cvec));
3496: if (v) *v = NULL;
3497: PetscFunctionReturn(PETSC_SUCCESS);
3498: }
3500: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3501: {
3502: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3504: PetscFunctionBegin;
3505: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3506: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3507: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3508: a->vecinuse = col + 1;
3509: PetscCall(MatDenseGetArrayRead(A, &a->ptrinuse));
3510: PetscCall(VecPlaceArray(a->cvec, PetscSafePointerPlusOffset(a->ptrinuse, (size_t)col * (size_t)a->lda)));
3511: PetscCall(VecLockReadPush(a->cvec));
3512: *v = a->cvec;
3513: PetscFunctionReturn(PETSC_SUCCESS);
3514: }
3516: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3517: {
3518: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3520: PetscFunctionBegin;
3521: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3522: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3523: VecCheckAssembled(a->cvec);
3524: a->vecinuse = 0;
3525: PetscCall(MatDenseRestoreArrayRead(A, &a->ptrinuse));
3526: PetscCall(VecLockReadPop(a->cvec));
3527: PetscCall(VecResetArray(a->cvec));
3528: if (v) *v = NULL;
3529: PetscFunctionReturn(PETSC_SUCCESS);
3530: }
3532: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3533: {
3534: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3536: PetscFunctionBegin;
3537: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3538: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3539: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3540: a->vecinuse = col + 1;
3541: PetscCall(MatDenseGetArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3542: PetscCall(VecPlaceArray(a->cvec, PetscSafePointerPlusOffset(a->ptrinuse, (size_t)col * (size_t)a->lda)));
3543: *v = a->cvec;
3544: PetscFunctionReturn(PETSC_SUCCESS);
3545: }
3547: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3548: {
3549: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3551: PetscFunctionBegin;
3552: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3553: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3554: VecCheckAssembled(a->cvec);
3555: a->vecinuse = 0;
3556: PetscCall(MatDenseRestoreArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3557: PetscCall(VecResetArray(a->cvec));
3558: if (v) *v = NULL;
3559: PetscFunctionReturn(PETSC_SUCCESS);
3560: }
3562: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3563: {
3564: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3565: PetscBool flg;
3567: PetscFunctionBegin;
3568: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3569: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3570: if (a->cmat && (cend - cbegin != a->cmat->cmap->N || rend - rbegin != a->cmat->rmap->N)) PetscCall(MatDestroy(&a->cmat));
3571: if (!a->cmat) {
3572: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), rend - rbegin, PETSC_DECIDE, rend - rbegin, cend - cbegin, PetscSafePointerPlusOffset(a->v, rbegin + (size_t)cbegin * a->lda), &a->cmat));
3573: /* A may be a MATSEQDENSECUDA or MATSEQDENSEHIP bound to the CPU, whose VecType is not valid for the MATSEQDENSE submatrix */
3574: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQDENSE, &flg));
3575: if (flg) PetscCall(MatSetVecType(a->cmat, A->defaultvectype));
3576: } else {
3577: PetscCall(MatDensePlaceArray(a->cmat, PetscSafePointerPlusOffset(a->v, rbegin + (size_t)cbegin * a->lda)));
3578: }
3579: PetscCall(MatDenseSetLDA(a->cmat, a->lda));
3580: a->matinuse = cbegin + 1;
3581: *v = a->cmat;
3582: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
3583: A->offloadmask = PETSC_OFFLOAD_CPU;
3584: #endif
3585: PetscFunctionReturn(PETSC_SUCCESS);
3586: }
3588: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A, Mat *v)
3589: {
3590: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3592: PetscFunctionBegin;
3593: PetscCheck(a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetSubMatrix() first");
3594: PetscCheck(a->cmat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column matrix");
3595: PetscCheck(*v == a->cmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not the matrix obtained from MatDenseGetSubMatrix()");
3596: a->matinuse = 0;
3597: PetscCall(MatDenseResetArray(a->cmat));
3598: *v = NULL;
3599: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
3600: A->offloadmask = PETSC_OFFLOAD_CPU;
3601: #endif
3602: PetscFunctionReturn(PETSC_SUCCESS);
3603: }
3605: static PetscErrorCode MatDenseUpdateColumnLayout_SeqDense(Mat A, PetscLayout clayout)
3606: {
3607: PetscFunctionBegin;
3608: PetscCall(PetscLayoutReference(clayout, &A->cmap));
3609: PetscFunctionReturn(PETSC_SUCCESS);
3610: }
3612: /*MC
3613: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
3615: Options Database Key:
3616: . -mat_type seqdense - sets the matrix type to `MATSEQDENSE` during a call to `MatSetFromOptions()`
3618: Level: beginner
3620: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreateSeqDense()`
3621: M*/
3622: PetscErrorCode MatCreate_SeqDense(Mat B)
3623: {
3624: Mat_SeqDense *b;
3625: PetscMPIInt size;
3627: PetscFunctionBegin;
3628: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
3629: PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1");
3631: PetscCall(PetscNew(&b));
3632: B->data = (void *)b;
3633: B->ops[0] = MatOps_Values;
3635: b->roworiented = PETSC_TRUE;
3637: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactor_C", MatQRFactor_SeqDense));
3638: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetLDA_C", MatDenseGetLDA_SeqDense));
3639: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseSetLDA_C", MatDenseSetLDA_SeqDense));
3640: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArray_C", MatDenseGetArray_SeqDense));
3641: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArray_C", MatDenseRestoreArray_SeqDense));
3642: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDensePlaceArray_C", MatDensePlaceArray_SeqDense));
3643: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseResetArray_C", MatDenseResetArray_SeqDense));
3644: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseReplaceArray_C", MatDenseReplaceArray_SeqDense));
3645: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayRead_C", MatDenseGetArray_SeqDense));
3646: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayRead_C", MatDenseRestoreArray_SeqDense));
3647: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayWrite_C", MatDenseGetArray_SeqDense));
3648: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayWrite_C", MatDenseRestoreArray_SeqDense));
3649: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqaij_C", MatConvert_SeqDense_SeqAIJ));
3650: #if PetscDefined(HAVE_ELEMENTAL)
3651: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_elemental_C", MatConvert_SeqDense_Elemental));
3652: #endif
3653: #if PetscDefined(HAVE_SCALAPACK) && (PetscDefined(USE_REAL_SINGLE) || PetscDefined(USE_REAL_DOUBLE))
3654: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_scalapack_C", MatConvert_Dense_ScaLAPACK));
3655: #endif
3656: #if PetscDefined(HAVE_CUDA)
3657: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensecuda_C", MatConvert_SeqDense_SeqDenseCUDA));
3658: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3659: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdense_C", MatProductSetFromOptions_SeqDense));
3660: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3661: #endif
3662: #if PetscDefined(HAVE_HIP)
3663: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensehip_C", MatConvert_SeqDense_SeqDenseHIP));
3664: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3665: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdense_C", MatProductSetFromOptions_SeqDense));
3666: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3667: #endif
3668: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqDenseSetPreallocation_C", MatSeqDenseSetPreallocation_SeqDense));
3669: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqdense_C", MatProductSetFromOptions_SeqAIJ_SeqDense));
3670: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdense_C", MatProductSetFromOptions_SeqDense));
3671: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));
3672: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqsbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));
3674: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumn_C", MatDenseGetColumn_SeqDense));
3675: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumn_C", MatDenseRestoreColumn_SeqDense));
3676: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVec_C", MatDenseGetColumnVec_SeqDense));
3677: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVec_C", MatDenseRestoreColumnVec_SeqDense));
3678: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecRead_C", MatDenseGetColumnVecRead_SeqDense));
3679: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecRead_C", MatDenseRestoreColumnVecRead_SeqDense));
3680: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecWrite_C", MatDenseGetColumnVecWrite_SeqDense));
3681: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecWrite_C", MatDenseRestoreColumnVecWrite_SeqDense));
3682: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetSubMatrix_C", MatDenseGetSubMatrix_SeqDense));
3683: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreSubMatrix_C", MatDenseRestoreSubMatrix_SeqDense));
3684: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultColumnRange_C", MatMultColumnRange_SeqDense));
3685: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultAddColumnRange_C", MatMultAddColumnRange_SeqDense));
3686: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultHermitianTransposeColumnRange_C", MatMultHermitianTransposeColumnRange_SeqDense));
3687: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultHermitianTransposeAddColumnRange_C", MatMultHermitianTransposeAddColumnRange_SeqDense));
3688: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseUpdateColumnLayout_C", MatDenseUpdateColumnLayout_SeqDense));
3689: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQDENSE));
3690: PetscFunctionReturn(PETSC_SUCCESS);
3691: }
3693: /*@
3694: 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.
3696: Not Collective
3698: Input Parameters:
3699: + A - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3700: - col - column index
3702: Output Parameter:
3703: . vals - pointer to the data
3705: Level: intermediate
3707: Note:
3708: Use `MatDenseGetColumnVec()` to get access to a column of a `MATDENSE` treated as a `Vec`
3710: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreColumn()`, `MatDenseGetColumnVec()`
3711: @*/
3712: PetscErrorCode MatDenseGetColumn(Mat A, PetscInt col, PetscScalar *vals[])
3713: {
3714: PetscFunctionBegin;
3717: PetscAssertPointer(vals, 3);
3718: PetscUseMethod(A, "MatDenseGetColumn_C", (Mat, PetscInt, PetscScalar **), (A, col, vals));
3719: PetscFunctionReturn(PETSC_SUCCESS);
3720: }
3722: /*@
3723: MatDenseRestoreColumn - returns access to a column of a `MATDENSE` matrix which is returned by `MatDenseGetColumn()`.
3725: Not Collective
3727: Input Parameters:
3728: + A - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3729: - vals - pointer to the data (may be `NULL`)
3731: Level: intermediate
3733: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetColumn()`
3734: @*/
3735: PetscErrorCode MatDenseRestoreColumn(Mat A, PetscScalar *vals[])
3736: {
3737: PetscFunctionBegin;
3739: PetscAssertPointer(vals, 2);
3740: PetscUseMethod(A, "MatDenseRestoreColumn_C", (Mat, PetscScalar **), (A, vals));
3741: PetscFunctionReturn(PETSC_SUCCESS);
3742: }
3744: /*@
3745: MatDenseGetColumnVec - Gives read-write access to a column of a `MATDENSE` matrix, represented as a `Vec`.
3747: Collective
3749: Input Parameters:
3750: + A - the `Mat` object
3751: - col - the column index
3753: Output Parameter:
3754: . v - the vector
3756: Level: intermediate
3758: Notes:
3759: The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVec()` when the vector is no longer needed.
3761: Use `MatDenseGetColumnVecRead()` to obtain read-only access or `MatDenseGetColumnVecWrite()` for write-only access.
3763: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`, `MatDenseGetColumn()`
3764: @*/
3765: PetscErrorCode MatDenseGetColumnVec(Mat A, PetscInt col, Vec *v)
3766: {
3767: PetscFunctionBegin;
3771: PetscAssertPointer(v, 3);
3772: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3773: 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);
3774: PetscUseMethod(A, "MatDenseGetColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3775: PetscFunctionReturn(PETSC_SUCCESS);
3776: }
3778: /*@
3779: MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVec()`.
3781: Collective
3783: Input Parameters:
3784: + A - the `Mat` object
3785: . col - the column index
3786: - v - the `Vec` object (may be `NULL`)
3788: Level: intermediate
3790: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3791: @*/
3792: PetscErrorCode MatDenseRestoreColumnVec(Mat A, PetscInt col, Vec *v)
3793: {
3794: PetscFunctionBegin;
3799: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3800: 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);
3801: PetscUseMethod(A, "MatDenseRestoreColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3802: PetscFunctionReturn(PETSC_SUCCESS);
3803: }
3805: /*@
3806: MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a `Vec`.
3808: Collective
3810: Input Parameters:
3811: + A - the `Mat` object
3812: - col - the column index
3814: Output Parameter:
3815: . v - the vector
3817: Level: intermediate
3819: Notes:
3820: The vector is owned by PETSc and users cannot modify it.
3822: Users need to call `MatDenseRestoreColumnVecRead()` when the vector is no longer needed.
3824: Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecWrite()` for write-only access.
3826: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3827: @*/
3828: PetscErrorCode MatDenseGetColumnVecRead(Mat A, PetscInt col, Vec *v)
3829: {
3830: PetscFunctionBegin;
3834: PetscAssertPointer(v, 3);
3835: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3836: 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);
3837: PetscUseMethod(A, "MatDenseGetColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3838: PetscFunctionReturn(PETSC_SUCCESS);
3839: }
3841: /*@
3842: MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVecRead()`.
3844: Collective
3846: Input Parameters:
3847: + A - the `Mat` object
3848: . col - the column index
3849: - v - the `Vec` object (may be `NULL`)
3851: Level: intermediate
3853: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecWrite()`
3854: @*/
3855: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A, PetscInt col, Vec *v)
3856: {
3857: PetscFunctionBegin;
3862: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3863: 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);
3864: PetscUseMethod(A, "MatDenseRestoreColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3865: PetscFunctionReturn(PETSC_SUCCESS);
3866: }
3868: /*@
3869: MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a `Vec`.
3871: Collective
3873: Input Parameters:
3874: + A - the `Mat` object
3875: - col - the column index
3877: Output Parameter:
3878: . v - the vector
3880: Level: intermediate
3882: Notes:
3883: The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVecWrite()` when the vector is no longer needed.
3885: Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecRead()` for read-only access.
3887: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3888: @*/
3889: PetscErrorCode MatDenseGetColumnVecWrite(Mat A, PetscInt col, Vec *v)
3890: {
3891: PetscFunctionBegin;
3895: PetscAssertPointer(v, 3);
3896: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3897: 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);
3898: PetscUseMethod(A, "MatDenseGetColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3899: PetscFunctionReturn(PETSC_SUCCESS);
3900: }
3902: /*@
3903: MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVecWrite()`.
3905: Collective
3907: Input Parameters:
3908: + A - the `Mat` object
3909: . col - the column index
3910: - v - the `Vec` object (may be `NULL`)
3912: Level: intermediate
3914: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`
3915: @*/
3916: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A, PetscInt col, Vec *v)
3917: {
3918: PetscFunctionBegin;
3923: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3924: 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);
3925: PetscUseMethod(A, "MatDenseRestoreColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3926: PetscFunctionReturn(PETSC_SUCCESS);
3927: }
3929: /*@
3930: MatDenseGetSubMatrix - Gives access to a block of rows and columns of a dense matrix, represented as a `Mat`.
3932: Collective
3934: Input Parameters:
3935: + A - the `Mat` object
3936: . rbegin - the first global row index in the block (if `PETSC_DECIDE`, is 0)
3937: . rend - the global row index past the last one in the block (if `PETSC_DECIDE`, is `M`)
3938: . cbegin - the first global column index in the block (if `PETSC_DECIDE`, is 0)
3939: - cend - the global column index past the last one in the block (if `PETSC_DECIDE`, is `N`)
3941: Output Parameter:
3942: . v - the matrix
3944: Level: intermediate
3946: Notes:
3947: The matrix is owned by PETSc. Users need to call `MatDenseRestoreSubMatrix()` when the matrix is no longer needed.
3949: The output matrix is not redistributed by PETSc, so depending on the values of `rbegin` and `rend`, some processes may have no local rows.
3951: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreSubMatrix()`
3952: @*/
3953: PetscErrorCode MatDenseGetSubMatrix(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3954: {
3955: PetscFunctionBegin;
3962: PetscAssertPointer(v, 6);
3963: if (rbegin == PETSC_DECIDE) rbegin = 0;
3964: if (rend == PETSC_DECIDE) rend = A->rmap->N;
3965: if (cbegin == PETSC_DECIDE) cbegin = 0;
3966: if (cend == PETSC_DECIDE) cend = A->cmap->N;
3967: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3968: 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);
3969: 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);
3970: 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);
3971: 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);
3972: PetscUseMethod(A, "MatDenseGetSubMatrix_C", (Mat, PetscInt, PetscInt, PetscInt, PetscInt, Mat *), (A, rbegin, rend, cbegin, cend, v));
3973: PetscFunctionReturn(PETSC_SUCCESS);
3974: }
3976: /*@
3977: MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from `MatDenseGetSubMatrix()`.
3979: Collective
3981: Input Parameters:
3982: + A - the `Mat` object
3983: - v - the `Mat` object (cannot be `NULL`)
3985: Level: intermediate
3987: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseGetSubMatrix()`
3988: @*/
3989: PetscErrorCode MatDenseRestoreSubMatrix(Mat A, Mat *v)
3990: {
3991: PetscFunctionBegin;
3994: PetscAssertPointer(v, 2);
3996: PetscUseMethod(A, "MatDenseRestoreSubMatrix_C", (Mat, Mat *), (A, v));
3997: PetscFunctionReturn(PETSC_SUCCESS);
3998: }
4000: /*@
4001: MatDenseUpdateColumnLayout - Update the column layout of the dense matrix.
4003: Collective
4005: Input Parameters:
4006: + A - the `Mat` object
4007: - clayout - the `PetscLayout` object (cannot be `NULL`)
4009: Level: advanced
4011: Notes:
4012: Because a dense matrix's storage is independent of its column layout, this routine can update the layout without modifying the underlying storage.
4014: It can be useful when users want to apply the operator, with `MatMult()` on a right vector with a different layout.
4016: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `PetscLayout`, `MatMult()`, `MatMultAdd()`
4017: @*/
4018: PetscErrorCode MatDenseUpdateColumnLayout(Mat A, PetscLayout clayout)
4019: {
4020: PetscMPIInt flag;
4021: MPI_Comm lcomm;
4023: PetscFunctionBegin;
4026: PetscCall(PetscLayoutGetComm(clayout, &lcomm));
4027: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)A), lcomm, &flag));
4028: PetscCheck(flag == MPI_CONGRUENT || flag == MPI_IDENT, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMECOMM, "Different communicators in the two objects: flag %d", flag);
4029: PetscCheck(A->cmap->N == clayout->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_SIZ, "Mat global dim %" PetscInt_FMT " does not match layout global dim %" PetscInt_FMT, A->cmap->N, clayout->N);
4030: PetscUseMethod(A, "MatDenseUpdateColumnLayout_C", (Mat, PetscLayout), (A, clayout));
4031: PetscFunctionReturn(PETSC_SUCCESS);
4032: }
4034: #include <petscblaslapack.h>
4035: #include <petsc/private/kernels/blockinvert.h>
4037: /*@
4038: MatSeqDenseInvert - Invert a small `MATSEQDENSE` matrix in place using a hard-coded kernel.
4040: Not Collective
4042: Input Parameter:
4043: . A - the `MATSEQDENSE` matrix
4045: Level: developer
4047: Note:
4048: Intended for small blocks; specialized kernels are used for sizes up to 7 and `LAPACK` is used for larger sizes.
4049: If the matrix is singular and `MatSetErrorIfFailure()` was called an error will be immediately generated, otherwise the factor error type
4050: in the matrix, which can be obtained with `MatFactorGetError()`, is set to `MAT_FACTOR_NUMERIC_ZEROPIVOT`.
4052: .seealso: `Mat`, `MATSEQDENSE`, `MatInvertBlockDiagonal()`, `MatLUFactor()`, `MatSetErrorIfFailure()`, `MatFactorGetError()`
4053: @*/
4054: PetscErrorCode MatSeqDenseInvert(Mat A)
4055: {
4056: PetscInt m;
4057: const PetscReal shift = 0.0;
4058: PetscBool allowzeropivot, zeropivotdetected = PETSC_FALSE;
4059: PetscScalar *values;
4061: PetscFunctionBegin;
4063: PetscCall(MatDenseGetArray(A, &values));
4064: PetscCall(MatGetLocalSize(A, &m, NULL));
4065: allowzeropivot = PetscNot(A->erroriffailure);
4066: /* factor and invert each block */
4067: switch (m) {
4068: case 1:
4069: values[0] = (PetscScalar)1.0 / (values[0] + shift);
4070: break;
4071: case 2:
4072: PetscCall(PetscKernel_A_gets_inverse_A_2(values, shift, allowzeropivot, &zeropivotdetected));
4073: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4074: break;
4075: case 3:
4076: PetscCall(PetscKernel_A_gets_inverse_A_3(values, shift, allowzeropivot, &zeropivotdetected));
4077: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4078: break;
4079: case 4:
4080: PetscCall(PetscKernel_A_gets_inverse_A_4(values, shift, allowzeropivot, &zeropivotdetected));
4081: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4082: break;
4083: case 5: {
4084: PetscScalar work[25];
4085: PetscInt ipvt[5];
4087: PetscCall(PetscKernel_A_gets_inverse_A_5(values, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
4088: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4089: } break;
4090: case 6:
4091: PetscCall(PetscKernel_A_gets_inverse_A_6(values, shift, allowzeropivot, &zeropivotdetected));
4092: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4093: break;
4094: case 7:
4095: PetscCall(PetscKernel_A_gets_inverse_A_7(values, shift, allowzeropivot, &zeropivotdetected));
4096: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4097: break;
4098: default: {
4099: PetscInt *v_pivots, *IJ, j;
4100: PetscScalar *v_work;
4102: PetscCall(PetscMalloc3(m, &v_work, m, &v_pivots, m, &IJ));
4103: for (j = 0; j < m; j++) IJ[j] = j;
4104: PetscCall(PetscKernel_A_gets_inverse_A(m, values, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
4105: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4106: PetscCall(PetscFree3(v_work, v_pivots, IJ));
4107: }
4108: }
4109: PetscCall(MatDenseRestoreArray(A, &values));
4110: PetscFunctionReturn(PETSC_SUCCESS);
4111: }
4113: /*@
4114: MatDenseReplaceArrayWithMemType - Allows one to replace the array in a `MATDENSE`, `MATDENSECUDA`, or `MATDENSEHIP`
4115: with an array provided by the user and a matching `PetscMemType`. This is useful to avoid copying an array into a matrix.
4117: Not Collective
4119: Input Parameters:
4120: + mat - the matrix
4121: . mtype - the `PetscMemType` of the array
4122: - array - the array in column major order
4124: Level: developer
4126: Note:
4127: Adding `const` to `array` was an oversight, see notes in `VecPlaceArray()`.
4129: This permanently replaces the GPU array and frees the memory associated with the old GPU
4130: array. The memory passed in CANNOT be freed by the user. It will be freed when the matrix is
4131: destroyed. The array should respect the matrix leading dimension.
4133: .seealso: `MatDenseReplaceArray()`, `MatDenseCUDAReplaceArray()`, `MatDenseHIPReplaceArray()`
4134: @*/
4135: PetscErrorCode MatDenseReplaceArrayWithMemType(Mat mat, PetscMemType mtype, const PetscScalar array[])
4136: {
4137: const char *type = PetscMemTypeToString(mtype) + 14; /* skip "PETSC_MEMTYPE_" */
4138: char buffer[256];
4140: PetscFunctionBegin;
4142: PetscAssertPointer(array, 3);
4143: PetscCall(PetscSNPrintf(buffer, sizeof(buffer), "MatDense%sReplaceArray_C", PetscMemTypeHost(mtype) ? "" : type));
4144: PetscUseMethod(mat, buffer, (Mat, const PetscScalar[]), (mat, array));
4145: PetscFunctionReturn(PETSC_SUCCESS);
4146: }