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