1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 typedef struct { 5 Mat A; 6 Mat D; /* local submatrix for diagonal part */ 7 Vec w, left, right, leftwork, rightwork; 8 PetscScalar scale; 9 } Mat_Normal; 10 11 PetscErrorCode MatScale_Normal(Mat inA, PetscScalar scale) { 12 Mat_Normal *a = (Mat_Normal *)inA->data; 13 14 PetscFunctionBegin; 15 a->scale *= scale; 16 PetscFunctionReturn(0); 17 } 18 19 PetscErrorCode MatDiagonalScale_Normal(Mat inA, Vec left, Vec right) { 20 Mat_Normal *a = (Mat_Normal *)inA->data; 21 22 PetscFunctionBegin; 23 if (left) { 24 if (!a->left) { 25 PetscCall(VecDuplicate(left, &a->left)); 26 PetscCall(VecCopy(left, a->left)); 27 } else { 28 PetscCall(VecPointwiseMult(a->left, left, a->left)); 29 } 30 } 31 if (right) { 32 if (!a->right) { 33 PetscCall(VecDuplicate(right, &a->right)); 34 PetscCall(VecCopy(right, a->right)); 35 } else { 36 PetscCall(VecPointwiseMult(a->right, right, a->right)); 37 } 38 } 39 PetscFunctionReturn(0); 40 } 41 42 PetscErrorCode MatIncreaseOverlap_Normal(Mat A, PetscInt is_max, IS is[], PetscInt ov) { 43 Mat_Normal *a = (Mat_Normal *)A->data; 44 Mat pattern; 45 46 PetscFunctionBegin; 47 PetscCheck(ov >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Negative overlap specified"); 48 PetscCall(MatProductCreate(a->A, a->A, NULL, &pattern)); 49 PetscCall(MatProductSetType(pattern, MATPRODUCT_AtB)); 50 PetscCall(MatProductSetFromOptions(pattern)); 51 PetscCall(MatProductSymbolic(pattern)); 52 PetscCall(MatIncreaseOverlap(pattern, is_max, is, ov)); 53 PetscCall(MatDestroy(&pattern)); 54 PetscFunctionReturn(0); 55 } 56 57 PetscErrorCode MatCreateSubMatrices_Normal(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[]) { 58 Mat_Normal *a = (Mat_Normal *)mat->data; 59 Mat B = a->A, *suba; 60 IS *row; 61 PetscInt M; 62 63 PetscFunctionBegin; 64 PetscCheck(!a->left && !a->right && irow == icol, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Not implemented"); 65 if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat)); 66 PetscCall(MatGetSize(B, &M, NULL)); 67 PetscCall(PetscMalloc1(n, &row)); 68 PetscCall(ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0])); 69 PetscCall(ISSetIdentity(row[0])); 70 for (M = 1; M < n; ++M) row[M] = row[0]; 71 PetscCall(MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba)); 72 for (M = 0; M < n; ++M) { 73 PetscCall(MatCreateNormal(suba[M], *submat + M)); 74 ((Mat_Normal *)(*submat)[M]->data)->scale = a->scale; 75 } 76 PetscCall(ISDestroy(&row[0])); 77 PetscCall(PetscFree(row)); 78 PetscCall(MatDestroySubMatrices(n, &suba)); 79 PetscFunctionReturn(0); 80 } 81 82 PetscErrorCode MatPermute_Normal(Mat A, IS rowp, IS colp, Mat *B) { 83 Mat_Normal *a = (Mat_Normal *)A->data; 84 Mat C, Aa = a->A; 85 IS row; 86 87 PetscFunctionBegin; 88 PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same"); 89 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row)); 90 PetscCall(ISSetIdentity(row)); 91 PetscCall(MatPermute(Aa, row, colp, &C)); 92 PetscCall(ISDestroy(&row)); 93 PetscCall(MatCreateNormal(C, B)); 94 PetscCall(MatDestroy(&C)); 95 PetscFunctionReturn(0); 96 } 97 98 PetscErrorCode MatDuplicate_Normal(Mat A, MatDuplicateOption op, Mat *B) { 99 Mat_Normal *a = (Mat_Normal *)A->data; 100 Mat C; 101 102 PetscFunctionBegin; 103 PetscCheck(!a->left && !a->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented"); 104 PetscCall(MatDuplicate(a->A, op, &C)); 105 PetscCall(MatCreateNormal(C, B)); 106 PetscCall(MatDestroy(&C)); 107 if (op == MAT_COPY_VALUES) ((Mat_Normal *)(*B)->data)->scale = a->scale; 108 PetscFunctionReturn(0); 109 } 110 111 PetscErrorCode MatCopy_Normal(Mat A, Mat B, MatStructure str) { 112 Mat_Normal *a = (Mat_Normal *)A->data, *b = (Mat_Normal *)B->data; 113 114 PetscFunctionBegin; 115 PetscCheck(!a->left && !a->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented"); 116 PetscCall(MatCopy(a->A, b->A, str)); 117 b->scale = a->scale; 118 PetscCall(VecDestroy(&b->left)); 119 PetscCall(VecDestroy(&b->right)); 120 PetscCall(VecDestroy(&b->leftwork)); 121 PetscCall(VecDestroy(&b->rightwork)); 122 PetscFunctionReturn(0); 123 } 124 125 PetscErrorCode MatMult_Normal(Mat N, Vec x, Vec y) { 126 Mat_Normal *Na = (Mat_Normal *)N->data; 127 Vec in; 128 129 PetscFunctionBegin; 130 in = x; 131 if (Na->right) { 132 if (!Na->rightwork) PetscCall(VecDuplicate(Na->right, &Na->rightwork)); 133 PetscCall(VecPointwiseMult(Na->rightwork, Na->right, in)); 134 in = Na->rightwork; 135 } 136 PetscCall(MatMult(Na->A, in, Na->w)); 137 PetscCall(MatMultTranspose(Na->A, Na->w, y)); 138 if (Na->left) PetscCall(VecPointwiseMult(y, Na->left, y)); 139 PetscCall(VecScale(y, Na->scale)); 140 PetscFunctionReturn(0); 141 } 142 143 PetscErrorCode MatMultAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3) { 144 Mat_Normal *Na = (Mat_Normal *)N->data; 145 Vec in; 146 147 PetscFunctionBegin; 148 in = v1; 149 if (Na->right) { 150 if (!Na->rightwork) PetscCall(VecDuplicate(Na->right, &Na->rightwork)); 151 PetscCall(VecPointwiseMult(Na->rightwork, Na->right, in)); 152 in = Na->rightwork; 153 } 154 PetscCall(MatMult(Na->A, in, Na->w)); 155 PetscCall(VecScale(Na->w, Na->scale)); 156 if (Na->left) { 157 PetscCall(MatMultTranspose(Na->A, Na->w, v3)); 158 PetscCall(VecPointwiseMult(v3, Na->left, v3)); 159 PetscCall(VecAXPY(v3, 1.0, v2)); 160 } else { 161 PetscCall(MatMultTransposeAdd(Na->A, Na->w, v2, v3)); 162 } 163 PetscFunctionReturn(0); 164 } 165 166 PetscErrorCode MatMultTranspose_Normal(Mat N, Vec x, Vec y) { 167 Mat_Normal *Na = (Mat_Normal *)N->data; 168 Vec in; 169 170 PetscFunctionBegin; 171 in = x; 172 if (Na->left) { 173 if (!Na->leftwork) PetscCall(VecDuplicate(Na->left, &Na->leftwork)); 174 PetscCall(VecPointwiseMult(Na->leftwork, Na->left, in)); 175 in = Na->leftwork; 176 } 177 PetscCall(MatMult(Na->A, in, Na->w)); 178 PetscCall(MatMultTranspose(Na->A, Na->w, y)); 179 if (Na->right) PetscCall(VecPointwiseMult(y, Na->right, y)); 180 PetscCall(VecScale(y, Na->scale)); 181 PetscFunctionReturn(0); 182 } 183 184 PetscErrorCode MatMultTransposeAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3) { 185 Mat_Normal *Na = (Mat_Normal *)N->data; 186 Vec in; 187 188 PetscFunctionBegin; 189 in = v1; 190 if (Na->left) { 191 if (!Na->leftwork) PetscCall(VecDuplicate(Na->left, &Na->leftwork)); 192 PetscCall(VecPointwiseMult(Na->leftwork, Na->left, in)); 193 in = Na->leftwork; 194 } 195 PetscCall(MatMult(Na->A, in, Na->w)); 196 PetscCall(VecScale(Na->w, Na->scale)); 197 if (Na->right) { 198 PetscCall(MatMultTranspose(Na->A, Na->w, v3)); 199 PetscCall(VecPointwiseMult(v3, Na->right, v3)); 200 PetscCall(VecAXPY(v3, 1.0, v2)); 201 } else { 202 PetscCall(MatMultTransposeAdd(Na->A, Na->w, v2, v3)); 203 } 204 PetscFunctionReturn(0); 205 } 206 207 PetscErrorCode MatDestroy_Normal(Mat N) { 208 Mat_Normal *Na = (Mat_Normal *)N->data; 209 210 PetscFunctionBegin; 211 PetscCall(MatDestroy(&Na->A)); 212 PetscCall(MatDestroy(&Na->D)); 213 PetscCall(VecDestroy(&Na->w)); 214 PetscCall(VecDestroy(&Na->left)); 215 PetscCall(VecDestroy(&Na->right)); 216 PetscCall(VecDestroy(&Na->leftwork)); 217 PetscCall(VecDestroy(&Na->rightwork)); 218 PetscCall(PetscFree(N->data)); 219 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMat_C", NULL)); 220 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normal_seqaij_C", NULL)); 221 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normal_mpiaij_C", NULL)); 222 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_seqdense_C", NULL)); 223 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_mpidense_C", NULL)); 224 PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_dense_C", NULL)); 225 PetscFunctionReturn(0); 226 } 227 228 /* 229 Slow, nonscalable version 230 */ 231 PetscErrorCode MatGetDiagonal_Normal(Mat N, Vec v) { 232 Mat_Normal *Na = (Mat_Normal *)N->data; 233 Mat A = Na->A; 234 PetscInt i, j, rstart, rend, nnz; 235 const PetscInt *cols; 236 PetscScalar *diag, *work, *values; 237 const PetscScalar *mvalues; 238 239 PetscFunctionBegin; 240 PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work)); 241 PetscCall(PetscArrayzero(work, A->cmap->N)); 242 PetscCall(MatGetOwnershipRange(A, &rstart, &rend)); 243 for (i = rstart; i < rend; i++) { 244 PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues)); 245 for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * mvalues[j]; 246 PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues)); 247 } 248 PetscCall(MPIU_Allreduce(work, diag, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N))); 249 rstart = N->cmap->rstart; 250 rend = N->cmap->rend; 251 PetscCall(VecGetArray(v, &values)); 252 PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart)); 253 PetscCall(VecRestoreArray(v, &values)); 254 PetscCall(PetscFree2(diag, work)); 255 PetscCall(VecScale(v, Na->scale)); 256 PetscFunctionReturn(0); 257 } 258 259 PetscErrorCode MatGetDiagonalBlock_Normal(Mat N, Mat *D) { 260 Mat_Normal *Na = (Mat_Normal *)N->data; 261 Mat M, A = Na->A; 262 263 PetscFunctionBegin; 264 PetscCall(MatGetDiagonalBlock(A, &M)); 265 PetscCall(MatCreateNormal(M, &Na->D)); 266 *D = Na->D; 267 PetscFunctionReturn(0); 268 } 269 270 PetscErrorCode MatNormalGetMat_Normal(Mat A, Mat *M) { 271 Mat_Normal *Aa = (Mat_Normal *)A->data; 272 273 PetscFunctionBegin; 274 *M = Aa->A; 275 PetscFunctionReturn(0); 276 } 277 278 /*@ 279 MatNormalGetMat - Gets the Mat object stored inside a MATNORMAL 280 281 Logically collective on Mat 282 283 Input Parameter: 284 . A - the MATNORMAL matrix 285 286 Output Parameter: 287 . M - the matrix object stored inside A 288 289 Level: intermediate 290 291 .seealso: `MatCreateNormal()` 292 293 @*/ 294 PetscErrorCode MatNormalGetMat(Mat A, Mat *M) { 295 PetscFunctionBegin; 296 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 297 PetscValidType(A, 1); 298 PetscValidPointer(M, 2); 299 PetscUseMethod(A, "MatNormalGetMat_C", (Mat, Mat *), (A, M)); 300 PetscFunctionReturn(0); 301 } 302 303 PetscErrorCode MatConvert_Normal_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) { 304 Mat_Normal *Aa = (Mat_Normal *)A->data; 305 Mat B; 306 PetscInt m, n, M, N; 307 308 PetscFunctionBegin; 309 PetscCall(MatGetSize(A, &M, &N)); 310 PetscCall(MatGetLocalSize(A, &m, &n)); 311 if (reuse == MAT_REUSE_MATRIX) { 312 B = *newmat; 313 PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B)); 314 } else { 315 PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B)); 316 PetscCall(MatProductSetType(B, MATPRODUCT_AtB)); 317 PetscCall(MatProductSetFromOptions(B)); 318 PetscCall(MatProductSymbolic(B)); 319 PetscCall(MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE)); 320 } 321 PetscCall(MatProductNumeric(B)); 322 if (reuse == MAT_INPLACE_MATRIX) { 323 PetscCall(MatHeaderReplace(A, &B)); 324 } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B; 325 PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat)); 326 PetscFunctionReturn(0); 327 } 328 329 typedef struct { 330 Mat work[2]; 331 } Normal_Dense; 332 333 PetscErrorCode MatProductNumeric_Normal_Dense(Mat C) { 334 Mat A, B; 335 Normal_Dense *contents; 336 Mat_Normal *a; 337 PetscScalar *array; 338 339 PetscFunctionBegin; 340 MatCheckProduct(C, 1); 341 A = C->product->A; 342 a = (Mat_Normal *)A->data; 343 B = C->product->B; 344 contents = (Normal_Dense *)C->product->data; 345 PetscCheck(contents, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty"); 346 if (a->right) { 347 PetscCall(MatCopy(B, C, SAME_NONZERO_PATTERN)); 348 PetscCall(MatDiagonalScale(C, a->right, NULL)); 349 } 350 PetscCall(MatProductNumeric(contents->work[0])); 351 PetscCall(MatDenseGetArrayWrite(C, &array)); 352 PetscCall(MatDensePlaceArray(contents->work[1], array)); 353 PetscCall(MatProductNumeric(contents->work[1])); 354 PetscCall(MatDenseRestoreArrayWrite(C, &array)); 355 PetscCall(MatDenseResetArray(contents->work[1])); 356 PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 357 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 358 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 359 PetscCall(MatScale(C, a->scale)); 360 PetscFunctionReturn(0); 361 } 362 363 PetscErrorCode MatNormal_DenseDestroy(void *ctx) { 364 Normal_Dense *contents = (Normal_Dense *)ctx; 365 366 PetscFunctionBegin; 367 PetscCall(MatDestroy(contents->work)); 368 PetscCall(MatDestroy(contents->work + 1)); 369 PetscCall(PetscFree(contents)); 370 PetscFunctionReturn(0); 371 } 372 373 PetscErrorCode MatProductSymbolic_Normal_Dense(Mat C) { 374 Mat A, B; 375 Normal_Dense *contents = NULL; 376 Mat_Normal *a; 377 PetscScalar *array; 378 PetscInt n, N, m, M; 379 380 PetscFunctionBegin; 381 MatCheckProduct(C, 1); 382 PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty"); 383 A = C->product->A; 384 a = (Mat_Normal *)A->data; 385 PetscCheck(!a->left, PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "Not implemented"); 386 B = C->product->B; 387 PetscCall(MatGetLocalSize(C, &m, &n)); 388 PetscCall(MatGetSize(C, &M, &N)); 389 if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) { 390 PetscCall(MatGetLocalSize(B, NULL, &n)); 391 PetscCall(MatGetSize(B, NULL, &N)); 392 PetscCall(MatGetLocalSize(A, &m, NULL)); 393 PetscCall(MatGetSize(A, &M, NULL)); 394 PetscCall(MatSetSizes(C, m, n, M, N)); 395 } 396 PetscCall(MatSetType(C, ((PetscObject)B)->type_name)); 397 PetscCall(MatSetUp(C)); 398 PetscCall(PetscNew(&contents)); 399 C->product->data = contents; 400 C->product->destroy = MatNormal_DenseDestroy; 401 if (a->right) { 402 PetscCall(MatProductCreate(a->A, C, NULL, contents->work)); 403 } else { 404 PetscCall(MatProductCreate(a->A, B, NULL, contents->work)); 405 } 406 PetscCall(MatProductSetType(contents->work[0], MATPRODUCT_AB)); 407 PetscCall(MatProductSetFromOptions(contents->work[0])); 408 PetscCall(MatProductSymbolic(contents->work[0])); 409 PetscCall(MatProductCreate(a->A, contents->work[0], NULL, contents->work + 1)); 410 PetscCall(MatProductSetType(contents->work[1], MATPRODUCT_AtB)); 411 PetscCall(MatProductSetFromOptions(contents->work[1])); 412 PetscCall(MatProductSymbolic(contents->work[1])); 413 PetscCall(MatDenseGetArrayWrite(C, &array)); 414 PetscCall(MatSeqDenseSetPreallocation(contents->work[1], array)); 415 PetscCall(MatMPIDenseSetPreallocation(contents->work[1], array)); 416 PetscCall(MatDenseRestoreArrayWrite(C, &array)); 417 C->ops->productnumeric = MatProductNumeric_Normal_Dense; 418 PetscFunctionReturn(0); 419 } 420 421 PetscErrorCode MatProductSetFromOptions_Normal_Dense_AB(Mat C) { 422 PetscFunctionBegin; 423 C->ops->productsymbolic = MatProductSymbolic_Normal_Dense; 424 PetscFunctionReturn(0); 425 } 426 427 PetscErrorCode MatProductSetFromOptions_Normal_Dense(Mat C) { 428 Mat_Product *product = C->product; 429 430 PetscFunctionBegin; 431 if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_Normal_Dense_AB(C)); 432 PetscFunctionReturn(0); 433 } 434 435 /*@ 436 MatCreateNormal - Creates a new matrix object that behaves like A'*A. 437 438 Collective on Mat 439 440 Input Parameter: 441 . A - the (possibly rectangular) matrix 442 443 Output Parameter: 444 . N - the matrix that represents A'*A 445 446 Level: intermediate 447 448 Notes: 449 The product A'*A is NOT actually formed! Rather the new matrix 450 object performs the matrix-vector product by first multiplying by 451 A and then A' 452 @*/ 453 PetscErrorCode MatCreateNormal(Mat A, Mat *N) { 454 PetscInt n, nn; 455 Mat_Normal *Na; 456 VecType vtype; 457 458 PetscFunctionBegin; 459 PetscCall(MatGetSize(A, NULL, &nn)); 460 PetscCall(MatGetLocalSize(A, NULL, &n)); 461 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N)); 462 PetscCall(MatSetSizes(*N, n, n, nn, nn)); 463 PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMAL)); 464 PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap)); 465 PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap)); 466 467 PetscCall(PetscNewLog(*N, &Na)); 468 (*N)->data = (void *)Na; 469 PetscCall(PetscObjectReference((PetscObject)A)); 470 Na->A = A; 471 Na->scale = 1.0; 472 473 PetscCall(MatCreateVecs(A, NULL, &Na->w)); 474 475 (*N)->ops->destroy = MatDestroy_Normal; 476 (*N)->ops->mult = MatMult_Normal; 477 (*N)->ops->multtranspose = MatMultTranspose_Normal; 478 (*N)->ops->multtransposeadd = MatMultTransposeAdd_Normal; 479 (*N)->ops->multadd = MatMultAdd_Normal; 480 (*N)->ops->getdiagonal = MatGetDiagonal_Normal; 481 (*N)->ops->getdiagonalblock = MatGetDiagonalBlock_Normal; 482 (*N)->ops->scale = MatScale_Normal; 483 (*N)->ops->diagonalscale = MatDiagonalScale_Normal; 484 (*N)->ops->increaseoverlap = MatIncreaseOverlap_Normal; 485 (*N)->ops->createsubmatrices = MatCreateSubMatrices_Normal; 486 (*N)->ops->permute = MatPermute_Normal; 487 (*N)->ops->duplicate = MatDuplicate_Normal; 488 (*N)->ops->copy = MatCopy_Normal; 489 (*N)->assembled = PETSC_TRUE; 490 (*N)->preallocated = PETSC_TRUE; 491 492 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatNormalGetMat_C", MatNormalGetMat_Normal)); 493 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normal_seqaij_C", MatConvert_Normal_AIJ)); 494 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normal_mpiaij_C", MatConvert_Normal_AIJ)); 495 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_seqdense_C", MatProductSetFromOptions_Normal_Dense)); 496 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_mpidense_C", MatProductSetFromOptions_Normal_Dense)); 497 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_dense_C", MatProductSetFromOptions_Normal_Dense)); 498 PetscCall(MatSetOption(*N, MAT_SYMMETRIC, PETSC_TRUE)); 499 PetscCall(MatGetVecType(A, &vtype)); 500 PetscCall(MatSetVecType(*N, vtype)); 501 #if defined(PETSC_HAVE_DEVICE) 502 PetscCall(MatBindToCPU(*N, A->boundtocpu)); 503 #endif 504 PetscFunctionReturn(0); 505 } 506