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(PETSC_SUCCESS); 32 } 33 34 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(PETSC_SUCCESS); 41 PetscCall(PetscBLASIntCast(A->cmap->n, &n)); 42 if (A->factortype == MAT_FACTOR_LU) { 43 PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present"); 44 if (!mat->fwork) { 45 mat->lfwork = n; 46 PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork)); 47 } 48 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 49 PetscCallBLAS("LAPACKgetri", LAPACKgetri_(&n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info)); 50 PetscCall(PetscFPTrapPop()); 51 PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0)); 52 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 53 if (A->spd == PETSC_BOOL3_TRUE) { 54 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 55 PetscCallBLAS("LAPACKpotri", LAPACKpotri_("L", &n, mat->v, &mat->lda, &info)); 56 PetscCall(PetscFPTrapPop()); 57 PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_TRUE)); 58 #if defined(PETSC_USE_COMPLEX) 59 } else if (A->hermitian == PETSC_BOOL3_TRUE) { 60 PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present"); 61 PetscCheck(mat->fwork, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Fwork not present"); 62 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 63 PetscCallBLAS("LAPACKhetri", LAPACKhetri_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &info)); 64 PetscCall(PetscFPTrapPop()); 65 PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_TRUE)); 66 #endif 67 } else { /* symmetric case */ 68 PetscCheck(mat->pivots, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Pivots not present"); 69 PetscCheck(mat->fwork, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Fwork not present"); 70 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 71 PetscCallBLAS("LAPACKsytri", LAPACKsytri_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &info)); 72 PetscCall(PetscFPTrapPop()); 73 PetscCall(MatSeqDenseSymmetrize_Private(A, PETSC_FALSE)); 74 } 75 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad Inversion: zero pivot in row %" PetscInt_FMT, (PetscInt)info - 1); 76 PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0)); 77 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix must be factored to solve"); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 828 } 829 830 /* COMMENT: I have chosen to hide row permutation in the pivots, 831 rather than put it in the Mat->row slot.*/ 832 PetscErrorCode MatLUFactor_SeqDense(Mat A, IS row, IS col, const MatFactorInfo *minfo) 833 { 834 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 835 PetscBLASInt n, m, info; 836 837 PetscFunctionBegin; 838 PetscCall(PetscBLASIntCast(A->cmap->n, &n)); 839 PetscCall(PetscBLASIntCast(A->rmap->n, &m)); 840 if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); } 841 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS); 842 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 843 PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&m, &n, mat->v, &mat->lda, mat->pivots, &info)); 844 PetscCall(PetscFPTrapPop()); 845 846 PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to LU factorization %d", (int)info); 847 PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Bad LU factorization %d", (int)info); 848 849 A->ops->solve = MatSolve_SeqDense_LU; 850 A->ops->matsolve = MatMatSolve_SeqDense_LU; 851 A->ops->solvetranspose = MatSolveTranspose_SeqDense_LU; 852 A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_LU; 853 A->factortype = MAT_FACTOR_LU; 854 855 PetscCall(PetscFree(A->solvertype)); 856 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype)); 857 858 PetscCall(PetscLogFlops((2.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3)); 859 PetscFunctionReturn(PETSC_SUCCESS); 860 } 861 862 static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info_dummy) 863 { 864 MatFactorInfo info; 865 866 PetscFunctionBegin; 867 PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES)); 868 PetscUseTypeMethod(fact, lufactor, NULL, NULL, &info); 869 PetscFunctionReturn(PETSC_SUCCESS); 870 } 871 872 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, IS col, const MatFactorInfo *info) 873 { 874 PetscFunctionBegin; 875 fact->preallocated = PETSC_TRUE; 876 fact->assembled = PETSC_TRUE; 877 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 878 PetscFunctionReturn(PETSC_SUCCESS); 879 } 880 881 /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */ 882 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A, IS perm, const MatFactorInfo *factinfo) 883 { 884 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 885 PetscBLASInt info, n; 886 887 PetscFunctionBegin; 888 PetscCall(PetscBLASIntCast(A->cmap->n, &n)); 889 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS); 890 if (A->spd == PETSC_BOOL3_TRUE) { 891 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 892 PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("L", &n, mat->v, &mat->lda, &info)); 893 PetscCall(PetscFPTrapPop()); 894 #if defined(PETSC_USE_COMPLEX) 895 } else if (A->hermitian == PETSC_BOOL3_TRUE) { 896 if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); } 897 if (!mat->fwork) { 898 PetscScalar dummy; 899 900 mat->lfwork = -1; 901 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 902 PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info)); 903 PetscCall(PetscFPTrapPop()); 904 mat->lfwork = (PetscInt)PetscRealPart(dummy); 905 PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork)); 906 } 907 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 908 PetscCallBLAS("LAPACKhetrf", LAPACKhetrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info)); 909 PetscCall(PetscFPTrapPop()); 910 #endif 911 } else { /* symmetric case */ 912 if (!mat->pivots) { PetscCall(PetscMalloc1(A->rmap->n, &mat->pivots)); } 913 if (!mat->fwork) { 914 PetscScalar dummy; 915 916 mat->lfwork = -1; 917 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 918 PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, &dummy, &mat->lfwork, &info)); 919 PetscCall(PetscFPTrapPop()); 920 mat->lfwork = (PetscInt)PetscRealPart(dummy); 921 PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork)); 922 } 923 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 924 PetscCallBLAS("LAPACKsytrf", LAPACKsytrf_("L", &n, mat->v, &mat->lda, mat->pivots, mat->fwork, &mat->lfwork, &info)); 925 PetscCall(PetscFPTrapPop()); 926 } 927 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT, "Bad factorization: zero pivot in row %" PetscInt_FMT, (PetscInt)info - 1); 928 929 A->ops->solve = MatSolve_SeqDense_Cholesky; 930 A->ops->matsolve = MatMatSolve_SeqDense_Cholesky; 931 A->ops->solvetranspose = MatSolveTranspose_SeqDense_Cholesky; 932 A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_Cholesky; 933 A->factortype = MAT_FACTOR_CHOLESKY; 934 935 PetscCall(PetscFree(A->solvertype)); 936 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype)); 937 938 PetscCall(PetscLogFlops((1.0 * A->cmap->n * A->cmap->n * A->cmap->n) / 3.0)); 939 PetscFunctionReturn(PETSC_SUCCESS); 940 } 941 942 static PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info_dummy) 943 { 944 MatFactorInfo info; 945 946 PetscFunctionBegin; 947 info.fill = 1.0; 948 949 PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES)); 950 PetscUseTypeMethod(fact, choleskyfactor, NULL, &info); 951 PetscFunctionReturn(PETSC_SUCCESS); 952 } 953 954 PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info) 955 { 956 PetscFunctionBegin; 957 fact->assembled = PETSC_TRUE; 958 fact->preallocated = PETSC_TRUE; 959 fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense; 960 PetscFunctionReturn(PETSC_SUCCESS); 961 } 962 963 PetscErrorCode MatQRFactor_SeqDense(Mat A, IS col, const MatFactorInfo *minfo) 964 { 965 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 966 PetscBLASInt n, m, info, min, max; 967 968 PetscFunctionBegin; 969 PetscCall(PetscBLASIntCast(A->cmap->n, &n)); 970 PetscCall(PetscBLASIntCast(A->rmap->n, &m)); 971 max = PetscMax(m, n); 972 min = PetscMin(m, n); 973 if (!mat->tau) { PetscCall(PetscMalloc1(min, &mat->tau)); } 974 if (!mat->pivots) { PetscCall(PetscMalloc1(n, &mat->pivots)); } 975 if (!mat->qrrhs) PetscCall(MatCreateVecs(A, NULL, &(mat->qrrhs))); 976 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS); 977 if (!mat->fwork) { 978 PetscScalar dummy; 979 980 mat->lfwork = -1; 981 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 982 PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, &dummy, &mat->lfwork, &info)); 983 PetscCall(PetscFPTrapPop()); 984 mat->lfwork = (PetscInt)PetscRealPart(dummy); 985 PetscCall(PetscMalloc1(mat->lfwork, &mat->fwork)); 986 } 987 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 988 PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&m, &n, mat->v, &mat->lda, mat->tau, mat->fwork, &mat->lfwork, &info)); 989 PetscCall(PetscFPTrapPop()); 990 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to QR factorization %d", (int)info); 991 // 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 992 mat->rank = min; 993 994 A->ops->solve = MatSolve_SeqDense_QR; 995 A->ops->matsolve = MatMatSolve_SeqDense_QR; 996 A->factortype = MAT_FACTOR_QR; 997 if (m == n) { 998 A->ops->solvetranspose = MatSolveTranspose_SeqDense_QR; 999 A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_QR; 1000 } 1001 1002 PetscCall(PetscFree(A->solvertype)); 1003 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &A->solvertype)); 1004 1005 PetscCall(PetscLogFlops(2.0 * min * min * (max - min / 3.0))); 1006 PetscFunctionReturn(PETSC_SUCCESS); 1007 } 1008 1009 static PetscErrorCode MatQRFactorNumeric_SeqDense(Mat fact, Mat A, const MatFactorInfo *info_dummy) 1010 { 1011 MatFactorInfo info; 1012 1013 PetscFunctionBegin; 1014 info.fill = 1.0; 1015 1016 PetscCall(MatDuplicateNoCreate_SeqDense(fact, A, MAT_COPY_VALUES)); 1017 PetscUseMethod(fact, "MatQRFactor_C", (Mat, IS, const MatFactorInfo *), (fact, NULL, &info)); 1018 PetscFunctionReturn(PETSC_SUCCESS); 1019 } 1020 1021 PetscErrorCode MatQRFactorSymbolic_SeqDense(Mat fact, Mat A, IS row, const MatFactorInfo *info) 1022 { 1023 PetscFunctionBegin; 1024 fact->assembled = PETSC_TRUE; 1025 fact->preallocated = PETSC_TRUE; 1026 PetscCall(PetscObjectComposeFunction((PetscObject)fact, "MatQRFactorNumeric_C", MatQRFactorNumeric_SeqDense)); 1027 PetscFunctionReturn(PETSC_SUCCESS); 1028 } 1029 1030 /* uses LAPACK */ 1031 PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A, MatFactorType ftype, Mat *fact) 1032 { 1033 PetscFunctionBegin; 1034 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), fact)); 1035 PetscCall(MatSetSizes(*fact, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n)); 1036 PetscCall(MatSetType(*fact, MATDENSE)); 1037 (*fact)->trivialsymbolic = PETSC_TRUE; 1038 if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU) { 1039 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 1040 (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense; 1041 } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1042 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 1043 } else if (ftype == MAT_FACTOR_QR) { 1044 PetscCall(PetscObjectComposeFunction((PetscObject)(*fact), "MatQRFactorSymbolic_C", MatQRFactorSymbolic_SeqDense)); 1045 } 1046 (*fact)->factortype = ftype; 1047 1048 PetscCall(PetscFree((*fact)->solvertype)); 1049 PetscCall(PetscStrallocpy(MATSOLVERPETSC, &(*fact)->solvertype)); 1050 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_LU])); 1051 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ILU])); 1052 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_CHOLESKY])); 1053 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&(*fact)->preferredordering[MAT_FACTOR_ICC])); 1054 PetscFunctionReturn(PETSC_SUCCESS); 1055 } 1056 1057 static PetscErrorCode MatSOR_SeqDense(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal shift, PetscInt its, PetscInt lits, Vec xx) 1058 { 1059 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 1060 PetscScalar *x, *v = mat->v, zero = 0.0, xt; 1061 const PetscScalar *b; 1062 PetscInt m = A->rmap->n, i; 1063 PetscBLASInt o = 1, bm = 0; 1064 1065 PetscFunctionBegin; 1066 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP) 1067 PetscCheck(A->offloadmask != PETSC_OFFLOAD_GPU, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented"); 1068 #endif 1069 if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */ 1070 PetscCall(PetscBLASIntCast(m, &bm)); 1071 if (flag & SOR_ZERO_INITIAL_GUESS) { 1072 /* this is a hack fix, should have another version without the second BLASdotu */ 1073 PetscCall(VecSet(xx, zero)); 1074 } 1075 PetscCall(VecGetArray(xx, &x)); 1076 PetscCall(VecGetArrayRead(bb, &b)); 1077 its = its * lits; 1078 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); 1079 while (its--) { 1080 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1081 for (i = 0; i < m; i++) { 1082 PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o)); 1083 x[i] = (1. - omega) * x[i] + omega * (xt + v[i + i * m] * x[i]) / (v[i + i * m] + shift); 1084 } 1085 } 1086 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1087 for (i = m - 1; i >= 0; i--) { 1088 PetscCallBLAS("BLASdotu", xt = b[i] - BLASdotu_(&bm, v + i, &bm, x, &o)); 1089 x[i] = (1. - omega) * x[i] + omega * (xt + v[i + i * m] * x[i]) / (v[i + i * m] + shift); 1090 } 1091 } 1092 } 1093 PetscCall(VecRestoreArrayRead(bb, &b)); 1094 PetscCall(VecRestoreArray(xx, &x)); 1095 PetscFunctionReturn(PETSC_SUCCESS); 1096 } 1097 1098 PetscErrorCode MatMultTranspose_SeqDense(Mat A, Vec xx, Vec yy) 1099 { 1100 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 1101 const PetscScalar *v = mat->v, *x; 1102 PetscScalar *y; 1103 PetscBLASInt m, n, _One = 1; 1104 PetscScalar _DOne = 1.0, _DZero = 0.0; 1105 1106 PetscFunctionBegin; 1107 PetscCall(PetscBLASIntCast(A->rmap->n, &m)); 1108 PetscCall(PetscBLASIntCast(A->cmap->n, &n)); 1109 PetscCall(VecGetArrayRead(xx, &x)); 1110 PetscCall(VecGetArrayWrite(yy, &y)); 1111 if (!A->rmap->n || !A->cmap->n) { 1112 PetscBLASInt i; 1113 for (i = 0; i < n; i++) y[i] = 0.0; 1114 } else { 1115 PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v, &mat->lda, x, &_One, &_DZero, y, &_One)); 1116 PetscCall(PetscLogFlops(2.0 * A->rmap->n * A->cmap->n - A->cmap->n)); 1117 } 1118 PetscCall(VecRestoreArrayRead(xx, &x)); 1119 PetscCall(VecRestoreArrayWrite(yy, &y)); 1120 PetscFunctionReturn(PETSC_SUCCESS); 1121 } 1122 1123 PetscErrorCode MatMult_SeqDense(Mat A, Vec xx, Vec yy) 1124 { 1125 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 1126 PetscScalar *y, _DOne = 1.0, _DZero = 0.0; 1127 PetscBLASInt m, n, _One = 1; 1128 const PetscScalar *v = mat->v, *x; 1129 1130 PetscFunctionBegin; 1131 PetscCall(PetscBLASIntCast(A->rmap->n, &m)); 1132 PetscCall(PetscBLASIntCast(A->cmap->n, &n)); 1133 PetscCall(VecGetArrayRead(xx, &x)); 1134 PetscCall(VecGetArrayWrite(yy, &y)); 1135 if (!A->rmap->n || !A->cmap->n) { 1136 PetscBLASInt i; 1137 for (i = 0; i < m; i++) y[i] = 0.0; 1138 } else { 1139 PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v, &(mat->lda), x, &_One, &_DZero, y, &_One)); 1140 PetscCall(PetscLogFlops(2.0 * A->rmap->n * A->cmap->n - A->rmap->n)); 1141 } 1142 PetscCall(VecRestoreArrayRead(xx, &x)); 1143 PetscCall(VecRestoreArrayWrite(yy, &y)); 1144 PetscFunctionReturn(PETSC_SUCCESS); 1145 } 1146 1147 PetscErrorCode MatMultAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy) 1148 { 1149 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 1150 const PetscScalar *v = mat->v, *x; 1151 PetscScalar *y, _DOne = 1.0; 1152 PetscBLASInt m, n, _One = 1; 1153 1154 PetscFunctionBegin; 1155 PetscCall(PetscBLASIntCast(A->rmap->n, &m)); 1156 PetscCall(PetscBLASIntCast(A->cmap->n, &n)); 1157 PetscCall(VecCopy(zz, yy)); 1158 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS); 1159 PetscCall(VecGetArrayRead(xx, &x)); 1160 PetscCall(VecGetArray(yy, &y)); 1161 PetscCallBLAS("BLASgemv", BLASgemv_("N", &m, &n, &_DOne, v, &(mat->lda), x, &_One, &_DOne, y, &_One)); 1162 PetscCall(VecRestoreArrayRead(xx, &x)); 1163 PetscCall(VecRestoreArray(yy, &y)); 1164 PetscCall(PetscLogFlops(2.0 * A->rmap->n * A->cmap->n)); 1165 PetscFunctionReturn(PETSC_SUCCESS); 1166 } 1167 1168 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A, Vec xx, Vec zz, Vec yy) 1169 { 1170 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 1171 const PetscScalar *v = mat->v, *x; 1172 PetscScalar *y; 1173 PetscBLASInt m, n, _One = 1; 1174 PetscScalar _DOne = 1.0; 1175 1176 PetscFunctionBegin; 1177 PetscCall(PetscBLASIntCast(A->rmap->n, &m)); 1178 PetscCall(PetscBLASIntCast(A->cmap->n, &n)); 1179 PetscCall(VecCopy(zz, yy)); 1180 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(PETSC_SUCCESS); 1181 PetscCall(VecGetArrayRead(xx, &x)); 1182 PetscCall(VecGetArray(yy, &y)); 1183 PetscCallBLAS("BLASgemv", BLASgemv_("T", &m, &n, &_DOne, v, &(mat->lda), x, &_One, &_DOne, y, &_One)); 1184 PetscCall(VecRestoreArrayRead(xx, &x)); 1185 PetscCall(VecRestoreArray(yy, &y)); 1186 PetscCall(PetscLogFlops(2.0 * A->rmap->n * A->cmap->n)); 1187 PetscFunctionReturn(PETSC_SUCCESS); 1188 } 1189 1190 static PetscErrorCode MatGetRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals) 1191 { 1192 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 1193 PetscInt i; 1194 1195 PetscFunctionBegin; 1196 *ncols = A->cmap->n; 1197 if (cols) { 1198 PetscCall(PetscMalloc1(A->cmap->n, cols)); 1199 for (i = 0; i < A->cmap->n; i++) (*cols)[i] = i; 1200 } 1201 if (vals) { 1202 const PetscScalar *v; 1203 1204 PetscCall(MatDenseGetArrayRead(A, &v)); 1205 PetscCall(PetscMalloc1(A->cmap->n, vals)); 1206 v += row; 1207 for (i = 0; i < A->cmap->n; i++) { 1208 (*vals)[i] = *v; 1209 v += mat->lda; 1210 } 1211 PetscCall(MatDenseRestoreArrayRead(A, &v)); 1212 } 1213 PetscFunctionReturn(PETSC_SUCCESS); 1214 } 1215 1216 static PetscErrorCode MatRestoreRow_SeqDense(Mat A, PetscInt row, PetscInt *ncols, PetscInt **cols, PetscScalar **vals) 1217 { 1218 PetscFunctionBegin; 1219 if (ncols) *ncols = 0; 1220 if (cols) PetscCall(PetscFree(*cols)); 1221 if (vals) PetscCall(PetscFree(*vals)); 1222 PetscFunctionReturn(PETSC_SUCCESS); 1223 } 1224 1225 static PetscErrorCode MatSetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], const PetscScalar v[], InsertMode addv) 1226 { 1227 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 1228 PetscScalar *av; 1229 PetscInt i, j, idx = 0; 1230 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP) 1231 PetscOffloadMask oldf; 1232 #endif 1233 1234 PetscFunctionBegin; 1235 PetscCall(MatDenseGetArray(A, &av)); 1236 if (!mat->roworiented) { 1237 if (addv == INSERT_VALUES) { 1238 for (j = 0; j < n; j++) { 1239 if (indexn[j] < 0) { 1240 idx += m; 1241 continue; 1242 } 1243 PetscCheck(indexn[j] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, indexn[j], A->cmap->n - 1); 1244 for (i = 0; i < m; i++) { 1245 if (indexm[i] < 0) { 1246 idx++; 1247 continue; 1248 } 1249 PetscCheck(indexm[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, indexm[i], A->rmap->n - 1); 1250 av[indexn[j] * mat->lda + indexm[i]] = v[idx++]; 1251 } 1252 } 1253 } else { 1254 for (j = 0; j < n; j++) { 1255 if (indexn[j] < 0) { 1256 idx += m; 1257 continue; 1258 } 1259 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); 1260 for (i = 0; i < m; i++) { 1261 if (indexm[i] < 0) { 1262 idx++; 1263 continue; 1264 } 1265 PetscCheck(indexm[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, indexm[i], A->rmap->n - 1); 1266 av[indexn[j] * mat->lda + indexm[i]] += v[idx++]; 1267 } 1268 } 1269 } 1270 } else { 1271 if (addv == INSERT_VALUES) { 1272 for (i = 0; i < m; i++) { 1273 if (indexm[i] < 0) { 1274 idx += n; 1275 continue; 1276 } 1277 PetscCheck(indexm[i] < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, indexm[i], A->rmap->n - 1); 1278 for (j = 0; j < n; j++) { 1279 if (indexn[j] < 0) { 1280 idx++; 1281 continue; 1282 } 1283 PetscCheck(indexn[j] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, indexn[j], A->cmap->n - 1); 1284 av[indexn[j] * mat->lda + indexm[i]] = v[idx++]; 1285 } 1286 } 1287 } else { 1288 for (i = 0; i < m; i++) { 1289 if (indexm[i] < 0) { 1290 idx += n; 1291 continue; 1292 } 1293 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); 1294 for (j = 0; j < n; j++) { 1295 if (indexn[j] < 0) { 1296 idx++; 1297 continue; 1298 } 1299 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); 1300 av[indexn[j] * mat->lda + indexm[i]] += v[idx++]; 1301 } 1302 } 1303 } 1304 } 1305 /* hack to prevent unneeded copy to the GPU while returning the array */ 1306 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP) 1307 oldf = A->offloadmask; 1308 A->offloadmask = PETSC_OFFLOAD_GPU; 1309 #endif 1310 PetscCall(MatDenseRestoreArray(A, &av)); 1311 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP) 1312 A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU); 1313 #endif 1314 PetscFunctionReturn(PETSC_SUCCESS); 1315 } 1316 1317 static PetscErrorCode MatGetValues_SeqDense(Mat A, PetscInt m, const PetscInt indexm[], PetscInt n, const PetscInt indexn[], PetscScalar v[]) 1318 { 1319 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 1320 const PetscScalar *vv; 1321 PetscInt i, j; 1322 1323 PetscFunctionBegin; 1324 PetscCall(MatDenseGetArrayRead(A, &vv)); 1325 /* row-oriented output */ 1326 for (i = 0; i < m; i++) { 1327 if (indexm[i] < 0) { 1328 v += n; 1329 continue; 1330 } 1331 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); 1332 for (j = 0; j < n; j++) { 1333 if (indexn[j] < 0) { 1334 v++; 1335 continue; 1336 } 1337 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); 1338 *v++ = vv[indexn[j] * mat->lda + indexm[i]]; 1339 } 1340 } 1341 PetscCall(MatDenseRestoreArrayRead(A, &vv)); 1342 PetscFunctionReturn(PETSC_SUCCESS); 1343 } 1344 1345 PetscErrorCode MatView_Dense_Binary(Mat mat, PetscViewer viewer) 1346 { 1347 PetscBool skipHeader; 1348 PetscViewerFormat format; 1349 PetscInt header[4], M, N, m, lda, i, j, k; 1350 const PetscScalar *v; 1351 PetscScalar *vwork; 1352 1353 PetscFunctionBegin; 1354 PetscCall(PetscViewerSetUp(viewer)); 1355 PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader)); 1356 PetscCall(PetscViewerGetFormat(viewer, &format)); 1357 if (skipHeader) format = PETSC_VIEWER_NATIVE; 1358 1359 PetscCall(MatGetSize(mat, &M, &N)); 1360 1361 /* write matrix header */ 1362 header[0] = MAT_FILE_CLASSID; 1363 header[1] = M; 1364 header[2] = N; 1365 header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M * N; 1366 if (!skipHeader) PetscCall(PetscViewerBinaryWrite(viewer, header, 4, PETSC_INT)); 1367 1368 PetscCall(MatGetLocalSize(mat, &m, NULL)); 1369 if (format != PETSC_VIEWER_NATIVE) { 1370 PetscInt nnz = m * N, *iwork; 1371 /* store row lengths for each row */ 1372 PetscCall(PetscMalloc1(nnz, &iwork)); 1373 for (i = 0; i < m; i++) iwork[i] = N; 1374 PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT)); 1375 /* store column indices (zero start index) */ 1376 for (k = 0, i = 0; i < m; i++) 1377 for (j = 0; j < N; j++, k++) iwork[k] = j; 1378 PetscCall(PetscViewerBinaryWriteAll(viewer, iwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT)); 1379 PetscCall(PetscFree(iwork)); 1380 } 1381 /* store matrix values as a dense matrix in row major order */ 1382 PetscCall(PetscMalloc1(m * N, &vwork)); 1383 PetscCall(MatDenseGetArrayRead(mat, &v)); 1384 PetscCall(MatDenseGetLDA(mat, &lda)); 1385 for (k = 0, i = 0; i < m; i++) 1386 for (j = 0; j < N; j++, k++) vwork[k] = v[i + lda * j]; 1387 PetscCall(MatDenseRestoreArrayRead(mat, &v)); 1388 PetscCall(PetscViewerBinaryWriteAll(viewer, vwork, m * N, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR)); 1389 PetscCall(PetscFree(vwork)); 1390 PetscFunctionReturn(PETSC_SUCCESS); 1391 } 1392 1393 PetscErrorCode MatLoad_Dense_Binary(Mat mat, PetscViewer viewer) 1394 { 1395 PetscBool skipHeader; 1396 PetscInt header[4], M, N, m, nz, lda, i, j, k; 1397 PetscInt rows, cols; 1398 PetscScalar *v, *vwork; 1399 1400 PetscFunctionBegin; 1401 PetscCall(PetscViewerSetUp(viewer)); 1402 PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &skipHeader)); 1403 1404 if (!skipHeader) { 1405 PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT)); 1406 PetscCheck(header[0] == MAT_FILE_CLASSID, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file"); 1407 M = header[1]; 1408 N = header[2]; 1409 PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M); 1410 PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N); 1411 nz = header[3]; 1412 PetscCheck(nz == MATRIX_BINARY_FORMAT_DENSE || nz >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Unknown matrix format %" PetscInt_FMT " in file", nz); 1413 } else { 1414 PetscCall(MatGetSize(mat, &M, &N)); 1415 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"); 1416 nz = MATRIX_BINARY_FORMAT_DENSE; 1417 } 1418 1419 /* setup global sizes if not set */ 1420 if (mat->rmap->N < 0) mat->rmap->N = M; 1421 if (mat->cmap->N < 0) mat->cmap->N = N; 1422 PetscCall(MatSetUp(mat)); 1423 /* check if global sizes are correct */ 1424 PetscCall(MatGetSize(mat, &rows, &cols)); 1425 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); 1426 1427 PetscCall(MatGetSize(mat, NULL, &N)); 1428 PetscCall(MatGetLocalSize(mat, &m, NULL)); 1429 PetscCall(MatDenseGetArray(mat, &v)); 1430 PetscCall(MatDenseGetLDA(mat, &lda)); 1431 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */ 1432 PetscInt nnz = m * N; 1433 /* read in matrix values */ 1434 PetscCall(PetscMalloc1(nnz, &vwork)); 1435 PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR)); 1436 /* store values in column major order */ 1437 for (j = 0; j < N; j++) 1438 for (i = 0; i < m; i++) v[i + lda * j] = vwork[i * N + j]; 1439 PetscCall(PetscFree(vwork)); 1440 } else { /* matrix in file is sparse format */ 1441 PetscInt nnz = 0, *rlens, *icols; 1442 /* read in row lengths */ 1443 PetscCall(PetscMalloc1(m, &rlens)); 1444 PetscCall(PetscViewerBinaryReadAll(viewer, rlens, m, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT)); 1445 for (i = 0; i < m; i++) nnz += rlens[i]; 1446 /* read in column indices and values */ 1447 PetscCall(PetscMalloc2(nnz, &icols, nnz, &vwork)); 1448 PetscCall(PetscViewerBinaryReadAll(viewer, icols, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_INT)); 1449 PetscCall(PetscViewerBinaryReadAll(viewer, vwork, nnz, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_SCALAR)); 1450 /* store values in column major order */ 1451 for (k = 0, i = 0; i < m; i++) 1452 for (j = 0; j < rlens[i]; j++, k++) v[i + lda * icols[k]] = vwork[k]; 1453 PetscCall(PetscFree(rlens)); 1454 PetscCall(PetscFree2(icols, vwork)); 1455 } 1456 PetscCall(MatDenseRestoreArray(mat, &v)); 1457 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 1458 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 1459 PetscFunctionReturn(PETSC_SUCCESS); 1460 } 1461 1462 PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer) 1463 { 1464 PetscBool isbinary, ishdf5; 1465 1466 PetscFunctionBegin; 1467 PetscValidHeaderSpecific(newMat, MAT_CLASSID, 1); 1468 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 1469 /* force binary viewer to load .info file if it has not yet done so */ 1470 PetscCall(PetscViewerSetUp(viewer)); 1471 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 1472 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5)); 1473 if (isbinary) { 1474 PetscCall(MatLoad_Dense_Binary(newMat, viewer)); 1475 } else if (ishdf5) { 1476 #if defined(PETSC_HAVE_HDF5) 1477 PetscCall(MatLoad_Dense_HDF5(newMat, viewer)); 1478 #else 1479 SETERRQ(PetscObjectComm((PetscObject)newMat), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 1480 #endif 1481 } else { 1482 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); 1483 } 1484 PetscFunctionReturn(PETSC_SUCCESS); 1485 } 1486 1487 static PetscErrorCode MatView_SeqDense_ASCII(Mat A, PetscViewer viewer) 1488 { 1489 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 1490 PetscInt i, j; 1491 const char *name; 1492 PetscScalar *v, *av; 1493 PetscViewerFormat format; 1494 #if defined(PETSC_USE_COMPLEX) 1495 PetscBool allreal = PETSC_TRUE; 1496 #endif 1497 1498 PetscFunctionBegin; 1499 PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&av)); 1500 PetscCall(PetscViewerGetFormat(viewer, &format)); 1501 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 1502 PetscFunctionReturn(PETSC_SUCCESS); /* do nothing for now */ 1503 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 1504 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 1505 for (i = 0; i < A->rmap->n; i++) { 1506 v = av + i; 1507 PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i)); 1508 for (j = 0; j < A->cmap->n; j++) { 1509 #if defined(PETSC_USE_COMPLEX) 1510 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 1511 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", j, (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v))); 1512 } else if (PetscRealPart(*v)) { 1513 PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)PetscRealPart(*v))); 1514 } 1515 #else 1516 if (*v) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", j, (double)*v)); 1517 #endif 1518 v += a->lda; 1519 } 1520 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 1521 } 1522 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 1523 } else { 1524 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE)); 1525 #if defined(PETSC_USE_COMPLEX) 1526 /* determine if matrix has all real values */ 1527 for (j = 0; j < A->cmap->n; j++) { 1528 v = av + j * a->lda; 1529 for (i = 0; i < A->rmap->n; i++) { 1530 if (PetscImaginaryPart(v[i])) { 1531 allreal = PETSC_FALSE; 1532 break; 1533 } 1534 } 1535 } 1536 #endif 1537 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1538 PetscCall(PetscObjectGetName((PetscObject)A, &name)); 1539 PetscCall(PetscViewerASCIIPrintf(viewer, "%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n", A->rmap->n, A->cmap->n)); 1540 PetscCall(PetscViewerASCIIPrintf(viewer, "%s = zeros(%" PetscInt_FMT ",%" PetscInt_FMT ");\n", name, A->rmap->n, A->cmap->n)); 1541 PetscCall(PetscViewerASCIIPrintf(viewer, "%s = [\n", name)); 1542 } 1543 1544 for (i = 0; i < A->rmap->n; i++) { 1545 v = av + i; 1546 for (j = 0; j < A->cmap->n; j++) { 1547 #if defined(PETSC_USE_COMPLEX) 1548 if (allreal) { 1549 PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)PetscRealPart(*v))); 1550 } else { 1551 PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e + %18.16ei ", (double)PetscRealPart(*v), (double)PetscImaginaryPart(*v))); 1552 } 1553 #else 1554 PetscCall(PetscViewerASCIIPrintf(viewer, "%18.16e ", (double)*v)); 1555 #endif 1556 v += a->lda; 1557 } 1558 PetscCall(PetscViewerASCIIPrintf(viewer, "\n")); 1559 } 1560 if (format == PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer, "];\n")); 1561 PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE)); 1562 } 1563 PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&av)); 1564 PetscCall(PetscViewerFlush(viewer)); 1565 PetscFunctionReturn(PETSC_SUCCESS); 1566 } 1567 1568 #include <petscdraw.h> 1569 static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw, void *Aa) 1570 { 1571 Mat A = (Mat)Aa; 1572 PetscInt m = A->rmap->n, n = A->cmap->n, i, j; 1573 int color = PETSC_DRAW_WHITE; 1574 const PetscScalar *v; 1575 PetscViewer viewer; 1576 PetscReal xl, yl, xr, yr, x_l, x_r, y_l, y_r; 1577 PetscViewerFormat format; 1578 1579 PetscFunctionBegin; 1580 PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer)); 1581 PetscCall(PetscViewerGetFormat(viewer, &format)); 1582 PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr)); 1583 1584 /* Loop over matrix elements drawing boxes */ 1585 PetscCall(MatDenseGetArrayRead(A, &v)); 1586 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1587 PetscDrawCollectiveBegin(draw); 1588 /* Blue for negative and Red for positive */ 1589 for (j = 0; j < n; j++) { 1590 x_l = j; 1591 x_r = x_l + 1.0; 1592 for (i = 0; i < m; i++) { 1593 y_l = m - i - 1.0; 1594 y_r = y_l + 1.0; 1595 if (PetscRealPart(v[j * m + i]) > 0.) color = PETSC_DRAW_RED; 1596 else if (PetscRealPart(v[j * m + i]) < 0.) color = PETSC_DRAW_BLUE; 1597 else continue; 1598 PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color)); 1599 } 1600 } 1601 PetscDrawCollectiveEnd(draw); 1602 } else { 1603 /* use contour shading to indicate magnitude of values */ 1604 /* first determine max of all nonzero values */ 1605 PetscReal minv = 0.0, maxv = 0.0; 1606 PetscDraw popup; 1607 1608 for (i = 0; i < m * n; i++) { 1609 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1610 } 1611 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1612 PetscCall(PetscDrawGetPopup(draw, &popup)); 1613 PetscCall(PetscDrawScalePopup(popup, minv, maxv)); 1614 1615 PetscDrawCollectiveBegin(draw); 1616 for (j = 0; j < n; j++) { 1617 x_l = j; 1618 x_r = x_l + 1.0; 1619 for (i = 0; i < m; i++) { 1620 y_l = m - i - 1.0; 1621 y_r = y_l + 1.0; 1622 color = PetscDrawRealToColor(PetscAbsScalar(v[j * m + i]), minv, maxv); 1623 PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color)); 1624 } 1625 } 1626 PetscDrawCollectiveEnd(draw); 1627 } 1628 PetscCall(MatDenseRestoreArrayRead(A, &v)); 1629 PetscFunctionReturn(PETSC_SUCCESS); 1630 } 1631 1632 static PetscErrorCode MatView_SeqDense_Draw(Mat A, PetscViewer viewer) 1633 { 1634 PetscDraw draw; 1635 PetscBool isnull; 1636 PetscReal xr, yr, xl, yl, h, w; 1637 1638 PetscFunctionBegin; 1639 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 1640 PetscCall(PetscDrawIsNull(draw, &isnull)); 1641 if (isnull) PetscFunctionReturn(PETSC_SUCCESS); 1642 1643 xr = A->cmap->n; 1644 yr = A->rmap->n; 1645 h = yr / 10.0; 1646 w = xr / 10.0; 1647 xr += w; 1648 yr += h; 1649 xl = -w; 1650 yl = -h; 1651 PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr)); 1652 PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer)); 1653 PetscCall(PetscDrawZoom(draw, MatView_SeqDense_Draw_Zoom, A)); 1654 PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL)); 1655 PetscCall(PetscDrawSave(draw)); 1656 PetscFunctionReturn(PETSC_SUCCESS); 1657 } 1658 1659 PetscErrorCode MatView_SeqDense(Mat A, PetscViewer viewer) 1660 { 1661 PetscBool iascii, isbinary, isdraw; 1662 1663 PetscFunctionBegin; 1664 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 1665 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 1666 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 1667 if (iascii) PetscCall(MatView_SeqDense_ASCII(A, viewer)); 1668 else if (isbinary) PetscCall(MatView_Dense_Binary(A, viewer)); 1669 else if (isdraw) PetscCall(MatView_SeqDense_Draw(A, viewer)); 1670 PetscFunctionReturn(PETSC_SUCCESS); 1671 } 1672 1673 static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A, const PetscScalar *array) 1674 { 1675 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 1676 1677 PetscFunctionBegin; 1678 PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first"); 1679 PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 1680 PetscCheck(!a->unplacedarray, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreArray() first"); 1681 a->unplacedarray = a->v; 1682 a->unplaced_user_alloc = a->user_alloc; 1683 a->v = (PetscScalar *)array; 1684 a->user_alloc = PETSC_TRUE; 1685 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP) 1686 A->offloadmask = PETSC_OFFLOAD_CPU; 1687 #endif 1688 PetscFunctionReturn(PETSC_SUCCESS); 1689 } 1690 1691 static PetscErrorCode MatDenseResetArray_SeqDense(Mat A) 1692 { 1693 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 1694 1695 PetscFunctionBegin; 1696 PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first"); 1697 PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 1698 a->v = a->unplacedarray; 1699 a->user_alloc = a->unplaced_user_alloc; 1700 a->unplacedarray = NULL; 1701 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP) 1702 A->offloadmask = PETSC_OFFLOAD_CPU; 1703 #endif 1704 PetscFunctionReturn(PETSC_SUCCESS); 1705 } 1706 1707 static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A, const PetscScalar *array) 1708 { 1709 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 1710 1711 PetscFunctionBegin; 1712 PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first"); 1713 PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 1714 if (!a->user_alloc) PetscCall(PetscFree(a->v)); 1715 a->v = (PetscScalar *)array; 1716 a->user_alloc = PETSC_FALSE; 1717 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP) 1718 A->offloadmask = PETSC_OFFLOAD_CPU; 1719 #endif 1720 PetscFunctionReturn(PETSC_SUCCESS); 1721 } 1722 1723 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1724 { 1725 Mat_SeqDense *l = (Mat_SeqDense *)mat->data; 1726 1727 PetscFunctionBegin; 1728 #if defined(PETSC_USE_LOG) 1729 PetscCall(PetscLogObjectState((PetscObject)mat, "Rows %" PetscInt_FMT " Cols %" PetscInt_FMT, mat->rmap->n, mat->cmap->n)); 1730 #endif 1731 PetscCall(VecDestroy(&(l->qrrhs))); 1732 PetscCall(PetscFree(l->tau)); 1733 PetscCall(PetscFree(l->pivots)); 1734 PetscCall(PetscFree(l->fwork)); 1735 PetscCall(MatDestroy(&l->ptapwork)); 1736 if (!l->user_alloc) PetscCall(PetscFree(l->v)); 1737 if (!l->unplaced_user_alloc) PetscCall(PetscFree(l->unplacedarray)); 1738 PetscCheck(!l->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first"); 1739 PetscCheck(!l->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 1740 PetscCall(VecDestroy(&l->cvec)); 1741 PetscCall(MatDestroy(&l->cmat)); 1742 PetscCall(PetscFree(mat->data)); 1743 1744 PetscCall(PetscObjectChangeTypeName((PetscObject)mat, NULL)); 1745 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactor_C", NULL)); 1746 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorSymbolic_C", NULL)); 1747 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatQRFactorNumeric_C", NULL)); 1748 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetLDA_C", NULL)); 1749 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseSetLDA_C", NULL)); 1750 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArray_C", NULL)); 1751 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArray_C", NULL)); 1752 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDensePlaceArray_C", NULL)); 1753 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseResetArray_C", NULL)); 1754 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseReplaceArray_C", NULL)); 1755 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayRead_C", NULL)); 1756 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayRead_C", NULL)); 1757 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetArrayWrite_C", NULL)); 1758 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreArrayWrite_C", NULL)); 1759 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqaij_C", NULL)); 1760 #if defined(PETSC_HAVE_ELEMENTAL) 1761 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_elemental_C", NULL)); 1762 #endif 1763 #if defined(PETSC_HAVE_SCALAPACK) 1764 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_scalapack_C", NULL)); 1765 #endif 1766 #if defined(PETSC_HAVE_CUDA) 1767 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensecuda_C", NULL)); 1768 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", NULL)); 1769 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensecuda_seqdense_C", NULL)); 1770 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensecuda_C", NULL)); 1771 #endif 1772 #if defined(PETSC_HAVE_HIP) 1773 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatConvert_seqdense_seqdensehip_C", NULL)); 1774 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", NULL)); 1775 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdensehip_seqdense_C", NULL)); 1776 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdensehip_C", NULL)); 1777 #endif 1778 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatSeqDenseSetPreallocation_C", NULL)); 1779 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqaij_seqdense_C", NULL)); 1780 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqdense_seqdense_C", NULL)); 1781 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqbaij_seqdense_C", NULL)); 1782 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatProductSetFromOptions_seqsbaij_seqdense_C", NULL)); 1783 1784 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumn_C", NULL)); 1785 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumn_C", NULL)); 1786 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVec_C", NULL)); 1787 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVec_C", NULL)); 1788 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecRead_C", NULL)); 1789 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecRead_C", NULL)); 1790 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetColumnVecWrite_C", NULL)); 1791 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreColumnVecWrite_C", NULL)); 1792 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseGetSubMatrix_C", NULL)); 1793 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatDenseRestoreSubMatrix_C", NULL)); 1794 PetscFunctionReturn(PETSC_SUCCESS); 1795 } 1796 1797 static PetscErrorCode MatTranspose_SeqDense(Mat A, MatReuse reuse, Mat *matout) 1798 { 1799 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 1800 PetscInt k, j, m = A->rmap->n, M = mat->lda, n = A->cmap->n; 1801 PetscScalar *v, tmp; 1802 1803 PetscFunctionBegin; 1804 if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *matout)); 1805 if (reuse == MAT_INPLACE_MATRIX) { 1806 if (m == n) { /* in place transpose */ 1807 PetscCall(MatDenseGetArray(A, &v)); 1808 for (j = 0; j < m; j++) { 1809 for (k = 0; k < j; k++) { 1810 tmp = v[j + k * M]; 1811 v[j + k * M] = v[k + j * M]; 1812 v[k + j * M] = tmp; 1813 } 1814 } 1815 PetscCall(MatDenseRestoreArray(A, &v)); 1816 } else { /* reuse memory, temporary allocates new memory */ 1817 PetscScalar *v2; 1818 PetscLayout tmplayout; 1819 1820 PetscCall(PetscMalloc1((size_t)m * n, &v2)); 1821 PetscCall(MatDenseGetArray(A, &v)); 1822 for (j = 0; j < n; j++) { 1823 for (k = 0; k < m; k++) v2[j + (size_t)k * n] = v[k + (size_t)j * M]; 1824 } 1825 PetscCall(PetscArraycpy(v, v2, (size_t)m * n)); 1826 PetscCall(PetscFree(v2)); 1827 PetscCall(MatDenseRestoreArray(A, &v)); 1828 /* cleanup size dependent quantities */ 1829 PetscCall(VecDestroy(&mat->cvec)); 1830 PetscCall(MatDestroy(&mat->cmat)); 1831 PetscCall(PetscFree(mat->pivots)); 1832 PetscCall(PetscFree(mat->fwork)); 1833 PetscCall(MatDestroy(&mat->ptapwork)); 1834 /* swap row/col layouts */ 1835 mat->lda = n; 1836 tmplayout = A->rmap; 1837 A->rmap = A->cmap; 1838 A->cmap = tmplayout; 1839 } 1840 } else { /* out-of-place transpose */ 1841 Mat tmat; 1842 Mat_SeqDense *tmatd; 1843 PetscScalar *v2; 1844 PetscInt M2; 1845 1846 if (reuse == MAT_INITIAL_MATRIX) { 1847 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &tmat)); 1848 PetscCall(MatSetSizes(tmat, A->cmap->n, A->rmap->n, A->cmap->n, A->rmap->n)); 1849 PetscCall(MatSetType(tmat, ((PetscObject)A)->type_name)); 1850 PetscCall(MatSeqDenseSetPreallocation(tmat, NULL)); 1851 } else tmat = *matout; 1852 1853 PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&v)); 1854 PetscCall(MatDenseGetArray(tmat, &v2)); 1855 tmatd = (Mat_SeqDense *)tmat->data; 1856 M2 = tmatd->lda; 1857 for (j = 0; j < n; j++) { 1858 for (k = 0; k < m; k++) v2[j + k * M2] = v[k + j * M]; 1859 } 1860 PetscCall(MatDenseRestoreArray(tmat, &v2)); 1861 PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&v)); 1862 PetscCall(MatAssemblyBegin(tmat, MAT_FINAL_ASSEMBLY)); 1863 PetscCall(MatAssemblyEnd(tmat, MAT_FINAL_ASSEMBLY)); 1864 *matout = tmat; 1865 } 1866 PetscFunctionReturn(PETSC_SUCCESS); 1867 } 1868 1869 static PetscErrorCode MatEqual_SeqDense(Mat A1, Mat A2, PetscBool *flg) 1870 { 1871 Mat_SeqDense *mat1 = (Mat_SeqDense *)A1->data; 1872 Mat_SeqDense *mat2 = (Mat_SeqDense *)A2->data; 1873 PetscInt i; 1874 const PetscScalar *v1, *v2; 1875 1876 PetscFunctionBegin; 1877 if (A1->rmap->n != A2->rmap->n) { 1878 *flg = PETSC_FALSE; 1879 PetscFunctionReturn(PETSC_SUCCESS); 1880 } 1881 if (A1->cmap->n != A2->cmap->n) { 1882 *flg = PETSC_FALSE; 1883 PetscFunctionReturn(PETSC_SUCCESS); 1884 } 1885 PetscCall(MatDenseGetArrayRead(A1, &v1)); 1886 PetscCall(MatDenseGetArrayRead(A2, &v2)); 1887 for (i = 0; i < A1->cmap->n; i++) { 1888 PetscCall(PetscArraycmp(v1, v2, A1->rmap->n, flg)); 1889 if (*flg == PETSC_FALSE) PetscFunctionReturn(PETSC_SUCCESS); 1890 v1 += mat1->lda; 1891 v2 += mat2->lda; 1892 } 1893 PetscCall(MatDenseRestoreArrayRead(A1, &v1)); 1894 PetscCall(MatDenseRestoreArrayRead(A2, &v2)); 1895 *flg = PETSC_TRUE; 1896 PetscFunctionReturn(PETSC_SUCCESS); 1897 } 1898 1899 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A, Vec v) 1900 { 1901 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 1902 PetscInt i, n, len; 1903 PetscScalar *x; 1904 const PetscScalar *vv; 1905 1906 PetscFunctionBegin; 1907 PetscCall(VecGetSize(v, &n)); 1908 PetscCall(VecGetArray(v, &x)); 1909 len = PetscMin(A->rmap->n, A->cmap->n); 1910 PetscCall(MatDenseGetArrayRead(A, &vv)); 1911 PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming mat and vec"); 1912 for (i = 0; i < len; i++) x[i] = vv[i * mat->lda + i]; 1913 PetscCall(MatDenseRestoreArrayRead(A, &vv)); 1914 PetscCall(VecRestoreArray(v, &x)); 1915 PetscFunctionReturn(PETSC_SUCCESS); 1916 } 1917 1918 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A, Vec ll, Vec rr) 1919 { 1920 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 1921 const PetscScalar *l, *r; 1922 PetscScalar x, *v, *vv; 1923 PetscInt i, j, m = A->rmap->n, n = A->cmap->n; 1924 1925 PetscFunctionBegin; 1926 PetscCall(MatDenseGetArray(A, &vv)); 1927 if (ll) { 1928 PetscCall(VecGetSize(ll, &m)); 1929 PetscCall(VecGetArrayRead(ll, &l)); 1930 PetscCheck(m == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vec wrong size"); 1931 for (i = 0; i < m; i++) { 1932 x = l[i]; 1933 v = vv + i; 1934 for (j = 0; j < n; j++) { 1935 (*v) *= x; 1936 v += mat->lda; 1937 } 1938 } 1939 PetscCall(VecRestoreArrayRead(ll, &l)); 1940 PetscCall(PetscLogFlops(1.0 * n * m)); 1941 } 1942 if (rr) { 1943 PetscCall(VecGetSize(rr, &n)); 1944 PetscCall(VecGetArrayRead(rr, &r)); 1945 PetscCheck(n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vec wrong size"); 1946 for (i = 0; i < n; i++) { 1947 x = r[i]; 1948 v = vv + i * mat->lda; 1949 for (j = 0; j < m; j++) (*v++) *= x; 1950 } 1951 PetscCall(VecRestoreArrayRead(rr, &r)); 1952 PetscCall(PetscLogFlops(1.0 * n * m)); 1953 } 1954 PetscCall(MatDenseRestoreArray(A, &vv)); 1955 PetscFunctionReturn(PETSC_SUCCESS); 1956 } 1957 1958 PetscErrorCode MatNorm_SeqDense(Mat A, NormType type, PetscReal *nrm) 1959 { 1960 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 1961 PetscScalar *v, *vv; 1962 PetscReal sum = 0.0; 1963 PetscInt lda, m = A->rmap->n, i, j; 1964 1965 PetscFunctionBegin; 1966 PetscCall(MatDenseGetArrayRead(A, (const PetscScalar **)&vv)); 1967 PetscCall(MatDenseGetLDA(A, &lda)); 1968 v = vv; 1969 if (type == NORM_FROBENIUS) { 1970 if (lda > m) { 1971 for (j = 0; j < A->cmap->n; j++) { 1972 v = vv + j * lda; 1973 for (i = 0; i < m; i++) { 1974 sum += PetscRealPart(PetscConj(*v) * (*v)); 1975 v++; 1976 } 1977 } 1978 } else { 1979 #if defined(PETSC_USE_REAL___FP16) 1980 PetscBLASInt one = 1, cnt = A->cmap->n * A->rmap->n; 1981 PetscCallBLAS("BLASnrm2", *nrm = BLASnrm2_(&cnt, v, &one)); 1982 } 1983 #else 1984 for (i = 0; i < A->cmap->n * A->rmap->n; i++) { 1985 sum += PetscRealPart(PetscConj(*v) * (*v)); 1986 v++; 1987 } 1988 } 1989 *nrm = PetscSqrtReal(sum); 1990 #endif 1991 PetscCall(PetscLogFlops(2.0 * A->cmap->n * A->rmap->n)); 1992 } else if (type == NORM_1) { 1993 *nrm = 0.0; 1994 for (j = 0; j < A->cmap->n; j++) { 1995 v = vv + j * mat->lda; 1996 sum = 0.0; 1997 for (i = 0; i < A->rmap->n; i++) { 1998 sum += PetscAbsScalar(*v); 1999 v++; 2000 } 2001 if (sum > *nrm) *nrm = sum; 2002 } 2003 PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n)); 2004 } else if (type == NORM_INFINITY) { 2005 *nrm = 0.0; 2006 for (j = 0; j < A->rmap->n; j++) { 2007 v = vv + j; 2008 sum = 0.0; 2009 for (i = 0; i < A->cmap->n; i++) { 2010 sum += PetscAbsScalar(*v); 2011 v += mat->lda; 2012 } 2013 if (sum > *nrm) *nrm = sum; 2014 } 2015 PetscCall(PetscLogFlops(1.0 * A->cmap->n * A->rmap->n)); 2016 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No two norm"); 2017 PetscCall(MatDenseRestoreArrayRead(A, (const PetscScalar **)&vv)); 2018 PetscFunctionReturn(PETSC_SUCCESS); 2019 } 2020 2021 static PetscErrorCode MatSetOption_SeqDense(Mat A, MatOption op, PetscBool flg) 2022 { 2023 Mat_SeqDense *aij = (Mat_SeqDense *)A->data; 2024 2025 PetscFunctionBegin; 2026 switch (op) { 2027 case MAT_ROW_ORIENTED: 2028 aij->roworiented = flg; 2029 break; 2030 case MAT_NEW_NONZERO_LOCATIONS: 2031 case MAT_NEW_NONZERO_LOCATION_ERR: 2032 case MAT_NEW_NONZERO_ALLOCATION_ERR: 2033 case MAT_FORCE_DIAGONAL_ENTRIES: 2034 case MAT_KEEP_NONZERO_PATTERN: 2035 case MAT_IGNORE_OFF_PROC_ENTRIES: 2036 case MAT_USE_HASH_TABLE: 2037 case MAT_IGNORE_ZERO_ENTRIES: 2038 case MAT_IGNORE_LOWER_TRIANGULAR: 2039 case MAT_SORTED_FULL: 2040 PetscCall(PetscInfo(A, "Option %s ignored\n", MatOptions[op])); 2041 break; 2042 case MAT_SPD: 2043 case MAT_SYMMETRIC: 2044 case MAT_STRUCTURALLY_SYMMETRIC: 2045 case MAT_HERMITIAN: 2046 case MAT_SYMMETRY_ETERNAL: 2047 case MAT_STRUCTURAL_SYMMETRY_ETERNAL: 2048 case MAT_SPD_ETERNAL: 2049 break; 2050 default: 2051 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unknown option %s", MatOptions[op]); 2052 } 2053 PetscFunctionReturn(PETSC_SUCCESS); 2054 } 2055 2056 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 2057 { 2058 Mat_SeqDense *l = (Mat_SeqDense *)A->data; 2059 PetscInt lda = l->lda, m = A->rmap->n, n = A->cmap->n, j; 2060 PetscScalar *v; 2061 2062 PetscFunctionBegin; 2063 PetscCall(MatDenseGetArrayWrite(A, &v)); 2064 if (lda > m) { 2065 for (j = 0; j < n; j++) PetscCall(PetscArrayzero(v + j * lda, m)); 2066 } else { 2067 PetscCall(PetscArrayzero(v, PetscInt64Mult(m, n))); 2068 } 2069 PetscCall(MatDenseRestoreArrayWrite(A, &v)); 2070 PetscFunctionReturn(PETSC_SUCCESS); 2071 } 2072 2073 static PetscErrorCode MatZeroRows_SeqDense(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b) 2074 { 2075 Mat_SeqDense *l = (Mat_SeqDense *)A->data; 2076 PetscInt m = l->lda, n = A->cmap->n, i, j; 2077 PetscScalar *slot, *bb, *v; 2078 const PetscScalar *xx; 2079 2080 PetscFunctionBegin; 2081 if (PetscDefined(USE_DEBUG)) { 2082 for (i = 0; i < N; i++) { 2083 PetscCheck(rows[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row requested to be zeroed"); 2084 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); 2085 } 2086 } 2087 if (!N) PetscFunctionReturn(PETSC_SUCCESS); 2088 2089 /* fix right hand side if needed */ 2090 if (x && b) { 2091 PetscCall(VecGetArrayRead(x, &xx)); 2092 PetscCall(VecGetArray(b, &bb)); 2093 for (i = 0; i < N; i++) bb[rows[i]] = diag * xx[rows[i]]; 2094 PetscCall(VecRestoreArrayRead(x, &xx)); 2095 PetscCall(VecRestoreArray(b, &bb)); 2096 } 2097 2098 PetscCall(MatDenseGetArray(A, &v)); 2099 for (i = 0; i < N; i++) { 2100 slot = v + rows[i]; 2101 for (j = 0; j < n; j++) { 2102 *slot = 0.0; 2103 slot += m; 2104 } 2105 } 2106 if (diag != 0.0) { 2107 PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only coded for square matrices"); 2108 for (i = 0; i < N; i++) { 2109 slot = v + (m + 1) * rows[i]; 2110 *slot = diag; 2111 } 2112 } 2113 PetscCall(MatDenseRestoreArray(A, &v)); 2114 PetscFunctionReturn(PETSC_SUCCESS); 2115 } 2116 2117 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A, PetscInt *lda) 2118 { 2119 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 2120 2121 PetscFunctionBegin; 2122 *lda = mat->lda; 2123 PetscFunctionReturn(PETSC_SUCCESS); 2124 } 2125 2126 PetscErrorCode MatDenseGetArray_SeqDense(Mat A, PetscScalar **array) 2127 { 2128 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 2129 2130 PetscFunctionBegin; 2131 PetscCheck(!mat->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 2132 *array = mat->v; 2133 PetscFunctionReturn(PETSC_SUCCESS); 2134 } 2135 2136 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A, PetscScalar **array) 2137 { 2138 PetscFunctionBegin; 2139 if (array) *array = NULL; 2140 PetscFunctionReturn(PETSC_SUCCESS); 2141 } 2142 2143 /*@ 2144 MatDenseGetLDA - gets the leading dimension of the array returned from `MatDenseGetArray()` 2145 2146 Not Collective 2147 2148 Input Parameter: 2149 . mat - a `MATDENSE` or `MATDENSECUDA` matrix 2150 2151 Output Parameter: 2152 . lda - the leading dimension 2153 2154 Level: intermediate 2155 2156 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseSetLDA()` 2157 @*/ 2158 PetscErrorCode MatDenseGetLDA(Mat A, PetscInt *lda) 2159 { 2160 PetscFunctionBegin; 2161 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2162 PetscValidIntPointer(lda, 2); 2163 MatCheckPreallocated(A, 1); 2164 PetscUseMethod(A, "MatDenseGetLDA_C", (Mat, PetscInt *), (A, lda)); 2165 PetscFunctionReturn(PETSC_SUCCESS); 2166 } 2167 2168 /*@ 2169 MatDenseSetLDA - Sets the leading dimension of the array used by the `MATDENSE` matrix 2170 2171 Not Collective 2172 2173 Input Parameters: 2174 + mat - a `MATDENSE` or `MATDENSECUDA` matrix 2175 - lda - the leading dimension 2176 2177 Level: intermediate 2178 2179 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetLDA()` 2180 @*/ 2181 PetscErrorCode MatDenseSetLDA(Mat A, PetscInt lda) 2182 { 2183 PetscFunctionBegin; 2184 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2185 PetscTryMethod(A, "MatDenseSetLDA_C", (Mat, PetscInt), (A, lda)); 2186 PetscFunctionReturn(PETSC_SUCCESS); 2187 } 2188 2189 /*@C 2190 MatDenseGetArray - gives read-write access to the array where the data for a `MATDENSE` matrix is stored 2191 2192 Logically Collective 2193 2194 Input Parameter: 2195 . mat - a dense matrix 2196 2197 Output Parameter: 2198 . array - pointer to the data 2199 2200 Level: intermediate 2201 2202 Fortran Note: 2203 `MatDenseGetArray()` Fortran binding is deprecated (since PETSc 3.19), use `MatDenseGetArrayF90()` 2204 2205 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2206 @*/ 2207 PetscErrorCode MatDenseGetArray(Mat A, PetscScalar **array) 2208 { 2209 PetscFunctionBegin; 2210 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2211 PetscValidPointer(array, 2); 2212 PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array)); 2213 PetscFunctionReturn(PETSC_SUCCESS); 2214 } 2215 2216 /*@C 2217 MatDenseRestoreArray - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArray()` 2218 2219 Logically Collective 2220 2221 Input Parameters: 2222 + mat - a dense matrix 2223 - array - pointer to the data (may be `NULL`) 2224 2225 Level: intermediate 2226 2227 Fortran Note: 2228 `MatDenseRestoreArray()` Fortran binding is deprecated (since PETSc 3.19), use `MatDenseRestoreArrayF90()` 2229 2230 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2231 @*/ 2232 PetscErrorCode MatDenseRestoreArray(Mat A, PetscScalar **array) 2233 { 2234 PetscFunctionBegin; 2235 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2236 PetscValidPointer(array, 2); 2237 PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array)); 2238 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 2239 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP) 2240 A->offloadmask = PETSC_OFFLOAD_CPU; 2241 #endif 2242 PetscFunctionReturn(PETSC_SUCCESS); 2243 } 2244 2245 /*@C 2246 MatDenseGetArrayRead - gives read-only access to the array where the data for a `MATDENSE` matrix is stored 2247 2248 Not Collective; No Fortran Support 2249 2250 Input Parameter: 2251 . mat - a dense matrix 2252 2253 Output Parameter: 2254 . array - pointer to the data 2255 2256 Level: intermediate 2257 2258 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2259 @*/ 2260 PetscErrorCode MatDenseGetArrayRead(Mat A, const PetscScalar **array) 2261 { 2262 PetscFunctionBegin; 2263 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2264 PetscValidPointer(array, 2); 2265 PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, const PetscScalar **), (A, array)); 2266 PetscFunctionReturn(PETSC_SUCCESS); 2267 } 2268 2269 /*@C 2270 MatDenseRestoreArrayRead - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayRead()` 2271 2272 Not Collective; No Fortran Support 2273 2274 Input Parameters: 2275 + mat - a dense matrix 2276 - array - pointer to the data (may be `NULL`) 2277 2278 Level: intermediate 2279 2280 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayRead()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2281 @*/ 2282 PetscErrorCode MatDenseRestoreArrayRead(Mat A, const PetscScalar **array) 2283 { 2284 PetscFunctionBegin; 2285 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2286 PetscValidPointer(array, 2); 2287 PetscUseMethod(A, "MatDenseRestoreArrayRead_C", (Mat, const PetscScalar **), (A, array)); 2288 PetscFunctionReturn(PETSC_SUCCESS); 2289 } 2290 2291 /*@C 2292 MatDenseGetArrayWrite - gives write-only access to the array where the data for a `MATDENSE` matrix is stored 2293 2294 Not Collective; No Fortran Support 2295 2296 Input Parameter: 2297 . mat - a dense matrix 2298 2299 Output Parameter: 2300 . array - pointer to the data 2301 2302 Level: intermediate 2303 2304 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()` 2305 @*/ 2306 PetscErrorCode MatDenseGetArrayWrite(Mat A, PetscScalar **array) 2307 { 2308 PetscFunctionBegin; 2309 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2310 PetscValidPointer(array, 2); 2311 PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array)); 2312 PetscFunctionReturn(PETSC_SUCCESS); 2313 } 2314 2315 /*@C 2316 MatDenseRestoreArrayWrite - returns access to the array where the data for a `MATDENSE` matrix is stored obtained by `MatDenseGetArrayWrite()` 2317 2318 Not Collective; No Fortran Support 2319 2320 Input Parameters: 2321 + mat - a dense matrix 2322 - array - pointer to the data (may be `NULL`) 2323 2324 Level: intermediate 2325 2326 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWrite()`, `MatDenseGetArray()`, `MatDenseRestoreArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()` 2327 @*/ 2328 PetscErrorCode MatDenseRestoreArrayWrite(Mat A, PetscScalar **array) 2329 { 2330 PetscFunctionBegin; 2331 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2332 PetscValidPointer(array, 2); 2333 PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array)); 2334 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 2335 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP) 2336 A->offloadmask = PETSC_OFFLOAD_CPU; 2337 #endif 2338 PetscFunctionReturn(PETSC_SUCCESS); 2339 } 2340 2341 /*@C 2342 MatDenseGetArrayAndMemType - gives read-write access to the array where the data for a `MATDENSE` matrix is stored 2343 2344 Logically Collective 2345 2346 Input Parameter: 2347 . mat - a dense matrix 2348 2349 Output Parameters: 2350 + array - pointer to the data 2351 - mtype - memory type of the returned pointer 2352 2353 Level: intermediate 2354 2355 Notes: 2356 If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc., 2357 an array on device is always returned and is guaranteed to contain the matrix's latest data. 2358 2359 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayAndMemType()`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArrayRead()`, 2360 `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()` 2361 @*/ 2362 PetscErrorCode MatDenseGetArrayAndMemType(Mat A, PetscScalar **array, PetscMemType *mtype) 2363 { 2364 PetscBool isMPI; 2365 2366 PetscFunctionBegin; 2367 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2368 PetscValidPointer(array, 2); 2369 PetscCall(MatBindToCPU(A, PETSC_FALSE)); /* We want device matrices to always return device arrays, so we unbind the matrix if it is bound to CPU */ 2370 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI)); 2371 if (isMPI) { 2372 /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */ 2373 PetscCall(MatDenseGetArrayAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype)); 2374 } else { 2375 PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *); 2376 2377 PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayAndMemType_C", &fptr)); 2378 if (fptr) { 2379 PetscCall((*fptr)(A, array, mtype)); 2380 } else { 2381 PetscUseMethod(A, "MatDenseGetArray_C", (Mat, PetscScalar **), (A, array)); 2382 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 2383 } 2384 } 2385 PetscFunctionReturn(PETSC_SUCCESS); 2386 } 2387 2388 /*@C 2389 MatDenseRestoreArrayAndMemType - returns access to the array that is obtained by `MatDenseGetArrayAndMemType()` 2390 2391 Logically Collective 2392 2393 Input Parameters: 2394 + mat - a dense matrix 2395 - array - pointer to the data 2396 2397 Level: intermediate 2398 2399 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2400 @*/ 2401 PetscErrorCode MatDenseRestoreArrayAndMemType(Mat A, PetscScalar **array) 2402 { 2403 PetscBool isMPI; 2404 2405 PetscFunctionBegin; 2406 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2407 PetscValidPointer(array, 2); 2408 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI)); 2409 if (isMPI) { 2410 PetscCall(MatDenseRestoreArrayAndMemType(((Mat_MPIDense *)A->data)->A, array)); 2411 } else { 2412 PetscErrorCode (*fptr)(Mat, PetscScalar **); 2413 2414 PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayAndMemType_C", &fptr)); 2415 if (fptr) { 2416 PetscCall((*fptr)(A, array)); 2417 } else { 2418 PetscUseMethod(A, "MatDenseRestoreArray_C", (Mat, PetscScalar **), (A, array)); 2419 } 2420 *array = NULL; 2421 } 2422 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 2423 PetscFunctionReturn(PETSC_SUCCESS); 2424 } 2425 2426 /*@C 2427 MatDenseGetArrayReadAndMemType - gives read-only access to the array where the data for a `MATDENSE` matrix is stored 2428 2429 Logically Collective 2430 2431 Input Parameter: 2432 . mat - a dense matrix 2433 2434 Output Parameters: 2435 + array - pointer to the data 2436 - mtype - memory type of the returned pointer 2437 2438 Level: intermediate 2439 2440 Notes: 2441 If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc., 2442 an array on device is always returned and is guaranteed to contain the matrix's latest data. 2443 2444 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreArrayReadAndMemType()`, `MatDenseGetArrayWriteAndMemType()`, 2445 `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()`, `MatSeqAIJGetCSRAndMemType()` 2446 @*/ 2447 PetscErrorCode MatDenseGetArrayReadAndMemType(Mat A, const PetscScalar **array, PetscMemType *mtype) 2448 { 2449 PetscBool isMPI; 2450 2451 PetscFunctionBegin; 2452 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2453 PetscValidPointer(array, 2); 2454 PetscCall(MatBindToCPU(A, PETSC_FALSE)); /* We want device matrices to always return device arrays, so we unbind the matrix if it is bound to CPU */ 2455 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI)); 2456 if (isMPI) { /* Dispatch here so that the code can be reused for all subclasses of MATDENSE */ 2457 PetscCall(MatDenseGetArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype)); 2458 } else { 2459 PetscErrorCode (*fptr)(Mat, const PetscScalar **, PetscMemType *); 2460 2461 PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayReadAndMemType_C", &fptr)); 2462 if (fptr) { 2463 PetscCall((*fptr)(A, array, mtype)); 2464 } else { 2465 PetscUseMethod(A, "MatDenseGetArrayRead_C", (Mat, const PetscScalar **), (A, array)); 2466 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 2467 } 2468 } 2469 PetscFunctionReturn(PETSC_SUCCESS); 2470 } 2471 2472 /*@C 2473 MatDenseRestoreArrayReadAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()` 2474 2475 Logically Collective 2476 2477 Input Parameters: 2478 + mat - a dense matrix 2479 - array - pointer to the data 2480 2481 Level: intermediate 2482 2483 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayReadAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2484 @*/ 2485 PetscErrorCode MatDenseRestoreArrayReadAndMemType(Mat A, const PetscScalar **array) 2486 { 2487 PetscBool isMPI; 2488 2489 PetscFunctionBegin; 2490 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2491 PetscValidPointer(array, 2); 2492 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI)); 2493 if (isMPI) { 2494 PetscCall(MatDenseRestoreArrayReadAndMemType(((Mat_MPIDense *)A->data)->A, array)); 2495 } else { 2496 PetscErrorCode (*fptr)(Mat, const PetscScalar **); 2497 2498 PetscCall(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(PETSC_SUCCESS); 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 Level: intermediate 2522 2523 Notes: 2524 If the matrix is of a device type such as `MATDENSECUDA`, `MATDENSEHIP`, etc., 2525 an array on device is always returned and is guaranteed to contain the matrix's latest data. 2526 2527 .seealso: [](chapter_matrices), `Mat`, `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(MatBindToCPU(A, PETSC_FALSE)); /* We want device matrices to always return device arrays, so we unbind the matrix if it is bound to CPU */ 2538 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI)); 2539 if (isMPI) { 2540 PetscCall(MatDenseGetArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array, mtype)); 2541 } else { 2542 PetscErrorCode (*fptr)(Mat, PetscScalar **, PetscMemType *); 2543 2544 PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseGetArrayWriteAndMemType_C", &fptr)); 2545 if (fptr) { 2546 PetscCall((*fptr)(A, array, mtype)); 2547 } else { 2548 PetscUseMethod(A, "MatDenseGetArrayWrite_C", (Mat, PetscScalar **), (A, array)); 2549 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 2550 } 2551 } 2552 PetscFunctionReturn(PETSC_SUCCESS); 2553 } 2554 2555 /*@C 2556 MatDenseRestoreArrayWriteAndMemType - returns access to the array that is obtained by `MatDenseGetArrayReadAndMemType()` 2557 2558 Logically Collective 2559 2560 Input Parameters: 2561 + mat - a dense matrix 2562 - array - pointer to the data 2563 2564 Level: intermediate 2565 2566 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MatDenseGetArrayWriteAndMemType()`, `MatDenseGetArray()`, `MatDenseGetArrayRead()`, `MatDenseRestoreArrayRead()`, `MatDenseGetArrayWrite()`, `MatDenseRestoreArrayWrite()` 2567 @*/ 2568 PetscErrorCode MatDenseRestoreArrayWriteAndMemType(Mat A, PetscScalar **array) 2569 { 2570 PetscBool isMPI; 2571 2572 PetscFunctionBegin; 2573 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 2574 PetscValidPointer(array, 2); 2575 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIDENSE, &isMPI)); 2576 if (isMPI) { 2577 PetscCall(MatDenseRestoreArrayWriteAndMemType(((Mat_MPIDense *)A->data)->A, array)); 2578 } else { 2579 PetscErrorCode (*fptr)(Mat, PetscScalar **); 2580 2581 PetscCall(PetscObjectQueryFunction((PetscObject)A, "MatDenseRestoreArrayWriteAndMemType_C", &fptr)); 2582 if (fptr) { 2583 PetscCall((*fptr)(A, array)); 2584 } else { 2585 PetscUseMethod(A, "MatDenseRestoreArrayWrite_C", (Mat, PetscScalar **), (A, array)); 2586 } 2587 *array = NULL; 2588 } 2589 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 2590 PetscFunctionReturn(PETSC_SUCCESS); 2591 } 2592 2593 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B) 2594 { 2595 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 2596 PetscInt i, j, nrows, ncols, ldb; 2597 const PetscInt *irow, *icol; 2598 PetscScalar *av, *bv, *v = mat->v; 2599 Mat newmat; 2600 2601 PetscFunctionBegin; 2602 PetscCall(ISGetIndices(isrow, &irow)); 2603 PetscCall(ISGetIndices(iscol, &icol)); 2604 PetscCall(ISGetLocalSize(isrow, &nrows)); 2605 PetscCall(ISGetLocalSize(iscol, &ncols)); 2606 2607 /* Check submatrixcall */ 2608 if (scall == MAT_REUSE_MATRIX) { 2609 PetscInt n_cols, n_rows; 2610 PetscCall(MatGetSize(*B, &n_rows, &n_cols)); 2611 if (n_rows != nrows || n_cols != ncols) { 2612 /* resize the result matrix to match number of requested rows/columns */ 2613 PetscCall(MatSetSizes(*B, nrows, ncols, nrows, ncols)); 2614 } 2615 newmat = *B; 2616 } else { 2617 /* Create and fill new matrix */ 2618 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &newmat)); 2619 PetscCall(MatSetSizes(newmat, nrows, ncols, nrows, ncols)); 2620 PetscCall(MatSetType(newmat, ((PetscObject)A)->type_name)); 2621 PetscCall(MatSeqDenseSetPreallocation(newmat, NULL)); 2622 } 2623 2624 /* Now extract the data pointers and do the copy,column at a time */ 2625 PetscCall(MatDenseGetArray(newmat, &bv)); 2626 PetscCall(MatDenseGetLDA(newmat, &ldb)); 2627 for (i = 0; i < ncols; i++) { 2628 av = v + mat->lda * icol[i]; 2629 for (j = 0; j < nrows; j++) bv[j] = av[irow[j]]; 2630 bv += ldb; 2631 } 2632 PetscCall(MatDenseRestoreArray(newmat, &bv)); 2633 2634 /* Assemble the matrices so that the correct flags are set */ 2635 PetscCall(MatAssemblyBegin(newmat, MAT_FINAL_ASSEMBLY)); 2636 PetscCall(MatAssemblyEnd(newmat, MAT_FINAL_ASSEMBLY)); 2637 2638 /* Free work space */ 2639 PetscCall(ISRestoreIndices(isrow, &irow)); 2640 PetscCall(ISRestoreIndices(iscol, &icol)); 2641 *B = newmat; 2642 PetscFunctionReturn(PETSC_SUCCESS); 2643 } 2644 2645 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[]) 2646 { 2647 PetscInt i; 2648 2649 PetscFunctionBegin; 2650 if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n, B)); 2651 2652 for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SeqDense(A, irow[i], icol[i], scall, &(*B)[i])); 2653 PetscFunctionReturn(PETSC_SUCCESS); 2654 } 2655 2656 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat, MatAssemblyType mode) 2657 { 2658 PetscFunctionBegin; 2659 PetscFunctionReturn(PETSC_SUCCESS); 2660 } 2661 2662 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat, MatAssemblyType mode) 2663 { 2664 PetscFunctionBegin; 2665 PetscFunctionReturn(PETSC_SUCCESS); 2666 } 2667 2668 PetscErrorCode MatCopy_SeqDense(Mat A, Mat B, MatStructure str) 2669 { 2670 Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data; 2671 const PetscScalar *va; 2672 PetscScalar *vb; 2673 PetscInt lda1 = a->lda, lda2 = b->lda, m = A->rmap->n, n = A->cmap->n, j; 2674 2675 PetscFunctionBegin; 2676 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2677 if (A->ops->copy != B->ops->copy) { 2678 PetscCall(MatCopy_Basic(A, B, str)); 2679 PetscFunctionReturn(PETSC_SUCCESS); 2680 } 2681 PetscCheck(m == B->rmap->n && n == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "size(B) != size(A)"); 2682 PetscCall(MatDenseGetArrayRead(A, &va)); 2683 PetscCall(MatDenseGetArray(B, &vb)); 2684 if (lda1 > m || lda2 > m) { 2685 for (j = 0; j < n; j++) PetscCall(PetscArraycpy(vb + j * lda2, va + j * lda1, m)); 2686 } else { 2687 PetscCall(PetscArraycpy(vb, va, A->rmap->n * A->cmap->n)); 2688 } 2689 PetscCall(MatDenseRestoreArray(B, &vb)); 2690 PetscCall(MatDenseRestoreArrayRead(A, &va)); 2691 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 2692 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 2693 PetscFunctionReturn(PETSC_SUCCESS); 2694 } 2695 2696 PetscErrorCode MatSetUp_SeqDense(Mat A) 2697 { 2698 PetscFunctionBegin; 2699 PetscCall(PetscLayoutSetUp(A->rmap)); 2700 PetscCall(PetscLayoutSetUp(A->cmap)); 2701 if (!A->preallocated) PetscCall(MatSeqDenseSetPreallocation(A, NULL)); 2702 PetscFunctionReturn(PETSC_SUCCESS); 2703 } 2704 2705 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2706 { 2707 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 2708 PetscInt i, j; 2709 PetscInt min = PetscMin(A->rmap->n, A->cmap->n); 2710 PetscScalar *aa; 2711 2712 PetscFunctionBegin; 2713 PetscCall(MatDenseGetArray(A, &aa)); 2714 for (j = 0; j < A->cmap->n; j++) { 2715 for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscConj(aa[i + j * mat->lda]); 2716 } 2717 PetscCall(MatDenseRestoreArray(A, &aa)); 2718 if (mat->tau) 2719 for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]); 2720 PetscFunctionReturn(PETSC_SUCCESS); 2721 } 2722 2723 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2724 { 2725 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 2726 PetscInt i, j; 2727 PetscScalar *aa; 2728 2729 PetscFunctionBegin; 2730 PetscCall(MatDenseGetArray(A, &aa)); 2731 for (j = 0; j < A->cmap->n; j++) { 2732 for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscRealPart(aa[i + j * mat->lda]); 2733 } 2734 PetscCall(MatDenseRestoreArray(A, &aa)); 2735 PetscFunctionReturn(PETSC_SUCCESS); 2736 } 2737 2738 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2739 { 2740 Mat_SeqDense *mat = (Mat_SeqDense *)A->data; 2741 PetscInt i, j; 2742 PetscScalar *aa; 2743 2744 PetscFunctionBegin; 2745 PetscCall(MatDenseGetArray(A, &aa)); 2746 for (j = 0; j < A->cmap->n; j++) { 2747 for (i = 0; i < A->rmap->n; i++) aa[i + j * mat->lda] = PetscImaginaryPart(aa[i + j * mat->lda]); 2748 } 2749 PetscCall(MatDenseRestoreArray(A, &aa)); 2750 PetscFunctionReturn(PETSC_SUCCESS); 2751 } 2752 2753 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C) 2754 { 2755 PetscInt m = A->rmap->n, n = B->cmap->n; 2756 PetscBool cisdense = PETSC_FALSE; 2757 2758 PetscFunctionBegin; 2759 PetscCall(MatSetSizes(C, m, n, m, n)); 2760 #if defined(PETSC_HAVE_CUDA) 2761 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, "")); 2762 #endif 2763 #if defined(PETSC_HAVE_HIP) 2764 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, "")); 2765 #endif 2766 if (!cisdense) { 2767 PetscBool flg; 2768 2769 PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg)); 2770 PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE)); 2771 } 2772 PetscCall(MatSetUp(C)); 2773 PetscFunctionReturn(PETSC_SUCCESS); 2774 } 2775 2776 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C) 2777 { 2778 Mat_SeqDense *a = (Mat_SeqDense *)A->data, *b = (Mat_SeqDense *)B->data, *c = (Mat_SeqDense *)C->data; 2779 PetscBLASInt m, n, k; 2780 const PetscScalar *av, *bv; 2781 PetscScalar *cv; 2782 PetscScalar _DOne = 1.0, _DZero = 0.0; 2783 2784 PetscFunctionBegin; 2785 PetscCall(PetscBLASIntCast(C->rmap->n, &m)); 2786 PetscCall(PetscBLASIntCast(C->cmap->n, &n)); 2787 PetscCall(PetscBLASIntCast(A->cmap->n, &k)); 2788 if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS); 2789 PetscCall(MatDenseGetArrayRead(A, &av)); 2790 PetscCall(MatDenseGetArrayRead(B, &bv)); 2791 PetscCall(MatDenseGetArrayWrite(C, &cv)); 2792 PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda)); 2793 PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1))); 2794 PetscCall(MatDenseRestoreArrayRead(A, &av)); 2795 PetscCall(MatDenseRestoreArrayRead(B, &bv)); 2796 PetscCall(MatDenseRestoreArrayWrite(C, &cv)); 2797 PetscFunctionReturn(PETSC_SUCCESS); 2798 } 2799 2800 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C) 2801 { 2802 PetscInt m = A->rmap->n, n = B->rmap->n; 2803 PetscBool cisdense = PETSC_FALSE; 2804 2805 PetscFunctionBegin; 2806 PetscCall(MatSetSizes(C, m, n, m, n)); 2807 #if defined(PETSC_HAVE_CUDA) 2808 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, "")); 2809 #endif 2810 #if defined(PETSC_HAVE_HIP) 2811 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, "")); 2812 #endif 2813 if (!cisdense) { 2814 PetscBool flg; 2815 2816 PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg)); 2817 PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE)); 2818 } 2819 PetscCall(MatSetUp(C)); 2820 PetscFunctionReturn(PETSC_SUCCESS); 2821 } 2822 2823 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C) 2824 { 2825 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 2826 Mat_SeqDense *b = (Mat_SeqDense *)B->data; 2827 Mat_SeqDense *c = (Mat_SeqDense *)C->data; 2828 const PetscScalar *av, *bv; 2829 PetscScalar *cv; 2830 PetscBLASInt m, n, k; 2831 PetscScalar _DOne = 1.0, _DZero = 0.0; 2832 2833 PetscFunctionBegin; 2834 PetscCall(PetscBLASIntCast(C->rmap->n, &m)); 2835 PetscCall(PetscBLASIntCast(C->cmap->n, &n)); 2836 PetscCall(PetscBLASIntCast(A->cmap->n, &k)); 2837 if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS); 2838 PetscCall(MatDenseGetArrayRead(A, &av)); 2839 PetscCall(MatDenseGetArrayRead(B, &bv)); 2840 PetscCall(MatDenseGetArrayWrite(C, &cv)); 2841 PetscCallBLAS("BLASgemm", BLASgemm_("N", "T", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda)); 2842 PetscCall(MatDenseRestoreArrayRead(A, &av)); 2843 PetscCall(MatDenseRestoreArrayRead(B, &bv)); 2844 PetscCall(MatDenseRestoreArrayWrite(C, &cv)); 2845 PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1))); 2846 PetscFunctionReturn(PETSC_SUCCESS); 2847 } 2848 2849 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A, Mat B, PetscReal fill, Mat C) 2850 { 2851 PetscInt m = A->cmap->n, n = B->cmap->n; 2852 PetscBool cisdense = PETSC_FALSE; 2853 2854 PetscFunctionBegin; 2855 PetscCall(MatSetSizes(C, m, n, m, n)); 2856 #if defined(PETSC_HAVE_CUDA) 2857 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSECUDA, "")); 2858 #endif 2859 #if defined(PETSC_HAVE_HIP) 2860 PetscCall(PetscObjectTypeCompareAny((PetscObject)C, &cisdense, MATSEQDENSE, MATSEQDENSEHIP, "")); 2861 #endif 2862 if (!cisdense) { 2863 PetscBool flg; 2864 2865 PetscCall(PetscObjectTypeCompare((PetscObject)B, ((PetscObject)A)->type_name, &flg)); 2866 PetscCall(MatSetType(C, flg ? ((PetscObject)A)->type_name : MATDENSE)); 2867 } 2868 PetscCall(MatSetUp(C)); 2869 PetscFunctionReturn(PETSC_SUCCESS); 2870 } 2871 2872 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A, Mat B, Mat C) 2873 { 2874 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 2875 Mat_SeqDense *b = (Mat_SeqDense *)B->data; 2876 Mat_SeqDense *c = (Mat_SeqDense *)C->data; 2877 const PetscScalar *av, *bv; 2878 PetscScalar *cv; 2879 PetscBLASInt m, n, k; 2880 PetscScalar _DOne = 1.0, _DZero = 0.0; 2881 2882 PetscFunctionBegin; 2883 PetscCall(PetscBLASIntCast(C->rmap->n, &m)); 2884 PetscCall(PetscBLASIntCast(C->cmap->n, &n)); 2885 PetscCall(PetscBLASIntCast(A->rmap->n, &k)); 2886 if (!m || !n || !k) PetscFunctionReturn(PETSC_SUCCESS); 2887 PetscCall(MatDenseGetArrayRead(A, &av)); 2888 PetscCall(MatDenseGetArrayRead(B, &bv)); 2889 PetscCall(MatDenseGetArrayWrite(C, &cv)); 2890 PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &m, &n, &k, &_DOne, av, &a->lda, bv, &b->lda, &_DZero, cv, &c->lda)); 2891 PetscCall(MatDenseRestoreArrayRead(A, &av)); 2892 PetscCall(MatDenseRestoreArrayRead(B, &bv)); 2893 PetscCall(MatDenseRestoreArrayWrite(C, &cv)); 2894 PetscCall(PetscLogFlops(1.0 * m * n * k + 1.0 * m * n * (k - 1))); 2895 PetscFunctionReturn(PETSC_SUCCESS); 2896 } 2897 2898 static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C) 2899 { 2900 PetscFunctionBegin; 2901 C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense; 2902 C->ops->productsymbolic = MatProductSymbolic_AB; 2903 PetscFunctionReturn(PETSC_SUCCESS); 2904 } 2905 2906 static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C) 2907 { 2908 PetscFunctionBegin; 2909 C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense; 2910 C->ops->productsymbolic = MatProductSymbolic_AtB; 2911 PetscFunctionReturn(PETSC_SUCCESS); 2912 } 2913 2914 static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C) 2915 { 2916 PetscFunctionBegin; 2917 C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense; 2918 C->ops->productsymbolic = MatProductSymbolic_ABt; 2919 PetscFunctionReturn(PETSC_SUCCESS); 2920 } 2921 2922 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C) 2923 { 2924 Mat_Product *product = C->product; 2925 2926 PetscFunctionBegin; 2927 switch (product->type) { 2928 case MATPRODUCT_AB: 2929 PetscCall(MatProductSetFromOptions_SeqDense_AB(C)); 2930 break; 2931 case MATPRODUCT_AtB: 2932 PetscCall(MatProductSetFromOptions_SeqDense_AtB(C)); 2933 break; 2934 case MATPRODUCT_ABt: 2935 PetscCall(MatProductSetFromOptions_SeqDense_ABt(C)); 2936 break; 2937 default: 2938 break; 2939 } 2940 PetscFunctionReturn(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 3102 } 3103 3104 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A, PetscBool *missing, PetscInt *d) 3105 { 3106 PetscFunctionBegin; 3107 *missing = PETSC_FALSE; 3108 PetscFunctionReturn(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 3130 } 3131 3132 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 3133 MatGetRow_SeqDense, 3134 MatRestoreRow_SeqDense, 3135 MatMult_SeqDense, 3136 /* 4*/ MatMultAdd_SeqDense, 3137 MatMultTranspose_SeqDense, 3138 MatMultTransposeAdd_SeqDense, 3139 NULL, 3140 NULL, 3141 NULL, 3142 /* 10*/ NULL, 3143 MatLUFactor_SeqDense, 3144 MatCholeskyFactor_SeqDense, 3145 MatSOR_SeqDense, 3146 MatTranspose_SeqDense, 3147 /* 15*/ MatGetInfo_SeqDense, 3148 MatEqual_SeqDense, 3149 MatGetDiagonal_SeqDense, 3150 MatDiagonalScale_SeqDense, 3151 MatNorm_SeqDense, 3152 /* 20*/ MatAssemblyBegin_SeqDense, 3153 MatAssemblyEnd_SeqDense, 3154 MatSetOption_SeqDense, 3155 MatZeroEntries_SeqDense, 3156 /* 24*/ MatZeroRows_SeqDense, 3157 NULL, 3158 NULL, 3159 NULL, 3160 NULL, 3161 /* 29*/ MatSetUp_SeqDense, 3162 NULL, 3163 NULL, 3164 NULL, 3165 NULL, 3166 /* 34*/ MatDuplicate_SeqDense, 3167 NULL, 3168 NULL, 3169 NULL, 3170 NULL, 3171 /* 39*/ MatAXPY_SeqDense, 3172 MatCreateSubMatrices_SeqDense, 3173 NULL, 3174 MatGetValues_SeqDense, 3175 MatCopy_SeqDense, 3176 /* 44*/ MatGetRowMax_SeqDense, 3177 MatScale_SeqDense, 3178 MatShift_SeqDense, 3179 NULL, 3180 MatZeroRowsColumns_SeqDense, 3181 /* 49*/ MatSetRandom_SeqDense, 3182 NULL, 3183 NULL, 3184 NULL, 3185 NULL, 3186 /* 54*/ NULL, 3187 NULL, 3188 NULL, 3189 NULL, 3190 NULL, 3191 /* 59*/ MatCreateSubMatrix_SeqDense, 3192 MatDestroy_SeqDense, 3193 MatView_SeqDense, 3194 NULL, 3195 NULL, 3196 /* 64*/ NULL, 3197 NULL, 3198 NULL, 3199 NULL, 3200 NULL, 3201 /* 69*/ MatGetRowMaxAbs_SeqDense, 3202 NULL, 3203 NULL, 3204 NULL, 3205 NULL, 3206 /* 74*/ NULL, 3207 NULL, 3208 NULL, 3209 NULL, 3210 NULL, 3211 /* 79*/ NULL, 3212 NULL, 3213 NULL, 3214 NULL, 3215 /* 83*/ MatLoad_SeqDense, 3216 MatIsSymmetric_SeqDense, 3217 MatIsHermitian_SeqDense, 3218 NULL, 3219 NULL, 3220 NULL, 3221 /* 89*/ NULL, 3222 NULL, 3223 MatMatMultNumeric_SeqDense_SeqDense, 3224 NULL, 3225 NULL, 3226 /* 94*/ NULL, 3227 NULL, 3228 NULL, 3229 MatMatTransposeMultNumeric_SeqDense_SeqDense, 3230 NULL, 3231 /* 99*/ MatProductSetFromOptions_SeqDense, 3232 NULL, 3233 NULL, 3234 MatConjugate_SeqDense, 3235 NULL, 3236 /*104*/ NULL, 3237 MatRealPart_SeqDense, 3238 MatImaginaryPart_SeqDense, 3239 NULL, 3240 NULL, 3241 /*109*/ NULL, 3242 NULL, 3243 MatGetRowMin_SeqDense, 3244 MatGetColumnVector_SeqDense, 3245 MatMissingDiagonal_SeqDense, 3246 /*114*/ NULL, 3247 NULL, 3248 NULL, 3249 NULL, 3250 NULL, 3251 /*119*/ NULL, 3252 NULL, 3253 NULL, 3254 NULL, 3255 NULL, 3256 /*124*/ NULL, 3257 MatGetColumnReductions_SeqDense, 3258 NULL, 3259 NULL, 3260 NULL, 3261 /*129*/ NULL, 3262 NULL, 3263 NULL, 3264 MatTransposeMatMultNumeric_SeqDense_SeqDense, 3265 NULL, 3266 /*134*/ NULL, 3267 NULL, 3268 NULL, 3269 NULL, 3270 NULL, 3271 /*139*/ NULL, 3272 NULL, 3273 NULL, 3274 NULL, 3275 NULL, 3276 MatCreateMPIMatConcatenateSeqMat_SeqDense, 3277 /*145*/ NULL, 3278 NULL, 3279 NULL, 3280 NULL, 3281 NULL, 3282 /*150*/ NULL, 3283 NULL}; 3284 3285 /*@C 3286 MatCreateSeqDense - Creates a `MATSEQDENSE` that 3287 is stored in column major order (the usual Fortran manner). Many 3288 of the matrix operations use the BLAS and LAPACK routines. 3289 3290 Collective 3291 3292 Input Parameters: 3293 + comm - MPI communicator, set to `PETSC_COMM_SELF` 3294 . m - number of rows 3295 . n - number of columns 3296 - data - optional location of matrix data in column major order. Use `NULL` for PETSc 3297 to control all matrix memory allocation. 3298 3299 Output Parameter: 3300 . A - the matrix 3301 3302 Level: intermediate 3303 3304 Note: 3305 The data input variable is intended primarily for Fortran programmers 3306 who wish to allocate their own matrix memory space. Most users should 3307 set `data` = `NULL`. 3308 3309 .seealso: [](chapter_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()` 3310 @*/ 3311 PetscErrorCode MatCreateSeqDense(MPI_Comm comm, PetscInt m, PetscInt n, PetscScalar *data, Mat *A) 3312 { 3313 PetscFunctionBegin; 3314 PetscCall(MatCreate(comm, A)); 3315 PetscCall(MatSetSizes(*A, m, n, m, n)); 3316 PetscCall(MatSetType(*A, MATSEQDENSE)); 3317 PetscCall(MatSeqDenseSetPreallocation(*A, data)); 3318 PetscFunctionReturn(PETSC_SUCCESS); 3319 } 3320 3321 /*@C 3322 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements of a `MATSEQDENSE` matrix 3323 3324 Collective 3325 3326 Input Parameters: 3327 + B - the matrix 3328 - data - the array (or `NULL`) 3329 3330 Level: intermediate 3331 3332 Note: 3333 The data input variable is intended primarily for Fortran programmers 3334 who wish to allocate their own matrix memory space. Most users should 3335 need not call this routine. 3336 3337 .seealso: [](chapter_matrices), `Mat`, `MATSEQDENSE`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`, `MatDenseSetLDA()` 3338 @*/ 3339 PetscErrorCode MatSeqDenseSetPreallocation(Mat B, PetscScalar data[]) 3340 { 3341 PetscFunctionBegin; 3342 PetscValidHeaderSpecific(B, MAT_CLASSID, 1); 3343 PetscTryMethod(B, "MatSeqDenseSetPreallocation_C", (Mat, PetscScalar[]), (B, data)); 3344 PetscFunctionReturn(PETSC_SUCCESS); 3345 } 3346 3347 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B, PetscScalar *data) 3348 { 3349 Mat_SeqDense *b = (Mat_SeqDense *)B->data; 3350 3351 PetscFunctionBegin; 3352 PetscCheck(!b->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 3353 B->preallocated = PETSC_TRUE; 3354 3355 PetscCall(PetscLayoutSetUp(B->rmap)); 3356 PetscCall(PetscLayoutSetUp(B->cmap)); 3357 3358 if (b->lda <= 0) b->lda = B->rmap->n; 3359 3360 if (!data) { /* petsc-allocated storage */ 3361 if (!b->user_alloc) PetscCall(PetscFree(b->v)); 3362 PetscCall(PetscCalloc1((size_t)b->lda * B->cmap->n, &b->v)); 3363 3364 b->user_alloc = PETSC_FALSE; 3365 } else { /* user-allocated storage */ 3366 if (!b->user_alloc) PetscCall(PetscFree(b->v)); 3367 b->v = data; 3368 b->user_alloc = PETSC_TRUE; 3369 } 3370 B->assembled = PETSC_TRUE; 3371 PetscFunctionReturn(PETSC_SUCCESS); 3372 } 3373 3374 #if defined(PETSC_HAVE_ELEMENTAL) 3375 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 3376 { 3377 Mat mat_elemental; 3378 const PetscScalar *array; 3379 PetscScalar *v_colwise; 3380 PetscInt M = A->rmap->N, N = A->cmap->N, i, j, k, *rows, *cols; 3381 3382 PetscFunctionBegin; 3383 PetscCall(PetscMalloc3(M * N, &v_colwise, M, &rows, N, &cols)); 3384 PetscCall(MatDenseGetArrayRead(A, &array)); 3385 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 3386 k = 0; 3387 for (j = 0; j < N; j++) { 3388 cols[j] = j; 3389 for (i = 0; i < M; i++) v_colwise[j * M + i] = array[k++]; 3390 } 3391 for (i = 0; i < M; i++) rows[i] = i; 3392 PetscCall(MatDenseRestoreArrayRead(A, &array)); 3393 3394 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental)); 3395 PetscCall(MatSetSizes(mat_elemental, PETSC_DECIDE, PETSC_DECIDE, M, N)); 3396 PetscCall(MatSetType(mat_elemental, MATELEMENTAL)); 3397 PetscCall(MatSetUp(mat_elemental)); 3398 3399 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 3400 PetscCall(MatSetValues(mat_elemental, M, rows, N, cols, v_colwise, ADD_VALUES)); 3401 PetscCall(MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY)); 3402 PetscCall(MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY)); 3403 PetscCall(PetscFree3(v_colwise, rows, cols)); 3404 3405 if (reuse == MAT_INPLACE_MATRIX) { 3406 PetscCall(MatHeaderReplace(A, &mat_elemental)); 3407 } else { 3408 *newmat = mat_elemental; 3409 } 3410 PetscFunctionReturn(PETSC_SUCCESS); 3411 } 3412 #endif 3413 3414 PetscErrorCode MatDenseSetLDA_SeqDense(Mat B, PetscInt lda) 3415 { 3416 Mat_SeqDense *b = (Mat_SeqDense *)B->data; 3417 PetscBool data; 3418 3419 PetscFunctionBegin; 3420 data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE); 3421 PetscCheck(b->user_alloc || !data || b->lda == lda, PETSC_COMM_SELF, PETSC_ERR_ORDER, "LDA cannot be changed after allocation of internal storage"); 3422 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); 3423 b->lda = lda; 3424 PetscFunctionReturn(PETSC_SUCCESS); 3425 } 3426 3427 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat) 3428 { 3429 PetscFunctionBegin; 3430 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIDense(comm, inmat, n, scall, outmat)); 3431 PetscFunctionReturn(PETSC_SUCCESS); 3432 } 3433 3434 PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A, PetscInt col, Vec *v) 3435 { 3436 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3437 3438 PetscFunctionBegin; 3439 PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first"); 3440 PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 3441 if (!a->cvec) { PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, &a->cvec)); } 3442 a->vecinuse = col + 1; 3443 PetscCall(MatDenseGetArray(A, (PetscScalar **)&a->ptrinuse)); 3444 PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda)); 3445 *v = a->cvec; 3446 PetscFunctionReturn(PETSC_SUCCESS); 3447 } 3448 3449 PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A, PetscInt col, Vec *v) 3450 { 3451 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3452 3453 PetscFunctionBegin; 3454 PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first"); 3455 PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector"); 3456 a->vecinuse = 0; 3457 PetscCall(MatDenseRestoreArray(A, (PetscScalar **)&a->ptrinuse)); 3458 PetscCall(VecResetArray(a->cvec)); 3459 if (v) *v = NULL; 3460 PetscFunctionReturn(PETSC_SUCCESS); 3461 } 3462 3463 PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v) 3464 { 3465 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3466 3467 PetscFunctionBegin; 3468 PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first"); 3469 PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 3470 if (!a->cvec) { PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, &a->cvec)); } 3471 a->vecinuse = col + 1; 3472 PetscCall(MatDenseGetArrayRead(A, &a->ptrinuse)); 3473 PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda)); 3474 PetscCall(VecLockReadPush(a->cvec)); 3475 *v = a->cvec; 3476 PetscFunctionReturn(PETSC_SUCCESS); 3477 } 3478 3479 PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A, PetscInt col, Vec *v) 3480 { 3481 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3482 3483 PetscFunctionBegin; 3484 PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first"); 3485 PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector"); 3486 a->vecinuse = 0; 3487 PetscCall(MatDenseRestoreArrayRead(A, &a->ptrinuse)); 3488 PetscCall(VecLockReadPop(a->cvec)); 3489 PetscCall(VecResetArray(a->cvec)); 3490 if (v) *v = NULL; 3491 PetscFunctionReturn(PETSC_SUCCESS); 3492 } 3493 3494 PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v) 3495 { 3496 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3497 3498 PetscFunctionBegin; 3499 PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first"); 3500 PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 3501 if (!a->cvec) PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)A), A->rmap->bs, A->rmap->n, NULL, &a->cvec)); 3502 a->vecinuse = col + 1; 3503 PetscCall(MatDenseGetArrayWrite(A, (PetscScalar **)&a->ptrinuse)); 3504 PetscCall(VecPlaceArray(a->cvec, a->ptrinuse + (size_t)col * (size_t)a->lda)); 3505 *v = a->cvec; 3506 PetscFunctionReturn(PETSC_SUCCESS); 3507 } 3508 3509 PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A, PetscInt col, Vec *v) 3510 { 3511 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3512 3513 PetscFunctionBegin; 3514 PetscCheck(a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetColumnVec() first"); 3515 PetscCheck(a->cvec, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column vector"); 3516 a->vecinuse = 0; 3517 PetscCall(MatDenseRestoreArrayWrite(A, (PetscScalar **)&a->ptrinuse)); 3518 PetscCall(VecResetArray(a->cvec)); 3519 if (v) *v = NULL; 3520 PetscFunctionReturn(PETSC_SUCCESS); 3521 } 3522 3523 PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v) 3524 { 3525 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3526 3527 PetscFunctionBegin; 3528 PetscCheck(!a->vecinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreColumnVec() first"); 3529 PetscCheck(!a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseRestoreSubMatrix() first"); 3530 if (a->cmat && (cend - cbegin != a->cmat->cmap->N || rend - rbegin != a->cmat->rmap->N)) PetscCall(MatDestroy(&a->cmat)); 3531 if (!a->cmat) { 3532 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), rend - rbegin, PETSC_DECIDE, rend - rbegin, cend - cbegin, a->v + rbegin + (size_t)cbegin * a->lda, &a->cmat)); 3533 } else { 3534 PetscCall(MatDensePlaceArray(a->cmat, a->v + rbegin + (size_t)cbegin * a->lda)); 3535 } 3536 PetscCall(MatDenseSetLDA(a->cmat, a->lda)); 3537 a->matinuse = cbegin + 1; 3538 *v = a->cmat; 3539 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP) 3540 A->offloadmask = PETSC_OFFLOAD_CPU; 3541 #endif 3542 PetscFunctionReturn(PETSC_SUCCESS); 3543 } 3544 3545 PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A, Mat *v) 3546 { 3547 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3548 3549 PetscFunctionBegin; 3550 PetscCheck(a->matinuse, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Need to call MatDenseGetSubMatrix() first"); 3551 PetscCheck(a->cmat, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing internal column matrix"); 3552 PetscCheck(*v == a->cmat, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not the matrix obtained from MatDenseGetSubMatrix()"); 3553 a->matinuse = 0; 3554 PetscCall(MatDenseResetArray(a->cmat)); 3555 if (v) *v = NULL; 3556 #if defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP) 3557 A->offloadmask = PETSC_OFFLOAD_CPU; 3558 #endif 3559 PetscFunctionReturn(PETSC_SUCCESS); 3560 } 3561 3562 /*MC 3563 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 3564 3565 Options Database Key: 3566 . -mat_type seqdense - sets the matrix type to `MATSEQDENSE` during a call to `MatSetFromOptions()` 3567 3568 Level: beginner 3569 3570 .seealso: [](chapter_matrices), `Mat`, `MATSEQDENSE`, `MatCreateSeqDense()` 3571 M*/ 3572 PetscErrorCode MatCreate_SeqDense(Mat B) 3573 { 3574 Mat_SeqDense *b; 3575 PetscMPIInt size; 3576 3577 PetscFunctionBegin; 3578 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size)); 3579 PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1"); 3580 3581 PetscCall(PetscNew(&b)); 3582 PetscCall(PetscMemcpy(B->ops, &MatOps_Values, sizeof(struct _MatOps))); 3583 B->data = (void *)b; 3584 3585 b->roworiented = PETSC_TRUE; 3586 3587 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatQRFactor_C", MatQRFactor_SeqDense)); 3588 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetLDA_C", MatDenseGetLDA_SeqDense)); 3589 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseSetLDA_C", MatDenseSetLDA_SeqDense)); 3590 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArray_C", MatDenseGetArray_SeqDense)); 3591 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArray_C", MatDenseRestoreArray_SeqDense)); 3592 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDensePlaceArray_C", MatDensePlaceArray_SeqDense)); 3593 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseResetArray_C", MatDenseResetArray_SeqDense)); 3594 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseReplaceArray_C", MatDenseReplaceArray_SeqDense)); 3595 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayRead_C", MatDenseGetArray_SeqDense)); 3596 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayRead_C", MatDenseRestoreArray_SeqDense)); 3597 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetArrayWrite_C", MatDenseGetArray_SeqDense)); 3598 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreArrayWrite_C", MatDenseRestoreArray_SeqDense)); 3599 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqaij_C", MatConvert_SeqDense_SeqAIJ)); 3600 #if defined(PETSC_HAVE_ELEMENTAL) 3601 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_elemental_C", MatConvert_SeqDense_Elemental)); 3602 #endif 3603 #if defined(PETSC_HAVE_SCALAPACK) 3604 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_scalapack_C", MatConvert_Dense_ScaLAPACK)); 3605 #endif 3606 #if defined(PETSC_HAVE_CUDA) 3607 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensecuda_C", MatConvert_SeqDense_SeqDenseCUDA)); 3608 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdensecuda_C", MatProductSetFromOptions_SeqDense)); 3609 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensecuda_seqdense_C", MatProductSetFromOptions_SeqDense)); 3610 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensecuda_C", MatProductSetFromOptions_SeqDense)); 3611 #endif 3612 #if defined(PETSC_HAVE_HIP) 3613 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqdense_seqdensehip_C", MatConvert_SeqDense_SeqDenseHIP)); 3614 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdensehip_C", MatProductSetFromOptions_SeqDense)); 3615 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdensehip_seqdense_C", MatProductSetFromOptions_SeqDense)); 3616 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdensehip_C", MatProductSetFromOptions_SeqDense)); 3617 #endif 3618 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqDenseSetPreallocation_C", MatSeqDenseSetPreallocation_SeqDense)); 3619 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_seqdense_C", MatProductSetFromOptions_SeqAIJ_SeqDense)); 3620 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqdense_seqdense_C", MatProductSetFromOptions_SeqDense)); 3621 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense)); 3622 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqsbaij_seqdense_C", MatProductSetFromOptions_SeqXBAIJ_SeqDense)); 3623 3624 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumn_C", MatDenseGetColumn_SeqDense)); 3625 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumn_C", MatDenseRestoreColumn_SeqDense)); 3626 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVec_C", MatDenseGetColumnVec_SeqDense)); 3627 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVec_C", MatDenseRestoreColumnVec_SeqDense)); 3628 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecRead_C", MatDenseGetColumnVecRead_SeqDense)); 3629 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecRead_C", MatDenseRestoreColumnVecRead_SeqDense)); 3630 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetColumnVecWrite_C", MatDenseGetColumnVecWrite_SeqDense)); 3631 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreColumnVecWrite_C", MatDenseRestoreColumnVecWrite_SeqDense)); 3632 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseGetSubMatrix_C", MatDenseGetSubMatrix_SeqDense)); 3633 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatDenseRestoreSubMatrix_C", MatDenseRestoreSubMatrix_SeqDense)); 3634 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQDENSE)); 3635 PetscFunctionReturn(PETSC_SUCCESS); 3636 } 3637 3638 /*@C 3639 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. 3640 3641 Not Collective 3642 3643 Input Parameters: 3644 + mat - a `MATSEQDENSE` or `MATMPIDENSE` matrix 3645 - col - column index 3646 3647 Output Parameter: 3648 . vals - pointer to the data 3649 3650 Level: intermediate 3651 3652 Note: 3653 Use `MatDenseGetColumnVec()` to get access to a column of a `MATDENSE` treated as a `Vec` 3654 3655 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MatDenseRestoreColumn()`, `MatDenseGetColumnVec()` 3656 @*/ 3657 PetscErrorCode MatDenseGetColumn(Mat A, PetscInt col, PetscScalar **vals) 3658 { 3659 PetscFunctionBegin; 3660 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3661 PetscValidLogicalCollectiveInt(A, col, 2); 3662 PetscValidPointer(vals, 3); 3663 PetscUseMethod(A, "MatDenseGetColumn_C", (Mat, PetscInt, PetscScalar **), (A, col, vals)); 3664 PetscFunctionReturn(PETSC_SUCCESS); 3665 } 3666 3667 /*@C 3668 MatDenseRestoreColumn - returns access to a column of a `MATDENSE` matrix which is returned by `MatDenseGetColumn()`. 3669 3670 Not Collective 3671 3672 Input Parameters: 3673 + mat - a `MATSEQDENSE` or `MATMPIDENSE` matrix 3674 - vals - pointer to the data (may be `NULL`) 3675 3676 Level: intermediate 3677 3678 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MatDenseGetColumn()` 3679 @*/ 3680 PetscErrorCode MatDenseRestoreColumn(Mat A, PetscScalar **vals) 3681 { 3682 PetscFunctionBegin; 3683 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3684 PetscValidPointer(vals, 2); 3685 PetscUseMethod(A, "MatDenseRestoreColumn_C", (Mat, PetscScalar **), (A, vals)); 3686 PetscFunctionReturn(PETSC_SUCCESS); 3687 } 3688 3689 /*@ 3690 MatDenseGetColumnVec - Gives read-write access to a column of a `MATDENSE` matrix, represented as a `Vec`. 3691 3692 Collective 3693 3694 Input Parameters: 3695 + mat - the `Mat` object 3696 - col - the column index 3697 3698 Output Parameter: 3699 . v - the vector 3700 3701 Level: intermediate 3702 3703 Notes: 3704 The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVec()` when the vector is no longer needed. 3705 3706 Use `MatDenseGetColumnVecRead()` to obtain read-only access or `MatDenseGetColumnVecWrite()` for write-only access. 3707 3708 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()`, `MatDenseGetColumn()` 3709 @*/ 3710 PetscErrorCode MatDenseGetColumnVec(Mat A, PetscInt col, Vec *v) 3711 { 3712 PetscFunctionBegin; 3713 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3714 PetscValidType(A, 1); 3715 PetscValidLogicalCollectiveInt(A, col, 2); 3716 PetscValidPointer(v, 3); 3717 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3718 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); 3719 PetscUseMethod(A, "MatDenseGetColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v)); 3720 PetscFunctionReturn(PETSC_SUCCESS); 3721 } 3722 3723 /*@ 3724 MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec(). 3725 3726 Collective 3727 3728 Input Parameters: 3729 + mat - the Mat object 3730 . col - the column index 3731 - v - the Vec object (may be `NULL`) 3732 3733 Level: intermediate 3734 3735 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3736 @*/ 3737 PetscErrorCode MatDenseRestoreColumnVec(Mat A, PetscInt col, Vec *v) 3738 { 3739 PetscFunctionBegin; 3740 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3741 PetscValidType(A, 1); 3742 PetscValidLogicalCollectiveInt(A, col, 2); 3743 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3744 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); 3745 PetscUseMethod(A, "MatDenseRestoreColumnVec_C", (Mat, PetscInt, Vec *), (A, col, v)); 3746 PetscFunctionReturn(PETSC_SUCCESS); 3747 } 3748 3749 /*@ 3750 MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec. 3751 3752 Collective 3753 3754 Input Parameters: 3755 + mat - the `Mat` object 3756 - col - the column index 3757 3758 Output Parameter: 3759 . v - the vector 3760 3761 Level: intermediate 3762 3763 Notes: 3764 The vector is owned by PETSc and users cannot modify it. 3765 3766 Users need to call `MatDenseRestoreColumnVecRead()` when the vector is no longer needed. 3767 3768 Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecWrite()` for write-only access. 3769 3770 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3771 @*/ 3772 PetscErrorCode MatDenseGetColumnVecRead(Mat A, PetscInt col, Vec *v) 3773 { 3774 PetscFunctionBegin; 3775 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3776 PetscValidType(A, 1); 3777 PetscValidLogicalCollectiveInt(A, col, 2); 3778 PetscValidPointer(v, 3); 3779 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3780 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); 3781 PetscUseMethod(A, "MatDenseGetColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v)); 3782 PetscFunctionReturn(PETSC_SUCCESS); 3783 } 3784 3785 /*@ 3786 MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead(). 3787 3788 Collective 3789 3790 Input Parameters: 3791 + mat - the `Mat` object 3792 . col - the column index 3793 - v - the Vec object (may be `NULL`) 3794 3795 Level: intermediate 3796 3797 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecWrite()` 3798 @*/ 3799 PetscErrorCode MatDenseRestoreColumnVecRead(Mat A, PetscInt col, Vec *v) 3800 { 3801 PetscFunctionBegin; 3802 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3803 PetscValidType(A, 1); 3804 PetscValidLogicalCollectiveInt(A, col, 2); 3805 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3806 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); 3807 PetscUseMethod(A, "MatDenseRestoreColumnVecRead_C", (Mat, PetscInt, Vec *), (A, col, v)); 3808 PetscFunctionReturn(PETSC_SUCCESS); 3809 } 3810 3811 /*@ 3812 MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec. 3813 3814 Collective 3815 3816 Input Parameters: 3817 + mat - the `Mat` object 3818 - col - the column index 3819 3820 Output Parameter: 3821 . v - the vector 3822 3823 Level: intermediate 3824 3825 Notes: 3826 The vector is owned by PETSc. Users need to call `MatDenseRestoreColumnVecWrite()` when the vector is no longer needed. 3827 3828 Use `MatDenseGetColumnVec()` to obtain read-write access or `MatDenseGetColumnVecRead()` for read-only access. 3829 3830 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()`, `MatDenseRestoreColumnVecWrite()` 3831 @*/ 3832 PetscErrorCode MatDenseGetColumnVecWrite(Mat A, PetscInt col, Vec *v) 3833 { 3834 PetscFunctionBegin; 3835 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3836 PetscValidType(A, 1); 3837 PetscValidLogicalCollectiveInt(A, col, 2); 3838 PetscValidPointer(v, 3); 3839 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3840 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); 3841 PetscUseMethod(A, "MatDenseGetColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v)); 3842 PetscFunctionReturn(PETSC_SUCCESS); 3843 } 3844 3845 /*@ 3846 MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite(). 3847 3848 Collective 3849 3850 Input Parameters: 3851 + mat - the `Mat` object 3852 . col - the column index 3853 - v - the `Vec` object (may be `NULL`) 3854 3855 Level: intermediate 3856 3857 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseGetColumnVecRead()`, `MatDenseGetColumnVecWrite()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreColumnVecRead()` 3858 @*/ 3859 PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A, PetscInt col, Vec *v) 3860 { 3861 PetscFunctionBegin; 3862 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3863 PetscValidType(A, 1); 3864 PetscValidLogicalCollectiveInt(A, col, 2); 3865 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3866 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); 3867 PetscUseMethod(A, "MatDenseRestoreColumnVecWrite_C", (Mat, PetscInt, Vec *), (A, col, v)); 3868 PetscFunctionReturn(PETSC_SUCCESS); 3869 } 3870 3871 /*@ 3872 MatDenseGetSubMatrix - Gives access to a block of rows and columns of a dense matrix, represented as a Mat. 3873 3874 Collective 3875 3876 Input Parameters: 3877 + mat - the Mat object 3878 . rbegin - the first global row index in the block (if `PETSC_DECIDE`, is 0) 3879 . rend - the global row index past the last one in the block (if `PETSC_DECIDE`, is `M`) 3880 . cbegin - the first global column index in the block (if `PETSC_DECIDE`, is 0) 3881 - cend - the global column index past the last one in the block (if `PETSC_DECIDE`, is `N`) 3882 3883 Output Parameter: 3884 . v - the matrix 3885 3886 Level: intermediate 3887 3888 Notes: 3889 The matrix is owned by PETSc. Users need to call `MatDenseRestoreSubMatrix()` when the matrix is no longer needed. 3890 3891 The output matrix is not redistributed by PETSc, so depending on the values of `rbegin` and `rend`, some processes may have no local rows. 3892 3893 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseRestoreSubMatrix()` 3894 @*/ 3895 PetscErrorCode MatDenseGetSubMatrix(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *v) 3896 { 3897 PetscFunctionBegin; 3898 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3899 PetscValidType(A, 1); 3900 PetscValidLogicalCollectiveInt(A, rbegin, 2); 3901 PetscValidLogicalCollectiveInt(A, rend, 3); 3902 PetscValidLogicalCollectiveInt(A, cbegin, 4); 3903 PetscValidLogicalCollectiveInt(A, cend, 5); 3904 PetscValidPointer(v, 6); 3905 if (rbegin == PETSC_DECIDE) rbegin = 0; 3906 if (rend == PETSC_DECIDE) rend = A->rmap->N; 3907 if (cbegin == PETSC_DECIDE) cbegin = 0; 3908 if (cend == PETSC_DECIDE) cend = A->cmap->N; 3909 PetscCheck(A->preallocated, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Matrix not preallocated"); 3910 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); 3911 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); 3912 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); 3913 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); 3914 PetscUseMethod(A, "MatDenseGetSubMatrix_C", (Mat, PetscInt, PetscInt, PetscInt, PetscInt, Mat *), (A, rbegin, rend, cbegin, cend, v)); 3915 PetscFunctionReturn(PETSC_SUCCESS); 3916 } 3917 3918 /*@ 3919 MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix(). 3920 3921 Collective 3922 3923 Input Parameters: 3924 + mat - the `Mat` object 3925 - v - the `Mat` object (may be `NULL`) 3926 3927 Level: intermediate 3928 3929 .seealso: [](chapter_matrices), `Mat`, `MATDENSE`, `MATDENSECUDA`, `MATDENSEHIP`, `MatDenseGetColumnVec()`, `MatDenseRestoreColumnVec()`, `MatDenseGetSubMatrix()` 3930 @*/ 3931 PetscErrorCode MatDenseRestoreSubMatrix(Mat A, Mat *v) 3932 { 3933 PetscFunctionBegin; 3934 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 3935 PetscValidType(A, 1); 3936 PetscValidPointer(v, 2); 3937 PetscUseMethod(A, "MatDenseRestoreSubMatrix_C", (Mat, Mat *), (A, v)); 3938 PetscFunctionReturn(PETSC_SUCCESS); 3939 } 3940 3941 #include <petscblaslapack.h> 3942 #include <petsc/private/kernels/blockinvert.h> 3943 3944 PetscErrorCode MatSeqDenseInvert(Mat A) 3945 { 3946 Mat_SeqDense *a = (Mat_SeqDense *)A->data; 3947 PetscInt bs = A->rmap->n; 3948 MatScalar *values = a->v; 3949 const PetscReal shift = 0.0; 3950 PetscBool allowzeropivot = PetscNot(A->erroriffailure), zeropivotdetected = PETSC_FALSE; 3951 3952 PetscFunctionBegin; 3953 /* factor and invert each block */ 3954 switch (bs) { 3955 case 1: 3956 values[0] = (PetscScalar)1.0 / (values[0] + shift); 3957 break; 3958 case 2: 3959 PetscCall(PetscKernel_A_gets_inverse_A_2(values, shift, allowzeropivot, &zeropivotdetected)); 3960 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3961 break; 3962 case 3: 3963 PetscCall(PetscKernel_A_gets_inverse_A_3(values, shift, allowzeropivot, &zeropivotdetected)); 3964 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3965 break; 3966 case 4: 3967 PetscCall(PetscKernel_A_gets_inverse_A_4(values, shift, allowzeropivot, &zeropivotdetected)); 3968 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3969 break; 3970 case 5: { 3971 PetscScalar work[25]; 3972 PetscInt ipvt[5]; 3973 3974 PetscCall(PetscKernel_A_gets_inverse_A_5(values, ipvt, work, shift, allowzeropivot, &zeropivotdetected)); 3975 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3976 } break; 3977 case 6: 3978 PetscCall(PetscKernel_A_gets_inverse_A_6(values, shift, allowzeropivot, &zeropivotdetected)); 3979 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3980 break; 3981 case 7: 3982 PetscCall(PetscKernel_A_gets_inverse_A_7(values, shift, allowzeropivot, &zeropivotdetected)); 3983 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3984 break; 3985 default: { 3986 PetscInt *v_pivots, *IJ, j; 3987 PetscScalar *v_work; 3988 3989 PetscCall(PetscMalloc3(bs, &v_work, bs, &v_pivots, bs, &IJ)); 3990 for (j = 0; j < bs; j++) IJ[j] = j; 3991 PetscCall(PetscKernel_A_gets_inverse_A(bs, values, v_pivots, v_work, allowzeropivot, &zeropivotdetected)); 3992 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3993 PetscCall(PetscFree3(v_work, v_pivots, IJ)); 3994 } 3995 } 3996 PetscFunctionReturn(PETSC_SUCCESS); 3997 } 3998