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