1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 const char *const MatCompositeMergeTypes[] = {"left","right","MatCompositeMergeType","MAT_COMPOSITE_",0}; 5 6 typedef struct _Mat_CompositeLink *Mat_CompositeLink; 7 struct _Mat_CompositeLink { 8 Mat mat; 9 Vec work; 10 Mat_CompositeLink next,prev; 11 }; 12 13 typedef struct { 14 MatCompositeType type; 15 Mat_CompositeLink head,tail; 16 Vec work; 17 PetscScalar scale; /* scale factor supplied with MatScale() */ 18 Vec left,right; /* left and right diagonal scaling provided with MatDiagonalScale() */ 19 Vec leftwork,rightwork,leftwork2,rightwork2; /* Two pairs of working vectors */ 20 PetscInt nmat; 21 PetscBool merge; 22 MatCompositeMergeType mergetype; 23 MatStructure structure; 24 25 PetscScalar *scalings; 26 PetscBool merge_mvctx; /* Whether need to merge mvctx of component matrices */ 27 Vec *lvecs; /* [nmat] Basically, they are Mvctx->lvec of each component matrix */ 28 PetscScalar *larray; /* [len] Data arrays of lvecs[] are stored consecutively in larray */ 29 PetscInt len; /* Length of larray[] */ 30 Vec gvec; /* Union of lvecs[] without duplicated entries */ 31 PetscInt *location; /* A map that maps entries in garray[] to larray[] */ 32 VecScatter Mvctx; 33 } Mat_Composite; 34 35 PetscErrorCode MatDestroy_Composite(Mat mat) 36 { 37 PetscErrorCode ierr; 38 Mat_Composite *shell = (Mat_Composite*)mat->data; 39 Mat_CompositeLink next = shell->head,oldnext; 40 PetscInt i; 41 42 PetscFunctionBegin; 43 while (next) { 44 ierr = MatDestroy(&next->mat);CHKERRQ(ierr); 45 if (next->work && (!next->next || next->work != next->next->work)) { 46 ierr = VecDestroy(&next->work);CHKERRQ(ierr); 47 } 48 oldnext = next; 49 next = next->next; 50 ierr = PetscFree(oldnext);CHKERRQ(ierr); 51 } 52 ierr = VecDestroy(&shell->work);CHKERRQ(ierr); 53 ierr = VecDestroy(&shell->left);CHKERRQ(ierr); 54 ierr = VecDestroy(&shell->right);CHKERRQ(ierr); 55 ierr = VecDestroy(&shell->leftwork);CHKERRQ(ierr); 56 ierr = VecDestroy(&shell->rightwork);CHKERRQ(ierr); 57 ierr = VecDestroy(&shell->leftwork2);CHKERRQ(ierr); 58 ierr = VecDestroy(&shell->rightwork2);CHKERRQ(ierr); 59 60 if (shell->Mvctx) { 61 for (i=0; i<shell->nmat; i++) {ierr = VecDestroy(&shell->lvecs[i]);CHKERRQ(ierr);} 62 ierr = PetscFree3(shell->location,shell->larray,shell->lvecs);CHKERRQ(ierr); 63 ierr = PetscFree(shell->larray);CHKERRQ(ierr); 64 ierr = VecDestroy(&shell->gvec);CHKERRQ(ierr); 65 ierr = VecScatterDestroy(&shell->Mvctx);CHKERRQ(ierr); 66 } 67 68 ierr = PetscFree(shell->scalings);CHKERRQ(ierr); 69 ierr = PetscFree(mat->data);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 73 PetscErrorCode MatMult_Composite_Multiplicative(Mat A,Vec x,Vec y) 74 { 75 Mat_Composite *shell = (Mat_Composite*)A->data; 76 Mat_CompositeLink next = shell->head; 77 PetscErrorCode ierr; 78 Vec in,out; 79 PetscScalar scale; 80 PetscInt i; 81 82 PetscFunctionBegin; 83 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 84 in = x; 85 if (shell->right) { 86 if (!shell->rightwork) { 87 ierr = VecDuplicate(shell->right,&shell->rightwork);CHKERRQ(ierr); 88 } 89 ierr = VecPointwiseMult(shell->rightwork,shell->right,in);CHKERRQ(ierr); 90 in = shell->rightwork; 91 } 92 while (next->next) { 93 if (!next->work) { /* should reuse previous work if the same size */ 94 ierr = MatCreateVecs(next->mat,NULL,&next->work);CHKERRQ(ierr); 95 } 96 out = next->work; 97 ierr = MatMult(next->mat,in,out);CHKERRQ(ierr); 98 in = out; 99 next = next->next; 100 } 101 ierr = MatMult(next->mat,in,y);CHKERRQ(ierr); 102 if (shell->left) { 103 ierr = VecPointwiseMult(y,shell->left,y);CHKERRQ(ierr); 104 } 105 scale = shell->scale; 106 if (shell->scalings) {for (i=0; i<shell->nmat; i++) scale *= shell->scalings[i];} 107 ierr = VecScale(y,scale);CHKERRQ(ierr); 108 PetscFunctionReturn(0); 109 } 110 111 PetscErrorCode MatMultTranspose_Composite_Multiplicative(Mat A,Vec x,Vec y) 112 { 113 Mat_Composite *shell = (Mat_Composite*)A->data; 114 Mat_CompositeLink tail = shell->tail; 115 PetscErrorCode ierr; 116 Vec in,out; 117 PetscScalar scale; 118 PetscInt i; 119 120 PetscFunctionBegin; 121 if (!tail) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 122 in = x; 123 if (shell->left) { 124 if (!shell->leftwork) { 125 ierr = VecDuplicate(shell->left,&shell->leftwork);CHKERRQ(ierr); 126 } 127 ierr = VecPointwiseMult(shell->leftwork,shell->left,in);CHKERRQ(ierr); 128 in = shell->leftwork; 129 } 130 while (tail->prev) { 131 if (!tail->prev->work) { /* should reuse previous work if the same size */ 132 ierr = MatCreateVecs(tail->mat,NULL,&tail->prev->work);CHKERRQ(ierr); 133 } 134 out = tail->prev->work; 135 ierr = MatMultTranspose(tail->mat,in,out);CHKERRQ(ierr); 136 in = out; 137 tail = tail->prev; 138 } 139 ierr = MatMultTranspose(tail->mat,in,y);CHKERRQ(ierr); 140 if (shell->right) { 141 ierr = VecPointwiseMult(y,shell->right,y);CHKERRQ(ierr); 142 } 143 144 scale = shell->scale; 145 if (shell->scalings) {for (i=0; i<shell->nmat; i++) scale *= shell->scalings[i];} 146 ierr = VecScale(y,scale);CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 } 149 150 PetscErrorCode MatMult_Composite(Mat mat,Vec x,Vec y) 151 { 152 PetscErrorCode ierr; 153 Mat_Composite *shell = (Mat_Composite*)mat->data; 154 Mat_CompositeLink cur = shell->head; 155 Vec in,y2,xin; 156 Mat A,B; 157 PetscInt i,j,k,n,nuniq,lo,hi,mid,*gindices,*buf,*tmp,tot; 158 const PetscScalar *vals; 159 const PetscInt *garray; 160 IS ix,iy; 161 PetscBool match; 162 163 PetscFunctionBegin; 164 if (!cur) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 165 in = x; 166 if (shell->right) { 167 if (!shell->rightwork) { 168 ierr = VecDuplicate(shell->right,&shell->rightwork);CHKERRQ(ierr); 169 } 170 ierr = VecPointwiseMult(shell->rightwork,shell->right,in);CHKERRQ(ierr); 171 in = shell->rightwork; 172 } 173 174 /* Try to merge Mvctx when instructed but not yet done. We did not do it in MatAssemblyEnd() since at that time 175 we did not know whether mat is ADDITIVE or MULTIPLICATIVE. Only now we are assured mat is ADDITIVE and 176 it is legal to merge Mvctx, because all component matrices have the same size. 177 */ 178 if (shell->merge_mvctx && !shell->Mvctx) { 179 /* Currently only implemented for MATMPIAIJ */ 180 for (cur=shell->head; cur; cur=cur->next) { 181 ierr = PetscObjectTypeCompare((PetscObject)cur->mat,MATMPIAIJ,&match);CHKERRQ(ierr); 182 if (!match) { 183 shell->merge_mvctx = PETSC_FALSE; 184 goto skip_merge_mvctx; 185 } 186 } 187 188 /* Go through matrices first time to count total number of nonzero off-diag columns (may have dups) */ 189 tot = 0; 190 for (cur=shell->head; cur; cur=cur->next) { 191 ierr = MatMPIAIJGetSeqAIJ(cur->mat,NULL,&B,NULL);CHKERRQ(ierr); 192 ierr = MatGetLocalSize(B,NULL,&n);CHKERRQ(ierr); 193 tot += n; 194 } 195 ierr = PetscMalloc3(tot,&shell->location,tot,&shell->larray,shell->nmat,&shell->lvecs);CHKERRQ(ierr); 196 shell->len = tot; 197 198 /* Go through matrices second time to sort off-diag columns and remove dups */ 199 ierr = PetscMalloc1(tot,&gindices);CHKERRQ(ierr); /* No Malloc2() since we will give one to petsc and free the other */ 200 ierr = PetscMalloc1(tot,&buf);CHKERRQ(ierr); 201 nuniq = 0; /* Number of unique nonzero columns */ 202 for (cur=shell->head; cur; cur=cur->next) { 203 ierr = MatMPIAIJGetSeqAIJ(cur->mat,NULL,&B,&garray);CHKERRQ(ierr); 204 ierr = MatGetLocalSize(B,NULL,&n);CHKERRQ(ierr); 205 /* Merge pre-sorted garray[0,n) and gindices[0,nuniq) to buf[] */ 206 i = j = k = 0; 207 while (i < n && j < nuniq) { 208 if (garray[i] < gindices[j]) buf[k++] = garray[i++]; 209 else if (garray[i] > gindices[j]) buf[k++] = gindices[j++]; 210 else {buf[k++] = garray[i++]; j++;} 211 } 212 /* Copy leftover in garray[] or gindices[] */ 213 if (i < n) { 214 ierr = PetscArraycpy(buf+k,garray+i,n-i);CHKERRQ(ierr); 215 nuniq = k + n-i; 216 } else if (j < nuniq) { 217 ierr = PetscArraycpy(buf+k,gindices+j,nuniq-j);CHKERRQ(ierr); 218 nuniq = k + nuniq-j; 219 } else nuniq = k; 220 /* Swap gindices and buf to merge garray of the next matrix */ 221 tmp = gindices; 222 gindices = buf; 223 buf = tmp; 224 } 225 ierr = PetscFree(buf);CHKERRQ(ierr); 226 227 /* Go through matrices third time to build a map from gindices[] to garray[] */ 228 tot = 0; 229 for (cur=shell->head,j=0; cur; cur=cur->next,j++) { /* j-th matrix */ 230 ierr = MatMPIAIJGetSeqAIJ(cur->mat,NULL,&B,&garray);CHKERRQ(ierr); 231 ierr = MatGetLocalSize(B,NULL,&n);CHKERRQ(ierr); 232 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,n,NULL,&shell->lvecs[j]);CHKERRQ(ierr); 233 /* This is an optimized PetscFindInt(garray[i],nuniq,gindices,&shell->location[tot+i]), using the fact that garray[] is also sorted */ 234 lo = 0; 235 for (i=0; i<n; i++) { 236 hi = nuniq; 237 while (hi - lo > 1) { 238 mid = lo + (hi - lo)/2; 239 if (garray[i] < gindices[mid]) hi = mid; 240 else lo = mid; 241 } 242 shell->location[tot+i] = lo; /* gindices[lo] = garray[i] */ 243 lo++; /* Since garray[i+1] > garray[i], we can safely advance lo */ 244 } 245 tot += n; 246 } 247 248 /* Build merged Mvctx */ 249 ierr = ISCreateGeneral(PETSC_COMM_SELF,nuniq,gindices,PETSC_OWN_POINTER,&ix);CHKERRQ(ierr); 250 ierr = ISCreateStride(PETSC_COMM_SELF,nuniq,0,1,&iy);CHKERRQ(ierr); 251 ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)mat),1,mat->cmap->n,mat->cmap->N,NULL,&xin);CHKERRQ(ierr); 252 ierr = VecCreateSeq(PETSC_COMM_SELF,nuniq,&shell->gvec);CHKERRQ(ierr); 253 ierr = VecScatterCreate(xin,ix,shell->gvec,iy,&shell->Mvctx);CHKERRQ(ierr); 254 ierr = VecDestroy(&xin);CHKERRQ(ierr); 255 ierr = ISDestroy(&ix);CHKERRQ(ierr); 256 ierr = ISDestroy(&iy);CHKERRQ(ierr); 257 } 258 259 skip_merge_mvctx: 260 ierr = VecSet(y,0);CHKERRQ(ierr); 261 if (!shell->leftwork2) {ierr = VecDuplicate(y,&shell->leftwork2);CHKERRQ(ierr);} 262 y2 = shell->leftwork2; 263 264 if (shell->Mvctx) { /* Have a merged Mvctx */ 265 /* Suppose we want to compute y = sMx, where s is the scaling factor and A, B are matrix M's diagonal/off-diagonal part. We could do 266 in y = s(Ax1 + Bx2) or y = sAx1 + sBx2. The former incurs less FLOPS than the latter, but the latter provides an oppertunity to 267 overlap communication/computation since we can do sAx1 while communicating x2. Here, we use the former approach. 268 */ 269 ierr = VecScatterBegin(shell->Mvctx,in,shell->gvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 270 ierr = VecScatterEnd(shell->Mvctx,in,shell->gvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 271 272 ierr = VecGetArrayRead(shell->gvec,&vals);CHKERRQ(ierr); 273 for (i=0; i<shell->len; i++) shell->larray[i] = vals[shell->location[i]]; 274 ierr = VecRestoreArrayRead(shell->gvec,&vals);CHKERRQ(ierr); 275 276 for (cur=shell->head,tot=i=0; cur; cur=cur->next,i++) { /* i-th matrix */ 277 ierr = MatMPIAIJGetSeqAIJ(cur->mat,&A,&B,NULL);CHKERRQ(ierr); 278 ierr = (*A->ops->mult)(A,in,y2);CHKERRQ(ierr); 279 ierr = MatGetLocalSize(B,NULL,&n);CHKERRQ(ierr); 280 ierr = VecPlaceArray(shell->lvecs[i],&shell->larray[tot]);CHKERRQ(ierr); 281 ierr = (*B->ops->multadd)(B,shell->lvecs[i],y2,y2);CHKERRQ(ierr); 282 ierr = VecResetArray(shell->lvecs[i]);CHKERRQ(ierr); 283 ierr = VecAXPY(y,(shell->scalings ? shell->scalings[i] : 1.0),y2);CHKERRQ(ierr); 284 tot += n; 285 } 286 } else { 287 if (shell->scalings) { 288 for (cur=shell->head,i=0; cur; cur=cur->next,i++) { 289 ierr = MatMult(cur->mat,in,y2);CHKERRQ(ierr); 290 ierr = VecAXPY(y,shell->scalings[i],y2);CHKERRQ(ierr); 291 } 292 } else { 293 for (cur=shell->head; cur; cur=cur->next) {ierr = MatMultAdd(cur->mat,in,y,y);CHKERRQ(ierr);} 294 } 295 } 296 297 if (shell->left) {ierr = VecPointwiseMult(y,shell->left,y);CHKERRQ(ierr);} 298 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 299 PetscFunctionReturn(0); 300 } 301 302 PetscErrorCode MatMultTranspose_Composite(Mat A,Vec x,Vec y) 303 { 304 Mat_Composite *shell = (Mat_Composite*)A->data; 305 Mat_CompositeLink next = shell->head; 306 PetscErrorCode ierr; 307 Vec in,y2 = NULL; 308 PetscInt i; 309 310 PetscFunctionBegin; 311 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 312 in = x; 313 if (shell->left) { 314 if (!shell->leftwork) { 315 ierr = VecDuplicate(shell->left,&shell->leftwork);CHKERRQ(ierr); 316 } 317 ierr = VecPointwiseMult(shell->leftwork,shell->left,in);CHKERRQ(ierr); 318 in = shell->leftwork; 319 } 320 321 ierr = MatMultTranspose(next->mat,in,y);CHKERRQ(ierr); 322 if (shell->scalings) { 323 ierr = VecScale(y,shell->scalings[0]);CHKERRQ(ierr); 324 if (!shell->rightwork2) {ierr = VecDuplicate(y,&shell->rightwork2);CHKERRQ(ierr);} 325 y2 = shell->rightwork2; 326 } 327 i = 1; 328 while ((next = next->next)) { 329 if (!shell->scalings) {ierr = MatMultTransposeAdd(next->mat,in,y,y);CHKERRQ(ierr);} 330 else { 331 ierr = MatMultTranspose(next->mat,in,y2);CHKERRQ(ierr); 332 ierr = VecAXPY(y,shell->scalings[i++],y2);CHKERRQ(ierr); 333 } 334 } 335 if (shell->right) { 336 ierr = VecPointwiseMult(y,shell->right,y);CHKERRQ(ierr); 337 } 338 ierr = VecScale(y,shell->scale);CHKERRQ(ierr); 339 PetscFunctionReturn(0); 340 } 341 342 PetscErrorCode MatMultAdd_Composite(Mat A,Vec x,Vec y,Vec z) 343 { 344 Mat_Composite *shell = (Mat_Composite*)A->data; 345 PetscErrorCode ierr; 346 347 PetscFunctionBegin; 348 if (y != z) { 349 ierr = MatMult(A,x,z);CHKERRQ(ierr); 350 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 351 } else { 352 if (!shell->leftwork) { 353 ierr = VecDuplicate(z,&shell->leftwork);CHKERRQ(ierr); 354 } 355 ierr = MatMult(A,x,shell->leftwork);CHKERRQ(ierr); 356 ierr = VecCopy(y,z);CHKERRQ(ierr); 357 ierr = VecAXPY(z,1.0,shell->leftwork);CHKERRQ(ierr); 358 } 359 PetscFunctionReturn(0); 360 } 361 362 PetscErrorCode MatMultTransposeAdd_Composite(Mat A,Vec x,Vec y, Vec z) 363 { 364 Mat_Composite *shell = (Mat_Composite*)A->data; 365 PetscErrorCode ierr; 366 367 PetscFunctionBegin; 368 if (y != z) { 369 ierr = MatMultTranspose(A,x,z);CHKERRQ(ierr); 370 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 371 } else { 372 if (!shell->rightwork) { 373 ierr = VecDuplicate(z,&shell->rightwork);CHKERRQ(ierr); 374 } 375 ierr = MatMultTranspose(A,x,shell->rightwork);CHKERRQ(ierr); 376 ierr = VecCopy(y,z);CHKERRQ(ierr); 377 ierr = VecAXPY(z,1.0,shell->rightwork);CHKERRQ(ierr); 378 } 379 PetscFunctionReturn(0); 380 } 381 382 PetscErrorCode MatGetDiagonal_Composite(Mat A,Vec v) 383 { 384 Mat_Composite *shell = (Mat_Composite*)A->data; 385 Mat_CompositeLink next = shell->head; 386 PetscErrorCode ierr; 387 PetscInt i; 388 389 PetscFunctionBegin; 390 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 391 if (shell->right || shell->left) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get diagonal if left or right scaling"); 392 393 ierr = MatGetDiagonal(next->mat,v);CHKERRQ(ierr); 394 if (shell->scalings) {ierr = VecScale(v,shell->scalings[0]);CHKERRQ(ierr);} 395 396 if (next->next && !shell->work) { 397 ierr = VecDuplicate(v,&shell->work);CHKERRQ(ierr); 398 } 399 i = 1; 400 while ((next = next->next)) { 401 ierr = MatGetDiagonal(next->mat,shell->work);CHKERRQ(ierr); 402 ierr = VecAXPY(v,(shell->scalings ? shell->scalings[i++] : 1.0),shell->work);CHKERRQ(ierr); 403 } 404 ierr = VecScale(v,shell->scale);CHKERRQ(ierr); 405 PetscFunctionReturn(0); 406 } 407 408 PetscErrorCode MatAssemblyEnd_Composite(Mat Y,MatAssemblyType t) 409 { 410 Mat_Composite *shell = (Mat_Composite*)Y->data; 411 PetscErrorCode ierr; 412 413 PetscFunctionBegin; 414 if (shell->merge) { 415 ierr = MatCompositeMerge(Y);CHKERRQ(ierr); 416 } 417 418 PetscFunctionReturn(0); 419 } 420 421 PetscErrorCode MatScale_Composite(Mat inA,PetscScalar alpha) 422 { 423 Mat_Composite *a = (Mat_Composite*)inA->data; 424 425 PetscFunctionBegin; 426 a->scale *= alpha; 427 PetscFunctionReturn(0); 428 } 429 430 PetscErrorCode MatDiagonalScale_Composite(Mat inA,Vec left,Vec right) 431 { 432 Mat_Composite *a = (Mat_Composite*)inA->data; 433 PetscErrorCode ierr; 434 435 PetscFunctionBegin; 436 if (left) { 437 if (!a->left) { 438 ierr = VecDuplicate(left,&a->left);CHKERRQ(ierr); 439 ierr = VecCopy(left,a->left);CHKERRQ(ierr); 440 } else { 441 ierr = VecPointwiseMult(a->left,left,a->left);CHKERRQ(ierr); 442 } 443 } 444 if (right) { 445 if (!a->right) { 446 ierr = VecDuplicate(right,&a->right);CHKERRQ(ierr); 447 ierr = VecCopy(right,a->right);CHKERRQ(ierr); 448 } else { 449 ierr = VecPointwiseMult(a->right,right,a->right);CHKERRQ(ierr); 450 } 451 } 452 PetscFunctionReturn(0); 453 } 454 455 PetscErrorCode MatSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,Mat A) 456 { 457 Mat_Composite *a = (Mat_Composite*)A->data; 458 PetscErrorCode ierr; 459 460 PetscFunctionBegin; 461 ierr = PetscOptionsHead(PetscOptionsObject,"MATCOMPOSITE options");CHKERRQ(ierr); 462 ierr = PetscOptionsBool("-mat_composite_merge","Merge at MatAssemblyEnd","MatCompositeMerge",a->merge,&a->merge,NULL);CHKERRQ(ierr); 463 ierr = PetscOptionsEnum("-mat_composite_merge_type","Set composite merge direction","MatCompositeSetMergeType",MatCompositeMergeTypes,(PetscEnum)a->mergetype,(PetscEnum*)&a->mergetype,NULL);CHKERRQ(ierr); 464 ierr = PetscOptionsBool("-mat_composite_merge_mvctx","Merge MatMult() vecscat contexts","MatCreateComposite",a->merge_mvctx,&a->merge_mvctx,NULL);CHKERRQ(ierr); 465 ierr = PetscOptionsTail();CHKERRQ(ierr); 466 PetscFunctionReturn(0); 467 } 468 469 /*@ 470 MatCreateComposite - Creates a matrix as the sum or product of one or more matrices 471 472 Collective 473 474 Input Parameters: 475 + comm - MPI communicator 476 . nmat - number of matrices to put in 477 - mats - the matrices 478 479 Output Parameter: 480 . mat - the matrix 481 482 Options Database Keys: 483 + -mat_composite_merge - merge in MatAssemblyEnd() 484 . -mat_composite_merge_mvctx - merge Mvctx of component matrices to optimize communication in MatMult() for ADDITIVE matrices 485 - -mat_composite_merge_type - set merge direction 486 487 Level: advanced 488 489 Notes: 490 Alternative construction 491 $ MatCreate(comm,&mat); 492 $ MatSetSizes(mat,m,n,M,N); 493 $ MatSetType(mat,MATCOMPOSITE); 494 $ MatCompositeAddMat(mat,mats[0]); 495 $ .... 496 $ MatCompositeAddMat(mat,mats[nmat-1]); 497 $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY); 498 $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY); 499 500 For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0] 501 502 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeGetMat(), MatCompositeMerge(), MatCompositeSetType(), MATCOMPOSITE 503 504 @*/ 505 PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat) 506 { 507 PetscErrorCode ierr; 508 PetscInt m,n,M,N,i; 509 510 PetscFunctionBegin; 511 if (nmat < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix"); 512 PetscValidPointer(mat,3); 513 514 ierr = MatGetLocalSize(mats[0],PETSC_IGNORE,&n);CHKERRQ(ierr); 515 ierr = MatGetLocalSize(mats[nmat-1],&m,PETSC_IGNORE);CHKERRQ(ierr); 516 ierr = MatGetSize(mats[0],PETSC_IGNORE,&N);CHKERRQ(ierr); 517 ierr = MatGetSize(mats[nmat-1],&M,PETSC_IGNORE);CHKERRQ(ierr); 518 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 519 ierr = MatSetSizes(*mat,m,n,M,N);CHKERRQ(ierr); 520 ierr = MatSetType(*mat,MATCOMPOSITE);CHKERRQ(ierr); 521 for (i=0; i<nmat; i++) { 522 ierr = MatCompositeAddMat(*mat,mats[i]);CHKERRQ(ierr); 523 } 524 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 525 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 526 PetscFunctionReturn(0); 527 } 528 529 530 static PetscErrorCode MatCompositeAddMat_Composite(Mat mat,Mat smat) 531 { 532 Mat_Composite *shell = (Mat_Composite*)mat->data; 533 Mat_CompositeLink ilink,next = shell->head; 534 PetscErrorCode ierr; 535 536 PetscFunctionBegin; 537 ierr = PetscNewLog(mat,&ilink);CHKERRQ(ierr); 538 ilink->next = 0; 539 ierr = PetscObjectReference((PetscObject)smat);CHKERRQ(ierr); 540 ilink->mat = smat; 541 542 if (!next) shell->head = ilink; 543 else { 544 while (next->next) { 545 next = next->next; 546 } 547 next->next = ilink; 548 ilink->prev = next; 549 } 550 shell->tail = ilink; 551 shell->nmat += 1; 552 553 /* Retain the old scalings (if any) and expand it with a 1.0 for the newly added matrix */ 554 if (shell->scalings) { 555 ierr = PetscRealloc(sizeof(PetscScalar)*shell->nmat,&shell->scalings);CHKERRQ(ierr); 556 shell->scalings[shell->nmat-1] = 1.0; 557 } 558 PetscFunctionReturn(0); 559 } 560 561 /*@ 562 MatCompositeAddMat - Add another matrix to a composite matrix. 563 564 Collective on Mat 565 566 Input Parameters: 567 + mat - the composite matrix 568 - smat - the partial matrix 569 570 Level: advanced 571 572 .seealso: MatCreateComposite(), MatCompositeGetMat(), MATCOMPOSITE 573 @*/ 574 PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat) 575 { 576 PetscErrorCode ierr; 577 578 PetscFunctionBegin; 579 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 580 PetscValidHeaderSpecific(smat,MAT_CLASSID,2); 581 ierr = PetscUseMethod(mat,"MatCompositeAddMat_C",(Mat,Mat),(mat,smat));CHKERRQ(ierr); 582 PetscFunctionReturn(0); 583 } 584 585 static PetscErrorCode MatCompositeSetType_Composite(Mat mat,MatCompositeType type) 586 { 587 Mat_Composite *b = (Mat_Composite*)mat->data; 588 589 PetscFunctionBegin; 590 b->type = type; 591 if (type == MAT_COMPOSITE_MULTIPLICATIVE) { 592 mat->ops->getdiagonal = 0; 593 mat->ops->mult = MatMult_Composite_Multiplicative; 594 mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative; 595 b->merge_mvctx = PETSC_FALSE; 596 } else { 597 mat->ops->getdiagonal = MatGetDiagonal_Composite; 598 mat->ops->mult = MatMult_Composite; 599 mat->ops->multtranspose = MatMultTranspose_Composite; 600 } 601 PetscFunctionReturn(0); 602 } 603 604 /*@ 605 MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product. 606 607 Logically Collective on Mat 608 609 Input Parameters: 610 . mat - the composite matrix 611 612 Level: advanced 613 614 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MatCompositeGetType(), MATCOMPOSITE 615 616 @*/ 617 PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type) 618 { 619 PetscErrorCode ierr; 620 621 PetscFunctionBegin; 622 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 623 PetscValidLogicalCollectiveEnum(mat,type,2); 624 ierr = PetscUseMethod(mat,"MatCompositeSetType_C",(Mat,MatCompositeType),(mat,type));CHKERRQ(ierr); 625 PetscFunctionReturn(0); 626 } 627 628 static PetscErrorCode MatCompositeGetType_Composite(Mat mat,MatCompositeType *type) 629 { 630 Mat_Composite *b = (Mat_Composite*)mat->data; 631 632 PetscFunctionBegin; 633 *type = b->type; 634 PetscFunctionReturn(0); 635 } 636 637 /*@ 638 MatCompositeGetType - Returns type of composite. 639 640 Not Collective 641 642 Input Parameter: 643 . mat - the composite matrix 644 645 Output Parameter: 646 . type - type of composite 647 648 Level: advanced 649 650 .seealso: MatCreateComposite(), MatCompositeSetType(), MATCOMPOSITE 651 652 @*/ 653 PetscErrorCode MatCompositeGetType(Mat mat,MatCompositeType *type) 654 { 655 PetscErrorCode ierr; 656 657 PetscFunctionBegin; 658 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 659 PetscValidPointer(type,2); 660 ierr = PetscUseMethod(mat,"MatCompositeGetType_C",(Mat,MatCompositeType*),(mat,type));CHKERRQ(ierr); 661 PetscFunctionReturn(0); 662 } 663 664 static PetscErrorCode MatCompositeSetMatStructure_Composite(Mat mat,MatStructure str) 665 { 666 Mat_Composite *b = (Mat_Composite*)mat->data; 667 668 PetscFunctionBegin; 669 b->structure = str; 670 PetscFunctionReturn(0); 671 } 672 673 /*@ 674 MatCompositeSetMatStructure - Indicates structure of matrices in the composite matrix. 675 676 Not Collective 677 678 Input Parameters: 679 + mat - the composite matrix 680 - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN (default) or SUBSET_NONZERO_PATTERN 681 682 Level: advanced 683 684 Notes: 685 Information about the matrices structure is used in MatCompositeMerge() for additive composite matrix. 686 687 .seealso: MatAXPY(), MatCreateComposite(), MatCompositeMerge() MatCompositeGetMatStructure(), MATCOMPOSITE 688 689 @*/ 690 PetscErrorCode MatCompositeSetMatStructure(Mat mat,MatStructure str) 691 { 692 PetscErrorCode ierr; 693 694 PetscFunctionBegin; 695 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 696 ierr = PetscUseMethod(mat,"MatCompositeSetMatStructure_C",(Mat,MatStructure),(mat,str));CHKERRQ(ierr); 697 PetscFunctionReturn(0); 698 } 699 700 static PetscErrorCode MatCompositeGetMatStructure_Composite(Mat mat,MatStructure *str) 701 { 702 Mat_Composite *b = (Mat_Composite*)mat->data; 703 704 PetscFunctionBegin; 705 *str = b->structure; 706 PetscFunctionReturn(0); 707 } 708 709 /*@ 710 MatCompositeGetMatStructure - Returns the structure of matrices in the composite matrix. 711 712 Not Collective 713 714 Input Parameter: 715 . mat - the composite matrix 716 717 Output Parameter: 718 . str - structure of the matrices 719 720 Level: advanced 721 722 .seealso: MatCreateComposite(), MatCompositeSetMatStructure(), MATCOMPOSITE 723 724 @*/ 725 PetscErrorCode MatCompositeGetMatStructure(Mat mat,MatStructure *str) 726 { 727 PetscErrorCode ierr; 728 729 PetscFunctionBegin; 730 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 731 PetscValidPointer(str,2); 732 ierr = PetscUseMethod(mat,"MatCompositeGetMatStructure_C",(Mat,MatStructure*),(mat,str));CHKERRQ(ierr); 733 PetscFunctionReturn(0); 734 } 735 736 static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat,MatCompositeMergeType type) 737 { 738 Mat_Composite *shell = (Mat_Composite*)mat->data; 739 740 PetscFunctionBegin; 741 shell->mergetype = type; 742 PetscFunctionReturn(0); 743 } 744 745 /*@ 746 MatCompositeSetMergeType - Sets order of MatCompositeMerge(). 747 748 Logically Collective on Mat 749 750 Input Parameter: 751 + mat - the composite matrix 752 - type - MAT_COMPOSITE_MERGE RIGHT (default) to start merge from right with the first added matrix (mat[0]), 753 MAT_COMPOSITE_MERGE_LEFT to start merge from left with the last added matrix (mat[nmat-1]) 754 755 Level: advanced 756 757 Notes: 758 The resulting matrix is the same regardles of the MergeType. Only the order of operation is changed. 759 If set to MAT_COMPOSITE_MERGE_RIGHT the order of the merge is mat[nmat-1]*(mat[nmat-2]*(...*(mat[1]*mat[0]))) 760 otherwise the order is (((mat[nmat-1]*mat[nmat-2])*mat[nmat-3])*...)*mat[0]. 761 762 .seealso: MatCreateComposite(), MatCompositeMerge(), MATCOMPOSITE 763 764 @*/ 765 PetscErrorCode MatCompositeSetMergeType(Mat mat,MatCompositeMergeType type) 766 { 767 PetscErrorCode ierr; 768 769 PetscFunctionBegin; 770 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 771 PetscValidLogicalCollectiveEnum(mat,type,2); 772 ierr = PetscUseMethod(mat,"MatCompositeSetMergeType_C",(Mat,MatCompositeMergeType),(mat,type));CHKERRQ(ierr); 773 PetscFunctionReturn(0); 774 } 775 776 static PetscErrorCode MatCompositeMerge_Composite(Mat mat) 777 { 778 Mat_Composite *shell = (Mat_Composite*)mat->data; 779 Mat_CompositeLink next = shell->head, prev = shell->tail; 780 PetscErrorCode ierr; 781 Mat tmat,newmat; 782 Vec left,right; 783 PetscScalar scale; 784 PetscInt i; 785 786 PetscFunctionBegin; 787 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 788 789 PetscFunctionBegin; 790 scale = shell->scale; 791 if (shell->type == MAT_COMPOSITE_ADDITIVE) { 792 if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) { 793 i = 0; 794 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 795 if (shell->scalings) {ierr = MatScale(tmat,shell->scalings[i++]);CHKERRQ(ierr);} 796 while ((next = next->next)) { 797 ierr = MatAXPY(tmat,(shell->scalings ? shell->scalings[i++] : 1.0),next->mat,shell->structure);CHKERRQ(ierr); 798 } 799 } else { 800 i = shell->nmat-1; 801 ierr = MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 802 if (shell->scalings) {ierr = MatScale(tmat,shell->scalings[i--]);CHKERRQ(ierr);} 803 while ((prev = prev->prev)) { 804 ierr = MatAXPY(tmat,(shell->scalings ? shell->scalings[i--] : 1.0),prev->mat,shell->structure);CHKERRQ(ierr); 805 } 806 } 807 } else { 808 if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) { 809 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 810 while ((next = next->next)) { 811 ierr = MatMatMult(next->mat,tmat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr); 812 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 813 tmat = newmat; 814 } 815 } else { 816 ierr = MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 817 while ((prev = prev->prev)) { 818 ierr = MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr); 819 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 820 tmat = newmat; 821 } 822 } 823 if (shell->scalings) {for (i=0; i<shell->nmat; i++) scale *= shell->scalings[i];} 824 } 825 826 if ((left = shell->left)) {ierr = PetscObjectReference((PetscObject)left);CHKERRQ(ierr);} 827 if ((right = shell->right)) {ierr = PetscObjectReference((PetscObject)right);CHKERRQ(ierr);} 828 829 ierr = MatHeaderReplace(mat,&tmat);CHKERRQ(ierr); 830 831 ierr = MatDiagonalScale(mat,left,right);CHKERRQ(ierr); 832 ierr = MatScale(mat,scale);CHKERRQ(ierr); 833 ierr = VecDestroy(&left);CHKERRQ(ierr); 834 ierr = VecDestroy(&right);CHKERRQ(ierr); 835 PetscFunctionReturn(0); 836 } 837 838 /*@ 839 MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix 840 by summing or computing the product of all the matrices inside the composite matrix. 841 842 Collective 843 844 Input Parameters: 845 . mat - the composite matrix 846 847 848 Options Database Keys: 849 + -mat_composite_merge - merge in MatAssemblyEnd() 850 - -mat_composite_merge_type - set merge direction 851 852 Level: advanced 853 854 Notes: 855 The MatType of the resulting matrix will be the same as the MatType of the FIRST 856 matrix in the composite matrix. 857 858 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MatCompositeSetMatStructure(), MatCompositeSetMergeType(), MATCOMPOSITE 859 860 @*/ 861 PetscErrorCode MatCompositeMerge(Mat mat) 862 { 863 PetscErrorCode ierr; 864 865 PetscFunctionBegin; 866 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 867 ierr = PetscUseMethod(mat,"MatCompositeMerge_C",(Mat),(mat));CHKERRQ(ierr); 868 PetscFunctionReturn(0); 869 } 870 871 static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat,PetscInt *nmat) 872 { 873 Mat_Composite *shell = (Mat_Composite*)mat->data; 874 875 PetscFunctionBegin; 876 *nmat = shell->nmat; 877 PetscFunctionReturn(0); 878 } 879 880 /*@ 881 MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix. 882 883 Not Collective 884 885 Input Parameter: 886 . mat - the composite matrix 887 888 Output Parameter: 889 . nmat - number of matrices in the composite matrix 890 891 Level: advanced 892 893 .seealso: MatCreateComposite(), MatCompositeGetMat(), MATCOMPOSITE 894 895 @*/ 896 PetscErrorCode MatCompositeGetNumberMat(Mat mat,PetscInt *nmat) 897 { 898 PetscErrorCode ierr; 899 900 PetscFunctionBegin; 901 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 902 PetscValidPointer(nmat,2); 903 ierr = PetscUseMethod(mat,"MatCompositeGetNumberMat_C",(Mat,PetscInt*),(mat,nmat));CHKERRQ(ierr); 904 PetscFunctionReturn(0); 905 } 906 907 static PetscErrorCode MatCompositeGetMat_Composite(Mat mat,PetscInt i,Mat *Ai) 908 { 909 Mat_Composite *shell = (Mat_Composite*)mat->data; 910 Mat_CompositeLink ilink; 911 PetscInt k; 912 913 PetscFunctionBegin; 914 if (i >= shell->nmat) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_OUTOFRANGE,"index out of range: %d >= %d",i,shell->nmat); 915 ilink = shell->head; 916 for (k=0; k<i; k++) { 917 ilink = ilink->next; 918 } 919 *Ai = ilink->mat; 920 PetscFunctionReturn(0); 921 } 922 923 /*@ 924 MatCompositeGetMat - Returns the ith matrix from the composite matrix. 925 926 Logically Collective on Mat 927 928 Input Parameter: 929 + mat - the composite matrix 930 - i - the number of requested matrix 931 932 Output Parameter: 933 . Ai - ith matrix in composite 934 935 Level: advanced 936 937 .seealso: MatCreateComposite(), MatCompositeGetNumberMat(), MatCompositeAddMat(), MATCOMPOSITE 938 939 @*/ 940 PetscErrorCode MatCompositeGetMat(Mat mat,PetscInt i,Mat *Ai) 941 { 942 PetscErrorCode ierr; 943 944 PetscFunctionBegin; 945 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 946 PetscValidLogicalCollectiveInt(mat,i,2); 947 PetscValidPointer(Ai,3); 948 ierr = PetscUseMethod(mat,"MatCompositeGetMat_C",(Mat,PetscInt,Mat*),(mat,i,Ai));CHKERRQ(ierr); 949 PetscFunctionReturn(0); 950 } 951 952 PetscErrorCode MatCompositeSetScalings_Composite(Mat mat,const PetscScalar *scalings) 953 { 954 PetscErrorCode ierr; 955 Mat_Composite *shell = (Mat_Composite*)mat->data; 956 PetscInt nmat; 957 958 PetscFunctionBegin; 959 ierr = MatCompositeGetNumberMat(mat,&nmat);CHKERRQ(ierr); 960 if (!shell->scalings) {ierr = PetscMalloc1(nmat,&shell->scalings);CHKERRQ(ierr);} 961 ierr = PetscArraycpy(shell->scalings,scalings,nmat);CHKERRQ(ierr); 962 PetscFunctionReturn(0); 963 } 964 965 /*@ 966 MatCompositeSetScalings - Sets separate scaling factors for component matrices. 967 968 Logically Collective on Mat 969 970 Input Parameter: 971 + mat - the composite matrix 972 - scalings - array of scaling factors with scalings[i] being factor of i-th matrix, for i in [0, nmat) 973 974 Level: advanced 975 976 .seealso: MatScale(), MatDiagonalScale(), MATCOMPOSITE 977 978 @*/ 979 PetscErrorCode MatCompositeSetScalings(Mat mat,const PetscScalar *scalings) 980 { 981 PetscErrorCode ierr; 982 983 PetscFunctionBegin; 984 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 985 PetscValidPointer(scalings,2); 986 PetscValidLogicalCollectiveScalar(mat,*scalings,2); 987 ierr = PetscUseMethod(mat,"MatCompositeSetScalings_C",(Mat,const PetscScalar*),(mat,scalings));CHKERRQ(ierr); 988 PetscFunctionReturn(0); 989 } 990 991 static struct _MatOps MatOps_Values = {0, 992 0, 993 0, 994 MatMult_Composite, 995 MatMultAdd_Composite, 996 /* 5*/ MatMultTranspose_Composite, 997 MatMultTransposeAdd_Composite, 998 0, 999 0, 1000 0, 1001 /* 10*/ 0, 1002 0, 1003 0, 1004 0, 1005 0, 1006 /* 15*/ 0, 1007 0, 1008 MatGetDiagonal_Composite, 1009 MatDiagonalScale_Composite, 1010 0, 1011 /* 20*/ 0, 1012 MatAssemblyEnd_Composite, 1013 0, 1014 0, 1015 /* 24*/ 0, 1016 0, 1017 0, 1018 0, 1019 0, 1020 /* 29*/ 0, 1021 0, 1022 0, 1023 0, 1024 0, 1025 /* 34*/ 0, 1026 0, 1027 0, 1028 0, 1029 0, 1030 /* 39*/ 0, 1031 0, 1032 0, 1033 0, 1034 0, 1035 /* 44*/ 0, 1036 MatScale_Composite, 1037 MatShift_Basic, 1038 0, 1039 0, 1040 /* 49*/ 0, 1041 0, 1042 0, 1043 0, 1044 0, 1045 /* 54*/ 0, 1046 0, 1047 0, 1048 0, 1049 0, 1050 /* 59*/ 0, 1051 MatDestroy_Composite, 1052 0, 1053 0, 1054 0, 1055 /* 64*/ 0, 1056 0, 1057 0, 1058 0, 1059 0, 1060 /* 69*/ 0, 1061 0, 1062 0, 1063 0, 1064 0, 1065 /* 74*/ 0, 1066 0, 1067 MatSetFromOptions_Composite, 1068 0, 1069 0, 1070 /* 79*/ 0, 1071 0, 1072 0, 1073 0, 1074 0, 1075 /* 84*/ 0, 1076 0, 1077 0, 1078 0, 1079 0, 1080 /* 89*/ 0, 1081 0, 1082 0, 1083 0, 1084 0, 1085 /* 94*/ 0, 1086 0, 1087 0, 1088 0, 1089 0, 1090 /*99*/ 0, 1091 0, 1092 0, 1093 0, 1094 0, 1095 /*104*/ 0, 1096 0, 1097 0, 1098 0, 1099 0, 1100 /*109*/ 0, 1101 0, 1102 0, 1103 0, 1104 0, 1105 /*114*/ 0, 1106 0, 1107 0, 1108 0, 1109 0, 1110 /*119*/ 0, 1111 0, 1112 0, 1113 0, 1114 0, 1115 /*124*/ 0, 1116 0, 1117 0, 1118 0, 1119 0, 1120 /*129*/ 0, 1121 0, 1122 0, 1123 0, 1124 0, 1125 /*134*/ 0, 1126 0, 1127 0, 1128 0, 1129 0, 1130 /*139*/ 0, 1131 0, 1132 0 1133 }; 1134 1135 /*MC 1136 MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices. 1137 The matrices need to have a correct size and parallel layout for the sum or product to be valid. 1138 1139 Notes: 1140 to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE); 1141 1142 Level: advanced 1143 1144 .seealso: MatCreateComposite(), MatCompositeSetScalings(), MatCompositeAddMat(), MatSetType(), MatCompositeSetType(), MatCompositeGetType(), MatCompositeSetMatStructure(), MatCompositeGetMatStructure(), MatCompositeMerge(), MatCompositeSetMergeType(), MatCompositeGetNumberMat(), MatCompositeGetMat() 1145 M*/ 1146 1147 PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A) 1148 { 1149 Mat_Composite *b; 1150 PetscErrorCode ierr; 1151 1152 PetscFunctionBegin; 1153 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 1154 A->data = (void*)b; 1155 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1156 1157 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1158 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1159 1160 A->assembled = PETSC_TRUE; 1161 A->preallocated = PETSC_TRUE; 1162 b->type = MAT_COMPOSITE_ADDITIVE; 1163 b->scale = 1.0; 1164 b->nmat = 0; 1165 b->merge = PETSC_FALSE; 1166 b->mergetype = MAT_COMPOSITE_MERGE_RIGHT; 1167 b->structure = DIFFERENT_NONZERO_PATTERN; 1168 b->merge_mvctx = PETSC_TRUE; 1169 1170 1171 1172 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeAddMat_C",MatCompositeAddMat_Composite);CHKERRQ(ierr); 1173 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetType_C",MatCompositeSetType_Composite);CHKERRQ(ierr); 1174 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetType_C",MatCompositeGetType_Composite);CHKERRQ(ierr); 1175 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetMergeType_C",MatCompositeSetMergeType_Composite);CHKERRQ(ierr); 1176 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetMatStructure_C",MatCompositeSetMatStructure_Composite);CHKERRQ(ierr); 1177 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMatStructure_C",MatCompositeGetMatStructure_Composite);CHKERRQ(ierr); 1178 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeMerge_C",MatCompositeMerge_Composite);CHKERRQ(ierr); 1179 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetNumberMat_C",MatCompositeGetNumberMat_Composite);CHKERRQ(ierr); 1180 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMat_C",MatCompositeGetMat_Composite);CHKERRQ(ierr); 1181 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetScalings_C",MatCompositeSetScalings_Composite);CHKERRQ(ierr); 1182 1183 ierr = PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);CHKERRQ(ierr); 1184 PetscFunctionReturn(0); 1185 } 1186 1187