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: PetscInt j;
254: for (j = 0; j < X->cmap->n; j++) PetscCallBLAS("BLASaxpy", BLASaxpy_(&m, &alpha, xv + j * ldax, &one, yv + j * lday, &one));
255: } else {
256: PetscCallBLAS("BLASaxpy", BLASaxpy_(&N, &alpha, xv, &one, yv, &one));
257: }
258: PetscCall(MatDenseRestoreArrayRead(X, &xv));
259: PetscCall(MatDenseRestoreArray(Y, &yv));
260: PetscCall(PetscLogFlops(PetscMax(2.0 * N - 1, 0)));
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: static PetscErrorCode MatGetInfo_SeqDense(Mat A, MatInfoType flag, MatInfo *info)
265: {
266: PetscLogDouble N = A->rmap->n * A->cmap->n;
268: PetscFunctionBegin;
269: info->block_size = 1.0;
270: info->nz_allocated = N;
271: info->nz_used = N;
272: info->nz_unneeded = 0;
273: info->assemblies = A->num_ass;
274: info->mallocs = 0;
275: info->memory = 0; /* REVIEW ME */
276: info->fill_ratio_given = 0;
277: info->fill_ratio_needed = 0;
278: info->factor_mallocs = 0;
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: PetscErrorCode MatScale_SeqDense(Mat A, PetscScalar alpha)
283: {
284: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
285: PetscScalar *v;
286: PetscBLASInt one = 1, j, nz, lda = 0;
288: PetscFunctionBegin;
289: PetscCall(MatDenseGetArray(A, &v));
290: PetscCall(PetscBLASIntCast(a->lda, &lda));
291: if (lda > A->rmap->n) {
292: PetscCall(PetscBLASIntCast(A->rmap->n, &nz));
293: for (j = 0; j < A->cmap->n; j++) PetscCallBLAS("BLASscal", BLASscal_(&nz, &alpha, v + j * lda, &one));
294: } else {
295: PetscCall(PetscBLASIntCast(A->rmap->n * A->cmap->n, &nz));
296: PetscCallBLAS("BLASscal", BLASscal_(&nz, &alpha, v, &one));
297: }
298: PetscCall(PetscLogFlops(A->rmap->n * A->cmap->n));
299: PetscCall(MatDenseRestoreArray(A, &v));
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: PetscErrorCode MatShift_SeqDense(Mat A, PetscScalar alpha)
304: {
305: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
306: PetscScalar *v;
307: PetscInt j, k;
309: PetscFunctionBegin;
310: PetscCall(MatDenseGetArray(A, &v));
311: k = PetscMin(A->rmap->n, A->cmap->n);
312: for (j = 0; j < k; j++) v[j + j * a->lda] += alpha;
313: PetscCall(PetscLogFlops(k));
314: PetscCall(MatDenseRestoreArray(A, &v));
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: static PetscErrorCode MatIsHermitian_SeqDense(Mat A, PetscReal rtol, PetscBool *fl)
319: {
320: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
321: PetscInt i, j, m = A->rmap->n, N = a->lda;
322: const PetscScalar *v;
324: PetscFunctionBegin;
325: *fl = PETSC_FALSE;
326: if (A->rmap->n != A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
327: PetscCall(MatDenseGetArrayRead(A, &v));
328: for (i = 0; i < m; i++) {
329: for (j = i; j < m; j++) {
330: if (PetscAbsScalar(v[i + j * N] - PetscConj(v[j + i * N])) > rtol) goto restore;
331: }
332: }
333: *fl = PETSC_TRUE;
334: restore:
335: PetscCall(MatDenseRestoreArrayRead(A, &v));
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: static PetscErrorCode MatIsSymmetric_SeqDense(Mat A, PetscReal rtol, PetscBool *fl)
340: {
341: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
342: PetscInt i, j, m = A->rmap->n, N = a->lda;
343: const PetscScalar *v;
345: PetscFunctionBegin;
346: *fl = PETSC_FALSE;
347: if (A->rmap->n != A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
348: PetscCall(MatDenseGetArrayRead(A, &v));
349: for (i = 0; i < m; i++) {
350: for (j = i; j < m; j++) {
351: if (PetscAbsScalar(v[i + j * N] - v[j + i * N]) > rtol) goto restore;
352: }
353: }
354: *fl = PETSC_TRUE;
355: restore:
356: PetscCall(MatDenseRestoreArrayRead(A, &v));
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi, Mat A, MatDuplicateOption cpvalues)
361: {
362: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
363: PetscInt lda = (PetscInt)mat->lda, j, m, nlda = lda;
364: PetscBool isdensecpu;
366: PetscFunctionBegin;
367: PetscCall(PetscLayoutReference(A->rmap, &newi->rmap));
368: PetscCall(PetscLayoutReference(A->cmap, &newi->cmap));
369: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */
370: PetscCall(MatDenseSetLDA(newi, lda));
371: }
372: PetscCall(PetscObjectTypeCompare((PetscObject)newi, MATSEQDENSE, &isdensecpu));
373: if (isdensecpu) PetscCall(MatSeqDenseSetPreallocation(newi, NULL));
374: if (cpvalues == MAT_COPY_VALUES) {
375: const PetscScalar *av;
376: PetscScalar *v;
378: PetscCall(MatDenseGetArrayRead(A, &av));
379: PetscCall(MatDenseGetArrayWrite(newi, &v));
380: PetscCall(MatDenseGetLDA(newi, &nlda));
381: m = A->rmap->n;
382: if (lda > m || nlda > m) {
383: for (j = 0; j < A->cmap->n; j++) PetscCall(PetscArraycpy(v + j * nlda, av + j * lda, m));
384: } else {
385: PetscCall(PetscArraycpy(v, av, A->rmap->n * A->cmap->n));
386: }
387: PetscCall(MatDenseRestoreArrayWrite(newi, &v));
388: PetscCall(MatDenseRestoreArrayRead(A, &av));
389: }
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: PetscErrorCode MatDuplicate_SeqDense(Mat A, MatDuplicateOption cpvalues, Mat *newmat)
394: {
395: PetscFunctionBegin;
396: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), newmat));
397: PetscCall(MatSetSizes(*newmat, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
398: PetscCall(MatSetType(*newmat, ((PetscObject)A)->type_name));
399: PetscCall(MatDuplicateNoCreate_SeqDense(*newmat, A, cpvalues));
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: static PetscErrorCode MatSolve_SeqDense_Internal_LU(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
404: {
405: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
406: PetscBLASInt info;
408: PetscFunctionBegin;
409: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
410: PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_(T ? "T" : "N", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
411: PetscCall(PetscFPTrapPop());
412: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "GETRS - Bad solve %d", (int)info);
413: PetscCall(PetscLogFlops(nrhs * (2.0 * m * m - m)));
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: static PetscErrorCode MatConjugate_SeqDense(Mat);
419: static PetscErrorCode MatSolve_SeqDense_Internal_Cholesky(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k, PetscBool T)
420: {
421: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
422: PetscBLASInt info;
424: PetscFunctionBegin;
425: if (A->spd == PETSC_BOOL3_TRUE) {
426: if (PetscDefined(USE_COMPLEX) && T) PetscCall(MatConjugate_SeqDense(A));
427: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
428: PetscCallBLAS("LAPACKpotrs", LAPACKpotrs_("L", &m, &nrhs, mat->v, &mat->lda, x, &m, &info));
429: PetscCall(PetscFPTrapPop());
430: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "POTRS Bad solve %d", (int)info);
431: if (PetscDefined(USE_COMPLEX) && T) PetscCall(MatConjugate_SeqDense(A));
432: #if defined(PETSC_USE_COMPLEX)
433: } else if (A->hermitian == PETSC_BOOL3_TRUE) {
434: if (T) PetscCall(MatConjugate_SeqDense(A));
435: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
436: PetscCallBLAS("LAPACKhetrs", LAPACKhetrs_("L", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
437: PetscCall(PetscFPTrapPop());
438: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "HETRS Bad solve %d", (int)info);
439: if (T) PetscCall(MatConjugate_SeqDense(A));
440: #endif
441: } else { /* symmetric case */
442: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
443: PetscCallBLAS("LAPACKsytrs", LAPACKsytrs_("L", &m, &nrhs, mat->v, &mat->lda, mat->pivots, x, &m, &info));
444: PetscCall(PetscFPTrapPop());
445: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "SYTRS Bad solve %d", (int)info);
446: }
447: PetscCall(PetscLogFlops(nrhs * (2.0 * m * m - m)));
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: static PetscErrorCode MatSolve_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
452: {
453: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
454: PetscBLASInt info;
455: char trans;
457: PetscFunctionBegin;
458: if (PetscDefined(USE_COMPLEX)) {
459: trans = 'C';
460: } else {
461: trans = 'T';
462: }
463: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
464: { /* lwork depends on the number of right-hand sides */
465: PetscBLASInt nlfwork, lfwork = -1;
466: PetscScalar fwork;
468: PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", &trans, &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, &fwork, &lfwork, &info));
469: nlfwork = (PetscBLASInt)PetscRealPart(fwork);
470: if (nlfwork > mat->lfwork) {
471: mat->lfwork = nlfwork;
472: PetscCall(PetscFree(mat->fwork));
473: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
474: }
475: }
476: PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", &trans, &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, mat->fwork, &mat->lfwork, &info));
477: PetscCall(PetscFPTrapPop());
478: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "ORMQR - Bad orthogonal transform %d", (int)info);
479: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
480: PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "N", "N", &mat->rank, &nrhs, mat->v, &mat->lda, x, &ldx, &info));
481: PetscCall(PetscFPTrapPop());
482: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "TRTRS - Bad triangular solve %d", (int)info);
483: for (PetscInt j = 0; j < nrhs; j++) {
484: for (PetscInt i = mat->rank; i < k; i++) x[j * ldx + i] = 0.;
485: }
486: PetscCall(PetscLogFlops(nrhs * (4.0 * m * mat->rank - PetscSqr(mat->rank))));
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: static PetscErrorCode MatSolveTranspose_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt ldx, PetscBLASInt m, PetscBLASInt nrhs, PetscBLASInt k)
491: {
492: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
493: PetscBLASInt info;
495: PetscFunctionBegin;
496: if (A->rmap->n == A->cmap->n && mat->rank == A->rmap->n) {
497: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
498: PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "T", "N", &m, &nrhs, mat->v, &mat->lda, x, &ldx, &info));
499: PetscCall(PetscFPTrapPop());
500: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "TRTRS - Bad triangular solve %d", (int)info);
501: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate_SeqDense(A));
502: { /* lwork depends on the number of right-hand sides */
503: PetscBLASInt nlfwork, lfwork = -1;
504: PetscScalar fwork;
506: PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", "N", &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, &fwork, &lfwork, &info));
507: nlfwork = (PetscBLASInt)PetscRealPart(fwork);
508: if (nlfwork > mat->lfwork) {
509: mat->lfwork = nlfwork;
510: PetscCall(PetscFree(mat->fwork));
511: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
512: }
513: }
514: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
515: PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", "N", &m, &nrhs, &mat->rank, mat->v, &mat->lda, mat->tau, x, &ldx, mat->fwork, &mat->lfwork, &info));
516: PetscCall(PetscFPTrapPop());
517: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "ORMQR - Bad orthogonal transform %d", (int)info);
518: if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate_SeqDense(A));
519: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "QR factored matrix cannot be used for transpose solve");
520: PetscCall(PetscLogFlops(nrhs * (4.0 * m * mat->rank - PetscSqr(mat->rank))));
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: static PetscErrorCode MatSolve_SeqDense_SetUp(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
525: {
526: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
527: PetscScalar *y;
528: PetscBLASInt m = 0, k = 0;
530: PetscFunctionBegin;
531: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
532: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
533: if (k < m) {
534: PetscCall(VecCopy(xx, mat->qrrhs));
535: PetscCall(VecGetArray(mat->qrrhs, &y));
536: } else {
537: PetscCall(VecCopy(xx, yy));
538: PetscCall(VecGetArray(yy, &y));
539: }
540: *_y = y;
541: *_k = k;
542: *_m = m;
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: static PetscErrorCode MatSolve_SeqDense_TearDown(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k)
547: {
548: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
549: PetscScalar *y = NULL;
550: PetscBLASInt m, k;
552: PetscFunctionBegin;
553: y = *_y;
554: *_y = NULL;
555: k = *_k;
556: m = *_m;
557: if (k < m) {
558: PetscScalar *yv;
559: PetscCall(VecGetArray(yy, &yv));
560: PetscCall(PetscArraycpy(yv, y, k));
561: PetscCall(VecRestoreArray(yy, &yv));
562: PetscCall(VecRestoreArray(mat->qrrhs, &y));
563: } else {
564: PetscCall(VecRestoreArray(yy, &y));
565: }
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: static PetscErrorCode MatSolve_SeqDense_LU(Mat A, Vec xx, Vec yy)
570: {
571: PetscScalar *y = NULL;
572: PetscBLASInt m = 0, k = 0;
574: PetscFunctionBegin;
575: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
576: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_FALSE));
577: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
581: static PetscErrorCode MatSolveTranspose_SeqDense_LU(Mat A, Vec xx, Vec yy)
582: {
583: PetscScalar *y = NULL;
584: PetscBLASInt m = 0, k = 0;
586: PetscFunctionBegin;
587: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
588: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_TRUE));
589: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: static PetscErrorCode MatSolve_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
594: {
595: PetscScalar *y = NULL;
596: PetscBLASInt m = 0, k = 0;
598: PetscFunctionBegin;
599: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
600: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_FALSE));
601: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: static PetscErrorCode MatSolveTranspose_SeqDense_Cholesky(Mat A, Vec xx, Vec yy)
606: {
607: PetscScalar *y = NULL;
608: PetscBLASInt m = 0, k = 0;
610: PetscFunctionBegin;
611: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
612: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_TRUE));
613: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: static PetscErrorCode MatSolve_SeqDense_QR(Mat A, Vec xx, Vec yy)
618: {
619: PetscScalar *y = NULL;
620: PetscBLASInt m = 0, k = 0;
622: PetscFunctionBegin;
623: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
624: PetscCall(MatSolve_SeqDense_Internal_QR(A, y, PetscMax(m, k), m, 1, k));
625: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
629: static PetscErrorCode MatSolveTranspose_SeqDense_QR(Mat A, Vec xx, Vec yy)
630: {
631: PetscScalar *y = NULL;
632: PetscBLASInt m = 0, k = 0;
634: PetscFunctionBegin;
635: PetscCall(MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k));
636: PetscCall(MatSolveTranspose_SeqDense_Internal_QR(A, y, PetscMax(m, k), m, 1, k));
637: PetscCall(MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k));
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: static PetscErrorCode MatMatSolve_SeqDense_SetUp(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
642: {
643: const PetscScalar *b;
644: PetscScalar *y;
645: PetscInt n, _ldb, _ldx;
646: PetscBLASInt nrhs = 0, m = 0, k = 0, ldb = 0, ldx = 0, ldy = 0;
648: PetscFunctionBegin;
649: *_ldy = 0;
650: *_m = 0;
651: *_nrhs = 0;
652: *_k = 0;
653: *_y = NULL;
654: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
655: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
656: PetscCall(MatGetSize(B, NULL, &n));
657: PetscCall(PetscBLASIntCast(n, &nrhs));
658: PetscCall(MatDenseGetLDA(B, &_ldb));
659: PetscCall(PetscBLASIntCast(_ldb, &ldb));
660: PetscCall(MatDenseGetLDA(X, &_ldx));
661: PetscCall(PetscBLASIntCast(_ldx, &ldx));
662: if (ldx < m) {
663: PetscCall(MatDenseGetArrayRead(B, &b));
664: PetscCall(PetscMalloc1(nrhs * m, &y));
665: if (ldb == m) {
666: PetscCall(PetscArraycpy(y, b, ldb * nrhs));
667: } else {
668: for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&y[j * m], &b[j * ldb], m));
669: }
670: ldy = m;
671: PetscCall(MatDenseRestoreArrayRead(B, &b));
672: } else {
673: if (ldb == ldx) {
674: PetscCall(MatCopy(B, X, SAME_NONZERO_PATTERN));
675: PetscCall(MatDenseGetArray(X, &y));
676: } else {
677: PetscCall(MatDenseGetArray(X, &y));
678: PetscCall(MatDenseGetArrayRead(B, &b));
679: for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&y[j * ldx], &b[j * ldb], m));
680: PetscCall(MatDenseRestoreArrayRead(B, &b));
681: }
682: ldy = ldx;
683: }
684: *_y = y;
685: *_ldy = ldy;
686: *_k = k;
687: *_m = m;
688: *_nrhs = nrhs;
689: PetscFunctionReturn(PETSC_SUCCESS);
690: }
692: static PetscErrorCode MatMatSolve_SeqDense_TearDown(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_ldy, PetscBLASInt *_m, PetscBLASInt *_nrhs, PetscBLASInt *_k)
693: {
694: PetscScalar *y;
695: PetscInt _ldx;
696: PetscBLASInt k, ldy, nrhs, ldx = 0;
698: PetscFunctionBegin;
699: y = *_y;
700: *_y = NULL;
701: k = *_k;
702: ldy = *_ldy;
703: nrhs = *_nrhs;
704: PetscCall(MatDenseGetLDA(X, &_ldx));
705: PetscCall(PetscBLASIntCast(_ldx, &ldx));
706: if (ldx != ldy) {
707: PetscScalar *xv;
708: PetscCall(MatDenseGetArray(X, &xv));
709: for (PetscInt j = 0; j < nrhs; j++) PetscCall(PetscArraycpy(&xv[j * ldx], &y[j * ldy], k));
710: PetscCall(MatDenseRestoreArray(X, &xv));
711: PetscCall(PetscFree(y));
712: } else {
713: PetscCall(MatDenseRestoreArray(X, &y));
714: }
715: PetscFunctionReturn(PETSC_SUCCESS);
716: }
718: static PetscErrorCode MatMatSolve_SeqDense_LU(Mat A, Mat B, Mat X)
719: {
720: PetscScalar *y;
721: PetscBLASInt m, k, ldy, nrhs;
723: PetscFunctionBegin;
724: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
725: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_FALSE));
726: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
727: PetscFunctionReturn(PETSC_SUCCESS);
728: }
730: static PetscErrorCode MatMatSolveTranspose_SeqDense_LU(Mat A, Mat B, Mat X)
731: {
732: PetscScalar *y;
733: PetscBLASInt m, k, ldy, nrhs;
735: PetscFunctionBegin;
736: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
737: PetscCall(MatSolve_SeqDense_Internal_LU(A, y, ldy, m, nrhs, k, PETSC_TRUE));
738: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
739: PetscFunctionReturn(PETSC_SUCCESS);
740: }
742: static PetscErrorCode MatMatSolve_SeqDense_Cholesky(Mat A, Mat B, Mat X)
743: {
744: PetscScalar *y;
745: PetscBLASInt m, k, ldy, nrhs;
747: PetscFunctionBegin;
748: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
749: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_FALSE));
750: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: static PetscErrorCode MatMatSolveTranspose_SeqDense_Cholesky(Mat A, Mat B, Mat X)
755: {
756: PetscScalar *y;
757: PetscBLASInt m, k, ldy, nrhs;
759: PetscFunctionBegin;
760: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
761: PetscCall(MatSolve_SeqDense_Internal_Cholesky(A, y, ldy, m, nrhs, k, PETSC_TRUE));
762: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
763: PetscFunctionReturn(PETSC_SUCCESS);
764: }
766: static PetscErrorCode MatMatSolve_SeqDense_QR(Mat A, Mat B, Mat X)
767: {
768: PetscScalar *y;
769: PetscBLASInt m, k, ldy, nrhs;
771: PetscFunctionBegin;
772: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
773: PetscCall(MatSolve_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k));
774: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
775: PetscFunctionReturn(PETSC_SUCCESS);
776: }
778: static PetscErrorCode MatMatSolveTranspose_SeqDense_QR(Mat A, Mat B, Mat X)
779: {
780: PetscScalar *y;
781: PetscBLASInt m, k, ldy, nrhs;
783: PetscFunctionBegin;
784: PetscCall(MatMatSolve_SeqDense_SetUp(A, B, X, &y, &ldy, &m, &nrhs, &k));
785: PetscCall(MatSolveTranspose_SeqDense_Internal_QR(A, y, ldy, m, nrhs, k));
786: PetscCall(MatMatSolve_SeqDense_TearDown(A, B, X, &y, &ldy, &m, &nrhs, &k));
787: PetscFunctionReturn(PETSC_SUCCESS);
788: }
790: /* COMMENT: I have chosen to hide row permutation in the pivots,
791: rather than put it in the Mat->row slot.*/
792: PetscErrorCode MatLUFactor_SeqDense(Mat A, IS row, IS col, const MatFactorInfo *minfo)
793: {
794: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
795: PetscBLASInt n, m, info;
797: PetscFunctionBegin;
798: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
799: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
800: if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); }
801: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
802: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
803: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&m, &n, mat->v, &mat->lda, mat->pivots, &info));
804: PetscCall(PetscFPTrapPop());
806: PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to LU factorization %d", (int)info);
807: PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Bad LU factorization %d", (int)info);
809: A->ops->solve = MatSolve_SeqDense_LU;
810: A->ops->matsolve = MatMatSolve_SeqDense_LU;
811: A->ops->solvetranspose = MatSolveTranspose_SeqDense_LU;
812: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_LU;
813: A->factortype = MAT_FACTOR_LU;
815: PetscCall(PetscFree(A->solvertype));
816: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
818: PetscCall(PetscLogFlops((2.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3));
819: PetscFunctionReturn(PETSC_SUCCESS);
820: }
822: static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info_dummy)
823: {
824: MatFactorInfo info;
826: PetscFunctionBegin;
827: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
828: PetscUseTypeMethod(fact, lufactor, NULL, NULL, &info);
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
832: PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, IS col, const MatFactorInfo *info)
833: {
834: PetscFunctionBegin;
835: fact->preallocated = PETSC_TRUE;
836: fact->assembled = PETSC_TRUE;
837: fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense;
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }
841: /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */
842: PetscErrorCode MatCholeskyFactor_SeqDense(Mat A, IS perm, const MatFactorInfo *factinfo)
843: {
844: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
845: PetscBLASInt info, n;
847: PetscFunctionBegin;
848: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
849: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
850: if (A->spd == PETSC_BOOL3_TRUE) {
851: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
852: PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("L", &n, mat->v, &mat->lda, &info));
853: PetscCall(PetscFPTrapPop());
854: #if defined(PETSC_USE_COMPLEX)
855: } else if (A->hermitian == PETSC_BOOL3_TRUE) {
856: if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); }
857: if (!mat->fwork) {
858: PetscScalar dummy;
860: mat->lfwork = -1;
861: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
862: PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info));
863: PetscCall(PetscFPTrapPop());
864: mat->lfwork = (PetscInt)PetscRealPart(dummy);
865: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
866: }
867: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
868: PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
869: PetscCall(PetscFPTrapPop());
870: #endif
871: } else { /* symmetric case */
872: if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); }
873: if (!mat->fwork) {
874: PetscScalar dummy;
876: mat->lfwork = -1;
877: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
878: PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info));
879: PetscCall(PetscFPTrapPop());
880: mat->lfwork = (PetscInt)PetscRealPart(dummy);
881: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
882: }
883: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
884: PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info));
885: PetscCall(PetscFPTrapPop());
886: }
887: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad factorization: zero pivot in row %" PetscInt_FMT, (PetscInt)info - 1);
889: A->ops->solve = MatSolve_SeqDense_Cholesky;
890: A->ops->matsolve = MatMatSolve_SeqDense_Cholesky;
891: A->ops->solvetranspose = MatSolveTranspose_SeqDense_Cholesky;
892: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky;
893: A->factortype = MAT_FACTOR_CHOLESKY;
895: PetscCall(PetscFree(A->solvertype));
896: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
898: PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0));
899: PetscFunctionReturn(PETSC_SUCCESS);
900: }
902: static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info_dummy)
903: {
904: MatFactorInfo info;
906: PetscFunctionBegin;
907: info.fill = 1.0;
909: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
910: PetscUseTypeMethod(fact, choleskyfactor, NULL, &info);
911: PetscFunctionReturn(PETSC_SUCCESS);
912: }
914: PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info)
915: {
916: PetscFunctionBegin;
917: fact->assembled = PETSC_TRUE;
918: fact->preallocated = PETSC_TRUE;
919: fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense;
920: PetscFunctionReturn(PETSC_SUCCESS);
921: }
923: PetscErrorCode MatQRFactor_SeqDense(Mat A, IS col, const MatFactorInfo *minfo)
924: {
925: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
926: PetscBLASInt n, m, info, min, max;
928: PetscFunctionBegin;
929: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
930: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
931: max = PetscMax(m, n);
932: min = PetscMin(m, n);
933: if (!mat->tau) { PetscCall(PetscMalloc1(min, &mat->tau)); }
934: if (!mat->pivots) { PetscCall(PetscMalloc1(n, &mat->pivots)); }
935: if (!mat->qrrhs) PetscCall(MatCreateVecs(A, NULL, &(mat->qrrhs)));
936: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
937: if (!mat->fwork) {
938: PetscScalar dummy;
940: mat->lfwork = -1;
941: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
942: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, &dummy, &mat->lfwork, &info));
943: PetscCall(PetscFPTrapPop());
944: mat->lfwork = (PetscInt)PetscRealPart(dummy);
945: PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork));
946: }
947: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
948: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, mat->fwork, &mat->lfwork, &info));
949: PetscCall(PetscFPTrapPop());
950: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to QR factorization %d", (int)info);
951: // 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
952: mat->rank = min;
954: A->ops->solve = MatSolve_SeqDense_QR;
955: A->ops->matsolve = MatMatSolve_SeqDense_QR;
956: A->factortype = MAT_FACTOR_QR;
957: if (m == n) {
958: A->ops->solvetranspose = MatSolveTranspose_SeqDense_QR;
959: A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR;
960: }
962: PetscCall(PetscFree(A->solvertype));
963: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype));
965: PetscCall(PetscLogFlops(2.0 * min * min * (max - min / 3.0)));
966: PetscFunctionReturn(PETSC_SUCCESS);
967: }
969: static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info_dummy)
970: {
971: MatFactorInfo info;
973: PetscFunctionBegin;
974: info.fill = 1.0;
976: PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES));
977: PetscUseMethod(fact, "MatQRFactor_C", (Mat, IS, const MatFactorInfo *), (fact, NULL, &info));
978: PetscFunctionReturn(PETSC_SUCCESS);
979: }
981: PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info)
982: {
983: PetscFunctionBegin;
984: fact->assembled = PETSC_TRUE;
985: fact->preallocated = PETSC_TRUE;
986: PetscCall(PetscObjectComposeFunction((PetscObject)fact, "MatQRFactorNumeric_C", MatQRFactorNumeric_SeqDense));
987: PetscFunctionReturn(PETSC_SUCCESS);
988: }
990: /* uses LAPACK */
991: PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A, MatFactorType ftype, Mat *fact)
992: {
993: PetscFunctionBegin;
994: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), fact));
995: PetscCall(MatSetSizes(*fact, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
996: PetscCall(MatSetType(*fact, MATDENSE));
997: (*fact)->trivialsymbolic = PETSC_TRUE;
998: if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) {
999: (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense;
1000: (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense;
1001: } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) {
1002: (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense;
1003: } else if (ftype == MAT_FACTOR_QR) {
1004: PetscCall(PetscObjectComposeFunction((PetscObject)(*fact), "MatQRFactorSymbolic_C", MatQRFactorSymbolic_SeqDense));
1005: }
1006: (*fact)->factortype = ftype;
1008: PetscCall(PetscFree((*fact)->solvertype));
1009: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*fact)->solvertype));
1010: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_LU]));
1011: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ILU]));
1012: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_CHOLESKY]));
1013: PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ICC]));
1014: PetscFunctionReturn(PETSC_SUCCESS);
1015: }
1017: static PetscErrorCode MatSOR_SeqDense(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal shift, PetscInt its, PetscInt lits, Vec xx)
1018: {
1019: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1020: PetscScalar *x, *v = mat->v, zero = 0.0, xt;
1021: const PetscScalar *b;
1022: PetscInt m = A->rmap->n, i;
1023: PetscBLASInt o = 1, bm = 0;
1025: PetscFunctionBegin;
1026: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1027: PetscCheck(A->offloadmask != PETSC_OFFLOAD_GPU, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented");
1028: #endif
1029: if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */
1030: PetscCall(PetscBLASIntCast(m, &bm));
1031: if (flag & SOR_ZERO_INITIAL_GUESS) {
1032: /* this is a hack fix, should have another version without the second BLASdotu */
1033: PetscCall(VecSet(xx, zero));
1034: }
1035: PetscCall(VecGetArray(xx, &x));
1036: PetscCall(VecGetArrayRead(bb, &b));
1037: its = its * lits;
1038: 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);
1039: while (its--) {
1040: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
1041: for (i = 0; i < m; i++) {
1042: PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o));
1043: x[i] = (1. - omega) * x[i] + omega * (xt + v[i + i * m] * x[i]) / (v[i + i * m] + shift);
1044: }
1045: }
1046: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1047: for (i = m - 1; i >= 0; i--) {
1048: PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o));
1049: x[i] = (1. - omega) * x[i] + omega * (xt + v[i + i * m] * x[i]) / (v[i + i * m] + shift);
1050: }
1051: }
1052: }
1053: PetscCall(VecRestoreArrayRead(bb, &b));
1054: PetscCall(VecRestoreArray(xx, &x));
1055: PetscFunctionReturn(PETSC_SUCCESS);
1056: }
1058: PetscErrorCode MatMultTranspose_SeqDense(Mat A, Vec xx, Vec yy)
1059: {
1060: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1061: const PetscScalar *v = mat->v, *x;
1062: PetscScalar *y;
1063: PetscBLASInt m, n, _One = 1;
1064: PetscScalar _DOne = 1.0, _DZero = 0.0;
1066: PetscFunctionBegin;
1067: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1068: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
1069: PetscCall(VecGetArrayRead(xx, &x));
1070: PetscCall(VecGetArrayWrite(yy, &y));
1071: if (!A->rmap->n || !A->cmap->n) {
1072: PetscBLASInt i;
1073: for (i = 0; i < n; i++) y[i] = 0.0;
1074: } else {
1075: PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v, &mat->lda, x, &_One, &_DZero, y, &_One));
1076: PetscCall(PetscLogFlops(2.0 * A->rmap->n * A->cmap->n - A->cmap->n));
1077: }
1078: PetscCall(VecRestoreArrayRead(xx, &x));
1079: PetscCall(VecRestoreArrayWrite(yy, &y));
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: PetscErrorCode MatMult_SeqDense(Mat A, Vec xx, Vec yy)
1084: {
1085: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1086: PetscScalar *y, _DOne = 1.0, _DZero = 0.0;
1087: PetscBLASInt m, n, _One = 1;
1088: const PetscScalar *v = mat->v, *x;
1090: PetscFunctionBegin;
1091: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1092: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
1093: PetscCall(VecGetArrayRead(xx, &x));
1094: PetscCall(VecGetArrayWrite(yy, &y));
1095: if (!A->rmap->n || !A->cmap->n) {
1096: PetscBLASInt i;
1097: for (i = 0; i < m; i++) y[i] = 0.0;
1098: } else {
1099: PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v, &(mat->lda), x, &_One, &_DZero, y, &_One));
1100: PetscCall(PetscLogFlops(2.0 * A->rmap->n * A->cmap->n - A->rmap->n));
1101: }
1102: PetscCall(VecRestoreArrayRead(xx, &x));
1103: PetscCall(VecRestoreArrayWrite(yy, &y));
1104: PetscFunctionReturn(PETSC_SUCCESS);
1105: }
1107: PetscErrorCode MatMultAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1108: {
1109: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1110: const PetscScalar *v = mat->v, *x;
1111: PetscScalar *y, _DOne = 1.0;
1112: PetscBLASInt m, n, _One = 1;
1114: PetscFunctionBegin;
1115: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1116: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
1117: PetscCall(VecCopy(zz, yy));
1118: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
1119: PetscCall(VecGetArrayRead(xx, &x));
1120: PetscCall(VecGetArray(yy, &y));
1121: PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v, &(mat->lda), x, &_One, &_DOne, y, &_One));
1122: PetscCall(VecRestoreArrayRead(xx, &x));
1123: PetscCall(VecRestoreArray(yy, &y));
1124: PetscCall(PetscLogFlops(2.0 * A->rmap->n * A->cmap->n));
1125: PetscFunctionReturn(PETSC_SUCCESS);
1126: }
1128: PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy)
1129: {
1130: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1131: const PetscScalar *v = mat->v, *x;
1132: PetscScalar *y;
1133: PetscBLASInt m, n, _One = 1;
1134: PetscScalar _DOne = 1.0;
1136: PetscFunctionBegin;
1137: PetscCall(PetscBLASIntCast(A->rmap->n, &m));
1138: PetscCall(PetscBLASIntCast(A->cmap->n, &n));
1139: PetscCall(VecCopy(zz, yy));
1140: if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS);
1141: PetscCall(VecGetArrayRead(xx, &x));
1142: PetscCall(VecGetArray(yy, &y));
1143: PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v, &(mat->lda), x, &_One, &_DOne, y, &_One));
1144: PetscCall(VecRestoreArrayRead(xx, &x));
1145: PetscCall(VecRestoreArray(yy, &y));
1146: PetscCall(PetscLogFlops(2.0 * A->rmap->n * A->cmap->n));
1147: PetscFunctionReturn(PETSC_SUCCESS);
1148: }
1150: static PetscErrorCode MatGetRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1151: {
1152: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1153: PetscInt i;
1155: PetscFunctionBegin;
1156: if (ncols) *ncols = A->cmap->n;
1157: if (cols) {
1158: PetscCall(PetscMalloc1(A->cmap->n, cols));
1159: for (i = 0; i < A->cmap->n; i++) (*cols)[i] = i;
1160: }
1161: if (vals) {
1162: const PetscScalar *v;
1164: PetscCall(MatDenseGetArrayRead(A, &v));
1165: PetscCall(PetscMalloc1(A->cmap->n, vals));
1166: v += row;
1167: for (i = 0; i < A->cmap->n; i++) {
1168: (*vals)[i] = *v;
1169: v += mat->lda;
1170: }
1171: PetscCall(MatDenseRestoreArrayRead(A, &v));
1172: }
1173: PetscFunctionReturn(PETSC_SUCCESS);
1174: }
1176: static PetscErrorCode MatRestoreRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals)
1177: {
1178: PetscFunctionBegin;
1179: if (cols) PetscCall(PetscFree(*cols));
1180: if (vals) PetscCall(PetscFree(*vals));
1181: PetscFunctionReturn(PETSC_SUCCESS);
1182: }
1184: static PetscErrorCode MatSetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], const PetscScalar v[], InsertMode addv)
1185: {
1186: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1187: PetscScalar *av;
1188: PetscInt i, j, idx = 0;
1189: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1190: PetscOffloadMask oldf;
1191: #endif
1193: PetscFunctionBegin;
1194: PetscCall(MatDenseGetArray(A, &av));
1195: if (!mat->roworiented) {
1196: if (addv == INSERT_VALUES) {
1197: for (j = 0; j < n; j++) {
1198: if (indexn[j] < 0) {
1199: idx += m;
1200: continue;
1201: }
1202: 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);
1203: for (i = 0; i < m; i++) {
1204: if (indexm[i] < 0) {
1205: idx++;
1206: continue;
1207: }
1208: 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);
1209: av[indexn[j] * mat->lda + indexm[i]] = v[idx++];
1210: }
1211: }
1212: } else {
1213: for (j = 0; j < n; j++) {
1214: if (indexn[j] < 0) {
1215: idx += m;
1216: continue;
1217: }
1218: 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);
1219: for (i = 0; i < m; i++) {
1220: if (indexm[i] < 0) {
1221: idx++;
1222: continue;
1223: }
1224: 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);
1225: av[indexn[j] * mat->lda + indexm[i]] += v[idx++];
1226: }
1227: }
1228: }
1229: } else {
1230: if (addv == INSERT_VALUES) {
1231: for (i = 0; i < m; i++) {
1232: if (indexm[i] < 0) {
1233: idx += n;
1234: continue;
1235: }
1236: 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);
1237: for (j = 0; j < n; j++) {
1238: if (indexn[j] < 0) {
1239: idx++;
1240: continue;
1241: }
1242: 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);
1243: av[indexn[j] * mat->lda + indexm[i]] = v[idx++];
1244: }
1245: }
1246: } else {
1247: for (i = 0; i < m; i++) {
1248: if (indexm[i] < 0) {
1249: idx += n;
1250: continue;
1251: }
1252: 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);
1253: for (j = 0; j < n; j++) {
1254: if (indexn[j] < 0) {
1255: idx++;
1256: continue;
1257: }
1258: 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);
1259: av[indexn[j] * mat->lda + indexm[i]] += v[idx++];
1260: }
1261: }
1262: }
1263: }
1264: /* hack to prevent unneeded copy to the GPU while returning the array */
1265: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1266: oldf = A->offloadmask;
1267: A->offloadmask = PETSC_OFFLOAD_GPU;
1268: #endif
1269: PetscCall(MatDenseRestoreArray(A, &av));
1270: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1271: A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU);
1272: #endif
1273: PetscFunctionReturn(PETSC_SUCCESS);
1274: }
1276: static PetscErrorCode MatGetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], PetscScalar v[])
1277: {
1278: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1279: const PetscScalar *vv;
1280: PetscInt i, j;
1282: PetscFunctionBegin;
1283: PetscCall(MatDenseGetArrayRead(A, &vv));
1284: /* row-oriented output */
1285: for (i = 0; i < m; i++) {
1286: if (indexm[i] < 0) {
1287: v += n;
1288: continue;
1289: }
1290: 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);
1291: for (j = 0; j < n; j++) {
1292: if (indexn[j] < 0) {
1293: v++;
1294: continue;
1295: }
1296: 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);
1297: *v++ = vv[indexn[j] * mat->lda + indexm[i]];
1298: }
1299: }
1300: PetscCall(MatDenseRestoreArrayRead(A, &vv));
1301: PetscFunctionReturn(PETSC_SUCCESS);
1302: }
1304: PetscErrorCode MatView_Dense_Binary(Mat mat, PetscViewer viewer)
1305: {
1306: PetscBool skipHeader;
1307: PetscViewerFormat format;
1308: PetscInt header[4], M, N, m, lda, i, j, k;
1309: const PetscScalar *v;
1310: PetscScalar *vwork;
1312: PetscFunctionBegin;
1313: PetscCall(PetscViewerSetUp(viewer));
1314: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
1315: PetscCall(PetscViewerGetFormat(viewer, &format));
1316: if (skipHeader) format = PETSC_VIEWER_NATIVE;
1318: PetscCall(MatGetSize(mat, &M, &N));
1320: /* write matrix header */
1321: header[0] = MAT_FILE_CLASSID;
1322: header[1] = M;
1323: header[2] = N;
1324: header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M * N;
1325: if (!skipHeader) PetscCall(PetscViewerBinaryWrite(viewer, header, 4, PETSC_INT));
1327: PetscCall(MatGetLocalSize(mat, &m, NULL));
1328: if (format != PETSC_VIEWER_NATIVE) {
1329: PetscInt nnz = m * N, *iwork;
1330: /* store row lengths for each row */
1331: PetscCall(PetscMalloc1(nnz, &iwork));
1332: for (i = 0; i < m; i++) iwork[i] = N;
1333: PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1334: /* store column indices (zero start index) */
1335: for (k = 0, i = 0; i < m; i++)
1336: for (j = 0; j < N; j++, k++) iwork[k] = j;
1337: PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1338: PetscCall(PetscFree(iwork));
1339: }
1340: /* store matrix values as a dense matrix in row major order */
1341: PetscCall(PetscMalloc1(m * N, &vwork));
1342: PetscCall(MatDenseGetArrayRead(mat, &v));
1343: PetscCall(MatDenseGetLDA(mat, &lda));
1344: for (k = 0, i = 0; i < m; i++)
1345: for (j = 0; j < N; j++, k++) vwork[k] = v[i + lda * j];
1346: PetscCall(MatDenseRestoreArrayRead(mat, &v));
1347: PetscCall(PetscViewerBinaryWriteAll(viewer, vwork, m * N, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1348: PetscCall(PetscFree(vwork));
1349: PetscFunctionReturn(PETSC_SUCCESS);
1350: }
1352: PetscErrorCode MatLoad_Dense_Binary(Mat mat, PetscViewer viewer)
1353: {
1354: PetscBool skipHeader;
1355: PetscInt header[4], M, N, m, nz, lda, i, j, k;
1356: PetscInt rows, cols;
1357: PetscScalar *v, *vwork;
1359: PetscFunctionBegin;
1360: PetscCall(PetscViewerSetUp(viewer));
1361: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader));
1363: if (!skipHeader) {
1364: PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT));
1365: PetscCheck(header[0] == MAT_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file");
1366: M = header[1];
1367: N = header[2];
1368: PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M);
1369: PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N);
1370: nz = header[3];
1371: PetscCheck(nz == MATRIX_BINARY_FORMAT_DENSE || nz >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Unknown matrix format %" PetscInt_FMT " in file", nz);
1372: } else {
1373: PetscCall(MatGetSize(mat, &M, &N));
1374: 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");
1375: nz = MATRIX_BINARY_FORMAT_DENSE;
1376: }
1378: /* setup global sizes if not set */
1379: if (mat->rmap->N < 0) mat->rmap->N = M;
1380: if (mat->cmap->N < 0) mat->cmap->N = N;
1381: PetscCall(MatSetUp(mat));
1382: /* check if global sizes are correct */
1383: PetscCall(MatGetSize(mat, &rows, &cols));
1384: 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);
1386: PetscCall(MatGetSize(mat, NULL, &N));
1387: PetscCall(MatGetLocalSize(mat, &m, NULL));
1388: PetscCall(MatDenseGetArray(mat, &v));
1389: PetscCall(MatDenseGetLDA(mat, &lda));
1390: if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */
1391: PetscInt nnz = m * N;
1392: /* read in matrix values */
1393: PetscCall(PetscMalloc1(nnz, &vwork));
1394: PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1395: /* store values in column major order */
1396: for (j = 0; j < N; j++)
1397: for (i = 0; i < m; i++) v[i + lda * j] = vwork[i * N + j];
1398: PetscCall(PetscFree(vwork));
1399: } else { /* matrix in file is sparse format */
1400: PetscInt nnz = 0, *rlens, *icols;
1401: /* read in row lengths */
1402: PetscCall(PetscMalloc1(m, &rlens));
1403: PetscCall(PetscViewerBinaryReadAll(viewer, rlens, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1404: for (i = 0; i < m; i++) nnz += rlens[i];
1405: /* read in column indices and values */
1406: PetscCall(PetscMalloc2(nnz, &icols, nnz, &vwork));
1407: PetscCall(PetscViewerBinaryReadAll(viewer, icols, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT));
1408: PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR));
1409: /* store values in column major order */
1410: for (k = 0, i = 0; i < m; i++)
1411: for (j = 0; j < rlens[i]; j++, k++) v[i + lda * icols[k]] = vwork[k];
1412: PetscCall(PetscFree(rlens));
1413: PetscCall(PetscFree2(icols, vwork));
1414: }
1415: PetscCall(MatDenseRestoreArray(mat, &v));
1416: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
1417: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
1418: PetscFunctionReturn(PETSC_SUCCESS);
1419: }
1421: static PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer)
1422: {
1423: PetscBool isbinary, ishdf5;
1425: PetscFunctionBegin;
1428: /* force binary viewer to load .info file if it has not yet done so */
1429: PetscCall(PetscViewerSetUp(viewer));
1430: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1431: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1432: if (isbinary) {
1433: PetscCall(MatLoad_Dense_Binary(newMat, viewer));
1434: } else if (ishdf5) {
1435: #if defined(PETSC_HAVE_HDF5)
1436: PetscCall(MatLoad_Dense_HDF5(newMat, viewer));
1437: #else
1438: SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1439: #endif
1440: } else {
1441: 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);
1442: }
1443: PetscFunctionReturn(PETSC_SUCCESS);
1444: }
1446: static PetscErrorCode MatView_SeqDense_ASCII(Mat A, PetscViewer viewer)
1447: {
1448: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1449: PetscInt i, j;
1450: const char *name;
1451: PetscScalar *v, *av;
1452: PetscViewerFormat format;
1453: #if defined(PETSC_USE_COMPLEX)
1454: PetscBool allreal = PETSC_TRUE;
1455: #endif
1457: PetscFunctionBegin;
1458: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&av));
1459: PetscCall(PetscViewerGetFormat(viewer, &format));
1460: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1461: PetscFunctionReturn(PETSC_SUCCESS); /* do nothing for now */
1462: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1463: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1464: for (i = 0; i < A->rmap->n; i++) {
1465: v = av + i;
1466: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i));
1467: for (j = 0; j < A->cmap->n; j++) {
1468: #if defined(PETSC_USE_COMPLEX)
1469: if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) {
1470: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", j, (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1471: } else if (PetscRealPart(*v)) {
1472: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)PetscRealPart(*v)));
1473: }
1474: #else
1475: if (*v) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)*v));
1476: #endif
1477: v += a->lda;
1478: }
1479: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1480: }
1481: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1482: } else {
1483: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1484: #if defined(PETSC_USE_COMPLEX)
1485: /* determine if matrix has all real values */
1486: for (j = 0; j < A->cmap->n; j++) {
1487: v = av + j * a->lda;
1488: for (i = 0; i < A->rmap->n; i++) {
1489: if (PetscImaginaryPart(v[i])) {
1490: allreal = PETSC_FALSE;
1491: break;
1492: }
1493: }
1494: }
1495: #endif
1496: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1497: PetscCall(PetscObjectGetName((PetscObject)A, &name));
1498: PetscCall(PetscViewerASCIIPrintf(viewer, "%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n", A->rmap->n, A->cmap->n));
1499: PetscCall(PetscViewerASCIIPrintf(viewer, "%s = zeros(%" PetscInt_FMT ",%" PetscInt_FMT ");\n", name, A->rmap->n, A->cmap->n));
1500: PetscCall(PetscViewerASCIIPrintf(viewer, "%s = [\n", name));
1501: }
1503: for (i = 0; i < A->rmap->n; i++) {
1504: v = av + i;
1505: for (j = 0; j < A->cmap->n; j++) {
1506: #if defined(PETSC_USE_COMPLEX)
1507: if (allreal) {
1508: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)PetscRealPart(*v)));
1509: } else {
1510: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei ", (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v)));
1511: }
1512: #else
1513: PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)*v));
1514: #endif
1515: v += a->lda;
1516: }
1517: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1518: }
1519: if (format == PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer, "];\n"));
1520: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1521: }
1522: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&av));
1523: PetscCall(PetscViewerFlush(viewer));
1524: PetscFunctionReturn(PETSC_SUCCESS);
1525: }
1527: #include <petscdraw.h>
1528: static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw, void *Aa)
1529: {
1530: Mat A = (Mat)Aa;
1531: PetscInt m = A->rmap->n, n = A->cmap->n, i, j;
1532: int color = PETSC_DRAW_WHITE;
1533: const PetscScalar *v;
1534: PetscViewer viewer;
1535: PetscReal xl, yl, xr, yr, x_l, x_r, y_l, y_r;
1536: PetscViewerFormat format;
1538: PetscFunctionBegin;
1539: PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer));
1540: PetscCall(PetscViewerGetFormat(viewer, &format));
1541: PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));
1543: /* Loop over matrix elements drawing boxes */
1544: PetscCall(MatDenseGetArrayRead(A, &v));
1545: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1546: PetscDrawCollectiveBegin(draw);
1547: /* Blue for negative and Red for positive */
1548: for (j = 0; j < n; j++) {
1549: x_l = j;
1550: x_r = x_l + 1.0;
1551: for (i = 0; i < m; i++) {
1552: y_l = m - i - 1.0;
1553: y_r = y_l + 1.0;
1554: if (PetscRealPart(v[j * m + i]) > 0.) color = PETSC_DRAW_RED;
1555: else if (PetscRealPart(v[j * m + i]) < 0.) color = PETSC_DRAW_BLUE;
1556: else continue;
1557: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1558: }
1559: }
1560: PetscDrawCollectiveEnd(draw);
1561: } else {
1562: /* use contour shading to indicate magnitude of values */
1563: /* first determine max of all nonzero values */
1564: PetscReal minv = 0.0, maxv = 0.0;
1565: PetscDraw popup;
1567: for (i = 0; i < m * n; i++) {
1568: if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]);
1569: }
1570: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1571: PetscCall(PetscDrawGetPopup(draw, &popup));
1572: PetscCall(PetscDrawScalePopup(popup, minv, maxv));
1574: PetscDrawCollectiveBegin(draw);
1575: for (j = 0; j < n; j++) {
1576: x_l = j;
1577: x_r = x_l + 1.0;
1578: for (i = 0; i < m; i++) {
1579: y_l = m - i - 1.0;
1580: y_r = y_l + 1.0;
1581: color = PetscDrawRealToColor(PetscAbsScalar(v[j * m + i]), minv, maxv);
1582: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1583: }
1584: }
1585: PetscDrawCollectiveEnd(draw);
1586: }
1587: PetscCall(MatDenseRestoreArrayRead(A, &v));
1588: PetscFunctionReturn(PETSC_SUCCESS);
1589: }
1591: static PetscErrorCode MatView_SeqDense_Draw(Mat A, PetscViewer viewer)
1592: {
1593: PetscDraw draw;
1594: PetscBool isnull;
1595: PetscReal xr, yr, xl, yl, h, w;
1597: PetscFunctionBegin;
1598: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1599: PetscCall(PetscDrawIsNull(draw, &isnull));
1600: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
1602: xr = A->cmap->n;
1603: yr = A->rmap->n;
1604: h = yr / 10.0;
1605: w = xr / 10.0;
1606: xr += w;
1607: yr += h;
1608: xl = -w;
1609: yl = -h;
1610: PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr));
1611: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer));
1612: PetscCall(PetscDrawZoom(draw, MatView_SeqDense_Draw_Zoom, A));
1613: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL));
1614: PetscCall(PetscDrawSave(draw));
1615: PetscFunctionReturn(PETSC_SUCCESS);
1616: }
1618: PetscErrorCode MatView_SeqDense(Mat A, PetscViewer viewer)
1619: {
1620: PetscBool iascii, isbinary, isdraw;
1622: PetscFunctionBegin;
1623: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1624: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1625: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1626: if (iascii) PetscCall(MatView_SeqDense_ASCII(A, viewer));
1627: else if (isbinary) PetscCall(MatView_Dense_Binary(A, viewer));
1628: else if (isdraw) PetscCall(MatView_SeqDense_Draw(A, viewer));
1629: PetscFunctionReturn(PETSC_SUCCESS);
1630: }
1632: static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A, const PetscScalar *array)
1633: {
1634: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1636: PetscFunctionBegin;
1637: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1638: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1639: PetscCheck(!a->unplacedarray, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreArray() first");
1640: a->unplacedarray = a->v;
1641: a->unplaced_user_alloc = a->user_alloc;
1642: a->v = (PetscScalar *)array;
1643: a->user_alloc = PETSC_TRUE;
1644: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1645: A->offloadmask = PETSC_OFFLOAD_CPU;
1646: #endif
1647: PetscFunctionReturn(PETSC_SUCCESS);
1648: }
1650: static PetscErrorCode MatDenseResetArray_SeqDense(Mat A)
1651: {
1652: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1654: PetscFunctionBegin;
1655: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1656: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1657: a->v = a->unplacedarray;
1658: a->user_alloc = a->unplaced_user_alloc;
1659: a->unplacedarray = NULL;
1660: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1661: A->offloadmask = PETSC_OFFLOAD_CPU;
1662: #endif
1663: PetscFunctionReturn(PETSC_SUCCESS);
1664: }
1666: static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A, const PetscScalar *array)
1667: {
1668: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
1670: PetscFunctionBegin;
1671: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1672: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1673: if (!a->user_alloc) PetscCall(PetscFree(a->v));
1674: a->v = (PetscScalar *)array;
1675: a->user_alloc = PETSC_FALSE;
1676: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1677: A->offloadmask = PETSC_OFFLOAD_CPU;
1678: #endif
1679: PetscFunctionReturn(PETSC_SUCCESS);
1680: }
1682: PetscErrorCode MatDestroy_SeqDense(Mat mat)
1683: {
1684: Mat_SeqDense *l = (Mat_SeqDense *)mat->data;
1686: PetscFunctionBegin;
1687: PetscCall(PetscLogObjectState((PetscObject)mat, "Rows %" PetscInt_FMT " Cols %" PetscInt_FMT, mat->rmap->n, mat->cmap->n));
1688: PetscCall(VecDestroy(&(l->qrrhs)));
1689: PetscCall(PetscFree(l->tau));
1690: PetscCall(PetscFree(l->pivots));
1691: PetscCall(PetscFree(l->fwork));
1692: if (!l->user_alloc) PetscCall(PetscFree(l->v));
1693: if (!l->unplaced_user_alloc) PetscCall(PetscFree(l->unplacedarray));
1694: PetscCheck(!l->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
1695: PetscCheck(!l->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
1696: PetscCall(VecDestroy(&l->cvec));
1697: PetscCall(MatDestroy(&l->cmat));
1698: PetscCall(PetscFree(mat->data));
1700: PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL));
1701: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactor_C", NULL));
1702: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorSymbolic_C", NULL));
1703: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorNumeric_C", NULL));
1704: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetLDA_C", NULL));
1705: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseSetLDA_C", NULL));
1706: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArray_C", NULL));
1707: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArray_C", NULL));
1708: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDensePlaceArray_C", NULL));
1709: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseResetArray_C", NULL));
1710: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseReplaceArray_C", NULL));
1711: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayRead_C", NULL));
1712: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayRead_C", NULL));
1713: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayWrite_C", NULL));
1714: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayWrite_C", NULL));
1715: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqaij_C", NULL));
1716: #if defined(PETSC_HAVE_ELEMENTAL)
1717: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_elemental_C", NULL));
1718: #endif
1719: #if defined(PETSC_HAVE_SCALAPACK)
1720: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_scalapack_C", NULL));
1721: #endif
1722: #if defined(PETSC_HAVE_CUDA)
1723: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensecuda_C", NULL));
1724: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", NULL));
1725: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdense_C", NULL));
1726: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensecuda_C", NULL));
1727: #endif
1728: #if defined(PETSC_HAVE_HIP)
1729: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensehip_C", NULL));
1730: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", NULL));
1731: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdense_C", NULL));
1732: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensehip_C", NULL));
1733: #endif
1734: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatSeqDenseSetPreallocation_C", NULL));
1735: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqaij_seqdense_C", NULL));
1736: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdense_C", NULL));
1737: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqbaij_seqdense_C", NULL));
1738: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqsbaij_seqdense_C", NULL));
1740: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumn_C", NULL));
1741: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumn_C", NULL));
1742: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVec_C", NULL));
1743: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVec_C", NULL));
1744: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecRead_C", NULL));
1745: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecRead_C", NULL));
1746: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecWrite_C", NULL));
1747: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecWrite_C", NULL));
1748: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetSubMatrix_C", NULL));
1749: PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreSubMatrix_C", NULL));
1750: PetscFunctionReturn(PETSC_SUCCESS);
1751: }
1753: static PetscErrorCode MatTranspose_SeqDense(Mat A, MatReuse reuse, Mat *matout)
1754: {
1755: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1756: PetscInt k, j, m = A->rmap->n, M = mat->lda, n = A->cmap->n;
1757: PetscScalar *v, tmp;
1759: PetscFunctionBegin;
1760: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *matout));
1761: if (reuse == MAT_INPLACE_MATRIX) {
1762: if (m == n) { /* in place transpose */
1763: PetscCall(MatDenseGetArray(A, &v));
1764: for (j = 0; j < m; j++) {
1765: for (k = 0; k < j; k++) {
1766: tmp = v[j + k * M];
1767: v[j + k * M] = v[k + j * M];
1768: v[k + j * M] = tmp;
1769: }
1770: }
1771: PetscCall(MatDenseRestoreArray(A, &v));
1772: } else { /* reuse memory, temporary allocates new memory */
1773: PetscScalar *v2;
1774: PetscLayout tmplayout;
1776: PetscCall(PetscMalloc1((size_t)m * n, &v2));
1777: PetscCall(MatDenseGetArray(A, &v));
1778: for (j = 0; j < n; j++) {
1779: for (k = 0; k < m; k++) v2[j + (size_t)k * n] = v[k + (size_t)j * M];
1780: }
1781: PetscCall(PetscArraycpy(v, v2, (size_t)m * n));
1782: PetscCall(PetscFree(v2));
1783: PetscCall(MatDenseRestoreArray(A, &v));
1784: /* cleanup size dependent quantities */
1785: PetscCall(VecDestroy(&mat->cvec));
1786: PetscCall(MatDestroy(&mat->cmat));
1787: PetscCall(PetscFree(mat->pivots));
1788: PetscCall(PetscFree(mat->fwork));
1789: /* swap row/col layouts */
1790: mat->lda = n;
1791: tmplayout = A->rmap;
1792: A->rmap = A->cmap;
1793: A->cmap = tmplayout;
1794: }
1795: } else { /* out-of-place transpose */
1796: Mat tmat;
1797: Mat_SeqDense *tmatd;
1798: PetscScalar *v2;
1799: PetscInt M2;
1801: if (reuse == MAT_INITIAL_MATRIX) {
1802: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &tmat));
1803: PetscCall(MatSetSizes(tmat, A->cmap->n, A->rmap->n, A->cmap->n, A->rmap->n));
1804: PetscCall(MatSetType(tmat, ((PetscObject)A)->type_name));
1805: PetscCall(MatSeqDenseSetPreallocation(tmat, NULL));
1806: } else tmat = *matout;
1808: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&v));
1809: PetscCall(MatDenseGetArray(tmat, &v2));
1810: tmatd = (Mat_SeqDense *)tmat->data;
1811: M2 = tmatd->lda;
1812: for (j = 0; j < n; j++) {
1813: for (k = 0; k < m; k++) v2[j + k * M2] = v[k + j * M];
1814: }
1815: PetscCall(MatDenseRestoreArray(tmat, &v2));
1816: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&v));
1817: PetscCall(MatAssemblyBegin(tmat, MAT_FINAL_ASSEMBLY));
1818: PetscCall(MatAssemblyEnd(tmat, MAT_FINAL_ASSEMBLY));
1819: *matout = tmat;
1820: }
1821: PetscFunctionReturn(PETSC_SUCCESS);
1822: }
1824: static PetscErrorCode MatEqual_SeqDense(Mat A1, Mat A2, PetscBool *flg)
1825: {
1826: Mat_SeqDense *mat1 = (Mat_SeqDense *)A1->data;
1827: Mat_SeqDense *mat2 = (Mat_SeqDense *)A2->data;
1828: PetscInt i;
1829: const PetscScalar *v1, *v2;
1831: PetscFunctionBegin;
1832: if (A1->rmap->n != A2->rmap->n) {
1833: *flg = PETSC_FALSE;
1834: PetscFunctionReturn(PETSC_SUCCESS);
1835: }
1836: if (A1->cmap->n != A2->cmap->n) {
1837: *flg = PETSC_FALSE;
1838: PetscFunctionReturn(PETSC_SUCCESS);
1839: }
1840: PetscCall(MatDenseGetArrayRead(A1, &v1));
1841: PetscCall(MatDenseGetArrayRead(A2, &v2));
1842: for (i = 0; i < A1->cmap->n; i++) {
1843: PetscCall(PetscArraycmp(v1, v2, A1->rmap->n, flg));
1844: if (*flg == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS);
1845: v1 += mat1->lda;
1846: v2 += mat2->lda;
1847: }
1848: PetscCall(MatDenseRestoreArrayRead(A1, &v1));
1849: PetscCall(MatDenseRestoreArrayRead(A2, &v2));
1850: *flg = PETSC_TRUE;
1851: PetscFunctionReturn(PETSC_SUCCESS);
1852: }
1854: PetscErrorCode MatGetDiagonal_SeqDense(Mat A, Vec v)
1855: {
1856: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1857: PetscInt i, n, len;
1858: PetscScalar *x;
1859: const PetscScalar *vv;
1861: PetscFunctionBegin;
1862: PetscCall(VecGetSize(v, &n));
1863: PetscCall(VecGetArray(v, &x));
1864: len = PetscMin(A->rmap->n, A->cmap->n);
1865: PetscCall(MatDenseGetArrayRead(A, &vv));
1866: PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming mat and vec");
1867: for (i = 0; i < len; i++) x[i] = vv[i * mat->lda + i];
1868: PetscCall(MatDenseRestoreArrayRead(A, &vv));
1869: PetscCall(VecRestoreArray(v, &x));
1870: PetscFunctionReturn(PETSC_SUCCESS);
1871: }
1873: static PetscErrorCode MatDiagonalScale_SeqDense(Mat A, Vec ll, Vec rr)
1874: {
1875: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1876: const PetscScalar *l, *r;
1877: PetscScalar x, *v, *vv;
1878: PetscInt i, j, m = A->rmap->n, n = A->cmap->n;
1880: PetscFunctionBegin;
1881: PetscCall(MatDenseGetArray(A, &vv));
1882: if (ll) {
1883: PetscCall(VecGetSize(ll, &m));
1884: PetscCall(VecGetArrayRead(ll, &l));
1885: PetscCheck(m == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vec wrong size");
1886: for (i = 0; i < m; i++) {
1887: x = l[i];
1888: v = vv + i;
1889: for (j = 0; j < n; j++) {
1890: (*v) *= x;
1891: v += mat->lda;
1892: }
1893: }
1894: PetscCall(VecRestoreArrayRead(ll, &l));
1895: PetscCall(PetscLogFlops(1.0 * n * m));
1896: }
1897: if (rr) {
1898: PetscCall(VecGetSize(rr, &n));
1899: PetscCall(VecGetArrayRead(rr, &r));
1900: PetscCheck(n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vec wrong size");
1901: for (i = 0; i < n; i++) {
1902: x = r[i];
1903: v = vv + i * mat->lda;
1904: for (j = 0; j < m; j++) (*v++) *= x;
1905: }
1906: PetscCall(VecRestoreArrayRead(rr, &r));
1907: PetscCall(PetscLogFlops(1.0 * n * m));
1908: }
1909: PetscCall(MatDenseRestoreArray(A, &vv));
1910: PetscFunctionReturn(PETSC_SUCCESS);
1911: }
1913: PetscErrorCode MatNorm_SeqDense(Mat A, NormType type, PetscReal *nrm)
1914: {
1915: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
1916: PetscScalar *v, *vv;
1917: PetscReal sum = 0.0;
1918: PetscInt lda, m = A->rmap->n, i, j;
1920: PetscFunctionBegin;
1921: PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&vv));
1922: PetscCall(MatDenseGetLDA(A, &lda));
1923: v = vv;
1924: if (type == NORM_FROBENIUS) {
1925: if (lda > m) {
1926: for (j = 0; j < A->cmap->n; j++) {
1927: v = vv + j * lda;
1928: for (i = 0; i < m; i++) {
1929: sum += PetscRealPart(PetscConj(*v) * (*v));
1930: v++;
1931: }
1932: }
1933: } else {
1934: #if defined(PETSC_USE_REAL___FP16)
1935: PetscBLASInt one = 1, cnt = A->cmap->n * A->rmap->n;
1936: PetscCallBLAS("BLASnrm2", *nrm = BLASnrm2_(&cnt, v, &one));
1937: }
1938: #else
1939: for (i = 0; i < A->cmap->n * A->rmap->n; i++) {
1940: sum += PetscRealPart(PetscConj(*v) * (*v));
1941: v++;
1942: }
1943: }
1944: *nrm = PetscSqrtReal(sum);
1945: #endif
1946: PetscCall(PetscLogFlops(2.0 * A->cmap->n * A->rmap->n));
1947: } else if (type == NORM_1) {
1948: *nrm = 0.0;
1949: for (j = 0; j < A->cmap->n; j++) {
1950: v = vv + j * mat->lda;
1951: sum = 0.0;
1952: for (i = 0; i < A->rmap->n; i++) {
1953: sum += PetscAbsScalar(*v);
1954: v++;
1955: }
1956: if (sum > *nrm) *nrm = sum;
1957: }
1958: PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
1959: } else if (type == NORM_INFINITY) {
1960: *nrm = 0.0;
1961: for (j = 0; j < A->rmap->n; j++) {
1962: v = vv + j;
1963: sum = 0.0;
1964: for (i = 0; i < A->cmap->n; i++) {
1965: sum += PetscAbsScalar(*v);
1966: v += mat->lda;
1967: }
1968: if (sum > *nrm) *nrm = sum;
1969: }
1970: PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n));
1971: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No two norm");
1972: PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&vv));
1973: PetscFunctionReturn(PETSC_SUCCESS);
1974: }
1976: static PetscErrorCode MatSetOption_SeqDense(Mat A, MatOption op, PetscBool flg)
1977: {
1978: Mat_SeqDense *aij = (Mat_SeqDense *)A->data;
1980: PetscFunctionBegin;
1981: switch (op) {
1982: case MAT_ROW_ORIENTED:
1983: aij->roworiented = flg;
1984: break;
1985: case MAT_NEW_NONZERO_LOCATIONS:
1986: case MAT_NEW_NONZERO_LOCATION_ERR:
1987: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1988: case MAT_FORCE_DIAGONAL_ENTRIES:
1989: case MAT_KEEP_NONZERO_PATTERN:
1990: case MAT_IGNORE_OFF_PROC_ENTRIES:
1991: case MAT_USE_HASH_TABLE:
1992: case MAT_IGNORE_ZERO_ENTRIES:
1993: case MAT_IGNORE_LOWER_TRIANGULAR:
1994: case MAT_SORTED_FULL:
1995: PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op]));
1996: break;
1997: case MAT_SPD:
1998: case MAT_SYMMETRIC:
1999: case MAT_STRUCTURALLY_SYMMETRIC:
2000: case MAT_HERMITIAN:
2001: case MAT_SYMMETRY_ETERNAL:
2002: case MAT_STRUCTURAL_SYMMETRY_ETERNAL:
2003: case MAT_SPD_ETERNAL:
2004: break;
2005: default:
2006: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unknown option %s", MatOptions[op]);
2007: }
2008: PetscFunctionReturn(PETSC_SUCCESS);
2009: }
2011: PetscErrorCode MatZeroEntries_SeqDense(Mat A)
2012: {
2013: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
2014: PetscInt lda = l->lda, m = A->rmap->n, n = A->cmap->n, j;
2015: PetscScalar *v;
2017: PetscFunctionBegin;
2018: PetscCall(MatDenseGetArrayWrite(A, &v));
2019: if (lda > m) {
2020: for (j = 0; j < n; j++) PetscCall(PetscArrayzero(v + j * lda, m));
2021: } else {
2022: PetscCall(PetscArrayzero(v, PetscInt64Mult(m, n)));
2023: }
2024: PetscCall(MatDenseRestoreArrayWrite(A, &v));
2025: PetscFunctionReturn(PETSC_SUCCESS);
2026: }
2028: static PetscErrorCode MatZeroRows_SeqDense(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2029: {
2030: Mat_SeqDense *l = (Mat_SeqDense *)A->data;
2031: PetscInt m = l->lda, n = A->cmap->n, i, j;
2032: PetscScalar *slot, *bb, *v;
2033: const PetscScalar *xx;
2035: PetscFunctionBegin;
2036: if (PetscDefined(USE_DEBUG)) {
2037: for (i = 0; i < N; i++) {
2038: PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed");
2039: 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);
2040: }
2041: }
2042: if (!N) PetscFunctionReturn(PETSC_SUCCESS);
2044: /* fix right hand side if needed */
2045: if (x && b) {
2046: PetscCall(VecGetArrayRead(x, &xx));
2047: PetscCall(VecGetArray(b, &bb));
2048: for (i = 0; i < N; i++) bb[rows[i]] = diag * xx[rows[i]];
2049: PetscCall(VecRestoreArrayRead(x, &xx));
2050: PetscCall(VecRestoreArray(b, &bb));
2051: }
2053: PetscCall(MatDenseGetArray(A, &v));
2054: for (i = 0; i < N; i++) {
2055: slot = v + rows[i];
2056: for (j = 0; j < n; j++) {
2057: *slot = 0.0;
2058: slot += m;
2059: }
2060: }
2061: if (diag != 0.0) {
2062: PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices");
2063: for (i = 0; i < N; i++) {
2064: slot = v + (m + 1) * rows[i];
2065: *slot = diag;
2066: }
2067: }
2068: PetscCall(MatDenseRestoreArray(A, &v));
2069: PetscFunctionReturn(PETSC_SUCCESS);
2070: }
2072: static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A, PetscInt *lda)
2073: {
2074: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2076: PetscFunctionBegin;
2077: *lda = mat->lda;
2078: PetscFunctionReturn(PETSC_SUCCESS);
2079: }
2081: PetscErrorCode MatDenseGetArray_SeqDense(Mat A, PetscScalar **array)
2082: {
2083: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2085: PetscFunctionBegin;
2086: PetscCheck(!mat->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
2087: *array = mat->v;
2088: PetscFunctionReturn(PETSC_SUCCESS);
2089: }
2091: PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A, PetscScalar **array)
2092: {
2093: PetscFunctionBegin;
2094: if (array) *array = NULL;
2095: PetscFunctionReturn(PETSC_SUCCESS);
2096: }
2098: /*@
2099: MatDenseGetLDA - gets the leading dimension of the array returned from `MatDenseGetArray()`
2101: Not Collective
2103: Input Parameter:
2104: . A - a `MATDENSE` or `MATDENSECUDA` matrix
2106: Output Parameter:
2107: . lda - the leading dimension
2109: Level: intermediate
2111: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseSetLDA()`
2112: @*/
2113: PetscErrorCode MatDenseGetLDA(Mat A, PetscInt *lda)
2114: {
2115: PetscFunctionBegin;
2117: PetscAssertPointer(lda, 2);
2118: MatCheckPreallocated(A, 1);
2119: PetscUseMethod(A, "MatDenseGetLDA_C", (Mat, PetscInt *), (A, lda));
2120: PetscFunctionReturn(PETSC_SUCCESS);
2121: }
2123: /*@
2124: MatDenseSetLDA - Sets the leading dimension of the array used by the `MATDENSE` matrix
2126: Collective if the matrix layouts have not yet been setup
2128: Input Parameters:
2129: + A - a `MATDENSE` or `MATDENSECUDA` matrix
2130: - lda - the leading dimension
2132: Level: intermediate
2134: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetLDA()`
2135: @*/
2136: PetscErrorCode MatDenseSetLDA(Mat A, PetscInt lda)
2137: {
2138: PetscFunctionBegin;
2140: PetscTryMethod(A, "MatDenseSetLDA_C", (Mat, PetscInt), (A, lda));
2141: PetscFunctionReturn(PETSC_SUCCESS);
2142: }
2144: /*@C
2145: MatDenseGetArray - gives read-write access to the array where the data for a `MATDENSE` matrix is stored
2147: Logically Collective
2149: Input Parameter:
2150: . A - a dense matrix
2152: Output Parameter:
2153: . array - pointer to the data
2155: Level: intermediate
2157: Fortran Notes:
2158: `MatDenseGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `MatDenseGetArrayF90()`
2160: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2161: @*/
2162: PetscErrorCode MatDenseGetArray(Mat A, PetscScalar **array)
2163: {
2164: PetscFunctionBegin;
2166: PetscAssertPointer(array, 2);
2167: PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2168: PetscFunctionReturn(PETSC_SUCCESS);
2169: }
2171: /*@C
2172: MatDenseRestoreArray - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArray()`
2174: Logically Collective
2176: Input Parameters:
2177: + A - a dense matrix
2178: - array - pointer to the data (may be `NULL`)
2180: Level: intermediate
2182: Fortran Notes:
2183: `MatDenseRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `MatDenseRestoreArrayF90()`
2185: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2186: @*/
2187: PetscErrorCode MatDenseRestoreArray(Mat A, PetscScalar **array)
2188: {
2189: PetscFunctionBegin;
2191: if (array) PetscAssertPointer(array, 2);
2192: PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2193: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2194: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2195: A->offloadmask = PETSC_OFFLOAD_CPU;
2196: #endif
2197: PetscFunctionReturn(PETSC_SUCCESS);
2198: }
2200: /*@C
2201: MatDenseGetArrayRead - gives read-only access to the array where the data for a `MATDENSE` matrix is stored
2203: Not Collective
2205: Input Parameter:
2206: . A - a dense matrix
2208: Output Parameter:
2209: . array - pointer to the data
2211: Level: intermediate
2213: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2214: @*/
2215: PetscErrorCode MatDenseGetArrayRead(Mat A, const PetscScalar **array)
2216: {
2217: PetscFunctionBegin;
2219: PetscAssertPointer(array, 2);
2220: PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2221: PetscFunctionReturn(PETSC_SUCCESS);
2222: }
2224: /*@C
2225: MatDenseRestoreArrayRead - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayRead()`
2227: Not Collective
2229: Input Parameters:
2230: + A - a dense matrix
2231: - array - pointer to the data (may be `NULL`)
2233: Level: intermediate
2235: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2236: @*/
2237: PetscErrorCode MatDenseRestoreArrayRead(Mat A, const PetscScalar **array)
2238: {
2239: PetscFunctionBegin;
2241: if (array) PetscAssertPointer(array, 2);
2242: PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2243: PetscFunctionReturn(PETSC_SUCCESS);
2244: }
2246: /*@C
2247: MatDenseGetArrayWrite - gives write-only access to the array where the data for a `MATDENSE` matrix is stored
2249: Not Collective
2251: Input Parameter:
2252: . A - a dense matrix
2254: Output Parameter:
2255: . array - pointer to the data
2257: Level: intermediate
2259: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2260: @*/
2261: PetscErrorCode MatDenseGetArrayWrite(Mat A, PetscScalar **array)
2262: {
2263: PetscFunctionBegin;
2265: PetscAssertPointer(array, 2);
2266: PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2267: PetscFunctionReturn(PETSC_SUCCESS);
2268: }
2270: /*@C
2271: MatDenseRestoreArrayWrite - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayWrite()`
2273: Not Collective
2275: Input Parameters:
2276: + A - a dense matrix
2277: - array - pointer to the data (may be `NULL`)
2279: Level: intermediate
2281: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`
2282: @*/
2283: PetscErrorCode MatDenseRestoreArrayWrite(Mat A, PetscScalar **array)
2284: {
2285: PetscFunctionBegin;
2287: if (array) PetscAssertPointer(array, 2);
2288: PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2289: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2290: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
2291: A->offloadmask = PETSC_OFFLOAD_CPU;
2292: #endif
2293: PetscFunctionReturn(PETSC_SUCCESS);
2294: }
2296: /*@C
2297: MatDenseGetArrayAndMemType - gives read-write access to the array where the data for a `MATDENSE` matrix is stored
2299: Logically Collective
2301: Input Parameter:
2302: . A - a dense matrix
2304: Output Parameters:
2305: + array - pointer to the data
2306: - mtype - memory type of the returned pointer
2308: Level: intermediate
2310: Note:
2311: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2312: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2314: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArrayRead()`,
2315: `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2316: @*/
2317: PetscErrorCode MatDenseGetArrayAndMemType(Mat A, PetscScalar **array, PetscMemType *mtype)
2318: {
2319: PetscBool isMPI;
2321: PetscFunctionBegin;
2323: PetscAssertPointer(array, 2);
2324: 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 */
2325: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2326: if (isMPI) {
2327: /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2328: PetscCall(MatDenseGetArrayAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2329: } else {
2330: PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);
2332: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayAndMemType_C", &fptr));
2333: if (fptr) {
2334: PetscCall((*fptr)(A, array, mtype));
2335: } else {
2336: PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array));
2337: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2338: }
2339: }
2340: PetscFunctionReturn(PETSC_SUCCESS);
2341: }
2343: /*@C
2344: MatDenseRestoreArrayAndMemType - returns access to the array that is obtained by `MatDenseGetArrayAndMemType()`
2346: Logically Collective
2348: Input Parameters:
2349: + A - a dense matrix
2350: - array - pointer to the data
2352: Level: intermediate
2354: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2355: @*/
2356: PetscErrorCode MatDenseRestoreArrayAndMemType(Mat A, PetscScalar **array)
2357: {
2358: PetscBool isMPI;
2360: PetscFunctionBegin;
2362: PetscAssertPointer(array, 2);
2363: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2364: if (isMPI) {
2365: PetscCall(MatDenseRestoreArrayAndMemType(((Mat_MPIDense *)A->data)->A, array));
2366: } else {
2367: PetscErrorCode (*fptr)(Mat, PetscScalar **);
2369: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayAndMemType_C", &fptr));
2370: if (fptr) {
2371: PetscCall((*fptr)(A, array));
2372: } else {
2373: PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array));
2374: }
2375: *array = NULL;
2376: }
2377: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2378: PetscFunctionReturn(PETSC_SUCCESS);
2379: }
2381: /*@C
2382: MatDenseGetArrayReadAndMemType - gives read-only access to the array where the data for a `MATDENSE` matrix is stored
2384: Logically Collective
2386: Input Parameter:
2387: . A - a dense matrix
2389: Output Parameters:
2390: + array - pointer to the data
2391: - mtype - memory type of the returned pointer
2393: Level: intermediate
2395: Note:
2396: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2397: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2399: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`,
2400: `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2401: @*/
2402: PetscErrorCode MatDenseGetArrayReadAndMemType(Mat A, const PetscScalar **array, PetscMemType *mtype)
2403: {
2404: PetscBool isMPI;
2406: PetscFunctionBegin;
2408: PetscAssertPointer(array, 2);
2409: 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 */
2410: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2411: if (isMPI) { /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */
2412: PetscCall(MatDenseGetArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2413: } else {
2414: PetscErrorCode (*fptr)(Mat, const PetscScalar **, PetscMemType *);
2416: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayReadAndMemType_C", &fptr));
2417: if (fptr) {
2418: PetscCall((*fptr)(A, array, mtype));
2419: } else {
2420: PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2421: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2422: }
2423: }
2424: PetscFunctionReturn(PETSC_SUCCESS);
2425: }
2427: /*@C
2428: MatDenseRestoreArrayReadAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()`
2430: Logically Collective
2432: Input Parameters:
2433: + A - a dense matrix
2434: - array - pointer to the data
2436: Level: intermediate
2438: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2439: @*/
2440: PetscErrorCode MatDenseRestoreArrayReadAndMemType(Mat A, const PetscScalar **array)
2441: {
2442: PetscBool isMPI;
2444: PetscFunctionBegin;
2446: PetscAssertPointer(array, 2);
2447: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2448: if (isMPI) {
2449: PetscCall(MatDenseRestoreArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array));
2450: } else {
2451: PetscErrorCode (*fptr)(Mat, const PetscScalar **);
2453: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayReadAndMemType_C", &fptr));
2454: if (fptr) {
2455: PetscCall((*fptr)(A, array));
2456: } else {
2457: PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, PetscScalar **), (A, (PetscScalar **)array));
2458: }
2459: *array = NULL;
2460: }
2461: PetscFunctionReturn(PETSC_SUCCESS);
2462: }
2464: /*@C
2465: MatDenseGetArrayWriteAndMemType - gives write-only access to the array where the data for a `MATDENSE` matrix is stored
2467: Logically Collective
2469: Input Parameter:
2470: . A - a dense matrix
2472: Output Parameters:
2473: + array - pointer to the data
2474: - mtype - memory type of the returned pointer
2476: Level: intermediate
2478: Note:
2479: If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc.,
2480: an array on device is always returned and is guaranteed to contain the matrix's latest data.
2482: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWriteAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayRead()`,
2483: `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()`
2484: @*/
2485: PetscErrorCode MatDenseGetArrayWriteAndMemType(Mat A, PetscScalar **array, PetscMemType *mtype)
2486: {
2487: PetscBool isMPI;
2489: PetscFunctionBegin;
2491: PetscAssertPointer(array, 2);
2492: 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 */
2493: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2494: if (isMPI) {
2495: PetscCall(MatDenseGetArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype));
2496: } else {
2497: PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *);
2499: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayWriteAndMemType_C", &fptr));
2500: if (fptr) {
2501: PetscCall((*fptr)(A, array, mtype));
2502: } else {
2503: PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array));
2504: if (mtype) *mtype = PETSC_MEMTYPE_HOST;
2505: }
2506: }
2507: PetscFunctionReturn(PETSC_SUCCESS);
2508: }
2510: /*@C
2511: MatDenseRestoreArrayWriteAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()`
2513: Logically Collective
2515: Input Parameters:
2516: + A - a dense matrix
2517: - array - pointer to the data
2519: Level: intermediate
2521: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`
2522: @*/
2523: PetscErrorCode MatDenseRestoreArrayWriteAndMemType(Mat A, PetscScalar **array)
2524: {
2525: PetscBool isMPI;
2527: PetscFunctionBegin;
2529: PetscAssertPointer(array, 2);
2530: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI));
2531: if (isMPI) {
2532: PetscCall(MatDenseRestoreArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array));
2533: } else {
2534: PetscErrorCode (*fptr)(Mat, PetscScalar **);
2536: PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayWriteAndMemType_C", &fptr));
2537: if (fptr) {
2538: PetscCall((*fptr)(A, array));
2539: } else {
2540: PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array));
2541: }
2542: *array = NULL;
2543: }
2544: PetscCall(PetscObjectStateIncrease((PetscObject)A));
2545: PetscFunctionReturn(PETSC_SUCCESS);
2546: }
2548: static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B)
2549: {
2550: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2551: PetscInt i, j, nrows, ncols, ldb;
2552: const PetscInt *irow, *icol;
2553: PetscScalar *av, *bv, *v = mat->v;
2554: Mat newmat;
2556: PetscFunctionBegin;
2557: PetscCall(ISGetIndices(isrow, &irow));
2558: PetscCall(ISGetIndices(iscol, &icol));
2559: PetscCall(ISGetLocalSize(isrow, &nrows));
2560: PetscCall(ISGetLocalSize(iscol, &ncols));
2562: /* Check submatrixcall */
2563: if (scall == MAT_REUSE_MATRIX) {
2564: PetscInt n_cols, n_rows;
2565: PetscCall(MatGetSize(*B, &n_rows, &n_cols));
2566: if (n_rows != nrows || n_cols != ncols) {
2567: /* resize the result matrix to match number of requested rows/columns */
2568: PetscCall(MatSetSizes(*B, nrows, ncols, nrows, ncols));
2569: }
2570: newmat = *B;
2571: } else {
2572: /* Create and fill new matrix */
2573: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &newmat));
2574: PetscCall(MatSetSizes(newmat, nrows, ncols, nrows, ncols));
2575: PetscCall(MatSetType(newmat, ((PetscObject)A)->type_name));
2576: PetscCall(MatSeqDenseSetPreallocation(newmat, NULL));
2577: }
2579: /* Now extract the data pointers and do the copy,column at a time */
2580: PetscCall(MatDenseGetArray(newmat, &bv));
2581: PetscCall(MatDenseGetLDA(newmat, &ldb));
2582: for (i = 0; i < ncols; i++) {
2583: av = v + mat->lda * icol[i];
2584: for (j = 0; j < nrows; j++) bv[j] = av[irow[j]];
2585: bv += ldb;
2586: }
2587: PetscCall(MatDenseRestoreArray(newmat, &bv));
2589: /* Assemble the matrices so that the correct flags are set */
2590: PetscCall(MatAssemblyBegin(newmat, MAT_FINAL_ASSEMBLY));
2591: PetscCall(MatAssemblyEnd(newmat, MAT_FINAL_ASSEMBLY));
2593: /* Free work space */
2594: PetscCall(ISRestoreIndices(isrow, &irow));
2595: PetscCall(ISRestoreIndices(iscol, &icol));
2596: *B = newmat;
2597: PetscFunctionReturn(PETSC_SUCCESS);
2598: }
2600: static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[])
2601: {
2602: PetscInt i;
2604: PetscFunctionBegin;
2605: if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n, B));
2607: for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SeqDense(A, irow[i], icol[i], scall, &(*B)[i]));
2608: PetscFunctionReturn(PETSC_SUCCESS);
2609: }
2611: static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat, MatAssemblyType mode)
2612: {
2613: PetscFunctionBegin;
2614: PetscFunctionReturn(PETSC_SUCCESS);
2615: }
2617: static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat, MatAssemblyType mode)
2618: {
2619: PetscFunctionBegin;
2620: PetscFunctionReturn(PETSC_SUCCESS);
2621: }
2623: PetscErrorCode MatCopy_SeqDense(Mat A, Mat B, MatStructure str)
2624: {
2625: Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data;
2626: const PetscScalar *va;
2627: PetscScalar *vb;
2628: PetscInt lda1 = a->lda, lda2 = b->lda, m = A->rmap->n, n = A->cmap->n, j;
2630: PetscFunctionBegin;
2631: /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */
2632: if (A->ops->copy != B->ops->copy) {
2633: PetscCall(MatCopy_Basic(A, B, str));
2634: PetscFunctionReturn(PETSC_SUCCESS);
2635: }
2636: PetscCheck(m == B->rmap->n && n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "size(B) != size(A)");
2637: PetscCall(MatDenseGetArrayRead(A, &va));
2638: PetscCall(MatDenseGetArray(B, &vb));
2639: if (lda1 > m || lda2 > m) {
2640: for (j = 0; j < n; j++) PetscCall(PetscArraycpy(vb + j * lda2, va + j * lda1, m));
2641: } else {
2642: PetscCall(PetscArraycpy(vb, va, A->rmap->n * A->cmap->n));
2643: }
2644: PetscCall(MatDenseRestoreArray(B, &vb));
2645: PetscCall(MatDenseRestoreArrayRead(A, &va));
2646: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2647: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2648: PetscFunctionReturn(PETSC_SUCCESS);
2649: }
2651: PetscErrorCode MatSetUp_SeqDense(Mat A)
2652: {
2653: PetscFunctionBegin;
2654: PetscCall(PetscLayoutSetUp(A->rmap));
2655: PetscCall(PetscLayoutSetUp(A->cmap));
2656: if (!A->preallocated) PetscCall(MatSeqDenseSetPreallocation(A, NULL));
2657: PetscFunctionReturn(PETSC_SUCCESS);
2658: }
2660: static PetscErrorCode MatConjugate_SeqDense(Mat A)
2661: {
2662: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2663: PetscInt i, j;
2664: PetscInt min = PetscMin(A->rmap->n, A->cmap->n);
2665: PetscScalar *aa;
2667: PetscFunctionBegin;
2668: PetscCall(MatDenseGetArray(A, &aa));
2669: for (j = 0; j < A->cmap->n; j++) {
2670: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscConj(aa[i + j * mat->lda]);
2671: }
2672: PetscCall(MatDenseRestoreArray(A, &aa));
2673: if (mat->tau)
2674: for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]);
2675: PetscFunctionReturn(PETSC_SUCCESS);
2676: }
2678: static PetscErrorCode MatRealPart_SeqDense(Mat A)
2679: {
2680: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2681: PetscInt i, j;
2682: PetscScalar *aa;
2684: PetscFunctionBegin;
2685: PetscCall(MatDenseGetArray(A, &aa));
2686: for (j = 0; j < A->cmap->n; j++) {
2687: for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscRealPart(aa[i + j * mat->lda]);
2688: }
2689: PetscCall(MatDenseRestoreArray(A, &aa));
2690: PetscFunctionReturn(PETSC_SUCCESS);
2691: }
2693: static PetscErrorCode MatImaginaryPart_SeqDense(Mat A)
2694: {
2695: Mat_SeqDense *mat = (Mat_SeqDense *)A->data;
2696: PetscInt i, j;
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] = PetscImaginaryPart(aa[i + j * mat->lda]);
2703: }
2704: PetscCall(MatDenseRestoreArray(A, &aa));
2705: PetscFunctionReturn(PETSC_SUCCESS);
2706: }
2708: PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2709: {
2710: PetscInt m = A->rmap->n, n = B->cmap->n;
2711: PetscBool cisdense = PETSC_FALSE;
2713: PetscFunctionBegin;
2714: PetscCall(MatSetSizes(C, m, n, m, n));
2715: #if defined(PETSC_HAVE_CUDA)
2716: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2717: #endif
2718: #if defined(PETSC_HAVE_HIP)
2719: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2720: #endif
2721: if (!cisdense) {
2722: PetscBool flg;
2724: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2725: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2726: }
2727: PetscCall(MatSetUp(C));
2728: PetscFunctionReturn(PETSC_SUCCESS);
2729: }
2731: PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2732: {
2733: Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data, *c = (Mat_SeqDense *)C->data;
2734: PetscBLASInt m, n, k;
2735: const PetscScalar *av, *bv;
2736: PetscScalar *cv;
2737: PetscScalar _DOne = 1.0, _DZero = 0.0;
2739: PetscFunctionBegin;
2740: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2741: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2742: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2743: if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2744: PetscCall(MatDenseGetArrayRead(A, &av));
2745: PetscCall(MatDenseGetArrayRead(B, &bv));
2746: PetscCall(MatDenseGetArrayWrite(C, &cv));
2747: PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2748: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2749: PetscCall(MatDenseRestoreArrayRead(A, &av));
2750: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2751: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2752: PetscFunctionReturn(PETSC_SUCCESS);
2753: }
2755: PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2756: {
2757: PetscInt m = A->rmap->n, n = B->rmap->n;
2758: PetscBool cisdense = PETSC_FALSE;
2760: PetscFunctionBegin;
2761: PetscCall(MatSetSizes(C, m, n, m, n));
2762: #if defined(PETSC_HAVE_CUDA)
2763: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2764: #endif
2765: #if defined(PETSC_HAVE_HIP)
2766: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2767: #endif
2768: if (!cisdense) {
2769: PetscBool flg;
2771: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2772: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2773: }
2774: PetscCall(MatSetUp(C));
2775: PetscFunctionReturn(PETSC_SUCCESS);
2776: }
2778: PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2779: {
2780: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2781: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
2782: Mat_SeqDense *c = (Mat_SeqDense *)C->data;
2783: const PetscScalar *av, *bv;
2784: PetscScalar *cv;
2785: PetscBLASInt m, n, k;
2786: PetscScalar _DOne = 1.0, _DZero = 0.0;
2788: PetscFunctionBegin;
2789: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2790: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2791: PetscCall(PetscBLASIntCast(A->cmap->n, &k));
2792: if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2793: PetscCall(MatDenseGetArrayRead(A, &av));
2794: PetscCall(MatDenseGetArrayRead(B, &bv));
2795: PetscCall(MatDenseGetArrayWrite(C, &cv));
2796: PetscCallBLAS("BLASgemm", BLASgemm_("N", "T", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2797: PetscCall(MatDenseRestoreArrayRead(A, &av));
2798: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2799: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2800: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2801: PetscFunctionReturn(PETSC_SUCCESS);
2802: }
2804: PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
2805: {
2806: PetscInt m = A->cmap->n, n = B->cmap->n;
2807: PetscBool cisdense = PETSC_FALSE;
2809: PetscFunctionBegin;
2810: PetscCall(MatSetSizes(C, m, n, m, n));
2811: #if defined(PETSC_HAVE_CUDA)
2812: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, ""));
2813: #endif
2814: #if defined(PETSC_HAVE_HIP)
2815: PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, ""));
2816: #endif
2817: if (!cisdense) {
2818: PetscBool flg;
2820: PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg));
2821: PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE));
2822: }
2823: PetscCall(MatSetUp(C));
2824: PetscFunctionReturn(PETSC_SUCCESS);
2825: }
2827: PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C)
2828: {
2829: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2830: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
2831: Mat_SeqDense *c = (Mat_SeqDense *)C->data;
2832: const PetscScalar *av, *bv;
2833: PetscScalar *cv;
2834: PetscBLASInt m, n, k;
2835: PetscScalar _DOne = 1.0, _DZero = 0.0;
2837: PetscFunctionBegin;
2838: PetscCall(PetscBLASIntCast(C->rmap->n, &m));
2839: PetscCall(PetscBLASIntCast(C->cmap->n, &n));
2840: PetscCall(PetscBLASIntCast(A->rmap->n, &k));
2841: if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS);
2842: PetscCall(MatDenseGetArrayRead(A, &av));
2843: PetscCall(MatDenseGetArrayRead(B, &bv));
2844: PetscCall(MatDenseGetArrayWrite(C, &cv));
2845: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda));
2846: PetscCall(MatDenseRestoreArrayRead(A, &av));
2847: PetscCall(MatDenseRestoreArrayRead(B, &bv));
2848: PetscCall(MatDenseRestoreArrayWrite(C, &cv));
2849: PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1)));
2850: PetscFunctionReturn(PETSC_SUCCESS);
2851: }
2853: static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C)
2854: {
2855: PetscFunctionBegin;
2856: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense;
2857: C->ops->productsymbolic = MatProductSymbolic_AB;
2858: PetscFunctionReturn(PETSC_SUCCESS);
2859: }
2861: static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C)
2862: {
2863: PetscFunctionBegin;
2864: C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense;
2865: C->ops->productsymbolic = MatProductSymbolic_AtB;
2866: PetscFunctionReturn(PETSC_SUCCESS);
2867: }
2869: static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C)
2870: {
2871: PetscFunctionBegin;
2872: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense;
2873: C->ops->productsymbolic = MatProductSymbolic_ABt;
2874: PetscFunctionReturn(PETSC_SUCCESS);
2875: }
2877: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C)
2878: {
2879: Mat_Product *product = C->product;
2881: PetscFunctionBegin;
2882: switch (product->type) {
2883: case MATPRODUCT_AB:
2884: PetscCall(MatProductSetFromOptions_SeqDense_AB(C));
2885: break;
2886: case MATPRODUCT_AtB:
2887: PetscCall(MatProductSetFromOptions_SeqDense_AtB(C));
2888: break;
2889: case MATPRODUCT_ABt:
2890: PetscCall(MatProductSetFromOptions_SeqDense_ABt(C));
2891: break;
2892: default:
2893: break;
2894: }
2895: PetscFunctionReturn(PETSC_SUCCESS);
2896: }
2898: static PetscErrorCode MatGetRowMax_SeqDense(Mat A, Vec v, PetscInt idx[])
2899: {
2900: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2901: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2902: PetscScalar *x;
2903: const PetscScalar *aa;
2905: PetscFunctionBegin;
2906: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2907: PetscCall(VecGetArray(v, &x));
2908: PetscCall(VecGetLocalSize(v, &p));
2909: PetscCall(MatDenseGetArrayRead(A, &aa));
2910: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2911: for (i = 0; i < m; i++) {
2912: x[i] = aa[i];
2913: if (idx) idx[i] = 0;
2914: for (j = 1; j < n; j++) {
2915: if (PetscRealPart(x[i]) < PetscRealPart(aa[i + a->lda * j])) {
2916: x[i] = aa[i + a->lda * j];
2917: if (idx) idx[i] = j;
2918: }
2919: }
2920: }
2921: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2922: PetscCall(VecRestoreArray(v, &x));
2923: PetscFunctionReturn(PETSC_SUCCESS);
2924: }
2926: static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A, Vec v, PetscInt idx[])
2927: {
2928: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2929: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2930: PetscScalar *x;
2931: PetscReal atmp;
2932: const PetscScalar *aa;
2934: PetscFunctionBegin;
2935: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2936: PetscCall(VecGetArray(v, &x));
2937: PetscCall(VecGetLocalSize(v, &p));
2938: PetscCall(MatDenseGetArrayRead(A, &aa));
2939: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2940: for (i = 0; i < m; i++) {
2941: x[i] = PetscAbsScalar(aa[i]);
2942: for (j = 1; j < n; j++) {
2943: atmp = PetscAbsScalar(aa[i + a->lda * j]);
2944: if (PetscAbsScalar(x[i]) < atmp) {
2945: x[i] = atmp;
2946: if (idx) idx[i] = j;
2947: }
2948: }
2949: }
2950: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2951: PetscCall(VecRestoreArray(v, &x));
2952: PetscFunctionReturn(PETSC_SUCCESS);
2953: }
2955: static PetscErrorCode MatGetRowMin_SeqDense(Mat A, Vec v, PetscInt idx[])
2956: {
2957: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2958: PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p;
2959: PetscScalar *x;
2960: const PetscScalar *aa;
2962: PetscFunctionBegin;
2963: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2964: PetscCall(MatDenseGetArrayRead(A, &aa));
2965: PetscCall(VecGetArray(v, &x));
2966: PetscCall(VecGetLocalSize(v, &p));
2967: PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2968: for (i = 0; i < m; i++) {
2969: x[i] = aa[i];
2970: if (idx) idx[i] = 0;
2971: for (j = 1; j < n; j++) {
2972: if (PetscRealPart(x[i]) > PetscRealPart(aa[i + a->lda * j])) {
2973: x[i] = aa[i + a->lda * j];
2974: if (idx) idx[i] = j;
2975: }
2976: }
2977: }
2978: PetscCall(VecRestoreArray(v, &x));
2979: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2980: PetscFunctionReturn(PETSC_SUCCESS);
2981: }
2983: PetscErrorCode MatGetColumnVector_SeqDense(Mat A, Vec v, PetscInt col)
2984: {
2985: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
2986: PetscScalar *x;
2987: const PetscScalar *aa;
2989: PetscFunctionBegin;
2990: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2991: PetscCall(MatDenseGetArrayRead(A, &aa));
2992: PetscCall(VecGetArray(v, &x));
2993: PetscCall(PetscArraycpy(x, aa + col * a->lda, A->rmap->n));
2994: PetscCall(VecRestoreArray(v, &x));
2995: PetscCall(MatDenseRestoreArrayRead(A, &aa));
2996: PetscFunctionReturn(PETSC_SUCCESS);
2997: }
2999: PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A, PetscInt type, PetscReal *reductions)
3000: {
3001: PetscInt i, j, m, n;
3002: const PetscScalar *a;
3004: PetscFunctionBegin;
3005: PetscCall(MatGetSize(A, &m, &n));
3006: PetscCall(PetscArrayzero(reductions, n));
3007: PetscCall(MatDenseGetArrayRead(A, &a));
3008: if (type == NORM_2) {
3009: for (i = 0; i < n; i++) {
3010: for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j] * a[j]);
3011: a += m;
3012: }
3013: } else if (type == NORM_1) {
3014: for (i = 0; i < n; i++) {
3015: for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j]);
3016: a += m;
3017: }
3018: } else if (type == NORM_INFINITY) {
3019: for (i = 0; i < n; i++) {
3020: for (j = 0; j < m; j++) reductions[i] = PetscMax(PetscAbsScalar(a[j]), reductions[i]);
3021: a += m;
3022: }
3023: } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
3024: for (i = 0; i < n; i++) {
3025: for (j = 0; j < m; j++) reductions[i] += PetscRealPart(a[j]);
3026: a += m;
3027: }
3028: } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3029: for (i = 0; i < n; i++) {
3030: for (j = 0; j < m; j++) reductions[i] += PetscImaginaryPart(a[j]);
3031: a += m;
3032: }
3033: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown reduction type");
3034: PetscCall(MatDenseRestoreArrayRead(A, &a));
3035: if (type == NORM_2) {
3036: for (i = 0; i < n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
3037: } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
3038: for (i = 0; i < n; i++) reductions[i] /= m;
3039: }
3040: PetscFunctionReturn(PETSC_SUCCESS);
3041: }
3043: PetscErrorCode MatSetRandom_SeqDense(Mat x, PetscRandom rctx)
3044: {
3045: PetscScalar *a;
3046: PetscInt lda, m, n, i, j;
3048: PetscFunctionBegin;
3049: PetscCall(MatGetSize(x, &m, &n));
3050: PetscCall(MatDenseGetLDA(x, &lda));
3051: PetscCall(MatDenseGetArrayWrite(x, &a));
3052: for (j = 0; j < n; j++) {
3053: for (i = 0; i < m; i++) PetscCall(PetscRandomGetValue(rctx, a + j * lda + i));
3054: }
3055: PetscCall(MatDenseRestoreArrayWrite(x, &a));
3056: PetscFunctionReturn(PETSC_SUCCESS);
3057: }
3059: static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A, PetscBool *missing, PetscInt *d)
3060: {
3061: PetscFunctionBegin;
3062: *missing = PETSC_FALSE;
3063: PetscFunctionReturn(PETSC_SUCCESS);
3064: }
3066: /* vals is not const */
3067: static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A, PetscInt col, PetscScalar **vals)
3068: {
3069: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3070: PetscScalar *v;
3072: PetscFunctionBegin;
3073: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3074: PetscCall(MatDenseGetArray(A, &v));
3075: *vals = v + col * a->lda;
3076: PetscCall(MatDenseRestoreArray(A, &v));
3077: PetscFunctionReturn(PETSC_SUCCESS);
3078: }
3080: static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A, PetscScalar **vals)
3081: {
3082: PetscFunctionBegin;
3083: if (vals) *vals = NULL; /* user cannot accidentally use the array later */
3084: PetscFunctionReturn(PETSC_SUCCESS);
3085: }
3087: static struct _MatOps MatOps_Values = {MatSetValues_SeqDense,
3088: MatGetRow_SeqDense,
3089: MatRestoreRow_SeqDense,
3090: MatMult_SeqDense,
3091: /* 4*/ MatMultAdd_SeqDense,
3092: MatMultTranspose_SeqDense,
3093: MatMultTransposeAdd_SeqDense,
3094: NULL,
3095: NULL,
3096: NULL,
3097: /* 10*/ NULL,
3098: MatLUFactor_SeqDense,
3099: MatCholeskyFactor_SeqDense,
3100: MatSOR_SeqDense,
3101: MatTranspose_SeqDense,
3102: /* 15*/ MatGetInfo_SeqDense,
3103: MatEqual_SeqDense,
3104: MatGetDiagonal_SeqDense,
3105: MatDiagonalScale_SeqDense,
3106: MatNorm_SeqDense,
3107: /* 20*/ MatAssemblyBegin_SeqDense,
3108: MatAssemblyEnd_SeqDense,
3109: MatSetOption_SeqDense,
3110: MatZeroEntries_SeqDense,
3111: /* 24*/ MatZeroRows_SeqDense,
3112: NULL,
3113: NULL,
3114: NULL,
3115: NULL,
3116: /* 29*/ MatSetUp_SeqDense,
3117: NULL,
3118: NULL,
3119: NULL,
3120: NULL,
3121: /* 34*/ MatDuplicate_SeqDense,
3122: NULL,
3123: NULL,
3124: NULL,
3125: NULL,
3126: /* 39*/ MatAXPY_SeqDense,
3127: MatCreateSubMatrices_SeqDense,
3128: NULL,
3129: MatGetValues_SeqDense,
3130: MatCopy_SeqDense,
3131: /* 44*/ MatGetRowMax_SeqDense,
3132: MatScale_SeqDense,
3133: MatShift_SeqDense,
3134: NULL,
3135: MatZeroRowsColumns_SeqDense,
3136: /* 49*/ MatSetRandom_SeqDense,
3137: NULL,
3138: NULL,
3139: NULL,
3140: NULL,
3141: /* 54*/ NULL,
3142: NULL,
3143: NULL,
3144: NULL,
3145: NULL,
3146: /* 59*/ MatCreateSubMatrix_SeqDense,
3147: MatDestroy_SeqDense,
3148: MatView_SeqDense,
3149: NULL,
3150: NULL,
3151: /* 64*/ NULL,
3152: NULL,
3153: NULL,
3154: NULL,
3155: NULL,
3156: /* 69*/ MatGetRowMaxAbs_SeqDense,
3157: NULL,
3158: NULL,
3159: NULL,
3160: NULL,
3161: /* 74*/ NULL,
3162: NULL,
3163: NULL,
3164: NULL,
3165: NULL,
3166: /* 79*/ NULL,
3167: NULL,
3168: NULL,
3169: NULL,
3170: /* 83*/ MatLoad_SeqDense,
3171: MatIsSymmetric_SeqDense,
3172: MatIsHermitian_SeqDense,
3173: NULL,
3174: NULL,
3175: NULL,
3176: /* 89*/ NULL,
3177: NULL,
3178: MatMatMultNumeric_SeqDense_SeqDense,
3179: NULL,
3180: NULL,
3181: /* 94*/ NULL,
3182: NULL,
3183: NULL,
3184: MatMatTransposeMultNumeric_SeqDense_SeqDense,
3185: NULL,
3186: /* 99*/ MatProductSetFromOptions_SeqDense,
3187: NULL,
3188: NULL,
3189: MatConjugate_SeqDense,
3190: NULL,
3191: /*104*/ NULL,
3192: MatRealPart_SeqDense,
3193: MatImaginaryPart_SeqDense,
3194: NULL,
3195: NULL,
3196: /*109*/ NULL,
3197: NULL,
3198: MatGetRowMin_SeqDense,
3199: MatGetColumnVector_SeqDense,
3200: MatMissingDiagonal_SeqDense,
3201: /*114*/ NULL,
3202: NULL,
3203: NULL,
3204: NULL,
3205: NULL,
3206: /*119*/ NULL,
3207: NULL,
3208: NULL,
3209: NULL,
3210: NULL,
3211: /*124*/ NULL,
3212: MatGetColumnReductions_SeqDense,
3213: NULL,
3214: NULL,
3215: NULL,
3216: /*129*/ NULL,
3217: NULL,
3218: NULL,
3219: MatTransposeMatMultNumeric_SeqDense_SeqDense,
3220: NULL,
3221: /*134*/ NULL,
3222: NULL,
3223: NULL,
3224: NULL,
3225: NULL,
3226: /*139*/ NULL,
3227: NULL,
3228: NULL,
3229: NULL,
3230: NULL,
3231: MatCreateMPIMatConcatenateSeqMat_SeqDense,
3232: /*145*/ NULL,
3233: NULL,
3234: NULL,
3235: NULL,
3236: NULL,
3237: /*150*/ NULL,
3238: NULL};
3240: /*@C
3241: MatCreateSeqDense - Creates a `MATSEQDENSE` that
3242: is stored in column major order (the usual Fortran format).
3244: Collective
3246: Input Parameters:
3247: + comm - MPI communicator, set to `PETSC_COMM_SELF`
3248: . m - number of rows
3249: . n - number of columns
3250: - data - optional location of matrix data in column major order. Use `NULL` for PETSc
3251: to control all matrix memory allocation.
3253: Output Parameter:
3254: . A - the matrix
3256: Level: intermediate
3258: Note:
3259: The data input variable is intended primarily for Fortran programmers
3260: who wish to allocate their own matrix memory space. Most users should
3261: set `data` = `NULL`.
3263: Developer Note:
3264: Many of the matrix operations for this variant use the BLAS and LAPACK routines.
3266: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`
3267: @*/
3268: PetscErrorCode MatCreateSeqDense(MPI_Comm comm, PetscInt m, PetscInt n, PetscScalar *data, Mat *A)
3269: {
3270: PetscFunctionBegin;
3271: PetscCall(MatCreate(comm, A));
3272: PetscCall(MatSetSizes(*A, m, n, m, n));
3273: PetscCall(MatSetType(*A, MATSEQDENSE));
3274: PetscCall(MatSeqDenseSetPreallocation(*A, data));
3275: PetscFunctionReturn(PETSC_SUCCESS);
3276: }
3278: /*@C
3279: MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements of a `MATSEQDENSE` matrix
3281: Collective
3283: Input Parameters:
3284: + B - the matrix
3285: - data - the array (or `NULL`)
3287: Level: intermediate
3289: Note:
3290: The data input variable is intended primarily for Fortran programmers
3291: who wish to allocate their own matrix memory space. Most users should
3292: need not call this routine.
3294: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`, `MatDenseSetLDA()`
3295: @*/
3296: PetscErrorCode MatSeqDenseSetPreallocation(Mat B, PetscScalar data[])
3297: {
3298: PetscFunctionBegin;
3300: PetscTryMethod(B, "MatSeqDenseSetPreallocation_C", (Mat, PetscScalar[]), (B, data));
3301: PetscFunctionReturn(PETSC_SUCCESS);
3302: }
3304: PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B, PetscScalar *data)
3305: {
3306: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
3308: PetscFunctionBegin;
3309: PetscCheck(!b->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3310: B->preallocated = PETSC_TRUE;
3312: PetscCall(PetscLayoutSetUp(B->rmap));
3313: PetscCall(PetscLayoutSetUp(B->cmap));
3315: if (b->lda <= 0) b->lda = B->rmap->n;
3317: if (!data) { /* petsc-allocated storage */
3318: if (!b->user_alloc) PetscCall(PetscFree(b->v));
3319: PetscCall(PetscCalloc1((size_t)b->lda * B->cmap->n, &b->v));
3321: b->user_alloc = PETSC_FALSE;
3322: } else { /* user-allocated storage */
3323: if (!b->user_alloc) PetscCall(PetscFree(b->v));
3324: b->v = data;
3325: b->user_alloc = PETSC_TRUE;
3326: }
3327: B->assembled = PETSC_TRUE;
3328: PetscFunctionReturn(PETSC_SUCCESS);
3329: }
3331: #if defined(PETSC_HAVE_ELEMENTAL)
3332: PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
3333: {
3334: Mat mat_elemental;
3335: const PetscScalar *array;
3336: PetscScalar *v_colwise;
3337: PetscInt M = A->rmap->N, N = A->cmap->N, i, j, k, *rows, *cols;
3339: PetscFunctionBegin;
3340: PetscCall(PetscMalloc3(M * N, &v_colwise, M, &rows, N, &cols));
3341: PetscCall(MatDenseGetArrayRead(A, &array));
3342: /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */
3343: k = 0;
3344: for (j = 0; j < N; j++) {
3345: cols[j] = j;
3346: for (i = 0; i < M; i++) v_colwise[j * M + i] = array[k++];
3347: }
3348: for (i = 0; i < M; i++) rows[i] = i;
3349: PetscCall(MatDenseRestoreArrayRead(A, &array));
3351: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental));
3352: PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N));
3353: PetscCall(MatSetType(mat_elemental, MATELEMENTAL));
3354: PetscCall(MatSetUp(mat_elemental));
3356: /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */
3357: PetscCall(MatSetValues(mat_elemental, M, rows, N, cols, v_colwise, ADD_VALUES));
3358: PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY));
3359: PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY));
3360: PetscCall(PetscFree3(v_colwise, rows, cols));
3362: if (reuse == MAT_INPLACE_MATRIX) {
3363: PetscCall(MatHeaderReplace(A, &mat_elemental));
3364: } else {
3365: *newmat = mat_elemental;
3366: }
3367: PetscFunctionReturn(PETSC_SUCCESS);
3368: }
3369: #endif
3371: PetscErrorCode MatDenseSetLDA_SeqDense(Mat B, PetscInt lda)
3372: {
3373: Mat_SeqDense *b = (Mat_SeqDense *)B->data;
3374: PetscBool data;
3376: PetscFunctionBegin;
3377: data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE);
3378: PetscCheck(b->user_alloc || !data || b->lda == lda, PETSC_COMM_SELF, PETSC_ERR_ORDER, "LDA cannot be changed after allocation of internal storage");
3379: 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);
3380: b->lda = lda;
3381: PetscFunctionReturn(PETSC_SUCCESS);
3382: }
3384: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
3385: {
3386: PetscFunctionBegin;
3387: PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIDense(comm, inmat, n, scall, outmat));
3388: PetscFunctionReturn(PETSC_SUCCESS);
3389: }
3391: PetscErrorCode MatDenseCreateColumnVec_Private(Mat A, Vec *v)
3392: {
3393: PetscBool isstd, iskok, iscuda, iship;
3394: PetscMPIInt size;
3395: #if PetscDefined(HAVE_CUDA) || PetscDefined(HAVE_HIP)
3396: /* we pass the data of A, to prevent allocating needless GPU memory the first time VecCUPMPlaceArray is called. */
3397: const PetscScalar *a;
3398: #endif
3400: PetscFunctionBegin;
3401: *v = NULL;
3402: PetscCall(PetscStrcmpAny(A->defaultvectype, &isstd, VECSTANDARD, VECSEQ, VECMPI, ""));
3403: PetscCall(PetscStrcmpAny(A->defaultvectype, &iskok, VECKOKKOS, VECSEQKOKKOS, VECMPIKOKKOS, ""));
3404: PetscCall(PetscStrcmpAny(A->defaultvectype, &iscuda, VECCUDA, VECSEQCUDA, VECMPICUDA, ""));
3405: PetscCall(PetscStrcmpAny(A->defaultvectype, &iship, VECHIP, VECSEQHIP, VECMPIHIP, ""));
3406: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
3407: if (isstd) {
3408: if (size > 1) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, v));
3409: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, v));
3410: } else if (iskok) {
3411: PetscCheck(PetscDefined(HAVE_KOKKOS_KERNELS), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using KOKKOS kernels support");
3412: #if PetscDefined(HAVE_KOKKOS_KERNELS)
3413: if (size > 1) PetscCall(VecCreateMPIKokkosWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, v));
3414: else PetscCall(VecCreateSeqKokkosWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, v));
3415: #endif
3416: } else if (iscuda) {
3417: PetscCheck(PetscDefined(HAVE_CUDA), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using CUDA support");
3418: #if PetscDefined(HAVE_CUDA)
3419: PetscCall(MatDenseCUDAGetArrayRead(A, &a));
3420: if (size > 1) PetscCall(VecCreateMPICUDAWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, a, v));
3421: else PetscCall(VecCreateSeqCUDAWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, a, v));
3422: #endif
3423: } else if (iship) {
3424: PetscCheck(PetscDefined(HAVE_HIP), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Reconfigure using HIP support");
3425: #if PetscDefined(HAVE_HIP)
3426: PetscCall(MatDenseHIPGetArrayRead(A, &a));
3427: if (size > 1) PetscCall(VecCreateMPIHIPWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, A->rmap->N, NULL, a, v));
3428: else PetscCall(VecCreateSeqHIPWithArrays(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, a, v));
3429: #endif
3430: }
3431: PetscCheck(*v, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not coded for type %s", A->defaultvectype);
3432: PetscFunctionReturn(PETSC_SUCCESS);
3433: }
3435: PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3436: {
3437: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3439: PetscFunctionBegin;
3440: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3441: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3442: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3443: a->vecinuse = col + 1;
3444: PetscCall(MatDenseGetArray(A, (PetscScalar **)&a->ptrinuse));
3445: PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda));
3446: *v = a->cvec;
3447: PetscFunctionReturn(PETSC_SUCCESS);
3448: }
3450: PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A, PetscInt col, Vec *v)
3451: {
3452: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3454: PetscFunctionBegin;
3455: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3456: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3457: a->vecinuse = 0;
3458: PetscCall(MatDenseRestoreArray(A, (PetscScalar **)&a->ptrinuse));
3459: PetscCall(VecResetArray(a->cvec));
3460: if (v) *v = NULL;
3461: PetscFunctionReturn(PETSC_SUCCESS);
3462: }
3464: PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3465: {
3466: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3468: PetscFunctionBegin;
3469: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3470: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3471: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3472: a->vecinuse = col + 1;
3473: PetscCall(MatDenseGetArrayRead(A, &a->ptrinuse));
3474: PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda));
3475: PetscCall(VecLockReadPush(a->cvec));
3476: *v = a->cvec;
3477: PetscFunctionReturn(PETSC_SUCCESS);
3478: }
3480: PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v)
3481: {
3482: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3484: PetscFunctionBegin;
3485: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3486: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3487: a->vecinuse = 0;
3488: PetscCall(MatDenseRestoreArrayRead(A, &a->ptrinuse));
3489: PetscCall(VecLockReadPop(a->cvec));
3490: PetscCall(VecResetArray(a->cvec));
3491: if (v) *v = NULL;
3492: PetscFunctionReturn(PETSC_SUCCESS);
3493: }
3495: PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3496: {
3497: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3499: PetscFunctionBegin;
3500: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3501: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3502: if (!a->cvec) PetscCall(MatDenseCreateColumnVec_Private(A, &a->cvec));
3503: a->vecinuse = col + 1;
3504: PetscCall(MatDenseGetArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3505: PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda));
3506: *v = a->cvec;
3507: PetscFunctionReturn(PETSC_SUCCESS);
3508: }
3510: PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v)
3511: {
3512: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3514: PetscFunctionBegin;
3515: PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first");
3516: PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector");
3517: a->vecinuse = 0;
3518: PetscCall(MatDenseRestoreArrayWrite(A, (PetscScalar **)&a->ptrinuse));
3519: PetscCall(VecResetArray(a->cvec));
3520: if (v) *v = NULL;
3521: PetscFunctionReturn(PETSC_SUCCESS);
3522: }
3524: PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3525: {
3526: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3528: PetscFunctionBegin;
3529: PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first");
3530: PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first");
3531: if (a->cmat && (cend - cbegin != a->cmat->cmap->N || rend - rbegin != a->cmat->rmap->N)) PetscCall(MatDestroy(&a->cmat));
3532: if (!a->cmat) {
3533: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), rend - rbegin, PETSC_DECIDE, rend - rbegin, cend - cbegin, a->v ? a->v + rbegin + (size_t)cbegin * a->lda : NULL, &a->cmat));
3534: } else {
3535: PetscCall(MatDensePlaceArray(a->cmat, a->v ? a->v + rbegin + (size_t)cbegin * a->lda : NULL));
3536: }
3537: PetscCall(MatDenseSetLDA(a->cmat, a->lda));
3538: a->matinuse = cbegin + 1;
3539: *v = a->cmat;
3540: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
3541: A->offloadmask = PETSC_OFFLOAD_CPU;
3542: #endif
3543: PetscFunctionReturn(PETSC_SUCCESS);
3544: }
3546: PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A, Mat *v)
3547: {
3548: Mat_SeqDense *a = (Mat_SeqDense *)A->data;
3550: PetscFunctionBegin;
3551: PetscCheck(a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetSubMatrix() first");
3552: PetscCheck(a->cmat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column matrix");
3553: PetscCheck(*v == a->cmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not the matrix obtained from MatDenseGetSubMatrix()");
3554: a->matinuse = 0;
3555: PetscCall(MatDenseResetArray(a->cmat));
3556: if (v) *v = NULL;
3557: #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
3558: A->offloadmask = PETSC_OFFLOAD_CPU;
3559: #endif
3560: PetscFunctionReturn(PETSC_SUCCESS);
3561: }
3563: /*MC
3564: MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices.
3566: Options Database Key:
3567: . -mat_type seqdense - sets the matrix type to `MATSEQDENSE` during a call to `MatSetFromOptions()`
3569: Level: beginner
3571: .seealso: [](ch_matrices), `Mat`, `MATSEQDENSE`, `MatCreateSeqDense()`
3572: M*/
3573: PetscErrorCode MatCreate_SeqDense(Mat B)
3574: {
3575: Mat_SeqDense *b;
3576: PetscMPIInt size;
3578: PetscFunctionBegin;
3579: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
3580: PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1");
3582: PetscCall(PetscNew(&b));
3583: B->data = (void *)b;
3584: B->ops[0] = MatOps_Values;
3586: b->roworiented = PETSC_TRUE;
3588: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactor_C", MatQRFactor_SeqDense));
3589: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetLDA_C", MatDenseGetLDA_SeqDense));
3590: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseSetLDA_C", MatDenseSetLDA_SeqDense));
3591: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArray_C", MatDenseGetArray_SeqDense));
3592: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArray_C", MatDenseRestoreArray_SeqDense));
3593: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDensePlaceArray_C", MatDensePlaceArray_SeqDense));
3594: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseResetArray_C", MatDenseResetArray_SeqDense));
3595: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseReplaceArray_C", MatDenseReplaceArray_SeqDense));
3596: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayRead_C", MatDenseGetArray_SeqDense));
3597: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayRead_C", MatDenseRestoreArray_SeqDense));
3598: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayWrite_C", MatDenseGetArray_SeqDense));
3599: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayWrite_C", MatDenseRestoreArray_SeqDense));
3600: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqaij_C", MatConvert_SeqDense_SeqAIJ));
3601: #if defined(PETSC_HAVE_ELEMENTAL)
3602: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_elemental_C", MatConvert_SeqDense_Elemental));
3603: #endif
3604: #if defined(PETSC_HAVE_SCALAPACK)
3605: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_scalapack_C", MatConvert_Dense_ScaLAPACK));
3606: #endif
3607: #if defined(PETSC_HAVE_CUDA)
3608: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensecuda_C", MatConvert_SeqDense_SeqDenseCUDA));
3609: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3610: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdense_C", MatProductSetFromOptions_SeqDense));
3611: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensecuda_C", MatProductSetFromOptions_SeqDense));
3612: #endif
3613: #if defined(PETSC_HAVE_HIP)
3614: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensehip_C", MatConvert_SeqDense_SeqDenseHIP));
3615: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3616: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdense_C", MatProductSetFromOptions_SeqDense));
3617: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensehip_C", MatProductSetFromOptions_SeqDense));
3618: #endif
3619: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqDenseSetPreallocation_C", MatSeqDenseSetPreallocation_SeqDense));
3620: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqdense_C", MatProductSetFromOptions_SeqAIJ_SeqDense));
3621: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdense_C", MatProductSetFromOptions_SeqDense));
3622: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));
3623: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqsbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense));
3625: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumn_C", MatDenseGetColumn_SeqDense));
3626: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumn_C", MatDenseRestoreColumn_SeqDense));
3627: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVec_C", MatDenseGetColumnVec_SeqDense));
3628: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVec_C", MatDenseRestoreColumnVec_SeqDense));
3629: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecRead_C", MatDenseGetColumnVecRead_SeqDense));
3630: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecRead_C", MatDenseRestoreColumnVecRead_SeqDense));
3631: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecWrite_C", MatDenseGetColumnVecWrite_SeqDense));
3632: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecWrite_C", MatDenseRestoreColumnVecWrite_SeqDense));
3633: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetSubMatrix_C", MatDenseGetSubMatrix_SeqDense));
3634: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreSubMatrix_C", MatDenseRestoreSubMatrix_SeqDense));
3635: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQDENSE));
3636: PetscFunctionReturn(PETSC_SUCCESS);
3637: }
3639: /*@C
3640: 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.
3642: Not Collective
3644: Input Parameters:
3645: + A - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3646: - col - column index
3648: Output Parameter:
3649: . vals - pointer to the data
3651: Level: intermediate
3653: Note:
3654: Use `MatDenseGetColumnVec()` to get access to a column of a `MATDENSE` treated as a `Vec`
3656: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreColumn()`, `MatDenseGetColumnVec()`
3657: @*/
3658: PetscErrorCode MatDenseGetColumn(Mat A, PetscInt col, PetscScalar **vals)
3659: {
3660: PetscFunctionBegin;
3663: PetscAssertPointer(vals, 3);
3664: PetscUseMethod(A, "MatDenseGetColumn_C", (Mat, PetscInt, PetscScalar **), (A, col, vals));
3665: PetscFunctionReturn(PETSC_SUCCESS);
3666: }
3668: /*@C
3669: MatDenseRestoreColumn - returns access to a column of a `MATDENSE` matrix which is returned by `MatDenseGetColumn()`.
3671: Not Collective
3673: Input Parameters:
3674: + A - a `MATSEQDENSE` or `MATMPIDENSE` matrix
3675: - vals - pointer to the data (may be `NULL`)
3677: Level: intermediate
3679: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MatDenseGetColumn()`
3680: @*/
3681: PetscErrorCode MatDenseRestoreColumn(Mat A, PetscScalar **vals)
3682: {
3683: PetscFunctionBegin;
3685: PetscAssertPointer(vals, 2);
3686: PetscUseMethod(A, "MatDenseRestoreColumn_C", (Mat, PetscScalar **), (A, vals));
3687: PetscFunctionReturn(PETSC_SUCCESS);
3688: }
3690: /*@
3691: MatDenseGetColumnVec - Gives read-write access to a column of a `MATDENSE` matrix, represented as a `Vec`.
3693: Collective
3695: Input Parameters:
3696: + A - the `Mat` object
3697: - col - the column index
3699: Output Parameter:
3700: . v - the vector
3702: Level: intermediate
3704: Notes:
3705: The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVec()` when the vector is no longer needed.
3707: Use `MatDenseGetColumnVecRead()` to obtain read-only access or `MatDenseGetColumnVecWrite()` for write-only access.
3709: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`, `MatDenseGetColumn()`
3710: @*/
3711: PetscErrorCode MatDenseGetColumnVec(Mat A, PetscInt col, Vec *v)
3712: {
3713: PetscFunctionBegin;
3717: PetscAssertPointer(v, 3);
3718: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3719: 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);
3720: PetscUseMethod(A, "MatDenseGetColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3721: PetscFunctionReturn(PETSC_SUCCESS);
3722: }
3724: /*@
3725: MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVec()`.
3727: Collective
3729: Input Parameters:
3730: + A - the `Mat` object
3731: . col - the column index
3732: - v - the `Vec` object (may be `NULL`)
3734: Level: intermediate
3736: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3737: @*/
3738: PetscErrorCode MatDenseRestoreColumnVec(Mat A, PetscInt col, Vec *v)
3739: {
3740: PetscFunctionBegin;
3744: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3745: 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);
3746: PetscUseMethod(A, "MatDenseRestoreColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v));
3747: PetscFunctionReturn(PETSC_SUCCESS);
3748: }
3750: /*@
3751: MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a `Vec`.
3753: Collective
3755: Input Parameters:
3756: + A - the `Mat` object
3757: - col - the column index
3759: Output Parameter:
3760: . v - the vector
3762: Level: intermediate
3764: Notes:
3765: The vector is owned by PETSc and users cannot modify it.
3767: Users need to call `MatDenseRestoreColumnVecRead()` when the vector is no longer needed.
3769: Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecWrite()` for write-only access.
3771: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3772: @*/
3773: PetscErrorCode MatDenseGetColumnVecRead(Mat A, PetscInt col, Vec *v)
3774: {
3775: PetscFunctionBegin;
3779: PetscAssertPointer(v, 3);
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, "MatDenseGetColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3783: PetscFunctionReturn(PETSC_SUCCESS);
3784: }
3786: /*@
3787: MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVecRead()`.
3789: Collective
3791: Input Parameters:
3792: + A - the `Mat` object
3793: . col - the column index
3794: - v - the `Vec` object (may be `NULL`)
3796: Level: intermediate
3798: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecWrite()`
3799: @*/
3800: PetscErrorCode MatDenseRestoreColumnVecRead(Mat A, PetscInt col, Vec *v)
3801: {
3802: PetscFunctionBegin;
3806: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3807: 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);
3808: PetscUseMethod(A, "MatDenseRestoreColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v));
3809: PetscFunctionReturn(PETSC_SUCCESS);
3810: }
3812: /*@
3813: MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a `Vec`.
3815: Collective
3817: Input Parameters:
3818: + A - the `Mat` object
3819: - col - the column index
3821: Output Parameter:
3822: . v - the vector
3824: Level: intermediate
3826: Notes:
3827: The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVecWrite()` when the vector is no longer needed.
3829: Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecRead()` for read-only access.
3831: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`
3832: @*/
3833: PetscErrorCode MatDenseGetColumnVecWrite(Mat A, PetscInt col, Vec *v)
3834: {
3835: PetscFunctionBegin;
3839: PetscAssertPointer(v, 3);
3840: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3841: 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);
3842: PetscUseMethod(A, "MatDenseGetColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3843: PetscFunctionReturn(PETSC_SUCCESS);
3844: }
3846: /*@
3847: MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from `MatDenseGetColumnVecWrite()`.
3849: Collective
3851: Input Parameters:
3852: + A - the `Mat` object
3853: . col - the column index
3854: - v - the `Vec` object (may be `NULL`)
3856: Level: intermediate
3858: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`
3859: @*/
3860: PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A, PetscInt col, Vec *v)
3861: {
3862: PetscFunctionBegin;
3866: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3867: 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);
3868: PetscUseMethod(A, "MatDenseRestoreColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v));
3869: PetscFunctionReturn(PETSC_SUCCESS);
3870: }
3872: /*@
3873: MatDenseGetSubMatrix - Gives access to a block of rows and columns of a dense matrix, represented as a `Mat`.
3875: Collective
3877: Input Parameters:
3878: + A - the `Mat` object
3879: . rbegin - the first global row index in the block (if `PETSC_DECIDE`, is 0)
3880: . rend - the global row index past the last one in the block (if `PETSC_DECIDE`, is `M`)
3881: . cbegin - the first global column index in the block (if `PETSC_DECIDE`, is 0)
3882: - cend - the global column index past the last one in the block (if `PETSC_DECIDE`, is `N`)
3884: Output Parameter:
3885: . v - the matrix
3887: Level: intermediate
3889: Notes:
3890: The matrix is owned by PETSc. Users need to call `MatDenseRestoreSubMatrix()` when the matrix is no longer needed.
3892: The output matrix is not redistributed by PETSc, so depending on the values of `rbegin` and `rend`, some processes may have no local rows.
3894: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreSubMatrix()`
3895: @*/
3896: PetscErrorCode MatDenseGetSubMatrix(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v)
3897: {
3898: PetscFunctionBegin;
3905: PetscAssertPointer(v, 6);
3906: if (rbegin == PETSC_DECIDE) rbegin = 0;
3907: if (rend == PETSC_DECIDE) rend = A->rmap->N;
3908: if (cbegin == PETSC_DECIDE) cbegin = 0;
3909: if (cend == PETSC_DECIDE) cend = A->cmap->N;
3910: PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated");
3911: 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);
3912: 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);
3913: 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);
3914: 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);
3915: PetscUseMethod(A, "MatDenseGetSubMatrix_C", (Mat, PetscInt, PetscInt, PetscInt, PetscInt, Mat *), (A, rbegin, rend, cbegin, cend, v));
3916: PetscFunctionReturn(PETSC_SUCCESS);
3917: }
3919: /*@
3920: MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from `MatDenseGetSubMatrix()`.
3922: Collective
3924: Input Parameters:
3925: + A - the `Mat` object
3926: - v - the `Mat` object (may be `NULL`)
3928: Level: intermediate
3930: .seealso: [](ch_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseGetSubMatrix()`
3931: @*/
3932: PetscErrorCode MatDenseRestoreSubMatrix(Mat A, Mat *v)
3933: {
3934: PetscFunctionBegin;
3937: PetscAssertPointer(v, 2);
3938: PetscUseMethod(A, "MatDenseRestoreSubMatrix_C", (Mat, Mat *), (A, v));
3939: PetscFunctionReturn(PETSC_SUCCESS);
3940: }
3942: #include <petscblaslapack.h>
3943: #include <petsc/private/kernels/blockinvert.h>
3945: PetscErrorCode MatSeqDenseInvert(Mat A)
3946: {
3947: PetscInt m;
3948: const PetscReal shift = 0.0;
3949: PetscBool allowzeropivot, zeropivotdetected = PETSC_FALSE;
3950: PetscScalar *values;
3952: PetscFunctionBegin;
3954: PetscCall(MatDenseGetArray(A, &values));
3955: PetscCall(MatGetLocalSize(A, &m, NULL));
3956: allowzeropivot = PetscNot(A->erroriffailure);
3957: /* factor and invert each block */
3958: switch (m) {
3959: case 1:
3960: values[0] = (PetscScalar)1.0 / (values[0] + shift);
3961: break;
3962: case 2:
3963: PetscCall(PetscKernel_A_gets_inverse_A_2(values, shift, allowzeropivot, &zeropivotdetected));
3964: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3965: break;
3966: case 3:
3967: PetscCall(PetscKernel_A_gets_inverse_A_3(values, shift, allowzeropivot, &zeropivotdetected));
3968: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3969: break;
3970: case 4:
3971: PetscCall(PetscKernel_A_gets_inverse_A_4(values, shift, allowzeropivot, &zeropivotdetected));
3972: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3973: break;
3974: case 5: {
3975: PetscScalar work[25];
3976: PetscInt ipvt[5];
3978: PetscCall(PetscKernel_A_gets_inverse_A_5(values, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
3979: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3980: } break;
3981: case 6:
3982: PetscCall(PetscKernel_A_gets_inverse_A_6(values, shift, allowzeropivot, &zeropivotdetected));
3983: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3984: break;
3985: case 7:
3986: PetscCall(PetscKernel_A_gets_inverse_A_7(values, shift, allowzeropivot, &zeropivotdetected));
3987: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3988: break;
3989: default: {
3990: PetscInt *v_pivots, *IJ, j;
3991: PetscScalar *v_work;
3993: PetscCall(PetscMalloc3(m, &v_work, m, &v_pivots, m, &IJ));
3994: for (j = 0; j < m; j++) IJ[j] = j;
3995: PetscCall(PetscKernel_A_gets_inverse_A(m, values, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
3996: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
3997: PetscCall(PetscFree3(v_work, v_pivots, IJ));
3998: }
3999: }
4000: PetscCall(MatDenseRestoreArray(A, &values));
4001: PetscFunctionReturn(PETSC_SUCCESS);
4002: }