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