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