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