1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 const char *const MatCompositeMergeTypes[] = {"left","right","MatCompositeMergeType","MAT_COMPOSITE_",NULL}; 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 PetscFunctionReturn(0); 418 } 419 420 PetscErrorCode MatScale_Composite(Mat inA,PetscScalar alpha) 421 { 422 Mat_Composite *a = (Mat_Composite*)inA->data; 423 424 PetscFunctionBegin; 425 a->scale *= alpha; 426 PetscFunctionReturn(0); 427 } 428 429 PetscErrorCode MatDiagonalScale_Composite(Mat inA,Vec left,Vec right) 430 { 431 Mat_Composite *a = (Mat_Composite*)inA->data; 432 PetscErrorCode ierr; 433 434 PetscFunctionBegin; 435 if (left) { 436 if (!a->left) { 437 ierr = VecDuplicate(left,&a->left);CHKERRQ(ierr); 438 ierr = VecCopy(left,a->left);CHKERRQ(ierr); 439 } else { 440 ierr = VecPointwiseMult(a->left,left,a->left);CHKERRQ(ierr); 441 } 442 } 443 if (right) { 444 if (!a->right) { 445 ierr = VecDuplicate(right,&a->right);CHKERRQ(ierr); 446 ierr = VecCopy(right,a->right);CHKERRQ(ierr); 447 } else { 448 ierr = VecPointwiseMult(a->right,right,a->right);CHKERRQ(ierr); 449 } 450 } 451 PetscFunctionReturn(0); 452 } 453 454 PetscErrorCode MatSetFromOptions_Composite(PetscOptionItems *PetscOptionsObject,Mat A) 455 { 456 Mat_Composite *a = (Mat_Composite*)A->data; 457 PetscErrorCode ierr; 458 459 PetscFunctionBegin; 460 ierr = PetscOptionsHead(PetscOptionsObject,"MATCOMPOSITE options");CHKERRQ(ierr); 461 ierr = PetscOptionsBool("-mat_composite_merge","Merge at MatAssemblyEnd","MatCompositeMerge",a->merge,&a->merge,NULL);CHKERRQ(ierr); 462 ierr = PetscOptionsEnum("-mat_composite_merge_type","Set composite merge direction","MatCompositeSetMergeType",MatCompositeMergeTypes,(PetscEnum)a->mergetype,(PetscEnum*)&a->mergetype,NULL);CHKERRQ(ierr); 463 ierr = PetscOptionsBool("-mat_composite_merge_mvctx","Merge MatMult() vecscat contexts","MatCreateComposite",a->merge_mvctx,&a->merge_mvctx,NULL);CHKERRQ(ierr); 464 ierr = PetscOptionsTail();CHKERRQ(ierr); 465 PetscFunctionReturn(0); 466 } 467 468 /*@ 469 MatCreateComposite - Creates a matrix as the sum or product of one or more matrices 470 471 Collective 472 473 Input Parameters: 474 + comm - MPI communicator 475 . nmat - number of matrices to put in 476 - mats - the matrices 477 478 Output Parameter: 479 . mat - the matrix 480 481 Options Database Keys: 482 + -mat_composite_merge - merge in MatAssemblyEnd() 483 . -mat_composite_merge_mvctx - merge Mvctx of component matrices to optimize communication in MatMult() for ADDITIVE matrices 484 - -mat_composite_merge_type - set merge direction 485 486 Level: advanced 487 488 Notes: 489 Alternative construction 490 $ MatCreate(comm,&mat); 491 $ MatSetSizes(mat,m,n,M,N); 492 $ MatSetType(mat,MATCOMPOSITE); 493 $ MatCompositeAddMat(mat,mats[0]); 494 $ .... 495 $ MatCompositeAddMat(mat,mats[nmat-1]); 496 $ MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY); 497 $ MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY); 498 499 For the multiplicative form the product is mat[nmat-1]*mat[nmat-2]*....*mat[0] 500 501 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCompositeGetMat(), MatCompositeMerge(), MatCompositeSetType(), MATCOMPOSITE 502 503 @*/ 504 PetscErrorCode MatCreateComposite(MPI_Comm comm,PetscInt nmat,const Mat *mats,Mat *mat) 505 { 506 PetscErrorCode ierr; 507 PetscInt m,n,M,N,i; 508 509 PetscFunctionBegin; 510 if (nmat < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must pass in at least one matrix"); 511 PetscValidPointer(mat,4); 512 513 ierr = MatGetLocalSize(mats[0],PETSC_IGNORE,&n);CHKERRQ(ierr); 514 ierr = MatGetLocalSize(mats[nmat-1],&m,PETSC_IGNORE);CHKERRQ(ierr); 515 ierr = MatGetSize(mats[0],PETSC_IGNORE,&N);CHKERRQ(ierr); 516 ierr = MatGetSize(mats[nmat-1],&M,PETSC_IGNORE);CHKERRQ(ierr); 517 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 518 ierr = MatSetSizes(*mat,m,n,M,N);CHKERRQ(ierr); 519 ierr = MatSetType(*mat,MATCOMPOSITE);CHKERRQ(ierr); 520 for (i=0; i<nmat; i++) { 521 ierr = MatCompositeAddMat(*mat,mats[i]);CHKERRQ(ierr); 522 } 523 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 524 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 525 PetscFunctionReturn(0); 526 } 527 528 static PetscErrorCode MatCompositeAddMat_Composite(Mat mat,Mat smat) 529 { 530 Mat_Composite *shell = (Mat_Composite*)mat->data; 531 Mat_CompositeLink ilink,next = shell->head; 532 PetscErrorCode ierr; 533 534 PetscFunctionBegin; 535 ierr = PetscNewLog(mat,&ilink);CHKERRQ(ierr); 536 ilink->next = NULL; 537 ierr = PetscObjectReference((PetscObject)smat);CHKERRQ(ierr); 538 ilink->mat = smat; 539 540 if (!next) shell->head = ilink; 541 else { 542 while (next->next) { 543 next = next->next; 544 } 545 next->next = ilink; 546 ilink->prev = next; 547 } 548 shell->tail = ilink; 549 shell->nmat += 1; 550 551 /* Retain the old scalings (if any) and expand it with a 1.0 for the newly added matrix */ 552 if (shell->scalings) { 553 ierr = PetscRealloc(sizeof(PetscScalar)*shell->nmat,&shell->scalings);CHKERRQ(ierr); 554 shell->scalings[shell->nmat-1] = 1.0; 555 } 556 PetscFunctionReturn(0); 557 } 558 559 /*@ 560 MatCompositeAddMat - Add another matrix to a composite matrix. 561 562 Collective on Mat 563 564 Input Parameters: 565 + mat - the composite matrix 566 - smat - the partial matrix 567 568 Level: advanced 569 570 .seealso: MatCreateComposite(), MatCompositeGetMat(), MATCOMPOSITE 571 @*/ 572 PetscErrorCode MatCompositeAddMat(Mat mat,Mat smat) 573 { 574 PetscErrorCode ierr; 575 576 PetscFunctionBegin; 577 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 578 PetscValidHeaderSpecific(smat,MAT_CLASSID,2); 579 ierr = PetscUseMethod(mat,"MatCompositeAddMat_C",(Mat,Mat),(mat,smat));CHKERRQ(ierr); 580 PetscFunctionReturn(0); 581 } 582 583 static PetscErrorCode MatCompositeSetType_Composite(Mat mat,MatCompositeType type) 584 { 585 Mat_Composite *b = (Mat_Composite*)mat->data; 586 587 PetscFunctionBegin; 588 b->type = type; 589 if (type == MAT_COMPOSITE_MULTIPLICATIVE) { 590 mat->ops->getdiagonal = NULL; 591 mat->ops->mult = MatMult_Composite_Multiplicative; 592 mat->ops->multtranspose = MatMultTranspose_Composite_Multiplicative; 593 b->merge_mvctx = PETSC_FALSE; 594 } else { 595 mat->ops->getdiagonal = MatGetDiagonal_Composite; 596 mat->ops->mult = MatMult_Composite; 597 mat->ops->multtranspose = MatMultTranspose_Composite; 598 } 599 PetscFunctionReturn(0); 600 } 601 602 /*@ 603 MatCompositeSetType - Indicates if the matrix is defined as the sum of a set of matrices or the product. 604 605 Logically Collective on Mat 606 607 Input Parameters: 608 . mat - the composite matrix 609 610 Level: advanced 611 612 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MatCompositeGetType(), MATCOMPOSITE 613 614 @*/ 615 PetscErrorCode MatCompositeSetType(Mat mat,MatCompositeType type) 616 { 617 PetscErrorCode ierr; 618 619 PetscFunctionBegin; 620 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 621 PetscValidLogicalCollectiveEnum(mat,type,2); 622 ierr = PetscUseMethod(mat,"MatCompositeSetType_C",(Mat,MatCompositeType),(mat,type));CHKERRQ(ierr); 623 PetscFunctionReturn(0); 624 } 625 626 static PetscErrorCode MatCompositeGetType_Composite(Mat mat,MatCompositeType *type) 627 { 628 Mat_Composite *b = (Mat_Composite*)mat->data; 629 630 PetscFunctionBegin; 631 *type = b->type; 632 PetscFunctionReturn(0); 633 } 634 635 /*@ 636 MatCompositeGetType - Returns type of composite. 637 638 Not Collective 639 640 Input Parameter: 641 . mat - the composite matrix 642 643 Output Parameter: 644 . type - type of composite 645 646 Level: advanced 647 648 .seealso: MatCreateComposite(), MatCompositeSetType(), MATCOMPOSITE 649 650 @*/ 651 PetscErrorCode MatCompositeGetType(Mat mat,MatCompositeType *type) 652 { 653 PetscErrorCode ierr; 654 655 PetscFunctionBegin; 656 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 657 PetscValidPointer(type,2); 658 ierr = PetscUseMethod(mat,"MatCompositeGetType_C",(Mat,MatCompositeType*),(mat,type));CHKERRQ(ierr); 659 PetscFunctionReturn(0); 660 } 661 662 static PetscErrorCode MatCompositeSetMatStructure_Composite(Mat mat,MatStructure str) 663 { 664 Mat_Composite *b = (Mat_Composite*)mat->data; 665 666 PetscFunctionBegin; 667 b->structure = str; 668 PetscFunctionReturn(0); 669 } 670 671 /*@ 672 MatCompositeSetMatStructure - Indicates structure of matrices in the composite matrix. 673 674 Not Collective 675 676 Input Parameters: 677 + mat - the composite matrix 678 - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN (default) or SUBSET_NONZERO_PATTERN 679 680 Level: advanced 681 682 Notes: 683 Information about the matrices structure is used in MatCompositeMerge() for additive composite matrix. 684 685 .seealso: MatAXPY(), MatCreateComposite(), MatCompositeMerge() MatCompositeGetMatStructure(), MATCOMPOSITE 686 687 @*/ 688 PetscErrorCode MatCompositeSetMatStructure(Mat mat,MatStructure str) 689 { 690 PetscErrorCode ierr; 691 692 PetscFunctionBegin; 693 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 694 ierr = PetscUseMethod(mat,"MatCompositeSetMatStructure_C",(Mat,MatStructure),(mat,str));CHKERRQ(ierr); 695 PetscFunctionReturn(0); 696 } 697 698 static PetscErrorCode MatCompositeGetMatStructure_Composite(Mat mat,MatStructure *str) 699 { 700 Mat_Composite *b = (Mat_Composite*)mat->data; 701 702 PetscFunctionBegin; 703 *str = b->structure; 704 PetscFunctionReturn(0); 705 } 706 707 /*@ 708 MatCompositeGetMatStructure - Returns the structure of matrices in the composite matrix. 709 710 Not Collective 711 712 Input Parameter: 713 . mat - the composite matrix 714 715 Output Parameter: 716 . str - structure of the matrices 717 718 Level: advanced 719 720 .seealso: MatCreateComposite(), MatCompositeSetMatStructure(), MATCOMPOSITE 721 722 @*/ 723 PetscErrorCode MatCompositeGetMatStructure(Mat mat,MatStructure *str) 724 { 725 PetscErrorCode ierr; 726 727 PetscFunctionBegin; 728 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 729 PetscValidPointer(str,2); 730 ierr = PetscUseMethod(mat,"MatCompositeGetMatStructure_C",(Mat,MatStructure*),(mat,str));CHKERRQ(ierr); 731 PetscFunctionReturn(0); 732 } 733 734 static PetscErrorCode MatCompositeSetMergeType_Composite(Mat mat,MatCompositeMergeType type) 735 { 736 Mat_Composite *shell = (Mat_Composite*)mat->data; 737 738 PetscFunctionBegin; 739 shell->mergetype = type; 740 PetscFunctionReturn(0); 741 } 742 743 /*@ 744 MatCompositeSetMergeType - Sets order of MatCompositeMerge(). 745 746 Logically Collective on Mat 747 748 Input Parameter: 749 + mat - the composite matrix 750 - type - MAT_COMPOSITE_MERGE RIGHT (default) to start merge from right with the first added matrix (mat[0]), 751 MAT_COMPOSITE_MERGE_LEFT to start merge from left with the last added matrix (mat[nmat-1]) 752 753 Level: advanced 754 755 Notes: 756 The resulting matrix is the same regardles of the MergeType. Only the order of operation is changed. 757 If set to MAT_COMPOSITE_MERGE_RIGHT the order of the merge is mat[nmat-1]*(mat[nmat-2]*(...*(mat[1]*mat[0]))) 758 otherwise the order is (((mat[nmat-1]*mat[nmat-2])*mat[nmat-3])*...)*mat[0]. 759 760 .seealso: MatCreateComposite(), MatCompositeMerge(), MATCOMPOSITE 761 762 @*/ 763 PetscErrorCode MatCompositeSetMergeType(Mat mat,MatCompositeMergeType type) 764 { 765 PetscErrorCode ierr; 766 767 PetscFunctionBegin; 768 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 769 PetscValidLogicalCollectiveEnum(mat,type,2); 770 ierr = PetscUseMethod(mat,"MatCompositeSetMergeType_C",(Mat,MatCompositeMergeType),(mat,type));CHKERRQ(ierr); 771 PetscFunctionReturn(0); 772 } 773 774 static PetscErrorCode MatCompositeMerge_Composite(Mat mat) 775 { 776 Mat_Composite *shell = (Mat_Composite*)mat->data; 777 Mat_CompositeLink next = shell->head, prev = shell->tail; 778 PetscErrorCode ierr; 779 Mat tmat,newmat; 780 Vec left,right; 781 PetscScalar scale; 782 PetscInt i; 783 784 PetscFunctionBegin; 785 if (!next) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must provide at least one matrix with MatCompositeAddMat()"); 786 scale = shell->scale; 787 if (shell->type == MAT_COMPOSITE_ADDITIVE) { 788 if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) { 789 i = 0; 790 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 791 if (shell->scalings) {ierr = MatScale(tmat,shell->scalings[i++]);CHKERRQ(ierr);} 792 while ((next = next->next)) { 793 ierr = MatAXPY(tmat,(shell->scalings ? shell->scalings[i++] : 1.0),next->mat,shell->structure);CHKERRQ(ierr); 794 } 795 } else { 796 i = shell->nmat-1; 797 ierr = MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 798 if (shell->scalings) {ierr = MatScale(tmat,shell->scalings[i--]);CHKERRQ(ierr);} 799 while ((prev = prev->prev)) { 800 ierr = MatAXPY(tmat,(shell->scalings ? shell->scalings[i--] : 1.0),prev->mat,shell->structure);CHKERRQ(ierr); 801 } 802 } 803 } else { 804 if (shell->mergetype == MAT_COMPOSITE_MERGE_RIGHT) { 805 ierr = MatDuplicate(next->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 806 while ((next = next->next)) { 807 ierr = MatMatMult(next->mat,tmat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr); 808 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 809 tmat = newmat; 810 } 811 } else { 812 ierr = MatDuplicate(prev->mat,MAT_COPY_VALUES,&tmat);CHKERRQ(ierr); 813 while ((prev = prev->prev)) { 814 ierr = MatMatMult(tmat,prev->mat,MAT_INITIAL_MATRIX,PETSC_DECIDE,&newmat);CHKERRQ(ierr); 815 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 816 tmat = newmat; 817 } 818 } 819 if (shell->scalings) {for (i=0; i<shell->nmat; i++) scale *= shell->scalings[i];} 820 } 821 822 if ((left = shell->left)) {ierr = PetscObjectReference((PetscObject)left);CHKERRQ(ierr);} 823 if ((right = shell->right)) {ierr = PetscObjectReference((PetscObject)right);CHKERRQ(ierr);} 824 825 ierr = MatHeaderReplace(mat,&tmat);CHKERRQ(ierr); 826 827 ierr = MatDiagonalScale(mat,left,right);CHKERRQ(ierr); 828 ierr = MatScale(mat,scale);CHKERRQ(ierr); 829 ierr = VecDestroy(&left);CHKERRQ(ierr); 830 ierr = VecDestroy(&right);CHKERRQ(ierr); 831 PetscFunctionReturn(0); 832 } 833 834 /*@ 835 MatCompositeMerge - Given a composite matrix, replaces it with a "regular" matrix 836 by summing or computing the product of all the matrices inside the composite matrix. 837 838 Collective 839 840 Input Parameters: 841 . mat - the composite matrix 842 843 Options Database Keys: 844 + -mat_composite_merge - merge in MatAssemblyEnd() 845 - -mat_composite_merge_type - set merge direction 846 847 Level: advanced 848 849 Notes: 850 The MatType of the resulting matrix will be the same as the MatType of the FIRST 851 matrix in the composite matrix. 852 853 .seealso: MatDestroy(), MatMult(), MatCompositeAddMat(), MatCreateComposite(), MatCompositeSetMatStructure(), MatCompositeSetMergeType(), MATCOMPOSITE 854 855 @*/ 856 PetscErrorCode MatCompositeMerge(Mat mat) 857 { 858 PetscErrorCode ierr; 859 860 PetscFunctionBegin; 861 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 862 ierr = PetscUseMethod(mat,"MatCompositeMerge_C",(Mat),(mat));CHKERRQ(ierr); 863 PetscFunctionReturn(0); 864 } 865 866 static PetscErrorCode MatCompositeGetNumberMat_Composite(Mat mat,PetscInt *nmat) 867 { 868 Mat_Composite *shell = (Mat_Composite*)mat->data; 869 870 PetscFunctionBegin; 871 *nmat = shell->nmat; 872 PetscFunctionReturn(0); 873 } 874 875 /*@ 876 MatCompositeGetNumberMat - Returns the number of matrices in the composite matrix. 877 878 Not Collective 879 880 Input Parameter: 881 . mat - the composite matrix 882 883 Output Parameter: 884 . nmat - number of matrices in the composite matrix 885 886 Level: advanced 887 888 .seealso: MatCreateComposite(), MatCompositeGetMat(), MATCOMPOSITE 889 890 @*/ 891 PetscErrorCode MatCompositeGetNumberMat(Mat mat,PetscInt *nmat) 892 { 893 PetscErrorCode ierr; 894 895 PetscFunctionBegin; 896 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 897 PetscValidPointer(nmat,2); 898 ierr = PetscUseMethod(mat,"MatCompositeGetNumberMat_C",(Mat,PetscInt*),(mat,nmat));CHKERRQ(ierr); 899 PetscFunctionReturn(0); 900 } 901 902 static PetscErrorCode MatCompositeGetMat_Composite(Mat mat,PetscInt i,Mat *Ai) 903 { 904 Mat_Composite *shell = (Mat_Composite*)mat->data; 905 Mat_CompositeLink ilink; 906 PetscInt k; 907 908 PetscFunctionBegin; 909 if (i >= shell->nmat) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_OUTOFRANGE,"index out of range: %d >= %d",i,shell->nmat); 910 ilink = shell->head; 911 for (k=0; k<i; k++) { 912 ilink = ilink->next; 913 } 914 *Ai = ilink->mat; 915 PetscFunctionReturn(0); 916 } 917 918 /*@ 919 MatCompositeGetMat - Returns the ith matrix from the composite matrix. 920 921 Logically Collective on Mat 922 923 Input Parameter: 924 + mat - the composite matrix 925 - i - the number of requested matrix 926 927 Output Parameter: 928 . Ai - ith matrix in composite 929 930 Level: advanced 931 932 .seealso: MatCreateComposite(), MatCompositeGetNumberMat(), MatCompositeAddMat(), MATCOMPOSITE 933 934 @*/ 935 PetscErrorCode MatCompositeGetMat(Mat mat,PetscInt i,Mat *Ai) 936 { 937 PetscErrorCode ierr; 938 939 PetscFunctionBegin; 940 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 941 PetscValidLogicalCollectiveInt(mat,i,2); 942 PetscValidPointer(Ai,3); 943 ierr = PetscUseMethod(mat,"MatCompositeGetMat_C",(Mat,PetscInt,Mat*),(mat,i,Ai));CHKERRQ(ierr); 944 PetscFunctionReturn(0); 945 } 946 947 PetscErrorCode MatCompositeSetScalings_Composite(Mat mat,const PetscScalar *scalings) 948 { 949 PetscErrorCode ierr; 950 Mat_Composite *shell = (Mat_Composite*)mat->data; 951 PetscInt nmat; 952 953 PetscFunctionBegin; 954 ierr = MatCompositeGetNumberMat(mat,&nmat);CHKERRQ(ierr); 955 if (!shell->scalings) {ierr = PetscMalloc1(nmat,&shell->scalings);CHKERRQ(ierr);} 956 ierr = PetscArraycpy(shell->scalings,scalings,nmat);CHKERRQ(ierr); 957 PetscFunctionReturn(0); 958 } 959 960 /*@ 961 MatCompositeSetScalings - Sets separate scaling factors for component matrices. 962 963 Logically Collective on Mat 964 965 Input Parameter: 966 + mat - the composite matrix 967 - scalings - array of scaling factors with scalings[i] being factor of i-th matrix, for i in [0, nmat) 968 969 Level: advanced 970 971 .seealso: MatScale(), MatDiagonalScale(), MATCOMPOSITE 972 973 @*/ 974 PetscErrorCode MatCompositeSetScalings(Mat mat,const PetscScalar *scalings) 975 { 976 PetscErrorCode ierr; 977 978 PetscFunctionBegin; 979 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 980 PetscValidPointer(scalings,2); 981 PetscValidLogicalCollectiveScalar(mat,*scalings,2); 982 ierr = PetscUseMethod(mat,"MatCompositeSetScalings_C",(Mat,const PetscScalar*),(mat,scalings));CHKERRQ(ierr); 983 PetscFunctionReturn(0); 984 } 985 986 static struct _MatOps MatOps_Values = {NULL, 987 NULL, 988 NULL, 989 MatMult_Composite, 990 MatMultAdd_Composite, 991 /* 5*/ MatMultTranspose_Composite, 992 MatMultTransposeAdd_Composite, 993 NULL, 994 NULL, 995 NULL, 996 /* 10*/ NULL, 997 NULL, 998 NULL, 999 NULL, 1000 NULL, 1001 /* 15*/ NULL, 1002 NULL, 1003 MatGetDiagonal_Composite, 1004 MatDiagonalScale_Composite, 1005 NULL, 1006 /* 20*/ NULL, 1007 MatAssemblyEnd_Composite, 1008 NULL, 1009 NULL, 1010 /* 24*/ NULL, 1011 NULL, 1012 NULL, 1013 NULL, 1014 NULL, 1015 /* 29*/ NULL, 1016 NULL, 1017 NULL, 1018 NULL, 1019 NULL, 1020 /* 34*/ NULL, 1021 NULL, 1022 NULL, 1023 NULL, 1024 NULL, 1025 /* 39*/ NULL, 1026 NULL, 1027 NULL, 1028 NULL, 1029 NULL, 1030 /* 44*/ NULL, 1031 MatScale_Composite, 1032 MatShift_Basic, 1033 NULL, 1034 NULL, 1035 /* 49*/ NULL, 1036 NULL, 1037 NULL, 1038 NULL, 1039 NULL, 1040 /* 54*/ NULL, 1041 NULL, 1042 NULL, 1043 NULL, 1044 NULL, 1045 /* 59*/ NULL, 1046 MatDestroy_Composite, 1047 NULL, 1048 NULL, 1049 NULL, 1050 /* 64*/ NULL, 1051 NULL, 1052 NULL, 1053 NULL, 1054 NULL, 1055 /* 69*/ NULL, 1056 NULL, 1057 NULL, 1058 NULL, 1059 NULL, 1060 /* 74*/ NULL, 1061 NULL, 1062 MatSetFromOptions_Composite, 1063 NULL, 1064 NULL, 1065 /* 79*/ NULL, 1066 NULL, 1067 NULL, 1068 NULL, 1069 NULL, 1070 /* 84*/ NULL, 1071 NULL, 1072 NULL, 1073 NULL, 1074 NULL, 1075 /* 89*/ NULL, 1076 NULL, 1077 NULL, 1078 NULL, 1079 NULL, 1080 /* 94*/ NULL, 1081 NULL, 1082 NULL, 1083 NULL, 1084 NULL, 1085 /*99*/ NULL, 1086 NULL, 1087 NULL, 1088 NULL, 1089 NULL, 1090 /*104*/ NULL, 1091 NULL, 1092 NULL, 1093 NULL, 1094 NULL, 1095 /*109*/ NULL, 1096 NULL, 1097 NULL, 1098 NULL, 1099 NULL, 1100 /*114*/ NULL, 1101 NULL, 1102 NULL, 1103 NULL, 1104 NULL, 1105 /*119*/ NULL, 1106 NULL, 1107 NULL, 1108 NULL, 1109 NULL, 1110 /*124*/ NULL, 1111 NULL, 1112 NULL, 1113 NULL, 1114 NULL, 1115 /*129*/ NULL, 1116 NULL, 1117 NULL, 1118 NULL, 1119 NULL, 1120 /*134*/ NULL, 1121 NULL, 1122 NULL, 1123 NULL, 1124 NULL, 1125 /*139*/ NULL, 1126 NULL, 1127 NULL 1128 }; 1129 1130 /*MC 1131 MATCOMPOSITE - A matrix defined by the sum (or product) of one or more matrices. 1132 The matrices need to have a correct size and parallel layout for the sum or product to be valid. 1133 1134 Notes: 1135 to use the product of the matrices call MatCompositeSetType(mat,MAT_COMPOSITE_MULTIPLICATIVE); 1136 1137 Level: advanced 1138 1139 .seealso: MatCreateComposite(), MatCompositeSetScalings(), MatCompositeAddMat(), MatSetType(), MatCompositeSetType(), MatCompositeGetType(), MatCompositeSetMatStructure(), MatCompositeGetMatStructure(), MatCompositeMerge(), MatCompositeSetMergeType(), MatCompositeGetNumberMat(), MatCompositeGetMat() 1140 M*/ 1141 1142 PETSC_EXTERN PetscErrorCode MatCreate_Composite(Mat A) 1143 { 1144 Mat_Composite *b; 1145 PetscErrorCode ierr; 1146 1147 PetscFunctionBegin; 1148 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 1149 A->data = (void*)b; 1150 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1151 1152 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1153 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1154 1155 A->assembled = PETSC_TRUE; 1156 A->preallocated = PETSC_TRUE; 1157 b->type = MAT_COMPOSITE_ADDITIVE; 1158 b->scale = 1.0; 1159 b->nmat = 0; 1160 b->merge = PETSC_FALSE; 1161 b->mergetype = MAT_COMPOSITE_MERGE_RIGHT; 1162 b->structure = DIFFERENT_NONZERO_PATTERN; 1163 b->merge_mvctx = PETSC_TRUE; 1164 1165 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeAddMat_C",MatCompositeAddMat_Composite);CHKERRQ(ierr); 1166 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetType_C",MatCompositeSetType_Composite);CHKERRQ(ierr); 1167 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetType_C",MatCompositeGetType_Composite);CHKERRQ(ierr); 1168 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetMergeType_C",MatCompositeSetMergeType_Composite);CHKERRQ(ierr); 1169 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetMatStructure_C",MatCompositeSetMatStructure_Composite);CHKERRQ(ierr); 1170 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMatStructure_C",MatCompositeGetMatStructure_Composite);CHKERRQ(ierr); 1171 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeMerge_C",MatCompositeMerge_Composite);CHKERRQ(ierr); 1172 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetNumberMat_C",MatCompositeGetNumberMat_Composite);CHKERRQ(ierr); 1173 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeGetMat_C",MatCompositeGetMat_Composite);CHKERRQ(ierr); 1174 ierr = PetscObjectComposeFunction((PetscObject)A,"MatCompositeSetScalings_C",MatCompositeSetScalings_Composite);CHKERRQ(ierr); 1175 1176 ierr = PetscObjectChangeTypeName((PetscObject)A,MATCOMPOSITE);CHKERRQ(ierr); 1177 PetscFunctionReturn(0); 1178 } 1179 1180