1 2 /* 3 Routines for matrix products. Calling procedure: 4 5 MatProductCreate(A,B,C,&D); or MatProductCreateWithMat(A,B,C,D); 6 MatProductSetType(D, MATPRODUCT_AB/AtB/ABt/PtAP/RARt/ABC); 7 MatProductSetAlgorithm(D, alg); 8 MatProductSetFill(D,fill); 9 MatProductSetFromOptions(D); 10 -> MatProductSetFromOptions_producttype(D): 11 # Check matrix global sizes 12 -> MatProductSetFromOptions_Atype_Btype_Ctype(D); 13 ->MatProductSetFromOptions_Atype_Btype_Ctype_productype(D): 14 # Check matrix local sizes for mpi matrices 15 # Set default algorithm 16 # Get runtime option 17 # Set D->ops->productsymbolic = MatProductSymbolic_productype_Atype_Btype_Ctype; 18 19 PetscLogEventBegin() 20 MatProductSymbolic(D): 21 # Call MatxxxSymbolic_Atype_Btype_Ctype(); 22 # Set D->ops->productnumeric = MatProductNumeric_productype_Atype_Btype_Ctype; 23 PetscLogEventEnd() 24 25 PetscLogEventBegin() 26 MatProductNumeric(D); 27 # Call (D->ops->matxxxnumeric)(); 28 PetscLogEventEnd() 29 */ 30 31 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 32 33 static PetscErrorCode MatProductNumeric_PtAP_Basic(Mat C) 34 { 35 PetscErrorCode ierr; 36 Mat_Product *product = C->product; 37 Mat P = product->B,AP = product->Dwork; 38 39 PetscFunctionBegin; 40 /* AP = A*P */ 41 ierr = MatProductNumeric(AP);CHKERRQ(ierr); 42 /* C = P^T*AP */ 43 ierr = (C->ops->transposematmultnumeric)(P,AP,C);CHKERRQ(ierr); 44 PetscFunctionReturn(0); 45 } 46 47 static PetscErrorCode MatProductSymbolic_PtAP_Basic(Mat C) 48 { 49 PetscErrorCode ierr; 50 Mat_Product *product = C->product; 51 Mat A=product->A,P=product->B,AP; 52 PetscReal fill=product->fill; 53 54 PetscFunctionBegin; 55 /* AP = A*P */ 56 ierr = MatProductCreate(A,P,NULL,&AP);CHKERRQ(ierr); 57 ierr = MatProductSetType(AP,MATPRODUCT_AB);CHKERRQ(ierr); 58 ierr = MatProductSetAlgorithm(AP,"default");CHKERRQ(ierr); 59 ierr = MatProductSetFill(AP,fill);CHKERRQ(ierr); 60 ierr = MatProductSetFromOptions(AP);CHKERRQ(ierr); 61 ierr = MatProductSymbolic(AP);CHKERRQ(ierr); 62 63 /* C = P^T*AP */ 64 ierr = MatProductSetType(C,MATPRODUCT_AtB);CHKERRQ(ierr); 65 product->alg = "default"; 66 product->A = P; 67 product->B = AP; 68 ierr = MatProductSetFromOptions(C);CHKERRQ(ierr); 69 ierr = MatProductSymbolic(C);CHKERRQ(ierr); 70 71 /* resume user's original input matrix setting for A and B */ 72 product->A = A; 73 product->B = P; 74 product->Dwork = AP; 75 76 C->ops->productnumeric = MatProductNumeric_PtAP_Basic; 77 PetscFunctionReturn(0); 78 } 79 80 static PetscErrorCode MatProductNumeric_RARt_Basic(Mat C) 81 { 82 PetscErrorCode ierr; 83 Mat_Product *product = C->product; 84 Mat R=product->B,RA=product->Dwork; 85 86 PetscFunctionBegin; 87 /* RA = R*A */ 88 ierr = MatProductNumeric(RA);CHKERRQ(ierr); 89 /* C = RA*R^T */ 90 ierr = (C->ops->mattransposemultnumeric)(RA,R,C);CHKERRQ(ierr); 91 PetscFunctionReturn(0); 92 } 93 94 static PetscErrorCode MatProductSymbolic_RARt_Basic(Mat C) 95 { 96 PetscErrorCode ierr; 97 Mat_Product *product = C->product; 98 Mat A=product->A,R=product->B,RA; 99 PetscReal fill=product->fill; 100 101 PetscFunctionBegin; 102 /* RA = R*A */ 103 ierr = MatProductCreate(R,A,NULL,&RA);CHKERRQ(ierr); 104 ierr = MatProductSetType(RA,MATPRODUCT_AB);CHKERRQ(ierr); 105 ierr = MatProductSetAlgorithm(RA,"default");CHKERRQ(ierr); 106 ierr = MatProductSetFill(RA,fill);CHKERRQ(ierr); 107 ierr = MatProductSetFromOptions(RA);CHKERRQ(ierr); 108 ierr = MatProductSymbolic(RA);CHKERRQ(ierr); 109 110 /* C = RA*R^T */ 111 ierr = MatProductSetType(C,MATPRODUCT_ABt);CHKERRQ(ierr); 112 product->alg = "default"; 113 product->A = RA; 114 ierr = MatProductSetFromOptions(C);CHKERRQ(ierr); 115 ierr = MatProductSymbolic(C);CHKERRQ(ierr); 116 117 /* resume user's original input matrix setting for A */ 118 product->A = A; 119 product->Dwork = RA; /* save here so it will be destroyed with product C */ 120 C->ops->productnumeric = MatProductNumeric_RARt_Basic; 121 PetscFunctionReturn(0); 122 } 123 124 static PetscErrorCode MatProductNumeric_ABC_Basic(Mat mat) 125 { 126 PetscErrorCode ierr; 127 Mat_Product *product = mat->product; 128 Mat A=product->A,BC=product->Dwork; 129 130 PetscFunctionBegin; 131 /* Numeric BC = B*C */ 132 ierr = MatProductNumeric(BC);CHKERRQ(ierr); 133 /* Numeric mat = A*BC */ 134 ierr = (mat->ops->matmultnumeric)(A,BC,mat);CHKERRQ(ierr); 135 PetscFunctionReturn(0); 136 } 137 138 static PetscErrorCode MatProductSymbolic_ABC_Basic(Mat mat) 139 { 140 PetscErrorCode ierr; 141 Mat_Product *product = mat->product; 142 Mat B=product->B,C=product->C,BC; 143 PetscReal fill=product->fill; 144 145 PetscFunctionBegin; 146 /* Symbolic BC = B*C */ 147 ierr = MatProductCreate(B,C,NULL,&BC);CHKERRQ(ierr); 148 ierr = MatProductSetType(BC,MATPRODUCT_AB);CHKERRQ(ierr); 149 ierr = MatProductSetAlgorithm(BC,"default");CHKERRQ(ierr); 150 ierr = MatProductSetFill(BC,fill);CHKERRQ(ierr); 151 ierr = MatProductSetFromOptions(BC);CHKERRQ(ierr); 152 ierr = MatProductSymbolic(BC);CHKERRQ(ierr); 153 154 /* Symbolic mat = A*BC */ 155 ierr = MatProductSetType(mat,MATPRODUCT_AB);CHKERRQ(ierr); 156 product->alg = "default"; 157 product->B = BC; 158 product->Dwork = BC; 159 ierr = MatProductSetFromOptions(mat);CHKERRQ(ierr); 160 ierr = MatProductSymbolic(mat);CHKERRQ(ierr); 161 162 /* resume user's original input matrix setting for B */ 163 product->B = B; 164 mat->ops->productnumeric = MatProductNumeric_ABC_Basic; 165 PetscFunctionReturn(0); 166 } 167 168 PetscErrorCode MatProductSymbolic_Basic(Mat mat) 169 { 170 PetscErrorCode ierr; 171 Mat_Product *product = mat->product; 172 173 PetscFunctionBegin; 174 switch (product->type) { 175 case MATPRODUCT_PtAP: 176 PetscInfo2((PetscObject)mat, "MatProduct_Basic_PtAP() for A %s, P %s is used",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name); 177 ierr = MatProductSymbolic_PtAP_Basic(mat);CHKERRQ(ierr); 178 break; 179 case MATPRODUCT_RARt: 180 PetscInfo2((PetscObject)mat, "MatProduct_Basic_RARt() for A %s, R %s is used",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name); 181 ierr = MatProductSymbolic_RARt_Basic(mat);CHKERRQ(ierr); 182 break; 183 case MATPRODUCT_ABC: 184 PetscInfo3((PetscObject)mat, "MatProduct_Basic_ABC() for A %s, B %s, C %s is used",((PetscObject)product->A)->type_name,((PetscObject)product->B)->type_name,((PetscObject)product->C)->type_name); 185 ierr = MatProductSymbolic_ABC_Basic(mat);CHKERRQ(ierr); 186 break; 187 default: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType is not supported"); 188 } 189 PetscFunctionReturn(0); 190 } 191 192 /* ----------------------------------------------- */ 193 /*@C 194 MatProductReplaceMats - Replace input matrices for a matrix product. 195 196 Collective on Mat 197 198 Input Parameters: 199 + A - the matrix or NULL if not being replaced 200 . B - the matrix or NULL if not being replaced 201 . C - the matrix or NULL if not being replaced 202 - D - the matrix product 203 204 Level: intermediate 205 206 Notes: 207 Input matrix must have exactly same data structure as replaced one. 208 209 .seealso: MatProductCreate() 210 @*/ 211 PetscErrorCode MatProductReplaceMats(Mat A,Mat B,Mat C,Mat D) 212 { 213 Mat_Product *product=D->product; 214 215 PetscFunctionBegin; 216 if (!product) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_NULL,"Mat D does not have struct 'product'. Call MatProductReplaceProduct(). \n"); 217 if (A) { 218 if (!product->Areplaced) { 219 product->A = A; 220 } else SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"Matrix A was changed by a PETSc internal routine, cannot be replaced"); 221 } 222 if (B) { 223 if (!product->Breplaced) { 224 product->B = B; 225 } else SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"Matrix B was changed by a PETSc internal routine, cannot be replaced"); 226 } 227 if (C) product->C = C; 228 PetscFunctionReturn(0); 229 } 230 231 /* ----------------------------------------------- */ 232 static PetscErrorCode MatProductSetFromOptions_AB(Mat mat) 233 { 234 PetscErrorCode ierr; 235 Mat_Product *product = mat->product; 236 Mat A=product->A,B=product->B; 237 PetscBool sametype; 238 PetscErrorCode (*fA)(Mat); 239 PetscErrorCode (*fB)(Mat); 240 PetscErrorCode (*f)(Mat)=NULL; 241 PetscBool A_istrans,B_istrans; 242 243 PetscFunctionBegin; 244 /* Check matrix global sizes */ 245 if (B->rmap->N!=A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->rmap->N,A->cmap->N); 246 247 fA = A->ops->productsetfromoptions; 248 fB = B->ops->productsetfromoptions; 249 250 ierr = PetscStrcmp(((PetscObject)A)->type_name,((PetscObject)B)->type_name,&sametype);CHKERRQ(ierr); 251 ierr = PetscObjectTypeCompare((PetscObject)A,MATTRANSPOSEMAT,&A_istrans);CHKERRQ(ierr); 252 ierr = PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&B_istrans);CHKERRQ(ierr); 253 254 if (fB == fA && sametype && (!A_istrans || !B_istrans)) { 255 f = fB; 256 } else { 257 char mtypes[256]; 258 PetscBool At_istrans=PETSC_TRUE,Bt_istrans=PETSC_TRUE; 259 Mat At = NULL,Bt = NULL; 260 261 if (A_istrans && !B_istrans) { 262 ierr = MatTransposeGetMat(A,&At);CHKERRQ(ierr); 263 ierr = PetscObjectTypeCompare((PetscObject)At,MATTRANSPOSEMAT,&At_istrans);CHKERRQ(ierr); 264 if (At_istrans) { /* mat = ATT * B */ 265 Mat Att = NULL; 266 ierr = MatTransposeGetMat(At,&Att);CHKERRQ(ierr); 267 A = Att; 268 product->A = Att; /* use Att for matproduct */ 269 product->Areplaced = PETSC_TRUE; /* Att = A, but has native matrix type */ 270 } else { /* !At_istrans: mat = At^T*B */ 271 A = At; 272 product->A = At; 273 product->Areplaced = PETSC_TRUE; 274 product->type = MATPRODUCT_AtB; 275 } 276 } else if (!A_istrans && B_istrans) { 277 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr); 278 ierr = PetscObjectTypeCompare((PetscObject)Bt,MATTRANSPOSEMAT,&Bt_istrans);CHKERRQ(ierr); 279 if (Bt_istrans) { /* mat = A * BTT */ 280 Mat Btt = NULL; 281 ierr = MatTransposeGetMat(Bt,&Btt);CHKERRQ(ierr); 282 B = Btt; 283 product->B = Btt; /* use Btt for matproduct */ 284 product->Breplaced = PETSC_TRUE; 285 } else { /* !Bt_istrans */ 286 /* mat = A*Bt^T */ 287 B = Bt; 288 product->B = Bt; 289 product->Breplaced = PETSC_TRUE; 290 product->type = MATPRODUCT_ABt; 291 } 292 } else if (A_istrans && B_istrans) { /* mat = At^T * Bt^T */ 293 ierr = MatTransposeGetMat(A,&At);CHKERRQ(ierr); 294 ierr = PetscObjectTypeCompare((PetscObject)At,MATTRANSPOSEMAT,&At_istrans);CHKERRQ(ierr); 295 ierr = MatTransposeGetMat(B,&Bt);CHKERRQ(ierr); 296 ierr = PetscObjectTypeCompare((PetscObject)Bt,MATTRANSPOSEMAT,&Bt_istrans);CHKERRQ(ierr); 297 if (At_istrans && Bt_istrans) { 298 Mat Att= NULL,Btt = NULL; 299 ierr = MatTransposeGetMat(At,&Att);CHKERRQ(ierr); 300 ierr = MatTransposeGetMat(Bt,&Btt);CHKERRQ(ierr); 301 A = Att; 302 product->A = Att; product->Areplaced = PETSC_TRUE; 303 B = Btt; 304 product->B = Btt; product->Breplaced = PETSC_TRUE; 305 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Not supported yet"); 306 } 307 308 /* query MatProductSetFromOptions_Atype_Btype */ 309 ierr = PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));CHKERRQ(ierr); 310 ierr = PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));CHKERRQ(ierr); 311 ierr = PetscStrlcat(mtypes,"_",sizeof(mtypes));CHKERRQ(ierr); 312 ierr = PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));CHKERRQ(ierr); 313 ierr = PetscStrlcat(mtypes,"_C",sizeof(mtypes));CHKERRQ(ierr); 314 ierr = PetscObjectQueryFunction((PetscObject)A,mtypes,&f);CHKERRQ(ierr); 315 if (!f) { 316 ierr = PetscObjectQueryFunction((PetscObject)B,mtypes,&f);CHKERRQ(ierr); 317 } 318 } 319 320 if (!f) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"MatProductSetFromOptions_AB for A %s and B %s is not supported",((PetscObject)A)->type_name,((PetscObject)B)->type_name); 321 ierr = (*f)(mat);CHKERRQ(ierr); 322 PetscFunctionReturn(0); 323 } 324 325 static PetscErrorCode MatProductSetFromOptions_AtB(Mat mat) 326 { 327 PetscErrorCode ierr; 328 Mat_Product *product = mat->product; 329 Mat A=product->A,B=product->B; 330 PetscBool sametype; 331 PetscErrorCode (*fA)(Mat); 332 PetscErrorCode (*fB)(Mat); 333 PetscErrorCode (*f)(Mat)=NULL; 334 335 PetscFunctionBegin; 336 /* Check matrix global sizes */ 337 if (B->rmap->N!=A->rmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->rmap->N,A->rmap->N); 338 339 fA = A->ops->productsetfromoptions; 340 fB = B->ops->productsetfromoptions; 341 342 ierr = PetscStrcmp(((PetscObject)A)->type_name,((PetscObject)B)->type_name,&sametype);CHKERRQ(ierr); 343 344 if (fB == fA && sametype) { 345 f = fB; 346 } else { 347 char mtypes[256]; 348 PetscBool istrans; 349 ierr = PetscObjectTypeCompare((PetscObject)A,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr); 350 if (!istrans) { 351 /* query MatProductSetFromOptions_Atype_Btype */ 352 ierr = PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));CHKERRQ(ierr); 353 ierr = PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));CHKERRQ(ierr); 354 ierr = PetscStrlcat(mtypes,"_",sizeof(mtypes));CHKERRQ(ierr); 355 ierr = PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));CHKERRQ(ierr); 356 ierr = PetscStrlcat(mtypes,"_C",sizeof(mtypes));CHKERRQ(ierr); 357 ierr = PetscObjectQueryFunction((PetscObject)B,mtypes,&f);CHKERRQ(ierr); 358 } else { 359 Mat T = NULL; 360 SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"MatProductSetFromOptions_AtB for A %s and B %s is not supported",((PetscObject)A)->type_name,((PetscObject)B)->type_name); 361 362 ierr = MatTransposeGetMat(A,&T);CHKERRQ(ierr); 363 ierr = PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));CHKERRQ(ierr); 364 ierr = PetscStrlcat(mtypes,((PetscObject)T)->type_name,sizeof(mtypes));CHKERRQ(ierr); 365 ierr = PetscStrlcat(mtypes,"_",sizeof(mtypes));CHKERRQ(ierr); 366 ierr = PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));CHKERRQ(ierr); 367 ierr = PetscStrlcat(mtypes,"_C",sizeof(mtypes));CHKERRQ(ierr); 368 369 product->type = MATPRODUCT_AtB; 370 ierr = PetscObjectQueryFunction((PetscObject)B,mtypes,&f);CHKERRQ(ierr); 371 } 372 373 if (!f) { 374 ierr = PetscObjectQueryFunction((PetscObject)A,mtypes,&f);CHKERRQ(ierr); 375 } 376 } 377 if (!f) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"MatProductSetFromOptions_AB for A %s and B %s is not supported",((PetscObject)A)->type_name,((PetscObject)B)->type_name); 378 379 ierr = (*f)(mat);CHKERRQ(ierr); 380 PetscFunctionReturn(0); 381 } 382 383 static PetscErrorCode MatProductSetFromOptions_ABt(Mat mat) 384 { 385 PetscErrorCode ierr; 386 Mat_Product *product = mat->product; 387 Mat A=product->A,B=product->B; 388 PetscBool sametype; 389 PetscErrorCode (*fA)(Mat); 390 PetscErrorCode (*fB)(Mat); 391 PetscErrorCode (*f)(Mat)=NULL; 392 393 PetscFunctionBegin; 394 /* Check matrix global sizes */ 395 if (B->cmap->N!=A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, AN %D != BN %D",A->cmap->N,B->cmap->N); 396 397 fA = A->ops->productsetfromoptions; 398 fB = B->ops->productsetfromoptions; 399 400 ierr = PetscStrcmp(((PetscObject)A)->type_name,((PetscObject)B)->type_name,&sametype);CHKERRQ(ierr); 401 402 if (fB == fA && sametype) { 403 f = fB; 404 } else { 405 char mtypes[256]; 406 PetscBool istrans; 407 ierr = PetscObjectTypeCompare((PetscObject)A,MATTRANSPOSEMAT,&istrans);CHKERRQ(ierr); 408 if (!istrans) { 409 /* query MatProductSetFromOptions_Atype_Btype */ 410 ierr = PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));CHKERRQ(ierr); 411 ierr = PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));CHKERRQ(ierr); 412 ierr = PetscStrlcat(mtypes,"_",sizeof(mtypes));CHKERRQ(ierr); 413 ierr = PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));CHKERRQ(ierr); 414 ierr = PetscStrlcat(mtypes,"_C",sizeof(mtypes));CHKERRQ(ierr); 415 ierr = PetscObjectQueryFunction((PetscObject)B,mtypes,&f);CHKERRQ(ierr); 416 } else { 417 Mat T = NULL; 418 SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"MatProductSetFromOptions_ABt for A %s and B %s is not supported",((PetscObject)A)->type_name,((PetscObject)B)->type_name); 419 420 ierr = MatTransposeGetMat(A,&T);CHKERRQ(ierr); 421 ierr = PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));CHKERRQ(ierr); 422 ierr = PetscStrlcat(mtypes,((PetscObject)T)->type_name,sizeof(mtypes));CHKERRQ(ierr); 423 ierr = PetscStrlcat(mtypes,"_",sizeof(mtypes));CHKERRQ(ierr); 424 ierr = PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));CHKERRQ(ierr); 425 ierr = PetscStrlcat(mtypes,"_C",sizeof(mtypes));CHKERRQ(ierr); 426 427 product->type = MATPRODUCT_ABt; 428 ierr = PetscObjectQueryFunction((PetscObject)B,mtypes,&f);CHKERRQ(ierr); 429 } 430 431 if (!f) { 432 ierr = PetscObjectQueryFunction((PetscObject)A,mtypes,&f);CHKERRQ(ierr); 433 } 434 } 435 if (!f) { 436 SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"MatProductSetFromOptions_AB for A %s and B %s is not supported",((PetscObject)A)->type_name,((PetscObject)B)->type_name); 437 } 438 439 ierr = (*f)(mat);CHKERRQ(ierr); 440 PetscFunctionReturn(0); 441 } 442 443 static PetscErrorCode MatProductSetFromOptions_PtAP(Mat mat) 444 { 445 PetscErrorCode ierr; 446 Mat_Product *product = mat->product; 447 Mat A=product->A,B=product->B; 448 PetscBool sametype; 449 PetscErrorCode (*fA)(Mat); 450 PetscErrorCode (*fB)(Mat); 451 PetscErrorCode (*f)(Mat)=NULL; 452 453 PetscFunctionBegin; 454 /* Check matrix global sizes */ 455 if (A->rmap->N != A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix A must be square, %D != %D",A->rmap->N,A->cmap->N); 456 if (B->rmap->N != A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->rmap->N,A->cmap->N); 457 458 fA = A->ops->productsetfromoptions; 459 fB = B->ops->productsetfromoptions; 460 461 ierr = PetscStrcmp(((PetscObject)A)->type_name,((PetscObject)B)->type_name,&sametype);CHKERRQ(ierr); 462 if (fB == fA && sametype) { 463 f = fB; 464 } else { 465 /* query MatProductSetFromOptions_Atype_Btype */ 466 char mtypes[256]; 467 ierr = PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));CHKERRQ(ierr); 468 ierr = PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));CHKERRQ(ierr); 469 ierr = PetscStrlcat(mtypes,"_",sizeof(mtypes));CHKERRQ(ierr); 470 ierr = PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));CHKERRQ(ierr); 471 ierr = PetscStrlcat(mtypes,"_C",sizeof(mtypes));CHKERRQ(ierr); 472 ierr = PetscObjectQueryFunction((PetscObject)B,mtypes,&f);CHKERRQ(ierr); 473 474 if (!f) { 475 ierr = PetscObjectQueryFunction((PetscObject)A,mtypes,&f);CHKERRQ(ierr); 476 } 477 } 478 479 if (f) { 480 ierr = (*f)(mat);CHKERRQ(ierr); 481 } else { 482 mat->ops->productsymbolic = MatProductSymbolic_Basic; 483 PetscInfo2((PetscObject)mat, "MatProductSetFromOptions_PtAP for A %s, P %s uses MatProduct_Basic() implementation",((PetscObject)A)->type_name,((PetscObject)B)->type_name); 484 } 485 PetscFunctionReturn(0); 486 } 487 488 static PetscErrorCode MatProductSetFromOptions_RARt(Mat mat) 489 { 490 PetscErrorCode ierr; 491 Mat_Product *product = mat->product; 492 Mat A=product->A,B=product->B; 493 PetscBool sametype; 494 PetscErrorCode (*fA)(Mat); 495 PetscErrorCode (*fB)(Mat); 496 PetscErrorCode (*f)(Mat)=NULL; 497 498 PetscFunctionBegin; 499 /* Check matrix global sizes */ 500 if (A->rmap->N != B->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix A must be square, %D != %D",A->rmap->N,A->cmap->N); 501 if (B->cmap->N != A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->cmap->N,A->cmap->N); 502 503 fA = A->ops->productsetfromoptions; 504 fB = B->ops->productsetfromoptions; 505 506 ierr = PetscStrcmp(((PetscObject)A)->type_name,((PetscObject)B)->type_name,&sametype);CHKERRQ(ierr); 507 if (fB == fA && sametype) { 508 f = fB; 509 } else { 510 /* query MatProductSetFromOptions_Atype_Btype */ 511 char mtypes[256]; 512 ierr = PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));CHKERRQ(ierr); 513 ierr = PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));CHKERRQ(ierr); 514 ierr = PetscStrlcat(mtypes,"_",sizeof(mtypes));CHKERRQ(ierr); 515 ierr = PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));CHKERRQ(ierr); 516 ierr = PetscStrlcat(mtypes,"_C",sizeof(mtypes));CHKERRQ(ierr); 517 ierr = PetscObjectQueryFunction((PetscObject)B,mtypes,&f);CHKERRQ(ierr); 518 519 if (!f) { 520 ierr = PetscObjectQueryFunction((PetscObject)A,mtypes,&f);CHKERRQ(ierr); 521 } 522 } 523 524 if (f) { 525 ierr = (*f)(mat);CHKERRQ(ierr); 526 } else { 527 mat->ops->productsymbolic = MatProductSymbolic_Basic; 528 } 529 PetscFunctionReturn(0); 530 } 531 532 static PetscErrorCode MatProductSetFromOptions_ABC(Mat mat) 533 { 534 PetscErrorCode ierr; 535 Mat_Product *product = mat->product; 536 Mat A=product->A,B=product->B,C=product->C; 537 PetscErrorCode (*fA)(Mat); 538 PetscErrorCode (*fB)(Mat); 539 PetscErrorCode (*fC)(Mat); 540 PetscErrorCode (*f)(Mat)=NULL; 541 542 PetscFunctionBegin; 543 /* Check matrix global sizes */ 544 if (B->rmap->N!= A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",B->rmap->N,A->cmap->N); 545 if (C->rmap->N!= B->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_SIZ,"Matrix dimensions are incompatible, %D != %D",C->rmap->N,B->cmap->N); 546 547 fA = A->ops->productsetfromoptions; 548 fB = B->ops->productsetfromoptions; 549 fC = C->ops->productsetfromoptions; 550 if (fA == fB && fA == fC && fA) { 551 f = fA; 552 } else { 553 /* query MatProductSetFromOptions_Atype_Btype_Ctype */ 554 char mtypes[256]; 555 ierr = PetscStrncpy(mtypes,"MatProductSetFromOptions_",sizeof(mtypes));CHKERRQ(ierr); 556 ierr = PetscStrlcat(mtypes,((PetscObject)A)->type_name,sizeof(mtypes));CHKERRQ(ierr); 557 ierr = PetscStrlcat(mtypes,"_",sizeof(mtypes));CHKERRQ(ierr); 558 ierr = PetscStrlcat(mtypes,((PetscObject)B)->type_name,sizeof(mtypes));CHKERRQ(ierr); 559 ierr = PetscStrlcat(mtypes,"_",sizeof(mtypes));CHKERRQ(ierr); 560 ierr = PetscStrlcat(mtypes,((PetscObject)C)->type_name,sizeof(mtypes));CHKERRQ(ierr); 561 ierr = PetscStrlcat(mtypes,"_C",sizeof(mtypes));CHKERRQ(ierr); 562 563 ierr = PetscObjectQueryFunction((PetscObject)A,mtypes,&f);CHKERRQ(ierr); 564 if (!f) { 565 ierr = PetscObjectQueryFunction((PetscObject)B,mtypes,&f);CHKERRQ(ierr); 566 } 567 if (!f) { 568 ierr = PetscObjectQueryFunction((PetscObject)C,mtypes,&f);CHKERRQ(ierr); 569 } 570 } 571 572 if (f) { 573 ierr = (*f)(mat);CHKERRQ(ierr); 574 } else { /* use MatProductSymbolic/Numeric_Basic() */ 575 mat->ops->productsymbolic = MatProductSymbolic_Basic; 576 } 577 PetscFunctionReturn(0); 578 } 579 580 /*@C 581 MatProductSetFromOptions - Creates a matrix product where the type, the algorithm etc are determined from the options database. 582 583 Logically Collective on Mat 584 585 Input Parameter: 586 . mat - the matrix 587 588 Level: beginner 589 590 .seealso: MatSetFromOptions() 591 @*/ 592 PetscErrorCode MatProductSetFromOptions(Mat mat) 593 { 594 PetscErrorCode ierr; 595 596 PetscFunctionBegin; 597 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 598 599 if (mat->ops->productsetfromoptions) { 600 ierr = (*mat->ops->productsetfromoptions)(mat);CHKERRQ(ierr); 601 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_NULL,"Call MatProductSetType() first"); 602 PetscFunctionReturn(0); 603 } 604 605 /* ----------------------------------------------- */ 606 PetscErrorCode MatProductNumeric_AB(Mat mat) 607 { 608 PetscErrorCode ierr; 609 Mat_Product *product = mat->product; 610 Mat A=product->A,B=product->B; 611 612 PetscFunctionBegin; 613 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 614 ierr = (mat->ops->matmultnumeric)(A,B,mat);CHKERRQ(ierr); 615 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 616 PetscFunctionReturn(0); 617 } 618 619 PetscErrorCode MatProductNumeric_AtB(Mat mat) 620 { 621 PetscErrorCode ierr; 622 Mat_Product *product = mat->product; 623 Mat A=product->A,B=product->B; 624 625 PetscFunctionBegin; 626 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 627 ierr = (mat->ops->transposematmultnumeric)(A,B,mat);CHKERRQ(ierr); 628 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 629 PetscFunctionReturn(0); 630 } 631 632 PetscErrorCode MatProductNumeric_ABt(Mat mat) 633 { 634 PetscErrorCode ierr; 635 Mat_Product *product = mat->product; 636 Mat A=product->A,B=product->B; 637 638 PetscFunctionBegin; 639 ierr = PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 640 ierr = (mat->ops->mattransposemultnumeric)(A,B,mat);CHKERRQ(ierr); 641 ierr = PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 642 PetscFunctionReturn(0); 643 } 644 645 PetscErrorCode MatProductNumeric_PtAP(Mat mat) 646 { 647 PetscErrorCode ierr; 648 Mat_Product *product = mat->product; 649 Mat A=product->A,B=product->B; 650 651 PetscFunctionBegin; 652 ierr = PetscLogEventBegin(MAT_PtAPNumeric,mat,0,0,0);CHKERRQ(ierr); 653 ierr = (mat->ops->ptapnumeric)(A,B,mat);CHKERRQ(ierr); 654 ierr = PetscLogEventEnd(MAT_PtAPNumeric,mat,0,0,0);CHKERRQ(ierr); 655 PetscFunctionReturn(0); 656 } 657 658 PetscErrorCode MatProductNumeric_RARt(Mat mat) 659 { 660 PetscErrorCode ierr; 661 Mat_Product *product = mat->product; 662 Mat A=product->A,B=product->B; 663 664 PetscFunctionBegin; 665 ierr = PetscLogEventBegin(MAT_RARtNumeric,A,B,0,0);CHKERRQ(ierr); 666 ierr = (mat->ops->rartnumeric)(A,B,mat);CHKERRQ(ierr); 667 ierr = PetscLogEventEnd(MAT_RARtNumeric,A,B,0,0);CHKERRQ(ierr); 668 PetscFunctionReturn(0); 669 } 670 671 PetscErrorCode MatProductNumeric_ABC(Mat mat) 672 { 673 PetscErrorCode ierr; 674 Mat_Product *product = mat->product; 675 Mat A=product->A,B=product->B,C=product->C; 676 677 PetscFunctionBegin; 678 ierr = PetscLogEventBegin(MAT_MatMatMultNumeric,A,B,C,0);CHKERRQ(ierr); 679 ierr = (mat->ops->matmatmultnumeric)(A,B,C,mat);CHKERRQ(ierr); 680 ierr = PetscLogEventEnd(MAT_MatMatMultNumeric,A,B,C,0);CHKERRQ(ierr); 681 PetscFunctionReturn(0); 682 } 683 684 /*@ 685 MatProductNumeric - Implement a matrix product with numerical values. 686 687 Collective on Mat 688 689 Input Parameters: 690 . mat - the matrix to hold a product 691 692 Output Parameters: 693 . mat - the matrix product 694 695 Level: intermediate 696 697 .seealso: MatProductCreate(), MatSetType() 698 @*/ 699 PetscErrorCode MatProductNumeric(Mat mat) 700 { 701 PetscErrorCode ierr; 702 703 PetscFunctionBegin; 704 PetscValidPointer(mat,1); 705 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 706 707 if (mat->ops->productnumeric) { 708 ierr = (*mat->ops->productnumeric)(mat);CHKERRQ(ierr); 709 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_USER,"Call MatProductSymbolic() first"); 710 PetscFunctionReturn(0); 711 } 712 713 /* ----------------------------------------------- */ 714 PetscErrorCode MatProductSymbolic_AB(Mat mat) 715 { 716 PetscErrorCode ierr; 717 Mat_Product *product = mat->product; 718 Mat A=product->A,B=product->B; 719 720 PetscFunctionBegin; 721 ierr = (mat->ops->matmultsymbolic)(A,B,product->fill,mat);CHKERRQ(ierr); 722 mat->ops->productnumeric = MatProductNumeric_AB; 723 PetscFunctionReturn(0); 724 } 725 726 PetscErrorCode MatProductSymbolic_AtB(Mat mat) 727 { 728 PetscErrorCode ierr; 729 Mat_Product *product = mat->product; 730 Mat A=product->A,B=product->B; 731 732 PetscFunctionBegin; 733 ierr = (mat->ops->transposematmultsymbolic)(A,B,product->fill,mat);CHKERRQ(ierr); 734 mat->ops->productnumeric = MatProductNumeric_AtB; 735 PetscFunctionReturn(0); 736 } 737 738 PetscErrorCode MatProductSymbolic_ABt(Mat mat) 739 { 740 PetscErrorCode ierr; 741 Mat_Product *product = mat->product; 742 Mat A=product->A,B=product->B; 743 744 PetscFunctionBegin; 745 ierr = (mat->ops->mattransposemultsymbolic)(A,B,product->fill,mat);CHKERRQ(ierr); 746 mat->ops->productnumeric = MatProductNumeric_ABt; 747 PetscFunctionReturn(0); 748 } 749 750 PetscErrorCode MatProductSymbolic_ABC(Mat mat) 751 { 752 PetscErrorCode ierr; 753 Mat_Product *product = mat->product; 754 Mat A=product->A,B=product->B,C=product->C; 755 756 PetscFunctionBegin; 757 ierr = (mat->ops->matmatmultsymbolic)(A,B,C,product->fill,mat);CHKERRQ(ierr); 758 mat->ops->productnumeric = MatProductNumeric_ABC; 759 PetscFunctionReturn(0); 760 } 761 762 /*@ 763 MatProductSymbolic - Perform the symbolic portion of a matrix product, this creates a data structure for use with the numerical produce. 764 765 Collective on Mat 766 767 Input Parameters: 768 . mat - the matrix to hold a product 769 770 Output Parameters: 771 . mat - the matrix product data structure 772 773 Level: intermediate 774 775 .seealso: MatProductCreate(), MatSetType(), MatProductNumeric(), MatProductType, MatProductAlgorithm 776 @*/ 777 PetscErrorCode MatProductSymbolic(Mat mat) 778 { 779 PetscErrorCode ierr; 780 Mat_Product *product = mat->product; 781 MatProductType productype = product->type; 782 PetscLogEvent eventtype=-1; 783 784 PetscFunctionBegin; 785 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 786 787 /* log event */ 788 switch (productype) { 789 case MATPRODUCT_AB: 790 eventtype = MAT_MatMultSymbolic; 791 break; 792 case MATPRODUCT_AtB: 793 eventtype = MAT_TransposeMatMultSymbolic; 794 break; 795 case MATPRODUCT_ABt: 796 eventtype = MAT_MatTransposeMultSymbolic; 797 break; 798 case MATPRODUCT_PtAP: 799 eventtype = MAT_PtAPSymbolic; 800 break; 801 case MATPRODUCT_RARt: 802 eventtype = MAT_RARtSymbolic; 803 break; 804 case MATPRODUCT_ABC: 805 eventtype = MAT_MatMatMultSymbolic; 806 break; 807 default: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"MATPRODUCT type is not supported"); 808 } 809 810 if (mat->ops->productsymbolic) { 811 ierr = PetscLogEventBegin(eventtype,mat,0,0,0);CHKERRQ(ierr); 812 ierr = (*mat->ops->productsymbolic)(mat);CHKERRQ(ierr); 813 ierr = PetscLogEventEnd(eventtype,mat,0,0,0);CHKERRQ(ierr); 814 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_USER,"Call MatProductSetFromOptions() first"); 815 PetscFunctionReturn(0); 816 } 817 818 /*@ 819 MatProductSetFill - Set an expected fill of the matrix product. 820 821 Collective on Mat 822 823 Input Parameters: 824 + mat - the matrix product 825 - fill - expected fill as ratio of nnz(mat)/(nnz(A) + nnz(B) + nnz(C)); use PETSC_DEFAULT if you do not have a good estimate. If the product is a dense matrix, this is irrelevent. 826 827 Level: intermediate 828 829 .seealso: MatProductSetType(), MatProductSetAlgorithm(), MatProductCreate() 830 @*/ 831 PetscErrorCode MatProductSetFill(Mat mat,PetscReal fill) 832 { 833 Mat_Product *product = mat->product; 834 835 PetscFunctionBegin; 836 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 837 838 if (!product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_NULL,"Data struc Mat_Product is not created, call MatProductCreate() first"); 839 if (fill == PETSC_DEFAULT || fill == PETSC_DECIDE) { 840 product->fill = 2.0; 841 } else product->fill = fill; 842 PetscFunctionReturn(0); 843 } 844 845 /*@ 846 MatProductSetAlgorithm - Requests a particular algorithm for a matrix product implementation. 847 848 Collective on Mat 849 850 Input Parameters: 851 + mat - the matrix product 852 - alg - particular implementation algorithm of the matrix product, e.g., MATPRODUCTALGORITHM_DEFAULT. 853 854 Level: intermediate 855 856 .seealso: MatProductSetType(), MatProductSetFill(), MatProductCreate() 857 @*/ 858 PetscErrorCode MatProductSetAlgorithm(Mat mat,MatProductAlgorithm alg) 859 { 860 Mat_Product *product = mat->product; 861 862 PetscFunctionBegin; 863 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 864 865 if (!product) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_NULL,"Data struc Mat_Product is not created, call MatProductCreate() first"); 866 product->alg = alg; 867 PetscFunctionReturn(0); 868 } 869 870 /*@ 871 MatProductSetType - Sets a particular matrix product type, for example Mat*Mat. 872 873 Collective on Mat 874 875 Input Parameters: 876 + mat - the matrix 877 - productype - matrix product type, e.g., MATPRODUCT_AB,MATPRODUCT_AtB,MATPRODUCT_ABt,MATPRODUCT_PtAP,MATPRODUCT_RARt,MATPRODUCT_ABC. 878 879 Level: intermediate 880 881 .seealso: MatProductCreate(), MatProductType, MatProductAlgorithm 882 @*/ 883 PetscErrorCode MatProductSetType(Mat mat,MatProductType productype) 884 { 885 PetscErrorCode ierr; 886 Mat_Product *product = mat->product; 887 MPI_Comm comm; 888 889 PetscFunctionBegin; 890 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 891 892 ierr = PetscObjectGetComm((PetscObject)mat,&comm);CHKERRQ(ierr); 893 if (!product) SETERRQ(comm,PETSC_ERR_ARG_NULL,"Data struc Mat_Product is not created, call MatProductCreate() first"); 894 product->type = productype; 895 896 switch (productype) { 897 case MATPRODUCT_AB: 898 mat->ops->productsetfromoptions = MatProductSetFromOptions_AB; 899 break; 900 case MATPRODUCT_AtB: 901 mat->ops->productsetfromoptions = MatProductSetFromOptions_AtB; 902 break; 903 case MATPRODUCT_ABt: 904 mat->ops->productsetfromoptions = MatProductSetFromOptions_ABt; 905 break; 906 case MATPRODUCT_PtAP: 907 mat->ops->productsetfromoptions = MatProductSetFromOptions_PtAP; 908 break; 909 case MATPRODUCT_RARt: 910 mat->ops->productsetfromoptions = MatProductSetFromOptions_RARt; 911 break; 912 case MATPRODUCT_ABC: 913 mat->ops->productsetfromoptions = MatProductSetFromOptions_ABC; 914 break; 915 default: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"ProductType is not supported\n"); 916 } 917 PetscFunctionReturn(0); 918 } 919 920 /* Create a supporting struct and attach it to the matrix product */ 921 static PetscErrorCode MatProductCreate_Private(Mat A,Mat B,Mat C,Mat D) 922 { 923 PetscErrorCode ierr; 924 Mat_Product *product=NULL; 925 926 PetscFunctionBegin; 927 ierr = PetscNewLog(D,&product);CHKERRQ(ierr); 928 product->A = A; 929 product->B = B; 930 product->C = C; 931 product->Dwork = NULL; 932 product->alg = MATPRODUCTALGORITHM_DEFAULT; 933 product->fill = 2.0; /* PETSC_DEFAULT */ 934 product->Areplaced = PETSC_FALSE; 935 product->Breplaced = PETSC_FALSE; 936 product->api_user = PETSC_FALSE; 937 D->product = product; 938 PetscFunctionReturn(0); 939 } 940 941 /*@ 942 MatProductCreateWithMat - Setup a given matrix as a matrix product. 943 944 Collective on Mat 945 946 Input Parameters: 947 + A - the first matrix 948 . B - the second matrix 949 . C - the third matrix (optional) 950 - D - the matrix which will be used as a product 951 952 Output Parameters: 953 . D - the product matrix 954 955 Level: intermediate 956 957 .seealso: MatProductCreate() 958 @*/ 959 PetscErrorCode MatProductCreateWithMat(Mat A,Mat B,Mat C,Mat D) 960 { 961 PetscErrorCode ierr; 962 963 PetscFunctionBegin; 964 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 965 PetscValidType(A,1); 966 MatCheckPreallocated(A,1); 967 if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 968 if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 969 970 PetscValidHeaderSpecific(B,MAT_CLASSID,2); 971 PetscValidType(B,2); 972 MatCheckPreallocated(B,2); 973 if (!B->assembled) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 974 if (B->factortype) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 975 976 if (C) { 977 PetscValidHeaderSpecific(C,MAT_CLASSID,3); 978 PetscValidType(C,3); 979 MatCheckPreallocated(C,3); 980 if (!C->assembled) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 981 if (C->factortype) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 982 } 983 984 PetscValidHeaderSpecific(D,MAT_CLASSID,4); 985 PetscValidType(D,4); 986 MatCheckPreallocated(D,4); 987 if (!D->assembled) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 988 if (D->factortype) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 989 990 /* Create a supporting struct and attach it to D */ 991 ierr = MatProductCreate_Private(A,B,C,D);CHKERRQ(ierr); 992 PetscFunctionReturn(0); 993 } 994 995 /*@ 996 MatProductCreate - create a matrix product object that can be used to compute various matrix times matrix operations. 997 998 Collective on Mat 999 1000 Input Parameters: 1001 + A - the first matrix 1002 . B - the second matrix 1003 - C - the third matrix (optional) 1004 1005 Output Parameters: 1006 . D - the product matrix 1007 1008 Level: intermediate 1009 1010 .seealso: MatProductCreateWithMat(), MatProductSetType(), MatProductSetAlgorithm() 1011 @*/ 1012 PetscErrorCode MatProductCreate(Mat A,Mat B,Mat C,Mat *D) 1013 { 1014 PetscErrorCode ierr; 1015 1016 PetscFunctionBegin; 1017 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1018 PetscValidType(A,1); 1019 MatCheckPreallocated(A,1); 1020 if (!A->assembled) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1021 if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1022 1023 PetscValidHeaderSpecific(B,MAT_CLASSID,2); 1024 PetscValidType(B,2); 1025 MatCheckPreallocated(B,2); 1026 if (!B->assembled) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1027 if (B->factortype) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1028 1029 if (C) { 1030 PetscValidHeaderSpecific(C,MAT_CLASSID,3); 1031 PetscValidType(C,3); 1032 MatCheckPreallocated(C,3); 1033 if (!C->assembled) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 1034 if (C->factortype) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1035 } 1036 1037 PetscValidPointer(D,4); 1038 1039 ierr = MatCreate(PetscObjectComm((PetscObject)A),D);CHKERRQ(ierr); 1040 ierr = MatProductCreate_Private(A,B,C,*D);CHKERRQ(ierr); 1041 PetscFunctionReturn(0); 1042 } 1043