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