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