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 A 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: `MATNORMAL`, `MATNORMALHERMITIAN`, `MatCreateNormal()` 292 @*/ 293 PetscErrorCode MatNormalGetMat(Mat A, Mat *M) { 294 PetscFunctionBegin; 295 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 296 PetscValidType(A, 1); 297 PetscValidPointer(M, 2); 298 PetscUseMethod(A, "MatNormalGetMat_C", (Mat, Mat *), (A, M)); 299 PetscFunctionReturn(0); 300 } 301 302 PetscErrorCode MatConvert_Normal_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) { 303 Mat_Normal *Aa = (Mat_Normal *)A->data; 304 Mat B; 305 PetscInt m, n, M, N; 306 307 PetscFunctionBegin; 308 PetscCall(MatGetSize(A, &M, &N)); 309 PetscCall(MatGetLocalSize(A, &m, &n)); 310 if (reuse == MAT_REUSE_MATRIX) { 311 B = *newmat; 312 PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B)); 313 } else { 314 PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B)); 315 PetscCall(MatProductSetType(B, MATPRODUCT_AtB)); 316 PetscCall(MatProductSetFromOptions(B)); 317 PetscCall(MatProductSymbolic(B)); 318 PetscCall(MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE)); 319 } 320 PetscCall(MatProductNumeric(B)); 321 if (reuse == MAT_INPLACE_MATRIX) { 322 PetscCall(MatHeaderReplace(A, &B)); 323 } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B; 324 PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat)); 325 PetscFunctionReturn(0); 326 } 327 328 typedef struct { 329 Mat work[2]; 330 } Normal_Dense; 331 332 PetscErrorCode MatProductNumeric_Normal_Dense(Mat C) { 333 Mat A, B; 334 Normal_Dense *contents; 335 Mat_Normal *a; 336 PetscScalar *array; 337 338 PetscFunctionBegin; 339 MatCheckProduct(C, 1); 340 A = C->product->A; 341 a = (Mat_Normal *)A->data; 342 B = C->product->B; 343 contents = (Normal_Dense *)C->product->data; 344 PetscCheck(contents, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty"); 345 if (a->right) { 346 PetscCall(MatCopy(B, C, SAME_NONZERO_PATTERN)); 347 PetscCall(MatDiagonalScale(C, a->right, NULL)); 348 } 349 PetscCall(MatProductNumeric(contents->work[0])); 350 PetscCall(MatDenseGetArrayWrite(C, &array)); 351 PetscCall(MatDensePlaceArray(contents->work[1], array)); 352 PetscCall(MatProductNumeric(contents->work[1])); 353 PetscCall(MatDenseRestoreArrayWrite(C, &array)); 354 PetscCall(MatDenseResetArray(contents->work[1])); 355 PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE)); 356 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY)); 357 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY)); 358 PetscCall(MatScale(C, a->scale)); 359 PetscFunctionReturn(0); 360 } 361 362 PetscErrorCode MatNormal_DenseDestroy(void *ctx) { 363 Normal_Dense *contents = (Normal_Dense *)ctx; 364 365 PetscFunctionBegin; 366 PetscCall(MatDestroy(contents->work)); 367 PetscCall(MatDestroy(contents->work + 1)); 368 PetscCall(PetscFree(contents)); 369 PetscFunctionReturn(0); 370 } 371 372 PetscErrorCode MatProductSymbolic_Normal_Dense(Mat C) { 373 Mat A, B; 374 Normal_Dense *contents = NULL; 375 Mat_Normal *a; 376 PetscScalar *array; 377 PetscInt n, N, m, M; 378 379 PetscFunctionBegin; 380 MatCheckProduct(C, 1); 381 PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty"); 382 A = C->product->A; 383 a = (Mat_Normal *)A->data; 384 PetscCheck(!a->left, PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "Not implemented"); 385 B = C->product->B; 386 PetscCall(MatGetLocalSize(C, &m, &n)); 387 PetscCall(MatGetSize(C, &M, &N)); 388 if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) { 389 PetscCall(MatGetLocalSize(B, NULL, &n)); 390 PetscCall(MatGetSize(B, NULL, &N)); 391 PetscCall(MatGetLocalSize(A, &m, NULL)); 392 PetscCall(MatGetSize(A, &M, NULL)); 393 PetscCall(MatSetSizes(C, m, n, M, N)); 394 } 395 PetscCall(MatSetType(C, ((PetscObject)B)->type_name)); 396 PetscCall(MatSetUp(C)); 397 PetscCall(PetscNew(&contents)); 398 C->product->data = contents; 399 C->product->destroy = MatNormal_DenseDestroy; 400 if (a->right) { 401 PetscCall(MatProductCreate(a->A, C, NULL, contents->work)); 402 } else { 403 PetscCall(MatProductCreate(a->A, B, NULL, contents->work)); 404 } 405 PetscCall(MatProductSetType(contents->work[0], MATPRODUCT_AB)); 406 PetscCall(MatProductSetFromOptions(contents->work[0])); 407 PetscCall(MatProductSymbolic(contents->work[0])); 408 PetscCall(MatProductCreate(a->A, contents->work[0], NULL, contents->work + 1)); 409 PetscCall(MatProductSetType(contents->work[1], MATPRODUCT_AtB)); 410 PetscCall(MatProductSetFromOptions(contents->work[1])); 411 PetscCall(MatProductSymbolic(contents->work[1])); 412 PetscCall(MatDenseGetArrayWrite(C, &array)); 413 PetscCall(MatSeqDenseSetPreallocation(contents->work[1], array)); 414 PetscCall(MatMPIDenseSetPreallocation(contents->work[1], array)); 415 PetscCall(MatDenseRestoreArrayWrite(C, &array)); 416 C->ops->productnumeric = MatProductNumeric_Normal_Dense; 417 PetscFunctionReturn(0); 418 } 419 420 PetscErrorCode MatProductSetFromOptions_Normal_Dense_AB(Mat C) { 421 PetscFunctionBegin; 422 C->ops->productsymbolic = MatProductSymbolic_Normal_Dense; 423 PetscFunctionReturn(0); 424 } 425 426 PetscErrorCode MatProductSetFromOptions_Normal_Dense(Mat C) { 427 Mat_Product *product = C->product; 428 429 PetscFunctionBegin; 430 if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_Normal_Dense_AB(C)); 431 PetscFunctionReturn(0); 432 } 433 434 /*@ 435 MatCreateNormal - Creates a new `MATNORMAL` matrix object that behaves like A'*A. 436 437 Collective on mat 438 439 Input Parameter: 440 . A - the (possibly rectangular) matrix 441 442 Output Parameter: 443 . N - the matrix that represents A'*A 444 445 Level: intermediate 446 447 Notes: 448 The product A'*A is NOT actually formed! Rather the new matrix 449 object performs the matrix-vector product, `MatMult()`, by first multiplying by 450 A and then A' 451 452 .seealso: `MATNORMAL`, `MatMult()`, `MatNormalGetMat()`, `MATNORMALHERMITIAN`, 453 @*/ 454 PetscErrorCode MatCreateNormal(Mat A, Mat *N) { 455 PetscInt n, nn; 456 Mat_Normal *Na; 457 VecType vtype; 458 459 PetscFunctionBegin; 460 PetscCall(MatGetSize(A, NULL, &nn)); 461 PetscCall(MatGetLocalSize(A, NULL, &n)); 462 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N)); 463 PetscCall(MatSetSizes(*N, n, n, nn, nn)); 464 PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMAL)); 465 PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap)); 466 PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap)); 467 468 PetscCall(PetscNew(&Na)); 469 (*N)->data = (void *)Na; 470 PetscCall(PetscObjectReference((PetscObject)A)); 471 Na->A = A; 472 Na->scale = 1.0; 473 474 PetscCall(MatCreateVecs(A, NULL, &Na->w)); 475 476 (*N)->ops->destroy = MatDestroy_Normal; 477 (*N)->ops->mult = MatMult_Normal; 478 (*N)->ops->multtranspose = MatMultTranspose_Normal; 479 (*N)->ops->multtransposeadd = MatMultTransposeAdd_Normal; 480 (*N)->ops->multadd = MatMultAdd_Normal; 481 (*N)->ops->getdiagonal = MatGetDiagonal_Normal; 482 (*N)->ops->getdiagonalblock = MatGetDiagonalBlock_Normal; 483 (*N)->ops->scale = MatScale_Normal; 484 (*N)->ops->diagonalscale = MatDiagonalScale_Normal; 485 (*N)->ops->increaseoverlap = MatIncreaseOverlap_Normal; 486 (*N)->ops->createsubmatrices = MatCreateSubMatrices_Normal; 487 (*N)->ops->permute = MatPermute_Normal; 488 (*N)->ops->duplicate = MatDuplicate_Normal; 489 (*N)->ops->copy = MatCopy_Normal; 490 (*N)->assembled = PETSC_TRUE; 491 (*N)->preallocated = PETSC_TRUE; 492 493 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatNormalGetMat_C", MatNormalGetMat_Normal)); 494 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normal_seqaij_C", MatConvert_Normal_AIJ)); 495 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normal_mpiaij_C", MatConvert_Normal_AIJ)); 496 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_seqdense_C", MatProductSetFromOptions_Normal_Dense)); 497 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_mpidense_C", MatProductSetFromOptions_Normal_Dense)); 498 PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_dense_C", MatProductSetFromOptions_Normal_Dense)); 499 PetscCall(MatSetOption(*N, MAT_SYMMETRIC, PETSC_TRUE)); 500 PetscCall(MatGetVecType(A, &vtype)); 501 PetscCall(MatSetVecType(*N, vtype)); 502 #if defined(PETSC_HAVE_DEVICE) 503 PetscCall(MatBindToCPU(*N, A->boundtocpu)); 504 #endif 505 PetscFunctionReturn(0); 506 } 507