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