Actual source code: dense.c
1: /*
2: Defines the basic matrix operations for sequential dense.
3: Portions of this code are under:
4: Copyright (c) 2022 Advanced Micro Devices, Inc. All rights reserved.
5: */
7: #include <../src/mat/impls/dense/seq/dense.h>
8: #include <../src/mat/impls/dense/mpi/mpidense.h>
9: #include <petscblaslapack.h>
10: #include <../src/mat/impls/aij/seq/aij.h>
11: #include <petsc/private/vecimpl.h>
13: PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian)
14: {
15: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
16: PetscInt j, k, n = A->rmap->n;
17: PetscScalar *v;
19: PetscFunctionBegin;
20: PetscCheck(A->rmap->n == A->cmap->n, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot symmetrize a rectangular matrix");
21: PetscCall(MatDenseGetArray(A, &v));
22: if (!hermitian) {
23: for (k = 0; k < n; k++) {
24: for (j = k; j < n; j++) v[j * mat->lda + k] = v[k * mat->lda + j];
25: }
26: } else {
27: for (k = 0; k < n; k++) {
28: for (j = k; j < n; j++) v[j * mat->lda + k] = PetscConj(v[k * mat->lda + j]);
29: }
30: }
31: PetscCall(MatDenseRestoreArray(A, &v));
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A)
36: {
37: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
38: PetscBLASInt info, n;
40: PetscFunctionBegin;
41: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
42: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
43: if (A->factortype == MAT_FACTOR_LU) {
44: PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present");
45: if (!mat->fwork) {
46: mat->lfwork = n;
47: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
48: }
49: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
50: PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
51: PetscCall(PetscFPTrapPop());
52: PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
53: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
54: if (A->spd == PETSC_BOOL3_TRUE) {
55: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
56: PetscCallBLAS("LAPACKpotri", LAPACKpotri_("L", &n, mat->v, &mat->lda, &info));
57: PetscCall(PetscFPTrapPop());
58: PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_TRUE));
59: #if defined(PETSC_USE_COMPLEX)
60: } else if (A->hermitian == PETSC_BOOL3_TRUE) {
61: PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present");
62: PetscCheck(mat->fwork, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Fwork not present");
63: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
64: PetscCallBLAS("LAPACKhetri", LAPACKhetri_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &info));
65: PetscCall(PetscFPTrapPop());
66: PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_TRUE));
67: #endif
68: } else { /* symmetric case */
69: PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present");
70: PetscCheck(mat->fwork, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Fwork not present");
71: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
72: PetscCallBLAS("LAPACKsytri", LAPACKsytri_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &info));
73: PetscCall(PetscFPTrapPop());
74: PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_FALSE));
75: }
76: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad Inversion: zero pivot in row %" PetscBLASInt_FMT, info - 1);
77: PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
78: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix must be factored to solve");
80: A->ops->solve = NULL;
81: A->ops->matsolve = NULL;
82: A->ops->solvetranspose = NULL;
83: A->ops->matsolvetranspose = NULL;
84: A->ops->solveadd = NULL;
85: A->ops->solvetransposeadd = NULL;
86: A->factortype = MAT_FACTOR_NONE;
87: PetscCall(PetscFree(A->solvertype));
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
91: static PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
92: {
93: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
94: PetscInt m = l->lda, n = A->cmap->n, r = A->rmap->n, i, j;
95: PetscScalar *slot, *bb, *v;
96: const PetscScalar *xx;
98: PetscFunctionBegin;
99: if (PetscDefined(USE_DEBUG)) {
100: for (i = 0; i < N; i++) {
101: PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
102: PetscCheck(rows[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " requested to be zeroed greater than or equal number of rows %" PetscInt_FMT, rows[i], A->rmap->n);
103: PetscCheck(rows[i] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Col %" PetscInt_FMT " requested to be zeroed greater than or equal number of cols %" PetscInt_FMT, rows[i], A->cmap->n);
104: }
105: }
106: if (!N) PetscFunctionReturn(PETSC_SUCCESS);
108: /* fix right-hand side if needed */
109: if (x && b) {
110: Vec xt;
112: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
113: PetscCall(VecDuplicate(x, &xt));
114: PetscCall(VecCopy(x, xt));
115: PetscCall(VecScale(xt, -1.0));
116: PetscCall(MatMultAdd(A, xt, b, b));
117: PetscCall(VecDestroy(&xt));
118: PetscCall(VecGetArrayRead(x, &xx));
119: PetscCall(VecGetArray(b, &bb));
120: for (i = 0; i < N; i++) bb[rows[i]] = diag * xx[rows[i]];
121: PetscCall(VecRestoreArrayRead(x, &xx));
122: PetscCall(VecRestoreArray(b, &bb));
123: }
125: PetscCall(MatDenseGetArray(A, &v));
126: for (i = 0; i < N; i++) {
127: slot = v + rows[i] * m;
128: PetscCall(PetscArrayzero(slot, r));
129: }
130: for (i = 0; i < N; i++) {
131: slot = v + rows[i];
132: for (j = 0; j < n; j++) {
133: *slot = 0.0;
134: slot += m;
135: }
136: }
137: if (diag != 0.0) {
138: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
139: for (i = 0; i < N; i++) {
140: slot = v + (m + 1) * rows[i];
141: *slot = diag;
142: }
143: }
144: PetscCall(MatDenseRestoreArray(A, &v));
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
149: {
150: Mat B = NULL;
151: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
152: Mat_SeqDense *b;
153: PetscInt *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->N, i;
154: const MatScalar *av;
155: PetscBool isseqdense;
157: PetscFunctionBegin;
158: if (reuse == MAT_REUSE_MATRIX) {
159: PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATSEQDENSE, &isseqdense));
160: PetscCheck(isseqdense, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_USER, "Cannot reuse matrix of type %s", ((PetscObject)*newmat)->type_name);
161: }
162: if (reuse != MAT_REUSE_MATRIX) {
163: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
164: PetscCall(MatSetSizes(B, m, n, m, n));
165: PetscCall(MatSetType(B, MATSEQDENSE));
166: PetscCall(MatSeqDenseSetPreallocation(B, NULL));
167: b = (Mat_SeqDense *)B->data;
168: } else {
169: b = (Mat_SeqDense *)((*newmat)->data);
170: for (i = 0; i < n; i++) PetscCall(PetscArrayzero(b->v + i * b->lda, m));
171: }
172: PetscCall(MatSeqAIJGetArrayRead(A, &av));
173: for (i = 0; i < m; i++) {
174: PetscInt j;
175: for (j = 0; j < ai[1] - ai[0]; j++) {
176: b->v[*aj * b->lda + i] = *av;
177: aj++;
178: av++;
179: }
180: ai++;
181: }
182: PetscCall(MatSeqAIJRestoreArrayRead(A, &av));
184: if (reuse == MAT_INPLACE_MATRIX) {
185: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
186: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
187: PetscCall(MatHeaderReplace(A, &B));
188: } else {
189: if (B) *newmat = B;
190: PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
191: PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
192: }
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
197: {
198: Mat B = NULL;
199: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
200: PetscInt i, j;
201: PetscInt *rows, *nnz;
202: MatScalar *aa = a->v, *vals;
204: PetscFunctionBegin;
205: PetscCall(PetscCalloc3(A->rmap->n, &rows, A->rmap->n, &nnz, A->rmap->n, &vals));
206: if (reuse != MAT_REUSE_MATRIX) {
207: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
208: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
209: PetscCall(MatSetType(B, MATSEQAIJ));
210: for (j = 0; j < A->cmap->n; j++) {
211: for (i = 0; i < A->rmap->n; i++)
212: if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) ++nnz[i];
213: aa += a->lda;
214: }
215: PetscCall(MatSeqAIJSetPreallocation(B, PETSC_DETERMINE, nnz));
216: } else B = *newmat;
217: aa = a->v;
218: for (j = 0; j < A->cmap->n; j++) {
219: PetscInt numRows = 0;
220: for (i = 0; i < A->rmap->n; i++)
221: if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) {
222: rows[numRows] = i;
223: vals[numRows++] = aa[i];
224: }
225: PetscCall(MatSetValues(B, numRows, rows, 1, &j, vals, INSERT_VALUES));
226: aa += a->lda;
227: }
228: PetscCall(PetscFree3(rows, nnz, vals));
229: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
230: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
232: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
233: else if (reuse != MAT_REUSE_MATRIX) *newmat = B;
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: PetscErrorCode MatAXPY_SeqDense(Mat Y, PetscScalar alpha, Mat X, MatStructure str)
238: {
239: Mat_SeqDense *x = (Mat_SeqDense *)X->data, *y = (Mat_SeqDense *)Y->data;
240: const PetscScalar *xv;
241: PetscScalar *yv;
242: PetscBLASInt N, m, ldax = 0, lday = 0, one = 1;
244: PetscFunctionBegin;
245: PetscCall(MatDenseGetArrayRead(X, &xv));
246: PetscCall(MatDenseGetArray(Y, &yv));
247: PetscCall(PetscBLASIntCast(X->rmap->n * X->cmap->n, &N));
248: PetscCall(PetscBLASIntCast(X->rmap->n, &m));
249: PetscCall(PetscBLASIntCast(x->lda, &ldax));
250: PetscCall(PetscBLASIntCast(y->lda, &lday));
251: if (ldax > m || lday > m) {
252: for (PetscInt j = 0; j < X->cmap->n; j++) PetscCallBLAS("BLASaxpy", BLASaxpy_(&m, &alpha, PetscSafePointerPlusOffset(xv, j * ldax), &one, PetscSafePointerPlusOffset(yv, j * lday), &one));
253: } else {
254: PetscCallBLAS("BLASaxpy", BLASaxpy_(&N, &alpha, xv, &one, yv, &one));
255: }
256: PetscCall(MatDenseRestoreArrayRead(X, &xv));
257: PetscCall(MatDenseRestoreArray(Y, &yv));
258: PetscCall(PetscLogFlops(PetscMax(2.0 * N - 1, 0)));
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: static PetscErrorCode MatGetInfo_SeqDense(Mat A, MatInfoType flag, MatInfo *info)
263: {
264: PetscLogDouble N = A->rmap->n * A->cmap->n;
266: PetscFunctionBegin;
267: info->block_size = 1.0;
268: info->nz_allocated = N;
269: info->nz_used = N;
270: info->nz_unneeded = 0;
271: info->assemblies = A->num_ass;
272: info->mallocs = 0;
273: info->memory = 0; /* REVIEW ME */
274: info->fill_ratio_given = 0;
275: info->fill_ratio_needed = 0;
276: info->factor_mallocs = 0;
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: PetscErrorCode MatScale_SeqDense(Mat A, PetscScalar alpha)
281: {
282: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
283: PetscScalar *v;
284: PetscBLASInt one = 1, j, nz, lda = 0;
286: PetscFunctionBegin;
287: PetscCall(MatDenseGetArray(A, &v));
288: PetscCall(PetscBLASIntCast(a->lda, &lda));
289: if (lda > A->rmap->n) {
290: PetscCall(PetscBLASIntCast(A->rmap->n, &nz));
291: for (j = 0; j < A->cmap->n; j++) PetscCallBLAS("BLASscal", BLASscal_(&nz, &alpha, v + j * lda, &one));
292: } else {
293: PetscCall(PetscBLASIntCast(A->rmap->n * A->cmap->n, &nz));
294: PetscCallBLAS("BLASscal", BLASscal_(&nz, &alpha, v, &one));
295: }
296: PetscCall(PetscLogFlops(A->rmap->n * A->cmap->n));
297: PetscCall(MatDenseRestoreArray(A, &v));
298: PetscFunctionReturn(PETSC_SUCCESS);
299: }
301: PetscErrorCode MatShift_SeqDense(Mat A, PetscScalar alpha)
302: {
303: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
304: PetscScalar *v;
305: PetscInt j, k;
307: PetscFunctionBegin;
308: PetscCall(MatDenseGetArray(A, &v));
309: k = PetscMin(A->rmap->n, A->cmap->n);
310: for (j = 0; j < k; j++) v[j + j * a->lda] += alpha;
311: PetscCall(PetscLogFlops(k));
312: PetscCall(MatDenseRestoreArray(A, &v));
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: static PetscErrorCode MatIsHermitian_SeqDense(Mat A, PetscReal rtol, PetscBool *fl)
317: {
318: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
319: PetscInt i, j, m = A->rmap->n, N = a->lda;
320: const PetscScalar *v;
322: PetscFunctionBegin;
323: *fl = PETSC_FALSE;
324: if (A->rmap->n != A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
325: PetscCall(MatDenseGetArrayRead(A, &v));
326: for (i = 0; i < m; i++) {
327: for (j = i; j < m; j++) {
328: if (PetscAbsScalar(v[i + j * N] - PetscConj(v[j + i * N])) > rtol) goto restore;
329: }
330: }
331: *fl = PETSC_TRUE;
332: restore:
333: PetscCall(MatDenseRestoreArrayRead(A, &v));
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: static PetscErrorCode MatIsSymmetric_SeqDense(Mat A, PetscReal rtol, PetscBool *fl)
338: {
339: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
340: PetscInt i, j, m = A->rmap->n, N = a->lda;
341: const PetscScalar *v;
343: PetscFunctionBegin;
344: *fl = PETSC_FALSE;
345: if (A->rmap->n != A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
346: PetscCall(MatDenseGetArrayRead(A, &v));
347: for (i = 0; i < m; i++) {
348: for (j = i; j < m; j++) {
349: if (PetscAbsScalar(v[i + j * N] - v[j + i * N]) > rtol) goto restore;
350: }
351: }
352: *fl = PETSC_TRUE;
353: restore:
354: PetscCall(MatDenseRestoreArrayRead(A, &v));
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
358: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi, Mat A, MatDuplicateOption cpvalues)
359: {
360: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
361: PetscInt lda = mat->lda, j, m, nlda = lda;
362: PetscBool isdensecpu;
364: PetscFunctionBegin;
365: PetscCall(PetscLayoutReference(A->rmap, &newi->rmap));
366: PetscCall(PetscLayoutReference(A->cmap, &newi->cmap));
367: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */
368: PetscCall(MatDenseSetLDA(newi, lda));
369: }
370: PetscCall(PetscObjectTypeCompare((PetscObject)newi, MATSEQDENSE, &isdensecpu));
371: if (isdensecpu) PetscCall(MatSeqDenseSetPreallocation(newi, NULL));
372: if (cpvalues == MAT_COPY_VALUES) {
373: const PetscScalar *av;
374: PetscScalar *v;
376: PetscCall(MatDenseGetArrayRead(A, &av));
377: PetscCall(MatDenseGetArrayWrite(newi, &v));
378: PetscCall(MatDenseGetLDA(newi, &nlda));
379: m = A->rmap->n;
380: if (lda > m || nlda > m) {
381: for (j = 0; j < A->cmap->n; j++) PetscCall(PetscArraycpy(PetscSafePointerPlusOffset(v, j * nlda), PetscSafePointerPlusOffset(av, j * lda), m));
382: } else {
383: PetscCall(PetscArraycpy(v, av, A->rmap->n * A->cmap->n));
384: }
385: PetscCall(MatDenseRestoreArrayWrite(newi, &v));
386: PetscCall(MatDenseRestoreArrayRead(A, &av));
387: PetscCall(MatPropagateSymmetryOptions(A, newi));
388: }
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: PetscErrorCode MatDuplicate_SeqDense(Mat A, MatDuplicateOption cpvalues, Mat *newmat)
393: {
394: PetscFunctionBegin;
395: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), newmat));
396: PetscCall(MatSetSizes(*newmat, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
397: PetscCall(MatSetType(*newmat, ((PetscObject)A)->type_name));
398: PetscCall(MatDuplicateNoCreate_SeqDense(*newmat, A, cpvalues));
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: static PetscErrorCode MatSolve_SeqDense_Internal_LU(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
403: {
404: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
405: PetscBLASInt info;
407: PetscFunctionBegin;
408: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
409: PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_(T ? "T" : "N", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
410: PetscCall(PetscFPTrapPop());
411: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "GETRS - Bad solve %" PetscBLASInt_FMT, info);
412: PetscCall(PetscLogFlops(nrhs * (2.0 * m * m - m)));
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: static PetscErrorCode MatSolve_SeqDense_Internal_Cholesky(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
417: {
418: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
419: PetscBLASInt info;
421: PetscFunctionBegin;
422: if (A->spd == PETSC_BOOL3_TRUE) {
423: if (PetscDefined(USE_COMPLEX) && T) PetscCall(MatConjugate_SeqDense(A));
424: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
425: PetscCallBLAS("LAPACKpotrs", LAPACKpotrs_("L", &m, &nrhs, mat->v, &mat->lda, x, &m, &info));
426: PetscCall(PetscFPTrapPop());
427: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "POTRS Bad solve %" PetscBLASInt_FMT, info);
428: if (PetscDefined(USE_COMPLEX) && T) PetscCall(MatConjugate_SeqDense(A));
429: #if defined(PETSC_USE_COMPLEX)
430: } else if (A->hermitian == PETSC_BOOL3_TRUE) {
431: if (T) PetscCall(MatConjugate_SeqDense(A));
432: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
433: PetscCallBLAS("LAPACKhetrs", LAPACKhetrs_("L", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
434: PetscCall(PetscFPTrapPop());
435: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "HETRS Bad solve %" PetscBLASInt_FMT, info);
436: if (T) PetscCall(MatConjugate_SeqDense(A));
437: #endif
438: } else { /* symmetric case */
439: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
440: PetscCallBLAS("LAPACKsytrs", LAPACKsytrs_("L", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
441: PetscCall(PetscFPTrapPop());
442: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "SYTRS Bad solve %" PetscBLASInt_FMT, info);
443: }
444: PetscCall(PetscLogFlops(nrhs * (2.0 * m * m - m)));
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: static PetscErrorCode MatSolve_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
449: {
450: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
451: PetscBLASInt info;
452: char trans;
454: PetscFunctionBegin;
455: if (PetscDefined(USE_COMPLEX)) {
456: trans = 'C';
457: } else {
458: trans = 'T';
459: }
460: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
461: { /* lwork depends on the number of right-hand sides */
462: PetscBLASInt nlfwork, lfwork = -1;
463: PetscScalar fwork;
465: PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", &trans, &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, &fwork, &lfwork, &info));
466: nlfwork = (PetscBLASInt)PetscRealPart(fwork);
467: if (nlfwork > mat->lfwork) {
468: mat->lfwork = nlfwork;
469: PetscCall(PetscFree(mat->fwork));
470: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
471: }
472: }
473: PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", &trans, &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, mat->fwork, &mat->lfwork, &info));
474: PetscCall(PetscFPTrapPop());
475: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "ORMQR - Bad orthogonal transform %" PetscBLASInt_FMT, info);
476: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
477: PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "N", "N", &mat->rank, &nrhs, mat->v, &mat->lda, x, &ldx, &info));
478: PetscCall(PetscFPTrapPop());
479: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "TRTRS - Bad triangular solve %" PetscBLASInt_FMT, info);
480: for (PetscInt j = 0; j < nrhs; j++) {
481: for (PetscInt i = mat->rank; i < k; i++) x[j * ldx + i] = 0.;
482: }
483: PetscCall(PetscLogFlops(nrhs * (4.0 * m * mat->rank - PetscSqr(mat->rank))));
484: PetscFunctionReturn(PETSC_SUCCESS);
485: }
487: static PetscErrorCode MatSolveTranspose_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
488: {
489: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
490: PetscBLASInt info;
492: PetscFunctionBegin;
493: if (A->rmap->n == A->cmap->n && mat->rank == A->rmap->n) {
494: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
495: PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "T", "N", &m, &nrhs, mat->v, &mat->lda, x, &ldx, &info));
496: PetscCall(PetscFPTrapPop());
497: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "TRTRS - Bad triangular solve %" PetscBLASInt_FMT, info);
498: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate_SeqDense(A));
499: { /* lwork depends on the number of right-hand sides */
500: PetscBLASInt nlfwork, lfwork = -1;
501: PetscScalar fwork;
503: PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", "N", &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, &fwork, &lfwork, &info));
504: nlfwork = (PetscBLASInt)PetscRealPart(fwork);
505: if (nlfwork > mat->lfwork) {
506: mat->lfwork = nlfwork;
507: PetscCall(PetscFree(mat->fwork));
508: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
509: }
510: }
511: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
512: PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", "N", &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, mat->fwork, &mat->lfwork, &info));
513: PetscCall(PetscFPTrapPop());
514: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "ORMQR - Bad orthogonal transform %" PetscBLASInt_FMT, info);
515: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate_SeqDense(A));
516: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "QR factored matrix cannot be used for transpose solve");
517: PetscCall(PetscLogFlops(nrhs * (4.0 * m * mat->rank - PetscSqr(mat->rank))));
518: PetscFunctionReturn(PETSC_SUCCESS);
519: }
521: static PetscErrorCode MatSolve_SeqDense_SetUp(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
522: {
523: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
524: PetscScalar *y;
525: PetscBLASInt m = 0, k = 0;
527: PetscFunctionBegin;
528: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
529: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
530: if (k < m) {
531: PetscCall(VecCopy(xx, mat->qrrhs));
532: PetscCall(VecGetArray(mat->qrrhs, &y));
533: } else {
534: PetscCall(VecCopy(xx, yy));
535: PetscCall(VecGetArray(yy, &y));
536: }
537: *_y = y;
538: *_k = k;
539: *_m = m;
540: PetscFunctionReturn(PETSC_SUCCESS);
541: }
543: static PetscErrorCode MatSolve_SeqDense_TearDown(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
544: {
545: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
546: PetscScalar *y = NULL;
547: PetscBLASInt m, k;
549: PetscFunctionBegin;
550: y = *_y;
551: *_y = NULL;
552: k = *_k;
553: m = *_m;
554: if (k < m) {
555: PetscScalar *yv;
556: PetscCall(VecGetArray(yy, &yv));
557: PetscCall(PetscArraycpy(yv, y, k));
558: PetscCall(VecRestoreArray(yy, &yv));
559: PetscCall(VecRestoreArray(mat->qrrhs, &y));
560: } else {
561: PetscCall(VecRestoreArray(yy, &y));
562: }
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: static PetscErrorCode MatSolve_SeqDense_LU(Mat A, Vec xx, Vec yy)
567: {
568: PetscScalar *y = NULL;
569: PetscBLASInt m = 0, k = 0;
571: PetscFunctionBegin;
572: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
573: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_FALSE));
574: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
578: static PetscErrorCode MatSolveTranspose_SeqDense_LU(Mat A, Vec xx, Vec yy)
579: {
580: PetscScalar *y = NULL;
581: PetscBLASInt m = 0, k = 0;
583: PetscFunctionBegin;
584: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
585: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_TRUE));
586: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
590: static PetscErrorCode MatSolve_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
591: {
592: PetscScalar *y = NULL;
593: PetscBLASInt m = 0, k = 0;
595: PetscFunctionBegin;
596: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
597: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_FALSE));
598: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: static PetscErrorCode MatSolveTranspose_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
603: {
604: PetscScalar *y = NULL;
605: PetscBLASInt m = 0, k = 0;
607: PetscFunctionBegin;
608: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
609: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_TRUE));
610: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
614: static PetscErrorCode MatSolve_SeqDense_QR(Mat A, Vec xx, Vec yy)
615: {
616: PetscScalar *y = NULL;
617: PetscBLASInt m = 0, k = 0;
619: PetscFunctionBegin;
620: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
621: PetscCall(MatSolve_SeqDense_Internal_QR(A, y, PetscMax(m, k), m, 1, k));
622: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
623: PetscFunctionReturn(PETSC_SUCCESS);
624: }
626: static PetscErrorCode MatSolveTranspose_SeqDense_QR(Mat A, Vec xx, Vec yy)
627: {
628: PetscScalar *y = NULL;
629: PetscBLASInt m = 0, k = 0;
631: PetscFunctionBegin;
632: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
633: PetscCall(MatSolveTranspose_SeqDense_Internal_QR(A, y, PetscMax(m, k), m, 1, k));
634: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: static PetscErrorCode MatMatSolve_SeqDense_SetUp(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
639: {
640: const PetscScalar *b;
641: PetscScalar *y;
642: PetscInt n, _ldb, _ldx;
643: PetscBLASInt nrhs = 0, m = 0, k = 0, ldb = 0, ldx = 0, ldy = 0;
645: PetscFunctionBegin;
646: *_ldy = 0;
647: *_m = 0;
648: *_nrhs = 0;
649: *_k = 0;
650: *_y = NULL;
651: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
652: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
653: PetscCall(MatGetSize(B, NULL, &n));
654: PetscCall(PetscBLASIntCast(n, &nrhs));
655: PetscCall(MatDenseGetLDA(B, &_ldb));
656: PetscCall(PetscBLASIntCast(_ldb, &ldb));
657: PetscCall(MatDenseGetLDA(X, &_ldx));
658: PetscCall(PetscBLASIntCast(_ldx, &ldx));
659: if (ldx < m) {
660: PetscCall(MatDenseGetArrayRead(B, &b));
661: PetscCall(PetscMalloc1(nrhs * m, &y));
662: if (ldb == m) {
663: PetscCall(PetscArraycpy(y, b, ldb * nrhs));
664: } else {
665: for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&y[j * m], &b[j * ldb], m));
666: }
667: ldy = m;
668: PetscCall(MatDenseRestoreArrayRead(B, &b));
669: } else {
670: if (ldb == ldx) {
671: PetscCall(MatCopy(B, X, SAME_NONZERO_PATTERN));
672: PetscCall(MatDenseGetArray(X, &y));
673: } else {
674: PetscCall(MatDenseGetArray(X, &y));
675: PetscCall(MatDenseGetArrayRead(B, &b));
676: for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&y[j * ldx], &b[j * ldb], m));
677: PetscCall(MatDenseRestoreArrayRead(B, &b));
678: }
679: ldy = ldx;
680: }
681: *_y = y;
682: *_ldy = ldy;
683: *_k = k;
684: *_m = m;
685: *_nrhs = nrhs;
686: PetscFunctionReturn(PETSC_SUCCESS);
687: }
689: static PetscErrorCode MatMatSolve_SeqDense_TearDown(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
690: {
691: PetscScalar *y;
692: PetscInt _ldx;
693: PetscBLASInt k, ldy, nrhs, ldx = 0;
695: PetscFunctionBegin;
696: y = *_y;
697: *_y = NULL;
698: k = *_k;
699: ldy = *_ldy;
700: nrhs = *_nrhs;
701: PetscCall(MatDenseGetLDA(X, &_ldx));
702: PetscCall(PetscBLASIntCast(_ldx, &ldx));
703: if (ldx != ldy) {
704: PetscScalar *xv;
705: PetscCall(MatDenseGetArray(X, &xv));
706: for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&xv[j * ldx], &y[j * ldy], k));
707: PetscCall(MatDenseRestoreArray(X, &xv));
708: PetscCall(PetscFree(y));
709: } else {
710: PetscCall(MatDenseRestoreArray(X, &y));
711: }
712: PetscFunctionReturn(PETSC_SUCCESS);
713: }
715: static PetscErrorCode MatMatSolve_SeqDense_LU(Mat A, Mat B, Mat X)
716: {
717: PetscScalar *y;
718: PetscBLASInt m, k, ldy, nrhs;
720: PetscFunctionBegin;
721: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
722: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_FALSE));
723: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: static PetscErrorCode MatMatSolveTranspose_SeqDense_LU(Mat A, Mat B, Mat X)
728: {
729: PetscScalar *y;
730: PetscBLASInt m, k, ldy, nrhs;
732: PetscFunctionBegin;
733: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
734: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_TRUE));
735: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
736: PetscFunctionReturn(PETSC_SUCCESS);
737: }
739: static PetscErrorCode MatMatSolve_SeqDense_Cholesky(Mat A, Mat B, Mat X)
740: {
741: PetscScalar *y;
742: PetscBLASInt m, k, ldy, nrhs;
744: PetscFunctionBegin;
745: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
746: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_FALSE));
747: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
748: PetscFunctionReturn(PETSC_SUCCESS);
749: }
751: static PetscErrorCode MatMatSolveTranspose_SeqDense_Cholesky(Mat A, Mat B, Mat X)
752: {
753: PetscScalar *y;
754: PetscBLASInt m, k, ldy, nrhs;
756: PetscFunctionBegin;
757: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
758: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_TRUE));
759: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
760: PetscFunctionReturn(PETSC_SUCCESS);
761: }
763: static PetscErrorCode MatMatSolve_SeqDense_QR(Mat A, Mat B, Mat X)
764: {
765: PetscScalar *y;
766: PetscBLASInt m, k, ldy, nrhs;
768: PetscFunctionBegin;
769: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
770: PetscCall(MatSolve_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k));
771: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
772: PetscFunctionReturn(PETSC_SUCCESS);
773: }
775: static PetscErrorCode MatMatSolveTranspose_SeqDense_QR(Mat A, Mat B, Mat X)
776: {
777: PetscScalar *y;
778: PetscBLASInt m, k, ldy, nrhs;
780: PetscFunctionBegin;
781: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
782: PetscCall(MatSolveTranspose_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k));
783: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
784: PetscFunctionReturn(PETSC_SUCCESS);
785: }
787: /* COMMENT: I have chosen to hide row permutation in the pivots,
788: rather than put it in the Mat->row slot.*/
789: PetscErrorCode MatLUFactor_SeqDense(Mat A, IS row, IS col, PETSC_UNUSED const MatFactorInfo *minfo)
790: {
791: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
792: PetscBLASInt n, m, info;
794: PetscFunctionBegin;
795: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
796: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
797: if (!mat->pivots) PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots));
798: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
799: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
800: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&m, &n, mat->v, &mat->lda, mat->pivots, &info));
801: PetscCall(PetscFPTrapPop());
803: PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to LU factorization %" PetscBLASInt_FMT, info);
804: PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Bad LU factorization %" PetscBLASInt_FMT, info);
806: A->ops->solve = MatSolve_SeqDense_LU;
807: A->ops->matsolve = MatMatSolve_SeqDense_LU;
808: A->ops->solvetranspose = MatSolveTranspose_SeqDense_LU;
809: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_LU;
810: A->factortype = MAT_FACTOR_LU;
812: PetscCall(PetscFree(A->solvertype));
813: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
815: PetscCall(PetscLogFlops((2.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3));
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }
819: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info)
820: {
821: PetscFunctionBegin;
822: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
823: PetscUseTypeMethod(fact, lufactor, NULL, NULL, info);
824: PetscFunctionReturn(PETSC_SUCCESS);
825: }
827: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, IS col, PETSC_UNUSED const MatFactorInfo *info)
828: {
829: PetscFunctionBegin;
830: fact->preallocated = PETSC_TRUE;
831: fact->assembled = PETSC_TRUE;
832: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
833: PetscFunctionReturn(PETSC_SUCCESS);
834: }
836: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
837: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A, IS perm, PETSC_UNUSED const MatFactorInfo *minfo)
838: {
839: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
840: PetscBLASInt info, n;
842: PetscFunctionBegin;
843: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
844: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
845: if (A->spd == PETSC_BOOL3_TRUE) {
846: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
847: PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("L", &n, mat->v, &mat->lda, &info));
848: PetscCall(PetscFPTrapPop());
849: #if defined(PETSC_USE_COMPLEX)
850: } else if (A->hermitian == PETSC_BOOL3_TRUE) {
851: if (!mat->pivots) PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots));
852: if (!mat->fwork) {
853: PetscScalar dummy;
855: mat->lfwork = -1;
856: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
857: PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info));
858: PetscCall(PetscFPTrapPop());
859: PetscCall(PetscBLASIntCast((PetscCount)(PetscRealPart(dummy)), &mat->lfwork));
860: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
861: }
862: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
863: PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
864: PetscCall(PetscFPTrapPop());
865: #endif
866: } else { /* symmetric case */
867: if (!mat->pivots) PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots));
868: if (!mat->fwork) {
869: PetscScalar dummy;
871: mat->lfwork = -1;
872: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
873: PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info));
874: PetscCall(PetscFPTrapPop());
875: PetscCall(PetscBLASIntCast((PetscCount)(PetscRealPart(dummy)), &mat->lfwork));
876: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
877: }
878: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
879: PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
880: PetscCall(PetscFPTrapPop());
881: }
882: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad factorization: zero pivot in row %" PetscBLASInt_FMT, info - 1);
884: A->ops->solve = MatSolve_SeqDense_Cholesky;
885: A->ops->matsolve = MatMatSolve_SeqDense_Cholesky;
886: A->ops->solvetranspose = MatSolveTranspose_SeqDense_Cholesky;
887: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky;
888: A->factortype = MAT_FACTOR_CHOLESKY;
890: PetscCall(PetscFree(A->solvertype));
891: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
893: PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
894: PetscFunctionReturn(PETSC_SUCCESS);
895: }
897: static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info)
898: {
899: PetscFunctionBegin;
900: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
901: PetscUseTypeMethod(fact, choleskyfactor, NULL, info);
902: PetscFunctionReturn(PETSC_SUCCESS);
903: }
905: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info)
906: {
907: PetscFunctionBegin;
908: fact->assembled = PETSC_TRUE;
909: fact->preallocated = PETSC_TRUE;
910: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
911: PetscFunctionReturn(PETSC_SUCCESS);
912: }
914: PetscErrorCode MatQRFactor_SeqDense(Mat A, IS col, PETSC_UNUSED const MatFactorInfo *minfo)
915: {
916: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
917: PetscBLASInt n, m, info, min, max;
919: PetscFunctionBegin;
920: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
921: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
922: max = PetscMax(m, n);
923: min = PetscMin(m, n);
924: if (!mat->tau) PetscCall(PetscMalloc1(min, &mat->tau));
925: if (!mat->pivots) PetscCall(PetscMalloc1(n, &mat->pivots));
926: if (!mat->qrrhs) PetscCall(MatCreateVecs(A, NULL, &mat->qrrhs));
927: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
928: if (!mat->fwork) {
929: PetscScalar dummy;
931: mat->lfwork = -1;
932: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
933: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, &dummy, &mat->lfwork, &info));
934: PetscCall(PetscFPTrapPop());
935: PetscCall(PetscBLASIntCast((PetscCount)(PetscRealPart(dummy)), &mat->lfwork));
936: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
937: }
938: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
939: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, mat->fwork, &mat->lfwork, &info));
940: PetscCall(PetscFPTrapPop());
941: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to QR factorization %" PetscBLASInt_FMT, info);
942: // 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
943: mat->rank = min;
945: A->ops->solve = MatSolve_SeqDense_QR;
946: A->ops->matsolve = MatMatSolve_SeqDense_QR;
947: A->factortype = MAT_FACTOR_QR;
948: if (m == n) {
949: A->ops->solvetranspose = MatSolveTranspose_SeqDense_QR;
950: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
951: }
953: PetscCall(PetscFree(A->solvertype));
954: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
956: PetscCall(PetscLogFlops(2.0 * min * min * (max - min / 3.0)));
957: PetscFunctionReturn(PETSC_SUCCESS);
958: }
960: static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info)
961: {
962: PetscFunctionBegin;
963: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
964: PetscUseMethod(fact, "MatQRFactor_C", (Mat, IS, const MatFactorInfo *), (fact, NULL, info));
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }
968: PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info)
969: {
970: PetscFunctionBegin;
971: fact->assembled = PETSC_TRUE;
972: fact->preallocated = PETSC_TRUE;
973: PetscCall(PetscObjectComposeFunction((PetscObject)fact, "MatQRFactorNumeric_C", MatQRFactorNumeric_SeqDense));
974: PetscFunctionReturn(PETSC_SUCCESS);
975: }
977: /* uses LAPACK */
978: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A, MatFactorType ftype, Mat *fact)
979: {
980: PetscFunctionBegin;
981: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), fact));
982: PetscCall(MatSetSizes(*fact, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
983: PetscCall(MatSetType(*fact, MATDENSE));
984: (*fact)->trivialsymbolic = PETSC_TRUE;
985: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
986: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
987: (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
988: } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
989: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
990: } else if (ftype == MAT_FACTOR_QR) {
991: PetscCall(PetscObjectComposeFunction((PetscObject)*fact, "MatQRFactorSymbolic_C", MatQRFactorSymbolic_SeqDense));
992: }
993: (*fact)->factortype = ftype;
995: PetscCall(PetscFree((*fact)->solvertype));
996: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*fact)->solvertype));
997: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_LU]));
998: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ILU]));
999: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_CHOLESKY]));
1000: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ICC]));
1001: PetscFunctionReturn(PETSC_SUCCESS);
1002: }
1004: static PetscErrorCode MatSOR_SeqDense(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal shift, PetscInt its, PetscInt lits, Vec xx)
1005: {
1006: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1007: PetscScalar *x, *v = mat->v, zero = 0.0, xt;
1008: const PetscScalar *b;
1009: PetscInt m = A->rmap->n, i;
1010: PetscBLASInt o = 1, bm = 0;
1012: PetscFunctionBegin;
1013: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1014: PetscCheck(A->offloadmask != PETSC_OFFLOAD_GPU, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented");
1015: #endif
1016: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
1017: PetscCall(PetscBLASIntCast(m, &bm));
1018: if (flag & SOR_ZERO_INITIAL_GUESS) {
1019: /* this is a hack fix, should have another version without the second BLASdotu */
1020: PetscCall(VecSet(xx, zero));
1021: }
1022: PetscCall(VecGetArray(xx, &x));
1023: PetscCall(VecGetArrayRead(bb, &b));
1024: its = its * lits;
1025: 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);
1026: while (its--) {
1027: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1028: for (i = 0; i < m; i++) {
1029: PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o));
1030: x[i] = (1. - omega) * x[i] + (xt + v[i + i * m] * x[i]) * omega / (v[i + i * m] + shift);
1031: }
1032: }
1033: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1034: for (i = m - 1; i >= 0; i--) {
1035: PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o));
1036: x[i] = (1. - omega) * x[i] + (xt + v[i + i * m] * x[i]) * omega / (v[i + i * m] + shift);
1037: }
1038: }
1039: }
1040: PetscCall(VecRestoreArrayRead(bb, &b));
1041: PetscCall(VecRestoreArray(xx, &x));
1042: PetscFunctionReturn(PETSC_SUCCESS);
1043: }
1045: PETSC_INTERN PetscErrorCode MatMultColumnRangeKernel_SeqDense(Mat A, Vec xx, Vec yy, PetscInt c_start, PetscInt c_end, PetscBool trans, PetscBool herm)
1046: {
1047: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1048: PetscScalar *y, _DOne = 1.0, _DZero = 0.0;
1049: PetscBLASInt m, n, _One = 1;
1050: const PetscScalar *v = mat->v, *x;
1052: PetscFunctionBegin;
1053: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1054: PetscCall(PetscBLASIntCast(c_end - c_start, &n));
1055: PetscCall(VecGetArrayRead(xx, &x));
1056: PetscCall(VecGetArrayWrite(yy, &y));
1057: if (!m || !n) {
1058: PetscBLASInt i;
1059: if (trans)
1060: for (i = 0; i < n; i++) y[i] = 0.0;
1061: else
1062: for (i = 0; i < m; i++) y[i] = 0.0;
1063: } else {
1064: if (trans) {
1065: if (herm) PetscCallBLAS("BLASgemv", BLASgemv_("C", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DZero, y + c_start, &_One));
1066: else PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DZero, y + c_start, &_One));
1067: } else {
1068: PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x + c_start, &_One, &_DZero, y, &_One));
1069: }
1070: PetscCall(PetscLogFlops(2.0 * m * n - n));
1071: }
1072: PetscCall(VecRestoreArrayRead(xx, &x));
1073: PetscCall(VecRestoreArrayWrite(yy, &y));
1074: PetscFunctionReturn(PETSC_SUCCESS);
1075: }
1077: PetscErrorCode MatMultHermitianTransposeColumnRange_SeqDense(Mat A, Vec xx, Vec yy, PetscInt c_start, PetscInt c_end)
1078: {
1079: PetscFunctionBegin;
1080: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, c_start, c_end, PETSC_TRUE, PETSC_TRUE));
1081: PetscFunctionReturn(PETSC_SUCCESS);
1082: }
1084: PetscErrorCode MatMult_SeqDense(Mat A, Vec xx, Vec yy)
1085: {
1086: PetscFunctionBegin;
1087: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_FALSE, PETSC_FALSE));
1088: PetscFunctionReturn(PETSC_SUCCESS);
1089: }
1091: PetscErrorCode MatMultTranspose_SeqDense(Mat A, Vec xx, Vec yy)
1092: {
1093: PetscFunctionBegin;
1094: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_FALSE));
1095: PetscFunctionReturn(PETSC_SUCCESS);
1096: }
1098: PetscErrorCode MatMultHermitianTranspose_SeqDense(Mat A, Vec xx, Vec yy)
1099: {
1100: PetscFunctionBegin;
1101: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_TRUE));
1102: PetscFunctionReturn(PETSC_SUCCESS);
1103: }
1105: PETSC_INTERN PetscErrorCode MatMultAddColumnRangeKernel_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end, PetscBool trans, PetscBool herm)
1106: {
1107: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1108: const PetscScalar *v = mat->v, *x;
1109: PetscScalar *y, _DOne = 1.0;
1110: PetscBLASInt m, n, _One = 1;
1112: PetscFunctionBegin;
1113: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1114: PetscCall(PetscBLASIntCast(c_end - c_start, &n));
1115: PetscCall(VecCopy(zz, yy));
1116: if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS);
1117: PetscCall(VecGetArray(yy, &y));
1118: PetscCall(VecGetArrayRead(xx, &x));
1119: if (trans) {
1120: if (herm) PetscCallBLAS("BLASgemv", BLASgemv_("C", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DOne, y + c_start, &_One));
1121: else PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DOne, y + c_start, &_One));
1122: } else {
1123: PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x + c_start, &_One, &_DOne, y, &_One));
1124: }
1125: PetscCall(VecRestoreArrayRead(xx, &x));
1126: PetscCall(VecRestoreArray(yy, &y));
1127: PetscCall(PetscLogFlops(2.0 * m * n));
1128: PetscFunctionReturn(PETSC_SUCCESS);
1129: }
1131: PetscErrorCode MatMultColumnRange_SeqDense(Mat A, Vec xx, Vec yy, PetscInt c_start, PetscInt c_end)
1132: {
1133: PetscFunctionBegin;
1134: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, c_start, c_end, PETSC_FALSE, PETSC_FALSE));
1135: PetscFunctionReturn(PETSC_SUCCESS);
1136: }
1138: PetscErrorCode MatMultAddColumnRange_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end)
1139: {
1140: PetscFunctionBegin;
1141: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, c_start, c_end, PETSC_FALSE, PETSC_FALSE));
1142: PetscFunctionReturn(PETSC_SUCCESS);
1143: }
1145: PetscErrorCode MatMultHermitianTransposeAddColumnRange_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end)
1146: {
1147: PetscFunctionBegin;
1148: PetscMPIInt rank;
1149: PetscCallMPI(MPI_Comm_rank(MPI_COMM_WORLD, &rank));
1150: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, c_start, c_end, PETSC_TRUE, PETSC_TRUE));
1151: PetscFunctionReturn(PETSC_SUCCESS);
1152: }
1154: PetscErrorCode MatMultAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1155: {
1156: PetscFunctionBegin;
1157: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_FALSE, PETSC_FALSE));
1158: PetscFunctionReturn(PETSC_SUCCESS);
1159: }
1161: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1162: {
1163: PetscFunctionBegin;
1164: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_FALSE));
1165: PetscFunctionReturn(PETSC_SUCCESS);
1166: }
1168: PetscErrorCode MatMultHermitianTransposeAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1169: {
1170: PetscFunctionBegin;
1171: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_TRUE));
1172: PetscFunctionReturn(PETSC_SUCCESS);
1173: }
1175: static PetscErrorCode MatGetRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1176: {
1177: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1178: PetscInt i;
1180: PetscFunctionBegin;
1181: if (ncols) *ncols = A->cmap->n;
1182: if (cols) {
1183: PetscCall(PetscMalloc1(A->cmap->n, cols));
1184: for (i = 0; i < A->cmap->n; i++) (*cols)[i] = i;
1185: }
1186: if (vals) {
1187: const PetscScalar *v;
1189: PetscCall(MatDenseGetArrayRead(A, &v));
1190: PetscCall(PetscMalloc1(A->cmap->n, vals));
1191: v += row;
1192: for (i = 0; i < A->cmap->n; i++) {
1193: (*vals)[i] = *v;
1194: v += mat->lda;
1195: }
1196: PetscCall(MatDenseRestoreArrayRead(A, &v));
1197: }
1198: PetscFunctionReturn(PETSC_SUCCESS);
1199: }
1201: static PetscErrorCode MatRestoreRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1202: {
1203: PetscFunctionBegin;
1204: if (cols) PetscCall(PetscFree(*cols));
1205: if (vals) PetscCall(PetscFree(*vals));
1206: PetscFunctionReturn(PETSC_SUCCESS);
1207: }
1209: static PetscErrorCode MatSetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], const PetscScalar v[], InsertMode addv)
1210: {
1211: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1212: PetscScalar *av;
1213: PetscInt i, j, idx = 0;
1214: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1215: PetscOffloadMask oldf;
1216: #endif
1218: PetscFunctionBegin;
1219: PetscCall(MatDenseGetArray(A, &av));
1220: if (!mat->roworiented) {
1221: if (addv == INSERT_VALUES) {
1222: for (j = 0; j < n; j++) {
1223: if (indexn[j] < 0) {
1224: idx += m;
1225: continue;
1226: }
1227: 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);
1228: for (i = 0; i < m; i++) {
1229: if (indexm[i] < 0) {
1230: idx++;
1231: continue;
1232: }
1233: 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);
1234: av[indexn[j] * mat->lda + indexm[i]] = v ? v[idx++] : (idx++, 0.0);
1235: }
1236: }
1237: } else {
1238: for (j = 0; j < n; j++) {
1239: if (indexn[j] < 0) {
1240: idx += m;
1241: continue;
1242: }
1243: PetscCheck(indexn[j] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, indexn[j], A->cmap->n - 1);
1244: for (i = 0; i < m; i++) {
1245: if (indexm[i] < 0) {
1246: idx++;
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: av[indexn[j] * mat->lda + indexm[i]] += v ? v[idx++] : (idx++, 0.0);
1251: }
1252: }
1253: }
1254: } else {
1255: if (addv == INSERT_VALUES) {
1256: for (i = 0; i < m; i++) {
1257: if (indexm[i] < 0) {
1258: idx += n;
1259: continue;
1260: }
1261: 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);
1262: for (j = 0; j < n; j++) {
1263: if (indexn[j] < 0) {
1264: idx++;
1265: continue;
1266: }
1267: 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);
1268: av[indexn[j] * mat->lda + indexm[i]] = v ? v[idx++] : (idx++, 0.0);
1269: }
1270: }
1271: } else {
1272: for (i = 0; i < m; i++) {
1273: if (indexm[i] < 0) {
1274: idx += n;
1275: continue;
1276: }
1277: 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);
1278: for (j = 0; j < n; j++) {
1279: if (indexn[j] < 0) {
1280: idx++;
1281: continue;
1282: }
1283: 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);
1284: av[indexn[j] * mat->lda + indexm[i]] += v ? v[idx++] : (idx++, 0.0);
1285: }
1286: }
1287: }
1288: }
1289: /* hack to prevent unneeded copy to the GPU while returning the array */
1290: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1291: oldf = A->offloadmask;
1292: A->offloadmask = PETSC_OFFLOAD_GPU;
1293: #endif
1294: PetscCall(MatDenseRestoreArray(A, &av));
1295: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1296: A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1297: #endif
1298: PetscFunctionReturn(PETSC_SUCCESS);
1299: }
1301: static PetscErrorCode MatGetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], PetscScalar v[])
1302: {
1303: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1304: const PetscScalar *vv;
1305: PetscInt i, j;
1307: PetscFunctionBegin;
1308: PetscCall(MatDenseGetArrayRead(A, &vv));
1309: /* row-oriented output */
1310: for (i = 0; i < m; i++) {
1311: if (indexm[i] < 0) {
1312: v += n;
1313: continue;
1314: }
1315: 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);
1316: for (j = 0; j < n; j++) {
1317: if (indexn[j] < 0) {
1318: v++;
1319: continue;
1320: }
1321: 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);
1322: *v++ = vv[indexn[j] * mat->lda + indexm[i]];
1323: }
1324: }
1325: PetscCall(MatDenseRestoreArrayRead(A, &vv));
1326: PetscFunctionReturn(PETSC_SUCCESS);
1327: }
1329: PetscErrorCode MatView_Dense_Binary(Mat mat, PetscViewer viewer)
1330: {
1331: PetscBool skipHeader;
1332: PetscViewerFormat format;
1333: PetscInt header[4], M, N, m, lda, i, j;
1334: PetscCount k;
1335: const PetscScalar *v;
1336: PetscScalar *vwork;
1338: PetscFunctionBegin;
1339: PetscCall(PetscViewerSetUp(viewer));
1340: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
1341: PetscCall(PetscViewerGetFormat(viewer, &format));
1342: if (skipHeader) format = PETSC_VIEWER_NATIVE;
1344: PetscCall(MatGetSize(mat, &M, &N));
1346: /* write matrix header */
1347: header[0] = MAT_FILE_CLASSID;
1348: header[1] = M;
1349: header[2] = N;
1350: header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M * N;
1351: if (!skipHeader) PetscCall(PetscViewerBinaryWrite(viewer, header, 4, PETSC_INT));
1353: PetscCall(MatGetLocalSize(mat, &m, NULL));
1354: if (format != PETSC_VIEWER_NATIVE) {
1355: PetscInt nnz = m * N, *iwork;
1356: /* store row lengths for each row */
1357: PetscCall(PetscMalloc1(nnz, &iwork));
1358: for (i = 0; i < m; i++) iwork[i] = N;
1359: PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1360: /* store column indices (zero start index) */
1361: for (k = 0, i = 0; i < m; i++)
1362: for (j = 0; j < N; j++, k++) iwork[k] = j;
1363: PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1364: PetscCall(PetscFree(iwork));
1365: }
1366: /* store matrix values as a dense matrix in row major order */
1367: PetscCall(PetscMalloc1(m * N, &vwork));
1368: PetscCall(MatDenseGetArrayRead(mat, &v));
1369: PetscCall(MatDenseGetLDA(mat, &lda));
1370: for (k = 0, i = 0; i < m; i++)
1371: for (j = 0; j < N; j++, k++) vwork[k] = v[i + (size_t)lda * j];
1372: PetscCall(MatDenseRestoreArrayRead(mat, &v));
1373: PetscCall(PetscViewerBinaryWriteAll(viewer, vwork, m * N, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1374: PetscCall(PetscFree(vwork));
1375: PetscFunctionReturn(PETSC_SUCCESS);
1376: }
1378: PetscErrorCode MatLoad_Dense_Binary(Mat mat, PetscViewer viewer)
1379: {
1380: PetscBool skipHeader;
1381: PetscInt header[4], M, N, m, nz, lda, i, j, k;
1382: PetscInt rows, cols;
1383: PetscScalar *v, *vwork;
1385: PetscFunctionBegin;
1386: PetscCall(PetscViewerSetUp(viewer));
1387: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
1389: if (!skipHeader) {
1390: PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT));
1391: PetscCheck(header[0] == MAT_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file");
1392: M = header[1];
1393: N = header[2];
1394: PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M);
1395: PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N);
1396: nz = header[3];
1397: PetscCheck(nz == MATRIX_BINARY_FORMAT_DENSE || nz >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Unknown matrix format %" PetscInt_FMT " in file", nz);
1398: } else {
1399: PetscCall(MatGetSize(mat, &M, &N));
1400: 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");
1401: nz = MATRIX_BINARY_FORMAT_DENSE;
1402: }
1404: /* setup global sizes if not set */
1405: if (mat->rmap->N < 0) mat->rmap->N = M;
1406: if (mat->cmap->N < 0) mat->cmap->N = N;
1407: PetscCall(MatSetUp(mat));
1408: /* check if global sizes are correct */
1409: PetscCall(MatGetSize(mat, &rows, &cols));
1410: 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);
1412: PetscCall(MatGetSize(mat, NULL, &N));
1413: PetscCall(MatGetLocalSize(mat, &m, NULL));
1414: PetscCall(MatDenseGetArray(mat, &v));
1415: PetscCall(MatDenseGetLDA(mat, &lda));
1416: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */
1417: PetscCount nnz = (size_t)m * N;
1418: /* read in matrix values */
1419: PetscCall(PetscMalloc1(nnz, &vwork));
1420: PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1421: /* store values in column major order */
1422: for (j = 0; j < N; j++)
1423: for (i = 0; i < m; i++) v[i + (size_t)lda * j] = vwork[(size_t)i * N + j];
1424: PetscCall(PetscFree(vwork));
1425: } else { /* matrix in file is sparse format */
1426: PetscInt nnz = 0, *rlens, *icols;
1427: /* read in row lengths */
1428: PetscCall(PetscMalloc1(m, &rlens));
1429: PetscCall(PetscViewerBinaryReadAll(viewer, rlens, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1430: for (i = 0; i < m; i++) nnz += rlens[i];
1431: /* read in column indices and values */
1432: PetscCall(PetscMalloc2(nnz, &icols, nnz, &vwork));
1433: PetscCall(PetscViewerBinaryReadAll(viewer, icols, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1434: PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1435: /* store values in column major order */
1436: for (k = 0, i = 0; i < m; i++)
1437: for (j = 0; j < rlens[i]; j++, k++) v[i + lda * icols[k]] = vwork[k];
1438: PetscCall(PetscFree(rlens));
1439: PetscCall(PetscFree2(icols, vwork));
1440: }
1441: PetscCall(MatDenseRestoreArray(mat, &v));
1442: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
1443: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
1444: PetscFunctionReturn(PETSC_SUCCESS);
1445: }
1447: static PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1448: {
1449: PetscBool isbinary, ishdf5;
1451: PetscFunctionBegin;
1454: /* force binary viewer to load .info file if it has not yet done so */
1455: PetscCall(PetscViewerSetUp(viewer));
1456: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1457: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1458: if (isbinary) {
1459: PetscCall(MatLoad_Dense_Binary(newMat, viewer));
1460: } else if (ishdf5) {
1461: #if defined(PETSC_HAVE_HDF5)
1462: PetscCall(MatLoad_Dense_HDF5(newMat, viewer));
1463: #else
1464: SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1465: #endif
1466: } else {
1467: 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);
1468: }
1469: PetscFunctionReturn(PETSC_SUCCESS);
1470: }
1472: static PetscErrorCode MatView_SeqDense_ASCII(Mat A, PetscViewer viewer)
1473: {
1474: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1475: PetscInt i, j;
1476: const char *name;
1477: PetscScalar *v, *av;
1478: PetscViewerFormat format;
1479: #if defined(PETSC_USE_COMPLEX)
1480: PetscBool allreal = PETSC_TRUE;
1481: #endif
1483: PetscFunctionBegin;
1484: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&av));
1485: PetscCall(PetscViewerGetFormat(viewer, &format));
1486: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1487: PetscFunctionReturn(PETSC_SUCCESS); /* do nothing for now */
1488: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1489: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1490: for (i = 0; i < A->rmap->n; i++) {
1491: v = av + i;
1492: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i));
1493: for (j = 0; j < A->cmap->n; j++) {
1494: #if defined(PETSC_USE_COMPLEX)
1495: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1496: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", j, (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1497: } else if (PetscRealPart(*v)) {
1498: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)PetscRealPart(*v)));
1499: }
1500: #else
1501: if (*v) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)*v));
1502: #endif
1503: v += a->lda;
1504: }
1505: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1506: }
1507: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1508: } else {
1509: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1510: #if defined(PETSC_USE_COMPLEX)
1511: /* determine if matrix has all real values */
1512: for (j = 0; j < A->cmap->n; j++) {
1513: v = av + j * a->lda;
1514: for (i = 0; i < A->rmap->n; i++) {
1515: if (PetscImaginaryPart(v[i])) {
1516: allreal = PETSC_FALSE;
1517: break;
1518: }
1519: }
1520: }
1521: #endif
1522: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1523: PetscCall(PetscObjectGetName((PetscObject)A, &name));
1524: PetscCall(PetscViewerASCIIPrintf(viewer, "%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n", A->rmap->n, A->cmap->n));
1525: PetscCall(PetscViewerASCIIPrintf(viewer, "%s = zeros(%" PetscInt_FMT ",%" PetscInt_FMT ");\n", name, A->rmap->n, A->cmap->n));
1526: PetscCall(PetscViewerASCIIPrintf(viewer, "%s = [\n", name));
1527: }
1529: for (i = 0; i < A->rmap->n; i++) {
1530: v = av + i;
1531: for (j = 0; j < A->cmap->n; j++) {
1532: #if defined(PETSC_USE_COMPLEX)
1533: if (allreal) {
1534: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)PetscRealPart(*v)));
1535: } else {
1536: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei ", (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1537: }
1538: #else
1539: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)*v));
1540: #endif
1541: v += a->lda;
1542: }
1543: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1544: }
1545: if (format == PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer, "];\n"));
1546: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1547: }
1548: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&av));
1549: PetscCall(PetscViewerFlush(viewer));
1550: PetscFunctionReturn(PETSC_SUCCESS);
1551: }
1553: #include <petscdraw.h>
1554: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw, void *Aa)
1555: {
1556: Mat A = (Mat)Aa;
1557: PetscInt m = A->rmap->n, n = A->cmap->n, i, j;
1558: int color = PETSC_DRAW_WHITE;
1559: const PetscScalar *v;
1560: PetscViewer viewer;
1561: PetscReal xl, yl, xr, yr, x_l, x_r, y_l, y_r;
1562: PetscViewerFormat format;
1564: PetscFunctionBegin;
1565: PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer));
1566: PetscCall(PetscViewerGetFormat(viewer, &format));
1567: PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));
1569: /* Loop over matrix elements drawing boxes */
1570: PetscCall(MatDenseGetArrayRead(A, &v));
1571: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1572: PetscDrawCollectiveBegin(draw);
1573: /* Blue for negative and Red for positive */
1574: for (j = 0; j < n; j++) {
1575: x_l = j;
1576: x_r = x_l + 1.0;
1577: for (i = 0; i < m; i++) {
1578: y_l = m - i - 1.0;
1579: y_r = y_l + 1.0;
1580: if (PetscRealPart(v[j * m + i]) > 0.) color = PETSC_DRAW_RED;
1581: else if (PetscRealPart(v[j * m + i]) < 0.) color = PETSC_DRAW_BLUE;
1582: else continue;
1583: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1584: }
1585: }
1586: PetscDrawCollectiveEnd(draw);
1587: } else {
1588: /* use contour shading to indicate magnitude of values */
1589: /* first determine max of all nonzero values */
1590: PetscReal minv = 0.0, maxv = 0.0;
1591: PetscDraw popup;
1593: for (i = 0; i < m * n; i++) {
1594: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1595: }
1596: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1597: PetscCall(PetscDrawGetPopup(draw, &popup));
1598: PetscCall(PetscDrawScalePopup(popup, minv, maxv));
1600: PetscDrawCollectiveBegin(draw);
1601: for (j = 0; j < n; j++) {
1602: x_l = j;
1603: x_r = x_l + 1.0;
1604: for (i = 0; i < m; i++) {
1605: y_l = m - i - 1.0;
1606: y_r = y_l + 1.0;
1607: color = PetscDrawRealToColor(PetscAbsScalar(v[j * m + i]), minv, maxv);
1608: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1609: }
1610: }
1611: PetscDrawCollectiveEnd(draw);
1612: }
1613: PetscCall(MatDenseRestoreArrayRead(A, &v));
1614: PetscFunctionReturn(PETSC_SUCCESS);
1615: }
1617: static PetscErrorCode MatView_SeqDense_Draw(Mat A, PetscViewer viewer)
1618: {
1619: PetscDraw draw;
1620: PetscBool isnull;
1621: PetscReal xr, yr, xl, yl, h, w;
1623: PetscFunctionBegin;
1624: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1625: PetscCall(PetscDrawIsNull(draw, &isnull));
1626: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
1628: xr = A->cmap->n;
1629: yr = A->rmap->n;
1630: h = yr / 10.0;
1631: w = xr / 10.0;
1632: xr += w;
1633: yr += h;
1634: xl = -w;
1635: yl = -h;
1636: PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr));
1637: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer));
1638: PetscCall(PetscDrawZoom(draw, MatView_SeqDense_Draw_Zoom, A));
1639: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL));
1640: PetscCall(PetscDrawSave(draw));
1641: PetscFunctionReturn(PETSC_SUCCESS);
1642: }
1644: PetscErrorCode MatView_SeqDense(Mat A, PetscViewer viewer)
1645: {
1646: PetscBool isascii, isbinary, isdraw;
1648: PetscFunctionBegin;
1649: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1650: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1651: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1652: if (isascii) PetscCall(MatView_SeqDense_ASCII(A, viewer));
1653: else if (isbinary) PetscCall(MatView_Dense_Binary(A, viewer));
1654: else if (isdraw) PetscCall(MatView_SeqDense_Draw(A, viewer));
1655: PetscFunctionReturn(PETSC_SUCCESS);
1656: }
1658: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A, const PetscScalar *array)
1659: {
1660: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1662: PetscFunctionBegin;
1663: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1664: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1665: PetscCheck(!a->unplacedarray, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseResetArray() first");
1666: a->unplacedarray = a->v;
1667: a->unplaced_user_alloc = a->user_alloc;
1668: a->v = (PetscScalar *)array;
1669: a->user_alloc = PETSC_TRUE;
1670: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1671: A->offloadmask = PETSC_OFFLOAD_CPU;
1672: #endif
1673: PetscFunctionReturn(PETSC_SUCCESS);
1674: }
1676: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1677: {
1678: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1680: PetscFunctionBegin;
1681: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1682: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1683: a->v = a->unplacedarray;
1684: a->user_alloc = a->unplaced_user_alloc;
1685: a->unplacedarray = NULL;
1686: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1687: A->offloadmask = PETSC_OFFLOAD_CPU;
1688: #endif
1689: PetscFunctionReturn(PETSC_SUCCESS);
1690: }
1692: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A, const PetscScalar *array)
1693: {
1694: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1696: PetscFunctionBegin;
1697: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1698: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1699: if (!a->user_alloc) PetscCall(PetscFree(a->v));
1700: a->v = (PetscScalar *)array;
1701: a->user_alloc = PETSC_FALSE;
1702: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1703: A->offloadmask = PETSC_OFFLOAD_CPU;
1704: #endif
1705: PetscFunctionReturn(PETSC_SUCCESS);
1706: }
1708: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1709: {
1710: Mat_SeqDense *l = (Mat_SeqDense *)mat->data;
1712: PetscFunctionBegin;
1713: PetscCall(PetscLogObjectState((PetscObject)mat, "Rows %" PetscInt_FMT " Cols %" PetscInt_FMT, mat->rmap->n, mat->cmap->n));
1714: PetscCall(VecDestroy(&l->qrrhs));
1715: PetscCall(PetscFree(l->tau));
1716: PetscCall(PetscFree(l->pivots));
1717: PetscCall(PetscFree(l->fwork));
1718: if (!l->user_alloc) PetscCall(PetscFree(l->v));
1719: if (!l->unplaced_user_alloc) PetscCall(PetscFree(l->unplacedarray));
1720: PetscCheck(!l->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1721: PetscCheck(!l->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1722: PetscCall(VecDestroy(&l->cvec));
1723: PetscCall(MatDestroy(&l->cmat));
1724: PetscCall(PetscFree(mat->data));
1726: PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL));
1727: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactor_C", NULL));
1728: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorSymbolic_C", NULL));
1729: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorNumeric_C", NULL));
1730: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetLDA_C", NULL));
1731: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseSetLDA_C", NULL));
1732: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArray_C", NULL));
1733: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArray_C", NULL));
1734: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDensePlaceArray_C", NULL));
1735: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseResetArray_C", NULL));
1736: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseReplaceArray_C", NULL));
1737: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayRead_C", NULL));
1738: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayRead_C", NULL));
1739: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayWrite_C", NULL));
1740: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayWrite_C", NULL));
1741: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqaij_C", NULL));
1742: #if defined(PETSC_HAVE_ELEMENTAL)
1743: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_elemental_C", NULL));
1744: #endif
1745: #if defined(PETSC_HAVE_SCALAPACK)
1746: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_scalapack_C", NULL));
1747: #endif
1748: #if defined(PETSC_HAVE_CUDA)
1749: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensecuda_C", NULL));
1750: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", NULL));
1751: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdense_C", NULL));
1752: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensecuda_C", NULL));
1753: #endif
1754: #if defined(PETSC_HAVE_HIP)
1755: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensehip_C", NULL));
1756: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", NULL));
1757: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdense_C", NULL));
1758: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensehip_C", NULL));
1759: #endif
1760: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatSeqDenseSetPreallocation_C", NULL));
1761: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqaij_seqdense_C", NULL));
1762: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdense_C", NULL));
1763: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqbaij_seqdense_C", NULL));
1764: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqsbaij_seqdense_C", NULL));
1766: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumn_C", NULL));
1767: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumn_C", NULL));
1768: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVec_C", NULL));
1769: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVec_C", NULL));
1770: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecRead_C", NULL));
1771: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecRead_C", NULL));
1772: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecWrite_C", NULL));
1773: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecWrite_C", NULL));
1774: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetSubMatrix_C", NULL));
1775: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreSubMatrix_C", NULL));
1776: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultColumnRange_C", NULL));
1777: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultAddColumnRange_C", NULL));
1778: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeColumnRange_C", NULL));
1779: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeAddColumnRange_C", NULL));
1780: PetscFunctionReturn(PETSC_SUCCESS);
1781: }
1783: static PetscErrorCode MatTranspose_SeqDense(Mat A, MatReuse reuse, Mat *matout)
1784: {
1785: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1786: PetscInt k, j, m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1787: PetscScalar *v, tmp;
1789: PetscFunctionBegin;
1790: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *matout));
1791: if (reuse == MAT_INPLACE_MATRIX) {
1792: if (m == n) { /* in place transpose */
1793: PetscCall(MatDenseGetArray(A, &v));
1794: for (j = 0; j < m; j++) {
1795: for (k = 0; k < j; k++) {
1796: tmp = v[j + k * M];
1797: v[j + k * M] = v[k + j * M];
1798: v[k + j * M] = tmp;
1799: }
1800: }
1801: PetscCall(MatDenseRestoreArray(A, &v));
1802: } else { /* reuse memory, temporary allocates new memory */
1803: PetscScalar *v2;
1804: PetscLayout tmplayout;
1806: PetscCall(PetscMalloc1((size_t)m * n, &v2));
1807: PetscCall(MatDenseGetArray(A, &v));
1808: for (j = 0; j < n; j++) {
1809: for (k = 0; k < m; k++) v2[j + (size_t)k * n] = v[k + (size_t)j * M];
1810: }
1811: PetscCall(PetscArraycpy(v, v2, (size_t)m * n));
1812: PetscCall(PetscFree(v2));
1813: PetscCall(MatDenseRestoreArray(A, &v));
1814: /* cleanup size dependent quantities */
1815: PetscCall(VecDestroy(&mat->cvec));
1816: PetscCall(MatDestroy(&mat->cmat));
1817: PetscCall(PetscFree(mat->pivots));
1818: PetscCall(PetscFree(mat->fwork));
1819: /* swap row/col layouts */
1820: PetscCall(PetscBLASIntCast(n, &mat->lda));
1821: tmplayout = A->rmap;
1822: A->rmap = A->cmap;
1823: A->cmap = tmplayout;
1824: }
1825: } else { /* out-of-place transpose */
1826: Mat tmat;
1827: Mat_SeqDense *tmatd;
1828: PetscScalar *v2;
1829: PetscInt M2;
1831: if (reuse == MAT_INITIAL_MATRIX) {
1832: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &tmat));
1833: PetscCall(MatSetSizes(tmat, A->cmap->n, A->rmap->n, A->cmap->n, A->rmap->n));
1834: PetscCall(MatSetType(tmat, ((PetscObject)A)->type_name));
1835: PetscCall(MatSeqDenseSetPreallocation(tmat, NULL));
1836: } else tmat = *matout;
1838: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&v));
1839: PetscCall(MatDenseGetArray(tmat, &v2));
1840: tmatd = (Mat_SeqDense *)tmat->data;
1841: M2 = tmatd->lda;
1842: for (j = 0; j < n; j++) {
1843: for (k = 0; k < m; k++) v2[j + k * M2] = v[k + j * M];
1844: }
1845: PetscCall(MatDenseRestoreArray(tmat, &v2));
1846: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&v));
1847: PetscCall(MatAssemblyBegin(tmat, MAT_FINAL_ASSEMBLY));
1848: PetscCall(MatAssemblyEnd(tmat, MAT_FINAL_ASSEMBLY));
1849: *matout = tmat;
1850: }
1851: PetscFunctionReturn(PETSC_SUCCESS);
1852: }
1854: static PetscErrorCode MatEqual_SeqDense(Mat A1, Mat A2, PetscBool *flg)
1855: {
1856: Mat_SeqDense *mat1 = (Mat_SeqDense *)A1->data;
1857: Mat_SeqDense *mat2 = (Mat_SeqDense *)A2->data;
1858: PetscInt i;
1859: const PetscScalar *v1, *v2;
1861: PetscFunctionBegin;
1862: if (A1->rmap->n != A2->rmap->n) {
1863: *flg = PETSC_FALSE;
1864: PetscFunctionReturn(PETSC_SUCCESS);
1865: }
1866: if (A1->cmap->n != A2->cmap->n) {
1867: *flg = PETSC_FALSE;
1868: PetscFunctionReturn(PETSC_SUCCESS);
1869: }
1870: PetscCall(MatDenseGetArrayRead(A1, &v1));
1871: PetscCall(MatDenseGetArrayRead(A2, &v2));
1872: for (i = 0; i < A1->cmap->n; i++) {
1873: PetscCall(PetscArraycmp(v1, v2, A1->rmap->n, flg));
1874: if (*flg == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS);
1875: v1 += mat1->lda;
1876: v2 += mat2->lda;
1877: }
1878: PetscCall(MatDenseRestoreArrayRead(A1, &v1));
1879: PetscCall(MatDenseRestoreArrayRead(A2, &v2));
1880: *flg = PETSC_TRUE;
1881: PetscFunctionReturn(PETSC_SUCCESS);
1882: }
1884: PetscErrorCode MatGetDiagonal_SeqDense(Mat A, Vec v)
1885: {
1886: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1887: PetscInt i, n, len;
1888: PetscScalar *x;
1889: const PetscScalar *vv;
1891: PetscFunctionBegin;
1892: PetscCall(VecGetSize(v, &n));
1893: PetscCall(VecGetArray(v, &x));
1894: len = PetscMin(A->rmap->n, A->cmap->n);
1895: PetscCall(MatDenseGetArrayRead(A, &vv));
1896: PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming mat and vec");
1897: for (i = 0; i < len; i++) x[i] = vv[i * mat->lda + i];
1898: PetscCall(MatDenseRestoreArrayRead(A, &vv));
1899: PetscCall(VecRestoreArray(v, &x));
1900: PetscFunctionReturn(PETSC_SUCCESS);
1901: }
1903: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A, Vec ll, Vec rr)
1904: {
1905: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1906: const PetscScalar *l, *r;
1907: PetscScalar x, *v, *vv;
1908: PetscInt i, j, m = A->rmap->n, n = A->cmap->n;
1910: PetscFunctionBegin;
1911: PetscCall(MatDenseGetArray(A, &vv));
1912: if (ll) {
1913: PetscCall(VecGetSize(ll, &m));
1914: PetscCall(VecGetArrayRead(ll, &l));
1915: PetscCheck(m == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vec wrong size");
1916: for (i = 0; i < m; i++) {
1917: x = l[i];
1918: v = vv + i;
1919: for (j = 0; j < n; j++) {
1920: (*v) *= x;
1921: v += mat->lda;
1922: }
1923: }
1924: PetscCall(VecRestoreArrayRead(ll, &l));
1925: PetscCall(PetscLogFlops(1.0 * n * m));
1926: }
1927: if (rr) {
1928: PetscCall(VecGetSize(rr, &n));
1929: PetscCall(VecGetArrayRead(rr, &r));
1930: PetscCheck(n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vec wrong size");
1931: for (i = 0; i < n; i++) {
1932: x = r[i];
1933: v = vv + i * mat->lda;
1934: for (j = 0; j < m; j++) (*v++) *= x;
1935: }
1936: PetscCall(VecRestoreArrayRead(rr, &r));
1937: PetscCall(PetscLogFlops(1.0 * n * m));
1938: }
1939: PetscCall(MatDenseRestoreArray(A, &vv));
1940: PetscFunctionReturn(PETSC_SUCCESS);
1941: }
1943: PetscErrorCode MatNorm_SeqDense(Mat A, NormType type, PetscReal *nrm)
1944: {
1945: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1946: PetscScalar *v, *vv;
1947: PetscReal sum = 0.0;
1948: PetscInt lda, m = A->rmap->n, i, j;
1950: PetscFunctionBegin;
1951: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&vv));
1952: PetscCall(MatDenseGetLDA(A, &lda));
1953: v = vv;
1954: if (type == NORM_FROBENIUS) {
1955: if (lda > m) {
1956: for (j = 0; j < A->cmap->n; j++) {
1957: v = vv + j * lda;
1958: for (i = 0; i < m; i++) {
1959: sum += PetscRealPart(PetscConj(*v) * (*v));
1960: v++;
1961: }
1962: }
1963: } else {
1964: #if defined(PETSC_USE_REAL___FP16)
1965: PetscBLASInt one = 1, cnt = A->cmap->n * A->rmap->n;
1966: PetscCallBLAS("BLASnrm2", *nrm = BLASnrm2_(&cnt, v, &one));
1967: }
1968: #else
1969: for (i = 0; i < A->cmap->n * A->rmap->n; i++) {
1970: sum += PetscRealPart(PetscConj(*v) * (*v));
1971: v++;
1972: }
1973: }
1974: *nrm = PetscSqrtReal(sum);
1975: #endif
1976: PetscCall(PetscLogFlops(2.0 * A->cmap->n * A->rmap->n));
1977: } else if (type == NORM_1) {
1978: *nrm = 0.0;
1979: for (j = 0; j < A->cmap->n; j++) {
1980: v = vv + j * mat->lda;
1981: sum = 0.0;
1982: for (i = 0; i < A->rmap->n; i++) {
1983: sum += PetscAbsScalar(*v);
1984: v++;
1985: }
1986: if (sum > *nrm) *nrm = sum;
1987: }
1988: PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
1989: } else if (type == NORM_INFINITY) {
1990: *nrm = 0.0;
1991: for (j = 0; j < A->rmap->n; j++) {
1992: v = vv + j;
1993: sum = 0.0;
1994: for (i = 0; i < A->cmap->n; i++) {
1995: sum += PetscAbsScalar(*v);
1996: v += mat->lda;
1997: }
1998: if (sum > *nrm) *nrm = sum;
1999: }
2000: PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
2001: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No two norm");
2002: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&vv));
2003: PetscFunctionReturn(PETSC_SUCCESS);
2004: }
2006: static PetscErrorCode MatSetOption_SeqDense(Mat A, MatOption op, PetscBool flg)
2007: {
2008: Mat_SeqDense *aij = (Mat_SeqDense *)A->data;
2010: PetscFunctionBegin;
2011: switch (op) {
2012: case MAT_ROW_ORIENTED:
2013: aij->roworiented = flg;
2014: break;
2015: default:
2016: break;
2017: }
2018: PetscFunctionReturn(PETSC_SUCCESS);
2019: }
2021: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
2022: {
2023: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
2024: PetscInt lda = l->lda, m = A->rmap->n, n = A->cmap->n, j;
2025: PetscScalar *v;
2027: PetscFunctionBegin;
2028: PetscCall(MatDenseGetArrayWrite(A, &v));
2029: if (lda > m) {
2030: for (j = 0; j < n; j++) PetscCall(PetscArrayzero(v + j * lda, m));
2031: } else {
2032: PetscCall(PetscArrayzero(v, PetscInt64Mult(m, n)));
2033: }
2034: PetscCall(MatDenseRestoreArrayWrite(A, &v));
2035: PetscFunctionReturn(PETSC_SUCCESS);
2036: }
2038: static PetscErrorCode MatZeroRows_SeqDense(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2039: {
2040: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
2041: PetscInt m = l->lda, n = A->cmap->n, i, j;
2042: PetscScalar *slot, *bb, *v;
2043: const PetscScalar *xx;
2045: PetscFunctionBegin;
2046: if (PetscDefined(USE_DEBUG)) {
2047: for (i = 0; i < N; i++) {
2048: PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
2049: 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);
2050: }
2051: }
2052: if (!N) PetscFunctionReturn(PETSC_SUCCESS);
2054: /* fix right-hand side if needed */
2055: if (x && b) {
2056: PetscCall(VecGetArrayRead(x, &xx));
2057: PetscCall(VecGetArray(b, &bb));
2058: for (i = 0; i < N; i++) bb[rows[i]] = diag * xx[rows[i]];
2059: PetscCall(VecRestoreArrayRead(x, &xx));
2060: PetscCall(VecRestoreArray(b, &bb));
2061: }
2063: PetscCall(MatDenseGetArray(A, &v));
2064: for (i = 0; i < N; i++) {
2065: slot = v + rows[i];
2066: for (j = 0; j < n; j++) {
2067: *slot = 0.0;
2068: slot += m;
2069: }
2070: }
2071: if (diag != 0.0) {
2072: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
2073: for (i = 0; i < N; i++) {
2074: slot = v + (m + 1) * rows[i];
2075: *slot = diag;
2076: }
2077: }
2078: PetscCall(MatDenseRestoreArray(A, &v));
2079: PetscFunctionReturn(PETSC_SUCCESS);
2080: }
2082: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A, PetscInt *lda)
2083: {
2084: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2086: PetscFunctionBegin;
2087: *lda = mat->lda;
2088: PetscFunctionReturn(PETSC_SUCCESS);
2089: }
2091: PetscErrorCode MatDenseGetArray_SeqDense(Mat A, PetscScalar **array)
2092: {
2093: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2095: PetscFunctionBegin;
2096: PetscCheck(!mat->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
2097: *array = mat->v;
2098: PetscFunctionReturn(PETSC_SUCCESS);
2099: }
2101: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A, PetscScalar **array)
2102: {
2103: PetscFunctionBegin;
2104: if (array) *array = NULL;
2105: PetscFunctionReturn(PETSC_SUCCESS);
2106: }
2108: /*@
2109: MatDenseGetLDA - gets the leading dimension of the array returned from `MatDenseGetArray()`
2111: Not Collective
2113: Input Parameter:
2114: . A - a `MATDENSE` or `MATDENSECUDA` matrix
2116: Output Parameter:
2117: . lda - the leading dimension
2119: Level: intermediate
2121: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseSetLDA()`
2122: @*/
2123: PetscErrorCode MatDenseGetLDA(Mat A, PetscInt *lda)
2124: {
2125: PetscFunctionBegin;
2127: PetscAssertPointer(lda, 2);
2128: MatCheckPreallocated(A, 1);
2129: PetscUseMethod(A, "MatDenseGetLDA_C", (Mat, PetscInt *), (A, lda));
2130: PetscFunctionReturn(PETSC_SUCCESS);
2131: }
2133: /*@
2134: MatDenseSetLDA - Sets the leading dimension of the array used by the `MATDENSE` matrix
2136: Collective if the matrix layouts have not yet been setup
2138: Input Parameters:
2139: + A - a `MATDENSE` or `MATDENSECUDA` matrix
2140: - lda - the leading dimension
2142: Level: intermediate
2144: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetLDA()`
2145: @*/
2146: PetscErrorCode MatDenseSetLDA(Mat A, PetscInt lda)
2147: {
2148: PetscFunctionBegin;
2150: PetscTryMethod(A, "MatDenseSetLDA_C", (Mat, PetscInt), (A, lda));
2151: PetscFunctionReturn(PETSC_SUCCESS);
2152: }
2154: /*@C
2155: MatDenseGetArray - gives read-write access to the array where the data for a `MATDENSE` matrix is stored
2157: Logically Collective
2159: Input Parameter:
2160: . A - a dense matrix
2162: Output Parameter:
2163: . array - pointer to the data
2165: Level: intermediate
2167: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2168: @*/
2169: PetscErrorCode MatDenseGetArray(Mat A, PetscScalar *array[]) PeNS
2170: {
2171: PetscFunctionBegin;
2173: PetscAssertPointer(array, 2);
2174: PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2175: PetscFunctionReturn(PETSC_SUCCESS);
2176: }
2178: /*@C
2179: MatDenseRestoreArray - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArray()`
2181: Logically Collective
2183: Input Parameters:
2184: + A - a dense matrix
2185: - array - pointer to the data (may be `NULL`)
2187: Level: intermediate
2189: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2190: @*/
2191: PetscErrorCode MatDenseRestoreArray(Mat A, PetscScalar *array[]) PeNS
2192: {
2193: PetscFunctionBegin;
2195: if (array) PetscAssertPointer(array, 2);
2196: PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2197: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2198: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2199: A->offloadmask = PETSC_OFFLOAD_CPU;
2200: #endif
2201: PetscFunctionReturn(PETSC_SUCCESS);
2202: }
2204: /*@C
2205: MatDenseGetArrayRead - gives read-only access to the array where the data for a `MATDENSE` matrix is stored
2207: Not Collective
2209: Input Parameter:
2210: . A - a dense matrix
2212: Output Parameter:
2213: . array - pointer to the data
2215: Level: intermediate
2217: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2218: @*/
2219: PetscErrorCode MatDenseGetArrayRead(Mat A, const PetscScalar *array[]) PeNS
2220: {
2221: PetscFunctionBegin;
2223: PetscAssertPointer(array, 2);
2224: PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2225: PetscFunctionReturn(PETSC_SUCCESS);
2226: }
2228: /*@C
2229: MatDenseRestoreArrayRead - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayRead()`
2231: Not Collective
2233: Input Parameters:
2234: + A - a dense matrix
2235: - array - pointer to the data (may be `NULL`)
2237: Level: intermediate
2239: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2240: @*/
2241: PetscErrorCode MatDenseRestoreArrayRead(Mat A, const PetscScalar *array[]) PeNS
2242: {
2243: PetscFunctionBegin;
2245: if (array) PetscAssertPointer(array, 2);
2246: PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2247: PetscFunctionReturn(PETSC_SUCCESS);
2248: }
2250: /*@C
2251: MatDenseGetArrayWrite - gives write-only access to the array where the data for a `MATDENSE` matrix is stored
2253: Not Collective
2255: Input Parameter:
2256: . A - a dense matrix
2258: Output Parameter:
2259: . array - pointer to the data
2261: Level: intermediate
2263: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2264: @*/
2265: PetscErrorCode MatDenseGetArrayWrite(Mat A, PetscScalar *array[]) PeNS
2266: {
2267: PetscFunctionBegin;
2269: PetscAssertPointer(array, 2);
2270: PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2271: PetscFunctionReturn(PETSC_SUCCESS);
2272: }
2274: /*@C
2275: MatDenseRestoreArrayWrite - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayWrite()`
2277: Not Collective
2279: Input Parameters:
2280: + A - a dense matrix
2281: - array - pointer to the data (may be `NULL`)
2283: Level: intermediate
2285: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2286: @*/
2287: PetscErrorCode MatDenseRestoreArrayWrite(Mat A, PetscScalar *array[]) PeNS
2288: {
2289: PetscFunctionBegin;
2291: if (array) PetscAssertPointer(array, 2);
2292: PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2293: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2294: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2295: A->offloadmask = PETSC_OFFLOAD_CPU;
2296: #endif
2297: PetscFunctionReturn(PETSC_SUCCESS);
2298: }
2300: /*@C
2301: MatDenseGetArrayAndMemType - gives read-write access to the array where the data for a `MATDENSE` matrix is stored
2303: Logically Collective
2305: Input Parameter:
2306: . A - a dense matrix
2308: Output Parameters:
2309: + array - pointer to the data
2310: - mtype - memory type of the returned pointer
2312: Level: intermediate
2314: Note:
2315: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2316: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2318: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArrayRead()`,
2319: `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2320: @*/
2321: PetscErrorCode MatDenseGetArrayAndMemType(Mat A, PetscScalar *array[], PetscMemType *mtype)
2322: {
2323: PetscBool isMPI;
2325: PetscFunctionBegin;
2327: PetscAssertPointer(array, 2);
2328: 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 */
2329: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2330: if (isMPI) {
2331: /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2332: PetscCall(MatDenseGetArrayAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2333: } else {
2334: PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);
2336: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayAndMemType_C", &fptr));
2337: if (fptr) {
2338: PetscCall((*fptr)(A, array, mtype));
2339: } else {
2340: PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2341: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2342: }
2343: }
2344: PetscFunctionReturn(PETSC_SUCCESS);
2345: }
2347: /*@C
2348: MatDenseRestoreArrayAndMemType - returns access to the array that is obtained by `MatDenseGetArrayAndMemType()`
2350: Logically Collective
2352: Input Parameters:
2353: + A - a dense matrix
2354: - array - pointer to the data
2356: Level: intermediate
2358: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2359: @*/
2360: PetscErrorCode MatDenseRestoreArrayAndMemType(Mat A, PetscScalar *array[])
2361: {
2362: PetscBool isMPI;
2364: PetscFunctionBegin;
2366: if (array) PetscAssertPointer(array, 2);
2367: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2368: if (isMPI) {
2369: PetscCall(MatDenseRestoreArrayAndMemType(((Mat_MPIDense *)A->data)->A, array));
2370: } else {
2371: PetscErrorCode (*fptr)(Mat, PetscScalar **);
2373: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayAndMemType_C", &fptr));
2374: if (fptr) {
2375: PetscCall((*fptr)(A, array));
2376: } else {
2377: PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2378: }
2379: if (array) *array = NULL;
2380: }
2381: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2382: PetscFunctionReturn(PETSC_SUCCESS);
2383: }
2385: /*@C
2386: MatDenseGetArrayReadAndMemType - gives read-only access to the array where the data for a `MATDENSE` matrix is stored
2388: Logically Collective
2390: Input Parameter:
2391: . A - a dense matrix
2393: Output Parameters:
2394: + array - pointer to the data
2395: - mtype - memory type of the returned pointer
2397: Level: intermediate
2399: Note:
2400: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2401: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2403: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`,
2404: `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2405: @*/
2406: PetscErrorCode MatDenseGetArrayReadAndMemType(Mat A, const PetscScalar *array[], PetscMemType *mtype)
2407: {
2408: PetscBool isMPI;
2410: PetscFunctionBegin;
2412: PetscAssertPointer(array, 2);
2413: 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 */
2414: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2415: if (isMPI) { /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2416: PetscCall(MatDenseGetArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2417: } else {
2418: PetscErrorCode (*fptr)(Mat, const PetscScalar **, PetscMemType *);
2420: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayReadAndMemType_C", &fptr));
2421: if (fptr) {
2422: PetscCall((*fptr)(A, array, mtype));
2423: } else {
2424: PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2425: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2426: }
2427: }
2428: PetscFunctionReturn(PETSC_SUCCESS);
2429: }
2431: /*@C
2432: MatDenseRestoreArrayReadAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()`
2434: Logically Collective
2436: Input Parameters:
2437: + A - a dense matrix
2438: - array - pointer to the data
2440: Level: intermediate
2442: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2443: @*/
2444: PetscErrorCode MatDenseRestoreArrayReadAndMemType(Mat A, const PetscScalar *array[])
2445: {
2446: PetscBool isMPI;
2448: PetscFunctionBegin;
2450: if (array) PetscAssertPointer(array, 2);
2451: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2452: if (isMPI) {
2453: PetscCall(MatDenseRestoreArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array));
2454: } else {
2455: PetscErrorCode (*fptr)(Mat, const PetscScalar **);
2457: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayReadAndMemType_C", &fptr));
2458: if (fptr) {
2459: PetscCall((*fptr)(A, array));
2460: } else {
2461: PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2462: }
2463: if (array) *array = NULL;
2464: }
2465: PetscFunctionReturn(PETSC_SUCCESS);
2466: }
2468: /*@C
2469: MatDenseGetArrayWriteAndMemType - gives write-only access to the array where the data for a `MATDENSE` matrix is stored
2471: Logically Collective
2473: Input Parameter:
2474: . A - a dense matrix
2476: Output Parameters:
2477: + array - pointer to the data
2478: - mtype - memory type of the returned pointer
2480: Level: intermediate
2482: Note:
2483: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2484: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2486: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWriteAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayRead()`,
2487: `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2488: @*/
2489: PetscErrorCode MatDenseGetArrayWriteAndMemType(Mat A, PetscScalar *array[], PetscMemType *mtype)
2490: {
2491: PetscBool isMPI;
2493: PetscFunctionBegin;
2495: PetscAssertPointer(array, 2);
2496: 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 */
2497: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2498: if (isMPI) {
2499: PetscCall(MatDenseGetArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2500: } else {
2501: PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);
2503: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayWriteAndMemType_C", &fptr));
2504: if (fptr) {
2505: PetscCall((*fptr)(A, array, mtype));
2506: } else {
2507: PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2508: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2509: }
2510: }
2511: PetscFunctionReturn(PETSC_SUCCESS);
2512: }
2514: /*@C
2515: MatDenseRestoreArrayWriteAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()`
2517: Logically Collective
2519: Input Parameters:
2520: + A - a dense matrix
2521: - array - pointer to the data
2523: Level: intermediate
2525: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2526: @*/
2527: PetscErrorCode MatDenseRestoreArrayWriteAndMemType(Mat A, PetscScalar *array[])
2528: {
2529: PetscBool isMPI;
2531: PetscFunctionBegin;
2533: if (array) PetscAssertPointer(array, 2);
2534: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2535: if (isMPI) {
2536: PetscCall(MatDenseRestoreArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array));
2537: } else {
2538: PetscErrorCode (*fptr)(Mat, PetscScalar **);
2540: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayWriteAndMemType_C", &fptr));
2541: if (fptr) {
2542: PetscCall((*fptr)(A, array));
2543: } else {
2544: PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2545: }
2546: if (array) *array = NULL;
2547: }
2548: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2549: PetscFunctionReturn(PETSC_SUCCESS);
2550: }
2552: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B)
2553: {
2554: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2555: PetscInt i, j, nrows, ncols, ldb;
2556: const PetscInt *irow, *icol;
2557: PetscScalar *av, *bv, *v = mat->v;
2558: Mat newmat;
2560: PetscFunctionBegin;
2561: PetscCall(ISGetIndices(isrow, &irow));
2562: PetscCall(ISGetIndices(iscol, &icol));
2563: PetscCall(ISGetLocalSize(isrow, &nrows));
2564: PetscCall(ISGetLocalSize(iscol, &ncols));
2566: /* Check submatrixcall */
2567: if (scall == MAT_REUSE_MATRIX) {
2568: PetscInt n_cols, n_rows;
2569: PetscCall(MatGetSize(*B, &n_rows, &n_cols));
2570: if (n_rows != nrows || n_cols != ncols) {
2571: /* resize the result matrix to match number of requested rows/columns */
2572: PetscCall(MatSetSizes(*B, nrows, ncols, nrows, ncols));
2573: }
2574: newmat = *B;
2575: } else {
2576: /* Create and fill new matrix */
2577: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &newmat));
2578: PetscCall(MatSetSizes(newmat, nrows, ncols, nrows, ncols));
2579: PetscCall(MatSetType(newmat, ((PetscObject)A)->type_name));
2580: PetscCall(MatSeqDenseSetPreallocation(newmat, NULL));
2581: }
2583: /* Now extract the data pointers and do the copy,column at a time */
2584: PetscCall(MatDenseGetArray(newmat, &bv));
2585: PetscCall(MatDenseGetLDA(newmat, &ldb));
2586: for (i = 0; i < ncols; i++) {
2587: av = v + mat->lda * icol[i];
2588: for (j = 0; j < nrows; j++) bv[j] = av[irow[j]];
2589: bv += ldb;
2590: }
2591: PetscCall(MatDenseRestoreArray(newmat, &bv));
2593: /* Assemble the matrices so that the correct flags are set */
2594: PetscCall(MatAssemblyBegin(newmat, MAT_FINAL_ASSEMBLY));
2595: PetscCall(MatAssemblyEnd(newmat, MAT_FINAL_ASSEMBLY));
2597: /* Free work space */
2598: PetscCall(ISRestoreIndices(isrow, &irow));
2599: PetscCall(ISRestoreIndices(iscol, &icol));
2600: *B = newmat;
2601: PetscFunctionReturn(PETSC_SUCCESS);
2602: }
2604: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[])
2605: {
2606: PetscInt i;
2608: PetscFunctionBegin;
2609: if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n, B));
2611: for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SeqDense(A, irow[i], icol[i], scall, &(*B)[i]));
2612: PetscFunctionReturn(PETSC_SUCCESS);
2613: }
2615: PetscErrorCode MatCopy_SeqDense(Mat A, Mat B, MatStructure str)
2616: {
2617: Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data;
2618: const PetscScalar *va;
2619: PetscScalar *vb;
2620: PetscInt lda1 = a->lda, lda2 = b->lda, m = A->rmap->n, n = A->cmap->n, j;
2622: PetscFunctionBegin;
2623: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2624: if (A->ops->copy != B->ops->copy) {
2625: PetscCall(MatCopy_Basic(A, B, str));
2626: PetscFunctionReturn(PETSC_SUCCESS);
2627: }
2628: PetscCheck(m == B->rmap->n && n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "size(B) != size(A)");
2629: PetscCall(MatDenseGetArrayRead(A, &va));
2630: PetscCall(MatDenseGetArray(B, &vb));
2631: if (lda1 > m || lda2 > m) {
2632: for (j = 0; j < n; j++) PetscCall(PetscArraycpy(vb + j * lda2, va + j * lda1, m));
2633: } else {
2634: PetscCall(PetscArraycpy(vb, va, A->rmap->n * A->cmap->n));
2635: }
2636: PetscCall(MatDenseRestoreArray(B, &vb));
2637: PetscCall(MatDenseRestoreArrayRead(A, &va));
2638: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2639: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2640: PetscFunctionReturn(PETSC_SUCCESS);
2641: }
2643: PetscErrorCode MatSetUp_SeqDense(Mat A)
2644: {
2645: PetscFunctionBegin;
2646: PetscCall(PetscLayoutSetUp(A->rmap));
2647: PetscCall(PetscLayoutSetUp(A->cmap));
2648: if (!A->preallocated) PetscCall(MatSeqDenseSetPreallocation(A, NULL));
2649: PetscFunctionReturn(PETSC_SUCCESS);
2650: }
2652: PetscErrorCode MatConjugate_SeqDense(Mat A)
2653: {
2654: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2655: PetscInt i, j;
2656: PetscInt min = PetscMin(A->rmap->n, A->cmap->n);
2657: PetscScalar *aa;
2659: PetscFunctionBegin;
2660: PetscCall(MatDenseGetArray(A, &aa));
2661: for (j = 0; j < A->cmap->n; j++) {
2662: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscConj(aa[i + j * mat->lda]);
2663: }
2664: PetscCall(MatDenseRestoreArray(A, &aa));
2665: if (mat->tau)
2666: for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
2667: PetscFunctionReturn(PETSC_SUCCESS);
2668: }
2670: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2671: {
2672: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2673: PetscInt i, j;
2674: PetscScalar *aa;
2676: PetscFunctionBegin;
2677: PetscCall(MatDenseGetArray(A, &aa));
2678: for (j = 0; j < A->cmap->n; j++) {
2679: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscRealPart(aa[i + j * mat->lda]);
2680: }
2681: PetscCall(MatDenseRestoreArray(A, &aa));
2682: PetscFunctionReturn(PETSC_SUCCESS);
2683: }
2685: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2686: {
2687: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2688: PetscInt i, j;
2689: PetscScalar *aa;
2691: PetscFunctionBegin;
2692: PetscCall(MatDenseGetArray(A, &aa));
2693: for (j = 0; j < A->cmap->n; j++) {
2694: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscImaginaryPart(aa[i + j * mat->lda]);
2695: }
2696: PetscCall(MatDenseRestoreArray(A, &aa));
2697: PetscFunctionReturn(PETSC_SUCCESS);
2698: }
2700: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2701: {
2702: PetscInt m = A->rmap->n, n = B->cmap->n;
2703: PetscBool cisdense = PETSC_FALSE;
2705: PetscFunctionBegin;
2706: PetscCall(MatSetSizes(C, m, n, m, n));
2707: #if defined(PETSC_HAVE_CUDA)
2708: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2709: #endif
2710: #if defined(PETSC_HAVE_HIP)
2711: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2712: #endif
2713: if (!cisdense) {
2714: PetscBool flg;
2716: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2717: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2718: }
2719: PetscCall(MatSetUp(C));
2720: PetscFunctionReturn(PETSC_SUCCESS);
2721: }
2723: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2724: {
2725: Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data, *c = (Mat_SeqDense *)C->data;
2726: PetscBLASInt m, n, k;
2727: const PetscScalar *av, *bv;
2728: PetscScalar *cv;
2729: PetscScalar _DOne = 1.0, _DZero = 0.0;
2731: PetscFunctionBegin;
2732: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2733: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2734: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2735: if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2736: PetscCall(MatDenseGetArrayRead(A, &av));
2737: PetscCall(MatDenseGetArrayRead(B, &bv));
2738: PetscCall(MatDenseGetArrayWrite(C, &cv));
2739: PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2740: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2741: PetscCall(MatDenseRestoreArrayRead(A, &av));
2742: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2743: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2744: PetscFunctionReturn(PETSC_SUCCESS);
2745: }
2747: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2748: {
2749: PetscInt m = A->rmap->n, n = B->rmap->n;
2750: PetscBool cisdense = PETSC_FALSE;
2752: PetscFunctionBegin;
2753: PetscCall(MatSetSizes(C, m, n, m, n));
2754: #if defined(PETSC_HAVE_CUDA)
2755: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2756: #endif
2757: #if defined(PETSC_HAVE_HIP)
2758: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2759: #endif
2760: if (!cisdense) {
2761: PetscBool flg;
2763: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2764: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2765: }
2766: PetscCall(MatSetUp(C));
2767: PetscFunctionReturn(PETSC_SUCCESS);
2768: }
2770: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2771: {
2772: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2773: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
2774: Mat_SeqDense *c = (Mat_SeqDense *)C->data;
2775: const PetscScalar *av, *bv;
2776: PetscScalar *cv;
2777: PetscBLASInt m, n, k;
2778: PetscScalar _DOne = 1.0, _DZero = 0.0;
2780: PetscFunctionBegin;
2781: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2782: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2783: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2784: if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2785: PetscCall(MatDenseGetArrayRead(A, &av));
2786: PetscCall(MatDenseGetArrayRead(B, &bv));
2787: PetscCall(MatDenseGetArrayWrite(C, &cv));
2788: PetscCallBLAS("BLASgemm", BLASgemm_("N", "T", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2789: PetscCall(MatDenseRestoreArrayRead(A, &av));
2790: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2791: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2792: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2793: PetscFunctionReturn(PETSC_SUCCESS);
2794: }
2796: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2797: {
2798: PetscInt m = A->cmap->n, n = B->cmap->n;
2799: PetscBool cisdense = PETSC_FALSE;
2801: PetscFunctionBegin;
2802: PetscCall(MatSetSizes(C, m, n, m, n));
2803: #if defined(PETSC_HAVE_CUDA)
2804: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2805: #endif
2806: #if defined(PETSC_HAVE_HIP)
2807: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2808: #endif
2809: if (!cisdense) {
2810: PetscBool flg;
2812: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2813: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2814: }
2815: PetscCall(MatSetUp(C));
2816: PetscFunctionReturn(PETSC_SUCCESS);
2817: }
2819: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2820: {
2821: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2822: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
2823: Mat_SeqDense *c = (Mat_SeqDense *)C->data;
2824: const PetscScalar *av, *bv;
2825: PetscScalar *cv;
2826: PetscBLASInt m, n, k;
2827: PetscScalar _DOne = 1.0, _DZero = 0.0;
2829: PetscFunctionBegin;
2830: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2831: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2832: PetscCall(PetscBLASIntCast(A->rmap->n, &k));
2833: if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2834: PetscCall(MatDenseGetArrayRead(A, &av));
2835: PetscCall(MatDenseGetArrayRead(B, &bv));
2836: PetscCall(MatDenseGetArrayWrite(C, &cv));
2837: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2838: PetscCall(MatDenseRestoreArrayRead(A, &av));
2839: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2840: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2841: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2842: PetscFunctionReturn(PETSC_SUCCESS);
2843: }
2845: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2846: {
2847: PetscFunctionBegin;
2848: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2849: C->ops->productsymbolic = MatProductSymbolic_AB;
2850: PetscFunctionReturn(PETSC_SUCCESS);
2851: }
2853: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2854: {
2855: PetscFunctionBegin;
2856: C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2857: C->ops->productsymbolic = MatProductSymbolic_AtB;
2858: PetscFunctionReturn(PETSC_SUCCESS);
2859: }
2861: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2862: {
2863: PetscFunctionBegin;
2864: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2865: C->ops->productsymbolic = MatProductSymbolic_ABt;
2866: PetscFunctionReturn(PETSC_SUCCESS);
2867: }
2869: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2870: {
2871: Mat_Product *product = C->product;
2873: PetscFunctionBegin;
2874: switch (product->type) {
2875: case MATPRODUCT_AB:
2876: PetscCall(MatProductSetFromOptions_SeqDense_AB(C));
2877: break;
2878: case MATPRODUCT_AtB:
2879: PetscCall(MatProductSetFromOptions_SeqDense_AtB(C));
2880: break;
2881: case MATPRODUCT_ABt:
2882: PetscCall(MatProductSetFromOptions_SeqDense_ABt(C));
2883: break;
2884: default:
2885: break;
2886: }
2887: PetscFunctionReturn(PETSC_SUCCESS);
2888: }
2890: static PetscErrorCode MatGetRowMax_SeqDense(Mat A, Vec v, PetscInt idx[])
2891: {
2892: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2893: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2894: PetscScalar *x;
2895: const PetscScalar *aa;
2897: PetscFunctionBegin;
2898: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2899: PetscCall(VecGetArray(v, &x));
2900: PetscCall(VecGetLocalSize(v, &p));
2901: PetscCall(MatDenseGetArrayRead(A, &aa));
2902: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2903: for (i = 0; i < m; i++) {
2904: x[i] = aa[i];
2905: if (idx) idx[i] = 0;
2906: for (j = 1; j < n; j++) {
2907: if (PetscRealPart(x[i]) < PetscRealPart(aa[i + a->lda * j])) {
2908: x[i] = aa[i + a->lda * j];
2909: if (idx) idx[i] = j;
2910: }
2911: }
2912: }
2913: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2914: PetscCall(VecRestoreArray(v, &x));
2915: PetscFunctionReturn(PETSC_SUCCESS);
2916: }
2918: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A, Vec v, PetscInt idx[])
2919: {
2920: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2921: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2922: PetscScalar *x;
2923: PetscReal atmp;
2924: const PetscScalar *aa;
2926: PetscFunctionBegin;
2927: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2928: PetscCall(VecGetArray(v, &x));
2929: PetscCall(VecGetLocalSize(v, &p));
2930: PetscCall(MatDenseGetArrayRead(A, &aa));
2931: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2932: for (i = 0; i < m; i++) {
2933: x[i] = PetscAbsScalar(aa[i]);
2934: for (j = 1; j < n; j++) {
2935: atmp = PetscAbsScalar(aa[i + a->lda * j]);
2936: if (PetscAbsScalar(x[i]) < atmp) {
2937: x[i] = atmp;
2938: if (idx) idx[i] = j;
2939: }
2940: }
2941: }
2942: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2943: PetscCall(VecRestoreArray(v, &x));
2944: PetscFunctionReturn(PETSC_SUCCESS);
2945: }
2947: static PetscErrorCode MatGetRowMin_SeqDense(Mat A, Vec v, PetscInt idx[])
2948: {
2949: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2950: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2951: PetscScalar *x;
2952: const PetscScalar *aa;
2954: PetscFunctionBegin;
2955: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2956: PetscCall(MatDenseGetArrayRead(A, &aa));
2957: PetscCall(VecGetArray(v, &x));
2958: PetscCall(VecGetLocalSize(v, &p));
2959: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2960: for (i = 0; i < m; i++) {
2961: x[i] = aa[i];
2962: if (idx) idx[i] = 0;
2963: for (j = 1; j < n; j++) {
2964: if (PetscRealPart(x[i]) > PetscRealPart(aa[i + a->lda * j])) {
2965: x[i] = aa[i + a->lda * j];
2966: if (idx) idx[i] = j;
2967: }
2968: }
2969: }
2970: PetscCall(VecRestoreArray(v, &x));
2971: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2972: PetscFunctionReturn(PETSC_SUCCESS);
2973: }
2975: PetscErrorCode MatGetColumnVector_SeqDense(Mat A, Vec v, PetscInt col)
2976: {
2977: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2978: PetscScalar *x;
2979: const PetscScalar *aa;
2981: PetscFunctionBegin;
2982: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2983: PetscCall(MatDenseGetArrayRead(A, &aa));
2984: PetscCall(VecGetArray(v, &x));
2985: PetscCall(PetscArraycpy(x, aa + col * a->lda, A->rmap->n));
2986: PetscCall(VecRestoreArray(v, &x));
2987: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2988: PetscFunctionReturn(PETSC_SUCCESS);
2989: }
2991: PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A, PetscInt type, PetscReal *reductions)
2992: {
2993: PetscInt i, j, m, n;
2994: const PetscScalar *a;
2996: PetscFunctionBegin;
2997: PetscCall(MatGetSize(A, &m, &n));
2998: PetscCall(PetscArrayzero(reductions, n));
2999: PetscCall(MatDenseGetArrayRead(A, &a));
3000: if (type == NORM_2) {
3001: for (i = 0; i < n; i++) {
3002: for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j] * a[j]);
3003: a = PetscSafePointerPlusOffset(a, m);
3004: }
3005: } else if (type == NORM_1) {
3006: for (i = 0; i < n; i++) {
3007: for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j]);
3008: a = PetscSafePointerPlusOffset(a, m);
3009: }
3010: } else if (type == NORM_INFINITY) {
3011: for (i = 0; i < n; i++) {
3012: for (j = 0; j < m; j++) reductions[i] = PetscMax(PetscAbsScalar(a[j]), reductions[i]);
3013: a = PetscSafePointerPlusOffset(a, m);
3014: }
3015: } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
3016: for (i = 0; i < n; i++) {
3017: for (j = 0; j < m; j++) reductions[i] += PetscRealPart(a[j]);
3018: a = PetscSafePointerPlusOffset(a, m);
3019: }
3020: } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3021: for (i = 0; i < n; i++) {
3022: for (j = 0; j < m; j++) reductions[i] += PetscImaginaryPart(a[j]);
3023: a = PetscSafePointerPlusOffset(a, m);
3024: }
3025: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown reduction type");
3026: PetscCall(MatDenseRestoreArrayRead(A, &a));
3027: if (type == NORM_2) {
3028: for (i = 0; i < n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
3029: } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3030: for (i = 0; i < n; i++) reductions[i] /= m;
3031: }
3032: PetscFunctionReturn(PETSC_SUCCESS);
3033: }
3035: PetscErrorCode MatSetRandom_SeqDense(Mat x, PetscRandom rctx)
3036: {
3037: PetscScalar *a;
3038: PetscInt lda, m, n, i, j;
3040: PetscFunctionBegin;
3041: PetscCall(MatGetSize(x, &m, &n));
3042: PetscCall(MatDenseGetLDA(x, &lda));
3043: PetscCall(MatDenseGetArrayWrite(x, &a));
3044: for (j = 0; j < n; j++) {
3045: for (i = 0; i < m; i++) PetscCall(PetscRandomGetValue(rctx, a + j * lda + i));
3046: }
3047: PetscCall(MatDenseRestoreArrayWrite(x, &a));
3048: PetscFunctionReturn(PETSC_SUCCESS);
3049: }
3051: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A, PetscBool *missing, PetscInt *d)
3052: {
3053: PetscFunctionBegin;
3054: *missing = PETSC_FALSE;
3055: PetscFunctionReturn(PETSC_SUCCESS);
3056: }
3058: /* vals is not const */
3059: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A, PetscInt col, PetscScalar **vals)
3060: {
3061: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3062: PetscScalar *v;
3064: PetscFunctionBegin;
3065: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3066: PetscCall(MatDenseGetArray(A, &v));
3067: *vals = v + col * a->lda;
3068: PetscCall(MatDenseRestoreArray(A, &v));
3069: PetscFunctionReturn(PETSC_SUCCESS);
3070: }
3072: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A, PetscScalar **vals)
3073: {
3074: PetscFunctionBegin;
3075: if (vals) *vals = NULL; /* user cannot accidentally use the array later */
3076: PetscFunctionReturn(PETSC_SUCCESS);
3077: }
3079: static struct _MatOps MatOps_Values = {MatSetValues_SeqDense,
3080: MatGetRow_SeqDense,
3081: MatRestoreRow_SeqDense,
3082: MatMult_SeqDense,
3083: /* 4*/ MatMultAdd_SeqDense,
3084: MatMultTranspose_SeqDense,
3085: MatMultTransposeAdd_SeqDense,
3086: NULL,
3087: NULL,
3088: NULL,
3089: /* 10*/ NULL,
3090: MatLUFactor_SeqDense,
3091: MatCholeskyFactor_SeqDense,
3092: MatSOR_SeqDense,
3093: MatTranspose_SeqDense,
3094: /* 15*/ MatGetInfo_SeqDense,
3095: MatEqual_SeqDense,
3096: MatGetDiagonal_SeqDense,
3097: MatDiagonalScale_SeqDense,
3098: MatNorm_SeqDense,
3099: /* 20*/ NULL,
3100: NULL,
3101: MatSetOption_SeqDense,
3102: MatZeroEntries_SeqDense,
3103: /* 24*/ MatZeroRows_SeqDense,
3104: NULL,
3105: NULL,
3106: NULL,
3107: NULL,
3108: /* 29*/ MatSetUp_SeqDense,
3109: NULL,
3110: NULL,
3111: NULL,
3112: NULL,
3113: /* 34*/ MatDuplicate_SeqDense,
3114: NULL,
3115: NULL,
3116: NULL,
3117: NULL,
3118: /* 39*/ MatAXPY_SeqDense,
3119: MatCreateSubMatrices_SeqDense,
3120: NULL,
3121: MatGetValues_SeqDense,
3122: MatCopy_SeqDense,
3123: /* 44*/ MatGetRowMax_SeqDense,
3124: MatScale_SeqDense,
3125: MatShift_SeqDense,
3126: NULL,
3127: MatZeroRowsColumns_SeqDense,
3128: /* 49*/ MatSetRandom_SeqDense,
3129: NULL,
3130: NULL,
3131: NULL,
3132: NULL,
3133: /* 54*/ NULL,
3134: NULL,
3135: NULL,
3136: NULL,
3137: NULL,
3138: /* 59*/ MatCreateSubMatrix_SeqDense,
3139: MatDestroy_SeqDense,
3140: MatView_SeqDense,
3141: NULL,
3142: NULL,
3143: /* 64*/ NULL,
3144: NULL,
3145: NULL,
3146: NULL,
3147: MatGetRowMaxAbs_SeqDense,
3148: /* 69*/ NULL,
3149: NULL,
3150: NULL,
3151: NULL,
3152: NULL,
3153: /* 74*/ NULL,
3154: NULL,
3155: NULL,
3156: NULL,
3157: MatLoad_SeqDense,
3158: /* 79*/ MatIsSymmetric_SeqDense,
3159: MatIsHermitian_SeqDense,
3160: NULL,
3161: NULL,
3162: NULL,
3163: /* 84*/ NULL,
3164: MatMatMultNumeric_SeqDense_SeqDense,
3165: NULL,
3166: NULL,
3167: MatMatTransposeMultNumeric_SeqDense_SeqDense,
3168: /* 89*/ NULL,
3169: MatProductSetFromOptions_SeqDense,
3170: NULL,
3171: NULL,
3172: MatConjugate_SeqDense,
3173: /* 94*/ NULL,
3174: NULL,
3175: MatRealPart_SeqDense,
3176: MatImaginaryPart_SeqDense,
3177: NULL,
3178: /* 99*/ NULL,
3179: NULL,
3180: NULL,
3181: MatGetRowMin_SeqDense,
3182: MatGetColumnVector_SeqDense,
3183: /*104*/ MatMissingDiagonal_SeqDense,
3184: NULL,
3185: NULL,
3186: NULL,
3187: NULL,
3188: /*109*/ NULL,
3189: NULL,
3190: NULL,
3191: MatMultHermitianTranspose_SeqDense,
3192: MatMultHermitianTransposeAdd_SeqDense,
3193: /*114*/ NULL,
3194: NULL,
3195: MatGetColumnReductions_SeqDense,
3196: NULL,
3197: NULL,
3198: /*119*/ NULL,
3199: NULL,
3200: NULL,
3201: MatTransposeMatMultNumeric_SeqDense_SeqDense,
3202: NULL,
3203: /*124*/ NULL,
3204: NULL,
3205: NULL,
3206: NULL,
3207: NULL,
3208: /*129*/ NULL,
3209: NULL,
3210: MatCreateMPIMatConcatenateSeqMat_SeqDense,
3211: NULL,
3212: NULL,
3213: /*134*/ NULL,
3214: NULL,
3215: NULL,
3216: NULL,
3217: NULL,
3218: /*139*/ NULL,
3219: NULL,
3220: NULL,
3221: NULL,
3222: NULL};
3224: /*@
3225: MatCreateSeqDense - Creates a `MATSEQDENSE` that
3226: is stored in column major order (the usual Fortran format).
3228: Collective
3230: Input Parameters:
3231: + comm - MPI communicator, set to `PETSC_COMM_SELF`
3232: . m - number of rows
3233: . n - number of columns
3234: - data - optional location of matrix data in column major order. Use `NULL` for PETSc
3235: to control all matrix memory allocation.
3237: Output Parameter:
3238: . A - the matrix
3240: Level: intermediate
3242: Note:
3243: The data input variable is intended primarily for Fortran programmers
3244: who wish to allocate their own matrix memory space. Most users should
3245: set `data` = `NULL`.
3247: Developer Note:
3248: Many of the matrix operations for this variant use the BLAS and LAPACK routines.
3250: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`
3251: @*/
3252: PetscErrorCode MatCreateSeqDense(MPI_Comm comm, PetscInt m, PetscInt n, PetscScalar data[], Mat *A)
3253: {
3254: PetscFunctionBegin;
3255: PetscCall(MatCreate(comm, A));
3256: PetscCall(MatSetSizes(*A, m, n, m, n));
3257: PetscCall(MatSetType(*A, MATSEQDENSE));
3258: PetscCall(MatSeqDenseSetPreallocation(*A, data));
3259: PetscFunctionReturn(PETSC_SUCCESS);
3260: }
3262: /*@
3263: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements of a `MATSEQDENSE` matrix
3265: Collective
3267: Input Parameters:
3268: + B - the matrix
3269: - data - the array (or `NULL`)
3271: Level: intermediate
3273: Note:
3274: The data input variable is intended primarily for Fortran programmers
3275: who wish to allocate their own matrix memory space. Most users should
3276: need not call this routine.
3278: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`, `MatDenseSetLDA()`
3279: @*/
3280: PetscErrorCode MatSeqDenseSetPreallocation(Mat B, PetscScalar data[])
3281: {
3282: PetscFunctionBegin;
3284: PetscTryMethod(B, "MatSeqDenseSetPreallocation_C", (Mat, PetscScalar[]), (B, data));
3285: PetscFunctionReturn(PETSC_SUCCESS);
3286: }
3288: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B, PetscScalar *data)
3289: {
3290: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
3292: PetscFunctionBegin;
3293: PetscCheck(!b->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3294: B->preallocated = PETSC_TRUE;
3296: PetscCall(PetscLayoutSetUp(B->rmap));
3297: PetscCall(PetscLayoutSetUp(B->cmap));
3299: if (b->lda <= 0) PetscCall(PetscBLASIntCast(B->rmap->n, &b->lda));
3301: if (!data) { /* petsc-allocated storage */
3302: if (!b->user_alloc) PetscCall(PetscFree(b->v));
3303: PetscCall(PetscCalloc1((size_t)b->lda * B->cmap->n, &b->v));
3305: b->user_alloc = PETSC_FALSE;
3306: } else { /* user-allocated storage */
3307: if (!b->user_alloc) PetscCall(PetscFree(b->v));
3308: b->v = data;
3309: b->user_alloc = PETSC_TRUE;
3310: }
3311: B->assembled = PETSC_TRUE;
3312: PetscFunctionReturn(PETSC_SUCCESS);
3313: }
3315: #if defined(PETSC_HAVE_ELEMENTAL)
3316: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
3317: {
3318: Mat mat_elemental;
3319: const PetscScalar *array;
3320: PetscScalar *v_colwise;
3321: PetscInt M = A->rmap->N, N = A->cmap->N, i, j, k, *rows, *cols;
3323: PetscFunctionBegin;
3324: PetscCall(PetscMalloc3(M * N, &v_colwise, M, &rows, N, &cols));
3325: PetscCall(MatDenseGetArrayRead(A, &array));
3326: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
3327: k = 0;
3328: for (j = 0; j < N; j++) {
3329: cols[j] = j;
3330: for (i = 0; i < M; i++) v_colwise[j * M + i] = array[k++];
3331: }
3332: for (i = 0; i < M; i++) rows[i] = i;
3333: PetscCall(MatDenseRestoreArrayRead(A, &array));
3335: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
3336: PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
3337: PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
3338: PetscCall(MatSetUp(mat_elemental));
3340: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
3341: PetscCall(MatSetValues(mat_elemental, M, rows, N, cols, v_colwise, ADD_VALUES));
3342: PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
3343: PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
3344: PetscCall(PetscFree3(v_colwise, rows, cols));
3346: if (reuse == MAT_INPLACE_MATRIX) {
3347: PetscCall(MatHeaderReplace(A, &mat_elemental));
3348: } else {
3349: *newmat = mat_elemental;
3350: }
3351: PetscFunctionReturn(PETSC_SUCCESS);
3352: }
3353: #endif
3355: PetscErrorCode MatDenseSetLDA_SeqDense(Mat B, PetscInt lda)
3356: {
3357: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
3358: PetscBool data;
3360: PetscFunctionBegin;
3361: data = (B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE;
3362: PetscCheck(b->user_alloc || !data || b->lda == lda, PETSC_COMM_SELF, PETSC_ERR_ORDER, "LDA cannot be changed after allocation of internal storage");
3363: 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);
3364: PetscCall(PetscBLASIntCast(lda, &b->lda));
3365: PetscFunctionReturn(PETSC_SUCCESS);
3366: }
3368: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
3369: {
3370: PetscFunctionBegin;
3371: PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIDense(comm, inmat, n, scall, outmat));
3372: PetscFunctionReturn(PETSC_SUCCESS);
3373: }
3375: PetscErrorCode MatDenseCreateColumnVec_Private(Mat A, Vec *v)
3376: {
3377: PetscBool isstd, iskok, iscuda, iship;
3378: PetscMPIInt size;
3379: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
3380: /* we pass the data of A, to prevent allocating needless GPU memory the first time VecCUPMPlaceArray is called. */
3381: const PetscScalar *a;
3382: #endif
3384: PetscFunctionBegin;
3385: *v = NULL;
3386: PetscCall(PetscStrcmpAny(A->defaultvectype, &isstd, VECSTANDARD, VECSEQ, VECMPI, ""));
3387: PetscCall(PetscStrcmpAny(A->defaultvectype, &iskok, VECKOKKOS, VECSEQKOKKOS, VECMPIKOKKOS, ""));
3388: PetscCall(PetscStrcmpAny(A->defaultvectype, &iscuda, VECCUDA, VECSEQCUDA, VECMPICUDA, ""));
3389: PetscCall(PetscStrcmpAny(A->defaultvectype, &iship, VECHIP, VECSEQHIP, VECMPIHIP, ""));
3390: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
3391: if (isstd) {
3392: if (size > 1) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, v));
3393: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, v));
3394: } else if (iskok) {
3395: PetscCheck(PetscDefined(HAVE_KOKKOS_KERNELS), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using KOKKOS kernels support");
3396: #if PetscDefined(HAVE_KOKKOS_KERNELS)
3397: if (size > 1) PetscCall(VecCreateMPIKokkosWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, v));
3398: else PetscCall(VecCreateSeqKokkosWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, v));
3399: #endif
3400: } else if (iscuda) {
3401: PetscCheck(PetscDefined(HAVE_CUDA), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using CUDA support");
3402: #if PetscDefined(HAVE_CUDA)
3403: PetscCall(MatDenseCUDAGetArrayRead(A, &a));
3404: if (size > 1) PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, a, v));
3405: else PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, a, v));
3406: #endif
3407: } else if (iship) {
3408: PetscCheck(PetscDefined(HAVE_HIP), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using HIP support");
3409: #if PetscDefined(HAVE_HIP)
3410: PetscCall(MatDenseHIPGetArrayRead(A, &a));
3411: if (size > 1) PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, a, v));
3412: else PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, a, v));
3413: #endif
3414: }
3415: PetscCheck(*v, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not coded for type %s", A->defaultvectype);
3416: PetscFunctionReturn(PETSC_SUCCESS);
3417: }
3419: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3420: {
3421: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3423: PetscFunctionBegin;
3424: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3425: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3426: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3427: a->vecinuse = col + 1;
3428: PetscCall(MatDenseGetArray(A, (PetscScalar **)&a->ptrinuse));
3429: PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda));
3430: *v = a->cvec;
3431: PetscFunctionReturn(PETSC_SUCCESS);
3432: }
3434: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3435: {
3436: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3438: PetscFunctionBegin;
3439: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3440: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3441: VecCheckAssembled(a->cvec);
3442: a->vecinuse = 0;
3443: PetscCall(MatDenseRestoreArray(A, (PetscScalar **)&a->ptrinuse));
3444: PetscCall(VecResetArray(a->cvec));
3445: if (v) *v = NULL;
3446: PetscFunctionReturn(PETSC_SUCCESS);
3447: }
3449: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3450: {
3451: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3453: PetscFunctionBegin;
3454: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3455: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3456: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3457: a->vecinuse = col + 1;
3458: PetscCall(MatDenseGetArrayRead(A, &a->ptrinuse));
3459: PetscCall(VecPlaceArray(a->cvec, PetscSafePointerPlusOffset(a->ptrinuse, (size_t)col * (size_t)a->lda)));
3460: PetscCall(VecLockReadPush(a->cvec));
3461: *v = a->cvec;
3462: PetscFunctionReturn(PETSC_SUCCESS);
3463: }
3465: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3466: {
3467: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3469: PetscFunctionBegin;
3470: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3471: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3472: VecCheckAssembled(a->cvec);
3473: a->vecinuse = 0;
3474: PetscCall(MatDenseRestoreArrayRead(A, &a->ptrinuse));
3475: PetscCall(VecLockReadPop(a->cvec));
3476: PetscCall(VecResetArray(a->cvec));
3477: if (v) *v = NULL;
3478: PetscFunctionReturn(PETSC_SUCCESS);
3479: }
3481: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3482: {
3483: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3485: PetscFunctionBegin;
3486: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3487: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3488: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3489: a->vecinuse = col + 1;
3490: PetscCall(MatDenseGetArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3491: PetscCall(VecPlaceArray(a->cvec, PetscSafePointerPlusOffset(a->ptrinuse, (size_t)col * (size_t)a->lda)));
3492: *v = a->cvec;
3493: PetscFunctionReturn(PETSC_SUCCESS);
3494: }
3496: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3497: {
3498: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3500: PetscFunctionBegin;
3501: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3502: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3503: VecCheckAssembled(a->cvec);
3504: a->vecinuse = 0;
3505: PetscCall(MatDenseRestoreArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3506: PetscCall(VecResetArray(a->cvec));
3507: if (v) *v = NULL;
3508: PetscFunctionReturn(PETSC_SUCCESS);
3509: }
3511: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3512: {
3513: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3515: PetscFunctionBegin;
3516: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3517: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3518: if (a->cmat && (cend - cbegin != a->cmat->cmap->N || rend - rbegin != a->cmat->rmap->N)) PetscCall(MatDestroy(&a->cmat));
3519: if (!a->cmat) {
3520: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), rend - rbegin, PETSC_DECIDE, rend - rbegin, cend - cbegin, PetscSafePointerPlusOffset(a->v, rbegin + (size_t)cbegin * a->lda), &a->cmat));
3521: } else {
3522: PetscCall(MatDensePlaceArray(a->cmat, PetscSafePointerPlusOffset(a->v, rbegin + (size_t)cbegin * a->lda)));
3523: }
3524: PetscCall(MatDenseSetLDA(a->cmat, a->lda));
3525: a->matinuse = cbegin + 1;
3526: *v = a->cmat;
3527: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
3528: A->offloadmask = PETSC_OFFLOAD_CPU;
3529: #endif
3530: PetscFunctionReturn(PETSC_SUCCESS);
3531: }
3533: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A, Mat *v)
3534: {
3535: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3537: PetscFunctionBegin;
3538: PetscCheck(a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetSubMatrix() first");
3539: PetscCheck(a->cmat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column matrix");
3540: PetscCheck(*v == a->cmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not the matrix obtained from MatDenseGetSubMatrix()");
3541: a->matinuse = 0;
3542: PetscCall(MatDenseResetArray(a->cmat));
3543: if (v) *v = NULL;
3544: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
3545: A->offloadmask = PETSC_OFFLOAD_CPU;
3546: #endif
3547: PetscFunctionReturn(PETSC_SUCCESS);
3548: }
3550: /*MC
3551: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
3553: Options Database Key:
3554: . -mat_type seqdense - sets the matrix type to `MATSEQDENSE` during a call to `MatSetFromOptions()`
3556: Level: beginner
3558: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreateSeqDense()`
3559: M*/
3560: PetscErrorCode MatCreate_SeqDense(Mat B)
3561: {
3562: Mat_SeqDense *b;
3563: PetscMPIInt size;
3565: PetscFunctionBegin;
3566: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
3567: PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1");
3569: PetscCall(PetscNew(&b));
3570: B->data = (void *)b;
3571: B->ops[0] = MatOps_Values;
3573: b->roworiented = PETSC_TRUE;
3575: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactor_C", MatQRFactor_SeqDense));
3576: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetLDA_C", MatDenseGetLDA_SeqDense));
3577: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseSetLDA_C", MatDenseSetLDA_SeqDense));
3578: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArray_C", MatDenseGetArray_SeqDense));
3579: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArray_C", MatDenseRestoreArray_SeqDense));
3580: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDensePlaceArray_C", MatDensePlaceArray_SeqDense));
3581: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseResetArray_C", MatDenseResetArray_SeqDense));
3582: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseReplaceArray_C", MatDenseReplaceArray_SeqDense));
3583: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayRead_C", MatDenseGetArray_SeqDense));
3584: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayRead_C", MatDenseRestoreArray_SeqDense));
3585: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayWrite_C", MatDenseGetArray_SeqDense));
3586: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayWrite_C", MatDenseRestoreArray_SeqDense));
3587: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqaij_C", MatConvert_SeqDense_SeqAIJ));
3588: #if defined(PETSC_HAVE_ELEMENTAL)
3589: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_elemental_C", MatConvert_SeqDense_Elemental));
3590: #endif
3591: #if defined(PETSC_HAVE_SCALAPACK)
3592: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_scalapack_C", MatConvert_Dense_ScaLAPACK));
3593: #endif
3594: #if defined(PETSC_HAVE_CUDA)
3595: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensecuda_C", MatConvert_SeqDense_SeqDenseCUDA));
3596: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3597: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdense_C", MatProductSetFromOptions_SeqDense));
3598: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3599: #endif
3600: #if defined(PETSC_HAVE_HIP)
3601: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensehip_C", MatConvert_SeqDense_SeqDenseHIP));
3602: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3603: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdense_C", MatProductSetFromOptions_SeqDense));
3604: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3605: #endif
3606: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqDenseSetPreallocation_C", MatSeqDenseSetPreallocation_SeqDense));
3607: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqdense_C", MatProductSetFromOptions_SeqAIJ_SeqDense));
3608: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdense_C", MatProductSetFromOptions_SeqDense));
3609: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));
3610: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqsbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));
3612: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumn_C", MatDenseGetColumn_SeqDense));
3613: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumn_C", MatDenseRestoreColumn_SeqDense));
3614: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVec_C", MatDenseGetColumnVec_SeqDense));
3615: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVec_C", MatDenseRestoreColumnVec_SeqDense));
3616: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecRead_C", MatDenseGetColumnVecRead_SeqDense));
3617: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecRead_C", MatDenseRestoreColumnVecRead_SeqDense));
3618: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecWrite_C", MatDenseGetColumnVecWrite_SeqDense));
3619: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecWrite_C", MatDenseRestoreColumnVecWrite_SeqDense));
3620: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetSubMatrix_C", MatDenseGetSubMatrix_SeqDense));
3621: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreSubMatrix_C", MatDenseRestoreSubMatrix_SeqDense));
3622: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultColumnRange_C", MatMultColumnRange_SeqDense));
3623: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultAddColumnRange_C", MatMultAddColumnRange_SeqDense));
3624: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultHermitianTransposeColumnRange_C", MatMultHermitianTransposeColumnRange_SeqDense));
3625: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultHermitianTransposeAddColumnRange_C", MatMultHermitianTransposeAddColumnRange_SeqDense));
3626: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQDENSE));
3627: PetscFunctionReturn(PETSC_SUCCESS);
3628: }
3630: /*@C
3631: 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.
3633: Not Collective
3635: Input Parameters:
3636: + A - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3637: - col - column index
3639: Output Parameter:
3640: . vals - pointer to the data
3642: Level: intermediate
3644: Note:
3645: Use `MatDenseGetColumnVec()` to get access to a column of a `MATDENSE` treated as a `Vec`
3647: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreColumn()`, `MatDenseGetColumnVec()`
3648: @*/
3649: PetscErrorCode MatDenseGetColumn(Mat A, PetscInt col, PetscScalar *vals[])
3650: {
3651: PetscFunctionBegin;
3654: PetscAssertPointer(vals, 3);
3655: PetscUseMethod(A, "MatDenseGetColumn_C", (Mat, PetscInt, PetscScalar **), (A, col, vals));
3656: PetscFunctionReturn(PETSC_SUCCESS);
3657: }
3659: /*@C
3660: MatDenseRestoreColumn - returns access to a column of a `MATDENSE` matrix which is returned by `MatDenseGetColumn()`.
3662: Not Collective
3664: Input Parameters:
3665: + A - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3666: - vals - pointer to the data (may be `NULL`)
3668: Level: intermediate
3670: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetColumn()`
3671: @*/
3672: PetscErrorCode MatDenseRestoreColumn(Mat A, PetscScalar *vals[])
3673: {
3674: PetscFunctionBegin;
3676: PetscAssertPointer(vals, 2);
3677: PetscUseMethod(A, "MatDenseRestoreColumn_C", (Mat, PetscScalar **), (A, vals));
3678: PetscFunctionReturn(PETSC_SUCCESS);
3679: }
3681: /*@
3682: MatDenseGetColumnVec - Gives read-write access to a column of a `MATDENSE` matrix, represented as a `Vec`.
3684: Collective
3686: Input Parameters:
3687: + A - the `Mat` object
3688: - col - the column index
3690: Output Parameter:
3691: . v - the vector
3693: Level: intermediate
3695: Notes:
3696: The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVec()` when the vector is no longer needed.
3698: Use `MatDenseGetColumnVecRead()` to obtain read-only access or `MatDenseGetColumnVecWrite()` for write-only access.
3700: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`, `MatDenseGetColumn()`
3701: @*/
3702: PetscErrorCode MatDenseGetColumnVec(Mat A, PetscInt col, Vec *v)
3703: {
3704: PetscFunctionBegin;
3708: PetscAssertPointer(v, 3);
3709: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3710: 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);
3711: PetscUseMethod(A, "MatDenseGetColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3712: PetscFunctionReturn(PETSC_SUCCESS);
3713: }
3715: /*@
3716: MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVec()`.
3718: Collective
3720: Input Parameters:
3721: + A - the `Mat` object
3722: . col - the column index
3723: - v - the `Vec` object (may be `NULL`)
3725: Level: intermediate
3727: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3728: @*/
3729: PetscErrorCode MatDenseRestoreColumnVec(Mat A, PetscInt col, Vec *v)
3730: {
3731: PetscFunctionBegin;
3735: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3736: 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);
3737: PetscUseMethod(A, "MatDenseRestoreColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3738: PetscFunctionReturn(PETSC_SUCCESS);
3739: }
3741: /*@
3742: MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a `Vec`.
3744: Collective
3746: Input Parameters:
3747: + A - the `Mat` object
3748: - col - the column index
3750: Output Parameter:
3751: . v - the vector
3753: Level: intermediate
3755: Notes:
3756: The vector is owned by PETSc and users cannot modify it.
3758: Users need to call `MatDenseRestoreColumnVecRead()` when the vector is no longer needed.
3760: Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecWrite()` for write-only access.
3762: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3763: @*/
3764: PetscErrorCode MatDenseGetColumnVecRead(Mat A, PetscInt col, Vec *v)
3765: {
3766: PetscFunctionBegin;
3770: PetscAssertPointer(v, 3);
3771: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3772: 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);
3773: PetscUseMethod(A, "MatDenseGetColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3774: PetscFunctionReturn(PETSC_SUCCESS);
3775: }
3777: /*@
3778: MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVecRead()`.
3780: Collective
3782: Input Parameters:
3783: + A - the `Mat` object
3784: . col - the column index
3785: - v - the `Vec` object (may be `NULL`)
3787: Level: intermediate
3789: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecWrite()`
3790: @*/
3791: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A, PetscInt col, Vec *v)
3792: {
3793: PetscFunctionBegin;
3797: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3798: PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N);
3799: PetscUseMethod(A, "MatDenseRestoreColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3800: PetscFunctionReturn(PETSC_SUCCESS);
3801: }
3803: /*@
3804: MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a `Vec`.
3806: Collective
3808: Input Parameters:
3809: + A - the `Mat` object
3810: - col - the column index
3812: Output Parameter:
3813: . v - the vector
3815: Level: intermediate
3817: Notes:
3818: The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVecWrite()` when the vector is no longer needed.
3820: Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecRead()` for read-only access.
3822: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3823: @*/
3824: PetscErrorCode MatDenseGetColumnVecWrite(Mat A, PetscInt col, Vec *v)
3825: {
3826: PetscFunctionBegin;
3830: PetscAssertPointer(v, 3);
3831: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3832: 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);
3833: PetscUseMethod(A, "MatDenseGetColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3834: PetscFunctionReturn(PETSC_SUCCESS);
3835: }
3837: /*@
3838: MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVecWrite()`.
3840: Collective
3842: Input Parameters:
3843: + A - the `Mat` object
3844: . col - the column index
3845: - v - the `Vec` object (may be `NULL`)
3847: Level: intermediate
3849: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`
3850: @*/
3851: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A, PetscInt col, Vec *v)
3852: {
3853: PetscFunctionBegin;
3857: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3858: 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);
3859: PetscUseMethod(A, "MatDenseRestoreColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3860: PetscFunctionReturn(PETSC_SUCCESS);
3861: }
3863: /*@
3864: MatDenseGetSubMatrix - Gives access to a block of rows and columns of a dense matrix, represented as a `Mat`.
3866: Collective
3868: Input Parameters:
3869: + A - the `Mat` object
3870: . rbegin - the first global row index in the block (if `PETSC_DECIDE`, is 0)
3871: . rend - the global row index past the last one in the block (if `PETSC_DECIDE`, is `M`)
3872: . cbegin - the first global column index in the block (if `PETSC_DECIDE`, is 0)
3873: - cend - the global column index past the last one in the block (if `PETSC_DECIDE`, is `N`)
3875: Output Parameter:
3876: . v - the matrix
3878: Level: intermediate
3880: Notes:
3881: The matrix is owned by PETSc. Users need to call `MatDenseRestoreSubMatrix()` when the matrix is no longer needed.
3883: The output matrix is not redistributed by PETSc, so depending on the values of `rbegin` and `rend`, some processes may have no local rows.
3885: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreSubMatrix()`
3886: @*/
3887: PetscErrorCode MatDenseGetSubMatrix(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3888: {
3889: PetscFunctionBegin;
3896: PetscAssertPointer(v, 6);
3897: if (rbegin == PETSC_DECIDE) rbegin = 0;
3898: if (rend == PETSC_DECIDE) rend = A->rmap->N;
3899: if (cbegin == PETSC_DECIDE) cbegin = 0;
3900: if (cend == PETSC_DECIDE) cend = A->cmap->N;
3901: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3902: 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);
3903: 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);
3904: 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);
3905: 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);
3906: PetscUseMethod(A, "MatDenseGetSubMatrix_C", (Mat, PetscInt, PetscInt, PetscInt, PetscInt, Mat *), (A, rbegin, rend, cbegin, cend, v));
3907: PetscFunctionReturn(PETSC_SUCCESS);
3908: }
3910: /*@
3911: MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from `MatDenseGetSubMatrix()`.
3913: Collective
3915: Input Parameters:
3916: + A - the `Mat` object
3917: - v - the `Mat` object (may be `NULL`)
3919: Level: intermediate
3921: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseGetSubMatrix()`
3922: @*/
3923: PetscErrorCode MatDenseRestoreSubMatrix(Mat A, Mat *v)
3924: {
3925: PetscFunctionBegin;
3928: PetscAssertPointer(v, 2);
3929: PetscUseMethod(A, "MatDenseRestoreSubMatrix_C", (Mat, Mat *), (A, v));
3930: PetscFunctionReturn(PETSC_SUCCESS);
3931: }
3933: #include <petscblaslapack.h>
3934: #include <petsc/private/kernels/blockinvert.h>
3936: PetscErrorCode MatSeqDenseInvert(Mat A)
3937: {
3938: PetscInt m;
3939: const PetscReal shift = 0.0;
3940: PetscBool allowzeropivot, zeropivotdetected = PETSC_FALSE;
3941: PetscScalar *values;
3943: PetscFunctionBegin;
3945: PetscCall(MatDenseGetArray(A, &values));
3946: PetscCall(MatGetLocalSize(A, &m, NULL));
3947: allowzeropivot = PetscNot(A->erroriffailure);
3948: /* factor and invert each block */
3949: switch (m) {
3950: case 1:
3951: values[0] = (PetscScalar)1.0 / (values[0] + shift);
3952: break;
3953: case 2:
3954: PetscCall(PetscKernel_A_gets_inverse_A_2(values, shift, allowzeropivot, &zeropivotdetected));
3955: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3956: break;
3957: case 3:
3958: PetscCall(PetscKernel_A_gets_inverse_A_3(values, shift, allowzeropivot, &zeropivotdetected));
3959: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3960: break;
3961: case 4:
3962: PetscCall(PetscKernel_A_gets_inverse_A_4(values, shift, allowzeropivot, &zeropivotdetected));
3963: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3964: break;
3965: case 5: {
3966: PetscScalar work[25];
3967: PetscInt ipvt[5];
3969: PetscCall(PetscKernel_A_gets_inverse_A_5(values, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
3970: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3971: } break;
3972: case 6:
3973: PetscCall(PetscKernel_A_gets_inverse_A_6(values, shift, allowzeropivot, &zeropivotdetected));
3974: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3975: break;
3976: case 7:
3977: PetscCall(PetscKernel_A_gets_inverse_A_7(values, shift, allowzeropivot, &zeropivotdetected));
3978: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3979: break;
3980: default: {
3981: PetscInt *v_pivots, *IJ, j;
3982: PetscScalar *v_work;
3984: PetscCall(PetscMalloc3(m, &v_work, m, &v_pivots, m, &IJ));
3985: for (j = 0; j < m; j++) IJ[j] = j;
3986: PetscCall(PetscKernel_A_gets_inverse_A(m, values, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
3987: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3988: PetscCall(PetscFree3(v_work, v_pivots, IJ));
3989: }
3990: }
3991: PetscCall(MatDenseRestoreArray(A, &values));
3992: PetscFunctionReturn(PETSC_SUCCESS);
3993: }