Actual source code: dense.c
1: /*
2: Defines the basic matrix operations for sequential dense.
3: Portions of this code are under:
4: Copyright (c) 2022 Advanced Micro Devices, Inc. All rights reserved.
5: */
7: #include <../src/mat/impls/dense/seq/dense.h>
8: #include <../src/mat/impls/dense/mpi/mpidense.h>
9: #include <petscblaslapack.h>
10: #include <../src/mat/impls/aij/seq/aij.h>
12: PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian)
13: {
14: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
15: PetscInt j, k, n = A->rmap->n;
16: PetscScalar *v;
18: PetscFunctionBegin;
19: PetscCheck(A->rmap->n == A->cmap->n, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot symmetrize a rectangular matrix");
20: PetscCall(MatDenseGetArray(A, &v));
21: if (!hermitian) {
22: for (k = 0; k < n; k++) {
23: for (j = k; j < n; j++) v[j * mat->lda + k] = v[k * mat->lda + j];
24: }
25: } else {
26: for (k = 0; k < n; k++) {
27: for (j = k; j < n; j++) v[j * mat->lda + k] = PetscConj(v[k * mat->lda + j]);
28: }
29: }
30: PetscCall(MatDenseRestoreArray(A, &v));
31: PetscFunctionReturn(PETSC_SUCCESS);
32: }
34: PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A)
35: {
36: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
37: PetscBLASInt info, n;
39: PetscFunctionBegin;
40: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
41: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
42: if (A->factortype == MAT_FACTOR_LU) {
43: PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present");
44: if (!mat->fwork) {
45: mat->lfwork = n;
46: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
47: }
48: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
49: PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
50: PetscCall(PetscFPTrapPop());
51: PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
52: } else if (A->factortype == MAT_FACTOR_CHOLESKY) {
53: if (A->spd == PETSC_BOOL3_TRUE) {
54: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
55: PetscCallBLAS("LAPACKpotri", LAPACKpotri_("L", &n, mat->v, &mat->lda, &info));
56: PetscCall(PetscFPTrapPop());
57: PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_TRUE));
58: #if defined(PETSC_USE_COMPLEX)
59: } else if (A->hermitian == PETSC_BOOL3_TRUE) {
60: PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present");
61: PetscCheck(mat->fwork, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Fwork not present");
62: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
63: PetscCallBLAS("LAPACKhetri", LAPACKhetri_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &info));
64: PetscCall(PetscFPTrapPop());
65: PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_TRUE));
66: #endif
67: } else { /* symmetric case */
68: PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present");
69: PetscCheck(mat->fwork, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Fwork not present");
70: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
71: PetscCallBLAS("LAPACKsytri", LAPACKsytri_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &info));
72: PetscCall(PetscFPTrapPop());
73: PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_FALSE));
74: }
75: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad Inversion: zero pivot in row %" PetscInt_FMT, (PetscInt)info - 1);
76: PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
77: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix must be factored to solve");
79: A->ops->solve = NULL;
80: A->ops->matsolve = NULL;
81: A->ops->solvetranspose = NULL;
82: A->ops->matsolvetranspose = NULL;
83: A->ops->solveadd = NULL;
84: A->ops->solvetransposeadd = NULL;
85: A->factortype = MAT_FACTOR_NONE;
86: PetscCall(PetscFree(A->solvertype));
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
90: static PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
91: {
92: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
93: PetscInt m = l->lda, n = A->cmap->n, r = A->rmap->n, i, j;
94: PetscScalar *slot, *bb, *v;
95: const PetscScalar *xx;
97: PetscFunctionBegin;
98: if (PetscDefined(USE_DEBUG)) {
99: for (i = 0; i < N; i++) {
100: PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
101: PetscCheck(rows[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " requested to be zeroed greater than or equal number of rows %" PetscInt_FMT, rows[i], A->rmap->n);
102: PetscCheck(rows[i] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Col %" PetscInt_FMT " requested to be zeroed greater than or equal number of cols %" PetscInt_FMT, rows[i], A->cmap->n);
103: }
104: }
105: if (!N) PetscFunctionReturn(PETSC_SUCCESS);
107: /* fix right-hand side if needed */
108: if (x && b) {
109: Vec xt;
111: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
112: PetscCall(VecDuplicate(x, &xt));
113: PetscCall(VecCopy(x, xt));
114: PetscCall(VecScale(xt, -1.0));
115: PetscCall(MatMultAdd(A, xt, b, b));
116: PetscCall(VecDestroy(&xt));
117: PetscCall(VecGetArrayRead(x, &xx));
118: PetscCall(VecGetArray(b, &bb));
119: for (i = 0; i < N; i++) bb[rows[i]] = diag * xx[rows[i]];
120: PetscCall(VecRestoreArrayRead(x, &xx));
121: PetscCall(VecRestoreArray(b, &bb));
122: }
124: PetscCall(MatDenseGetArray(A, &v));
125: for (i = 0; i < N; i++) {
126: slot = v + rows[i] * m;
127: PetscCall(PetscArrayzero(slot, r));
128: }
129: for (i = 0; i < N; i++) {
130: slot = v + rows[i];
131: for (j = 0; j < n; j++) {
132: *slot = 0.0;
133: slot += m;
134: }
135: }
136: if (diag != 0.0) {
137: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
138: for (i = 0; i < N; i++) {
139: slot = v + (m + 1) * rows[i];
140: *slot = diag;
141: }
142: }
143: PetscCall(MatDenseRestoreArray(A, &v));
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
148: {
149: Mat B = NULL;
150: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
151: Mat_SeqDense *b;
152: PetscInt *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->N, i;
153: const MatScalar *av;
154: PetscBool isseqdense;
156: PetscFunctionBegin;
157: if (reuse == MAT_REUSE_MATRIX) {
158: PetscCall(PetscObjectTypeCompare((PetscObject)*newmat, MATSEQDENSE, &isseqdense));
159: PetscCheck(isseqdense, PetscObjectComm((PetscObject)*newmat), PETSC_ERR_USER, "Cannot reuse matrix of type %s", ((PetscObject)*newmat)->type_name);
160: }
161: if (reuse != MAT_REUSE_MATRIX) {
162: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
163: PetscCall(MatSetSizes(B, m, n, m, n));
164: PetscCall(MatSetType(B, MATSEQDENSE));
165: PetscCall(MatSeqDenseSetPreallocation(B, NULL));
166: b = (Mat_SeqDense *)B->data;
167: } else {
168: b = (Mat_SeqDense *)((*newmat)->data);
169: for (i = 0; i < n; i++) PetscCall(PetscArrayzero(b->v + i * b->lda, m));
170: }
171: PetscCall(MatSeqAIJGetArrayRead(A, &av));
172: for (i = 0; i < m; i++) {
173: PetscInt j;
174: for (j = 0; j < ai[1] - ai[0]; j++) {
175: b->v[*aj * b->lda + i] = *av;
176: aj++;
177: av++;
178: }
179: ai++;
180: }
181: PetscCall(MatSeqAIJRestoreArrayRead(A, &av));
183: if (reuse == MAT_INPLACE_MATRIX) {
184: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
185: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
186: PetscCall(MatHeaderReplace(A, &B));
187: } else {
188: if (B) *newmat = B;
189: PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
190: PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
191: }
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
196: {
197: Mat B = NULL;
198: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
199: PetscInt i, j;
200: PetscInt *rows, *nnz;
201: MatScalar *aa = a->v, *vals;
203: PetscFunctionBegin;
204: PetscCall(PetscCalloc3(A->rmap->n, &rows, A->rmap->n, &nnz, A->rmap->n, &vals));
205: if (reuse != MAT_REUSE_MATRIX) {
206: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
207: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
208: PetscCall(MatSetType(B, MATSEQAIJ));
209: for (j = 0; j < A->cmap->n; j++) {
210: for (i = 0; i < A->rmap->n; i++)
211: if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) ++nnz[i];
212: aa += a->lda;
213: }
214: PetscCall(MatSeqAIJSetPreallocation(B, PETSC_DETERMINE, nnz));
215: } else B = *newmat;
216: aa = a->v;
217: for (j = 0; j < A->cmap->n; j++) {
218: PetscInt numRows = 0;
219: for (i = 0; i < A->rmap->n; i++)
220: if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) {
221: rows[numRows] = i;
222: vals[numRows++] = aa[i];
223: }
224: PetscCall(MatSetValues(B, numRows, rows, 1, &j, vals, INSERT_VALUES));
225: aa += a->lda;
226: }
227: PetscCall(PetscFree3(rows, nnz, vals));
228: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
229: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
231: if (reuse == MAT_INPLACE_MATRIX) {
232: PetscCall(MatHeaderReplace(A, &B));
233: } else if (reuse != MAT_REUSE_MATRIX) *newmat = B;
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: PetscErrorCode MatAXPY_SeqDense(Mat Y, PetscScalar alpha, Mat X, MatStructure str)
238: {
239: Mat_SeqDense *x = (Mat_SeqDense *)X->data, *y = (Mat_SeqDense *)Y->data;
240: const PetscScalar *xv;
241: PetscScalar *yv;
242: PetscBLASInt N, m, ldax = 0, lday = 0, one = 1;
244: PetscFunctionBegin;
245: PetscCall(MatDenseGetArrayRead(X, &xv));
246: PetscCall(MatDenseGetArray(Y, &yv));
247: PetscCall(PetscBLASIntCast(X->rmap->n * X->cmap->n, &N));
248: PetscCall(PetscBLASIntCast(X->rmap->n, &m));
249: PetscCall(PetscBLASIntCast(x->lda, &ldax));
250: PetscCall(PetscBLASIntCast(y->lda, &lday));
251: if (ldax > m || lday > m) {
252: 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 = (PetscInt)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: }
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: PetscErrorCode MatDuplicate_SeqDense(Mat A, MatDuplicateOption cpvalues, Mat *newmat)
392: {
393: PetscFunctionBegin;
394: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), newmat));
395: PetscCall(MatSetSizes(*newmat, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
396: PetscCall(MatSetType(*newmat, ((PetscObject)A)->type_name));
397: PetscCall(MatDuplicateNoCreate_SeqDense(*newmat, A, cpvalues));
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
401: static PetscErrorCode MatSolve_SeqDense_Internal_LU(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
402: {
403: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
404: PetscBLASInt info;
406: PetscFunctionBegin;
407: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
408: PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_(T ? "T" : "N", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
409: PetscCall(PetscFPTrapPop());
410: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "GETRS - Bad solve %d", (int)info);
411: PetscCall(PetscLogFlops(nrhs * (2.0 * m * m - m)));
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: static PetscErrorCode MatConjugate_SeqDense(Mat);
417: static PetscErrorCode MatSolve_SeqDense_Internal_Cholesky(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
418: {
419: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
420: PetscBLASInt info;
422: PetscFunctionBegin;
423: if (A->spd == PETSC_BOOL3_TRUE) {
424: if (PetscDefined(USE_COMPLEX) && T) PetscCall(MatConjugate_SeqDense(A));
425: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
426: PetscCallBLAS("LAPACKpotrs", LAPACKpotrs_("L", &m, &nrhs, mat->v, &mat->lda, x, &m, &info));
427: PetscCall(PetscFPTrapPop());
428: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "POTRS Bad solve %d", (int)info);
429: if (PetscDefined(USE_COMPLEX) && T) PetscCall(MatConjugate_SeqDense(A));
430: #if defined(PETSC_USE_COMPLEX)
431: } else if (A->hermitian == PETSC_BOOL3_TRUE) {
432: if (T) PetscCall(MatConjugate_SeqDense(A));
433: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
434: PetscCallBLAS("LAPACKhetrs", LAPACKhetrs_("L", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
435: PetscCall(PetscFPTrapPop());
436: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "HETRS Bad solve %d", (int)info);
437: if (T) PetscCall(MatConjugate_SeqDense(A));
438: #endif
439: } else { /* symmetric case */
440: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
441: PetscCallBLAS("LAPACKsytrs", LAPACKsytrs_("L", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
442: PetscCall(PetscFPTrapPop());
443: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "SYTRS Bad solve %d", (int)info);
444: }
445: PetscCall(PetscLogFlops(nrhs * (2.0 * m * m - m)));
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: static PetscErrorCode MatSolve_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
450: {
451: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
452: PetscBLASInt info;
453: char trans;
455: PetscFunctionBegin;
456: if (PetscDefined(USE_COMPLEX)) {
457: trans = 'C';
458: } else {
459: trans = 'T';
460: }
461: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
462: { /* lwork depends on the number of right-hand sides */
463: PetscBLASInt nlfwork, lfwork = -1;
464: PetscScalar fwork;
466: PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", &trans, &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, &fwork, &lfwork, &info));
467: nlfwork = (PetscBLASInt)PetscRealPart(fwork);
468: if (nlfwork > mat->lfwork) {
469: mat->lfwork = nlfwork;
470: PetscCall(PetscFree(mat->fwork));
471: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
472: }
473: }
474: PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", &trans, &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, mat->fwork, &mat->lfwork, &info));
475: PetscCall(PetscFPTrapPop());
476: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "ORMQR - Bad orthogonal transform %d", (int)info);
477: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
478: PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "N", "N", &mat->rank, &nrhs, mat->v, &mat->lda, x, &ldx, &info));
479: PetscCall(PetscFPTrapPop());
480: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "TRTRS - Bad triangular solve %d", (int)info);
481: for (PetscInt j = 0; j < nrhs; j++) {
482: for (PetscInt i = mat->rank; i < k; i++) x[j * ldx + i] = 0.;
483: }
484: PetscCall(PetscLogFlops(nrhs * (4.0 * m * mat->rank - PetscSqr(mat->rank))));
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: static PetscErrorCode MatSolveTranspose_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
489: {
490: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
491: PetscBLASInt info;
493: PetscFunctionBegin;
494: if (A->rmap->n == A->cmap->n && mat->rank == A->rmap->n) {
495: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
496: PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "T", "N", &m, &nrhs, mat->v, &mat->lda, x, &ldx, &info));
497: PetscCall(PetscFPTrapPop());
498: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "TRTRS - Bad triangular solve %d", (int)info);
499: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate_SeqDense(A));
500: { /* lwork depends on the number of right-hand sides */
501: PetscBLASInt nlfwork, lfwork = -1;
502: PetscScalar fwork;
504: PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", "N", &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, &fwork, &lfwork, &info));
505: nlfwork = (PetscBLASInt)PetscRealPart(fwork);
506: if (nlfwork > mat->lfwork) {
507: mat->lfwork = nlfwork;
508: PetscCall(PetscFree(mat->fwork));
509: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
510: }
511: }
512: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
513: PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", "N", &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, mat->fwork, &mat->lfwork, &info));
514: PetscCall(PetscFPTrapPop());
515: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "ORMQR - Bad orthogonal transform %d", (int)info);
516: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate_SeqDense(A));
517: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "QR factored matrix cannot be used for transpose solve");
518: PetscCall(PetscLogFlops(nrhs * (4.0 * m * mat->rank - PetscSqr(mat->rank))));
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: static PetscErrorCode MatSolve_SeqDense_SetUp(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
523: {
524: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
525: PetscScalar *y;
526: PetscBLASInt m = 0, k = 0;
528: PetscFunctionBegin;
529: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
530: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
531: if (k < m) {
532: PetscCall(VecCopy(xx, mat->qrrhs));
533: PetscCall(VecGetArray(mat->qrrhs, &y));
534: } else {
535: PetscCall(VecCopy(xx, yy));
536: PetscCall(VecGetArray(yy, &y));
537: }
538: *_y = y;
539: *_k = k;
540: *_m = m;
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: static PetscErrorCode MatSolve_SeqDense_TearDown(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
545: {
546: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
547: PetscScalar *y = NULL;
548: PetscBLASInt m, k;
550: PetscFunctionBegin;
551: y = *_y;
552: *_y = NULL;
553: k = *_k;
554: m = *_m;
555: if (k < m) {
556: PetscScalar *yv;
557: PetscCall(VecGetArray(yy, &yv));
558: PetscCall(PetscArraycpy(yv, y, k));
559: PetscCall(VecRestoreArray(yy, &yv));
560: PetscCall(VecRestoreArray(mat->qrrhs, &y));
561: } else {
562: PetscCall(VecRestoreArray(yy, &y));
563: }
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: static PetscErrorCode MatSolve_SeqDense_LU(Mat A, Vec xx, Vec yy)
568: {
569: PetscScalar *y = NULL;
570: PetscBLASInt m = 0, k = 0;
572: PetscFunctionBegin;
573: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
574: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_FALSE));
575: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: static PetscErrorCode MatSolveTranspose_SeqDense_LU(Mat A, Vec xx, Vec yy)
580: {
581: PetscScalar *y = NULL;
582: PetscBLASInt m = 0, k = 0;
584: PetscFunctionBegin;
585: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
586: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_TRUE));
587: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
591: static PetscErrorCode MatSolve_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
592: {
593: PetscScalar *y = NULL;
594: PetscBLASInt m = 0, k = 0;
596: PetscFunctionBegin;
597: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
598: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_FALSE));
599: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
600: PetscFunctionReturn(PETSC_SUCCESS);
601: }
603: static PetscErrorCode MatSolveTranspose_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
604: {
605: PetscScalar *y = NULL;
606: PetscBLASInt m = 0, k = 0;
608: PetscFunctionBegin;
609: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
610: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_TRUE));
611: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: static PetscErrorCode MatSolve_SeqDense_QR(Mat A, Vec xx, Vec yy)
616: {
617: PetscScalar *y = NULL;
618: PetscBLASInt m = 0, k = 0;
620: PetscFunctionBegin;
621: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
622: PetscCall(MatSolve_SeqDense_Internal_QR(A, y, PetscMax(m, k), m, 1, k));
623: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
624: PetscFunctionReturn(PETSC_SUCCESS);
625: }
627: static PetscErrorCode MatSolveTranspose_SeqDense_QR(Mat A, Vec xx, Vec yy)
628: {
629: PetscScalar *y = NULL;
630: PetscBLASInt m = 0, k = 0;
632: PetscFunctionBegin;
633: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
634: PetscCall(MatSolveTranspose_SeqDense_Internal_QR(A, y, PetscMax(m, k), m, 1, k));
635: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
636: PetscFunctionReturn(PETSC_SUCCESS);
637: }
639: static PetscErrorCode MatMatSolve_SeqDense_SetUp(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
640: {
641: const PetscScalar *b;
642: PetscScalar *y;
643: PetscInt n, _ldb, _ldx;
644: PetscBLASInt nrhs = 0, m = 0, k = 0, ldb = 0, ldx = 0, ldy = 0;
646: PetscFunctionBegin;
647: *_ldy = 0;
648: *_m = 0;
649: *_nrhs = 0;
650: *_k = 0;
651: *_y = NULL;
652: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
653: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
654: PetscCall(MatGetSize(B, NULL, &n));
655: PetscCall(PetscBLASIntCast(n, &nrhs));
656: PetscCall(MatDenseGetLDA(B, &_ldb));
657: PetscCall(PetscBLASIntCast(_ldb, &ldb));
658: PetscCall(MatDenseGetLDA(X, &_ldx));
659: PetscCall(PetscBLASIntCast(_ldx, &ldx));
660: if (ldx < m) {
661: PetscCall(MatDenseGetArrayRead(B, &b));
662: PetscCall(PetscMalloc1(nrhs * m, &y));
663: if (ldb == m) {
664: PetscCall(PetscArraycpy(y, b, ldb * nrhs));
665: } else {
666: for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&y[j * m], &b[j * ldb], m));
667: }
668: ldy = m;
669: PetscCall(MatDenseRestoreArrayRead(B, &b));
670: } else {
671: if (ldb == ldx) {
672: PetscCall(MatCopy(B, X, SAME_NONZERO_PATTERN));
673: PetscCall(MatDenseGetArray(X, &y));
674: } else {
675: PetscCall(MatDenseGetArray(X, &y));
676: PetscCall(MatDenseGetArrayRead(B, &b));
677: for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&y[j * ldx], &b[j * ldb], m));
678: PetscCall(MatDenseRestoreArrayRead(B, &b));
679: }
680: ldy = ldx;
681: }
682: *_y = y;
683: *_ldy = ldy;
684: *_k = k;
685: *_m = m;
686: *_nrhs = nrhs;
687: PetscFunctionReturn(PETSC_SUCCESS);
688: }
690: static PetscErrorCode MatMatSolve_SeqDense_TearDown(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
691: {
692: PetscScalar *y;
693: PetscInt _ldx;
694: PetscBLASInt k, ldy, nrhs, ldx = 0;
696: PetscFunctionBegin;
697: y = *_y;
698: *_y = NULL;
699: k = *_k;
700: ldy = *_ldy;
701: nrhs = *_nrhs;
702: PetscCall(MatDenseGetLDA(X, &_ldx));
703: PetscCall(PetscBLASIntCast(_ldx, &ldx));
704: if (ldx != ldy) {
705: PetscScalar *xv;
706: PetscCall(MatDenseGetArray(X, &xv));
707: for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&xv[j * ldx], &y[j * ldy], k));
708: PetscCall(MatDenseRestoreArray(X, &xv));
709: PetscCall(PetscFree(y));
710: } else {
711: PetscCall(MatDenseRestoreArray(X, &y));
712: }
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }
716: static PetscErrorCode MatMatSolve_SeqDense_LU(Mat A, Mat B, Mat X)
717: {
718: PetscScalar *y;
719: PetscBLASInt m, k, ldy, nrhs;
721: PetscFunctionBegin;
722: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
723: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_FALSE));
724: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
725: PetscFunctionReturn(PETSC_SUCCESS);
726: }
728: static PetscErrorCode MatMatSolveTranspose_SeqDense_LU(Mat A, Mat B, Mat X)
729: {
730: PetscScalar *y;
731: PetscBLASInt m, k, ldy, nrhs;
733: PetscFunctionBegin;
734: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
735: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_TRUE));
736: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
737: PetscFunctionReturn(PETSC_SUCCESS);
738: }
740: static PetscErrorCode MatMatSolve_SeqDense_Cholesky(Mat A, Mat B, Mat X)
741: {
742: PetscScalar *y;
743: PetscBLASInt m, k, ldy, nrhs;
745: PetscFunctionBegin;
746: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
747: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_FALSE));
748: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
752: static PetscErrorCode MatMatSolveTranspose_SeqDense_Cholesky(Mat A, Mat B, Mat X)
753: {
754: PetscScalar *y;
755: PetscBLASInt m, k, ldy, nrhs;
757: PetscFunctionBegin;
758: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
759: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_TRUE));
760: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
761: PetscFunctionReturn(PETSC_SUCCESS);
762: }
764: static PetscErrorCode MatMatSolve_SeqDense_QR(Mat A, Mat B, Mat X)
765: {
766: PetscScalar *y;
767: PetscBLASInt m, k, ldy, nrhs;
769: PetscFunctionBegin;
770: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
771: PetscCall(MatSolve_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k));
772: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
773: PetscFunctionReturn(PETSC_SUCCESS);
774: }
776: static PetscErrorCode MatMatSolveTranspose_SeqDense_QR(Mat A, Mat B, Mat X)
777: {
778: PetscScalar *y;
779: PetscBLASInt m, k, ldy, nrhs;
781: PetscFunctionBegin;
782: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
783: PetscCall(MatSolveTranspose_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k));
784: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
785: PetscFunctionReturn(PETSC_SUCCESS);
786: }
788: /* COMMENT: I have chosen to hide row permutation in the pivots,
789: rather than put it in the Mat->row slot.*/
790: PetscErrorCode MatLUFactor_SeqDense(Mat A, IS row, IS col, const MatFactorInfo *minfo)
791: {
792: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
793: PetscBLASInt n, m, info;
795: PetscFunctionBegin;
796: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
797: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
798: if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); }
799: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
800: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
801: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&m, &n, mat->v, &mat->lda, mat->pivots, &info));
802: PetscCall(PetscFPTrapPop());
804: PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to LU factorization %d", (int)info);
805: PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Bad LU factorization %d", (int)info);
807: A->ops->solve = MatSolve_SeqDense_LU;
808: A->ops->matsolve = MatMatSolve_SeqDense_LU;
809: A->ops->solvetranspose = MatSolveTranspose_SeqDense_LU;
810: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_LU;
811: A->factortype = MAT_FACTOR_LU;
813: PetscCall(PetscFree(A->solvertype));
814: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
816: PetscCall(PetscLogFlops((2.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3));
817: PetscFunctionReturn(PETSC_SUCCESS);
818: }
820: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info_dummy)
821: {
822: MatFactorInfo info;
824: PetscFunctionBegin;
825: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
826: PetscUseTypeMethod(fact, lufactor, NULL, NULL, &info);
827: PetscFunctionReturn(PETSC_SUCCESS);
828: }
830: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, IS col, const MatFactorInfo *info)
831: {
832: PetscFunctionBegin;
833: fact->preallocated = PETSC_TRUE;
834: fact->assembled = PETSC_TRUE;
835: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
836: PetscFunctionReturn(PETSC_SUCCESS);
837: }
839: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
840: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A, IS perm, const MatFactorInfo *factinfo)
841: {
842: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
843: PetscBLASInt info, n;
845: PetscFunctionBegin;
846: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
847: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
848: if (A->spd == PETSC_BOOL3_TRUE) {
849: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
850: PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("L", &n, mat->v, &mat->lda, &info));
851: PetscCall(PetscFPTrapPop());
852: #if defined(PETSC_USE_COMPLEX)
853: } else if (A->hermitian == PETSC_BOOL3_TRUE) {
854: if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); }
855: if (!mat->fwork) {
856: PetscScalar dummy;
858: mat->lfwork = -1;
859: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
860: PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info));
861: PetscCall(PetscFPTrapPop());
862: mat->lfwork = (PetscInt)PetscRealPart(dummy);
863: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
864: }
865: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
866: PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
867: PetscCall(PetscFPTrapPop());
868: #endif
869: } else { /* symmetric case */
870: if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); }
871: if (!mat->fwork) {
872: PetscScalar dummy;
874: mat->lfwork = -1;
875: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
876: PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info));
877: PetscCall(PetscFPTrapPop());
878: mat->lfwork = (PetscInt)PetscRealPart(dummy);
879: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
880: }
881: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
882: PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
883: PetscCall(PetscFPTrapPop());
884: }
885: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad factorization: zero pivot in row %" PetscInt_FMT, (PetscInt)info - 1);
887: A->ops->solve = MatSolve_SeqDense_Cholesky;
888: A->ops->matsolve = MatMatSolve_SeqDense_Cholesky;
889: A->ops->solvetranspose = MatSolveTranspose_SeqDense_Cholesky;
890: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky;
891: A->factortype = MAT_FACTOR_CHOLESKY;
893: PetscCall(PetscFree(A->solvertype));
894: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
896: PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info_dummy)
901: {
902: MatFactorInfo info;
904: PetscFunctionBegin;
905: info.fill = 1.0;
907: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
908: PetscUseTypeMethod(fact, choleskyfactor, NULL, &info);
909: PetscFunctionReturn(PETSC_SUCCESS);
910: }
912: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info)
913: {
914: PetscFunctionBegin;
915: fact->assembled = PETSC_TRUE;
916: fact->preallocated = PETSC_TRUE;
917: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
918: PetscFunctionReturn(PETSC_SUCCESS);
919: }
921: PetscErrorCode MatQRFactor_SeqDense(Mat A, IS col, const MatFactorInfo *minfo)
922: {
923: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
924: PetscBLASInt n, m, info, min, max;
926: PetscFunctionBegin;
927: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
928: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
929: max = PetscMax(m, n);
930: min = PetscMin(m, n);
931: if (!mat->tau) { PetscCall(PetscMalloc1(min, &mat->tau)); }
932: if (!mat->pivots) { PetscCall(PetscMalloc1(n, &mat->pivots)); }
933: if (!mat->qrrhs) PetscCall(MatCreateVecs(A, NULL, &mat->qrrhs));
934: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
935: if (!mat->fwork) {
936: PetscScalar dummy;
938: mat->lfwork = -1;
939: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
940: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, &dummy, &mat->lfwork, &info));
941: PetscCall(PetscFPTrapPop());
942: mat->lfwork = (PetscInt)PetscRealPart(dummy);
943: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
944: }
945: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
946: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, mat->fwork, &mat->lfwork, &info));
947: PetscCall(PetscFPTrapPop());
948: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to QR factorization %d", (int)info);
949: // 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
950: mat->rank = min;
952: A->ops->solve = MatSolve_SeqDense_QR;
953: A->ops->matsolve = MatMatSolve_SeqDense_QR;
954: A->factortype = MAT_FACTOR_QR;
955: if (m == n) {
956: A->ops->solvetranspose = MatSolveTranspose_SeqDense_QR;
957: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
958: }
960: PetscCall(PetscFree(A->solvertype));
961: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
963: PetscCall(PetscLogFlops(2.0 * min * min * (max - min / 3.0)));
964: PetscFunctionReturn(PETSC_SUCCESS);
965: }
967: static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info_dummy)
968: {
969: MatFactorInfo info;
971: PetscFunctionBegin;
972: info.fill = 1.0;
974: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
975: PetscUseMethod(fact, "MatQRFactor_C", (Mat, IS, const MatFactorInfo *), (fact, NULL, &info));
976: PetscFunctionReturn(PETSC_SUCCESS);
977: }
979: PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info)
980: {
981: PetscFunctionBegin;
982: fact->assembled = PETSC_TRUE;
983: fact->preallocated = PETSC_TRUE;
984: PetscCall(PetscObjectComposeFunction((PetscObject)fact, "MatQRFactorNumeric_C", MatQRFactorNumeric_SeqDense));
985: PetscFunctionReturn(PETSC_SUCCESS);
986: }
988: /* uses LAPACK */
989: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A, MatFactorType ftype, Mat *fact)
990: {
991: PetscFunctionBegin;
992: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), fact));
993: PetscCall(MatSetSizes(*fact, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
994: PetscCall(MatSetType(*fact, MATDENSE));
995: (*fact)->trivialsymbolic = PETSC_TRUE;
996: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
997: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
998: (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
999: } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
1000: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
1001: } else if (ftype == MAT_FACTOR_QR) {
1002: PetscCall(PetscObjectComposeFunction((PetscObject)*fact, "MatQRFactorSymbolic_C", MatQRFactorSymbolic_SeqDense));
1003: }
1004: (*fact)->factortype = ftype;
1006: PetscCall(PetscFree((*fact)->solvertype));
1007: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*fact)->solvertype));
1008: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_LU]));
1009: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ILU]));
1010: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_CHOLESKY]));
1011: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ICC]));
1012: PetscFunctionReturn(PETSC_SUCCESS);
1013: }
1015: static PetscErrorCode MatSOR_SeqDense(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal shift, PetscInt its, PetscInt lits, Vec xx)
1016: {
1017: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1018: PetscScalar *x, *v = mat->v, zero = 0.0, xt;
1019: const PetscScalar *b;
1020: PetscInt m = A->rmap->n, i;
1021: PetscBLASInt o = 1, bm = 0;
1023: PetscFunctionBegin;
1024: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1025: PetscCheck(A->offloadmask != PETSC_OFFLOAD_GPU, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented");
1026: #endif
1027: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
1028: PetscCall(PetscBLASIntCast(m, &bm));
1029: if (flag & SOR_ZERO_INITIAL_GUESS) {
1030: /* this is a hack fix, should have another version without the second BLASdotu */
1031: PetscCall(VecSet(xx, zero));
1032: }
1033: PetscCall(VecGetArray(xx, &x));
1034: PetscCall(VecGetArrayRead(bb, &b));
1035: its = its * lits;
1036: 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);
1037: while (its--) {
1038: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1039: for (i = 0; i < m; i++) {
1040: PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o));
1041: x[i] = (1. - omega) * x[i] + omega * (xt + v[i + i * m] * x[i]) / (v[i + i * m] + shift);
1042: }
1043: }
1044: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1045: for (i = m - 1; i >= 0; i--) {
1046: PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o));
1047: x[i] = (1. - omega) * x[i] + omega * (xt + v[i + i * m] * x[i]) / (v[i + i * m] + shift);
1048: }
1049: }
1050: }
1051: PetscCall(VecRestoreArrayRead(bb, &b));
1052: PetscCall(VecRestoreArray(xx, &x));
1053: PetscFunctionReturn(PETSC_SUCCESS);
1054: }
1056: static PetscErrorCode MatMultColumnRangeKernel_SeqDense(Mat A, Vec xx, Vec yy, PetscInt c_start, PetscInt c_end, PetscBool trans, PetscBool herm)
1057: {
1058: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1059: PetscScalar *y, _DOne = 1.0, _DZero = 0.0;
1060: PetscBLASInt m, n, _One = 1;
1061: const PetscScalar *v = mat->v, *x;
1063: PetscFunctionBegin;
1064: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1065: PetscCall(PetscBLASIntCast(c_end - c_start, &n));
1066: PetscCall(VecGetArrayRead(xx, &x));
1067: PetscCall(VecGetArrayWrite(yy, &y));
1068: if (!m || !n) {
1069: PetscBLASInt i;
1070: if (trans)
1071: for (i = 0; i < n; i++) y[i] = 0.0;
1072: else
1073: for (i = 0; i < m; i++) y[i] = 0.0;
1074: } else {
1075: if (trans) {
1076: if (herm) PetscCallBLAS("BLASgemv", BLASgemv_("C", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DZero, y + c_start, &_One));
1077: else PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DZero, y + c_start, &_One));
1078: } else {
1079: PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x + c_start, &_One, &_DZero, y, &_One));
1080: }
1081: PetscCall(PetscLogFlops(2.0 * m * n - n));
1082: }
1083: PetscCall(VecRestoreArrayRead(xx, &x));
1084: PetscCall(VecRestoreArrayWrite(yy, &y));
1085: PetscFunctionReturn(PETSC_SUCCESS);
1086: }
1088: PetscErrorCode MatMultHermitianTransposeColumnRange_SeqDense(Mat A, Vec xx, Vec yy, PetscInt c_start, PetscInt c_end)
1089: {
1090: PetscFunctionBegin;
1091: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, c_start, c_end, PETSC_TRUE, PETSC_TRUE));
1092: PetscFunctionReturn(PETSC_SUCCESS);
1093: }
1095: PetscErrorCode MatMult_SeqDense(Mat A, Vec xx, Vec yy)
1096: {
1097: PetscFunctionBegin;
1098: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_FALSE, PETSC_FALSE));
1099: PetscFunctionReturn(PETSC_SUCCESS);
1100: }
1102: PetscErrorCode MatMultTranspose_SeqDense(Mat A, Vec xx, Vec yy)
1103: {
1104: PetscFunctionBegin;
1105: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_FALSE));
1106: PetscFunctionReturn(PETSC_SUCCESS);
1107: }
1109: PetscErrorCode MatMultHermitianTranspose_SeqDense(Mat A, Vec xx, Vec yy)
1110: {
1111: PetscFunctionBegin;
1112: PetscCall(MatMultColumnRangeKernel_SeqDense(A, xx, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_TRUE));
1113: PetscFunctionReturn(PETSC_SUCCESS);
1114: }
1116: static PetscErrorCode MatMultAddColumnRangeKernel_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end, PetscBool trans, PetscBool herm)
1117: {
1118: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1119: const PetscScalar *v = mat->v, *x;
1120: PetscScalar *y, _DOne = 1.0;
1121: PetscBLASInt m, n, _One = 1;
1123: PetscFunctionBegin;
1124: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1125: PetscCall(PetscBLASIntCast(c_end - c_start, &n));
1126: PetscCall(VecCopy(zz, yy));
1127: if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS);
1128: PetscCall(VecGetArray(yy, &y));
1129: PetscCall(VecGetArrayRead(xx, &x));
1130: if (trans) {
1131: if (herm) PetscCallBLAS("BLASgemv", BLASgemv_("C", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DOne, y + c_start, &_One));
1132: else PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x, &_One, &_DOne, y + c_start, &_One));
1133: } else {
1134: PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v + c_start * mat->lda, &mat->lda, x + c_start, &_One, &_DOne, y, &_One));
1135: }
1136: PetscCall(VecRestoreArrayRead(xx, &x));
1137: PetscCall(VecRestoreArray(yy, &y));
1138: PetscCall(PetscLogFlops(2.0 * m * n));
1139: PetscFunctionReturn(PETSC_SUCCESS);
1140: }
1142: PetscErrorCode MatMultAddColumnRange_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end)
1143: {
1144: PetscFunctionBegin;
1145: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, c_start, c_end, PETSC_FALSE, PETSC_FALSE));
1146: PetscFunctionReturn(PETSC_SUCCESS);
1147: }
1149: PetscErrorCode MatMultHermitianTransposeAddColumnRange_SeqDense(Mat A, Vec xx, Vec zz, Vec yy, PetscInt c_start, PetscInt c_end)
1150: {
1151: PetscFunctionBegin;
1152: PetscMPIInt rank;
1153: PetscCallMPI(MPI_Comm_rank(MPI_COMM_WORLD, &rank));
1154: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, c_start, c_end, PETSC_TRUE, PETSC_TRUE));
1155: PetscFunctionReturn(PETSC_SUCCESS);
1156: }
1158: PetscErrorCode MatMultAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1159: {
1160: PetscFunctionBegin;
1161: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_FALSE, PETSC_FALSE));
1162: PetscFunctionReturn(PETSC_SUCCESS);
1163: }
1165: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1166: {
1167: PetscFunctionBegin;
1168: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_FALSE));
1169: PetscFunctionReturn(PETSC_SUCCESS);
1170: }
1172: PetscErrorCode MatMultHermitianTransposeAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1173: {
1174: PetscFunctionBegin;
1175: PetscCall(MatMultAddColumnRangeKernel_SeqDense(A, xx, zz, yy, 0, A->cmap->n, PETSC_TRUE, PETSC_TRUE));
1176: PetscFunctionReturn(PETSC_SUCCESS);
1177: }
1179: static PetscErrorCode MatGetRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1180: {
1181: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1182: PetscInt i;
1184: PetscFunctionBegin;
1185: if (ncols) *ncols = A->cmap->n;
1186: if (cols) {
1187: PetscCall(PetscMalloc1(A->cmap->n, cols));
1188: for (i = 0; i < A->cmap->n; i++) (*cols)[i] = i;
1189: }
1190: if (vals) {
1191: const PetscScalar *v;
1193: PetscCall(MatDenseGetArrayRead(A, &v));
1194: PetscCall(PetscMalloc1(A->cmap->n, vals));
1195: v += row;
1196: for (i = 0; i < A->cmap->n; i++) {
1197: (*vals)[i] = *v;
1198: v += mat->lda;
1199: }
1200: PetscCall(MatDenseRestoreArrayRead(A, &v));
1201: }
1202: PetscFunctionReturn(PETSC_SUCCESS);
1203: }
1205: static PetscErrorCode MatRestoreRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1206: {
1207: PetscFunctionBegin;
1208: if (cols) PetscCall(PetscFree(*cols));
1209: if (vals) PetscCall(PetscFree(*vals));
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: static PetscErrorCode MatSetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], const PetscScalar v[], InsertMode addv)
1214: {
1215: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1216: PetscScalar *av;
1217: PetscInt i, j, idx = 0;
1218: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1219: PetscOffloadMask oldf;
1220: #endif
1222: PetscFunctionBegin;
1223: PetscCall(MatDenseGetArray(A, &av));
1224: if (!mat->roworiented) {
1225: if (addv == INSERT_VALUES) {
1226: for (j = 0; j < n; j++) {
1227: if (indexn[j] < 0) {
1228: idx += m;
1229: continue;
1230: }
1231: PetscCheck(indexn[j] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, indexn[j], A->cmap->n - 1);
1232: for (i = 0; i < m; i++) {
1233: if (indexm[i] < 0) {
1234: idx++;
1235: continue;
1236: }
1237: PetscCheck(indexm[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, indexm[i], A->rmap->n - 1);
1238: av[indexn[j] * mat->lda + indexm[i]] = v[idx++];
1239: }
1240: }
1241: } else {
1242: for (j = 0; j < n; j++) {
1243: if (indexn[j] < 0) {
1244: idx += m;
1245: continue;
1246: }
1247: 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);
1248: for (i = 0; i < m; i++) {
1249: if (indexm[i] < 0) {
1250: idx++;
1251: continue;
1252: }
1253: PetscCheck(indexm[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, indexm[i], A->rmap->n - 1);
1254: av[indexn[j] * mat->lda + indexm[i]] += v[idx++];
1255: }
1256: }
1257: }
1258: } else {
1259: if (addv == INSERT_VALUES) {
1260: for (i = 0; i < m; i++) {
1261: if (indexm[i] < 0) {
1262: idx += n;
1263: continue;
1264: }
1265: PetscCheck(indexm[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, indexm[i], A->rmap->n - 1);
1266: for (j = 0; j < n; j++) {
1267: if (indexn[j] < 0) {
1268: idx++;
1269: continue;
1270: }
1271: PetscCheck(indexn[j] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, indexn[j], A->cmap->n - 1);
1272: av[indexn[j] * mat->lda + indexm[i]] = v[idx++];
1273: }
1274: }
1275: } else {
1276: for (i = 0; i < m; i++) {
1277: if (indexm[i] < 0) {
1278: idx += n;
1279: continue;
1280: }
1281: 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);
1282: for (j = 0; j < n; j++) {
1283: if (indexn[j] < 0) {
1284: idx++;
1285: continue;
1286: }
1287: 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);
1288: av[indexn[j] * mat->lda + indexm[i]] += v[idx++];
1289: }
1290: }
1291: }
1292: }
1293: /* hack to prevent unneeded copy to the GPU while returning the array */
1294: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1295: oldf = A->offloadmask;
1296: A->offloadmask = PETSC_OFFLOAD_GPU;
1297: #endif
1298: PetscCall(MatDenseRestoreArray(A, &av));
1299: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1300: A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1301: #endif
1302: PetscFunctionReturn(PETSC_SUCCESS);
1303: }
1305: static PetscErrorCode MatGetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], PetscScalar v[])
1306: {
1307: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1308: const PetscScalar *vv;
1309: PetscInt i, j;
1311: PetscFunctionBegin;
1312: PetscCall(MatDenseGetArrayRead(A, &vv));
1313: /* row-oriented output */
1314: for (i = 0; i < m; i++) {
1315: if (indexm[i] < 0) {
1316: v += n;
1317: continue;
1318: }
1319: 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);
1320: for (j = 0; j < n; j++) {
1321: if (indexn[j] < 0) {
1322: v++;
1323: continue;
1324: }
1325: 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);
1326: *v++ = vv[indexn[j] * mat->lda + indexm[i]];
1327: }
1328: }
1329: PetscCall(MatDenseRestoreArrayRead(A, &vv));
1330: PetscFunctionReturn(PETSC_SUCCESS);
1331: }
1333: PetscErrorCode MatView_Dense_Binary(Mat mat, PetscViewer viewer)
1334: {
1335: PetscBool skipHeader;
1336: PetscViewerFormat format;
1337: PetscInt header[4], M, N, m, lda, i, j, k;
1338: const PetscScalar *v;
1339: PetscScalar *vwork;
1341: PetscFunctionBegin;
1342: PetscCall(PetscViewerSetUp(viewer));
1343: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
1344: PetscCall(PetscViewerGetFormat(viewer, &format));
1345: if (skipHeader) format = PETSC_VIEWER_NATIVE;
1347: PetscCall(MatGetSize(mat, &M, &N));
1349: /* write matrix header */
1350: header[0] = MAT_FILE_CLASSID;
1351: header[1] = M;
1352: header[2] = N;
1353: header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M * N;
1354: if (!skipHeader) PetscCall(PetscViewerBinaryWrite(viewer, header, 4, PETSC_INT));
1356: PetscCall(MatGetLocalSize(mat, &m, NULL));
1357: if (format != PETSC_VIEWER_NATIVE) {
1358: PetscInt nnz = m * N, *iwork;
1359: /* store row lengths for each row */
1360: PetscCall(PetscMalloc1(nnz, &iwork));
1361: for (i = 0; i < m; i++) iwork[i] = N;
1362: PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1363: /* store column indices (zero start index) */
1364: for (k = 0, i = 0; i < m; i++)
1365: for (j = 0; j < N; j++, k++) iwork[k] = j;
1366: PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1367: PetscCall(PetscFree(iwork));
1368: }
1369: /* store matrix values as a dense matrix in row major order */
1370: PetscCall(PetscMalloc1(m * N, &vwork));
1371: PetscCall(MatDenseGetArrayRead(mat, &v));
1372: PetscCall(MatDenseGetLDA(mat, &lda));
1373: for (k = 0, i = 0; i < m; i++)
1374: for (j = 0; j < N; j++, k++) vwork[k] = v[i + lda * j];
1375: PetscCall(MatDenseRestoreArrayRead(mat, &v));
1376: PetscCall(PetscViewerBinaryWriteAll(viewer, vwork, m * N, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1377: PetscCall(PetscFree(vwork));
1378: PetscFunctionReturn(PETSC_SUCCESS);
1379: }
1381: PetscErrorCode MatLoad_Dense_Binary(Mat mat, PetscViewer viewer)
1382: {
1383: PetscBool skipHeader;
1384: PetscInt header[4], M, N, m, nz, lda, i, j, k;
1385: PetscInt rows, cols;
1386: PetscScalar *v, *vwork;
1388: PetscFunctionBegin;
1389: PetscCall(PetscViewerSetUp(viewer));
1390: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
1392: if (!skipHeader) {
1393: PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT));
1394: PetscCheck(header[0] == MAT_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file");
1395: M = header[1];
1396: N = header[2];
1397: PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M);
1398: PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N);
1399: nz = header[3];
1400: PetscCheck(nz == MATRIX_BINARY_FORMAT_DENSE || nz >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Unknown matrix format %" PetscInt_FMT " in file", nz);
1401: } else {
1402: PetscCall(MatGetSize(mat, &M, &N));
1403: 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");
1404: nz = MATRIX_BINARY_FORMAT_DENSE;
1405: }
1407: /* setup global sizes if not set */
1408: if (mat->rmap->N < 0) mat->rmap->N = M;
1409: if (mat->cmap->N < 0) mat->cmap->N = N;
1410: PetscCall(MatSetUp(mat));
1411: /* check if global sizes are correct */
1412: PetscCall(MatGetSize(mat, &rows, &cols));
1413: 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);
1415: PetscCall(MatGetSize(mat, NULL, &N));
1416: PetscCall(MatGetLocalSize(mat, &m, NULL));
1417: PetscCall(MatDenseGetArray(mat, &v));
1418: PetscCall(MatDenseGetLDA(mat, &lda));
1419: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */
1420: PetscInt nnz = m * N;
1421: /* read in matrix values */
1422: PetscCall(PetscMalloc1(nnz, &vwork));
1423: PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1424: /* store values in column major order */
1425: for (j = 0; j < N; j++)
1426: for (i = 0; i < m; i++) v[i + lda * j] = vwork[i * N + j];
1427: PetscCall(PetscFree(vwork));
1428: } else { /* matrix in file is sparse format */
1429: PetscInt nnz = 0, *rlens, *icols;
1430: /* read in row lengths */
1431: PetscCall(PetscMalloc1(m, &rlens));
1432: PetscCall(PetscViewerBinaryReadAll(viewer, rlens, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1433: for (i = 0; i < m; i++) nnz += rlens[i];
1434: /* read in column indices and values */
1435: PetscCall(PetscMalloc2(nnz, &icols, nnz, &vwork));
1436: PetscCall(PetscViewerBinaryReadAll(viewer, icols, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1437: PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1438: /* store values in column major order */
1439: for (k = 0, i = 0; i < m; i++)
1440: for (j = 0; j < rlens[i]; j++, k++) v[i + lda * icols[k]] = vwork[k];
1441: PetscCall(PetscFree(rlens));
1442: PetscCall(PetscFree2(icols, vwork));
1443: }
1444: PetscCall(MatDenseRestoreArray(mat, &v));
1445: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
1446: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
1447: PetscFunctionReturn(PETSC_SUCCESS);
1448: }
1450: static PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1451: {
1452: PetscBool isbinary, ishdf5;
1454: PetscFunctionBegin;
1457: /* force binary viewer to load .info file if it has not yet done so */
1458: PetscCall(PetscViewerSetUp(viewer));
1459: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1460: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1461: if (isbinary) {
1462: PetscCall(MatLoad_Dense_Binary(newMat, viewer));
1463: } else if (ishdf5) {
1464: #if defined(PETSC_HAVE_HDF5)
1465: PetscCall(MatLoad_Dense_HDF5(newMat, viewer));
1466: #else
1467: SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1468: #endif
1469: } else {
1470: 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);
1471: }
1472: PetscFunctionReturn(PETSC_SUCCESS);
1473: }
1475: static PetscErrorCode MatView_SeqDense_ASCII(Mat A, PetscViewer viewer)
1476: {
1477: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1478: PetscInt i, j;
1479: const char *name;
1480: PetscScalar *v, *av;
1481: PetscViewerFormat format;
1482: #if defined(PETSC_USE_COMPLEX)
1483: PetscBool allreal = PETSC_TRUE;
1484: #endif
1486: PetscFunctionBegin;
1487: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&av));
1488: PetscCall(PetscViewerGetFormat(viewer, &format));
1489: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1490: PetscFunctionReturn(PETSC_SUCCESS); /* do nothing for now */
1491: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1492: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1493: for (i = 0; i < A->rmap->n; i++) {
1494: v = av + i;
1495: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i));
1496: for (j = 0; j < A->cmap->n; j++) {
1497: #if defined(PETSC_USE_COMPLEX)
1498: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1499: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", j, (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1500: } else if (PetscRealPart(*v)) {
1501: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)PetscRealPart(*v)));
1502: }
1503: #else
1504: if (*v) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)*v));
1505: #endif
1506: v += a->lda;
1507: }
1508: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1509: }
1510: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1511: } else {
1512: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1513: #if defined(PETSC_USE_COMPLEX)
1514: /* determine if matrix has all real values */
1515: for (j = 0; j < A->cmap->n; j++) {
1516: v = av + j * a->lda;
1517: for (i = 0; i < A->rmap->n; i++) {
1518: if (PetscImaginaryPart(v[i])) {
1519: allreal = PETSC_FALSE;
1520: break;
1521: }
1522: }
1523: }
1524: #endif
1525: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1526: PetscCall(PetscObjectGetName((PetscObject)A, &name));
1527: PetscCall(PetscViewerASCIIPrintf(viewer, "%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n", A->rmap->n, A->cmap->n));
1528: PetscCall(PetscViewerASCIIPrintf(viewer, "%s = zeros(%" PetscInt_FMT ",%" PetscInt_FMT ");\n", name, A->rmap->n, A->cmap->n));
1529: PetscCall(PetscViewerASCIIPrintf(viewer, "%s = [\n", name));
1530: }
1532: for (i = 0; i < A->rmap->n; i++) {
1533: v = av + i;
1534: for (j = 0; j < A->cmap->n; j++) {
1535: #if defined(PETSC_USE_COMPLEX)
1536: if (allreal) {
1537: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)PetscRealPart(*v)));
1538: } else {
1539: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei ", (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1540: }
1541: #else
1542: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)*v));
1543: #endif
1544: v += a->lda;
1545: }
1546: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1547: }
1548: if (format == PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer, "];\n"));
1549: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1550: }
1551: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&av));
1552: PetscCall(PetscViewerFlush(viewer));
1553: PetscFunctionReturn(PETSC_SUCCESS);
1554: }
1556: #include <petscdraw.h>
1557: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw, void *Aa)
1558: {
1559: Mat A = (Mat)Aa;
1560: PetscInt m = A->rmap->n, n = A->cmap->n, i, j;
1561: int color = PETSC_DRAW_WHITE;
1562: const PetscScalar *v;
1563: PetscViewer viewer;
1564: PetscReal xl, yl, xr, yr, x_l, x_r, y_l, y_r;
1565: PetscViewerFormat format;
1567: PetscFunctionBegin;
1568: PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer));
1569: PetscCall(PetscViewerGetFormat(viewer, &format));
1570: PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));
1572: /* Loop over matrix elements drawing boxes */
1573: PetscCall(MatDenseGetArrayRead(A, &v));
1574: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1575: PetscDrawCollectiveBegin(draw);
1576: /* Blue for negative and Red for positive */
1577: for (j = 0; j < n; j++) {
1578: x_l = j;
1579: x_r = x_l + 1.0;
1580: for (i = 0; i < m; i++) {
1581: y_l = m - i - 1.0;
1582: y_r = y_l + 1.0;
1583: if (PetscRealPart(v[j * m + i]) > 0.) color = PETSC_DRAW_RED;
1584: else if (PetscRealPart(v[j * m + i]) < 0.) color = PETSC_DRAW_BLUE;
1585: else continue;
1586: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1587: }
1588: }
1589: PetscDrawCollectiveEnd(draw);
1590: } else {
1591: /* use contour shading to indicate magnitude of values */
1592: /* first determine max of all nonzero values */
1593: PetscReal minv = 0.0, maxv = 0.0;
1594: PetscDraw popup;
1596: for (i = 0; i < m * n; i++) {
1597: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1598: }
1599: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1600: PetscCall(PetscDrawGetPopup(draw, &popup));
1601: PetscCall(PetscDrawScalePopup(popup, minv, maxv));
1603: PetscDrawCollectiveBegin(draw);
1604: for (j = 0; j < n; j++) {
1605: x_l = j;
1606: x_r = x_l + 1.0;
1607: for (i = 0; i < m; i++) {
1608: y_l = m - i - 1.0;
1609: y_r = y_l + 1.0;
1610: color = PetscDrawRealToColor(PetscAbsScalar(v[j * m + i]), minv, maxv);
1611: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1612: }
1613: }
1614: PetscDrawCollectiveEnd(draw);
1615: }
1616: PetscCall(MatDenseRestoreArrayRead(A, &v));
1617: PetscFunctionReturn(PETSC_SUCCESS);
1618: }
1620: static PetscErrorCode MatView_SeqDense_Draw(Mat A, PetscViewer viewer)
1621: {
1622: PetscDraw draw;
1623: PetscBool isnull;
1624: PetscReal xr, yr, xl, yl, h, w;
1626: PetscFunctionBegin;
1627: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1628: PetscCall(PetscDrawIsNull(draw, &isnull));
1629: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
1631: xr = A->cmap->n;
1632: yr = A->rmap->n;
1633: h = yr / 10.0;
1634: w = xr / 10.0;
1635: xr += w;
1636: yr += h;
1637: xl = -w;
1638: yl = -h;
1639: PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr));
1640: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer));
1641: PetscCall(PetscDrawZoom(draw, MatView_SeqDense_Draw_Zoom, A));
1642: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL));
1643: PetscCall(PetscDrawSave(draw));
1644: PetscFunctionReturn(PETSC_SUCCESS);
1645: }
1647: PetscErrorCode MatView_SeqDense(Mat A, PetscViewer viewer)
1648: {
1649: PetscBool iascii, isbinary, isdraw;
1651: PetscFunctionBegin;
1652: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1653: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1654: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1655: if (iascii) PetscCall(MatView_SeqDense_ASCII(A, viewer));
1656: else if (isbinary) PetscCall(MatView_Dense_Binary(A, viewer));
1657: else if (isdraw) PetscCall(MatView_SeqDense_Draw(A, viewer));
1658: PetscFunctionReturn(PETSC_SUCCESS);
1659: }
1661: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A, const PetscScalar *array)
1662: {
1663: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1665: PetscFunctionBegin;
1666: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1667: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1668: PetscCheck(!a->unplacedarray, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreArray() first");
1669: a->unplacedarray = a->v;
1670: a->unplaced_user_alloc = a->user_alloc;
1671: a->v = (PetscScalar *)array;
1672: a->user_alloc = PETSC_TRUE;
1673: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1674: A->offloadmask = PETSC_OFFLOAD_CPU;
1675: #endif
1676: PetscFunctionReturn(PETSC_SUCCESS);
1677: }
1679: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1680: {
1681: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1683: PetscFunctionBegin;
1684: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1685: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1686: a->v = a->unplacedarray;
1687: a->user_alloc = a->unplaced_user_alloc;
1688: a->unplacedarray = NULL;
1689: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1690: A->offloadmask = PETSC_OFFLOAD_CPU;
1691: #endif
1692: PetscFunctionReturn(PETSC_SUCCESS);
1693: }
1695: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A, const PetscScalar *array)
1696: {
1697: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1699: PetscFunctionBegin;
1700: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1701: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1702: if (!a->user_alloc) PetscCall(PetscFree(a->v));
1703: a->v = (PetscScalar *)array;
1704: a->user_alloc = PETSC_FALSE;
1705: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1706: A->offloadmask = PETSC_OFFLOAD_CPU;
1707: #endif
1708: PetscFunctionReturn(PETSC_SUCCESS);
1709: }
1711: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1712: {
1713: Mat_SeqDense *l = (Mat_SeqDense *)mat->data;
1715: PetscFunctionBegin;
1716: PetscCall(PetscLogObjectState((PetscObject)mat, "Rows %" PetscInt_FMT " Cols %" PetscInt_FMT, mat->rmap->n, mat->cmap->n));
1717: PetscCall(VecDestroy(&l->qrrhs));
1718: PetscCall(PetscFree(l->tau));
1719: PetscCall(PetscFree(l->pivots));
1720: PetscCall(PetscFree(l->fwork));
1721: if (!l->user_alloc) PetscCall(PetscFree(l->v));
1722: if (!l->unplaced_user_alloc) PetscCall(PetscFree(l->unplacedarray));
1723: PetscCheck(!l->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1724: PetscCheck(!l->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1725: PetscCall(VecDestroy(&l->cvec));
1726: PetscCall(MatDestroy(&l->cmat));
1727: PetscCall(PetscFree(mat->data));
1729: PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL));
1730: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactor_C", NULL));
1731: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorSymbolic_C", NULL));
1732: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorNumeric_C", NULL));
1733: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetLDA_C", NULL));
1734: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseSetLDA_C", NULL));
1735: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArray_C", NULL));
1736: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArray_C", NULL));
1737: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDensePlaceArray_C", NULL));
1738: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseResetArray_C", NULL));
1739: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseReplaceArray_C", NULL));
1740: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayRead_C", NULL));
1741: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayRead_C", NULL));
1742: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayWrite_C", NULL));
1743: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayWrite_C", NULL));
1744: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqaij_C", NULL));
1745: #if defined(PETSC_HAVE_ELEMENTAL)
1746: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_elemental_C", NULL));
1747: #endif
1748: #if defined(PETSC_HAVE_SCALAPACK)
1749: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_scalapack_C", NULL));
1750: #endif
1751: #if defined(PETSC_HAVE_CUDA)
1752: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensecuda_C", NULL));
1753: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", NULL));
1754: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdense_C", NULL));
1755: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensecuda_C", NULL));
1756: #endif
1757: #if defined(PETSC_HAVE_HIP)
1758: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensehip_C", NULL));
1759: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", NULL));
1760: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdense_C", NULL));
1761: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensehip_C", NULL));
1762: #endif
1763: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatSeqDenseSetPreallocation_C", NULL));
1764: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqaij_seqdense_C", NULL));
1765: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdense_C", NULL));
1766: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqbaij_seqdense_C", NULL));
1767: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqsbaij_seqdense_C", NULL));
1769: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumn_C", NULL));
1770: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumn_C", NULL));
1771: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVec_C", NULL));
1772: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVec_C", NULL));
1773: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecRead_C", NULL));
1774: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecRead_C", NULL));
1775: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecWrite_C", NULL));
1776: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecWrite_C", NULL));
1777: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetSubMatrix_C", NULL));
1778: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreSubMatrix_C", NULL));
1779: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultAddColumnRange_C", NULL));
1780: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeColumnRange_C", NULL));
1781: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatMultHermitianTransposeAddColumnRange_C", NULL));
1782: PetscFunctionReturn(PETSC_SUCCESS);
1783: }
1785: static PetscErrorCode MatTranspose_SeqDense(Mat A, MatReuse reuse, Mat *matout)
1786: {
1787: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1788: PetscInt k, j, m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1789: PetscScalar *v, tmp;
1791: PetscFunctionBegin;
1792: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *matout));
1793: if (reuse == MAT_INPLACE_MATRIX) {
1794: if (m == n) { /* in place transpose */
1795: PetscCall(MatDenseGetArray(A, &v));
1796: for (j = 0; j < m; j++) {
1797: for (k = 0; k < j; k++) {
1798: tmp = v[j + k * M];
1799: v[j + k * M] = v[k + j * M];
1800: v[k + j * M] = tmp;
1801: }
1802: }
1803: PetscCall(MatDenseRestoreArray(A, &v));
1804: } else { /* reuse memory, temporary allocates new memory */
1805: PetscScalar *v2;
1806: PetscLayout tmplayout;
1808: PetscCall(PetscMalloc1((size_t)m * n, &v2));
1809: PetscCall(MatDenseGetArray(A, &v));
1810: for (j = 0; j < n; j++) {
1811: for (k = 0; k < m; k++) v2[j + (size_t)k * n] = v[k + (size_t)j * M];
1812: }
1813: PetscCall(PetscArraycpy(v, v2, (size_t)m * n));
1814: PetscCall(PetscFree(v2));
1815: PetscCall(MatDenseRestoreArray(A, &v));
1816: /* cleanup size dependent quantities */
1817: PetscCall(VecDestroy(&mat->cvec));
1818: PetscCall(MatDestroy(&mat->cmat));
1819: PetscCall(PetscFree(mat->pivots));
1820: PetscCall(PetscFree(mat->fwork));
1821: /* swap row/col layouts */
1822: mat->lda = n;
1823: tmplayout = A->rmap;
1824: A->rmap = A->cmap;
1825: A->cmap = tmplayout;
1826: }
1827: } else { /* out-of-place transpose */
1828: Mat tmat;
1829: Mat_SeqDense *tmatd;
1830: PetscScalar *v2;
1831: PetscInt M2;
1833: if (reuse == MAT_INITIAL_MATRIX) {
1834: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &tmat));
1835: PetscCall(MatSetSizes(tmat, A->cmap->n, A->rmap->n, A->cmap->n, A->rmap->n));
1836: PetscCall(MatSetType(tmat, ((PetscObject)A)->type_name));
1837: PetscCall(MatSeqDenseSetPreallocation(tmat, NULL));
1838: } else tmat = *matout;
1840: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&v));
1841: PetscCall(MatDenseGetArray(tmat, &v2));
1842: tmatd = (Mat_SeqDense *)tmat->data;
1843: M2 = tmatd->lda;
1844: for (j = 0; j < n; j++) {
1845: for (k = 0; k < m; k++) v2[j + k * M2] = v[k + j * M];
1846: }
1847: PetscCall(MatDenseRestoreArray(tmat, &v2));
1848: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&v));
1849: PetscCall(MatAssemblyBegin(tmat, MAT_FINAL_ASSEMBLY));
1850: PetscCall(MatAssemblyEnd(tmat, MAT_FINAL_ASSEMBLY));
1851: *matout = tmat;
1852: }
1853: PetscFunctionReturn(PETSC_SUCCESS);
1854: }
1856: static PetscErrorCode MatEqual_SeqDense(Mat A1, Mat A2, PetscBool *flg)
1857: {
1858: Mat_SeqDense *mat1 = (Mat_SeqDense *)A1->data;
1859: Mat_SeqDense *mat2 = (Mat_SeqDense *)A2->data;
1860: PetscInt i;
1861: const PetscScalar *v1, *v2;
1863: PetscFunctionBegin;
1864: if (A1->rmap->n != A2->rmap->n) {
1865: *flg = PETSC_FALSE;
1866: PetscFunctionReturn(PETSC_SUCCESS);
1867: }
1868: if (A1->cmap->n != A2->cmap->n) {
1869: *flg = PETSC_FALSE;
1870: PetscFunctionReturn(PETSC_SUCCESS);
1871: }
1872: PetscCall(MatDenseGetArrayRead(A1, &v1));
1873: PetscCall(MatDenseGetArrayRead(A2, &v2));
1874: for (i = 0; i < A1->cmap->n; i++) {
1875: PetscCall(PetscArraycmp(v1, v2, A1->rmap->n, flg));
1876: if (*flg == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS);
1877: v1 += mat1->lda;
1878: v2 += mat2->lda;
1879: }
1880: PetscCall(MatDenseRestoreArrayRead(A1, &v1));
1881: PetscCall(MatDenseRestoreArrayRead(A2, &v2));
1882: *flg = PETSC_TRUE;
1883: PetscFunctionReturn(PETSC_SUCCESS);
1884: }
1886: PetscErrorCode MatGetDiagonal_SeqDense(Mat A, Vec v)
1887: {
1888: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1889: PetscInt i, n, len;
1890: PetscScalar *x;
1891: const PetscScalar *vv;
1893: PetscFunctionBegin;
1894: PetscCall(VecGetSize(v, &n));
1895: PetscCall(VecGetArray(v, &x));
1896: len = PetscMin(A->rmap->n, A->cmap->n);
1897: PetscCall(MatDenseGetArrayRead(A, &vv));
1898: PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming mat and vec");
1899: for (i = 0; i < len; i++) x[i] = vv[i * mat->lda + i];
1900: PetscCall(MatDenseRestoreArrayRead(A, &vv));
1901: PetscCall(VecRestoreArray(v, &x));
1902: PetscFunctionReturn(PETSC_SUCCESS);
1903: }
1905: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A, Vec ll, Vec rr)
1906: {
1907: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1908: const PetscScalar *l, *r;
1909: PetscScalar x, *v, *vv;
1910: PetscInt i, j, m = A->rmap->n, n = A->cmap->n;
1912: PetscFunctionBegin;
1913: PetscCall(MatDenseGetArray(A, &vv));
1914: if (ll) {
1915: PetscCall(VecGetSize(ll, &m));
1916: PetscCall(VecGetArrayRead(ll, &l));
1917: PetscCheck(m == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vec wrong size");
1918: for (i = 0; i < m; i++) {
1919: x = l[i];
1920: v = vv + i;
1921: for (j = 0; j < n; j++) {
1922: (*v) *= x;
1923: v += mat->lda;
1924: }
1925: }
1926: PetscCall(VecRestoreArrayRead(ll, &l));
1927: PetscCall(PetscLogFlops(1.0 * n * m));
1928: }
1929: if (rr) {
1930: PetscCall(VecGetSize(rr, &n));
1931: PetscCall(VecGetArrayRead(rr, &r));
1932: PetscCheck(n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vec wrong size");
1933: for (i = 0; i < n; i++) {
1934: x = r[i];
1935: v = vv + i * mat->lda;
1936: for (j = 0; j < m; j++) (*v++) *= x;
1937: }
1938: PetscCall(VecRestoreArrayRead(rr, &r));
1939: PetscCall(PetscLogFlops(1.0 * n * m));
1940: }
1941: PetscCall(MatDenseRestoreArray(A, &vv));
1942: PetscFunctionReturn(PETSC_SUCCESS);
1943: }
1945: PetscErrorCode MatNorm_SeqDense(Mat A, NormType type, PetscReal *nrm)
1946: {
1947: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1948: PetscScalar *v, *vv;
1949: PetscReal sum = 0.0;
1950: PetscInt lda, m = A->rmap->n, i, j;
1952: PetscFunctionBegin;
1953: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&vv));
1954: PetscCall(MatDenseGetLDA(A, &lda));
1955: v = vv;
1956: if (type == NORM_FROBENIUS) {
1957: if (lda > m) {
1958: for (j = 0; j < A->cmap->n; j++) {
1959: v = vv + j * lda;
1960: for (i = 0; i < m; i++) {
1961: sum += PetscRealPart(PetscConj(*v) * (*v));
1962: v++;
1963: }
1964: }
1965: } else {
1966: #if defined(PETSC_USE_REAL___FP16)
1967: PetscBLASInt one = 1, cnt = A->cmap->n * A->rmap->n;
1968: PetscCallBLAS("BLASnrm2", *nrm = BLASnrm2_(&cnt, v, &one));
1969: }
1970: #else
1971: for (i = 0; i < A->cmap->n * A->rmap->n; i++) {
1972: sum += PetscRealPart(PetscConj(*v) * (*v));
1973: v++;
1974: }
1975: }
1976: *nrm = PetscSqrtReal(sum);
1977: #endif
1978: PetscCall(PetscLogFlops(2.0 * A->cmap->n * A->rmap->n));
1979: } else if (type == NORM_1) {
1980: *nrm = 0.0;
1981: for (j = 0; j < A->cmap->n; j++) {
1982: v = vv + j * mat->lda;
1983: sum = 0.0;
1984: for (i = 0; i < A->rmap->n; i++) {
1985: sum += PetscAbsScalar(*v);
1986: v++;
1987: }
1988: if (sum > *nrm) *nrm = sum;
1989: }
1990: PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
1991: } else if (type == NORM_INFINITY) {
1992: *nrm = 0.0;
1993: for (j = 0; j < A->rmap->n; j++) {
1994: v = vv + j;
1995: sum = 0.0;
1996: for (i = 0; i < A->cmap->n; i++) {
1997: sum += PetscAbsScalar(*v);
1998: v += mat->lda;
1999: }
2000: if (sum > *nrm) *nrm = sum;
2001: }
2002: PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
2003: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No two norm");
2004: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&vv));
2005: PetscFunctionReturn(PETSC_SUCCESS);
2006: }
2008: static PetscErrorCode MatSetOption_SeqDense(Mat A, MatOption op, PetscBool flg)
2009: {
2010: Mat_SeqDense *aij = (Mat_SeqDense *)A->data;
2012: PetscFunctionBegin;
2013: switch (op) {
2014: case MAT_ROW_ORIENTED:
2015: aij->roworiented = flg;
2016: break;
2017: case MAT_NEW_NONZERO_LOCATIONS:
2018: case MAT_NEW_NONZERO_LOCATION_ERR:
2019: case MAT_NEW_NONZERO_ALLOCATION_ERR:
2020: case MAT_FORCE_DIAGONAL_ENTRIES:
2021: case MAT_KEEP_NONZERO_PATTERN:
2022: case MAT_IGNORE_OFF_PROC_ENTRIES:
2023: case MAT_USE_HASH_TABLE:
2024: case MAT_IGNORE_ZERO_ENTRIES:
2025: case MAT_IGNORE_LOWER_TRIANGULAR:
2026: case MAT_SORTED_FULL:
2027: PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op]));
2028: break;
2029: case MAT_SPD:
2030: case MAT_SYMMETRIC:
2031: case MAT_STRUCTURALLY_SYMMETRIC:
2032: case MAT_HERMITIAN:
2033: case MAT_SYMMETRY_ETERNAL:
2034: case MAT_STRUCTURAL_SYMMETRY_ETERNAL:
2035: case MAT_SPD_ETERNAL:
2036: break;
2037: default:
2038: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unknown option %s", MatOptions[op]);
2039: }
2040: PetscFunctionReturn(PETSC_SUCCESS);
2041: }
2043: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
2044: {
2045: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
2046: PetscInt lda = l->lda, m = A->rmap->n, n = A->cmap->n, j;
2047: PetscScalar *v;
2049: PetscFunctionBegin;
2050: PetscCall(MatDenseGetArrayWrite(A, &v));
2051: if (lda > m) {
2052: for (j = 0; j < n; j++) PetscCall(PetscArrayzero(v + j * lda, m));
2053: } else {
2054: PetscCall(PetscArrayzero(v, PetscInt64Mult(m, n)));
2055: }
2056: PetscCall(MatDenseRestoreArrayWrite(A, &v));
2057: PetscFunctionReturn(PETSC_SUCCESS);
2058: }
2060: static PetscErrorCode MatZeroRows_SeqDense(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2061: {
2062: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
2063: PetscInt m = l->lda, n = A->cmap->n, i, j;
2064: PetscScalar *slot, *bb, *v;
2065: const PetscScalar *xx;
2067: PetscFunctionBegin;
2068: if (PetscDefined(USE_DEBUG)) {
2069: for (i = 0; i < N; i++) {
2070: PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
2071: 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);
2072: }
2073: }
2074: if (!N) PetscFunctionReturn(PETSC_SUCCESS);
2076: /* fix right-hand side if needed */
2077: if (x && b) {
2078: PetscCall(VecGetArrayRead(x, &xx));
2079: PetscCall(VecGetArray(b, &bb));
2080: for (i = 0; i < N; i++) bb[rows[i]] = diag * xx[rows[i]];
2081: PetscCall(VecRestoreArrayRead(x, &xx));
2082: PetscCall(VecRestoreArray(b, &bb));
2083: }
2085: PetscCall(MatDenseGetArray(A, &v));
2086: for (i = 0; i < N; i++) {
2087: slot = v + rows[i];
2088: for (j = 0; j < n; j++) {
2089: *slot = 0.0;
2090: slot += m;
2091: }
2092: }
2093: if (diag != 0.0) {
2094: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
2095: for (i = 0; i < N; i++) {
2096: slot = v + (m + 1) * rows[i];
2097: *slot = diag;
2098: }
2099: }
2100: PetscCall(MatDenseRestoreArray(A, &v));
2101: PetscFunctionReturn(PETSC_SUCCESS);
2102: }
2104: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A, PetscInt *lda)
2105: {
2106: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2108: PetscFunctionBegin;
2109: *lda = mat->lda;
2110: PetscFunctionReturn(PETSC_SUCCESS);
2111: }
2113: PetscErrorCode MatDenseGetArray_SeqDense(Mat A, PetscScalar **array)
2114: {
2115: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2117: PetscFunctionBegin;
2118: PetscCheck(!mat->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
2119: *array = mat->v;
2120: PetscFunctionReturn(PETSC_SUCCESS);
2121: }
2123: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A, PetscScalar **array)
2124: {
2125: PetscFunctionBegin;
2126: if (array) *array = NULL;
2127: PetscFunctionReturn(PETSC_SUCCESS);
2128: }
2130: /*@
2131: MatDenseGetLDA - gets the leading dimension of the array returned from `MatDenseGetArray()`
2133: Not Collective
2135: Input Parameter:
2136: . A - a `MATDENSE` or `MATDENSECUDA` matrix
2138: Output Parameter:
2139: . lda - the leading dimension
2141: Level: intermediate
2143: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseSetLDA()`
2144: @*/
2145: PetscErrorCode MatDenseGetLDA(Mat A, PetscInt *lda)
2146: {
2147: PetscFunctionBegin;
2149: PetscAssertPointer(lda, 2);
2150: MatCheckPreallocated(A, 1);
2151: PetscUseMethod(A, "MatDenseGetLDA_C", (Mat, PetscInt *), (A, lda));
2152: PetscFunctionReturn(PETSC_SUCCESS);
2153: }
2155: /*@
2156: MatDenseSetLDA - Sets the leading dimension of the array used by the `MATDENSE` matrix
2158: Collective if the matrix layouts have not yet been setup
2160: Input Parameters:
2161: + A - a `MATDENSE` or `MATDENSECUDA` matrix
2162: - lda - the leading dimension
2164: Level: intermediate
2166: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetLDA()`
2167: @*/
2168: PetscErrorCode MatDenseSetLDA(Mat A, PetscInt lda)
2169: {
2170: PetscFunctionBegin;
2172: PetscTryMethod(A, "MatDenseSetLDA_C", (Mat, PetscInt), (A, lda));
2173: PetscFunctionReturn(PETSC_SUCCESS);
2174: }
2176: /*@C
2177: MatDenseGetArray - gives read-write access to the array where the data for a `MATDENSE` matrix is stored
2179: Logically Collective
2181: Input Parameter:
2182: . A - a dense matrix
2184: Output Parameter:
2185: . array - pointer to the data
2187: Level: intermediate
2189: Fortran Notes:
2190: `MatDenseGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `MatDenseGetArrayF90()`
2192: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2193: @*/
2194: PetscErrorCode MatDenseGetArray(Mat A, PetscScalar **array)
2195: {
2196: PetscFunctionBegin;
2198: PetscAssertPointer(array, 2);
2199: PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2200: PetscFunctionReturn(PETSC_SUCCESS);
2201: }
2203: /*@C
2204: MatDenseRestoreArray - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArray()`
2206: Logically Collective
2208: Input Parameters:
2209: + A - a dense matrix
2210: - array - pointer to the data (may be `NULL`)
2212: Level: intermediate
2214: Fortran Notes:
2215: `MatDenseRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `MatDenseRestoreArrayF90()`
2217: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2218: @*/
2219: PetscErrorCode MatDenseRestoreArray(Mat A, PetscScalar **array)
2220: {
2221: PetscFunctionBegin;
2223: if (array) PetscAssertPointer(array, 2);
2224: PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2225: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2226: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2227: A->offloadmask = PETSC_OFFLOAD_CPU;
2228: #endif
2229: PetscFunctionReturn(PETSC_SUCCESS);
2230: }
2232: /*@C
2233: MatDenseGetArrayRead - gives read-only access to the array where the data for a `MATDENSE` matrix is stored
2235: Not Collective
2237: Input Parameter:
2238: . A - a dense matrix
2240: Output Parameter:
2241: . array - pointer to the data
2243: Level: intermediate
2245: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2246: @*/
2247: PetscErrorCode MatDenseGetArrayRead(Mat A, const PetscScalar **array)
2248: {
2249: PetscFunctionBegin;
2251: PetscAssertPointer(array, 2);
2252: PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2253: PetscFunctionReturn(PETSC_SUCCESS);
2254: }
2256: /*@C
2257: MatDenseRestoreArrayRead - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayRead()`
2259: Not Collective
2261: Input Parameters:
2262: + A - a dense matrix
2263: - array - pointer to the data (may be `NULL`)
2265: Level: intermediate
2267: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2268: @*/
2269: PetscErrorCode MatDenseRestoreArrayRead(Mat A, const PetscScalar **array)
2270: {
2271: PetscFunctionBegin;
2273: if (array) PetscAssertPointer(array, 2);
2274: PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2275: PetscFunctionReturn(PETSC_SUCCESS);
2276: }
2278: /*@C
2279: MatDenseGetArrayWrite - gives write-only access to the array where the data for a `MATDENSE` matrix is stored
2281: Not Collective
2283: Input Parameter:
2284: . A - a dense matrix
2286: Output Parameter:
2287: . array - pointer to the data
2289: Level: intermediate
2291: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2292: @*/
2293: PetscErrorCode MatDenseGetArrayWrite(Mat A, PetscScalar **array)
2294: {
2295: PetscFunctionBegin;
2297: PetscAssertPointer(array, 2);
2298: PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2299: PetscFunctionReturn(PETSC_SUCCESS);
2300: }
2302: /*@C
2303: MatDenseRestoreArrayWrite - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayWrite()`
2305: Not Collective
2307: Input Parameters:
2308: + A - a dense matrix
2309: - array - pointer to the data (may be `NULL`)
2311: Level: intermediate
2313: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2314: @*/
2315: PetscErrorCode MatDenseRestoreArrayWrite(Mat A, PetscScalar **array)
2316: {
2317: PetscFunctionBegin;
2319: if (array) PetscAssertPointer(array, 2);
2320: PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2321: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2322: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2323: A->offloadmask = PETSC_OFFLOAD_CPU;
2324: #endif
2325: PetscFunctionReturn(PETSC_SUCCESS);
2326: }
2328: /*@C
2329: MatDenseGetArrayAndMemType - gives read-write access to the array where the data for a `MATDENSE` matrix is stored
2331: Logically Collective
2333: Input Parameter:
2334: . A - a dense matrix
2336: Output Parameters:
2337: + array - pointer to the data
2338: - mtype - memory type of the returned pointer
2340: Level: intermediate
2342: Note:
2343: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2344: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2346: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArrayRead()`,
2347: `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2348: @*/
2349: PetscErrorCode MatDenseGetArrayAndMemType(Mat A, PetscScalar **array, PetscMemType *mtype)
2350: {
2351: PetscBool isMPI;
2353: PetscFunctionBegin;
2355: PetscAssertPointer(array, 2);
2356: 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 */
2357: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2358: if (isMPI) {
2359: /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2360: PetscCall(MatDenseGetArrayAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2361: } else {
2362: PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);
2364: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayAndMemType_C", &fptr));
2365: if (fptr) {
2366: PetscCall((*fptr)(A, array, mtype));
2367: } else {
2368: PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2369: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2370: }
2371: }
2372: PetscFunctionReturn(PETSC_SUCCESS);
2373: }
2375: /*@C
2376: MatDenseRestoreArrayAndMemType - returns access to the array that is obtained by `MatDenseGetArrayAndMemType()`
2378: Logically Collective
2380: Input Parameters:
2381: + A - a dense matrix
2382: - array - pointer to the data
2384: Level: intermediate
2386: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2387: @*/
2388: PetscErrorCode MatDenseRestoreArrayAndMemType(Mat A, PetscScalar **array)
2389: {
2390: PetscBool isMPI;
2392: PetscFunctionBegin;
2394: PetscAssertPointer(array, 2);
2395: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2396: if (isMPI) {
2397: PetscCall(MatDenseRestoreArrayAndMemType(((Mat_MPIDense *)A->data)->A, array));
2398: } else {
2399: PetscErrorCode (*fptr)(Mat, PetscScalar **);
2401: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayAndMemType_C", &fptr));
2402: if (fptr) {
2403: PetscCall((*fptr)(A, array));
2404: } else {
2405: PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2406: }
2407: *array = NULL;
2408: }
2409: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2410: PetscFunctionReturn(PETSC_SUCCESS);
2411: }
2413: /*@C
2414: MatDenseGetArrayReadAndMemType - gives read-only access to the array where the data for a `MATDENSE` matrix is stored
2416: Logically Collective
2418: Input Parameter:
2419: . A - a dense matrix
2421: Output Parameters:
2422: + array - pointer to the data
2423: - mtype - memory type of the returned pointer
2425: Level: intermediate
2427: Note:
2428: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2429: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2431: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`,
2432: `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2433: @*/
2434: PetscErrorCode MatDenseGetArrayReadAndMemType(Mat A, const PetscScalar **array, PetscMemType *mtype)
2435: {
2436: PetscBool isMPI;
2438: PetscFunctionBegin;
2440: PetscAssertPointer(array, 2);
2441: 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 */
2442: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2443: if (isMPI) { /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2444: PetscCall(MatDenseGetArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2445: } else {
2446: PetscErrorCode (*fptr)(Mat, const PetscScalar **, PetscMemType *);
2448: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayReadAndMemType_C", &fptr));
2449: if (fptr) {
2450: PetscCall((*fptr)(A, array, mtype));
2451: } else {
2452: PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2453: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2454: }
2455: }
2456: PetscFunctionReturn(PETSC_SUCCESS);
2457: }
2459: /*@C
2460: MatDenseRestoreArrayReadAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()`
2462: Logically Collective
2464: Input Parameters:
2465: + A - a dense matrix
2466: - array - pointer to the data
2468: Level: intermediate
2470: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2471: @*/
2472: PetscErrorCode MatDenseRestoreArrayReadAndMemType(Mat A, const PetscScalar **array)
2473: {
2474: PetscBool isMPI;
2476: PetscFunctionBegin;
2478: PetscAssertPointer(array, 2);
2479: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2480: if (isMPI) {
2481: PetscCall(MatDenseRestoreArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array));
2482: } else {
2483: PetscErrorCode (*fptr)(Mat, const PetscScalar **);
2485: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayReadAndMemType_C", &fptr));
2486: if (fptr) {
2487: PetscCall((*fptr)(A, array));
2488: } else {
2489: PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2490: }
2491: *array = NULL;
2492: }
2493: PetscFunctionReturn(PETSC_SUCCESS);
2494: }
2496: /*@C
2497: MatDenseGetArrayWriteAndMemType - gives write-only access to the array where the data for a `MATDENSE` matrix is stored
2499: Logically Collective
2501: Input Parameter:
2502: . A - a dense matrix
2504: Output Parameters:
2505: + array - pointer to the data
2506: - mtype - memory type of the returned pointer
2508: Level: intermediate
2510: Note:
2511: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2512: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2514: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWriteAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayRead()`,
2515: `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2516: @*/
2517: PetscErrorCode MatDenseGetArrayWriteAndMemType(Mat A, PetscScalar **array, PetscMemType *mtype)
2518: {
2519: PetscBool isMPI;
2521: PetscFunctionBegin;
2523: PetscAssertPointer(array, 2);
2524: 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 */
2525: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2526: if (isMPI) {
2527: PetscCall(MatDenseGetArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2528: } else {
2529: PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);
2531: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayWriteAndMemType_C", &fptr));
2532: if (fptr) {
2533: PetscCall((*fptr)(A, array, mtype));
2534: } else {
2535: PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2536: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2537: }
2538: }
2539: PetscFunctionReturn(PETSC_SUCCESS);
2540: }
2542: /*@C
2543: MatDenseRestoreArrayWriteAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()`
2545: Logically Collective
2547: Input Parameters:
2548: + A - a dense matrix
2549: - array - pointer to the data
2551: Level: intermediate
2553: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2554: @*/
2555: PetscErrorCode MatDenseRestoreArrayWriteAndMemType(Mat A, PetscScalar **array)
2556: {
2557: PetscBool isMPI;
2559: PetscFunctionBegin;
2561: PetscAssertPointer(array, 2);
2562: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2563: if (isMPI) {
2564: PetscCall(MatDenseRestoreArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array));
2565: } else {
2566: PetscErrorCode (*fptr)(Mat, PetscScalar **);
2568: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayWriteAndMemType_C", &fptr));
2569: if (fptr) {
2570: PetscCall((*fptr)(A, array));
2571: } else {
2572: PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2573: }
2574: *array = NULL;
2575: }
2576: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2577: PetscFunctionReturn(PETSC_SUCCESS);
2578: }
2580: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B)
2581: {
2582: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2583: PetscInt i, j, nrows, ncols, ldb;
2584: const PetscInt *irow, *icol;
2585: PetscScalar *av, *bv, *v = mat->v;
2586: Mat newmat;
2588: PetscFunctionBegin;
2589: PetscCall(ISGetIndices(isrow, &irow));
2590: PetscCall(ISGetIndices(iscol, &icol));
2591: PetscCall(ISGetLocalSize(isrow, &nrows));
2592: PetscCall(ISGetLocalSize(iscol, &ncols));
2594: /* Check submatrixcall */
2595: if (scall == MAT_REUSE_MATRIX) {
2596: PetscInt n_cols, n_rows;
2597: PetscCall(MatGetSize(*B, &n_rows, &n_cols));
2598: if (n_rows != nrows || n_cols != ncols) {
2599: /* resize the result matrix to match number of requested rows/columns */
2600: PetscCall(MatSetSizes(*B, nrows, ncols, nrows, ncols));
2601: }
2602: newmat = *B;
2603: } else {
2604: /* Create and fill new matrix */
2605: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &newmat));
2606: PetscCall(MatSetSizes(newmat, nrows, ncols, nrows, ncols));
2607: PetscCall(MatSetType(newmat, ((PetscObject)A)->type_name));
2608: PetscCall(MatSeqDenseSetPreallocation(newmat, NULL));
2609: }
2611: /* Now extract the data pointers and do the copy,column at a time */
2612: PetscCall(MatDenseGetArray(newmat, &bv));
2613: PetscCall(MatDenseGetLDA(newmat, &ldb));
2614: for (i = 0; i < ncols; i++) {
2615: av = v + mat->lda * icol[i];
2616: for (j = 0; j < nrows; j++) bv[j] = av[irow[j]];
2617: bv += ldb;
2618: }
2619: PetscCall(MatDenseRestoreArray(newmat, &bv));
2621: /* Assemble the matrices so that the correct flags are set */
2622: PetscCall(MatAssemblyBegin(newmat, MAT_FINAL_ASSEMBLY));
2623: PetscCall(MatAssemblyEnd(newmat, MAT_FINAL_ASSEMBLY));
2625: /* Free work space */
2626: PetscCall(ISRestoreIndices(isrow, &irow));
2627: PetscCall(ISRestoreIndices(iscol, &icol));
2628: *B = newmat;
2629: PetscFunctionReturn(PETSC_SUCCESS);
2630: }
2632: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[])
2633: {
2634: PetscInt i;
2636: PetscFunctionBegin;
2637: if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n, B));
2639: for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SeqDense(A, irow[i], icol[i], scall, &(*B)[i]));
2640: PetscFunctionReturn(PETSC_SUCCESS);
2641: }
2643: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat, MatAssemblyType mode)
2644: {
2645: PetscFunctionBegin;
2646: PetscFunctionReturn(PETSC_SUCCESS);
2647: }
2649: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat, MatAssemblyType mode)
2650: {
2651: PetscFunctionBegin;
2652: PetscFunctionReturn(PETSC_SUCCESS);
2653: }
2655: PetscErrorCode MatCopy_SeqDense(Mat A, Mat B, MatStructure str)
2656: {
2657: Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data;
2658: const PetscScalar *va;
2659: PetscScalar *vb;
2660: PetscInt lda1 = a->lda, lda2 = b->lda, m = A->rmap->n, n = A->cmap->n, j;
2662: PetscFunctionBegin;
2663: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2664: if (A->ops->copy != B->ops->copy) {
2665: PetscCall(MatCopy_Basic(A, B, str));
2666: PetscFunctionReturn(PETSC_SUCCESS);
2667: }
2668: PetscCheck(m == B->rmap->n && n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "size(B) != size(A)");
2669: PetscCall(MatDenseGetArrayRead(A, &va));
2670: PetscCall(MatDenseGetArray(B, &vb));
2671: if (lda1 > m || lda2 > m) {
2672: for (j = 0; j < n; j++) PetscCall(PetscArraycpy(vb + j * lda2, va + j * lda1, m));
2673: } else {
2674: PetscCall(PetscArraycpy(vb, va, A->rmap->n * A->cmap->n));
2675: }
2676: PetscCall(MatDenseRestoreArray(B, &vb));
2677: PetscCall(MatDenseRestoreArrayRead(A, &va));
2678: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2679: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2680: PetscFunctionReturn(PETSC_SUCCESS);
2681: }
2683: PetscErrorCode MatSetUp_SeqDense(Mat A)
2684: {
2685: PetscFunctionBegin;
2686: PetscCall(PetscLayoutSetUp(A->rmap));
2687: PetscCall(PetscLayoutSetUp(A->cmap));
2688: if (!A->preallocated) PetscCall(MatSeqDenseSetPreallocation(A, NULL));
2689: PetscFunctionReturn(PETSC_SUCCESS);
2690: }
2692: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2693: {
2694: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2695: PetscInt i, j;
2696: PetscInt min = PetscMin(A->rmap->n, A->cmap->n);
2697: PetscScalar *aa;
2699: PetscFunctionBegin;
2700: PetscCall(MatDenseGetArray(A, &aa));
2701: for (j = 0; j < A->cmap->n; j++) {
2702: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscConj(aa[i + j * mat->lda]);
2703: }
2704: PetscCall(MatDenseRestoreArray(A, &aa));
2705: if (mat->tau)
2706: for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
2707: PetscFunctionReturn(PETSC_SUCCESS);
2708: }
2710: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2711: {
2712: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2713: PetscInt i, j;
2714: PetscScalar *aa;
2716: PetscFunctionBegin;
2717: PetscCall(MatDenseGetArray(A, &aa));
2718: for (j = 0; j < A->cmap->n; j++) {
2719: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscRealPart(aa[i + j * mat->lda]);
2720: }
2721: PetscCall(MatDenseRestoreArray(A, &aa));
2722: PetscFunctionReturn(PETSC_SUCCESS);
2723: }
2725: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2726: {
2727: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2728: PetscInt i, j;
2729: PetscScalar *aa;
2731: PetscFunctionBegin;
2732: PetscCall(MatDenseGetArray(A, &aa));
2733: for (j = 0; j < A->cmap->n; j++) {
2734: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscImaginaryPart(aa[i + j * mat->lda]);
2735: }
2736: PetscCall(MatDenseRestoreArray(A, &aa));
2737: PetscFunctionReturn(PETSC_SUCCESS);
2738: }
2740: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2741: {
2742: PetscInt m = A->rmap->n, n = B->cmap->n;
2743: PetscBool cisdense = PETSC_FALSE;
2745: PetscFunctionBegin;
2746: PetscCall(MatSetSizes(C, m, n, m, n));
2747: #if defined(PETSC_HAVE_CUDA)
2748: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2749: #endif
2750: #if defined(PETSC_HAVE_HIP)
2751: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2752: #endif
2753: if (!cisdense) {
2754: PetscBool flg;
2756: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2757: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2758: }
2759: PetscCall(MatSetUp(C));
2760: PetscFunctionReturn(PETSC_SUCCESS);
2761: }
2763: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2764: {
2765: Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data, *c = (Mat_SeqDense *)C->data;
2766: PetscBLASInt m, n, k;
2767: const PetscScalar *av, *bv;
2768: PetscScalar *cv;
2769: PetscScalar _DOne = 1.0, _DZero = 0.0;
2771: PetscFunctionBegin;
2772: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2773: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2774: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2775: if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2776: PetscCall(MatDenseGetArrayRead(A, &av));
2777: PetscCall(MatDenseGetArrayRead(B, &bv));
2778: PetscCall(MatDenseGetArrayWrite(C, &cv));
2779: PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2780: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2781: PetscCall(MatDenseRestoreArrayRead(A, &av));
2782: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2783: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2784: PetscFunctionReturn(PETSC_SUCCESS);
2785: }
2787: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2788: {
2789: PetscInt m = A->rmap->n, n = B->rmap->n;
2790: PetscBool cisdense = PETSC_FALSE;
2792: PetscFunctionBegin;
2793: PetscCall(MatSetSizes(C, m, n, m, n));
2794: #if defined(PETSC_HAVE_CUDA)
2795: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2796: #endif
2797: #if defined(PETSC_HAVE_HIP)
2798: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2799: #endif
2800: if (!cisdense) {
2801: PetscBool flg;
2803: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2804: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2805: }
2806: PetscCall(MatSetUp(C));
2807: PetscFunctionReturn(PETSC_SUCCESS);
2808: }
2810: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2811: {
2812: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2813: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
2814: Mat_SeqDense *c = (Mat_SeqDense *)C->data;
2815: const PetscScalar *av, *bv;
2816: PetscScalar *cv;
2817: PetscBLASInt m, n, k;
2818: PetscScalar _DOne = 1.0, _DZero = 0.0;
2820: PetscFunctionBegin;
2821: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2822: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2823: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2824: if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2825: PetscCall(MatDenseGetArrayRead(A, &av));
2826: PetscCall(MatDenseGetArrayRead(B, &bv));
2827: PetscCall(MatDenseGetArrayWrite(C, &cv));
2828: PetscCallBLAS("BLASgemm", BLASgemm_("N", "T", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2829: PetscCall(MatDenseRestoreArrayRead(A, &av));
2830: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2831: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2832: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2833: PetscFunctionReturn(PETSC_SUCCESS);
2834: }
2836: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2837: {
2838: PetscInt m = A->cmap->n, n = B->cmap->n;
2839: PetscBool cisdense = PETSC_FALSE;
2841: PetscFunctionBegin;
2842: PetscCall(MatSetSizes(C, m, n, m, n));
2843: #if defined(PETSC_HAVE_CUDA)
2844: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2845: #endif
2846: #if defined(PETSC_HAVE_HIP)
2847: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2848: #endif
2849: if (!cisdense) {
2850: PetscBool flg;
2852: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2853: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2854: }
2855: PetscCall(MatSetUp(C));
2856: PetscFunctionReturn(PETSC_SUCCESS);
2857: }
2859: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2860: {
2861: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2862: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
2863: Mat_SeqDense *c = (Mat_SeqDense *)C->data;
2864: const PetscScalar *av, *bv;
2865: PetscScalar *cv;
2866: PetscBLASInt m, n, k;
2867: PetscScalar _DOne = 1.0, _DZero = 0.0;
2869: PetscFunctionBegin;
2870: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2871: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2872: PetscCall(PetscBLASIntCast(A->rmap->n, &k));
2873: if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2874: PetscCall(MatDenseGetArrayRead(A, &av));
2875: PetscCall(MatDenseGetArrayRead(B, &bv));
2876: PetscCall(MatDenseGetArrayWrite(C, &cv));
2877: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2878: PetscCall(MatDenseRestoreArrayRead(A, &av));
2879: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2880: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2881: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2882: PetscFunctionReturn(PETSC_SUCCESS);
2883: }
2885: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2886: {
2887: PetscFunctionBegin;
2888: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2889: C->ops->productsymbolic = MatProductSymbolic_AB;
2890: PetscFunctionReturn(PETSC_SUCCESS);
2891: }
2893: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2894: {
2895: PetscFunctionBegin;
2896: C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2897: C->ops->productsymbolic = MatProductSymbolic_AtB;
2898: PetscFunctionReturn(PETSC_SUCCESS);
2899: }
2901: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2902: {
2903: PetscFunctionBegin;
2904: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2905: C->ops->productsymbolic = MatProductSymbolic_ABt;
2906: PetscFunctionReturn(PETSC_SUCCESS);
2907: }
2909: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2910: {
2911: Mat_Product *product = C->product;
2913: PetscFunctionBegin;
2914: switch (product->type) {
2915: case MATPRODUCT_AB:
2916: PetscCall(MatProductSetFromOptions_SeqDense_AB(C));
2917: break;
2918: case MATPRODUCT_AtB:
2919: PetscCall(MatProductSetFromOptions_SeqDense_AtB(C));
2920: break;
2921: case MATPRODUCT_ABt:
2922: PetscCall(MatProductSetFromOptions_SeqDense_ABt(C));
2923: break;
2924: default:
2925: break;
2926: }
2927: PetscFunctionReturn(PETSC_SUCCESS);
2928: }
2930: static PetscErrorCode MatGetRowMax_SeqDense(Mat A, Vec v, PetscInt idx[])
2931: {
2932: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2933: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2934: PetscScalar *x;
2935: const PetscScalar *aa;
2937: PetscFunctionBegin;
2938: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2939: PetscCall(VecGetArray(v, &x));
2940: PetscCall(VecGetLocalSize(v, &p));
2941: PetscCall(MatDenseGetArrayRead(A, &aa));
2942: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2943: for (i = 0; i < m; i++) {
2944: x[i] = aa[i];
2945: if (idx) idx[i] = 0;
2946: for (j = 1; j < n; j++) {
2947: if (PetscRealPart(x[i]) < PetscRealPart(aa[i + a->lda * j])) {
2948: x[i] = aa[i + a->lda * j];
2949: if (idx) idx[i] = j;
2950: }
2951: }
2952: }
2953: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2954: PetscCall(VecRestoreArray(v, &x));
2955: PetscFunctionReturn(PETSC_SUCCESS);
2956: }
2958: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A, Vec v, PetscInt idx[])
2959: {
2960: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2961: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2962: PetscScalar *x;
2963: PetscReal atmp;
2964: const PetscScalar *aa;
2966: PetscFunctionBegin;
2967: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2968: PetscCall(VecGetArray(v, &x));
2969: PetscCall(VecGetLocalSize(v, &p));
2970: PetscCall(MatDenseGetArrayRead(A, &aa));
2971: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2972: for (i = 0; i < m; i++) {
2973: x[i] = PetscAbsScalar(aa[i]);
2974: for (j = 1; j < n; j++) {
2975: atmp = PetscAbsScalar(aa[i + a->lda * j]);
2976: if (PetscAbsScalar(x[i]) < atmp) {
2977: x[i] = atmp;
2978: if (idx) idx[i] = j;
2979: }
2980: }
2981: }
2982: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2983: PetscCall(VecRestoreArray(v, &x));
2984: PetscFunctionReturn(PETSC_SUCCESS);
2985: }
2987: static PetscErrorCode MatGetRowMin_SeqDense(Mat A, Vec v, PetscInt idx[])
2988: {
2989: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2990: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2991: PetscScalar *x;
2992: const PetscScalar *aa;
2994: PetscFunctionBegin;
2995: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2996: PetscCall(MatDenseGetArrayRead(A, &aa));
2997: PetscCall(VecGetArray(v, &x));
2998: PetscCall(VecGetLocalSize(v, &p));
2999: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
3000: for (i = 0; i < m; i++) {
3001: x[i] = aa[i];
3002: if (idx) idx[i] = 0;
3003: for (j = 1; j < n; j++) {
3004: if (PetscRealPart(x[i]) > PetscRealPart(aa[i + a->lda * j])) {
3005: x[i] = aa[i + a->lda * j];
3006: if (idx) idx[i] = j;
3007: }
3008: }
3009: }
3010: PetscCall(VecRestoreArray(v, &x));
3011: PetscCall(MatDenseRestoreArrayRead(A, &aa));
3012: PetscFunctionReturn(PETSC_SUCCESS);
3013: }
3015: PetscErrorCode MatGetColumnVector_SeqDense(Mat A, Vec v, PetscInt col)
3016: {
3017: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3018: PetscScalar *x;
3019: const PetscScalar *aa;
3021: PetscFunctionBegin;
3022: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3023: PetscCall(MatDenseGetArrayRead(A, &aa));
3024: PetscCall(VecGetArray(v, &x));
3025: PetscCall(PetscArraycpy(x, aa + col * a->lda, A->rmap->n));
3026: PetscCall(VecRestoreArray(v, &x));
3027: PetscCall(MatDenseRestoreArrayRead(A, &aa));
3028: PetscFunctionReturn(PETSC_SUCCESS);
3029: }
3031: PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A, PetscInt type, PetscReal *reductions)
3032: {
3033: PetscInt i, j, m, n;
3034: const PetscScalar *a;
3036: PetscFunctionBegin;
3037: PetscCall(MatGetSize(A, &m, &n));
3038: PetscCall(PetscArrayzero(reductions, n));
3039: PetscCall(MatDenseGetArrayRead(A, &a));
3040: if (type == NORM_2) {
3041: for (i = 0; i < n; i++) {
3042: for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j] * a[j]);
3043: a = PetscSafePointerPlusOffset(a, m);
3044: }
3045: } else if (type == NORM_1) {
3046: for (i = 0; i < n; i++) {
3047: for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j]);
3048: a = PetscSafePointerPlusOffset(a, m);
3049: }
3050: } else if (type == NORM_INFINITY) {
3051: for (i = 0; i < n; i++) {
3052: for (j = 0; j < m; j++) reductions[i] = PetscMax(PetscAbsScalar(a[j]), reductions[i]);
3053: a = PetscSafePointerPlusOffset(a, m);
3054: }
3055: } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
3056: for (i = 0; i < n; i++) {
3057: for (j = 0; j < m; j++) reductions[i] += PetscRealPart(a[j]);
3058: a = PetscSafePointerPlusOffset(a, m);
3059: }
3060: } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3061: for (i = 0; i < n; i++) {
3062: for (j = 0; j < m; j++) reductions[i] += PetscImaginaryPart(a[j]);
3063: a = PetscSafePointerPlusOffset(a, m);
3064: }
3065: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown reduction type");
3066: PetscCall(MatDenseRestoreArrayRead(A, &a));
3067: if (type == NORM_2) {
3068: for (i = 0; i < n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
3069: } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3070: for (i = 0; i < n; i++) reductions[i] /= m;
3071: }
3072: PetscFunctionReturn(PETSC_SUCCESS);
3073: }
3075: PetscErrorCode MatSetRandom_SeqDense(Mat x, PetscRandom rctx)
3076: {
3077: PetscScalar *a;
3078: PetscInt lda, m, n, i, j;
3080: PetscFunctionBegin;
3081: PetscCall(MatGetSize(x, &m, &n));
3082: PetscCall(MatDenseGetLDA(x, &lda));
3083: PetscCall(MatDenseGetArrayWrite(x, &a));
3084: for (j = 0; j < n; j++) {
3085: for (i = 0; i < m; i++) PetscCall(PetscRandomGetValue(rctx, a + j * lda + i));
3086: }
3087: PetscCall(MatDenseRestoreArrayWrite(x, &a));
3088: PetscFunctionReturn(PETSC_SUCCESS);
3089: }
3091: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A, PetscBool *missing, PetscInt *d)
3092: {
3093: PetscFunctionBegin;
3094: *missing = PETSC_FALSE;
3095: PetscFunctionReturn(PETSC_SUCCESS);
3096: }
3098: /* vals is not const */
3099: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A, PetscInt col, PetscScalar **vals)
3100: {
3101: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3102: PetscScalar *v;
3104: PetscFunctionBegin;
3105: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3106: PetscCall(MatDenseGetArray(A, &v));
3107: *vals = v + col * a->lda;
3108: PetscCall(MatDenseRestoreArray(A, &v));
3109: PetscFunctionReturn(PETSC_SUCCESS);
3110: }
3112: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A, PetscScalar **vals)
3113: {
3114: PetscFunctionBegin;
3115: if (vals) *vals = NULL; /* user cannot accidentally use the array later */
3116: PetscFunctionReturn(PETSC_SUCCESS);
3117: }
3119: static struct _MatOps MatOps_Values = {MatSetValues_SeqDense,
3120: MatGetRow_SeqDense,
3121: MatRestoreRow_SeqDense,
3122: MatMult_SeqDense,
3123: /* 4*/ MatMultAdd_SeqDense,
3124: MatMultTranspose_SeqDense,
3125: MatMultTransposeAdd_SeqDense,
3126: NULL,
3127: NULL,
3128: NULL,
3129: /* 10*/ NULL,
3130: MatLUFactor_SeqDense,
3131: MatCholeskyFactor_SeqDense,
3132: MatSOR_SeqDense,
3133: MatTranspose_SeqDense,
3134: /* 15*/ MatGetInfo_SeqDense,
3135: MatEqual_SeqDense,
3136: MatGetDiagonal_SeqDense,
3137: MatDiagonalScale_SeqDense,
3138: MatNorm_SeqDense,
3139: /* 20*/ MatAssemblyBegin_SeqDense,
3140: MatAssemblyEnd_SeqDense,
3141: MatSetOption_SeqDense,
3142: MatZeroEntries_SeqDense,
3143: /* 24*/ MatZeroRows_SeqDense,
3144: NULL,
3145: NULL,
3146: NULL,
3147: NULL,
3148: /* 29*/ MatSetUp_SeqDense,
3149: NULL,
3150: NULL,
3151: NULL,
3152: NULL,
3153: /* 34*/ MatDuplicate_SeqDense,
3154: NULL,
3155: NULL,
3156: NULL,
3157: NULL,
3158: /* 39*/ MatAXPY_SeqDense,
3159: MatCreateSubMatrices_SeqDense,
3160: NULL,
3161: MatGetValues_SeqDense,
3162: MatCopy_SeqDense,
3163: /* 44*/ MatGetRowMax_SeqDense,
3164: MatScale_SeqDense,
3165: MatShift_SeqDense,
3166: NULL,
3167: MatZeroRowsColumns_SeqDense,
3168: /* 49*/ MatSetRandom_SeqDense,
3169: NULL,
3170: NULL,
3171: NULL,
3172: NULL,
3173: /* 54*/ NULL,
3174: NULL,
3175: NULL,
3176: NULL,
3177: NULL,
3178: /* 59*/ MatCreateSubMatrix_SeqDense,
3179: MatDestroy_SeqDense,
3180: MatView_SeqDense,
3181: NULL,
3182: NULL,
3183: /* 64*/ NULL,
3184: NULL,
3185: NULL,
3186: NULL,
3187: NULL,
3188: /* 69*/ MatGetRowMaxAbs_SeqDense,
3189: NULL,
3190: NULL,
3191: NULL,
3192: NULL,
3193: /* 74*/ NULL,
3194: NULL,
3195: NULL,
3196: NULL,
3197: NULL,
3198: /* 79*/ NULL,
3199: NULL,
3200: NULL,
3201: NULL,
3202: /* 83*/ MatLoad_SeqDense,
3203: MatIsSymmetric_SeqDense,
3204: MatIsHermitian_SeqDense,
3205: NULL,
3206: NULL,
3207: NULL,
3208: /* 89*/ NULL,
3209: NULL,
3210: MatMatMultNumeric_SeqDense_SeqDense,
3211: NULL,
3212: NULL,
3213: /* 94*/ NULL,
3214: NULL,
3215: NULL,
3216: MatMatTransposeMultNumeric_SeqDense_SeqDense,
3217: NULL,
3218: /* 99*/ MatProductSetFromOptions_SeqDense,
3219: NULL,
3220: NULL,
3221: MatConjugate_SeqDense,
3222: NULL,
3223: /*104*/ NULL,
3224: MatRealPart_SeqDense,
3225: MatImaginaryPart_SeqDense,
3226: NULL,
3227: NULL,
3228: /*109*/ NULL,
3229: NULL,
3230: MatGetRowMin_SeqDense,
3231: MatGetColumnVector_SeqDense,
3232: MatMissingDiagonal_SeqDense,
3233: /*114*/ NULL,
3234: NULL,
3235: NULL,
3236: NULL,
3237: NULL,
3238: /*119*/ NULL,
3239: NULL,
3240: MatMultHermitianTranspose_SeqDense,
3241: MatMultHermitianTransposeAdd_SeqDense,
3242: NULL,
3243: /*124*/ NULL,
3244: MatGetColumnReductions_SeqDense,
3245: NULL,
3246: NULL,
3247: NULL,
3248: /*129*/ NULL,
3249: NULL,
3250: NULL,
3251: MatTransposeMatMultNumeric_SeqDense_SeqDense,
3252: NULL,
3253: /*134*/ NULL,
3254: NULL,
3255: NULL,
3256: NULL,
3257: NULL,
3258: /*139*/ NULL,
3259: NULL,
3260: NULL,
3261: NULL,
3262: NULL,
3263: MatCreateMPIMatConcatenateSeqMat_SeqDense,
3264: /*145*/ NULL,
3265: NULL,
3266: NULL,
3267: NULL,
3268: NULL,
3269: /*150*/ NULL,
3270: NULL,
3271: NULL};
3273: /*@C
3274: MatCreateSeqDense - Creates a `MATSEQDENSE` that
3275: is stored in column major order (the usual Fortran format).
3277: Collective
3279: Input Parameters:
3280: + comm - MPI communicator, set to `PETSC_COMM_SELF`
3281: . m - number of rows
3282: . n - number of columns
3283: - data - optional location of matrix data in column major order. Use `NULL` for PETSc
3284: to control all matrix memory allocation.
3286: Output Parameter:
3287: . A - the matrix
3289: Level: intermediate
3291: Note:
3292: The data input variable is intended primarily for Fortran programmers
3293: who wish to allocate their own matrix memory space. Most users should
3294: set `data` = `NULL`.
3296: Developer Note:
3297: Many of the matrix operations for this variant use the BLAS and LAPACK routines.
3299: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`
3300: @*/
3301: PetscErrorCode MatCreateSeqDense(MPI_Comm comm, PetscInt m, PetscInt n, PetscScalar *data, Mat *A)
3302: {
3303: PetscFunctionBegin;
3304: PetscCall(MatCreate(comm, A));
3305: PetscCall(MatSetSizes(*A, m, n, m, n));
3306: PetscCall(MatSetType(*A, MATSEQDENSE));
3307: PetscCall(MatSeqDenseSetPreallocation(*A, data));
3308: PetscFunctionReturn(PETSC_SUCCESS);
3309: }
3311: /*@C
3312: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements of a `MATSEQDENSE` matrix
3314: Collective
3316: Input Parameters:
3317: + B - the matrix
3318: - data - the array (or `NULL`)
3320: Level: intermediate
3322: Note:
3323: The data input variable is intended primarily for Fortran programmers
3324: who wish to allocate their own matrix memory space. Most users should
3325: need not call this routine.
3327: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`, `MatDenseSetLDA()`
3328: @*/
3329: PetscErrorCode MatSeqDenseSetPreallocation(Mat B, PetscScalar data[])
3330: {
3331: PetscFunctionBegin;
3333: PetscTryMethod(B, "MatSeqDenseSetPreallocation_C", (Mat, PetscScalar[]), (B, data));
3334: PetscFunctionReturn(PETSC_SUCCESS);
3335: }
3337: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B, PetscScalar *data)
3338: {
3339: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
3341: PetscFunctionBegin;
3342: PetscCheck(!b->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3343: B->preallocated = PETSC_TRUE;
3345: PetscCall(PetscLayoutSetUp(B->rmap));
3346: PetscCall(PetscLayoutSetUp(B->cmap));
3348: if (b->lda <= 0) b->lda = B->rmap->n;
3350: if (!data) { /* petsc-allocated storage */
3351: if (!b->user_alloc) PetscCall(PetscFree(b->v));
3352: PetscCall(PetscCalloc1((size_t)b->lda * B->cmap->n, &b->v));
3354: b->user_alloc = PETSC_FALSE;
3355: } else { /* user-allocated storage */
3356: if (!b->user_alloc) PetscCall(PetscFree(b->v));
3357: b->v = data;
3358: b->user_alloc = PETSC_TRUE;
3359: }
3360: B->assembled = PETSC_TRUE;
3361: PetscFunctionReturn(PETSC_SUCCESS);
3362: }
3364: #if defined(PETSC_HAVE_ELEMENTAL)
3365: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
3366: {
3367: Mat mat_elemental;
3368: const PetscScalar *array;
3369: PetscScalar *v_colwise;
3370: PetscInt M = A->rmap->N, N = A->cmap->N, i, j, k, *rows, *cols;
3372: PetscFunctionBegin;
3373: PetscCall(PetscMalloc3(M * N, &v_colwise, M, &rows, N, &cols));
3374: PetscCall(MatDenseGetArrayRead(A, &array));
3375: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
3376: k = 0;
3377: for (j = 0; j < N; j++) {
3378: cols[j] = j;
3379: for (i = 0; i < M; i++) v_colwise[j * M + i] = array[k++];
3380: }
3381: for (i = 0; i < M; i++) rows[i] = i;
3382: PetscCall(MatDenseRestoreArrayRead(A, &array));
3384: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
3385: PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
3386: PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
3387: PetscCall(MatSetUp(mat_elemental));
3389: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
3390: PetscCall(MatSetValues(mat_elemental, M, rows, N, cols, v_colwise, ADD_VALUES));
3391: PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
3392: PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
3393: PetscCall(PetscFree3(v_colwise, rows, cols));
3395: if (reuse == MAT_INPLACE_MATRIX) {
3396: PetscCall(MatHeaderReplace(A, &mat_elemental));
3397: } else {
3398: *newmat = mat_elemental;
3399: }
3400: PetscFunctionReturn(PETSC_SUCCESS);
3401: }
3402: #endif
3404: PetscErrorCode MatDenseSetLDA_SeqDense(Mat B, PetscInt lda)
3405: {
3406: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
3407: PetscBool data;
3409: PetscFunctionBegin;
3410: data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
3411: PetscCheck(b->user_alloc || !data || b->lda == lda, PETSC_COMM_SELF, PETSC_ERR_ORDER, "LDA cannot be changed after allocation of internal storage");
3412: 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);
3413: b->lda = lda;
3414: PetscFunctionReturn(PETSC_SUCCESS);
3415: }
3417: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
3418: {
3419: PetscFunctionBegin;
3420: PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIDense(comm, inmat, n, scall, outmat));
3421: PetscFunctionReturn(PETSC_SUCCESS);
3422: }
3424: PetscErrorCode MatDenseCreateColumnVec_Private(Mat A, Vec *v)
3425: {
3426: PetscBool isstd, iskok, iscuda, iship;
3427: PetscMPIInt size;
3428: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
3429: /* we pass the data of A, to prevent allocating needless GPU memory the first time VecCUPMPlaceArray is called. */
3430: const PetscScalar *a;
3431: #endif
3433: PetscFunctionBegin;
3434: *v = NULL;
3435: PetscCall(PetscStrcmpAny(A->defaultvectype, &isstd, VECSTANDARD, VECSEQ, VECMPI, ""));
3436: PetscCall(PetscStrcmpAny(A->defaultvectype, &iskok, VECKOKKOS, VECSEQKOKKOS, VECMPIKOKKOS, ""));
3437: PetscCall(PetscStrcmpAny(A->defaultvectype, &iscuda, VECCUDA, VECSEQCUDA, VECMPICUDA, ""));
3438: PetscCall(PetscStrcmpAny(A->defaultvectype, &iship, VECHIP, VECSEQHIP, VECMPIHIP, ""));
3439: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
3440: if (isstd) {
3441: if (size > 1) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, v));
3442: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, v));
3443: } else if (iskok) {
3444: PetscCheck(PetscDefined(HAVE_KOKKOS_KERNELS), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using KOKKOS kernels support");
3445: #if PetscDefined(HAVE_KOKKOS_KERNELS)
3446: if (size > 1) PetscCall(VecCreateMPIKokkosWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, v));
3447: else PetscCall(VecCreateSeqKokkosWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, v));
3448: #endif
3449: } else if (iscuda) {
3450: PetscCheck(PetscDefined(HAVE_CUDA), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using CUDA support");
3451: #if PetscDefined(HAVE_CUDA)
3452: PetscCall(MatDenseCUDAGetArrayRead(A, &a));
3453: if (size > 1) PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, a, v));
3454: else PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, a, v));
3455: #endif
3456: } else if (iship) {
3457: PetscCheck(PetscDefined(HAVE_HIP), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using HIP support");
3458: #if PetscDefined(HAVE_HIP)
3459: PetscCall(MatDenseHIPGetArrayRead(A, &a));
3460: if (size > 1) PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, a, v));
3461: else PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, a, v));
3462: #endif
3463: }
3464: PetscCheck(*v, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not coded for type %s", A->defaultvectype);
3465: PetscFunctionReturn(PETSC_SUCCESS);
3466: }
3468: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3469: {
3470: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3472: PetscFunctionBegin;
3473: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3474: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3475: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3476: a->vecinuse = col + 1;
3477: PetscCall(MatDenseGetArray(A, (PetscScalar **)&a->ptrinuse));
3478: PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda));
3479: *v = a->cvec;
3480: PetscFunctionReturn(PETSC_SUCCESS);
3481: }
3483: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3484: {
3485: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3487: PetscFunctionBegin;
3488: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3489: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3490: a->vecinuse = 0;
3491: PetscCall(MatDenseRestoreArray(A, (PetscScalar **)&a->ptrinuse));
3492: PetscCall(VecResetArray(a->cvec));
3493: if (v) *v = NULL;
3494: PetscFunctionReturn(PETSC_SUCCESS);
3495: }
3497: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3498: {
3499: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3501: PetscFunctionBegin;
3502: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3503: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3504: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3505: a->vecinuse = col + 1;
3506: PetscCall(MatDenseGetArrayRead(A, &a->ptrinuse));
3507: PetscCall(VecPlaceArray(a->cvec, PetscSafePointerPlusOffset(a->ptrinuse, (size_t)col * (size_t)a->lda)));
3508: PetscCall(VecLockReadPush(a->cvec));
3509: *v = a->cvec;
3510: PetscFunctionReturn(PETSC_SUCCESS);
3511: }
3513: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3514: {
3515: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3517: PetscFunctionBegin;
3518: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3519: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3520: a->vecinuse = 0;
3521: PetscCall(MatDenseRestoreArrayRead(A, &a->ptrinuse));
3522: PetscCall(VecLockReadPop(a->cvec));
3523: PetscCall(VecResetArray(a->cvec));
3524: if (v) *v = NULL;
3525: PetscFunctionReturn(PETSC_SUCCESS);
3526: }
3528: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3529: {
3530: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3532: PetscFunctionBegin;
3533: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3534: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3535: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3536: a->vecinuse = col + 1;
3537: PetscCall(MatDenseGetArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3538: PetscCall(VecPlaceArray(a->cvec, PetscSafePointerPlusOffset(a->ptrinuse, (size_t)col * (size_t)a->lda)));
3539: *v = a->cvec;
3540: PetscFunctionReturn(PETSC_SUCCESS);
3541: }
3543: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3544: {
3545: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3547: PetscFunctionBegin;
3548: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3549: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3550: a->vecinuse = 0;
3551: PetscCall(MatDenseRestoreArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3552: PetscCall(VecResetArray(a->cvec));
3553: if (v) *v = NULL;
3554: PetscFunctionReturn(PETSC_SUCCESS);
3555: }
3557: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3558: {
3559: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3561: PetscFunctionBegin;
3562: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3563: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3564: if (a->cmat && (cend - cbegin != a->cmat->cmap->N || rend - rbegin != a->cmat->rmap->N)) PetscCall(MatDestroy(&a->cmat));
3565: if (!a->cmat) {
3566: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), rend - rbegin, PETSC_DECIDE, rend - rbegin, cend - cbegin, PetscSafePointerPlusOffset(a->v, rbegin + (size_t)cbegin * a->lda), &a->cmat));
3567: } else {
3568: PetscCall(MatDensePlaceArray(a->cmat, PetscSafePointerPlusOffset(a->v, rbegin + (size_t)cbegin * a->lda)));
3569: }
3570: PetscCall(MatDenseSetLDA(a->cmat, a->lda));
3571: a->matinuse = cbegin + 1;
3572: *v = a->cmat;
3573: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
3574: A->offloadmask = PETSC_OFFLOAD_CPU;
3575: #endif
3576: PetscFunctionReturn(PETSC_SUCCESS);
3577: }
3579: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A, Mat *v)
3580: {
3581: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3583: PetscFunctionBegin;
3584: PetscCheck(a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetSubMatrix() first");
3585: PetscCheck(a->cmat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column matrix");
3586: PetscCheck(*v == a->cmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not the matrix obtained from MatDenseGetSubMatrix()");
3587: a->matinuse = 0;
3588: PetscCall(MatDenseResetArray(a->cmat));
3589: if (v) *v = NULL;
3590: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
3591: A->offloadmask = PETSC_OFFLOAD_CPU;
3592: #endif
3593: PetscFunctionReturn(PETSC_SUCCESS);
3594: }
3596: /*MC
3597: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
3599: Options Database Key:
3600: . -mat_type seqdense - sets the matrix type to `MATSEQDENSE` during a call to `MatSetFromOptions()`
3602: Level: beginner
3604: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreateSeqDense()`
3605: M*/
3606: PetscErrorCode MatCreate_SeqDense(Mat B)
3607: {
3608: Mat_SeqDense *b;
3609: PetscMPIInt size;
3611: PetscFunctionBegin;
3612: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
3613: PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1");
3615: PetscCall(PetscNew(&b));
3616: B->data = (void *)b;
3617: B->ops[0] = MatOps_Values;
3619: b->roworiented = PETSC_TRUE;
3621: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactor_C", MatQRFactor_SeqDense));
3622: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetLDA_C", MatDenseGetLDA_SeqDense));
3623: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseSetLDA_C", MatDenseSetLDA_SeqDense));
3624: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArray_C", MatDenseGetArray_SeqDense));
3625: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArray_C", MatDenseRestoreArray_SeqDense));
3626: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDensePlaceArray_C", MatDensePlaceArray_SeqDense));
3627: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseResetArray_C", MatDenseResetArray_SeqDense));
3628: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseReplaceArray_C", MatDenseReplaceArray_SeqDense));
3629: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayRead_C", MatDenseGetArray_SeqDense));
3630: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayRead_C", MatDenseRestoreArray_SeqDense));
3631: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayWrite_C", MatDenseGetArray_SeqDense));
3632: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayWrite_C", MatDenseRestoreArray_SeqDense));
3633: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqaij_C", MatConvert_SeqDense_SeqAIJ));
3634: #if defined(PETSC_HAVE_ELEMENTAL)
3635: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_elemental_C", MatConvert_SeqDense_Elemental));
3636: #endif
3637: #if defined(PETSC_HAVE_SCALAPACK)
3638: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_scalapack_C", MatConvert_Dense_ScaLAPACK));
3639: #endif
3640: #if defined(PETSC_HAVE_CUDA)
3641: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensecuda_C", MatConvert_SeqDense_SeqDenseCUDA));
3642: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3643: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdense_C", MatProductSetFromOptions_SeqDense));
3644: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3645: #endif
3646: #if defined(PETSC_HAVE_HIP)
3647: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensehip_C", MatConvert_SeqDense_SeqDenseHIP));
3648: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3649: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdense_C", MatProductSetFromOptions_SeqDense));
3650: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3651: #endif
3652: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqDenseSetPreallocation_C", MatSeqDenseSetPreallocation_SeqDense));
3653: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqdense_C", MatProductSetFromOptions_SeqAIJ_SeqDense));
3654: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdense_C", MatProductSetFromOptions_SeqDense));
3655: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));
3656: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqsbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));
3658: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumn_C", MatDenseGetColumn_SeqDense));
3659: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumn_C", MatDenseRestoreColumn_SeqDense));
3660: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVec_C", MatDenseGetColumnVec_SeqDense));
3661: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVec_C", MatDenseRestoreColumnVec_SeqDense));
3662: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecRead_C", MatDenseGetColumnVecRead_SeqDense));
3663: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecRead_C", MatDenseRestoreColumnVecRead_SeqDense));
3664: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecWrite_C", MatDenseGetColumnVecWrite_SeqDense));
3665: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecWrite_C", MatDenseRestoreColumnVecWrite_SeqDense));
3666: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetSubMatrix_C", MatDenseGetSubMatrix_SeqDense));
3667: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreSubMatrix_C", MatDenseRestoreSubMatrix_SeqDense));
3668: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultAddColumnRange_C", MatMultAddColumnRange_SeqDense));
3669: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultHermitianTransposeColumnRange_C", MatMultHermitianTransposeColumnRange_SeqDense));
3670: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatMultHermitianTransposeAddColumnRange_C", MatMultHermitianTransposeAddColumnRange_SeqDense));
3671: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQDENSE));
3672: PetscFunctionReturn(PETSC_SUCCESS);
3673: }
3675: /*@C
3676: 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.
3678: Not Collective
3680: Input Parameters:
3681: + A - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3682: - col - column index
3684: Output Parameter:
3685: . vals - pointer to the data
3687: Level: intermediate
3689: Note:
3690: Use `MatDenseGetColumnVec()` to get access to a column of a `MATDENSE` treated as a `Vec`
3692: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreColumn()`, `MatDenseGetColumnVec()`
3693: @*/
3694: PetscErrorCode MatDenseGetColumn(Mat A, PetscInt col, PetscScalar **vals)
3695: {
3696: PetscFunctionBegin;
3699: PetscAssertPointer(vals, 3);
3700: PetscUseMethod(A, "MatDenseGetColumn_C", (Mat, PetscInt, PetscScalar **), (A, col, vals));
3701: PetscFunctionReturn(PETSC_SUCCESS);
3702: }
3704: /*@C
3705: MatDenseRestoreColumn - returns access to a column of a `MATDENSE` matrix which is returned by `MatDenseGetColumn()`.
3707: Not Collective
3709: Input Parameters:
3710: + A - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3711: - vals - pointer to the data (may be `NULL`)
3713: Level: intermediate
3715: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetColumn()`
3716: @*/
3717: PetscErrorCode MatDenseRestoreColumn(Mat A, PetscScalar **vals)
3718: {
3719: PetscFunctionBegin;
3721: PetscAssertPointer(vals, 2);
3722: PetscUseMethod(A, "MatDenseRestoreColumn_C", (Mat, PetscScalar **), (A, vals));
3723: PetscFunctionReturn(PETSC_SUCCESS);
3724: }
3726: /*@
3727: MatDenseGetColumnVec - Gives read-write access to a column of a `MATDENSE` matrix, represented as a `Vec`.
3729: Collective
3731: Input Parameters:
3732: + A - the `Mat` object
3733: - col - the column index
3735: Output Parameter:
3736: . v - the vector
3738: Level: intermediate
3740: Notes:
3741: The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVec()` when the vector is no longer needed.
3743: Use `MatDenseGetColumnVecRead()` to obtain read-only access or `MatDenseGetColumnVecWrite()` for write-only access.
3745: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`, `MatDenseGetColumn()`
3746: @*/
3747: PetscErrorCode MatDenseGetColumnVec(Mat A, PetscInt col, Vec *v)
3748: {
3749: PetscFunctionBegin;
3753: PetscAssertPointer(v, 3);
3754: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3755: 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);
3756: PetscUseMethod(A, "MatDenseGetColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3757: PetscFunctionReturn(PETSC_SUCCESS);
3758: }
3760: /*@
3761: MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVec()`.
3763: Collective
3765: Input Parameters:
3766: + A - the `Mat` object
3767: . col - the column index
3768: - v - the `Vec` object (may be `NULL`)
3770: Level: intermediate
3772: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3773: @*/
3774: PetscErrorCode MatDenseRestoreColumnVec(Mat A, PetscInt col, Vec *v)
3775: {
3776: PetscFunctionBegin;
3780: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3781: 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);
3782: PetscUseMethod(A, "MatDenseRestoreColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3783: PetscFunctionReturn(PETSC_SUCCESS);
3784: }
3786: /*@
3787: MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a `Vec`.
3789: Collective
3791: Input Parameters:
3792: + A - the `Mat` object
3793: - col - the column index
3795: Output Parameter:
3796: . v - the vector
3798: Level: intermediate
3800: Notes:
3801: The vector is owned by PETSc and users cannot modify it.
3803: Users need to call `MatDenseRestoreColumnVecRead()` when the vector is no longer needed.
3805: Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecWrite()` for write-only access.
3807: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3808: @*/
3809: PetscErrorCode MatDenseGetColumnVecRead(Mat A, PetscInt col, Vec *v)
3810: {
3811: PetscFunctionBegin;
3815: PetscAssertPointer(v, 3);
3816: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3817: PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N);
3818: PetscUseMethod(A, "MatDenseGetColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3819: PetscFunctionReturn(PETSC_SUCCESS);
3820: }
3822: /*@
3823: MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVecRead()`.
3825: Collective
3827: Input Parameters:
3828: + A - the `Mat` object
3829: . col - the column index
3830: - v - the `Vec` object (may be `NULL`)
3832: Level: intermediate
3834: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecWrite()`
3835: @*/
3836: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A, PetscInt col, Vec *v)
3837: {
3838: PetscFunctionBegin;
3842: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3843: 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);
3844: PetscUseMethod(A, "MatDenseRestoreColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3845: PetscFunctionReturn(PETSC_SUCCESS);
3846: }
3848: /*@
3849: MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a `Vec`.
3851: Collective
3853: Input Parameters:
3854: + A - the `Mat` object
3855: - col - the column index
3857: Output Parameter:
3858: . v - the vector
3860: Level: intermediate
3862: Notes:
3863: The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVecWrite()` when the vector is no longer needed.
3865: Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecRead()` for read-only access.
3867: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3868: @*/
3869: PetscErrorCode MatDenseGetColumnVecWrite(Mat A, PetscInt col, Vec *v)
3870: {
3871: PetscFunctionBegin;
3875: PetscAssertPointer(v, 3);
3876: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3877: 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);
3878: PetscUseMethod(A, "MatDenseGetColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3879: PetscFunctionReturn(PETSC_SUCCESS);
3880: }
3882: /*@
3883: MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVecWrite()`.
3885: Collective
3887: Input Parameters:
3888: + A - the `Mat` object
3889: . col - the column index
3890: - v - the `Vec` object (may be `NULL`)
3892: Level: intermediate
3894: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`
3895: @*/
3896: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A, PetscInt col, Vec *v)
3897: {
3898: PetscFunctionBegin;
3902: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3903: 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);
3904: PetscUseMethod(A, "MatDenseRestoreColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3905: PetscFunctionReturn(PETSC_SUCCESS);
3906: }
3908: /*@
3909: MatDenseGetSubMatrix - Gives access to a block of rows and columns of a dense matrix, represented as a `Mat`.
3911: Collective
3913: Input Parameters:
3914: + A - the `Mat` object
3915: . rbegin - the first global row index in the block (if `PETSC_DECIDE`, is 0)
3916: . rend - the global row index past the last one in the block (if `PETSC_DECIDE`, is `M`)
3917: . cbegin - the first global column index in the block (if `PETSC_DECIDE`, is 0)
3918: - cend - the global column index past the last one in the block (if `PETSC_DECIDE`, is `N`)
3920: Output Parameter:
3921: . v - the matrix
3923: Level: intermediate
3925: Notes:
3926: The matrix is owned by PETSc. Users need to call `MatDenseRestoreSubMatrix()` when the matrix is no longer needed.
3928: The output matrix is not redistributed by PETSc, so depending on the values of `rbegin` and `rend`, some processes may have no local rows.
3930: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreSubMatrix()`
3931: @*/
3932: PetscErrorCode MatDenseGetSubMatrix(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3933: {
3934: PetscFunctionBegin;
3941: PetscAssertPointer(v, 6);
3942: if (rbegin == PETSC_DECIDE) rbegin = 0;
3943: if (rend == PETSC_DECIDE) rend = A->rmap->N;
3944: if (cbegin == PETSC_DECIDE) cbegin = 0;
3945: if (cend == PETSC_DECIDE) cend = A->cmap->N;
3946: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3947: 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);
3948: 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);
3949: 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);
3950: 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);
3951: PetscUseMethod(A, "MatDenseGetSubMatrix_C", (Mat, PetscInt, PetscInt, PetscInt, PetscInt, Mat *), (A, rbegin, rend, cbegin, cend, v));
3952: PetscFunctionReturn(PETSC_SUCCESS);
3953: }
3955: /*@
3956: MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from `MatDenseGetSubMatrix()`.
3958: Collective
3960: Input Parameters:
3961: + A - the `Mat` object
3962: - v - the `Mat` object (may be `NULL`)
3964: Level: intermediate
3966: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseGetSubMatrix()`
3967: @*/
3968: PetscErrorCode MatDenseRestoreSubMatrix(Mat A, Mat *v)
3969: {
3970: PetscFunctionBegin;
3973: PetscAssertPointer(v, 2);
3974: PetscUseMethod(A, "MatDenseRestoreSubMatrix_C", (Mat, Mat *), (A, v));
3975: PetscFunctionReturn(PETSC_SUCCESS);
3976: }
3978: #include <petscblaslapack.h>
3979: #include <petsc/private/kernels/blockinvert.h>
3981: PetscErrorCode MatSeqDenseInvert(Mat A)
3982: {
3983: PetscInt m;
3984: const PetscReal shift = 0.0;
3985: PetscBool allowzeropivot, zeropivotdetected = PETSC_FALSE;
3986: PetscScalar *values;
3988: PetscFunctionBegin;
3990: PetscCall(MatDenseGetArray(A, &values));
3991: PetscCall(MatGetLocalSize(A, &m, NULL));
3992: allowzeropivot = PetscNot(A->erroriffailure);
3993: /* factor and invert each block */
3994: switch (m) {
3995: case 1:
3996: values[0] = (PetscScalar)1.0 / (values[0] + shift);
3997: break;
3998: case 2:
3999: PetscCall(PetscKernel_A_gets_inverse_A_2(values, shift, allowzeropivot, &zeropivotdetected));
4000: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4001: break;
4002: case 3:
4003: PetscCall(PetscKernel_A_gets_inverse_A_3(values, shift, allowzeropivot, &zeropivotdetected));
4004: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4005: break;
4006: case 4:
4007: PetscCall(PetscKernel_A_gets_inverse_A_4(values, shift, allowzeropivot, &zeropivotdetected));
4008: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4009: break;
4010: case 5: {
4011: PetscScalar work[25];
4012: PetscInt ipvt[5];
4014: PetscCall(PetscKernel_A_gets_inverse_A_5(values, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
4015: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4016: } break;
4017: case 6:
4018: PetscCall(PetscKernel_A_gets_inverse_A_6(values, shift, allowzeropivot, &zeropivotdetected));
4019: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4020: break;
4021: case 7:
4022: PetscCall(PetscKernel_A_gets_inverse_A_7(values, shift, allowzeropivot, &zeropivotdetected));
4023: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4024: break;
4025: default: {
4026: PetscInt *v_pivots, *IJ, j;
4027: PetscScalar *v_work;
4029: PetscCall(PetscMalloc3(m, &v_work, m, &v_pivots, m, &IJ));
4030: for (j = 0; j < m; j++) IJ[j] = j;
4031: PetscCall(PetscKernel_A_gets_inverse_A(m, values, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
4032: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
4033: PetscCall(PetscFree3(v_work, v_pivots, IJ));
4034: }
4035: }
4036: PetscCall(MatDenseRestoreArray(A, &values));
4037: PetscFunctionReturn(PETSC_SUCCESS);
4038: }