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., unless it is currently bound to CPU, 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(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI)); 2376 if (isMPI) { 2377 /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */ 2378 PetscCall(MatDenseGetArrayAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype)); 2379 } else { 2380 PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *); 2381 PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayAndMemType_C", &fptr); 2382 if (fptr) { 2383 PetscCall((*fptr)(A, array, mtype)); 2384 } else { 2385 PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array)); 2386 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 2387 } 2388 } 2389 PetscFunctionReturn(0); 2390 } 2391 2392 /*@C 2393 MatDenseRestoreArrayAndMemType - returns access to the array that is obtained by `MatDenseGetArrayAndMemType()` 2394 2395 Logically Collective 2396 2397 Input Parameters: 2398 + mat - a dense matrix 2399 - array - pointer to the data 2400 2401 Level: intermediate 2402 2403 .seealso: `MATDENSE`, `MatDenseGetArrayAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2404 @*/ 2405 PetscErrorCode MatDenseRestoreArrayAndMemType(Mat A, PetscScalar **array) 2406 { 2407 PetscBool isMPI; 2408 2409 PetscFunctionBegin; 2410 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2411 PetscValidPointer(array, 2); 2412 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI)); 2413 if (isMPI) { 2414 PetscCall(MatDenseRestoreArrayAndMemType(((Mat_MPIDense *)A->data)->A, array)); 2415 } else { 2416 PetscErrorCode (*fptr)(Mat, PetscScalar **); 2417 PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayAndMemType_C", &fptr); 2418 if (fptr) { 2419 PetscCall((*fptr)(A, array)); 2420 } else { 2421 PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array)); 2422 } 2423 *array = NULL; 2424 } 2425 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 2426 PetscFunctionReturn(0); 2427 } 2428 2429 /*@C 2430 MatDenseGetArrayReadAndMemType - gives read-only access to the array where the data for a `MATDENSE` matrix is stored 2431 2432 Logically Collective 2433 2434 Input Parameter: 2435 . mat - a dense matrix 2436 2437 Output Parameters: 2438 + array - pointer to the data 2439 - mtype - memory type of the returned pointer 2440 2441 Notes: 2442 If the matrix is of a device type such as MATDENSECUDA, MATDENSEHIP, etc., unless it is currently bound to CPU, 2443 an array on device is always returned and is guaranteed to contain the matrix's latest data. 2444 2445 Level: intermediate 2446 2447 .seealso: `MATDENSE`, `MatDenseRestoreArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`, 2448 `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()` 2449 @*/ 2450 PetscErrorCode MatDenseGetArrayReadAndMemType(Mat A, const PetscScalar **array, PetscMemType *mtype) 2451 { 2452 PetscBool isMPI; 2453 2454 PetscFunctionBegin; 2455 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2456 PetscValidPointer(array, 2); 2457 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI)); 2458 if (isMPI) { /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */ 2459 PetscCall(MatDenseGetArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype)); 2460 } else { 2461 PetscErrorCode (*fptr)(Mat, const PetscScalar **, PetscMemType *); 2462 PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayReadAndMemType_C", &fptr); 2463 if (fptr) { 2464 PetscCall((*fptr)(A, array, mtype)); 2465 } else { 2466 PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, const PetscScalar **), (A, array)); 2467 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 2468 } 2469 } 2470 PetscFunctionReturn(0); 2471 } 2472 2473 /*@C 2474 MatDenseRestoreArrayReadAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()` 2475 2476 Logically Collective 2477 2478 Input Parameters: 2479 + mat - a dense matrix 2480 - array - pointer to the data 2481 2482 Level: intermediate 2483 2484 .seealso: `MATDENSE`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2485 @*/ 2486 PetscErrorCode MatDenseRestoreArrayReadAndMemType(Mat A, const PetscScalar **array) 2487 { 2488 PetscBool isMPI; 2489 2490 PetscFunctionBegin; 2491 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2492 PetscValidPointer(array, 2); 2493 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI)); 2494 if (isMPI) { 2495 PetscCall(MatDenseRestoreArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array)); 2496 } else { 2497 PetscErrorCode (*fptr)(Mat, const PetscScalar **); 2498 PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayReadAndMemType_C", &fptr); 2499 if (fptr) { 2500 PetscCall((*fptr)(A, array)); 2501 } else { 2502 PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, const PetscScalar **), (A, array)); 2503 } 2504 *array = NULL; 2505 } 2506 PetscFunctionReturn(0); 2507 } 2508 2509 /*@C 2510 MatDenseGetArrayWriteAndMemType - gives write-only access to the array where the data for a `MATDENSE` matrix is stored 2511 2512 Logically Collective 2513 2514 Input Parameter: 2515 . mat - a dense matrix 2516 2517 Output Parameters: 2518 + array - pointer to the data 2519 - mtype - memory type of the returned pointer 2520 2521 Notes: 2522 If the matrix is of a device type such as MATDENSECUDA, MATDENSEHIP, etc., unless it is currently bound to CPU, 2523 an array on device is always returned and is guaranteed to contain the matrix's latest data. 2524 2525 Level: intermediate 2526 2527 .seealso: `MATDENSE`, `MatDenseRestoreArrayWriteAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayRead()`, 2528 `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()` 2529 @*/ 2530 PetscErrorCode MatDenseGetArrayWriteAndMemType(Mat A, PetscScalar **array, PetscMemType *mtype) 2531 { 2532 PetscBool isMPI; 2533 2534 PetscFunctionBegin; 2535 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2536 PetscValidPointer(array, 2); 2537 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI)); 2538 if (isMPI) { 2539 PetscCall(MatDenseGetArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype)); 2540 } else { 2541 PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *); 2542 PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayWriteAndMemType_C", &fptr); 2543 if (fptr) { 2544 PetscCall((*fptr)(A, array, mtype)); 2545 } else { 2546 PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array)); 2547 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 2548 } 2549 } 2550 PetscFunctionReturn(0); 2551 } 2552 2553 /*@C 2554 MatDenseRestoreArrayWriteAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()` 2555 2556 Logically Collective 2557 2558 Input Parameters: 2559 + mat - a dense matrix 2560 - array - pointer to the data 2561 2562 Level: intermediate 2563 2564 .seealso: `MATDENSE`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2565 @*/ 2566 PetscErrorCode MatDenseRestoreArrayWriteAndMemType(Mat A, PetscScalar **array) 2567 { 2568 PetscBool isMPI; 2569 2570 PetscFunctionBegin; 2571 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2572 PetscValidPointer(array, 2); 2573 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI)); 2574 if (isMPI) { 2575 PetscCall(MatDenseRestoreArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array)); 2576 } else { 2577 PetscErrorCode (*fptr)(Mat, PetscScalar **); 2578 PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayWriteAndMemType_C", &fptr); 2579 if (fptr) { 2580 PetscCall((*fptr)(A, array)); 2581 } else { 2582 PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array)); 2583 } 2584 *array = NULL; 2585 } 2586 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 2587 PetscFunctionReturn(0); 2588 } 2589 2590 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B) 2591 { 2592 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 2593 PetscInt i, j, nrows, ncols, ldb; 2594 const PetscInt *irow, *icol; 2595 PetscScalar *av, *bv, *v = mat->v; 2596 Mat newmat; 2597 2598 PetscFunctionBegin; 2599 PetscCall(ISGetIndices(isrow, &irow)); 2600 PetscCall(ISGetIndices(iscol, &icol)); 2601 PetscCall(ISGetLocalSize(isrow, &nrows)); 2602 PetscCall(ISGetLocalSize(iscol, &ncols)); 2603 2604 /* Check submatrixcall */ 2605 if (scall == MAT_REUSE_MATRIX) { 2606 PetscInt n_cols, n_rows; 2607 PetscCall(MatGetSize(*B, &n_rows, &n_cols)); 2608 if (n_rows != nrows || n_cols != ncols) { 2609 /* resize the result matrix to match number of requested rows/columns */ 2610 PetscCall(MatSetSizes(*B, nrows, ncols, nrows, ncols)); 2611 } 2612 newmat = *B; 2613 } else { 2614 /* Create and fill new matrix */ 2615 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &newmat)); 2616 PetscCall(MatSetSizes(newmat, nrows, ncols, nrows, ncols)); 2617 PetscCall(MatSetType(newmat, ((PetscObject)A)->type_name)); 2618 PetscCall(MatSeqDenseSetPreallocation(newmat, NULL)); 2619 } 2620 2621 /* Now extract the data pointers and do the copy,column at a time */ 2622 PetscCall(MatDenseGetArray(newmat, &bv)); 2623 PetscCall(MatDenseGetLDA(newmat, &ldb)); 2624 for (i = 0; i < ncols; i++) { 2625 av = v + mat->lda * icol[i]; 2626 for (j = 0; j < nrows; j++) bv[j] = av[irow[j]]; 2627 bv += ldb; 2628 } 2629 PetscCall(MatDenseRestoreArray(newmat, &bv)); 2630 2631 /* Assemble the matrices so that the correct flags are set */ 2632 PetscCall(MatAssemblyBegin(newmat, MAT_FINAL_ASSEMBLY)); 2633 PetscCall(MatAssemblyEnd(newmat, MAT_FINAL_ASSEMBLY)); 2634 2635 /* Free work space */ 2636 PetscCall(ISRestoreIndices(isrow, &irow)); 2637 PetscCall(ISRestoreIndices(iscol, &icol)); 2638 *B = newmat; 2639 PetscFunctionReturn(0); 2640 } 2641 2642 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[]) 2643 { 2644 PetscInt i; 2645 2646 PetscFunctionBegin; 2647 if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n, B)); 2648 2649 for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SeqDense(A, irow[i], icol[i], scall, &(*B)[i])); 2650 PetscFunctionReturn(0); 2651 } 2652 2653 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat, MatAssemblyType mode) 2654 { 2655 PetscFunctionBegin; 2656 PetscFunctionReturn(0); 2657 } 2658 2659 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat, MatAssemblyType mode) 2660 { 2661 PetscFunctionBegin; 2662 PetscFunctionReturn(0); 2663 } 2664 2665 PetscErrorCode MatCopy_SeqDense(Mat A, Mat B, MatStructure str) 2666 { 2667 Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data; 2668 const PetscScalar *va; 2669 PetscScalar *vb; 2670 PetscInt lda1 = a->lda, lda2 = b->lda, m = A->rmap->n, n = A->cmap->n, j; 2671 2672 PetscFunctionBegin; 2673 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2674 if (A->ops->copy != B->ops->copy) { 2675 PetscCall(MatCopy_Basic(A, B, str)); 2676 PetscFunctionReturn(0); 2677 } 2678 PetscCheck(m == B->rmap->n && n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "size(B) != size(A)"); 2679 PetscCall(MatDenseGetArrayRead(A, &va)); 2680 PetscCall(MatDenseGetArray(B, &vb)); 2681 if (lda1 > m || lda2 > m) { 2682 for (j = 0; j < n; j++) PetscCall(PetscArraycpy(vb + j * lda2, va + j * lda1, m)); 2683 } else { 2684 PetscCall(PetscArraycpy(vb, va, A->rmap->n * A->cmap->n)); 2685 } 2686 PetscCall(MatDenseRestoreArray(B, &vb)); 2687 PetscCall(MatDenseRestoreArrayRead(A, &va)); 2688 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 2689 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 2690 PetscFunctionReturn(0); 2691 } 2692 2693 PetscErrorCode MatSetUp_SeqDense(Mat A) 2694 { 2695 PetscFunctionBegin; 2696 PetscCall(PetscLayoutSetUp(A->rmap)); 2697 PetscCall(PetscLayoutSetUp(A->cmap)); 2698 if (!A->preallocated) PetscCall(MatSeqDenseSetPreallocation(A, NULL)); 2699 PetscFunctionReturn(0); 2700 } 2701 2702 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2703 { 2704 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 2705 PetscInt i, j; 2706 PetscInt min = PetscMin(A->rmap->n, A->cmap->n); 2707 PetscScalar *aa; 2708 2709 PetscFunctionBegin; 2710 PetscCall(MatDenseGetArray(A, &aa)); 2711 for (j = 0; j < A->cmap->n; j++) { 2712 for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscConj(aa[i + j * mat->lda]); 2713 } 2714 PetscCall(MatDenseRestoreArray(A, &aa)); 2715 if (mat->tau) 2716 for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]); 2717 PetscFunctionReturn(0); 2718 } 2719 2720 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2721 { 2722 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 2723 PetscInt i, j; 2724 PetscScalar *aa; 2725 2726 PetscFunctionBegin; 2727 PetscCall(MatDenseGetArray(A, &aa)); 2728 for (j = 0; j < A->cmap->n; j++) { 2729 for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscRealPart(aa[i + j * mat->lda]); 2730 } 2731 PetscCall(MatDenseRestoreArray(A, &aa)); 2732 PetscFunctionReturn(0); 2733 } 2734 2735 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2736 { 2737 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 2738 PetscInt i, j; 2739 PetscScalar *aa; 2740 2741 PetscFunctionBegin; 2742 PetscCall(MatDenseGetArray(A, &aa)); 2743 for (j = 0; j < A->cmap->n; j++) { 2744 for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscImaginaryPart(aa[i + j * mat->lda]); 2745 } 2746 PetscCall(MatDenseRestoreArray(A, &aa)); 2747 PetscFunctionReturn(0); 2748 } 2749 2750 /* ----------------------------------------------------------------*/ 2751 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C) 2752 { 2753 PetscInt m = A->rmap->n, n = B->cmap->n; 2754 PetscBool cisdense = PETSC_FALSE; 2755 2756 PetscFunctionBegin; 2757 PetscCall(MatSetSizes(C, m, n, m, n)); 2758 #if defined(PETSC_HAVE_CUDA) 2759 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, "")); 2760 #endif 2761 #if defined(PETSC_HAVE_HIP) 2762 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, "")); 2763 #endif 2764 if (!cisdense) { 2765 PetscBool flg; 2766 2767 PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg)); 2768 PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE)); 2769 } 2770 PetscCall(MatSetUp(C)); 2771 PetscFunctionReturn(0); 2772 } 2773 2774 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C) 2775 { 2776 Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data, *c = (Mat_SeqDense *)C->data; 2777 PetscBLASInt m, n, k; 2778 const PetscScalar *av, *bv; 2779 PetscScalar *cv; 2780 PetscScalar _DOne = 1.0, _DZero = 0.0; 2781 2782 PetscFunctionBegin; 2783 PetscCall(PetscBLASIntCast(C->rmap->n, &m)); 2784 PetscCall(PetscBLASIntCast(C->cmap->n, &n)); 2785 PetscCall(PetscBLASIntCast(A->cmap->n, &k)); 2786 if (!m || !n || !k) PetscFunctionReturn(0); 2787 PetscCall(MatDenseGetArrayRead(A, &av)); 2788 PetscCall(MatDenseGetArrayRead(B, &bv)); 2789 PetscCall(MatDenseGetArrayWrite(C, &cv)); 2790 PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda)); 2791 PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1))); 2792 PetscCall(MatDenseRestoreArrayRead(A, &av)); 2793 PetscCall(MatDenseRestoreArrayRead(B, &bv)); 2794 PetscCall(MatDenseRestoreArrayWrite(C, &cv)); 2795 PetscFunctionReturn(0); 2796 } 2797 2798 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C) 2799 { 2800 PetscInt m = A->rmap->n, n = B->rmap->n; 2801 PetscBool cisdense = PETSC_FALSE; 2802 2803 PetscFunctionBegin; 2804 PetscCall(MatSetSizes(C, m, n, m, n)); 2805 #if defined(PETSC_HAVE_CUDA) 2806 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, "")); 2807 #endif 2808 #if defined(PETSC_HAVE_HIP) 2809 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, "")); 2810 #endif 2811 if (!cisdense) { 2812 PetscBool flg; 2813 2814 PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg)); 2815 PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE)); 2816 } 2817 PetscCall(MatSetUp(C)); 2818 PetscFunctionReturn(0); 2819 } 2820 2821 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C) 2822 { 2823 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 2824 Mat_SeqDense *b = (Mat_SeqDense *)B->data; 2825 Mat_SeqDense *c = (Mat_SeqDense *)C->data; 2826 const PetscScalar *av, *bv; 2827 PetscScalar *cv; 2828 PetscBLASInt m, n, k; 2829 PetscScalar _DOne = 1.0, _DZero = 0.0; 2830 2831 PetscFunctionBegin; 2832 PetscCall(PetscBLASIntCast(C->rmap->n, &m)); 2833 PetscCall(PetscBLASIntCast(C->cmap->n, &n)); 2834 PetscCall(PetscBLASIntCast(A->cmap->n, &k)); 2835 if (!m || !n || !k) PetscFunctionReturn(0); 2836 PetscCall(MatDenseGetArrayRead(A, &av)); 2837 PetscCall(MatDenseGetArrayRead(B, &bv)); 2838 PetscCall(MatDenseGetArrayWrite(C, &cv)); 2839 PetscCallBLAS("BLASgemm", BLASgemm_("N", "T", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda)); 2840 PetscCall(MatDenseRestoreArrayRead(A, &av)); 2841 PetscCall(MatDenseRestoreArrayRead(B, &bv)); 2842 PetscCall(MatDenseRestoreArrayWrite(C, &cv)); 2843 PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1))); 2844 PetscFunctionReturn(0); 2845 } 2846 2847 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C) 2848 { 2849 PetscInt m = A->cmap->n, n = B->cmap->n; 2850 PetscBool cisdense = PETSC_FALSE; 2851 2852 PetscFunctionBegin; 2853 PetscCall(MatSetSizes(C, m, n, m, n)); 2854 #if defined(PETSC_HAVE_CUDA) 2855 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, "")); 2856 #endif 2857 #if defined(PETSC_HAVE_HIP) 2858 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, "")); 2859 #endif 2860 if (!cisdense) { 2861 PetscBool flg; 2862 2863 PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg)); 2864 PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE)); 2865 } 2866 PetscCall(MatSetUp(C)); 2867 PetscFunctionReturn(0); 2868 } 2869 2870 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C) 2871 { 2872 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 2873 Mat_SeqDense *b = (Mat_SeqDense *)B->data; 2874 Mat_SeqDense *c = (Mat_SeqDense *)C->data; 2875 const PetscScalar *av, *bv; 2876 PetscScalar *cv; 2877 PetscBLASInt m, n, k; 2878 PetscScalar _DOne = 1.0, _DZero = 0.0; 2879 2880 PetscFunctionBegin; 2881 PetscCall(PetscBLASIntCast(C->rmap->n, &m)); 2882 PetscCall(PetscBLASIntCast(C->cmap->n, &n)); 2883 PetscCall(PetscBLASIntCast(A->rmap->n, &k)); 2884 if (!m || !n || !k) PetscFunctionReturn(0); 2885 PetscCall(MatDenseGetArrayRead(A, &av)); 2886 PetscCall(MatDenseGetArrayRead(B, &bv)); 2887 PetscCall(MatDenseGetArrayWrite(C, &cv)); 2888 PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda)); 2889 PetscCall(MatDenseRestoreArrayRead(A, &av)); 2890 PetscCall(MatDenseRestoreArrayRead(B, &bv)); 2891 PetscCall(MatDenseRestoreArrayWrite(C, &cv)); 2892 PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1))); 2893 PetscFunctionReturn(0); 2894 } 2895 2896 /* ----------------------------------------------- */ 2897 static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C) 2898 { 2899 PetscFunctionBegin; 2900 C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense; 2901 C->ops->productsymbolic = MatProductSymbolic_AB; 2902 PetscFunctionReturn(0); 2903 } 2904 2905 static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C) 2906 { 2907 PetscFunctionBegin; 2908 C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense; 2909 C->ops->productsymbolic = MatProductSymbolic_AtB; 2910 PetscFunctionReturn(0); 2911 } 2912 2913 static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C) 2914 { 2915 PetscFunctionBegin; 2916 C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense; 2917 C->ops->productsymbolic = MatProductSymbolic_ABt; 2918 PetscFunctionReturn(0); 2919 } 2920 2921 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C) 2922 { 2923 Mat_Product *product = C->product; 2924 2925 PetscFunctionBegin; 2926 switch (product->type) { 2927 case MATPRODUCT_AB: 2928 PetscCall(MatProductSetFromOptions_SeqDense_AB(C)); 2929 break; 2930 case MATPRODUCT_AtB: 2931 PetscCall(MatProductSetFromOptions_SeqDense_AtB(C)); 2932 break; 2933 case MATPRODUCT_ABt: 2934 PetscCall(MatProductSetFromOptions_SeqDense_ABt(C)); 2935 break; 2936 default: 2937 break; 2938 } 2939 PetscFunctionReturn(0); 2940 } 2941 /* ----------------------------------------------- */ 2942 2943 static PetscErrorCode MatGetRowMax_SeqDense(Mat A, Vec v, PetscInt idx[]) 2944 { 2945 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 2946 PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p; 2947 PetscScalar *x; 2948 const PetscScalar *aa; 2949 2950 PetscFunctionBegin; 2951 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 2952 PetscCall(VecGetArray(v, &x)); 2953 PetscCall(VecGetLocalSize(v, &p)); 2954 PetscCall(MatDenseGetArrayRead(A, &aa)); 2955 PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector"); 2956 for (i = 0; i < m; i++) { 2957 x[i] = aa[i]; 2958 if (idx) idx[i] = 0; 2959 for (j = 1; j < n; j++) { 2960 if (PetscRealPart(x[i]) < PetscRealPart(aa[i + a->lda * j])) { 2961 x[i] = aa[i + a->lda * j]; 2962 if (idx) idx[i] = j; 2963 } 2964 } 2965 } 2966 PetscCall(MatDenseRestoreArrayRead(A, &aa)); 2967 PetscCall(VecRestoreArray(v, &x)); 2968 PetscFunctionReturn(0); 2969 } 2970 2971 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A, Vec v, PetscInt idx[]) 2972 { 2973 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 2974 PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p; 2975 PetscScalar *x; 2976 PetscReal atmp; 2977 const PetscScalar *aa; 2978 2979 PetscFunctionBegin; 2980 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 2981 PetscCall(VecGetArray(v, &x)); 2982 PetscCall(VecGetLocalSize(v, &p)); 2983 PetscCall(MatDenseGetArrayRead(A, &aa)); 2984 PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector"); 2985 for (i = 0; i < m; i++) { 2986 x[i] = PetscAbsScalar(aa[i]); 2987 for (j = 1; j < n; j++) { 2988 atmp = PetscAbsScalar(aa[i + a->lda * j]); 2989 if (PetscAbsScalar(x[i]) < atmp) { 2990 x[i] = atmp; 2991 if (idx) idx[i] = j; 2992 } 2993 } 2994 } 2995 PetscCall(MatDenseRestoreArrayRead(A, &aa)); 2996 PetscCall(VecRestoreArray(v, &x)); 2997 PetscFunctionReturn(0); 2998 } 2999 3000 static PetscErrorCode MatGetRowMin_SeqDense(Mat A, Vec v, PetscInt idx[]) 3001 { 3002 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3003 PetscInt i, j, m = A->rmap->n, n = A->cmap->n, p; 3004 PetscScalar *x; 3005 const PetscScalar *aa; 3006 3007 PetscFunctionBegin; 3008 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3009 PetscCall(MatDenseGetArrayRead(A, &aa)); 3010 PetscCall(VecGetArray(v, &x)); 3011 PetscCall(VecGetLocalSize(v, &p)); 3012 PetscCheck(p == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector"); 3013 for (i = 0; i < m; i++) { 3014 x[i] = aa[i]; 3015 if (idx) idx[i] = 0; 3016 for (j = 1; j < n; j++) { 3017 if (PetscRealPart(x[i]) > PetscRealPart(aa[i + a->lda * j])) { 3018 x[i] = aa[i + a->lda * j]; 3019 if (idx) idx[i] = j; 3020 } 3021 } 3022 } 3023 PetscCall(VecRestoreArray(v, &x)); 3024 PetscCall(MatDenseRestoreArrayRead(A, &aa)); 3025 PetscFunctionReturn(0); 3026 } 3027 3028 PetscErrorCode MatGetColumnVector_SeqDense(Mat A, Vec v, PetscInt col) 3029 { 3030 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3031 PetscScalar *x; 3032 const PetscScalar *aa; 3033 3034 PetscFunctionBegin; 3035 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3036 PetscCall(MatDenseGetArrayRead(A, &aa)); 3037 PetscCall(VecGetArray(v, &x)); 3038 PetscCall(PetscArraycpy(x, aa + col * a->lda, A->rmap->n)); 3039 PetscCall(VecRestoreArray(v, &x)); 3040 PetscCall(MatDenseRestoreArrayRead(A, &aa)); 3041 PetscFunctionReturn(0); 3042 } 3043 3044 PETSC_INTERN PetscErrorCode MatGetColumnReductions_SeqDense(Mat A, PetscInt type, PetscReal *reductions) 3045 { 3046 PetscInt i, j, m, n; 3047 const PetscScalar *a; 3048 3049 PetscFunctionBegin; 3050 PetscCall(MatGetSize(A, &m, &n)); 3051 PetscCall(PetscArrayzero(reductions, n)); 3052 PetscCall(MatDenseGetArrayRead(A, &a)); 3053 if (type == NORM_2) { 3054 for (i = 0; i < n; i++) { 3055 for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j] * a[j]); 3056 a += m; 3057 } 3058 } else if (type == NORM_1) { 3059 for (i = 0; i < n; i++) { 3060 for (j = 0; j < m; j++) reductions[i] += PetscAbsScalar(a[j]); 3061 a += m; 3062 } 3063 } else if (type == NORM_INFINITY) { 3064 for (i = 0; i < n; i++) { 3065 for (j = 0; j < m; j++) reductions[i] = PetscMax(PetscAbsScalar(a[j]), reductions[i]); 3066 a += m; 3067 } 3068 } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) { 3069 for (i = 0; i < n; i++) { 3070 for (j = 0; j < m; j++) reductions[i] += PetscRealPart(a[j]); 3071 a += m; 3072 } 3073 } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) { 3074 for (i = 0; i < n; i++) { 3075 for (j = 0; j < m; j++) reductions[i] += PetscImaginaryPart(a[j]); 3076 a += m; 3077 } 3078 } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown reduction type"); 3079 PetscCall(MatDenseRestoreArrayRead(A, &a)); 3080 if (type == NORM_2) { 3081 for (i = 0; i < n; i++) reductions[i] = PetscSqrtReal(reductions[i]); 3082 } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) { 3083 for (i = 0; i < n; i++) reductions[i] /= m; 3084 } 3085 PetscFunctionReturn(0); 3086 } 3087 3088 PetscErrorCode MatSetRandom_SeqDense(Mat x, PetscRandom rctx) 3089 { 3090 PetscScalar *a; 3091 PetscInt lda, m, n, i, j; 3092 3093 PetscFunctionBegin; 3094 PetscCall(MatGetSize(x, &m, &n)); 3095 PetscCall(MatDenseGetLDA(x, &lda)); 3096 PetscCall(MatDenseGetArrayWrite(x, &a)); 3097 for (j = 0; j < n; j++) { 3098 for (i = 0; i < m; i++) PetscCall(PetscRandomGetValue(rctx, a + j * lda + i)); 3099 } 3100 PetscCall(MatDenseRestoreArrayWrite(x, &a)); 3101 PetscFunctionReturn(0); 3102 } 3103 3104 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A, PetscBool *missing, PetscInt *d) 3105 { 3106 PetscFunctionBegin; 3107 *missing = PETSC_FALSE; 3108 PetscFunctionReturn(0); 3109 } 3110 3111 /* vals is not const */ 3112 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A, PetscInt col, PetscScalar **vals) 3113 { 3114 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3115 PetscScalar *v; 3116 3117 PetscFunctionBegin; 3118 PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix"); 3119 PetscCall(MatDenseGetArray(A, &v)); 3120 *vals = v + col * a->lda; 3121 PetscCall(MatDenseRestoreArray(A, &v)); 3122 PetscFunctionReturn(0); 3123 } 3124 3125 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A, PetscScalar **vals) 3126 { 3127 PetscFunctionBegin; 3128 if (vals) *vals = NULL; /* user cannot accidentally use the array later */ 3129 PetscFunctionReturn(0); 3130 } 3131 3132 /* -------------------------------------------------------------------*/ 3133 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 3134 MatGetRow_SeqDense, 3135 MatRestoreRow_SeqDense, 3136 MatMult_SeqDense, 3137 /* 4*/ MatMultAdd_SeqDense, 3138 MatMultTranspose_SeqDense, 3139 MatMultTransposeAdd_SeqDense, 3140 NULL, 3141 NULL, 3142 NULL, 3143 /* 10*/ NULL, 3144 MatLUFactor_SeqDense, 3145 MatCholeskyFactor_SeqDense, 3146 MatSOR_SeqDense, 3147 MatTranspose_SeqDense, 3148 /* 15*/ MatGetInfo_SeqDense, 3149 MatEqual_SeqDense, 3150 MatGetDiagonal_SeqDense, 3151 MatDiagonalScale_SeqDense, 3152 MatNorm_SeqDense, 3153 /* 20*/ MatAssemblyBegin_SeqDense, 3154 MatAssemblyEnd_SeqDense, 3155 MatSetOption_SeqDense, 3156 MatZeroEntries_SeqDense, 3157 /* 24*/ MatZeroRows_SeqDense, 3158 NULL, 3159 NULL, 3160 NULL, 3161 NULL, 3162 /* 29*/ MatSetUp_SeqDense, 3163 NULL, 3164 NULL, 3165 NULL, 3166 NULL, 3167 /* 34*/ MatDuplicate_SeqDense, 3168 NULL, 3169 NULL, 3170 NULL, 3171 NULL, 3172 /* 39*/ MatAXPY_SeqDense, 3173 MatCreateSubMatrices_SeqDense, 3174 NULL, 3175 MatGetValues_SeqDense, 3176 MatCopy_SeqDense, 3177 /* 44*/ MatGetRowMax_SeqDense, 3178 MatScale_SeqDense, 3179 MatShift_SeqDense, 3180 NULL, 3181 MatZeroRowsColumns_SeqDense, 3182 /* 49*/ MatSetRandom_SeqDense, 3183 NULL, 3184 NULL, 3185 NULL, 3186 NULL, 3187 /* 54*/ NULL, 3188 NULL, 3189 NULL, 3190 NULL, 3191 NULL, 3192 /* 59*/ MatCreateSubMatrix_SeqDense, 3193 MatDestroy_SeqDense, 3194 MatView_SeqDense, 3195 NULL, 3196 NULL, 3197 /* 64*/ NULL, 3198 NULL, 3199 NULL, 3200 NULL, 3201 NULL, 3202 /* 69*/ MatGetRowMaxAbs_SeqDense, 3203 NULL, 3204 NULL, 3205 NULL, 3206 NULL, 3207 /* 74*/ NULL, 3208 NULL, 3209 NULL, 3210 NULL, 3211 NULL, 3212 /* 79*/ NULL, 3213 NULL, 3214 NULL, 3215 NULL, 3216 /* 83*/ MatLoad_SeqDense, 3217 MatIsSymmetric_SeqDense, 3218 MatIsHermitian_SeqDense, 3219 NULL, 3220 NULL, 3221 NULL, 3222 /* 89*/ NULL, 3223 NULL, 3224 MatMatMultNumeric_SeqDense_SeqDense, 3225 NULL, 3226 NULL, 3227 /* 94*/ NULL, 3228 NULL, 3229 NULL, 3230 MatMatTransposeMultNumeric_SeqDense_SeqDense, 3231 NULL, 3232 /* 99*/ MatProductSetFromOptions_SeqDense, 3233 NULL, 3234 NULL, 3235 MatConjugate_SeqDense, 3236 NULL, 3237 /*104*/ NULL, 3238 MatRealPart_SeqDense, 3239 MatImaginaryPart_SeqDense, 3240 NULL, 3241 NULL, 3242 /*109*/ NULL, 3243 NULL, 3244 MatGetRowMin_SeqDense, 3245 MatGetColumnVector_SeqDense, 3246 MatMissingDiagonal_SeqDense, 3247 /*114*/ NULL, 3248 NULL, 3249 NULL, 3250 NULL, 3251 NULL, 3252 /*119*/ NULL, 3253 NULL, 3254 NULL, 3255 NULL, 3256 NULL, 3257 /*124*/ NULL, 3258 MatGetColumnReductions_SeqDense, 3259 NULL, 3260 NULL, 3261 NULL, 3262 /*129*/ NULL, 3263 NULL, 3264 NULL, 3265 MatTransposeMatMultNumeric_SeqDense_SeqDense, 3266 NULL, 3267 /*134*/ NULL, 3268 NULL, 3269 NULL, 3270 NULL, 3271 NULL, 3272 /*139*/ NULL, 3273 NULL, 3274 NULL, 3275 NULL, 3276 NULL, 3277 MatCreateMPIMatConcatenateSeqMat_SeqDense, 3278 /*145*/ NULL, 3279 NULL, 3280 NULL, 3281 NULL, 3282 NULL, 3283 /*150*/ NULL, 3284 NULL}; 3285 3286 /*@C 3287 MatCreateSeqDense - Creates a `MATSEQDENSE` that 3288 is stored in column major order (the usual Fortran 77 manner). Many 3289 of the matrix operations use the BLAS and LAPACK routines. 3290 3291 Collective 3292 3293 Input Parameters: 3294 + comm - MPI communicator, set to `PETSC_COMM_SELF` 3295 . m - number of rows 3296 . n - number of columns 3297 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 3298 to control all matrix memory allocation. 3299 3300 Output Parameter: 3301 . A - the matrix 3302 3303 Note: 3304 The data input variable is intended primarily for Fortran programmers 3305 who wish to allocate their own matrix memory space. Most users should 3306 set data=NULL. 3307 3308 Level: intermediate 3309 3310 .seealso: `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()` 3311 @*/ 3312 PetscErrorCode MatCreateSeqDense(MPI_Comm comm, PetscInt m, PetscInt n, PetscScalar *data, Mat *A) 3313 { 3314 PetscFunctionBegin; 3315 PetscCall(MatCreate(comm, A)); 3316 PetscCall(MatSetSizes(*A, m, n, m, n)); 3317 PetscCall(MatSetType(*A, MATSEQDENSE)); 3318 PetscCall(MatSeqDenseSetPreallocation(*A, data)); 3319 PetscFunctionReturn(0); 3320 } 3321 3322 /*@C 3323 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements of a `MATSEQDENSE` matrix 3324 3325 Collective 3326 3327 Input Parameters: 3328 + B - the matrix 3329 - data - the array (or NULL) 3330 3331 Note: 3332 The data input variable is intended primarily for Fortran programmers 3333 who wish to allocate their own matrix memory space. Most users should 3334 need not call this routine. 3335 3336 Level: intermediate 3337 3338 .seealso: `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`, `MatDenseSetLDA()` 3339 @*/ 3340 PetscErrorCode MatSeqDenseSetPreallocation(Mat B, PetscScalar data[]) 3341 { 3342 PetscFunctionBegin; 3343 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 3344 PetscTryMethod(B, "MatSeqDenseSetPreallocation_C", (Mat, PetscScalar[]), (B, data)); 3345 PetscFunctionReturn(0); 3346 } 3347 3348 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B, PetscScalar *data) 3349 { 3350 Mat_SeqDense *b = (Mat_SeqDense *)B->data; 3351 3352 PetscFunctionBegin; 3353 PetscCheck(!b->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 3354 B->preallocated = PETSC_TRUE; 3355 3356 PetscCall(PetscLayoutSetUp(B->rmap)); 3357 PetscCall(PetscLayoutSetUp(B->cmap)); 3358 3359 if (b->lda <= 0) b->lda = B->rmap->n; 3360 3361 if (!data) { /* petsc-allocated storage */ 3362 if (!b->user_alloc) PetscCall(PetscFree(b->v)); 3363 PetscCall(PetscCalloc1((size_t)b->lda * B->cmap->n, &b->v)); 3364 3365 b->user_alloc = PETSC_FALSE; 3366 } else { /* user-allocated storage */ 3367 if (!b->user_alloc) PetscCall(PetscFree(b->v)); 3368 b->v = data; 3369 b->user_alloc = PETSC_TRUE; 3370 } 3371 B->assembled = PETSC_TRUE; 3372 PetscFunctionReturn(0); 3373 } 3374 3375 #if defined(PETSC_HAVE_ELEMENTAL) 3376 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 3377 { 3378 Mat mat_elemental; 3379 const PetscScalar *array; 3380 PetscScalar *v_colwise; 3381 PetscInt M = A->rmap->N, N = A->cmap->N, i, j, k, *rows, *cols; 3382 3383 PetscFunctionBegin; 3384 PetscCall(PetscMalloc3(M * N, &v_colwise, M, &rows, N, &cols)); 3385 PetscCall(MatDenseGetArrayRead(A, &array)); 3386 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 3387 k = 0; 3388 for (j = 0; j < N; j++) { 3389 cols[j] = j; 3390 for (i = 0; i < M; i++) v_colwise[j * M + i] = array[k++]; 3391 } 3392 for (i = 0; i < M; i++) rows[i] = i; 3393 PetscCall(MatDenseRestoreArrayRead(A, &array)); 3394 3395 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental)); 3396 PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N)); 3397 PetscCall(MatSetType(mat_elemental, MATELEMENTAL)); 3398 PetscCall(MatSetUp(mat_elemental)); 3399 3400 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 3401 PetscCall(MatSetValues(mat_elemental, M, rows, N, cols, v_colwise, ADD_VALUES)); 3402 PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY)); 3403 PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY)); 3404 PetscCall(PetscFree3(v_colwise, rows, cols)); 3405 3406 if (reuse == MAT_INPLACE_MATRIX) { 3407 PetscCall(MatHeaderReplace(A, &mat_elemental)); 3408 } else { 3409 *newmat = mat_elemental; 3410 } 3411 PetscFunctionReturn(0); 3412 } 3413 #endif 3414 3415 PetscErrorCode MatDenseSetLDA_SeqDense(Mat B, PetscInt lda) 3416 { 3417 Mat_SeqDense *b = (Mat_SeqDense *)B->data; 3418 PetscBool data; 3419 3420 PetscFunctionBegin; 3421 data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE); 3422 PetscCheck(b->user_alloc || !data || b->lda == lda, PETSC_COMM_SELF, PETSC_ERR_ORDER, "LDA cannot be changed after allocation of internal storage"); 3423 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); 3424 b->lda = lda; 3425 PetscFunctionReturn(0); 3426 } 3427 3428 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat) 3429 { 3430 PetscFunctionBegin; 3431 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIDense(comm, inmat, n, scall, outmat)); 3432 PetscFunctionReturn(0); 3433 } 3434 3435 PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A, PetscInt col, Vec *v) 3436 { 3437 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3438 3439 PetscFunctionBegin; 3440 PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first"); 3441 PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 3442 if (!a->cvec) { PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, &a->cvec)); } 3443 a->vecinuse = col + 1; 3444 PetscCall(MatDenseGetArray(A, (PetscScalar **)&a->ptrinuse)); 3445 PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda)); 3446 *v = a->cvec; 3447 PetscFunctionReturn(0); 3448 } 3449 3450 PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A, PetscInt col, Vec *v) 3451 { 3452 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3453 3454 PetscFunctionBegin; 3455 PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first"); 3456 PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector"); 3457 a->vecinuse = 0; 3458 PetscCall(MatDenseRestoreArray(A, (PetscScalar **)&a->ptrinuse)); 3459 PetscCall(VecResetArray(a->cvec)); 3460 if (v) *v = NULL; 3461 PetscFunctionReturn(0); 3462 } 3463 3464 PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v) 3465 { 3466 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3467 3468 PetscFunctionBegin; 3469 PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first"); 3470 PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 3471 if (!a->cvec) { PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, &a->cvec)); } 3472 a->vecinuse = col + 1; 3473 PetscCall(MatDenseGetArrayRead(A, &a->ptrinuse)); 3474 PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda)); 3475 PetscCall(VecLockReadPush(a->cvec)); 3476 *v = a->cvec; 3477 PetscFunctionReturn(0); 3478 } 3479 3480 PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v) 3481 { 3482 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3483 3484 PetscFunctionBegin; 3485 PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first"); 3486 PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector"); 3487 a->vecinuse = 0; 3488 PetscCall(MatDenseRestoreArrayRead(A, &a->ptrinuse)); 3489 PetscCall(VecLockReadPop(a->cvec)); 3490 PetscCall(VecResetArray(a->cvec)); 3491 if (v) *v = NULL; 3492 PetscFunctionReturn(0); 3493 } 3494 3495 PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v) 3496 { 3497 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3498 3499 PetscFunctionBegin; 3500 PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first"); 3501 PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 3502 if (!a->cvec) { PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, &a->cvec)); } 3503 a->vecinuse = col + 1; 3504 PetscCall(MatDenseGetArrayWrite(A, (PetscScalar **)&a->ptrinuse)); 3505 PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda)); 3506 *v = a->cvec; 3507 PetscFunctionReturn(0); 3508 } 3509 3510 PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v) 3511 { 3512 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3513 3514 PetscFunctionBegin; 3515 PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first"); 3516 PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector"); 3517 a->vecinuse = 0; 3518 PetscCall(MatDenseRestoreArrayWrite(A, (PetscScalar **)&a->ptrinuse)); 3519 PetscCall(VecResetArray(a->cvec)); 3520 if (v) *v = NULL; 3521 PetscFunctionReturn(0); 3522 } 3523 3524 PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v) 3525 { 3526 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3527 3528 PetscFunctionBegin; 3529 PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first"); 3530 PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 3531 if (a->cmat && (cend - cbegin != a->cmat->cmap->N || rend - rbegin != a->cmat->rmap->N)) PetscCall(MatDestroy(&a->cmat)); 3532 if (!a->cmat) { 3533 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), rend - rbegin, PETSC_DECIDE, rend - rbegin, cend - cbegin, a->v + rbegin + (size_t)cbegin * a->lda, &a->cmat)); 3534 } else { 3535 PetscCall(MatDensePlaceArray(a->cmat, a->v + rbegin + (size_t)cbegin * a->lda)); 3536 } 3537 PetscCall(MatDenseSetLDA(a->cmat, a->lda)); 3538 a->matinuse = cbegin + 1; 3539 *v = a->cmat; 3540 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP) 3541 A->offloadmask = PETSC_OFFLOAD_CPU; 3542 #endif 3543 PetscFunctionReturn(0); 3544 } 3545 3546 PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A, Mat *v) 3547 { 3548 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3549 3550 PetscFunctionBegin; 3551 PetscCheck(a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetSubMatrix() first"); 3552 PetscCheck(a->cmat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column matrix"); 3553 PetscCheck(*v == a->cmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not the matrix obtained from MatDenseGetSubMatrix()"); 3554 a->matinuse = 0; 3555 PetscCall(MatDenseResetArray(a->cmat)); 3556 if (v) *v = NULL; 3557 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP) 3558 A->offloadmask = PETSC_OFFLOAD_CPU; 3559 #endif 3560 PetscFunctionReturn(0); 3561 } 3562 3563 /*MC 3564 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 3565 3566 Options Database Keys: 3567 . -mat_type seqdense - sets the matrix type to `MATSEQDENSE` during a call to `MatSetFromOptions()` 3568 3569 Level: beginner 3570 3571 .seealso: `MATSEQDENSE`, `MatCreateSeqDense()` 3572 M*/ 3573 PetscErrorCode MatCreate_SeqDense(Mat B) 3574 { 3575 Mat_SeqDense *b; 3576 PetscMPIInt size; 3577 3578 PetscFunctionBegin; 3579 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size)); 3580 PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1"); 3581 3582 PetscCall(PetscNew(&b)); 3583 PetscCall(PetscMemcpy(B->ops, &MatOps_Values, sizeof(struct _MatOps))); 3584 B->data = (void *)b; 3585 3586 b->roworiented = PETSC_TRUE; 3587 3588 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactor_C", MatQRFactor_SeqDense)); 3589 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetLDA_C", MatDenseGetLDA_SeqDense)); 3590 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseSetLDA_C", MatDenseSetLDA_SeqDense)); 3591 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArray_C", MatDenseGetArray_SeqDense)); 3592 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArray_C", MatDenseRestoreArray_SeqDense)); 3593 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDensePlaceArray_C", MatDensePlaceArray_SeqDense)); 3594 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseResetArray_C", MatDenseResetArray_SeqDense)); 3595 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseReplaceArray_C", MatDenseReplaceArray_SeqDense)); 3596 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayRead_C", MatDenseGetArray_SeqDense)); 3597 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayRead_C", MatDenseRestoreArray_SeqDense)); 3598 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayWrite_C", MatDenseGetArray_SeqDense)); 3599 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayWrite_C", MatDenseRestoreArray_SeqDense)); 3600 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqaij_C", MatConvert_SeqDense_SeqAIJ)); 3601 #if defined(PETSC_HAVE_ELEMENTAL) 3602 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_elemental_C", MatConvert_SeqDense_Elemental)); 3603 #endif 3604 #if defined(PETSC_HAVE_SCALAPACK) 3605 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_scalapack_C", MatConvert_Dense_ScaLAPACK)); 3606 #endif 3607 #if defined(PETSC_HAVE_CUDA) 3608 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensecuda_C", MatConvert_SeqDense_SeqDenseCUDA)); 3609 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", MatProductSetFromOptions_SeqDense)); 3610 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdense_C", MatProductSetFromOptions_SeqDense)); 3611 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensecuda_C", MatProductSetFromOptions_SeqDense)); 3612 #endif 3613 #if defined(PETSC_HAVE_HIP) 3614 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensehip_C", MatConvert_SeqDense_SeqDenseHIP)); 3615 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", MatProductSetFromOptions_SeqDense)); 3616 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdense_C", MatProductSetFromOptions_SeqDense)); 3617 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensehip_C", MatProductSetFromOptions_SeqDense)); 3618 #endif 3619 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqDenseSetPreallocation_C", MatSeqDenseSetPreallocation_SeqDense)); 3620 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqdense_C", MatProductSetFromOptions_SeqAIJ_SeqDense)); 3621 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdense_C", MatProductSetFromOptions_SeqDense)); 3622 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense)); 3623 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqsbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense)); 3624 3625 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumn_C", MatDenseGetColumn_SeqDense)); 3626 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumn_C", MatDenseRestoreColumn_SeqDense)); 3627 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVec_C", MatDenseGetColumnVec_SeqDense)); 3628 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVec_C", MatDenseRestoreColumnVec_SeqDense)); 3629 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecRead_C", MatDenseGetColumnVecRead_SeqDense)); 3630 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecRead_C", MatDenseRestoreColumnVecRead_SeqDense)); 3631 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecWrite_C", MatDenseGetColumnVecWrite_SeqDense)); 3632 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecWrite_C", MatDenseRestoreColumnVecWrite_SeqDense)); 3633 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetSubMatrix_C", MatDenseGetSubMatrix_SeqDense)); 3634 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreSubMatrix_C", MatDenseRestoreSubMatrix_SeqDense)); 3635 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQDENSE)); 3636 PetscFunctionReturn(0); 3637 } 3638 3639 /*@C 3640 MatDenseGetColumn - gives access to a column of a dense matrix. This is only the local part of the column. You MUST call `MatDenseRestoreColumn()` to avoid memory bleeding. 3641 3642 Not Collective 3643 3644 Input Parameters: 3645 + mat - a `MATSEQDENSE` or `MATMPIDENSE` matrix 3646 - col - column index 3647 3648 Output Parameter: 3649 . vals - pointer to the data 3650 3651 Level: intermediate 3652 3653 Note: 3654 Use `MatDenseGetColumnVec()` to get access to a column of a `MATDENSE` treated as a `Vec` 3655 3656 .seealso: `MATDENSE`, `MatDenseRestoreColumn()`, `MatDenseGetColumnVec()` 3657 @*/ 3658 PetscErrorCode MatDenseGetColumn(Mat A, PetscInt col, PetscScalar **vals) 3659 { 3660 PetscFunctionBegin; 3661 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3662 PetscValidLogicalCollectiveInt(A, col, 2); 3663 PetscValidPointer(vals, 3); 3664 PetscUseMethod(A, "MatDenseGetColumn_C", (Mat, PetscInt, PetscScalar **), (A, col, vals)); 3665 PetscFunctionReturn(0); 3666 } 3667 3668 /*@C 3669 MatDenseRestoreColumn - returns access to a column of a `MATDENSE` matrix which is returned by `MatDenseGetColumn()`. 3670 3671 Not Collective 3672 3673 Input Parameters: 3674 + mat - a `MATSEQDENSE` or `MATMPIDENSE` matrix 3675 - vals - pointer to the data (may be NULL) 3676 3677 Level: intermediate 3678 3679 .seealso: `MATDENSE`, `MatDenseGetColumn()` 3680 @*/ 3681 PetscErrorCode MatDenseRestoreColumn(Mat A, PetscScalar **vals) 3682 { 3683 PetscFunctionBegin; 3684 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3685 PetscValidPointer(vals, 2); 3686 PetscUseMethod(A, "MatDenseRestoreColumn_C", (Mat, PetscScalar **), (A, vals)); 3687 PetscFunctionReturn(0); 3688 } 3689 3690 /*@ 3691 MatDenseGetColumnVec - Gives read-write access to a column of a `MATDENSE` matrix, represented as a `Vec`. 3692 3693 Collective 3694 3695 Input Parameters: 3696 + mat - the `Mat` object 3697 - col - the column index 3698 3699 Output Parameter: 3700 . v - the vector 3701 3702 Notes: 3703 The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVec()` when the vector is no longer needed. 3704 3705 Use `MatDenseGetColumnVecRead()` to obtain read-only access or `MatDenseGetColumnVecWrite()` for write-only access. 3706 3707 Level: intermediate 3708 3709 .seealso: `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`, `MatDenseGetColumn()` 3710 @*/ 3711 PetscErrorCode MatDenseGetColumnVec(Mat A, PetscInt col, Vec *v) 3712 { 3713 PetscFunctionBegin; 3714 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3715 PetscValidType(A, 1); 3716 PetscValidLogicalCollectiveInt(A, col, 2); 3717 PetscValidPointer(v, 3); 3718 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3719 PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N); 3720 PetscUseMethod(A, "MatDenseGetColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v)); 3721 PetscFunctionReturn(0); 3722 } 3723 3724 /*@ 3725 MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec(). 3726 3727 Collective 3728 3729 Input Parameters: 3730 + mat - the Mat object 3731 . col - the column index 3732 - v - the Vec object (may be NULL) 3733 3734 Level: intermediate 3735 3736 .seealso: `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3737 @*/ 3738 PetscErrorCode MatDenseRestoreColumnVec(Mat A, PetscInt col, Vec *v) 3739 { 3740 PetscFunctionBegin; 3741 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3742 PetscValidType(A, 1); 3743 PetscValidLogicalCollectiveInt(A, col, 2); 3744 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3745 PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N); 3746 PetscUseMethod(A, "MatDenseRestoreColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v)); 3747 PetscFunctionReturn(0); 3748 } 3749 3750 /*@ 3751 MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec. 3752 3753 Collective 3754 3755 Input Parameters: 3756 + mat - the Mat object 3757 - col - the column index 3758 3759 Output Parameter: 3760 . v - the vector 3761 3762 Notes: 3763 The vector is owned by PETSc and users cannot modify it. 3764 3765 Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed. 3766 3767 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access. 3768 3769 Level: intermediate 3770 3771 .seealso: `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3772 @*/ 3773 PetscErrorCode MatDenseGetColumnVecRead(Mat A, PetscInt col, Vec *v) 3774 { 3775 PetscFunctionBegin; 3776 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3777 PetscValidType(A, 1); 3778 PetscValidLogicalCollectiveInt(A, col, 2); 3779 PetscValidPointer(v, 3); 3780 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3781 PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N); 3782 PetscUseMethod(A, "MatDenseGetColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v)); 3783 PetscFunctionReturn(0); 3784 } 3785 3786 /*@ 3787 MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead(). 3788 3789 Collective 3790 3791 Input Parameters: 3792 + mat - the Mat object 3793 . col - the column index 3794 - v - the Vec object (may be NULL) 3795 3796 Level: intermediate 3797 3798 .seealso: `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecWrite()` 3799 @*/ 3800 PetscErrorCode MatDenseRestoreColumnVecRead(Mat A, PetscInt col, Vec *v) 3801 { 3802 PetscFunctionBegin; 3803 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3804 PetscValidType(A, 1); 3805 PetscValidLogicalCollectiveInt(A, col, 2); 3806 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3807 PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N); 3808 PetscUseMethod(A, "MatDenseRestoreColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v)); 3809 PetscFunctionReturn(0); 3810 } 3811 3812 /*@ 3813 MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec. 3814 3815 Collective 3816 3817 Input Parameters: 3818 + mat - the Mat object 3819 - col - the column index 3820 3821 Output Parameter: 3822 . v - the vector 3823 3824 Notes: 3825 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed. 3826 3827 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access. 3828 3829 Level: intermediate 3830 3831 .seealso: `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3832 @*/ 3833 PetscErrorCode MatDenseGetColumnVecWrite(Mat A, PetscInt col, Vec *v) 3834 { 3835 PetscFunctionBegin; 3836 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3837 PetscValidType(A, 1); 3838 PetscValidLogicalCollectiveInt(A, col, 2); 3839 PetscValidPointer(v, 3); 3840 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3841 PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N); 3842 PetscUseMethod(A, "MatDenseGetColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v)); 3843 PetscFunctionReturn(0); 3844 } 3845 3846 /*@ 3847 MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite(). 3848 3849 Collective 3850 3851 Input Parameters: 3852 + mat - the Mat object 3853 . col - the column index 3854 - v - the Vec object (may be NULL) 3855 3856 Level: intermediate 3857 3858 .seealso: `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()` 3859 @*/ 3860 PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A, PetscInt col, Vec *v) 3861 { 3862 PetscFunctionBegin; 3863 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3864 PetscValidType(A, 1); 3865 PetscValidLogicalCollectiveInt(A, col, 2); 3866 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3867 PetscCheck(col >= 0 && col < A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid col %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT ")", col, A->cmap->N); 3868 PetscUseMethod(A, "MatDenseRestoreColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v)); 3869 PetscFunctionReturn(0); 3870 } 3871 3872 /*@ 3873 MatDenseGetSubMatrix - Gives access to a block of rows and columns of a dense matrix, represented as a Mat. 3874 3875 Collective 3876 3877 Input Parameters: 3878 + mat - the Mat object 3879 . rbegin - the first global row index in the block (if PETSC_DECIDE, is 0) 3880 . rend - the global row index past the last one in the block (if PETSC_DECIDE, is M) 3881 . cbegin - the first global column index in the block (if PETSC_DECIDE, is 0) 3882 - cend - the global column index past the last one in the block (if PETSC_DECIDE, is N) 3883 3884 Output Parameter: 3885 . v - the matrix 3886 3887 Notes: 3888 The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed. 3889 3890 The output matrix is not redistributed by PETSc, so depending on the values of rbegin and rend, some processes may have no local rows. 3891 3892 Level: intermediate 3893 3894 .seealso: `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreSubMatrix()` 3895 @*/ 3896 PetscErrorCode MatDenseGetSubMatrix(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v) 3897 { 3898 PetscFunctionBegin; 3899 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3900 PetscValidType(A, 1); 3901 PetscValidLogicalCollectiveInt(A, rbegin, 2); 3902 PetscValidLogicalCollectiveInt(A, rend, 3); 3903 PetscValidLogicalCollectiveInt(A, cbegin, 4); 3904 PetscValidLogicalCollectiveInt(A, cend, 5); 3905 PetscValidPointer(v, 6); 3906 if (rbegin == PETSC_DECIDE) rbegin = 0; 3907 if (rend == PETSC_DECIDE) rend = A->rmap->N; 3908 if (cbegin == PETSC_DECIDE) cbegin = 0; 3909 if (cend == PETSC_DECIDE) cend = A->cmap->N; 3910 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3911 PetscCheck(rbegin >= 0 && rbegin <= A->rmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid rbegin %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT "]", rbegin, A->rmap->N); 3912 PetscCheck(rend >= rbegin && rend <= A->rmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid rend %" PetscInt_FMT ", should be in [%" PetscInt_FMT ",%" PetscInt_FMT "]", rend, rbegin, A->rmap->N); 3913 PetscCheck(cbegin >= 0 && cbegin <= A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid cbegin %" PetscInt_FMT ", should be in [0,%" PetscInt_FMT "]", cbegin, A->cmap->N); 3914 PetscCheck(cend >= cbegin && cend <= A->cmap->N, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid cend %" PetscInt_FMT ", should be in [%" PetscInt_FMT ",%" PetscInt_FMT "]", cend, cbegin, A->cmap->N); 3915 PetscUseMethod(A, "MatDenseGetSubMatrix_C", (Mat, PetscInt, PetscInt, PetscInt, PetscInt, Mat *), (A, rbegin, rend, cbegin, cend, v)); 3916 PetscFunctionReturn(0); 3917 } 3918 3919 /*@ 3920 MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix(). 3921 3922 Collective 3923 3924 Input Parameters: 3925 + mat - the Mat object 3926 - v - the Mat object (may be NULL) 3927 3928 Level: intermediate 3929 3930 .seealso: `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseGetSubMatrix()` 3931 @*/ 3932 PetscErrorCode MatDenseRestoreSubMatrix(Mat A, Mat *v) 3933 { 3934 PetscFunctionBegin; 3935 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3936 PetscValidType(A, 1); 3937 PetscValidPointer(v, 2); 3938 PetscUseMethod(A, "MatDenseRestoreSubMatrix_C", (Mat, Mat *), (A, v)); 3939 PetscFunctionReturn(0); 3940 } 3941 3942 #include <petscblaslapack.h> 3943 #include <petsc/private/kernels/blockinvert.h> 3944 3945 PetscErrorCode MatSeqDenseInvert(Mat A) 3946 { 3947 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3948 PetscInt bs = A->rmap->n; 3949 MatScalar *values = a->v; 3950 const PetscReal shift = 0.0; 3951 PetscBool allowzeropivot = PetscNot(A->erroriffailure), zeropivotdetected = PETSC_FALSE; 3952 3953 PetscFunctionBegin; 3954 /* factor and invert each block */ 3955 switch (bs) { 3956 case 1: 3957 values[0] = (PetscScalar)1.0 / (values[0] + shift); 3958 break; 3959 case 2: 3960 PetscCall(PetscKernel_A_gets_inverse_A_2(values, shift, allowzeropivot, &zeropivotdetected)); 3961 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3962 break; 3963 case 3: 3964 PetscCall(PetscKernel_A_gets_inverse_A_3(values, shift, allowzeropivot, &zeropivotdetected)); 3965 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3966 break; 3967 case 4: 3968 PetscCall(PetscKernel_A_gets_inverse_A_4(values, shift, allowzeropivot, &zeropivotdetected)); 3969 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3970 break; 3971 case 5: { 3972 PetscScalar work[25]; 3973 PetscInt ipvt[5]; 3974 3975 PetscCall(PetscKernel_A_gets_inverse_A_5(values, ipvt, work, shift, allowzeropivot, &zeropivotdetected)); 3976 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3977 } break; 3978 case 6: 3979 PetscCall(PetscKernel_A_gets_inverse_A_6(values, shift, allowzeropivot, &zeropivotdetected)); 3980 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3981 break; 3982 case 7: 3983 PetscCall(PetscKernel_A_gets_inverse_A_7(values, shift, allowzeropivot, &zeropivotdetected)); 3984 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3985 break; 3986 default: { 3987 PetscInt *v_pivots, *IJ, j; 3988 PetscScalar *v_work; 3989 3990 PetscCall(PetscMalloc3(bs, &v_work, bs, &v_pivots, bs, &IJ)); 3991 for (j = 0; j < bs; j++) IJ[j] = j; 3992 PetscCall(PetscKernel_A_gets_inverse_A(bs, values, v_pivots, v_work, allowzeropivot, &zeropivotdetected)); 3993 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3994 PetscCall(PetscFree3(v_work, v_pivots, IJ)); 3995 } 3996 } 3997 PetscFunctionReturn(0); 3998 } 3999