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