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