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