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