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