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