1 /* 2 Creates a matrix class for using the Neumann-Neumann type preconditioners. 3 This stores the matrices in globally unassembled form. Each processor 4 assembles only its local Neumann problem and the parallel matrix vector 5 product is handled "implicitly". 6 7 Currently this allows for only one subdomain per processor. 8 */ 9 10 #include <../src/mat/impls/is/matis.h> /*I "petscmat.h" I*/ 11 #include <petsc/private/sfimpl.h> 12 #include <petsc/private/vecimpl.h> 13 #include <petsc/private/hashseti.h> 14 15 #define MATIS_MAX_ENTRIES_INSERTION 2048 16 static PetscErrorCode MatSetValuesLocal_IS(Mat,PetscInt,const PetscInt*,PetscInt,const PetscInt*,const PetscScalar*,InsertMode); 17 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat,PetscInt,const PetscInt*,PetscInt,const PetscInt*,const PetscScalar*,InsertMode); 18 static PetscErrorCode MatISSetUpScatters_Private(Mat); 19 20 static PetscErrorCode MatISContainerDestroyPtAP_Private(void *ptr) 21 { 22 MatISPtAP ptap = (MatISPtAP)ptr; 23 24 PetscFunctionBegin; 25 PetscCall(MatDestroySubMatrices(ptap->ris1 ? 2 : 1,&ptap->lP)); 26 PetscCall(ISDestroy(&ptap->cis0)); 27 PetscCall(ISDestroy(&ptap->cis1)); 28 PetscCall(ISDestroy(&ptap->ris0)); 29 PetscCall(ISDestroy(&ptap->ris1)); 30 PetscCall(PetscFree(ptap)); 31 PetscFunctionReturn(0); 32 } 33 34 static PetscErrorCode MatPtAPNumeric_IS_XAIJ(Mat A, Mat P, Mat C) 35 { 36 MatISPtAP ptap; 37 Mat_IS *matis = (Mat_IS*)A->data; 38 Mat lA,lC; 39 MatReuse reuse; 40 IS ris[2],cis[2]; 41 PetscContainer c; 42 PetscInt n; 43 44 PetscFunctionBegin; 45 PetscCall(PetscObjectQuery((PetscObject)C,"_MatIS_PtAP",(PetscObject*)&c)); 46 PetscCheck(c,PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing PtAP information"); 47 PetscCall(PetscContainerGetPointer(c,(void**)&ptap)); 48 ris[0] = ptap->ris0; 49 ris[1] = ptap->ris1; 50 cis[0] = ptap->cis0; 51 cis[1] = ptap->cis1; 52 n = ptap->ris1 ? 2 : 1; 53 reuse = ptap->lP ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX; 54 PetscCall(MatCreateSubMatrices(P,n,ris,cis,reuse,&ptap->lP)); 55 56 PetscCall(MatISGetLocalMat(A,&lA)); 57 PetscCall(MatISGetLocalMat(C,&lC)); 58 if (ptap->ris1) { /* unsymmetric A mapping */ 59 Mat lPt; 60 61 PetscCall(MatTranspose(ptap->lP[1],MAT_INITIAL_MATRIX,&lPt)); 62 PetscCall(MatMatMatMult(lPt,lA,ptap->lP[0],reuse,ptap->fill,&lC)); 63 if (matis->storel2l) { 64 PetscCall(PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",(PetscObject)lPt)); 65 } 66 PetscCall(MatDestroy(&lPt)); 67 } else { 68 PetscCall(MatPtAP(lA,ptap->lP[0],reuse,ptap->fill,&lC)); 69 if (matis->storel2l) { 70 PetscCall(PetscObjectCompose((PetscObject)C,"_MatIS_PtAP_l2l",(PetscObject)ptap->lP[0])); 71 } 72 } 73 if (reuse == MAT_INITIAL_MATRIX) { 74 PetscCall(MatISSetLocalMat(C,lC)); 75 PetscCall(MatDestroy(&lC)); 76 } 77 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 78 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 79 PetscFunctionReturn(0); 80 } 81 82 static PetscErrorCode MatGetNonzeroColumnsLocal_Private(Mat PT,IS *cis) 83 { 84 Mat Po,Pd; 85 IS zd,zo; 86 const PetscInt *garray; 87 PetscInt *aux,i,bs; 88 PetscInt dc,stc,oc,ctd,cto; 89 PetscBool ismpiaij,ismpibaij,isseqaij,isseqbaij; 90 MPI_Comm comm; 91 92 PetscFunctionBegin; 93 PetscValidHeaderSpecific(PT,MAT_CLASSID,1); 94 PetscValidPointer(cis,2); 95 PetscCall(PetscObjectGetComm((PetscObject)PT,&comm)); 96 bs = 1; 97 PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT,MATMPIAIJ,&ismpiaij)); 98 PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT,MATMPIBAIJ,&ismpibaij)); 99 PetscCall(PetscObjectBaseTypeCompare((PetscObject)PT,MATSEQAIJ,&isseqaij)); 100 PetscCall(PetscObjectTypeCompare((PetscObject)PT,MATSEQBAIJ,&isseqbaij)); 101 if (isseqaij || isseqbaij) { 102 Pd = PT; 103 Po = NULL; 104 garray = NULL; 105 } else if (ismpiaij) { 106 PetscCall(MatMPIAIJGetSeqAIJ(PT,&Pd,&Po,&garray)); 107 } else if (ismpibaij) { 108 PetscCall(MatMPIBAIJGetSeqBAIJ(PT,&Pd,&Po,&garray)); 109 PetscCall(MatGetBlockSize(PT,&bs)); 110 } else SETERRQ(comm,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(PT))->type_name); 111 112 /* identify any null columns in Pd or Po */ 113 /* We use a tolerance comparison since it may happen that, with geometric multigrid, 114 some of the columns are not really zero, but very close to */ 115 zo = zd = NULL; 116 if (Po) { 117 PetscCall(MatFindNonzeroRowsOrCols_Basic(Po,PETSC_TRUE,PETSC_SMALL,&zo)); 118 } 119 PetscCall(MatFindNonzeroRowsOrCols_Basic(Pd,PETSC_TRUE,PETSC_SMALL,&zd)); 120 121 PetscCall(MatGetLocalSize(PT,NULL,&dc)); 122 PetscCall(MatGetOwnershipRangeColumn(PT,&stc,NULL)); 123 if (Po) PetscCall(MatGetLocalSize(Po,NULL,&oc)); 124 else oc = 0; 125 PetscCall(PetscMalloc1((dc+oc)/bs,&aux)); 126 if (zd) { 127 const PetscInt *idxs; 128 PetscInt nz; 129 130 /* this will throw an error if bs is not valid */ 131 PetscCall(ISSetBlockSize(zd,bs)); 132 PetscCall(ISGetLocalSize(zd,&nz)); 133 PetscCall(ISGetIndices(zd,&idxs)); 134 ctd = nz/bs; 135 for (i=0; i<ctd; i++) aux[i] = (idxs[bs*i]+stc)/bs; 136 PetscCall(ISRestoreIndices(zd,&idxs)); 137 } else { 138 ctd = dc/bs; 139 for (i=0; i<ctd; i++) aux[i] = i+stc/bs; 140 } 141 if (zo) { 142 const PetscInt *idxs; 143 PetscInt nz; 144 145 /* this will throw an error if bs is not valid */ 146 PetscCall(ISSetBlockSize(zo,bs)); 147 PetscCall(ISGetLocalSize(zo,&nz)); 148 PetscCall(ISGetIndices(zo,&idxs)); 149 cto = nz/bs; 150 for (i=0; i<cto; i++) aux[i+ctd] = garray[idxs[bs*i]/bs]; 151 PetscCall(ISRestoreIndices(zo,&idxs)); 152 } else { 153 cto = oc/bs; 154 for (i=0; i<cto; i++) aux[i+ctd] = garray[i]; 155 } 156 PetscCall(ISCreateBlock(comm,bs,ctd+cto,aux,PETSC_OWN_POINTER,cis)); 157 PetscCall(ISDestroy(&zd)); 158 PetscCall(ISDestroy(&zo)); 159 PetscFunctionReturn(0); 160 } 161 162 static PetscErrorCode MatPtAPSymbolic_IS_XAIJ(Mat A,Mat P,PetscReal fill,Mat C) 163 { 164 Mat PT,lA; 165 MatISPtAP ptap; 166 ISLocalToGlobalMapping Crl2g,Ccl2g,rl2g,cl2g; 167 PetscContainer c; 168 MatType lmtype; 169 const PetscInt *garray; 170 PetscInt ibs,N,dc; 171 MPI_Comm comm; 172 173 PetscFunctionBegin; 174 PetscCall(PetscObjectGetComm((PetscObject)A,&comm)); 175 PetscCall(MatSetType(C,MATIS)); 176 PetscCall(MatISGetLocalMat(A,&lA)); 177 PetscCall(MatGetType(lA,&lmtype)); 178 PetscCall(MatISSetLocalMatType(C,lmtype)); 179 PetscCall(MatGetSize(P,NULL,&N)); 180 PetscCall(MatGetLocalSize(P,NULL,&dc)); 181 PetscCall(MatSetSizes(C,dc,dc,N,N)); 182 /* Not sure about this 183 PetscCall(MatGetBlockSizes(P,NULL,&ibs)); 184 PetscCall(MatSetBlockSize(*C,ibs)); 185 */ 186 187 PetscCall(PetscNew(&ptap)); 188 PetscCall(PetscContainerCreate(PETSC_COMM_SELF,&c)); 189 PetscCall(PetscContainerSetPointer(c,ptap)); 190 PetscCall(PetscContainerSetUserDestroy(c,MatISContainerDestroyPtAP_Private)); 191 PetscCall(PetscObjectCompose((PetscObject)C,"_MatIS_PtAP",(PetscObject)c)); 192 PetscCall(PetscContainerDestroy(&c)); 193 ptap->fill = fill; 194 195 PetscCall(MatISGetLocalToGlobalMapping(A,&rl2g,&cl2g)); 196 197 PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g,&ibs)); 198 PetscCall(ISLocalToGlobalMappingGetSize(cl2g,&N)); 199 PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g,&garray)); 200 PetscCall(ISCreateBlock(comm,ibs,N/ibs,garray,PETSC_COPY_VALUES,&ptap->ris0)); 201 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g,&garray)); 202 203 PetscCall(MatCreateSubMatrix(P,ptap->ris0,NULL,MAT_INITIAL_MATRIX,&PT)); 204 PetscCall(MatGetNonzeroColumnsLocal_Private(PT,&ptap->cis0)); 205 PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis0,&Ccl2g)); 206 PetscCall(MatDestroy(&PT)); 207 208 Crl2g = NULL; 209 if (rl2g != cl2g) { /* unsymmetric A mapping */ 210 PetscBool same,lsame = PETSC_FALSE; 211 PetscInt N1,ibs1; 212 213 PetscCall(ISLocalToGlobalMappingGetSize(rl2g,&N1)); 214 PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g,&ibs1)); 215 PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g,&garray)); 216 PetscCall(ISCreateBlock(comm,ibs,N/ibs,garray,PETSC_COPY_VALUES,&ptap->ris1)); 217 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g,&garray)); 218 if (ibs1 == ibs && N1 == N) { /* check if the l2gmaps are the same */ 219 const PetscInt *i1,*i2; 220 221 PetscCall(ISBlockGetIndices(ptap->ris0,&i1)); 222 PetscCall(ISBlockGetIndices(ptap->ris1,&i2)); 223 PetscCall(PetscArraycmp(i1,i2,N,&lsame)); 224 } 225 PetscCall(MPIU_Allreduce(&lsame,&same,1,MPIU_BOOL,MPI_LAND,comm)); 226 if (same) { 227 PetscCall(ISDestroy(&ptap->ris1)); 228 } else { 229 PetscCall(MatCreateSubMatrix(P,ptap->ris1,NULL,MAT_INITIAL_MATRIX,&PT)); 230 PetscCall(MatGetNonzeroColumnsLocal_Private(PT,&ptap->cis1)); 231 PetscCall(ISLocalToGlobalMappingCreateIS(ptap->cis1,&Crl2g)); 232 PetscCall(MatDestroy(&PT)); 233 } 234 } 235 /* Not sure about this 236 if (!Crl2g) { 237 PetscCall(MatGetBlockSize(C,&ibs)); 238 PetscCall(ISLocalToGlobalMappingSetBlockSize(Ccl2g,ibs)); 239 } 240 */ 241 PetscCall(MatSetLocalToGlobalMapping(C,Crl2g ? Crl2g : Ccl2g,Ccl2g)); 242 PetscCall(ISLocalToGlobalMappingDestroy(&Crl2g)); 243 PetscCall(ISLocalToGlobalMappingDestroy(&Ccl2g)); 244 245 C->ops->ptapnumeric = MatPtAPNumeric_IS_XAIJ; 246 PetscFunctionReturn(0); 247 } 248 249 /* ----------------------------------------- */ 250 static PetscErrorCode MatProductSymbolic_PtAP_IS_XAIJ(Mat C) 251 { 252 Mat_Product *product = C->product; 253 Mat A=product->A,P=product->B; 254 PetscReal fill=product->fill; 255 256 PetscFunctionBegin; 257 PetscCall(MatPtAPSymbolic_IS_XAIJ(A,P,fill,C)); 258 C->ops->productnumeric = MatProductNumeric_PtAP; 259 PetscFunctionReturn(0); 260 } 261 262 static PetscErrorCode MatProductSetFromOptions_IS_XAIJ_PtAP(Mat C) 263 { 264 PetscFunctionBegin; 265 C->ops->productsymbolic = MatProductSymbolic_PtAP_IS_XAIJ; 266 PetscFunctionReturn(0); 267 } 268 269 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat C) 270 { 271 Mat_Product *product = C->product; 272 273 PetscFunctionBegin; 274 if (product->type == MATPRODUCT_PtAP) { 275 PetscCall(MatProductSetFromOptions_IS_XAIJ_PtAP(C)); 276 } 277 PetscFunctionReturn(0); 278 } 279 280 /* ----------------------------------------- */ 281 static PetscErrorCode MatISContainerDestroyFields_Private(void *ptr) 282 { 283 MatISLocalFields lf = (MatISLocalFields)ptr; 284 PetscInt i; 285 286 PetscFunctionBegin; 287 for (i=0;i<lf->nr;i++) { 288 PetscCall(ISDestroy(&lf->rf[i])); 289 } 290 for (i=0;i<lf->nc;i++) { 291 PetscCall(ISDestroy(&lf->cf[i])); 292 } 293 PetscCall(PetscFree2(lf->rf,lf->cf)); 294 PetscCall(PetscFree(lf)); 295 PetscFunctionReturn(0); 296 } 297 298 static PetscErrorCode MatConvert_SeqXAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat) 299 { 300 Mat B,lB; 301 302 PetscFunctionBegin; 303 if (reuse != MAT_REUSE_MATRIX) { 304 ISLocalToGlobalMapping rl2g,cl2g; 305 PetscInt bs; 306 IS is; 307 308 PetscCall(MatGetBlockSize(A,&bs)); 309 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),A->rmap->n/bs,0,1,&is)); 310 if (bs > 1) { 311 IS is2; 312 PetscInt i,*aux; 313 314 PetscCall(ISGetLocalSize(is,&i)); 315 PetscCall(ISGetIndices(is,(const PetscInt**)&aux)); 316 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A),bs,i,aux,PETSC_COPY_VALUES,&is2)); 317 PetscCall(ISRestoreIndices(is,(const PetscInt**)&aux)); 318 PetscCall(ISDestroy(&is)); 319 is = is2; 320 } 321 PetscCall(ISSetIdentity(is)); 322 PetscCall(ISLocalToGlobalMappingCreateIS(is,&rl2g)); 323 PetscCall(ISDestroy(&is)); 324 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),A->cmap->n/bs,0,1,&is)); 325 if (bs > 1) { 326 IS is2; 327 PetscInt i,*aux; 328 329 PetscCall(ISGetLocalSize(is,&i)); 330 PetscCall(ISGetIndices(is,(const PetscInt**)&aux)); 331 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A),bs,i,aux,PETSC_COPY_VALUES,&is2)); 332 PetscCall(ISRestoreIndices(is,(const PetscInt**)&aux)); 333 PetscCall(ISDestroy(&is)); 334 is = is2; 335 } 336 PetscCall(ISSetIdentity(is)); 337 PetscCall(ISLocalToGlobalMappingCreateIS(is,&cl2g)); 338 PetscCall(ISDestroy(&is)); 339 PetscCall(MatCreateIS(PetscObjectComm((PetscObject)A),bs,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N,rl2g,cl2g,&B)); 340 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 341 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 342 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&lB)); 343 if (reuse == MAT_INITIAL_MATRIX) *newmat = B; 344 } else { 345 B = *newmat; 346 PetscCall(PetscObjectReference((PetscObject)A)); 347 lB = A; 348 } 349 PetscCall(MatISSetLocalMat(B,lB)); 350 PetscCall(MatDestroy(&lB)); 351 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 352 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 353 if (reuse == MAT_INPLACE_MATRIX) { 354 PetscCall(MatHeaderReplace(A,&B)); 355 } 356 PetscFunctionReturn(0); 357 } 358 359 static PetscErrorCode MatISScaleDisassembling_Private(Mat A) 360 { 361 Mat_IS *matis = (Mat_IS*)(A->data); 362 PetscScalar *aa; 363 const PetscInt *ii,*jj; 364 PetscInt i,n,m; 365 PetscInt *ecount,**eneighs; 366 PetscBool flg; 367 368 PetscFunctionBegin; 369 PetscCall(MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&m,&ii,&jj,&flg)); 370 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure"); 371 PetscCall(ISLocalToGlobalMappingGetNodeInfo(matis->rmapping,&n,&ecount,&eneighs)); 372 PetscCheck(m == n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected %" PetscInt_FMT " != %" PetscInt_FMT,m,n); 373 PetscCall(MatSeqAIJGetArray(matis->A,&aa)); 374 for (i=0;i<n;i++) { 375 if (ecount[i] > 1) { 376 PetscInt j; 377 378 for (j=ii[i];j<ii[i+1];j++) { 379 PetscInt i2 = jj[j],p,p2; 380 PetscReal scal = 0.0; 381 382 for (p=0;p<ecount[i];p++) { 383 for (p2=0;p2<ecount[i2];p2++) { 384 if (eneighs[i][p] == eneighs[i2][p2]) { scal += 1.0; break; } 385 } 386 } 387 if (scal) aa[j] /= scal; 388 } 389 } 390 } 391 PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(matis->rmapping,&n,&ecount,&eneighs)); 392 PetscCall(MatSeqAIJRestoreArray(matis->A,&aa)); 393 PetscCall(MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&m,&ii,&jj,&flg)); 394 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure"); 395 PetscFunctionReturn(0); 396 } 397 398 typedef enum {MAT_IS_DISASSEMBLE_L2G_NATURAL,MAT_IS_DISASSEMBLE_L2G_MAT, MAT_IS_DISASSEMBLE_L2G_ND} MatISDisassemblel2gType; 399 400 static PetscErrorCode MatMPIXAIJComputeLocalToGlobalMapping_Private(Mat A, ISLocalToGlobalMapping *l2g) 401 { 402 Mat Ad,Ao; 403 IS is,ndmap,ndsub; 404 MPI_Comm comm; 405 const PetscInt *garray,*ndmapi; 406 PetscInt bs,i,cnt,nl,*ncount,*ndmapc; 407 PetscBool ismpiaij,ismpibaij; 408 const char *const MatISDisassemblel2gTypes[] = {"NATURAL","MAT","ND","MatISDisassemblel2gType","MAT_IS_DISASSEMBLE_L2G_",NULL}; 409 MatISDisassemblel2gType mode = MAT_IS_DISASSEMBLE_L2G_NATURAL; 410 MatPartitioning part; 411 PetscSF sf; 412 PetscObject dm; 413 414 PetscFunctionBegin; 415 PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatIS l2g disassembling options","Mat"); 416 PetscCall(PetscOptionsEnum("-mat_is_disassemble_l2g_type","Type of local-to-global mapping to be used for disassembling","MatISDisassemblel2gType",MatISDisassemblel2gTypes,(PetscEnum)mode,(PetscEnum*)&mode,NULL)); 417 PetscOptionsEnd(); 418 if (mode == MAT_IS_DISASSEMBLE_L2G_MAT) { 419 PetscCall(MatGetLocalToGlobalMapping(A,l2g,NULL)); 420 PetscFunctionReturn(0); 421 } 422 PetscCall(PetscObjectGetComm((PetscObject)A,&comm)); 423 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ ,&ismpiaij)); 424 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A,MATMPIBAIJ,&ismpibaij)); 425 PetscCall(MatGetBlockSize(A,&bs)); 426 switch (mode) { 427 case MAT_IS_DISASSEMBLE_L2G_ND: 428 PetscCall(MatPartitioningCreate(comm,&part)); 429 PetscCall(MatPartitioningSetAdjacency(part,A)); 430 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)part,((PetscObject)A)->prefix)); 431 PetscCall(MatPartitioningSetFromOptions(part)); 432 PetscCall(MatPartitioningApplyND(part,&ndmap)); 433 PetscCall(MatPartitioningDestroy(&part)); 434 PetscCall(ISBuildTwoSided(ndmap,NULL,&ndsub)); 435 PetscCall(MatMPIAIJSetUseScalableIncreaseOverlap(A,PETSC_TRUE)); 436 PetscCall(MatIncreaseOverlap(A,1,&ndsub,1)); 437 PetscCall(ISLocalToGlobalMappingCreateIS(ndsub,l2g)); 438 439 /* it may happen that a separator node is not properly shared */ 440 PetscCall(ISLocalToGlobalMappingGetNodeInfo(*l2g,&nl,&ncount,NULL)); 441 PetscCall(PetscSFCreate(comm,&sf)); 442 PetscCall(ISLocalToGlobalMappingGetIndices(*l2g,&garray)); 443 PetscCall(PetscSFSetGraphLayout(sf,A->rmap,nl,NULL,PETSC_OWN_POINTER,garray)); 444 PetscCall(ISLocalToGlobalMappingRestoreIndices(*l2g,&garray)); 445 PetscCall(PetscCalloc1(A->rmap->n,&ndmapc)); 446 PetscCall(PetscSFReduceBegin(sf,MPIU_INT,ncount,ndmapc,MPI_REPLACE)); 447 PetscCall(PetscSFReduceEnd(sf,MPIU_INT,ncount,ndmapc,MPI_REPLACE)); 448 PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(*l2g,NULL,&ncount,NULL)); 449 PetscCall(ISGetIndices(ndmap,&ndmapi)); 450 for (i = 0, cnt = 0; i < A->rmap->n; i++) 451 if (ndmapi[i] < 0 && ndmapc[i] < 2) 452 cnt++; 453 454 PetscCall(MPIU_Allreduce(&cnt,&i,1,MPIU_INT,MPI_MAX,comm)); 455 if (i) { /* we detected isolated separator nodes */ 456 Mat A2,A3; 457 IS *workis,is2; 458 PetscScalar *vals; 459 PetscInt gcnt = i,*dnz,*onz,j,*lndmapi; 460 ISLocalToGlobalMapping ll2g; 461 PetscBool flg; 462 const PetscInt *ii,*jj; 463 464 /* communicate global id of separators */ 465 MatPreallocateBegin(comm,A->rmap->n,A->cmap->n,dnz,onz); 466 for (i = 0, cnt = 0; i < A->rmap->n; i++) 467 dnz[i] = ndmapi[i] < 0 ? i + A->rmap->rstart : -1; 468 469 PetscCall(PetscMalloc1(nl,&lndmapi)); 470 PetscCall(PetscSFBcastBegin(sf,MPIU_INT,dnz,lndmapi,MPI_REPLACE)); 471 472 /* compute adjacency of isolated separators node */ 473 PetscCall(PetscMalloc1(gcnt,&workis)); 474 for (i = 0, cnt = 0; i < A->rmap->n; i++) { 475 if (ndmapi[i] < 0 && ndmapc[i] < 2) { 476 PetscCall(ISCreateStride(comm,1,i+A->rmap->rstart,1,&workis[cnt++])); 477 } 478 } 479 for (i = cnt; i < gcnt; i++) { 480 PetscCall(ISCreateStride(comm,0,0,1,&workis[i])); 481 } 482 for (i = 0; i < gcnt; i++) { 483 PetscCall(PetscObjectSetName((PetscObject)workis[i],"ISOLATED")); 484 PetscCall(ISViewFromOptions(workis[i],NULL,"-view_isolated_separators")); 485 } 486 487 /* no communications since all the ISes correspond to locally owned rows */ 488 PetscCall(MatIncreaseOverlap(A,gcnt,workis,1)); 489 490 /* end communicate global id of separators */ 491 PetscCall(PetscSFBcastEnd(sf,MPIU_INT,dnz,lndmapi,MPI_REPLACE)); 492 493 /* communicate new layers : create a matrix and transpose it */ 494 PetscCall(PetscArrayzero(dnz,A->rmap->n)); 495 PetscCall(PetscArrayzero(onz,A->rmap->n)); 496 for (i = 0, j = 0; i < A->rmap->n; i++) { 497 if (ndmapi[i] < 0 && ndmapc[i] < 2) { 498 const PetscInt* idxs; 499 PetscInt s; 500 501 PetscCall(ISGetLocalSize(workis[j],&s)); 502 PetscCall(ISGetIndices(workis[j],&idxs)); 503 PetscCall(MatPreallocateSet(i+A->rmap->rstart,s,idxs,dnz,onz)); 504 j++; 505 } 506 } 507 PetscCheck(j == cnt,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT,j,cnt); 508 509 for (i = 0; i < gcnt; i++) { 510 PetscCall(PetscObjectSetName((PetscObject)workis[i],"EXTENDED")); 511 PetscCall(ISViewFromOptions(workis[i],NULL,"-view_isolated_separators")); 512 } 513 514 for (i = 0, j = 0; i < A->rmap->n; i++) j = PetscMax(j,dnz[i]+onz[i]); 515 PetscCall(PetscMalloc1(j,&vals)); 516 for (i = 0; i < j; i++) vals[i] = 1.0; 517 518 PetscCall(MatCreate(comm,&A2)); 519 PetscCall(MatSetType(A2,MATMPIAIJ)); 520 PetscCall(MatSetSizes(A2,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 521 PetscCall(MatMPIAIJSetPreallocation(A2,0,dnz,0,onz)); 522 PetscCall(MatSetOption(A2,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE)); 523 for (i = 0, j = 0; i < A2->rmap->n; i++) { 524 PetscInt row = i+A2->rmap->rstart,s = dnz[i] + onz[i]; 525 const PetscInt* idxs; 526 527 if (s) { 528 PetscCall(ISGetIndices(workis[j],&idxs)); 529 PetscCall(MatSetValues(A2,1,&row,s,idxs,vals,INSERT_VALUES)); 530 PetscCall(ISRestoreIndices(workis[j],&idxs)); 531 j++; 532 } 533 } 534 PetscCheck(j == cnt,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected local count %" PetscInt_FMT " != %" PetscInt_FMT,j,cnt); 535 PetscCall(PetscFree(vals)); 536 PetscCall(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 537 PetscCall(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 538 PetscCall(MatTranspose(A2,MAT_INPLACE_MATRIX,&A2)); 539 540 /* extract submatrix corresponding to the coupling "owned separators" x "isolated separators" */ 541 for (i = 0, j = 0; i < nl; i++) 542 if (lndmapi[i] >= 0) lndmapi[j++] = lndmapi[i]; 543 PetscCall(ISCreateGeneral(comm,j,lndmapi,PETSC_USE_POINTER,&is)); 544 PetscCall(MatMPIAIJGetLocalMatCondensed(A2,MAT_INITIAL_MATRIX,&is,NULL,&A3)); 545 PetscCall(ISDestroy(&is)); 546 PetscCall(MatDestroy(&A2)); 547 548 /* extend local to global map to include connected isolated separators */ 549 PetscCall(PetscObjectQuery((PetscObject)A3,"_petsc_GetLocalMatCondensed_iscol",(PetscObject*)&is)); 550 PetscCheck(is,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing column map"); 551 PetscCall(ISLocalToGlobalMappingCreateIS(is,&ll2g)); 552 PetscCall(MatGetRowIJ(A3,0,PETSC_FALSE,PETSC_FALSE,&i,&ii,&jj,&flg)); 553 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure"); 554 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,ii[i],jj,PETSC_COPY_VALUES,&is)); 555 PetscCall(MatRestoreRowIJ(A3,0,PETSC_FALSE,PETSC_FALSE,&i,&ii,&jj,&flg)); 556 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure"); 557 PetscCall(ISLocalToGlobalMappingApplyIS(ll2g,is,&is2)); 558 PetscCall(ISDestroy(&is)); 559 PetscCall(ISLocalToGlobalMappingDestroy(&ll2g)); 560 561 /* add new nodes to the local-to-global map */ 562 PetscCall(ISLocalToGlobalMappingDestroy(l2g)); 563 PetscCall(ISExpand(ndsub,is2,&is)); 564 PetscCall(ISDestroy(&is2)); 565 PetscCall(ISLocalToGlobalMappingCreateIS(is,l2g)); 566 PetscCall(ISDestroy(&is)); 567 568 PetscCall(MatDestroy(&A3)); 569 PetscCall(PetscFree(lndmapi)); 570 MatPreallocateEnd(dnz,onz); 571 for (i = 0; i < gcnt; i++) { 572 PetscCall(ISDestroy(&workis[i])); 573 } 574 PetscCall(PetscFree(workis)); 575 } 576 PetscCall(ISRestoreIndices(ndmap,&ndmapi)); 577 PetscCall(PetscSFDestroy(&sf)); 578 PetscCall(PetscFree(ndmapc)); 579 PetscCall(ISDestroy(&ndmap)); 580 PetscCall(ISDestroy(&ndsub)); 581 PetscCall(ISLocalToGlobalMappingSetBlockSize(*l2g,bs)); 582 PetscCall(ISLocalToGlobalMappingViewFromOptions(*l2g,NULL,"-matis_nd_l2g_view")); 583 break; 584 case MAT_IS_DISASSEMBLE_L2G_NATURAL: 585 PetscCall(PetscObjectQuery((PetscObject)A,"__PETSc_dm",(PetscObject*)&dm)); 586 if (dm) { /* if a matrix comes from a DM, most likely we can use the l2gmap if any */ 587 PetscCall(MatGetLocalToGlobalMapping(A,l2g,NULL)); 588 PetscCall(PetscObjectReference((PetscObject)*l2g)); 589 if (*l2g) PetscFunctionReturn(0); 590 } 591 if (ismpiaij) { 592 PetscCall(MatMPIAIJGetSeqAIJ(A,&Ad,&Ao,&garray)); 593 } else if (ismpibaij) { 594 PetscCall(MatMPIBAIJGetSeqBAIJ(A,&Ad,&Ao,&garray)); 595 } else SETERRQ(comm,PETSC_ERR_SUP,"Type %s",((PetscObject)A)->type_name); 596 PetscCheck(garray,comm,PETSC_ERR_ARG_WRONGSTATE,"garray not present"); 597 if (A->rmap->n) { 598 PetscInt dc,oc,stc,*aux; 599 600 PetscCall(MatGetLocalSize(Ad,NULL,&dc)); 601 PetscCall(MatGetLocalSize(Ao,NULL,&oc)); 602 PetscCall(MatGetOwnershipRangeColumn(A,&stc,NULL)); 603 PetscCall(PetscMalloc1((dc+oc)/bs,&aux)); 604 for (i=0; i<dc/bs; i++) aux[i] = i+stc/bs; 605 for (i=0; i<oc/bs; i++) aux[i+dc/bs] = (ismpiaij ? garray[i*bs]/bs : garray[i]); 606 PetscCall(ISCreateBlock(comm,bs,(dc+oc)/bs,aux,PETSC_OWN_POINTER,&is)); 607 } else { 608 PetscCall(ISCreateBlock(comm,1,0,NULL,PETSC_OWN_POINTER,&is)); 609 } 610 PetscCall(ISLocalToGlobalMappingCreateIS(is,l2g)); 611 PetscCall(ISDestroy(&is)); 612 break; 613 default: 614 SETERRQ(comm,PETSC_ERR_ARG_WRONG,"Unsupported l2g disassembling type %d",mode); 615 } 616 PetscFunctionReturn(0); 617 } 618 619 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat) 620 { 621 Mat lA,Ad,Ao,B = NULL; 622 ISLocalToGlobalMapping rl2g,cl2g; 623 IS is; 624 MPI_Comm comm; 625 void *ptrs[2]; 626 const char *names[2] = {"_convert_csr_aux","_convert_csr_data"}; 627 const PetscInt *garray; 628 PetscScalar *dd,*od,*aa,*data; 629 const PetscInt *di,*dj,*oi,*oj; 630 const PetscInt *odi,*odj,*ooi,*ooj; 631 PetscInt *aux,*ii,*jj; 632 PetscInt bs,lc,dr,dc,oc,str,stc,nnz,i,jd,jo,cum; 633 PetscBool flg,ismpiaij,ismpibaij,was_inplace = PETSC_FALSE; 634 PetscMPIInt size; 635 636 PetscFunctionBegin; 637 PetscCall(PetscObjectGetComm((PetscObject)A,&comm)); 638 PetscCallMPI(MPI_Comm_size(comm,&size)); 639 if (size == 1) { 640 PetscCall(MatConvert_SeqXAIJ_IS(A,type,reuse,newmat)); 641 PetscFunctionReturn(0); 642 } 643 if (reuse != MAT_REUSE_MATRIX && A->cmap->N == A->rmap->N) { 644 PetscCall(MatMPIXAIJComputeLocalToGlobalMapping_Private(A,&rl2g)); 645 PetscCall(MatCreate(comm,&B)); 646 PetscCall(MatSetType(B,MATIS)); 647 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 648 PetscCall(MatSetLocalToGlobalMapping(B,rl2g,rl2g)); 649 PetscCall(MatGetBlockSize(A,&bs)); 650 PetscCall(MatSetBlockSize(B,bs)); 651 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 652 if (reuse == MAT_INPLACE_MATRIX) was_inplace = PETSC_TRUE; 653 reuse = MAT_REUSE_MATRIX; 654 } 655 if (reuse == MAT_REUSE_MATRIX) { 656 Mat *newlA, lA; 657 IS rows, cols; 658 const PetscInt *ridx, *cidx; 659 PetscInt rbs, cbs, nr, nc; 660 661 if (!B) B = *newmat; 662 PetscCall(MatISGetLocalToGlobalMapping(B,&rl2g,&cl2g)); 663 PetscCall(ISLocalToGlobalMappingGetBlockIndices(rl2g,&ridx)); 664 PetscCall(ISLocalToGlobalMappingGetBlockIndices(cl2g,&cidx)); 665 PetscCall(ISLocalToGlobalMappingGetSize(rl2g,&nr)); 666 PetscCall(ISLocalToGlobalMappingGetSize(cl2g,&nc)); 667 PetscCall(ISLocalToGlobalMappingGetBlockSize(rl2g,&rbs)); 668 PetscCall(ISLocalToGlobalMappingGetBlockSize(cl2g,&cbs)); 669 PetscCall(ISCreateBlock(comm,rbs,nr/rbs,ridx,PETSC_USE_POINTER,&rows)); 670 if (rl2g != cl2g) { 671 PetscCall(ISCreateBlock(comm,cbs,nc/cbs,cidx,PETSC_USE_POINTER,&cols)); 672 } else { 673 PetscCall(PetscObjectReference((PetscObject)rows)); 674 cols = rows; 675 } 676 PetscCall(MatISGetLocalMat(B,&lA)); 677 PetscCall(MatCreateSubMatrices(A,1,&rows,&cols,MAT_INITIAL_MATRIX,&newlA)); 678 PetscCall(MatConvert(newlA[0],MATSEQAIJ,MAT_INPLACE_MATRIX,&newlA[0])); 679 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(rl2g,&ridx)); 680 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(cl2g,&cidx)); 681 PetscCall(ISDestroy(&rows)); 682 PetscCall(ISDestroy(&cols)); 683 if (!lA->preallocated) { /* first time */ 684 PetscCall(MatDuplicate(newlA[0],MAT_COPY_VALUES,&lA)); 685 PetscCall(MatISSetLocalMat(B,lA)); 686 PetscCall(PetscObjectDereference((PetscObject)lA)); 687 } 688 PetscCall(MatCopy(newlA[0],lA,SAME_NONZERO_PATTERN)); 689 PetscCall(MatDestroySubMatrices(1,&newlA)); 690 PetscCall(MatISScaleDisassembling_Private(B)); 691 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 692 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 693 if (was_inplace) PetscCall(MatHeaderReplace(A,&B)); 694 else *newmat = B; 695 PetscFunctionReturn(0); 696 } 697 /* rectangular case, just compress out the column space */ 698 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ ,&ismpiaij)); 699 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A,MATMPIBAIJ,&ismpibaij)); 700 if (ismpiaij) { 701 bs = 1; 702 PetscCall(MatMPIAIJGetSeqAIJ(A,&Ad,&Ao,&garray)); 703 } else if (ismpibaij) { 704 PetscCall(MatGetBlockSize(A,&bs)); 705 PetscCall(MatMPIBAIJGetSeqBAIJ(A,&Ad,&Ao,&garray)); 706 PetscCall(MatConvert(Ad,MATSEQAIJ,MAT_INITIAL_MATRIX,&Ad)); 707 PetscCall(MatConvert(Ao,MATSEQAIJ,MAT_INITIAL_MATRIX,&Ao)); 708 } else SETERRQ(comm,PETSC_ERR_SUP,"Type %s",((PetscObject)A)->type_name); 709 PetscCall(MatSeqAIJGetArray(Ad,&dd)); 710 PetscCall(MatSeqAIJGetArray(Ao,&od)); 711 PetscCheck(garray,comm,PETSC_ERR_ARG_WRONGSTATE,"garray not present"); 712 713 /* access relevant information from MPIAIJ */ 714 PetscCall(MatGetOwnershipRange(A,&str,NULL)); 715 PetscCall(MatGetOwnershipRangeColumn(A,&stc,NULL)); 716 PetscCall(MatGetLocalSize(A,&dr,&dc)); 717 PetscCall(MatGetLocalSize(Ao,NULL,&oc)); 718 PetscCall(MatGetRowIJ(Ad,0,PETSC_FALSE,PETSC_FALSE,&i,&di,&dj,&flg)); 719 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure"); 720 PetscCall(MatGetRowIJ(Ao,0,PETSC_FALSE,PETSC_FALSE,&i,&oi,&oj,&flg)); 721 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get IJ structure"); 722 nnz = di[dr] + oi[dr]; 723 /* store original pointers to be restored later */ 724 odi = di; odj = dj; ooi = oi; ooj = oj; 725 726 /* generate l2g maps for rows and cols */ 727 PetscCall(ISCreateStride(comm,dr/bs,str/bs,1,&is)); 728 if (bs > 1) { 729 IS is2; 730 731 PetscCall(ISGetLocalSize(is,&i)); 732 PetscCall(ISGetIndices(is,(const PetscInt**)&aux)); 733 PetscCall(ISCreateBlock(comm,bs,i,aux,PETSC_COPY_VALUES,&is2)); 734 PetscCall(ISRestoreIndices(is,(const PetscInt**)&aux)); 735 PetscCall(ISDestroy(&is)); 736 is = is2; 737 } 738 PetscCall(ISLocalToGlobalMappingCreateIS(is,&rl2g)); 739 PetscCall(ISDestroy(&is)); 740 if (dr) { 741 PetscCall(PetscMalloc1((dc+oc)/bs,&aux)); 742 for (i=0; i<dc/bs; i++) aux[i] = i+stc/bs; 743 for (i=0; i<oc/bs; i++) aux[i+dc/bs] = garray[i]; 744 PetscCall(ISCreateBlock(comm,bs,(dc+oc)/bs,aux,PETSC_OWN_POINTER,&is)); 745 lc = dc+oc; 746 } else { 747 PetscCall(ISCreateBlock(comm,bs,0,NULL,PETSC_OWN_POINTER,&is)); 748 lc = 0; 749 } 750 PetscCall(ISLocalToGlobalMappingCreateIS(is,&cl2g)); 751 PetscCall(ISDestroy(&is)); 752 753 /* create MATIS object */ 754 PetscCall(MatCreate(comm,&B)); 755 PetscCall(MatSetSizes(B,dr,dc,PETSC_DECIDE,PETSC_DECIDE)); 756 PetscCall(MatSetType(B,MATIS)); 757 PetscCall(MatSetBlockSize(B,bs)); 758 PetscCall(MatSetLocalToGlobalMapping(B,rl2g,cl2g)); 759 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 760 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 761 762 /* merge local matrices */ 763 PetscCall(PetscMalloc1(nnz+dr+1,&aux)); 764 PetscCall(PetscMalloc1(nnz,&data)); 765 ii = aux; 766 jj = aux+dr+1; 767 aa = data; 768 *ii = *(di++) + *(oi++); 769 for (jd=0,jo=0,cum=0;*ii<nnz;cum++) 770 { 771 for (;jd<*di;jd++) { *jj++ = *dj++; *aa++ = *dd++; } 772 for (;jo<*oi;jo++) { *jj++ = *oj++ + dc; *aa++ = *od++; } 773 *(++ii) = *(di++) + *(oi++); 774 } 775 for (;cum<dr;cum++) *(++ii) = nnz; 776 777 PetscCall(MatRestoreRowIJ(Ad,0,PETSC_FALSE,PETSC_FALSE,&i,&odi,&odj,&flg)); 778 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure"); 779 PetscCall(MatRestoreRowIJ(Ao,0,PETSC_FALSE,PETSC_FALSE,&i,&ooi,&ooj,&flg)); 780 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot restore IJ structure"); 781 PetscCall(MatSeqAIJRestoreArray(Ad,&dd)); 782 PetscCall(MatSeqAIJRestoreArray(Ao,&od)); 783 784 ii = aux; 785 jj = aux+dr+1; 786 aa = data; 787 PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,lc,ii,jj,aa,&lA)); 788 789 /* create containers to destroy the data */ 790 ptrs[0] = aux; 791 ptrs[1] = data; 792 for (i=0; i<2; i++) { 793 PetscContainer c; 794 795 PetscCall(PetscContainerCreate(PETSC_COMM_SELF,&c)); 796 PetscCall(PetscContainerSetPointer(c,ptrs[i])); 797 PetscCall(PetscContainerSetUserDestroy(c,PetscContainerUserDestroyDefault)); 798 PetscCall(PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c)); 799 PetscCall(PetscContainerDestroy(&c)); 800 } 801 if (ismpibaij) { /* destroy converted local matrices */ 802 PetscCall(MatDestroy(&Ad)); 803 PetscCall(MatDestroy(&Ao)); 804 } 805 806 /* finalize matrix */ 807 PetscCall(MatISSetLocalMat(B,lA)); 808 PetscCall(MatDestroy(&lA)); 809 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 810 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 811 if (reuse == MAT_INPLACE_MATRIX) { 812 PetscCall(MatHeaderReplace(A,&B)); 813 } else *newmat = B; 814 PetscFunctionReturn(0); 815 } 816 817 PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat) 818 { 819 Mat **nest,*snest,**rnest,lA,B; 820 IS *iscol,*isrow,*islrow,*islcol; 821 ISLocalToGlobalMapping rl2g,cl2g; 822 MPI_Comm comm; 823 PetscInt *lr,*lc,*l2gidxs; 824 PetscInt i,j,nr,nc,rbs,cbs; 825 PetscBool convert,lreuse,*istrans; 826 827 PetscFunctionBegin; 828 PetscCall(MatNestGetSubMats(A,&nr,&nc,&nest)); 829 lreuse = PETSC_FALSE; 830 rnest = NULL; 831 if (reuse == MAT_REUSE_MATRIX) { 832 PetscBool ismatis,isnest; 833 834 PetscCall(PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis)); 835 PetscCheck(ismatis,PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name); 836 PetscCall(MatISGetLocalMat(*newmat,&lA)); 837 PetscCall(PetscObjectTypeCompare((PetscObject)lA,MATNEST,&isnest)); 838 if (isnest) { 839 PetscCall(MatNestGetSubMats(lA,&i,&j,&rnest)); 840 lreuse = (PetscBool)(i == nr && j == nc); 841 if (!lreuse) rnest = NULL; 842 } 843 } 844 PetscCall(PetscObjectGetComm((PetscObject)A,&comm)); 845 PetscCall(PetscCalloc2(nr,&lr,nc,&lc)); 846 PetscCall(PetscCalloc6(nr,&isrow,nc,&iscol,nr,&islrow,nc,&islcol,nr*nc,&snest,nr*nc,&istrans)); 847 PetscCall(MatNestGetISs(A,isrow,iscol)); 848 for (i=0;i<nr;i++) { 849 for (j=0;j<nc;j++) { 850 PetscBool ismatis; 851 PetscInt l1,l2,lb1,lb2,ij=i*nc+j; 852 853 /* Null matrix pointers are allowed in MATNEST */ 854 if (!nest[i][j]) continue; 855 856 /* Nested matrices should be of type MATIS */ 857 PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j],MATTRANSPOSEMAT,&istrans[ij])); 858 if (istrans[ij]) { 859 Mat T,lT; 860 PetscCall(MatTransposeGetMat(nest[i][j],&T)); 861 PetscCall(PetscObjectTypeCompare((PetscObject)T,MATIS,&ismatis)); 862 PetscCheck(ismatis,comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") (transposed) is not of type MATIS",i,j); 863 PetscCall(MatISGetLocalMat(T,&lT)); 864 PetscCall(MatCreateTranspose(lT,&snest[ij])); 865 } else { 866 PetscCall(PetscObjectTypeCompare((PetscObject)nest[i][j],MATIS,&ismatis)); 867 PetscCheck(ismatis,comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") is not of type MATIS",i,j); 868 PetscCall(MatISGetLocalMat(nest[i][j],&snest[ij])); 869 } 870 871 /* Check compatibility of local sizes */ 872 PetscCall(MatGetSize(snest[ij],&l1,&l2)); 873 PetscCall(MatGetBlockSizes(snest[ij],&lb1,&lb2)); 874 if (!l1 || !l2) continue; 875 PetscCheck(!lr[i] || l1 == lr[i],PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT,i,j,lr[i],l1); 876 PetscCheck(!lc[j] || l2 == lc[j],PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid local size %" PetscInt_FMT " != %" PetscInt_FMT,i,j,lc[j],l2); 877 lr[i] = l1; 878 lc[j] = l2; 879 880 /* check compatibilty for local matrix reusage */ 881 if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE; 882 } 883 } 884 885 if (PetscDefined (USE_DEBUG)) { 886 /* Check compatibility of l2g maps for rows */ 887 for (i=0;i<nr;i++) { 888 rl2g = NULL; 889 for (j=0;j<nc;j++) { 890 PetscInt n1,n2; 891 892 if (!nest[i][j]) continue; 893 if (istrans[i*nc+j]) { 894 Mat T; 895 896 PetscCall(MatTransposeGetMat(nest[i][j],&T)); 897 PetscCall(MatISGetLocalToGlobalMapping(T,NULL,&cl2g)); 898 } else { 899 PetscCall(MatISGetLocalToGlobalMapping(nest[i][j],&cl2g,NULL)); 900 } 901 PetscCall(ISLocalToGlobalMappingGetSize(cl2g,&n1)); 902 if (!n1) continue; 903 if (!rl2g) { 904 rl2g = cl2g; 905 } else { 906 const PetscInt *idxs1,*idxs2; 907 PetscBool same; 908 909 PetscCall(ISLocalToGlobalMappingGetSize(rl2g,&n2)); 910 PetscCheck(n1 == n2,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid row l2gmap size %" PetscInt_FMT " != %" PetscInt_FMT,i,j,n1,n2); 911 PetscCall(ISLocalToGlobalMappingGetIndices(cl2g,&idxs1)); 912 PetscCall(ISLocalToGlobalMappingGetIndices(rl2g,&idxs2)); 913 PetscCall(PetscArraycmp(idxs1,idxs2,n1,&same)); 914 PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1)); 915 PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2)); 916 PetscCheck(same,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid row l2gmap",i,j); 917 } 918 } 919 } 920 /* Check compatibility of l2g maps for columns */ 921 for (i=0;i<nc;i++) { 922 rl2g = NULL; 923 for (j=0;j<nr;j++) { 924 PetscInt n1,n2; 925 926 if (!nest[j][i]) continue; 927 if (istrans[j*nc+i]) { 928 Mat T; 929 930 PetscCall(MatTransposeGetMat(nest[j][i],&T)); 931 PetscCall(MatISGetLocalToGlobalMapping(T,&cl2g,NULL)); 932 } else { 933 PetscCall(MatISGetLocalToGlobalMapping(nest[j][i],NULL,&cl2g)); 934 } 935 PetscCall(ISLocalToGlobalMappingGetSize(cl2g,&n1)); 936 if (!n1) continue; 937 if (!rl2g) { 938 rl2g = cl2g; 939 } else { 940 const PetscInt *idxs1,*idxs2; 941 PetscBool same; 942 943 PetscCall(ISLocalToGlobalMappingGetSize(rl2g,&n2)); 944 PetscCheck(n1 == n2,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap size %" PetscInt_FMT " != %" PetscInt_FMT,j,i,n1,n2); 945 PetscCall(ISLocalToGlobalMappingGetIndices(cl2g,&idxs1)); 946 PetscCall(ISLocalToGlobalMappingGetIndices(rl2g,&idxs2)); 947 PetscCall(PetscArraycmp(idxs1,idxs2,n1,&same)); 948 PetscCall(ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1)); 949 PetscCall(ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2)); 950 PetscCheck(same,PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%" PetscInt_FMT ",%" PetscInt_FMT ") has invalid column l2gmap",j,i); 951 } 952 } 953 } 954 } 955 956 B = NULL; 957 if (reuse != MAT_REUSE_MATRIX) { 958 PetscInt stl; 959 960 /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */ 961 for (i=0,stl=0;i<nr;i++) stl += lr[i]; 962 PetscCall(PetscMalloc1(stl,&l2gidxs)); 963 for (i=0,stl=0;i<nr;i++) { 964 Mat usedmat; 965 Mat_IS *matis; 966 const PetscInt *idxs; 967 968 /* local IS for local NEST */ 969 PetscCall(ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i])); 970 971 /* l2gmap */ 972 j = 0; 973 usedmat = nest[i][j]; 974 while (!usedmat && j < nc-1) usedmat = nest[i][++j]; 975 PetscCheck(usedmat,comm,PETSC_ERR_SUP,"Cannot find valid row mat"); 976 977 if (istrans[i*nc+j]) { 978 Mat T; 979 PetscCall(MatTransposeGetMat(usedmat,&T)); 980 usedmat = T; 981 } 982 matis = (Mat_IS*)(usedmat->data); 983 PetscCall(ISGetIndices(isrow[i],&idxs)); 984 if (istrans[i*nc+j]) { 985 PetscCall(PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 986 PetscCall(PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 987 } else { 988 PetscCall(PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 989 PetscCall(PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 990 } 991 PetscCall(ISRestoreIndices(isrow[i],&idxs)); 992 stl += lr[i]; 993 } 994 PetscCall(ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&rl2g)); 995 996 /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */ 997 for (i=0,stl=0;i<nc;i++) stl += lc[i]; 998 PetscCall(PetscMalloc1(stl,&l2gidxs)); 999 for (i=0,stl=0;i<nc;i++) { 1000 Mat usedmat; 1001 Mat_IS *matis; 1002 const PetscInt *idxs; 1003 1004 /* local IS for local NEST */ 1005 PetscCall(ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i])); 1006 1007 /* l2gmap */ 1008 j = 0; 1009 usedmat = nest[j][i]; 1010 while (!usedmat && j < nr-1) usedmat = nest[++j][i]; 1011 PetscCheck(usedmat,comm,PETSC_ERR_SUP,"Cannot find valid column mat"); 1012 if (istrans[j*nc+i]) { 1013 Mat T; 1014 PetscCall(MatTransposeGetMat(usedmat,&T)); 1015 usedmat = T; 1016 } 1017 matis = (Mat_IS*)(usedmat->data); 1018 PetscCall(ISGetIndices(iscol[i],&idxs)); 1019 if (istrans[j*nc+i]) { 1020 PetscCall(PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 1021 PetscCall(PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 1022 } else { 1023 PetscCall(PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 1024 PetscCall(PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 1025 } 1026 PetscCall(ISRestoreIndices(iscol[i],&idxs)); 1027 stl += lc[i]; 1028 } 1029 PetscCall(ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&cl2g)); 1030 1031 /* Create MATIS */ 1032 PetscCall(MatCreate(comm,&B)); 1033 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 1034 PetscCall(MatGetBlockSizes(A,&rbs,&cbs)); 1035 PetscCall(MatSetBlockSizes(B,rbs,cbs)); 1036 PetscCall(MatSetType(B,MATIS)); 1037 PetscCall(MatISSetLocalMatType(B,MATNEST)); 1038 { /* hack : avoid setup of scatters */ 1039 Mat_IS *matis = (Mat_IS*)(B->data); 1040 matis->islocalref = PETSC_TRUE; 1041 } 1042 PetscCall(MatSetLocalToGlobalMapping(B,rl2g,cl2g)); 1043 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 1044 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 1045 PetscCall(MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA)); 1046 PetscCall(MatNestSetVecType(lA,VECNEST)); 1047 for (i=0;i<nr*nc;i++) { 1048 if (istrans[i]) { 1049 PetscCall(MatDestroy(&snest[i])); 1050 } 1051 } 1052 PetscCall(MatISSetLocalMat(B,lA)); 1053 PetscCall(MatDestroy(&lA)); 1054 { /* hack : setup of scatters done here */ 1055 Mat_IS *matis = (Mat_IS*)(B->data); 1056 1057 matis->islocalref = PETSC_FALSE; 1058 PetscCall(MatISSetUpScatters_Private(B)); 1059 } 1060 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1061 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1062 if (reuse == MAT_INPLACE_MATRIX) { 1063 PetscCall(MatHeaderReplace(A,&B)); 1064 } else { 1065 *newmat = B; 1066 } 1067 } else { 1068 if (lreuse) { 1069 PetscCall(MatISGetLocalMat(*newmat,&lA)); 1070 for (i=0;i<nr;i++) { 1071 for (j=0;j<nc;j++) { 1072 if (snest[i*nc+j]) { 1073 PetscCall(MatNestSetSubMat(lA,i,j,snest[i*nc+j])); 1074 if (istrans[i*nc+j]) { 1075 PetscCall(MatDestroy(&snest[i*nc+j])); 1076 } 1077 } 1078 } 1079 } 1080 } else { 1081 PetscInt stl; 1082 for (i=0,stl=0;i<nr;i++) { 1083 PetscCall(ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i])); 1084 stl += lr[i]; 1085 } 1086 for (i=0,stl=0;i<nc;i++) { 1087 PetscCall(ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i])); 1088 stl += lc[i]; 1089 } 1090 PetscCall(MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA)); 1091 for (i=0;i<nr*nc;i++) { 1092 if (istrans[i]) { 1093 PetscCall(MatDestroy(&snest[i])); 1094 } 1095 } 1096 PetscCall(MatISSetLocalMat(*newmat,lA)); 1097 PetscCall(MatDestroy(&lA)); 1098 } 1099 PetscCall(MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY)); 1100 PetscCall(MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY)); 1101 } 1102 1103 /* Create local matrix in MATNEST format */ 1104 convert = PETSC_FALSE; 1105 PetscCall(PetscOptionsGetBool(NULL,((PetscObject)A)->prefix,"-matis_convert_local_nest",&convert,NULL)); 1106 if (convert) { 1107 Mat M; 1108 MatISLocalFields lf; 1109 PetscContainer c; 1110 1111 PetscCall(MatISGetLocalMat(*newmat,&lA)); 1112 PetscCall(MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M)); 1113 PetscCall(MatISSetLocalMat(*newmat,M)); 1114 PetscCall(MatDestroy(&M)); 1115 1116 /* attach local fields to the matrix */ 1117 PetscCall(PetscNew(&lf)); 1118 PetscCall(PetscMalloc2(nr,&lf->rf,nc,&lf->cf)); 1119 for (i=0;i<nr;i++) { 1120 PetscInt n,st; 1121 1122 PetscCall(ISGetLocalSize(islrow[i],&n)); 1123 PetscCall(ISStrideGetInfo(islrow[i],&st,NULL)); 1124 PetscCall(ISCreateStride(comm,n,st,1,&lf->rf[i])); 1125 } 1126 for (i=0;i<nc;i++) { 1127 PetscInt n,st; 1128 1129 PetscCall(ISGetLocalSize(islcol[i],&n)); 1130 PetscCall(ISStrideGetInfo(islcol[i],&st,NULL)); 1131 PetscCall(ISCreateStride(comm,n,st,1,&lf->cf[i])); 1132 } 1133 lf->nr = nr; 1134 lf->nc = nc; 1135 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)),&c)); 1136 PetscCall(PetscContainerSetPointer(c,lf)); 1137 PetscCall(PetscContainerSetUserDestroy(c,MatISContainerDestroyFields_Private)); 1138 PetscCall(PetscObjectCompose((PetscObject)(*newmat),"_convert_nest_lfields",(PetscObject)c)); 1139 PetscCall(PetscContainerDestroy(&c)); 1140 } 1141 1142 /* Free workspace */ 1143 for (i=0;i<nr;i++) { 1144 PetscCall(ISDestroy(&islrow[i])); 1145 } 1146 for (i=0;i<nc;i++) { 1147 PetscCall(ISDestroy(&islcol[i])); 1148 } 1149 PetscCall(PetscFree6(isrow,iscol,islrow,islcol,snest,istrans)); 1150 PetscCall(PetscFree2(lr,lc)); 1151 PetscFunctionReturn(0); 1152 } 1153 1154 static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r) 1155 { 1156 Mat_IS *matis = (Mat_IS*)A->data; 1157 Vec ll,rr; 1158 const PetscScalar *Y,*X; 1159 PetscScalar *x,*y; 1160 1161 PetscFunctionBegin; 1162 if (l) { 1163 ll = matis->y; 1164 PetscCall(VecGetArrayRead(l,&Y)); 1165 PetscCall(VecGetArray(ll,&y)); 1166 PetscCall(PetscSFBcastBegin(matis->sf,MPIU_SCALAR,Y,y,MPI_REPLACE)); 1167 } else { 1168 ll = NULL; 1169 } 1170 if (r) { 1171 rr = matis->x; 1172 PetscCall(VecGetArrayRead(r,&X)); 1173 PetscCall(VecGetArray(rr,&x)); 1174 PetscCall(PetscSFBcastBegin(matis->csf,MPIU_SCALAR,X,x,MPI_REPLACE)); 1175 } else { 1176 rr = NULL; 1177 } 1178 if (ll) { 1179 PetscCall(PetscSFBcastEnd(matis->sf,MPIU_SCALAR,Y,y,MPI_REPLACE)); 1180 PetscCall(VecRestoreArrayRead(l,&Y)); 1181 PetscCall(VecRestoreArray(ll,&y)); 1182 } 1183 if (rr) { 1184 PetscCall(PetscSFBcastEnd(matis->csf,MPIU_SCALAR,X,x,MPI_REPLACE)); 1185 PetscCall(VecRestoreArrayRead(r,&X)); 1186 PetscCall(VecRestoreArray(rr,&x)); 1187 } 1188 PetscCall(MatDiagonalScale(matis->A,ll,rr)); 1189 PetscFunctionReturn(0); 1190 } 1191 1192 static PetscErrorCode MatGetInfo_IS(Mat A,MatInfoType flag,MatInfo *ginfo) 1193 { 1194 Mat_IS *matis = (Mat_IS*)A->data; 1195 MatInfo info; 1196 PetscLogDouble isend[6],irecv[6]; 1197 PetscInt bs; 1198 1199 PetscFunctionBegin; 1200 PetscCall(MatGetBlockSize(A,&bs)); 1201 if (matis->A->ops->getinfo) { 1202 PetscCall(MatGetInfo(matis->A,MAT_LOCAL,&info)); 1203 isend[0] = info.nz_used; 1204 isend[1] = info.nz_allocated; 1205 isend[2] = info.nz_unneeded; 1206 isend[3] = info.memory; 1207 isend[4] = info.mallocs; 1208 } else { 1209 isend[0] = 0.; 1210 isend[1] = 0.; 1211 isend[2] = 0.; 1212 isend[3] = 0.; 1213 isend[4] = 0.; 1214 } 1215 isend[5] = matis->A->num_ass; 1216 if (flag == MAT_LOCAL) { 1217 ginfo->nz_used = isend[0]; 1218 ginfo->nz_allocated = isend[1]; 1219 ginfo->nz_unneeded = isend[2]; 1220 ginfo->memory = isend[3]; 1221 ginfo->mallocs = isend[4]; 1222 ginfo->assemblies = isend[5]; 1223 } else if (flag == MAT_GLOBAL_MAX) { 1224 PetscCall(MPIU_Allreduce(isend,irecv,6,MPIU_PETSCLOGDOUBLE,MPI_MAX,PetscObjectComm((PetscObject)A))); 1225 1226 ginfo->nz_used = irecv[0]; 1227 ginfo->nz_allocated = irecv[1]; 1228 ginfo->nz_unneeded = irecv[2]; 1229 ginfo->memory = irecv[3]; 1230 ginfo->mallocs = irecv[4]; 1231 ginfo->assemblies = irecv[5]; 1232 } else if (flag == MAT_GLOBAL_SUM) { 1233 PetscCall(MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_SUM,PetscObjectComm((PetscObject)A))); 1234 1235 ginfo->nz_used = irecv[0]; 1236 ginfo->nz_allocated = irecv[1]; 1237 ginfo->nz_unneeded = irecv[2]; 1238 ginfo->memory = irecv[3]; 1239 ginfo->mallocs = irecv[4]; 1240 ginfo->assemblies = A->num_ass; 1241 } 1242 ginfo->block_size = bs; 1243 ginfo->fill_ratio_given = 0; 1244 ginfo->fill_ratio_needed = 0; 1245 ginfo->factor_mallocs = 0; 1246 PetscFunctionReturn(0); 1247 } 1248 1249 static PetscErrorCode MatTranspose_IS(Mat A,MatReuse reuse,Mat *B) 1250 { 1251 Mat C,lC,lA; 1252 1253 PetscFunctionBegin; 1254 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) { 1255 ISLocalToGlobalMapping rl2g,cl2g; 1256 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 1257 PetscCall(MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N)); 1258 PetscCall(MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs))); 1259 PetscCall(MatSetType(C,MATIS)); 1260 PetscCall(MatGetLocalToGlobalMapping(A,&rl2g,&cl2g)); 1261 PetscCall(MatSetLocalToGlobalMapping(C,cl2g,rl2g)); 1262 } else C = *B; 1263 1264 /* perform local transposition */ 1265 PetscCall(MatISGetLocalMat(A,&lA)); 1266 PetscCall(MatTranspose(lA,MAT_INITIAL_MATRIX,&lC)); 1267 PetscCall(MatSetLocalToGlobalMapping(lC,lA->cmap->mapping,lA->rmap->mapping)); 1268 PetscCall(MatISSetLocalMat(C,lC)); 1269 PetscCall(MatDestroy(&lC)); 1270 1271 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) { 1272 *B = C; 1273 } else { 1274 PetscCall(MatHeaderMerge(A,&C)); 1275 } 1276 PetscCall(MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY)); 1277 PetscCall(MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY)); 1278 PetscFunctionReturn(0); 1279 } 1280 1281 static PetscErrorCode MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode) 1282 { 1283 Mat_IS *is = (Mat_IS*)A->data; 1284 1285 PetscFunctionBegin; 1286 if (D) { /* MatShift_IS pass D = NULL */ 1287 PetscCall(VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD)); 1288 PetscCall(VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD)); 1289 } 1290 PetscCall(VecPointwiseDivide(is->y,is->y,is->counter)); 1291 PetscCall(MatDiagonalSet(is->A,is->y,insmode)); 1292 PetscFunctionReturn(0); 1293 } 1294 1295 static PetscErrorCode MatShift_IS(Mat A,PetscScalar a) 1296 { 1297 Mat_IS *is = (Mat_IS*)A->data; 1298 1299 PetscFunctionBegin; 1300 PetscCall(VecSet(is->y,a)); 1301 PetscCall(MatDiagonalSet_IS(A,NULL,ADD_VALUES)); 1302 PetscFunctionReturn(0); 1303 } 1304 1305 static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 1306 { 1307 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1308 1309 PetscFunctionBegin; 1310 PetscCheck(m <= MATIS_MAX_ENTRIES_INSERTION && n <= MATIS_MAX_ENTRIES_INSERTION,PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column indices must be <= %d: they are %" PetscInt_FMT " %" PetscInt_FMT,MATIS_MAX_ENTRIES_INSERTION,m,n); 1311 PetscCall(ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l)); 1312 PetscCall(ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l)); 1313 PetscCall(MatSetValuesLocal_IS(A,m,rows_l,n,cols_l,values,addv)); 1314 PetscFunctionReturn(0); 1315 } 1316 1317 static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 1318 { 1319 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1320 1321 PetscFunctionBegin; 1322 PetscCheck(m <= MATIS_MAX_ENTRIES_INSERTION && n <= MATIS_MAX_ENTRIES_INSERTION,PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column block indices must be <= %d: they are %" PetscInt_FMT " %" PetscInt_FMT,MATIS_MAX_ENTRIES_INSERTION,m,n); 1323 PetscCall(ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l)); 1324 PetscCall(ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l)); 1325 PetscCall(MatSetValuesBlockedLocal_IS(A,m,rows_l,n,cols_l,values,addv)); 1326 PetscFunctionReturn(0); 1327 } 1328 1329 static PetscErrorCode MatCreateSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat) 1330 { 1331 Mat locmat,newlocmat; 1332 Mat_IS *newmatis; 1333 const PetscInt *idxs; 1334 PetscInt i,m,n; 1335 1336 PetscFunctionBegin; 1337 if (scall == MAT_REUSE_MATRIX) { 1338 PetscBool ismatis; 1339 1340 PetscCall(PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis)); 1341 PetscCheck(ismatis,PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Not of MATIS type"); 1342 newmatis = (Mat_IS*)(*newmat)->data; 1343 PetscCheck(newmatis->getsub_ris,PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local row IS"); 1344 PetscCheck(newmatis->getsub_cis,PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local col IS"); 1345 } 1346 /* irow and icol may not have duplicate entries */ 1347 if (PetscDefined(USE_DEBUG)) { 1348 Vec rtest,ltest; 1349 const PetscScalar *array; 1350 1351 PetscCall(MatCreateVecs(mat,<est,&rtest)); 1352 PetscCall(ISGetLocalSize(irow,&n)); 1353 PetscCall(ISGetIndices(irow,&idxs)); 1354 for (i=0;i<n;i++) { 1355 PetscCall(VecSetValue(rtest,idxs[i],1.0,ADD_VALUES)); 1356 } 1357 PetscCall(VecAssemblyBegin(rtest)); 1358 PetscCall(VecAssemblyEnd(rtest)); 1359 PetscCall(VecGetLocalSize(rtest,&n)); 1360 PetscCall(VecGetOwnershipRange(rtest,&m,NULL)); 1361 PetscCall(VecGetArrayRead(rtest,&array)); 1362 for (i=0;i<n;i++) PetscCheck(array[i] == 0. || array[i] == 1.,PETSC_COMM_SELF,PETSC_ERR_SUP,"Index %" PetscInt_FMT " counted %" PetscInt_FMT " times! Irow may not have duplicate entries",i+m,(PetscInt)PetscRealPart(array[i])); 1363 PetscCall(VecRestoreArrayRead(rtest,&array)); 1364 PetscCall(ISRestoreIndices(irow,&idxs)); 1365 PetscCall(ISGetLocalSize(icol,&n)); 1366 PetscCall(ISGetIndices(icol,&idxs)); 1367 for (i=0;i<n;i++) { 1368 PetscCall(VecSetValue(ltest,idxs[i],1.0,ADD_VALUES)); 1369 } 1370 PetscCall(VecAssemblyBegin(ltest)); 1371 PetscCall(VecAssemblyEnd(ltest)); 1372 PetscCall(VecGetLocalSize(ltest,&n)); 1373 PetscCall(VecGetOwnershipRange(ltest,&m,NULL)); 1374 PetscCall(VecGetArrayRead(ltest,&array)); 1375 for (i=0;i<n;i++) PetscCheck(array[i] == 0. || array[i] == 1.,PETSC_COMM_SELF,PETSC_ERR_SUP,"Index %" PetscInt_FMT " counted %" PetscInt_FMT " times! Icol may not have duplicate entries",i+m,(PetscInt)PetscRealPart(array[i])); 1376 PetscCall(VecRestoreArrayRead(ltest,&array)); 1377 PetscCall(ISRestoreIndices(icol,&idxs)); 1378 PetscCall(VecDestroy(&rtest)); 1379 PetscCall(VecDestroy(<est)); 1380 } 1381 if (scall == MAT_INITIAL_MATRIX) { 1382 Mat_IS *matis = (Mat_IS*)mat->data; 1383 ISLocalToGlobalMapping rl2g; 1384 IS is; 1385 PetscInt *lidxs,*lgidxs,*newgidxs; 1386 PetscInt ll,newloc,irbs,icbs,arbs,acbs,rbs,cbs; 1387 PetscBool cong; 1388 MPI_Comm comm; 1389 1390 PetscCall(PetscObjectGetComm((PetscObject)mat,&comm)); 1391 PetscCall(MatGetBlockSizes(mat,&arbs,&acbs)); 1392 PetscCall(ISGetBlockSize(irow,&irbs)); 1393 PetscCall(ISGetBlockSize(icol,&icbs)); 1394 rbs = arbs == irbs ? irbs : 1; 1395 cbs = acbs == icbs ? icbs : 1; 1396 PetscCall(ISGetLocalSize(irow,&m)); 1397 PetscCall(ISGetLocalSize(icol,&n)); 1398 PetscCall(MatCreate(comm,newmat)); 1399 PetscCall(MatSetType(*newmat,MATIS)); 1400 PetscCall(MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE)); 1401 PetscCall(MatSetBlockSizes(*newmat,rbs,cbs)); 1402 /* communicate irow to their owners in the layout */ 1403 PetscCall(ISGetIndices(irow,&idxs)); 1404 PetscCall(PetscLayoutMapLocal(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs)); 1405 PetscCall(ISRestoreIndices(irow,&idxs)); 1406 PetscCall(PetscArrayzero(matis->sf_rootdata,matis->sf->nroots)); 1407 for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1; 1408 PetscCall(PetscFree(lidxs)); 1409 PetscCall(PetscFree(lgidxs)); 1410 PetscCall(PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE)); 1411 PetscCall(PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE)); 1412 for (i=0,newloc=0;i<matis->sf->nleaves;i++) if (matis->sf_leafdata[i]) newloc++; 1413 PetscCall(PetscMalloc1(newloc,&newgidxs)); 1414 PetscCall(PetscMalloc1(newloc,&lidxs)); 1415 for (i=0,newloc=0;i<matis->sf->nleaves;i++) 1416 if (matis->sf_leafdata[i]) { 1417 lidxs[newloc] = i; 1418 newgidxs[newloc++] = matis->sf_leafdata[i]-1; 1419 } 1420 PetscCall(ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is)); 1421 PetscCall(ISLocalToGlobalMappingCreateIS(is,&rl2g)); 1422 PetscCall(ISLocalToGlobalMappingSetBlockSize(rl2g,rbs)); 1423 PetscCall(ISDestroy(&is)); 1424 /* local is to extract local submatrix */ 1425 newmatis = (Mat_IS*)(*newmat)->data; 1426 PetscCall(ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris)); 1427 PetscCall(MatHasCongruentLayouts(mat,&cong)); 1428 if (cong && irow == icol && matis->csf == matis->sf) { 1429 PetscCall(MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g)); 1430 PetscCall(PetscObjectReference((PetscObject)newmatis->getsub_ris)); 1431 newmatis->getsub_cis = newmatis->getsub_ris; 1432 } else { 1433 ISLocalToGlobalMapping cl2g; 1434 1435 /* communicate icol to their owners in the layout */ 1436 PetscCall(ISGetIndices(icol,&idxs)); 1437 PetscCall(PetscLayoutMapLocal(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs)); 1438 PetscCall(ISRestoreIndices(icol,&idxs)); 1439 PetscCall(PetscArrayzero(matis->csf_rootdata,matis->csf->nroots)); 1440 for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1; 1441 PetscCall(PetscFree(lidxs)); 1442 PetscCall(PetscFree(lgidxs)); 1443 PetscCall(PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata,MPI_REPLACE)); 1444 PetscCall(PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata,MPI_REPLACE)); 1445 for (i=0,newloc=0;i<matis->csf->nleaves;i++) if (matis->csf_leafdata[i]) newloc++; 1446 PetscCall(PetscMalloc1(newloc,&newgidxs)); 1447 PetscCall(PetscMalloc1(newloc,&lidxs)); 1448 for (i=0,newloc=0;i<matis->csf->nleaves;i++) 1449 if (matis->csf_leafdata[i]) { 1450 lidxs[newloc] = i; 1451 newgidxs[newloc++] = matis->csf_leafdata[i]-1; 1452 } 1453 PetscCall(ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is)); 1454 PetscCall(ISLocalToGlobalMappingCreateIS(is,&cl2g)); 1455 PetscCall(ISLocalToGlobalMappingSetBlockSize(cl2g,cbs)); 1456 PetscCall(ISDestroy(&is)); 1457 /* local is to extract local submatrix */ 1458 PetscCall(ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis)); 1459 PetscCall(MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g)); 1460 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 1461 } 1462 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 1463 } else { 1464 PetscCall(MatISGetLocalMat(*newmat,&newlocmat)); 1465 } 1466 PetscCall(MatISGetLocalMat(mat,&locmat)); 1467 newmatis = (Mat_IS*)(*newmat)->data; 1468 PetscCall(MatCreateSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat)); 1469 if (scall == MAT_INITIAL_MATRIX) { 1470 PetscCall(MatISSetLocalMat(*newmat,newlocmat)); 1471 PetscCall(MatDestroy(&newlocmat)); 1472 } 1473 PetscCall(MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY)); 1474 PetscCall(MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY)); 1475 PetscFunctionReturn(0); 1476 } 1477 1478 static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str) 1479 { 1480 Mat_IS *a = (Mat_IS*)A->data,*b; 1481 PetscBool ismatis; 1482 1483 PetscFunctionBegin; 1484 PetscCall(PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis)); 1485 PetscCheck(ismatis,PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Need to be implemented"); 1486 b = (Mat_IS*)B->data; 1487 PetscCall(MatCopy(a->A,b->A,str)); 1488 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 1489 PetscFunctionReturn(0); 1490 } 1491 1492 static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool *missing,PetscInt *d) 1493 { 1494 Vec v; 1495 const PetscScalar *array; 1496 PetscInt i,n; 1497 1498 PetscFunctionBegin; 1499 *missing = PETSC_FALSE; 1500 PetscCall(MatCreateVecs(A,NULL,&v)); 1501 PetscCall(MatGetDiagonal(A,v)); 1502 PetscCall(VecGetLocalSize(v,&n)); 1503 PetscCall(VecGetArrayRead(v,&array)); 1504 for (i=0;i<n;i++) if (array[i] == 0.) break; 1505 PetscCall(VecRestoreArrayRead(v,&array)); 1506 PetscCall(VecDestroy(&v)); 1507 if (i != n) *missing = PETSC_TRUE; 1508 if (d) { 1509 *d = -1; 1510 if (*missing) { 1511 PetscInt rstart; 1512 PetscCall(MatGetOwnershipRange(A,&rstart,NULL)); 1513 *d = i+rstart; 1514 } 1515 } 1516 PetscFunctionReturn(0); 1517 } 1518 1519 static PetscErrorCode MatISSetUpSF_IS(Mat B) 1520 { 1521 Mat_IS *matis = (Mat_IS*)(B->data); 1522 const PetscInt *gidxs; 1523 PetscInt nleaves; 1524 1525 PetscFunctionBegin; 1526 if (matis->sf) PetscFunctionReturn(0); 1527 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf)); 1528 PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping,&gidxs)); 1529 PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping,&nleaves)); 1530 PetscCall(PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs)); 1531 PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping,&gidxs)); 1532 PetscCall(PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata)); 1533 if (matis->rmapping != matis->cmapping) { /* setup SF for columns */ 1534 PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping,&nleaves)); 1535 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf)); 1536 PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping,&gidxs)); 1537 PetscCall(PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs)); 1538 PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping,&gidxs)); 1539 PetscCall(PetscMalloc2(matis->csf->nroots,&matis->csf_rootdata,matis->csf->nleaves,&matis->csf_leafdata)); 1540 } else { 1541 matis->csf = matis->sf; 1542 matis->csf_leafdata = matis->sf_leafdata; 1543 matis->csf_rootdata = matis->sf_rootdata; 1544 } 1545 PetscFunctionReturn(0); 1546 } 1547 1548 /*@ 1549 MatISStoreL2L - Store local-to-local operators during the Galerkin process of MatPtAP. 1550 1551 Collective 1552 1553 Input Parameters: 1554 + A - the matrix 1555 - store - the boolean flag 1556 1557 Level: advanced 1558 1559 Notes: 1560 1561 .seealso: `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatPtAP()` 1562 @*/ 1563 PetscErrorCode MatISStoreL2L(Mat A, PetscBool store) 1564 { 1565 PetscFunctionBegin; 1566 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1567 PetscValidType(A,1); 1568 PetscValidLogicalCollectiveBool(A,store,2); 1569 PetscTryMethod(A,"MatISStoreL2L_C",(Mat,PetscBool),(A,store)); 1570 PetscFunctionReturn(0); 1571 } 1572 1573 static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store) 1574 { 1575 Mat_IS *matis = (Mat_IS*)(A->data); 1576 1577 PetscFunctionBegin; 1578 matis->storel2l = store; 1579 if (!store) { 1580 PetscCall(PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",NULL)); 1581 } 1582 PetscFunctionReturn(0); 1583 } 1584 1585 /*@ 1586 MatISFixLocalEmpty - Compress out zero local rows from the local matrices 1587 1588 Collective 1589 1590 Input Parameters: 1591 + A - the matrix 1592 - fix - the boolean flag 1593 1594 Level: advanced 1595 1596 Notes: When fix is true, new local matrices and l2g maps are generated during the final assembly process. 1597 1598 .seealso: `MatCreate()`, `MatCreateIS()`, `MatISSetPreallocation()`, `MatAssemblyEnd()`, `MAT_FINAL_ASSEMBLY` 1599 @*/ 1600 PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix) 1601 { 1602 PetscFunctionBegin; 1603 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1604 PetscValidType(A,1); 1605 PetscValidLogicalCollectiveBool(A,fix,2); 1606 PetscTryMethod(A,"MatISFixLocalEmpty_C",(Mat,PetscBool),(A,fix)); 1607 PetscFunctionReturn(0); 1608 } 1609 1610 static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix) 1611 { 1612 Mat_IS *matis = (Mat_IS*)(A->data); 1613 1614 PetscFunctionBegin; 1615 matis->locempty = fix; 1616 PetscFunctionReturn(0); 1617 } 1618 1619 /*@ 1620 MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix. 1621 1622 Collective 1623 1624 Input Parameters: 1625 + B - the matrix 1626 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1627 (same value is used for all local rows) 1628 . d_nnz - array containing the number of nonzeros in the various rows of the 1629 DIAGONAL portion of the local submatrix (possibly different for each row) 1630 or NULL, if d_nz is used to specify the nonzero structure. 1631 The size of this array is equal to the number of local rows, i.e 'm'. 1632 For matrices that will be factored, you must leave room for (and set) 1633 the diagonal entry even if it is zero. 1634 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1635 submatrix (same value is used for all local rows). 1636 - o_nnz - array containing the number of nonzeros in the various rows of the 1637 OFF-DIAGONAL portion of the local submatrix (possibly different for 1638 each row) or NULL, if o_nz is used to specify the nonzero 1639 structure. The size of this array is equal to the number 1640 of local rows, i.e 'm'. 1641 1642 If the *_nnz parameter is given then the *_nz parameter is ignored 1643 1644 Level: intermediate 1645 1646 Notes: 1647 This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition 1648 from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local 1649 matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects. 1650 1651 .seealso: `MatCreate()`, `MatCreateIS()`, `MatMPIAIJSetPreallocation()`, `MatISGetLocalMat()`, `MATIS` 1652 @*/ 1653 PetscErrorCode MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1654 { 1655 PetscFunctionBegin; 1656 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1657 PetscValidType(B,1); 1658 PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz)); 1659 PetscFunctionReturn(0); 1660 } 1661 1662 /* this is used by DMDA */ 1663 PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1664 { 1665 Mat_IS *matis = (Mat_IS*)(B->data); 1666 PetscInt bs,i,nlocalcols; 1667 1668 PetscFunctionBegin; 1669 PetscCall(MatSetUp(B)); 1670 if (!d_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nz; 1671 else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nnz[i]; 1672 1673 if (!o_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nz; 1674 else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nnz[i]; 1675 1676 PetscCall(PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE)); 1677 PetscCall(MatGetSize(matis->A,NULL,&nlocalcols)); 1678 PetscCall(MatGetBlockSize(matis->A,&bs)); 1679 PetscCall(PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE)); 1680 1681 for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols); 1682 PetscCall(MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata)); 1683 #if defined(PETSC_HAVE_HYPRE) 1684 PetscCall(MatHYPRESetPreallocation(matis->A,0,matis->sf_leafdata,0,NULL)); 1685 #endif 1686 1687 for (i=0;i<matis->sf->nleaves/bs;i++) { 1688 PetscInt b; 1689 1690 matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs; 1691 for (b=1;b<bs;b++) { 1692 matis->sf_leafdata[i] = PetscMax(matis->sf_leafdata[i],matis->sf_leafdata[i*bs+b]/bs); 1693 } 1694 } 1695 PetscCall(MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata)); 1696 1697 nlocalcols /= bs; 1698 for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols - i); 1699 PetscCall(MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata)); 1700 1701 /* for other matrix types */ 1702 PetscCall(MatSetUp(matis->A)); 1703 PetscFunctionReturn(0); 1704 } 1705 1706 PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce) 1707 { 1708 Mat_IS *matis = (Mat_IS*)(A->data); 1709 PetscInt *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership; 1710 const PetscInt *global_indices_r,*global_indices_c; 1711 PetscInt i,j,bs,rows,cols; 1712 PetscInt lrows,lcols; 1713 PetscInt local_rows,local_cols; 1714 PetscMPIInt size; 1715 PetscBool isdense,issbaij; 1716 1717 PetscFunctionBegin; 1718 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A),&size)); 1719 PetscCall(MatGetSize(A,&rows,&cols)); 1720 PetscCall(MatGetBlockSize(A,&bs)); 1721 PetscCall(MatGetSize(matis->A,&local_rows,&local_cols)); 1722 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense)); 1723 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij)); 1724 PetscCall(ISLocalToGlobalMappingGetIndices(matis->rmapping,&global_indices_r)); 1725 if (matis->rmapping != matis->cmapping) { 1726 PetscCall(ISLocalToGlobalMappingGetIndices(matis->cmapping,&global_indices_c)); 1727 } else global_indices_c = global_indices_r; 1728 1729 if (issbaij) PetscCall(MatGetRowUpperTriangular(matis->A)); 1730 /* 1731 An SF reduce is needed to sum up properly on shared rows. 1732 Note that generally preallocation is not exact, since it overestimates nonzeros 1733 */ 1734 PetscCall(MatGetLocalSize(A,&lrows,&lcols)); 1735 MatPreallocateBegin(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz); 1736 /* All processes need to compute entire row ownership */ 1737 PetscCall(PetscMalloc1(rows,&row_ownership)); 1738 PetscCall(MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges)); 1739 for (i=0;i<size;i++) { 1740 for (j=mat_ranges[i];j<mat_ranges[i+1];j++) row_ownership[j] = i; 1741 } 1742 PetscCall(MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges)); 1743 1744 /* 1745 my_dnz and my_onz contains exact contribution to preallocation from each local mat 1746 then, they will be summed up properly. This way, preallocation is always sufficient 1747 */ 1748 PetscCall(PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz)); 1749 /* preallocation as a MATAIJ */ 1750 if (isdense) { /* special case for dense local matrices */ 1751 for (i=0;i<local_rows;i++) { 1752 PetscInt owner = row_ownership[global_indices_r[i]]; 1753 for (j=0;j<local_cols;j++) { 1754 PetscInt index_col = global_indices_c[j]; 1755 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1]) { /* diag block */ 1756 my_dnz[i] += 1; 1757 } else { /* offdiag block */ 1758 my_onz[i] += 1; 1759 } 1760 } 1761 } 1762 } else if (matis->A->ops->getrowij) { 1763 const PetscInt *ii,*jj,*jptr; 1764 PetscBool done; 1765 PetscCall(MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done)); 1766 PetscCheck(done,PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ"); 1767 jptr = jj; 1768 for (i=0;i<local_rows;i++) { 1769 PetscInt index_row = global_indices_r[i]; 1770 for (j=0;j<ii[i+1]-ii[i];j++,jptr++) { 1771 PetscInt owner = row_ownership[index_row]; 1772 PetscInt index_col = global_indices_c[*jptr]; 1773 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1]) { /* diag block */ 1774 my_dnz[i] += 1; 1775 } else { /* offdiag block */ 1776 my_onz[i] += 1; 1777 } 1778 /* same as before, interchanging rows and cols */ 1779 if (issbaij && index_col != index_row) { 1780 owner = row_ownership[index_col]; 1781 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1]) { 1782 my_dnz[*jptr] += 1; 1783 } else { 1784 my_onz[*jptr] += 1; 1785 } 1786 } 1787 } 1788 } 1789 PetscCall(MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done)); 1790 PetscCheck(done,PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ"); 1791 } else { /* loop over rows and use MatGetRow */ 1792 for (i=0;i<local_rows;i++) { 1793 const PetscInt *cols; 1794 PetscInt ncols,index_row = global_indices_r[i]; 1795 PetscCall(MatGetRow(matis->A,i,&ncols,&cols,NULL)); 1796 for (j=0;j<ncols;j++) { 1797 PetscInt owner = row_ownership[index_row]; 1798 PetscInt index_col = global_indices_c[cols[j]]; 1799 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1]) { /* diag block */ 1800 my_dnz[i] += 1; 1801 } else { /* offdiag block */ 1802 my_onz[i] += 1; 1803 } 1804 /* same as before, interchanging rows and cols */ 1805 if (issbaij && index_col != index_row) { 1806 owner = row_ownership[index_col]; 1807 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1]) { 1808 my_dnz[cols[j]] += 1; 1809 } else { 1810 my_onz[cols[j]] += 1; 1811 } 1812 } 1813 } 1814 PetscCall(MatRestoreRow(matis->A,i,&ncols,&cols,NULL)); 1815 } 1816 } 1817 if (global_indices_c != global_indices_r) { 1818 PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->cmapping,&global_indices_c)); 1819 } 1820 PetscCall(ISLocalToGlobalMappingRestoreIndices(matis->rmapping,&global_indices_r)); 1821 PetscCall(PetscFree(row_ownership)); 1822 1823 /* Reduce my_dnz and my_onz */ 1824 if (maxreduce) { 1825 PetscCall(PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX)); 1826 PetscCall(PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX)); 1827 PetscCall(PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX)); 1828 PetscCall(PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX)); 1829 } else { 1830 PetscCall(PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM)); 1831 PetscCall(PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM)); 1832 PetscCall(PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM)); 1833 PetscCall(PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM)); 1834 } 1835 PetscCall(PetscFree2(my_dnz,my_onz)); 1836 1837 /* Resize preallocation if overestimated */ 1838 for (i=0;i<lrows;i++) { 1839 dnz[i] = PetscMin(dnz[i],lcols); 1840 onz[i] = PetscMin(onz[i],cols-lcols); 1841 } 1842 1843 /* Set preallocation */ 1844 PetscCall(MatSetBlockSizesFromMats(B,A,A)); 1845 PetscCall(MatSeqAIJSetPreallocation(B,0,dnz)); 1846 PetscCall(MatMPIAIJSetPreallocation(B,0,dnz,0,onz)); 1847 for (i=0;i<lrows;i+=bs) { 1848 PetscInt b, d = dnz[i],o = onz[i]; 1849 1850 for (b=1;b<bs;b++) { 1851 d = PetscMax(d,dnz[i+b]); 1852 o = PetscMax(o,onz[i+b]); 1853 } 1854 dnz[i/bs] = PetscMin(d/bs + d%bs,lcols/bs); 1855 onz[i/bs] = PetscMin(o/bs + o%bs,(cols-lcols)/bs); 1856 } 1857 PetscCall(MatSeqBAIJSetPreallocation(B,bs,0,dnz)); 1858 PetscCall(MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz)); 1859 PetscCall(MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz)); 1860 MatPreallocateEnd(dnz,onz); 1861 if (issbaij) PetscCall(MatRestoreRowUpperTriangular(matis->A)); 1862 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE)); 1863 PetscFunctionReturn(0); 1864 } 1865 1866 PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M) 1867 { 1868 Mat_IS *matis = (Mat_IS*)(mat->data); 1869 Mat local_mat,MT; 1870 PetscInt rbs,cbs,rows,cols,lrows,lcols; 1871 PetscInt local_rows,local_cols; 1872 PetscBool isseqdense,isseqsbaij,isseqaij,isseqbaij; 1873 PetscMPIInt size; 1874 const PetscScalar *array; 1875 1876 PetscFunctionBegin; 1877 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size)); 1878 if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) { 1879 Mat B; 1880 IS irows = NULL,icols = NULL; 1881 PetscInt rbs,cbs; 1882 1883 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping,&rbs)); 1884 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping,&cbs)); 1885 if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */ 1886 IS rows,cols; 1887 const PetscInt *ridxs,*cidxs; 1888 PetscInt i,nw,*work; 1889 1890 PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->rmapping,&ridxs)); 1891 PetscCall(ISLocalToGlobalMappingGetSize(matis->rmapping,&nw)); 1892 nw = nw/rbs; 1893 PetscCall(PetscCalloc1(nw,&work)); 1894 for (i=0;i<nw;i++) work[ridxs[i]] += 1; 1895 for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break; 1896 if (i == nw) { 1897 PetscCall(ISCreateBlock(PETSC_COMM_SELF,rbs,nw,ridxs,PETSC_USE_POINTER,&rows)); 1898 PetscCall(ISSetPermutation(rows)); 1899 PetscCall(ISInvertPermutation(rows,PETSC_DECIDE,&irows)); 1900 PetscCall(ISDestroy(&rows)); 1901 } 1902 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->rmapping,&ridxs)); 1903 PetscCall(PetscFree(work)); 1904 if (irows && matis->rmapping != matis->cmapping) { 1905 PetscCall(ISLocalToGlobalMappingGetBlockIndices(matis->cmapping,&cidxs)); 1906 PetscCall(ISLocalToGlobalMappingGetSize(matis->cmapping,&nw)); 1907 nw = nw/cbs; 1908 PetscCall(PetscCalloc1(nw,&work)); 1909 for (i=0;i<nw;i++) work[cidxs[i]] += 1; 1910 for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break; 1911 if (i == nw) { 1912 PetscCall(ISCreateBlock(PETSC_COMM_SELF,cbs,nw,cidxs,PETSC_USE_POINTER,&cols)); 1913 PetscCall(ISSetPermutation(cols)); 1914 PetscCall(ISInvertPermutation(cols,PETSC_DECIDE,&icols)); 1915 PetscCall(ISDestroy(&cols)); 1916 } 1917 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(matis->cmapping,&cidxs)); 1918 PetscCall(PetscFree(work)); 1919 } else if (irows) { 1920 PetscCall(PetscObjectReference((PetscObject)irows)); 1921 icols = irows; 1922 } 1923 } else { 1924 PetscCall(PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject*)&irows)); 1925 PetscCall(PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject*)&icols)); 1926 if (irows) PetscCall(PetscObjectReference((PetscObject)irows)); 1927 if (icols) PetscCall(PetscObjectReference((PetscObject)icols)); 1928 } 1929 if (!irows || !icols) { 1930 PetscCall(ISDestroy(&icols)); 1931 PetscCall(ISDestroy(&irows)); 1932 goto general_assembly; 1933 } 1934 PetscCall(MatConvert(matis->A,mtype,MAT_INITIAL_MATRIX,&B)); 1935 if (reuse != MAT_INPLACE_MATRIX) { 1936 PetscCall(MatCreateSubMatrix(B,irows,icols,reuse,M)); 1937 PetscCall(PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject)irows)); 1938 PetscCall(PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject)icols)); 1939 } else { 1940 Mat C; 1941 1942 PetscCall(MatCreateSubMatrix(B,irows,icols,MAT_INITIAL_MATRIX,&C)); 1943 PetscCall(MatHeaderReplace(mat,&C)); 1944 } 1945 PetscCall(MatDestroy(&B)); 1946 PetscCall(ISDestroy(&icols)); 1947 PetscCall(ISDestroy(&irows)); 1948 PetscFunctionReturn(0); 1949 } 1950 general_assembly: 1951 PetscCall(MatGetSize(mat,&rows,&cols)); 1952 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->rmapping,&rbs)); 1953 PetscCall(ISLocalToGlobalMappingGetBlockSize(matis->cmapping,&cbs)); 1954 PetscCall(MatGetLocalSize(mat,&lrows,&lcols)); 1955 PetscCall(MatGetSize(matis->A,&local_rows,&local_cols)); 1956 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense)); 1957 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij)); 1958 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij)); 1959 PetscCall(PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij)); 1960 PetscCheck(isseqdense || isseqaij || isseqbaij || isseqsbaij,PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name); 1961 if (PetscDefined (USE_DEBUG)) { 1962 PetscBool lb[4],bb[4]; 1963 1964 lb[0] = isseqdense; 1965 lb[1] = isseqaij; 1966 lb[2] = isseqbaij; 1967 lb[3] = isseqsbaij; 1968 PetscCall(MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat))); 1969 PetscCheck(bb[0] || bb[1] || bb[2] || bb[3],PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type"); 1970 } 1971 1972 if (reuse != MAT_REUSE_MATRIX) { 1973 PetscCall(MatCreate(PetscObjectComm((PetscObject)mat),&MT)); 1974 PetscCall(MatSetSizes(MT,lrows,lcols,rows,cols)); 1975 PetscCall(MatSetType(MT,mtype)); 1976 PetscCall(MatSetBlockSizes(MT,rbs,cbs)); 1977 PetscCall(MatISSetMPIXAIJPreallocation_Private(mat,MT,PETSC_FALSE)); 1978 } else { 1979 PetscInt mrbs,mcbs,mrows,mcols,mlrows,mlcols; 1980 1981 /* some checks */ 1982 MT = *M; 1983 PetscCall(MatGetBlockSizes(MT,&mrbs,&mcbs)); 1984 PetscCall(MatGetSize(MT,&mrows,&mcols)); 1985 PetscCall(MatGetLocalSize(MT,&mlrows,&mlcols)); 1986 PetscCheck(mrows == rows,PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%" PetscInt_FMT " != %" PetscInt_FMT ")",rows,mrows); 1987 PetscCheck(mcols == cols,PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%" PetscInt_FMT " != %" PetscInt_FMT ")",cols,mcols); 1988 PetscCheck(mlrows == lrows,PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%" PetscInt_FMT " != %" PetscInt_FMT ")",lrows,mlrows); 1989 PetscCheck(mlcols == lcols,PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%" PetscInt_FMT " != %" PetscInt_FMT ")",lcols,mlcols); 1990 PetscCheck(mrbs == rbs,PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong row block size (%" PetscInt_FMT " != %" PetscInt_FMT ")",rbs,mrbs); 1991 PetscCheck(mcbs == cbs,PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong col block size (%" PetscInt_FMT " != %" PetscInt_FMT ")",cbs,mcbs); 1992 PetscCall(MatZeroEntries(MT)); 1993 } 1994 1995 if (isseqsbaij || isseqbaij) { 1996 PetscCall(MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&local_mat)); 1997 isseqaij = PETSC_TRUE; 1998 } else { 1999 PetscCall(PetscObjectReference((PetscObject)matis->A)); 2000 local_mat = matis->A; 2001 } 2002 2003 /* Set values */ 2004 PetscCall(MatSetLocalToGlobalMapping(MT,matis->rmapping,matis->cmapping)); 2005 if (isseqdense) { /* special case for dense local matrices */ 2006 PetscInt i,*dummy; 2007 2008 PetscCall(PetscMalloc1(PetscMax(local_rows,local_cols),&dummy)); 2009 for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i; 2010 PetscCall(MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_FALSE)); 2011 PetscCall(MatDenseGetArrayRead(local_mat,&array)); 2012 PetscCall(MatSetValuesLocal(MT,local_rows,dummy,local_cols,dummy,array,ADD_VALUES)); 2013 PetscCall(MatDenseRestoreArrayRead(local_mat,&array)); 2014 PetscCall(PetscFree(dummy)); 2015 } else if (isseqaij) { 2016 const PetscInt *blocks; 2017 PetscInt i,nvtxs,*xadj,*adjncy, nb; 2018 PetscBool done; 2019 PetscScalar *sarray; 2020 2021 PetscCall(MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done)); 2022 PetscCheck(done,PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ"); 2023 PetscCall(MatSeqAIJGetArray(local_mat,&sarray)); 2024 PetscCall(MatGetVariableBlockSizes(local_mat,&nb,&blocks)); 2025 if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */ 2026 PetscInt sum; 2027 2028 for (i=0,sum=0;i<nb;i++) sum += blocks[i]; 2029 if (sum == nvtxs) { 2030 PetscInt r; 2031 2032 for (i=0,r=0;i<nb;i++) { 2033 PetscAssert(blocks[i] == xadj[r+1] - xadj[r],PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid block sizes prescribed for block %" PetscInt_FMT ": expected %" PetscInt_FMT ", got %" PetscInt_FMT,i,blocks[i],xadj[r+1] - xadj[r]); 2034 PetscCall(MatSetValuesLocal(MT,blocks[i],adjncy+xadj[r],blocks[i],adjncy+xadj[r],sarray+xadj[r],ADD_VALUES)); 2035 r += blocks[i]; 2036 } 2037 } else { 2038 for (i=0;i<nvtxs;i++) { 2039 PetscCall(MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],sarray+xadj[i],ADD_VALUES)); 2040 } 2041 } 2042 } else { 2043 for (i=0;i<nvtxs;i++) { 2044 PetscCall(MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],sarray+xadj[i],ADD_VALUES)); 2045 } 2046 } 2047 PetscCall(MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done)); 2048 PetscCheck(done,PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ"); 2049 PetscCall(MatSeqAIJRestoreArray(local_mat,&sarray)); 2050 } else { /* very basic values insertion for all other matrix types */ 2051 PetscInt i; 2052 2053 for (i=0;i<local_rows;i++) { 2054 PetscInt j; 2055 const PetscInt *local_indices_cols; 2056 2057 PetscCall(MatGetRow(local_mat,i,&j,&local_indices_cols,&array)); 2058 PetscCall(MatSetValuesLocal(MT,1,&i,j,local_indices_cols,array,ADD_VALUES)); 2059 PetscCall(MatRestoreRow(local_mat,i,&j,&local_indices_cols,&array)); 2060 } 2061 } 2062 PetscCall(MatAssemblyBegin(MT,MAT_FINAL_ASSEMBLY)); 2063 PetscCall(MatDestroy(&local_mat)); 2064 PetscCall(MatAssemblyEnd(MT,MAT_FINAL_ASSEMBLY)); 2065 if (isseqdense) { 2066 PetscCall(MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_TRUE)); 2067 } 2068 if (reuse == MAT_INPLACE_MATRIX) { 2069 PetscCall(MatHeaderReplace(mat,&MT)); 2070 } else if (reuse == MAT_INITIAL_MATRIX) { 2071 *M = MT; 2072 } 2073 PetscFunctionReturn(0); 2074 } 2075 2076 /*@ 2077 MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format 2078 2079 Input Parameters: 2080 + mat - the matrix (should be of type MATIS) 2081 - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 2082 2083 Output Parameter: 2084 . newmat - the matrix in AIJ format 2085 2086 Level: developer 2087 2088 Notes: 2089 This function has been deprecated and it will be removed in future releases. Update your code to use the MatConvert() interface. 2090 2091 .seealso: `MATIS`, `MatConvert()` 2092 @*/ 2093 PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat) 2094 { 2095 PetscFunctionBegin; 2096 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2097 PetscValidLogicalCollectiveEnum(mat,reuse,2); 2098 PetscValidPointer(newmat,3); 2099 if (reuse == MAT_REUSE_MATRIX) { 2100 PetscValidHeaderSpecific(*newmat,MAT_CLASSID,3); 2101 PetscCheckSameComm(mat,1,*newmat,3); 2102 PetscCheck(mat != *newmat,PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix"); 2103 } 2104 PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatType,MatReuse,Mat*),(mat,MATAIJ,reuse,newmat)); 2105 PetscFunctionReturn(0); 2106 } 2107 2108 static PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat) 2109 { 2110 Mat_IS *matis = (Mat_IS*)(mat->data); 2111 PetscInt rbs,cbs,m,n,M,N; 2112 Mat B,localmat; 2113 2114 PetscFunctionBegin; 2115 PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs)); 2116 PetscCall(ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs)); 2117 PetscCall(MatGetSize(mat,&M,&N)); 2118 PetscCall(MatGetLocalSize(mat,&m,&n)); 2119 PetscCall(MatCreate(PetscObjectComm((PetscObject)mat),&B)); 2120 PetscCall(MatSetSizes(B,m,n,M,N)); 2121 PetscCall(MatSetBlockSize(B,rbs == cbs ? rbs : 1)); 2122 PetscCall(MatSetType(B,MATIS)); 2123 PetscCall(MatISSetLocalMatType(B,matis->lmattype)); 2124 PetscCall(MatSetLocalToGlobalMapping(B,mat->rmap->mapping,mat->cmap->mapping)); 2125 PetscCall(MatDuplicate(matis->A,op,&localmat)); 2126 PetscCall(MatSetLocalToGlobalMapping(localmat,matis->A->rmap->mapping,matis->A->cmap->mapping)); 2127 PetscCall(MatISSetLocalMat(B,localmat)); 2128 PetscCall(MatDestroy(&localmat)); 2129 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 2130 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 2131 *newmat = B; 2132 PetscFunctionReturn(0); 2133 } 2134 2135 static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg) 2136 { 2137 Mat_IS *matis = (Mat_IS*)A->data; 2138 PetscBool local_sym; 2139 2140 PetscFunctionBegin; 2141 PetscCall(MatIsHermitian(matis->A,tol,&local_sym)); 2142 PetscCall(MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A))); 2143 PetscFunctionReturn(0); 2144 } 2145 2146 static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg) 2147 { 2148 Mat_IS *matis = (Mat_IS*)A->data; 2149 PetscBool local_sym; 2150 2151 PetscFunctionBegin; 2152 if (matis->rmapping != matis->cmapping) { 2153 *flg = PETSC_FALSE; 2154 PetscFunctionReturn(0); 2155 } 2156 PetscCall(MatIsSymmetric(matis->A,tol,&local_sym)); 2157 PetscCall(MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A))); 2158 PetscFunctionReturn(0); 2159 } 2160 2161 static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool *flg) 2162 { 2163 Mat_IS *matis = (Mat_IS*)A->data; 2164 PetscBool local_sym; 2165 2166 PetscFunctionBegin; 2167 if (matis->rmapping != matis->cmapping) { 2168 *flg = PETSC_FALSE; 2169 PetscFunctionReturn(0); 2170 } 2171 PetscCall(MatIsStructurallySymmetric(matis->A,&local_sym)); 2172 PetscCall(MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A))); 2173 PetscFunctionReturn(0); 2174 } 2175 2176 static PetscErrorCode MatDestroy_IS(Mat A) 2177 { 2178 Mat_IS *b = (Mat_IS*)A->data; 2179 2180 PetscFunctionBegin; 2181 PetscCall(PetscFree(b->bdiag)); 2182 PetscCall(PetscFree(b->lmattype)); 2183 PetscCall(MatDestroy(&b->A)); 2184 PetscCall(VecScatterDestroy(&b->cctx)); 2185 PetscCall(VecScatterDestroy(&b->rctx)); 2186 PetscCall(VecDestroy(&b->x)); 2187 PetscCall(VecDestroy(&b->y)); 2188 PetscCall(VecDestroy(&b->counter)); 2189 PetscCall(ISDestroy(&b->getsub_ris)); 2190 PetscCall(ISDestroy(&b->getsub_cis)); 2191 if (b->sf != b->csf) { 2192 PetscCall(PetscSFDestroy(&b->csf)); 2193 PetscCall(PetscFree2(b->csf_rootdata,b->csf_leafdata)); 2194 } else b->csf = NULL; 2195 PetscCall(PetscSFDestroy(&b->sf)); 2196 PetscCall(PetscFree2(b->sf_rootdata,b->sf_leafdata)); 2197 PetscCall(ISLocalToGlobalMappingDestroy(&b->rmapping)); 2198 PetscCall(ISLocalToGlobalMappingDestroy(&b->cmapping)); 2199 PetscCall(PetscFree(A->data)); 2200 PetscCall(PetscObjectChangeTypeName((PetscObject)A,NULL)); 2201 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",NULL)); 2202 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL)); 2203 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL)); 2204 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL)); 2205 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL)); 2206 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL)); 2207 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",NULL)); 2208 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalToGlobalMapping_C",NULL)); 2209 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL)); 2210 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL)); 2211 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL)); 2212 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL)); 2213 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL)); 2214 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL)); 2215 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL)); 2216 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOOLocal_C",NULL)); 2217 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",NULL)); 2218 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",NULL)); 2219 PetscFunctionReturn(0); 2220 } 2221 2222 static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y) 2223 { 2224 Mat_IS *is = (Mat_IS*)A->data; 2225 PetscScalar zero = 0.0; 2226 2227 PetscFunctionBegin; 2228 /* scatter the global vector x into the local work vector */ 2229 PetscCall(VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD)); 2230 PetscCall(VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD)); 2231 2232 /* multiply the local matrix */ 2233 PetscCall(MatMult(is->A,is->x,is->y)); 2234 2235 /* scatter product back into global memory */ 2236 PetscCall(VecSet(y,zero)); 2237 PetscCall(VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE)); 2238 PetscCall(VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE)); 2239 PetscFunctionReturn(0); 2240 } 2241 2242 static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 2243 { 2244 Vec temp_vec; 2245 2246 PetscFunctionBegin; /* v3 = v2 + A * v1.*/ 2247 if (v3 != v2) { 2248 PetscCall(MatMult(A,v1,v3)); 2249 PetscCall(VecAXPY(v3,1.0,v2)); 2250 } else { 2251 PetscCall(VecDuplicate(v2,&temp_vec)); 2252 PetscCall(MatMult(A,v1,temp_vec)); 2253 PetscCall(VecAXPY(temp_vec,1.0,v2)); 2254 PetscCall(VecCopy(temp_vec,v3)); 2255 PetscCall(VecDestroy(&temp_vec)); 2256 } 2257 PetscFunctionReturn(0); 2258 } 2259 2260 static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x) 2261 { 2262 Mat_IS *is = (Mat_IS*)A->data; 2263 2264 PetscFunctionBegin; 2265 /* scatter the global vector x into the local work vector */ 2266 PetscCall(VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD)); 2267 PetscCall(VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD)); 2268 2269 /* multiply the local matrix */ 2270 PetscCall(MatMultTranspose(is->A,is->y,is->x)); 2271 2272 /* scatter product back into global vector */ 2273 PetscCall(VecSet(x,0)); 2274 PetscCall(VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE)); 2275 PetscCall(VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE)); 2276 PetscFunctionReturn(0); 2277 } 2278 2279 static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 2280 { 2281 Vec temp_vec; 2282 2283 PetscFunctionBegin; /* v3 = v2 + A' * v1.*/ 2284 if (v3 != v2) { 2285 PetscCall(MatMultTranspose(A,v1,v3)); 2286 PetscCall(VecAXPY(v3,1.0,v2)); 2287 } else { 2288 PetscCall(VecDuplicate(v2,&temp_vec)); 2289 PetscCall(MatMultTranspose(A,v1,temp_vec)); 2290 PetscCall(VecAXPY(temp_vec,1.0,v2)); 2291 PetscCall(VecCopy(temp_vec,v3)); 2292 PetscCall(VecDestroy(&temp_vec)); 2293 } 2294 PetscFunctionReturn(0); 2295 } 2296 2297 static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer) 2298 { 2299 Mat_IS *a = (Mat_IS*)A->data; 2300 PetscViewer sviewer; 2301 PetscBool isascii,view = PETSC_TRUE; 2302 2303 PetscFunctionBegin; 2304 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); 2305 if (isascii) { 2306 PetscViewerFormat format; 2307 2308 PetscCall(PetscViewerGetFormat(viewer,&format)); 2309 if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE; 2310 } 2311 if (!view) PetscFunctionReturn(0); 2312 PetscCall(PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); 2313 PetscCall(MatView(a->A,sviewer)); 2314 PetscCall(PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); 2315 PetscCall(PetscViewerFlush(viewer)); 2316 PetscFunctionReturn(0); 2317 } 2318 2319 static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat,const PetscScalar **values) 2320 { 2321 Mat_IS *is = (Mat_IS*)mat->data; 2322 MPI_Datatype nodeType; 2323 const PetscScalar *lv; 2324 PetscInt bs; 2325 2326 PetscFunctionBegin; 2327 PetscCall(MatGetBlockSize(mat,&bs)); 2328 PetscCall(MatSetBlockSize(is->A,bs)); 2329 PetscCall(MatInvertBlockDiagonal(is->A,&lv)); 2330 if (!is->bdiag) { 2331 PetscCall(PetscMalloc1(bs*mat->rmap->n,&is->bdiag)); 2332 } 2333 PetscCallMPI(MPI_Type_contiguous(bs,MPIU_SCALAR,&nodeType)); 2334 PetscCallMPI(MPI_Type_commit(&nodeType)); 2335 PetscCall(PetscSFReduceBegin(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE)); 2336 PetscCall(PetscSFReduceEnd(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE)); 2337 PetscCallMPI(MPI_Type_free(&nodeType)); 2338 if (values) *values = is->bdiag; 2339 PetscFunctionReturn(0); 2340 } 2341 2342 static PetscErrorCode MatISSetUpScatters_Private(Mat A) 2343 { 2344 Vec cglobal,rglobal; 2345 IS from; 2346 Mat_IS *is = (Mat_IS*)A->data; 2347 PetscScalar sum; 2348 const PetscInt *garray; 2349 PetscInt nr,rbs,nc,cbs; 2350 VecType rtype; 2351 2352 PetscFunctionBegin; 2353 PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping,&nr)); 2354 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping,&rbs)); 2355 PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping,&nc)); 2356 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping,&cbs)); 2357 PetscCall(VecDestroy(&is->x)); 2358 PetscCall(VecDestroy(&is->y)); 2359 PetscCall(VecDestroy(&is->counter)); 2360 PetscCall(VecScatterDestroy(&is->rctx)); 2361 PetscCall(VecScatterDestroy(&is->cctx)); 2362 PetscCall(MatCreateVecs(is->A,&is->x,&is->y)); 2363 PetscCall(VecBindToCPU(is->y,PETSC_TRUE)); 2364 PetscCall(VecGetRootType_Private(is->y,&rtype)); 2365 PetscCall(PetscFree(A->defaultvectype)); 2366 PetscCall(PetscStrallocpy(rtype,&A->defaultvectype)); 2367 PetscCall(MatCreateVecs(A,&cglobal,&rglobal)); 2368 PetscCall(VecBindToCPU(rglobal,PETSC_TRUE)); 2369 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping,&garray)); 2370 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A),rbs,nr/rbs,garray,PETSC_USE_POINTER,&from)); 2371 PetscCall(VecScatterCreate(rglobal,from,is->y,NULL,&is->rctx)); 2372 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping,&garray)); 2373 PetscCall(ISDestroy(&from)); 2374 if (is->rmapping != is->cmapping) { 2375 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping,&garray)); 2376 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A),cbs,nc/cbs,garray,PETSC_USE_POINTER,&from)); 2377 PetscCall(VecScatterCreate(cglobal,from,is->x,NULL,&is->cctx)); 2378 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping,&garray)); 2379 PetscCall(ISDestroy(&from)); 2380 } else { 2381 PetscCall(PetscObjectReference((PetscObject)is->rctx)); 2382 is->cctx = is->rctx; 2383 } 2384 PetscCall(VecDestroy(&cglobal)); 2385 2386 /* interface counter vector (local) */ 2387 PetscCall(VecDuplicate(is->y,&is->counter)); 2388 PetscCall(VecBindToCPU(is->counter,PETSC_TRUE)); 2389 PetscCall(VecSet(is->y,1.)); 2390 PetscCall(VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE)); 2391 PetscCall(VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE)); 2392 PetscCall(VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD)); 2393 PetscCall(VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD)); 2394 PetscCall(VecBindToCPU(is->y,PETSC_FALSE)); 2395 PetscCall(VecBindToCPU(is->counter,PETSC_FALSE)); 2396 2397 /* special functions for block-diagonal matrices */ 2398 PetscCall(VecSum(rglobal,&sum)); 2399 A->ops->invertblockdiagonal = NULL; 2400 if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && is->rmapping == is->cmapping) A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS; 2401 PetscCall(VecDestroy(&rglobal)); 2402 2403 /* setup SF for general purpose shared indices based communications */ 2404 PetscCall(MatISSetUpSF_IS(A)); 2405 PetscFunctionReturn(0); 2406 } 2407 2408 static PetscErrorCode MatISFilterL2GMap(Mat A, ISLocalToGlobalMapping map, ISLocalToGlobalMapping *nmap, ISLocalToGlobalMapping *lmap) 2409 { 2410 IS is; 2411 ISLocalToGlobalMappingType l2gtype; 2412 const PetscInt *idxs; 2413 PetscHSetI ht; 2414 PetscInt *nidxs; 2415 PetscInt i,n,bs,c; 2416 PetscBool flg[] = {PETSC_FALSE,PETSC_FALSE}; 2417 2418 PetscFunctionBegin; 2419 PetscCall(ISLocalToGlobalMappingGetSize(map,&n)); 2420 PetscCall(ISLocalToGlobalMappingGetBlockSize(map,&bs)); 2421 PetscCall(ISLocalToGlobalMappingGetBlockIndices(map,&idxs)); 2422 PetscCall(PetscHSetICreate(&ht)); 2423 PetscCall(PetscMalloc1(n/bs,&nidxs)); 2424 for (i=0,c=0;i<n/bs;i++) { 2425 PetscBool missing; 2426 if (idxs[i] < 0) { flg[0] = PETSC_TRUE; continue; } 2427 PetscCall(PetscHSetIQueryAdd(ht,idxs[i],&missing)); 2428 if (!missing) flg[1] = PETSC_TRUE; 2429 else nidxs[c++] = idxs[i]; 2430 } 2431 PetscCall(PetscHSetIDestroy(&ht)); 2432 PetscCall(MPIU_Allreduce(MPI_IN_PLACE,flg,2,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A))); 2433 if (!flg[0] && !flg[1]) { /* Entries are all non negative and unique */ 2434 *nmap = NULL; 2435 *lmap = NULL; 2436 PetscCall(PetscFree(nidxs)); 2437 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map,&idxs)); 2438 PetscFunctionReturn(0); 2439 } 2440 2441 /* New l2g map without negative or repeated indices */ 2442 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A),bs,c,nidxs,PETSC_USE_POINTER,&is)); 2443 PetscCall(ISLocalToGlobalMappingCreateIS(is,nmap)); 2444 PetscCall(ISDestroy(&is)); 2445 PetscCall(ISLocalToGlobalMappingGetType(map,&l2gtype)); 2446 PetscCall(ISLocalToGlobalMappingSetType(*nmap,l2gtype)); 2447 2448 /* New local l2g map for repeated indices */ 2449 PetscCall(ISGlobalToLocalMappingApplyBlock(*nmap,IS_GTOLM_MASK,n/bs,idxs,NULL,nidxs)); 2450 PetscCall(ISCreateBlock(PETSC_COMM_SELF,bs,n/bs,nidxs,PETSC_USE_POINTER,&is)); 2451 PetscCall(ISLocalToGlobalMappingCreateIS(is,lmap)); 2452 PetscCall(ISDestroy(&is)); 2453 2454 PetscCall(PetscFree(nidxs)); 2455 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map,&idxs)); 2456 PetscFunctionReturn(0); 2457 } 2458 2459 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping) 2460 { 2461 Mat_IS *is = (Mat_IS*)A->data; 2462 ISLocalToGlobalMapping localrmapping = NULL, localcmapping = NULL; 2463 PetscBool cong, freem[] = { PETSC_FALSE, PETSC_FALSE }; 2464 PetscInt nr,rbs,nc,cbs; 2465 2466 PetscFunctionBegin; 2467 if (rmapping) PetscCheckSameComm(A,1,rmapping,2); 2468 if (cmapping) PetscCheckSameComm(A,1,cmapping,3); 2469 2470 PetscCall(ISLocalToGlobalMappingDestroy(&is->rmapping)); 2471 PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping)); 2472 PetscCall(PetscLayoutSetUp(A->rmap)); 2473 PetscCall(PetscLayoutSetUp(A->cmap)); 2474 PetscCall(MatHasCongruentLayouts(A,&cong)); 2475 2476 /* If NULL, local space matches global space */ 2477 if (!rmapping) { 2478 IS is; 2479 2480 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),A->rmap->N,0,1,&is)); 2481 PetscCall(ISLocalToGlobalMappingCreateIS(is,&rmapping)); 2482 if (A->rmap->bs > 0) PetscCall(ISLocalToGlobalMappingSetBlockSize(rmapping,A->rmap->bs)); 2483 PetscCall(ISDestroy(&is)); 2484 freem[0] = PETSC_TRUE; 2485 if (!cmapping && cong && A->rmap->bs == A->cmap->bs) cmapping = rmapping; 2486 } else if (!is->islocalref) { /* check if the l2g map has negative or repeated entries */ 2487 PetscCall(MatISFilterL2GMap(A,rmapping,&is->rmapping,&localrmapping)); 2488 if (rmapping == cmapping) { 2489 PetscCall(PetscObjectReference((PetscObject)is->rmapping)); 2490 is->cmapping = is->rmapping; 2491 PetscCall(PetscObjectReference((PetscObject)localrmapping)); 2492 localcmapping = localrmapping; 2493 } 2494 } 2495 if (!cmapping) { 2496 IS is; 2497 2498 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),A->cmap->N,0,1,&is)); 2499 PetscCall(ISLocalToGlobalMappingCreateIS(is,&cmapping)); 2500 if (A->cmap->bs > 0) PetscCall(ISLocalToGlobalMappingSetBlockSize(cmapping,A->cmap->bs)); 2501 PetscCall(ISDestroy(&is)); 2502 freem[1] = PETSC_TRUE; 2503 } else if (cmapping != rmapping && !is->islocalref) { /* check if the l2g map has negative or repeated entries */ 2504 PetscCall(MatISFilterL2GMap(A,cmapping,&is->cmapping,&localcmapping)); 2505 } 2506 if (!is->rmapping) { 2507 PetscCall(PetscObjectReference((PetscObject)rmapping)); 2508 is->rmapping = rmapping; 2509 } 2510 if (!is->cmapping) { 2511 PetscCall(PetscObjectReference((PetscObject)cmapping)); 2512 is->cmapping = cmapping; 2513 } 2514 2515 /* Clean up */ 2516 PetscCall(MatDestroy(&is->A)); 2517 if (is->csf != is->sf) { 2518 PetscCall(PetscSFDestroy(&is->csf)); 2519 PetscCall(PetscFree2(is->csf_rootdata,is->csf_leafdata)); 2520 } else is->csf = NULL; 2521 PetscCall(PetscSFDestroy(&is->sf)); 2522 PetscCall(PetscFree2(is->sf_rootdata,is->sf_leafdata)); 2523 PetscCall(PetscFree(is->bdiag)); 2524 2525 /* check if the two mappings are actually the same for square matrices since MATIS has some optimization for this case 2526 (DOLFIN passes 2 different objects) */ 2527 PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping,&nr)); 2528 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping,&rbs)); 2529 PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping,&nc)); 2530 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping,&cbs)); 2531 if (is->rmapping != is->cmapping && cong) { 2532 PetscBool same = PETSC_FALSE; 2533 if (nr == nc && cbs == rbs) { 2534 const PetscInt *idxs1,*idxs2; 2535 2536 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping,&idxs1)); 2537 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping,&idxs2)); 2538 PetscCall(PetscArraycmp(idxs1,idxs2,nr/rbs,&same)); 2539 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping,&idxs1)); 2540 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping,&idxs2)); 2541 } 2542 PetscCall(MPIU_Allreduce(MPI_IN_PLACE,&same,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A))); 2543 if (same) { 2544 PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping)); 2545 PetscCall(PetscObjectReference((PetscObject)is->rmapping)); 2546 is->cmapping = is->rmapping; 2547 } 2548 } 2549 PetscCall(PetscLayoutSetBlockSize(A->rmap,rbs)); 2550 PetscCall(PetscLayoutSetBlockSize(A->cmap,cbs)); 2551 /* Pass the user defined maps to the layout */ 2552 PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping)); 2553 PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping)); 2554 if (freem[0]) PetscCall(ISLocalToGlobalMappingDestroy(&rmapping)); 2555 if (freem[1]) PetscCall(ISLocalToGlobalMappingDestroy(&cmapping)); 2556 2557 /* Create the local matrix A */ 2558 PetscCall(MatCreate(PETSC_COMM_SELF,&is->A)); 2559 PetscCall(MatSetType(is->A,is->lmattype)); 2560 PetscCall(MatSetSizes(is->A,nr,nc,nr,nc)); 2561 PetscCall(MatSetBlockSizes(is->A,rbs,cbs)); 2562 PetscCall(MatSetOptionsPrefix(is->A,"is_")); 2563 PetscCall(MatAppendOptionsPrefix(is->A,((PetscObject)A)->prefix)); 2564 PetscCall(PetscLayoutSetUp(is->A->rmap)); 2565 PetscCall(PetscLayoutSetUp(is->A->cmap)); 2566 PetscCall(MatSetLocalToGlobalMapping(is->A,localrmapping,localcmapping)); 2567 PetscCall(ISLocalToGlobalMappingDestroy(&localrmapping)); 2568 PetscCall(ISLocalToGlobalMappingDestroy(&localcmapping)); 2569 2570 /* setup scatters and local vectors for MatMult */ 2571 if (!is->islocalref) PetscCall(MatISSetUpScatters_Private(A)); 2572 A->preallocated = PETSC_TRUE; 2573 PetscFunctionReturn(0); 2574 } 2575 2576 static PetscErrorCode MatSetUp_IS(Mat A) 2577 { 2578 ISLocalToGlobalMapping rmap, cmap; 2579 2580 PetscFunctionBegin; 2581 PetscCall(MatGetLocalToGlobalMapping(A,&rmap,&cmap)); 2582 if (!rmap && !cmap) { 2583 PetscCall(MatSetLocalToGlobalMapping(A,NULL,NULL)); 2584 } 2585 PetscFunctionReturn(0); 2586 } 2587 2588 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2589 { 2590 Mat_IS *is = (Mat_IS*)mat->data; 2591 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2592 2593 PetscFunctionBegin; 2594 PetscCall(ISGlobalToLocalMappingApply(is->rmapping,IS_GTOLM_MASK,m,rows,&m,rows_l)); 2595 if (m != n || rows != cols || is->cmapping != is->rmapping) { 2596 PetscCall(ISGlobalToLocalMappingApply(is->cmapping,IS_GTOLM_MASK,n,cols,&n,cols_l)); 2597 PetscCall(MatSetValues(is->A,m,rows_l,n,cols_l,values,addv)); 2598 } else { 2599 PetscCall(MatSetValues(is->A,m,rows_l,m,rows_l,values,addv)); 2600 } 2601 PetscFunctionReturn(0); 2602 } 2603 2604 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2605 { 2606 Mat_IS *is = (Mat_IS*)mat->data; 2607 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2608 2609 PetscFunctionBegin; 2610 PetscCall(ISGlobalToLocalMappingApplyBlock(is->rmapping,IS_GTOLM_MASK,m,rows,&m,rows_l)); 2611 if (m != n || rows != cols || is->cmapping != is->rmapping) { 2612 PetscCall(ISGlobalToLocalMappingApply(is->cmapping,IS_GTOLM_MASK,n,cols,&n,cols_l)); 2613 PetscCall(MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv)); 2614 } else { 2615 PetscCall(MatSetValuesBlocked(is->A,m,rows_l,n,rows_l,values,addv)); 2616 } 2617 PetscFunctionReturn(0); 2618 } 2619 2620 static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2621 { 2622 Mat_IS *is = (Mat_IS*)A->data; 2623 2624 PetscFunctionBegin; 2625 if (is->A->rmap->mapping || is->A->cmap->mapping) { 2626 PetscCall(MatSetValuesLocal(is->A,m,rows,n,cols,values,addv)); 2627 } else { 2628 PetscCall(MatSetValues(is->A,m,rows,n,cols,values,addv)); 2629 } 2630 PetscFunctionReturn(0); 2631 } 2632 2633 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2634 { 2635 Mat_IS *is = (Mat_IS*)A->data; 2636 2637 PetscFunctionBegin; 2638 if (is->A->rmap->mapping || is->A->cmap->mapping) { 2639 PetscCall(MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv)); 2640 } else { 2641 PetscCall(MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv)); 2642 } 2643 PetscFunctionReturn(0); 2644 } 2645 2646 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns) 2647 { 2648 Mat_IS *is = (Mat_IS*)A->data; 2649 2650 PetscFunctionBegin; 2651 if (!n) { 2652 is->pure_neumann = PETSC_TRUE; 2653 } else { 2654 PetscInt i; 2655 is->pure_neumann = PETSC_FALSE; 2656 2657 if (columns) { 2658 PetscCall(MatZeroRowsColumns(is->A,n,rows,diag,NULL,NULL)); 2659 } else { 2660 PetscCall(MatZeroRows(is->A,n,rows,diag,NULL,NULL)); 2661 } 2662 if (diag != 0.) { 2663 const PetscScalar *array; 2664 PetscCall(VecGetArrayRead(is->counter,&array)); 2665 for (i=0; i<n; i++) { 2666 PetscCall(MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES)); 2667 } 2668 PetscCall(VecRestoreArrayRead(is->counter,&array)); 2669 } 2670 PetscCall(MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY)); 2671 PetscCall(MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY)); 2672 } 2673 PetscFunctionReturn(0); 2674 } 2675 2676 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns) 2677 { 2678 Mat_IS *matis = (Mat_IS*)A->data; 2679 PetscInt nr,nl,len,i; 2680 PetscInt *lrows; 2681 2682 PetscFunctionBegin; 2683 if (PetscUnlikelyDebug(columns || diag != 0. || (x && b))) { 2684 PetscBool cong; 2685 2686 PetscCall(PetscLayoutCompare(A->rmap,A->cmap,&cong)); 2687 cong = (PetscBool)(cong && matis->sf == matis->csf); 2688 PetscCheck(cong || !columns,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Columns can be zeroed if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS"); 2689 PetscCheck(cong || diag == 0.,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Nonzero diagonal value supported if and only if A->rmap and A->cmap are congruent and the l2g maps are the same for MATIS"); 2690 PetscCheck(cong || !x || !b,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"A->rmap and A->cmap need to be congruent, and the l2g maps be the same"); 2691 } 2692 /* get locally owned rows */ 2693 PetscCall(PetscLayoutMapLocal(A->rmap,n,rows,&len,&lrows,NULL)); 2694 /* fix right hand side if needed */ 2695 if (x && b) { 2696 const PetscScalar *xx; 2697 PetscScalar *bb; 2698 2699 PetscCall(VecGetArrayRead(x, &xx)); 2700 PetscCall(VecGetArray(b, &bb)); 2701 for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]]; 2702 PetscCall(VecRestoreArrayRead(x, &xx)); 2703 PetscCall(VecRestoreArray(b, &bb)); 2704 } 2705 /* get rows associated to the local matrices */ 2706 PetscCall(MatGetSize(matis->A,&nl,NULL)); 2707 PetscCall(PetscArrayzero(matis->sf_leafdata,nl)); 2708 PetscCall(PetscArrayzero(matis->sf_rootdata,A->rmap->n)); 2709 for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1; 2710 PetscCall(PetscFree(lrows)); 2711 PetscCall(PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE)); 2712 PetscCall(PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE)); 2713 PetscCall(PetscMalloc1(nl,&lrows)); 2714 for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i; 2715 PetscCall(MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns)); 2716 PetscCall(PetscFree(lrows)); 2717 PetscFunctionReturn(0); 2718 } 2719 2720 static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2721 { 2722 PetscFunctionBegin; 2723 PetscCall(MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE)); 2724 PetscFunctionReturn(0); 2725 } 2726 2727 static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2728 { 2729 PetscFunctionBegin; 2730 PetscCall(MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE)); 2731 PetscFunctionReturn(0); 2732 } 2733 2734 static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type) 2735 { 2736 Mat_IS *is = (Mat_IS*)A->data; 2737 2738 PetscFunctionBegin; 2739 PetscCall(MatAssemblyBegin(is->A,type)); 2740 PetscFunctionReturn(0); 2741 } 2742 2743 static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type) 2744 { 2745 Mat_IS *is = (Mat_IS*)A->data; 2746 2747 PetscFunctionBegin; 2748 PetscCall(MatAssemblyEnd(is->A,type)); 2749 /* fix for local empty rows/cols */ 2750 if (is->locempty && type == MAT_FINAL_ASSEMBLY) { 2751 Mat newlA; 2752 ISLocalToGlobalMapping rl2g,cl2g; 2753 IS nzr,nzc; 2754 PetscInt nr,nc,nnzr,nnzc; 2755 PetscBool lnewl2g,newl2g; 2756 2757 PetscCall(MatGetSize(is->A,&nr,&nc)); 2758 PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_FALSE,PETSC_SMALL,&nzr)); 2759 if (!nzr) { 2760 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&nzr)); 2761 } 2762 PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_TRUE,PETSC_SMALL,&nzc)); 2763 if (!nzc) { 2764 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&nzc)); 2765 } 2766 PetscCall(ISGetSize(nzr,&nnzr)); 2767 PetscCall(ISGetSize(nzc,&nnzc)); 2768 if (nnzr != nr || nnzc != nc) { /* need new global l2g map */ 2769 lnewl2g = PETSC_TRUE; 2770 PetscCallMPI(MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A))); 2771 2772 /* extract valid submatrix */ 2773 PetscCall(MatCreateSubMatrix(is->A,nzr,nzc,MAT_INITIAL_MATRIX,&newlA)); 2774 } else { /* local matrix fully populated */ 2775 lnewl2g = PETSC_FALSE; 2776 PetscCallMPI(MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A))); 2777 PetscCall(PetscObjectReference((PetscObject)is->A)); 2778 newlA = is->A; 2779 } 2780 2781 /* attach new global l2g map if needed */ 2782 if (newl2g) { 2783 IS zr,zc; 2784 const PetscInt *ridxs,*cidxs,*zridxs,*zcidxs; 2785 PetscInt *nidxs,i; 2786 2787 PetscCall(ISComplement(nzr,0,nr,&zr)); 2788 PetscCall(ISComplement(nzc,0,nc,&zc)); 2789 PetscCall(PetscMalloc1(PetscMax(nr,nc),&nidxs)); 2790 PetscCall(ISLocalToGlobalMappingGetIndices(is->rmapping,&ridxs)); 2791 PetscCall(ISLocalToGlobalMappingGetIndices(is->cmapping,&cidxs)); 2792 PetscCall(ISGetIndices(zr,&zridxs)); 2793 PetscCall(ISGetIndices(zc,&zcidxs)); 2794 PetscCall(ISGetLocalSize(zr,&nnzr)); 2795 PetscCall(ISGetLocalSize(zc,&nnzc)); 2796 2797 PetscCall(PetscArraycpy(nidxs,ridxs,nr)); 2798 for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1; 2799 PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nr,nidxs,PETSC_COPY_VALUES,&rl2g)); 2800 PetscCall(PetscArraycpy(nidxs,cidxs,nc)); 2801 for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1; 2802 PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nc,nidxs,PETSC_COPY_VALUES,&cl2g)); 2803 2804 PetscCall(ISRestoreIndices(zr,&zridxs)); 2805 PetscCall(ISRestoreIndices(zc,&zcidxs)); 2806 PetscCall(ISLocalToGlobalMappingRestoreIndices(is->rmapping,&ridxs)); 2807 PetscCall(ISLocalToGlobalMappingRestoreIndices(is->cmapping,&cidxs)); 2808 PetscCall(ISDestroy(&nzr)); 2809 PetscCall(ISDestroy(&nzc)); 2810 PetscCall(ISDestroy(&zr)); 2811 PetscCall(ISDestroy(&zc)); 2812 PetscCall(PetscFree(nidxs)); 2813 PetscCall(MatSetLocalToGlobalMapping(A,rl2g,cl2g)); 2814 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 2815 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 2816 } 2817 PetscCall(MatISSetLocalMat(A,newlA)); 2818 PetscCall(MatDestroy(&newlA)); 2819 PetscCall(ISDestroy(&nzr)); 2820 PetscCall(ISDestroy(&nzc)); 2821 is->locempty = PETSC_FALSE; 2822 } 2823 PetscFunctionReturn(0); 2824 } 2825 2826 static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local) 2827 { 2828 Mat_IS *is = (Mat_IS*)mat->data; 2829 2830 PetscFunctionBegin; 2831 *local = is->A; 2832 PetscFunctionReturn(0); 2833 } 2834 2835 static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local) 2836 { 2837 PetscFunctionBegin; 2838 *local = NULL; 2839 PetscFunctionReturn(0); 2840 } 2841 2842 /*@ 2843 MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix. 2844 2845 Input Parameter: 2846 . mat - the matrix 2847 2848 Output Parameter: 2849 . local - the local matrix 2850 2851 Level: advanced 2852 2853 Notes: 2854 This can be called if you have precomputed the nonzero structure of the 2855 matrix and want to provide it to the inner matrix object to improve the performance 2856 of the MatSetValues() operation. 2857 2858 Call MatISRestoreLocalMat() when finished with the local matrix. 2859 2860 .seealso: `MATIS` 2861 @*/ 2862 PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local) 2863 { 2864 PetscFunctionBegin; 2865 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2866 PetscValidPointer(local,2); 2867 PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local)); 2868 PetscFunctionReturn(0); 2869 } 2870 2871 /*@ 2872 MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat() 2873 2874 Input Parameter: 2875 . mat - the matrix 2876 2877 Output Parameter: 2878 . local - the local matrix 2879 2880 Level: advanced 2881 2882 .seealso: `MATIS` 2883 @*/ 2884 PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local) 2885 { 2886 PetscFunctionBegin; 2887 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2888 PetscValidPointer(local,2); 2889 PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local)); 2890 PetscFunctionReturn(0); 2891 } 2892 2893 static PetscErrorCode MatISSetLocalMatType_IS(Mat mat,MatType mtype) 2894 { 2895 Mat_IS *is = (Mat_IS*)mat->data; 2896 2897 PetscFunctionBegin; 2898 if (is->A) { 2899 PetscCall(MatSetType(is->A,mtype)); 2900 } 2901 PetscCall(PetscFree(is->lmattype)); 2902 PetscCall(PetscStrallocpy(mtype,&is->lmattype)); 2903 PetscFunctionReturn(0); 2904 } 2905 2906 /*@ 2907 MatISSetLocalMatType - Specifies the type of local matrix 2908 2909 Input Parameters: 2910 + mat - the matrix 2911 - mtype - the local matrix type 2912 2913 Output Parameter: 2914 2915 Level: advanced 2916 2917 .seealso: `MATIS`, `MatSetType()`, `MatType` 2918 @*/ 2919 PetscErrorCode MatISSetLocalMatType(Mat mat,MatType mtype) 2920 { 2921 PetscFunctionBegin; 2922 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2923 PetscUseMethod(mat,"MatISSetLocalMatType_C",(Mat,MatType),(mat,mtype)); 2924 PetscFunctionReturn(0); 2925 } 2926 2927 static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local) 2928 { 2929 Mat_IS *is = (Mat_IS*)mat->data; 2930 PetscInt nrows,ncols,orows,ocols; 2931 MatType mtype,otype; 2932 PetscBool sametype = PETSC_TRUE; 2933 2934 PetscFunctionBegin; 2935 if (is->A && !is->islocalref) { 2936 PetscCall(MatGetSize(is->A,&orows,&ocols)); 2937 PetscCall(MatGetSize(local,&nrows,&ncols)); 2938 PetscCheck(orows == nrows && ocols == ncols,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local MATIS matrix should be of size %" PetscInt_FMT "x%" PetscInt_FMT " (you passed a %" PetscInt_FMT "x%" PetscInt_FMT " matrix)",orows,ocols,nrows,ncols); 2939 PetscCall(MatGetType(local,&mtype)); 2940 PetscCall(MatGetType(is->A,&otype)); 2941 PetscCall(PetscStrcmp(mtype,otype,&sametype)); 2942 } 2943 PetscCall(PetscObjectReference((PetscObject)local)); 2944 PetscCall(MatDestroy(&is->A)); 2945 is->A = local; 2946 PetscCall(MatGetType(is->A,&mtype)); 2947 PetscCall(MatISSetLocalMatType(mat,mtype)); 2948 if (!sametype && !is->islocalref) { 2949 PetscCall(MatISSetUpScatters_Private(mat)); 2950 } 2951 PetscFunctionReturn(0); 2952 } 2953 2954 /*@ 2955 MatISSetLocalMat - Replace the local matrix stored inside a MATIS object. 2956 2957 Collective on Mat 2958 2959 Input Parameters: 2960 + mat - the matrix 2961 - local - the local matrix 2962 2963 Output Parameter: 2964 2965 Level: advanced 2966 2967 Notes: 2968 This can be called if you have precomputed the local matrix and 2969 want to provide it to the matrix object MATIS. 2970 2971 .seealso: `MATIS` 2972 @*/ 2973 PetscErrorCode MatISSetLocalMat(Mat mat,Mat local) 2974 { 2975 PetscFunctionBegin; 2976 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2977 PetscValidHeaderSpecific(local,MAT_CLASSID,2); 2978 PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local)); 2979 PetscFunctionReturn(0); 2980 } 2981 2982 static PetscErrorCode MatZeroEntries_IS(Mat A) 2983 { 2984 Mat_IS *a = (Mat_IS*)A->data; 2985 2986 PetscFunctionBegin; 2987 PetscCall(MatZeroEntries(a->A)); 2988 PetscFunctionReturn(0); 2989 } 2990 2991 static PetscErrorCode MatScale_IS(Mat A,PetscScalar a) 2992 { 2993 Mat_IS *is = (Mat_IS*)A->data; 2994 2995 PetscFunctionBegin; 2996 PetscCall(MatScale(is->A,a)); 2997 PetscFunctionReturn(0); 2998 } 2999 3000 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v) 3001 { 3002 Mat_IS *is = (Mat_IS*)A->data; 3003 3004 PetscFunctionBegin; 3005 /* get diagonal of the local matrix */ 3006 PetscCall(MatGetDiagonal(is->A,is->y)); 3007 3008 /* scatter diagonal back into global vector */ 3009 PetscCall(VecSet(v,0)); 3010 PetscCall(VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE)); 3011 PetscCall(VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE)); 3012 PetscFunctionReturn(0); 3013 } 3014 3015 static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg) 3016 { 3017 Mat_IS *a = (Mat_IS*)A->data; 3018 3019 PetscFunctionBegin; 3020 PetscCall(MatSetOption(a->A,op,flg)); 3021 PetscFunctionReturn(0); 3022 } 3023 3024 static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str) 3025 { 3026 Mat_IS *y = (Mat_IS*)Y->data; 3027 Mat_IS *x; 3028 3029 PetscFunctionBegin; 3030 if (PetscDefined(USE_DEBUG)) { 3031 PetscBool ismatis; 3032 PetscCall(PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis)); 3033 PetscCheck(ismatis,PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS"); 3034 } 3035 x = (Mat_IS*)X->data; 3036 PetscCall(MatAXPY(y->A,a,x->A,str)); 3037 PetscFunctionReturn(0); 3038 } 3039 3040 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat) 3041 { 3042 Mat lA; 3043 Mat_IS *matis = (Mat_IS*)(A->data); 3044 ISLocalToGlobalMapping rl2g,cl2g; 3045 IS is; 3046 const PetscInt *rg,*rl; 3047 PetscInt nrg; 3048 PetscInt N,M,nrl,i,*idxs; 3049 3050 PetscFunctionBegin; 3051 PetscCall(ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg)); 3052 PetscCall(ISGetLocalSize(row,&nrl)); 3053 PetscCall(ISGetIndices(row,&rl)); 3054 PetscCall(ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg)); 3055 if (PetscDefined(USE_DEBUG)) { 3056 for (i=0; i<nrl; i++) PetscCheck(rl[i]<nrg,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row index %" PetscInt_FMT " -> %" PetscInt_FMT " greater then maximum possible %" PetscInt_FMT,i,rl[i],nrg); 3057 } 3058 PetscCall(PetscMalloc1(nrg,&idxs)); 3059 /* map from [0,nrl) to row */ 3060 for (i=0;i<nrl;i++) idxs[i] = rl[i]; 3061 for (i=nrl;i<nrg;i++) idxs[i] = -1; 3062 PetscCall(ISRestoreIndices(row,&rl)); 3063 PetscCall(ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg)); 3064 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is)); 3065 PetscCall(ISLocalToGlobalMappingCreateIS(is,&rl2g)); 3066 PetscCall(ISDestroy(&is)); 3067 /* compute new l2g map for columns */ 3068 if (col != row || matis->rmapping != matis->cmapping || matis->A->rmap->mapping != matis->A->cmap->mapping) { 3069 const PetscInt *cg,*cl; 3070 PetscInt ncg; 3071 PetscInt ncl; 3072 3073 PetscCall(ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg)); 3074 PetscCall(ISGetLocalSize(col,&ncl)); 3075 PetscCall(ISGetIndices(col,&cl)); 3076 PetscCall(ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg)); 3077 if (PetscDefined(USE_DEBUG)) { 3078 for (i=0; i<ncl; i++) PetscCheck(cl[i]<ncg,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local column index %" PetscInt_FMT " -> %" PetscInt_FMT " greater then maximum possible %" PetscInt_FMT,i,cl[i],ncg); 3079 } 3080 PetscCall(PetscMalloc1(ncg,&idxs)); 3081 /* map from [0,ncl) to col */ 3082 for (i=0;i<ncl;i++) idxs[i] = cl[i]; 3083 for (i=ncl;i<ncg;i++) idxs[i] = -1; 3084 PetscCall(ISRestoreIndices(col,&cl)); 3085 PetscCall(ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg)); 3086 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is)); 3087 PetscCall(ISLocalToGlobalMappingCreateIS(is,&cl2g)); 3088 PetscCall(ISDestroy(&is)); 3089 } else { 3090 PetscCall(PetscObjectReference((PetscObject)rl2g)); 3091 cl2g = rl2g; 3092 } 3093 /* create the MATIS submatrix */ 3094 PetscCall(MatGetSize(A,&M,&N)); 3095 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),submat)); 3096 PetscCall(MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N)); 3097 PetscCall(MatSetType(*submat,MATIS)); 3098 matis = (Mat_IS*)((*submat)->data); 3099 matis->islocalref = PETSC_TRUE; 3100 PetscCall(MatSetLocalToGlobalMapping(*submat,rl2g,cl2g)); 3101 PetscCall(MatISGetLocalMat(A,&lA)); 3102 PetscCall(MatISSetLocalMat(*submat,lA)); 3103 PetscCall(MatSetUp(*submat)); 3104 PetscCall(MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY)); 3105 PetscCall(MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY)); 3106 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 3107 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 3108 3109 /* remove unsupported ops */ 3110 PetscCall(PetscMemzero((*submat)->ops,sizeof(struct _MatOps))); 3111 (*submat)->ops->destroy = MatDestroy_IS; 3112 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS; 3113 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS; 3114 (*submat)->ops->assemblybegin = MatAssemblyBegin_IS; 3115 (*submat)->ops->assemblyend = MatAssemblyEnd_IS; 3116 PetscFunctionReturn(0); 3117 } 3118 3119 static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A) 3120 { 3121 Mat_IS *a = (Mat_IS*)A->data; 3122 char type[256]; 3123 PetscBool flg; 3124 3125 PetscFunctionBegin; 3126 PetscOptionsHeadBegin(PetscOptionsObject,"MATIS options"); 3127 PetscCall(PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns","MatISFixLocalEmpty",a->locempty,&a->locempty,NULL)); 3128 PetscCall(PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL)); 3129 PetscCall(PetscOptionsFList("-matis_localmat_type","Matrix type","MatISSetLocalMatType",MatList,a->lmattype,type,256,&flg)); 3130 if (flg) { 3131 PetscCall(MatISSetLocalMatType(A,type)); 3132 } 3133 if (a->A) { 3134 PetscCall(MatSetFromOptions(a->A)); 3135 } 3136 PetscOptionsHeadEnd(); 3137 PetscFunctionReturn(0); 3138 } 3139 3140 /*@ 3141 MatCreateIS - Creates a "process" unassembled matrix, assembled on each 3142 process but not across processes. 3143 3144 Input Parameters: 3145 + comm - MPI communicator that will share the matrix 3146 . bs - block size of the matrix 3147 . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products 3148 . rmap - local to global map for rows 3149 - cmap - local to global map for cols 3150 3151 Output Parameter: 3152 . A - the resulting matrix 3153 3154 Level: advanced 3155 3156 Notes: 3157 See MATIS for more details. 3158 m and n are NOT related to the size of the map; they represent the size of the local parts of the distributed vectors 3159 used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices. 3160 If rmap (cmap) is NULL, then the local row (column) spaces matches the global space. 3161 3162 .seealso: `MATIS`, `MatSetLocalToGlobalMapping()` 3163 @*/ 3164 PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A) 3165 { 3166 PetscFunctionBegin; 3167 PetscCall(MatCreate(comm,A)); 3168 PetscCall(MatSetSizes(*A,m,n,M,N)); 3169 if (bs > 0) { 3170 PetscCall(MatSetBlockSize(*A,bs)); 3171 } 3172 PetscCall(MatSetType(*A,MATIS)); 3173 PetscCall(MatSetLocalToGlobalMapping(*A,rmap,cmap)); 3174 PetscFunctionReturn(0); 3175 } 3176 3177 static PetscErrorCode MatHasOperation_IS(Mat A, MatOperation op, PetscBool *has) 3178 { 3179 Mat_IS *a = (Mat_IS*)A->data; 3180 3181 PetscFunctionBegin; 3182 *has = PETSC_FALSE; 3183 if (!((void**)A->ops)[op]) PetscFunctionReturn(0); 3184 PetscCall(MatHasOperation(a->A,op,has)); 3185 PetscFunctionReturn(0); 3186 } 3187 3188 static PetscErrorCode MatSetValuesCOO_IS(Mat A,const PetscScalar v[],InsertMode imode) 3189 { 3190 Mat_IS *a = (Mat_IS*)A->data; 3191 3192 PetscFunctionBegin; 3193 PetscCall(MatSetValuesCOO(a->A,v,imode)); 3194 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 3195 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 3196 PetscFunctionReturn(0); 3197 } 3198 3199 static PetscErrorCode MatSetPreallocationCOOLocal_IS(Mat A,PetscCount ncoo,PetscInt coo_i[],PetscInt coo_j[]) 3200 { 3201 Mat_IS *a = (Mat_IS*)A->data; 3202 3203 PetscFunctionBegin; 3204 PetscCheck(a->A,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Need to provide l2g map first via MatSetLocalToGlobalMapping"); 3205 if (a->A->rmap->mapping || a->A->cmap->mapping) { 3206 PetscCall(MatSetPreallocationCOOLocal(a->A,ncoo,coo_i,coo_j)); 3207 } else { 3208 PetscCall(MatSetPreallocationCOO(a->A,ncoo,coo_i,coo_j)); 3209 } 3210 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",MatSetValuesCOO_IS)); 3211 A->preallocated = PETSC_TRUE; 3212 PetscFunctionReturn(0); 3213 } 3214 3215 static PetscErrorCode MatSetPreallocationCOO_IS(Mat A,PetscCount ncoo,const PetscInt coo_i[],const PetscInt coo_j[]) 3216 { 3217 Mat_IS *a = (Mat_IS*)A->data; 3218 PetscInt *coo_il, *coo_jl, incoo; 3219 3220 PetscFunctionBegin; 3221 PetscCheck(a->A,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Need to provide l2g map first via MatSetLocalToGlobalMapping"); 3222 PetscCheck(ncoo <= PETSC_MAX_INT,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ncoo %" PetscCount_FMT " overflowed PetscInt; configure --with-64-bit-indices or request support",ncoo); 3223 PetscCall(PetscMalloc2(ncoo,&coo_il,ncoo,&coo_jl)); 3224 PetscCall(ISGlobalToLocalMappingApply(a->rmapping,IS_GTOLM_MASK,ncoo,coo_i,&incoo,coo_il)); 3225 PetscCall(ISGlobalToLocalMappingApply(a->cmapping,IS_GTOLM_MASK,ncoo,coo_j,&incoo,coo_jl)); 3226 PetscCall(MatSetPreallocationCOO(a->A,ncoo,coo_il,coo_jl)); 3227 PetscCall(PetscFree2(coo_il,coo_jl)); 3228 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",MatSetValuesCOO_IS)); 3229 A->preallocated = PETSC_TRUE; 3230 PetscFunctionReturn(0); 3231 } 3232 3233 /*@ 3234 MatISGetLocalToGlobalMapping - Gets the local-to-global numbering of the MATIS object 3235 3236 Not Collective 3237 3238 Input Parameter: 3239 . A - the matrix 3240 3241 Output Parameters: 3242 + rmapping - row mapping 3243 - cmapping - column mapping 3244 3245 Notes: The returned map can be different from the one used to construct the MATIS object, since it will not contain negative or repeated indices. 3246 3247 Level: advanced 3248 3249 .seealso: `MatSetLocalToGlobalMapping()` 3250 @*/ 3251 PetscErrorCode MatISGetLocalToGlobalMapping(Mat A,ISLocalToGlobalMapping *rmapping,ISLocalToGlobalMapping *cmapping) 3252 { 3253 PetscFunctionBegin; 3254 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3255 PetscValidType(A,1); 3256 if (rmapping) PetscValidPointer(rmapping,2); 3257 if (cmapping) PetscValidPointer(cmapping,3); 3258 PetscUseMethod(A,"MatISGetLocalToGlobalMapping_C",(Mat,ISLocalToGlobalMapping*,ISLocalToGlobalMapping*),(A,rmapping,cmapping)); 3259 PetscFunctionReturn(0); 3260 } 3261 3262 static PetscErrorCode MatISGetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping *r, ISLocalToGlobalMapping *c) 3263 { 3264 Mat_IS *a = (Mat_IS*)A->data; 3265 3266 PetscFunctionBegin; 3267 if (r) *r = a->rmapping; 3268 if (c) *c = a->cmapping; 3269 PetscFunctionReturn(0); 3270 } 3271 3272 /*MC 3273 MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP). 3274 This stores the matrices in globally unassembled form. Each processor assembles only its local Neumann problem and the parallel matrix vector 3275 product is handled "implicitly". 3276 3277 Options Database Keys: 3278 + -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions() 3279 . -matis_fixempty - Fixes local matrices in case of empty local rows/columns. 3280 - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of MatPtAP(). 3281 3282 Notes: 3283 Options prefix for the inner matrix are given by -is_mat_xxx 3284 3285 You must call MatSetLocalToGlobalMapping() before using this matrix type. 3286 3287 You can do matrix preallocation on the local matrix after you obtain it with 3288 MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation() 3289 3290 Level: advanced 3291 3292 .seealso: `Mat`, `MatISGetLocalMat()`, `MatSetLocalToGlobalMapping()`, `MatISSetPreallocation()`, `MatCreateIS()`, `PCBDDC`, `KSPFETIDP` 3293 3294 M*/ 3295 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A) 3296 { 3297 Mat_IS *a; 3298 3299 PetscFunctionBegin; 3300 PetscCall(PetscNewLog(A,&a)); 3301 PetscCall(PetscStrallocpy(MATAIJ,&a->lmattype)); 3302 A->data = (void*)a; 3303 3304 /* matrix ops */ 3305 PetscCall(PetscMemzero(A->ops,sizeof(struct _MatOps))); 3306 A->ops->mult = MatMult_IS; 3307 A->ops->multadd = MatMultAdd_IS; 3308 A->ops->multtranspose = MatMultTranspose_IS; 3309 A->ops->multtransposeadd = MatMultTransposeAdd_IS; 3310 A->ops->destroy = MatDestroy_IS; 3311 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS; 3312 A->ops->setvalues = MatSetValues_IS; 3313 A->ops->setvaluesblocked = MatSetValuesBlocked_IS; 3314 A->ops->setvalueslocal = MatSetValuesLocal_IS; 3315 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS; 3316 A->ops->zerorows = MatZeroRows_IS; 3317 A->ops->zerorowscolumns = MatZeroRowsColumns_IS; 3318 A->ops->assemblybegin = MatAssemblyBegin_IS; 3319 A->ops->assemblyend = MatAssemblyEnd_IS; 3320 A->ops->view = MatView_IS; 3321 A->ops->zeroentries = MatZeroEntries_IS; 3322 A->ops->scale = MatScale_IS; 3323 A->ops->getdiagonal = MatGetDiagonal_IS; 3324 A->ops->setoption = MatSetOption_IS; 3325 A->ops->ishermitian = MatIsHermitian_IS; 3326 A->ops->issymmetric = MatIsSymmetric_IS; 3327 A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS; 3328 A->ops->duplicate = MatDuplicate_IS; 3329 A->ops->missingdiagonal = MatMissingDiagonal_IS; 3330 A->ops->copy = MatCopy_IS; 3331 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 3332 A->ops->createsubmatrix = MatCreateSubMatrix_IS; 3333 A->ops->axpy = MatAXPY_IS; 3334 A->ops->diagonalset = MatDiagonalSet_IS; 3335 A->ops->shift = MatShift_IS; 3336 A->ops->transpose = MatTranspose_IS; 3337 A->ops->getinfo = MatGetInfo_IS; 3338 A->ops->diagonalscale = MatDiagonalScale_IS; 3339 A->ops->setfromoptions = MatSetFromOptions_IS; 3340 A->ops->setup = MatSetUp_IS; 3341 A->ops->hasoperation = MatHasOperation_IS; 3342 3343 /* special MATIS functions */ 3344 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",MatISSetLocalMatType_IS)); 3345 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS)); 3346 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS)); 3347 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS)); 3348 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ)); 3349 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS)); 3350 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS)); 3351 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",MatISFixLocalEmpty_IS)); 3352 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalToGlobalMapping_C",MatISGetLocalToGlobalMapping_IS)); 3353 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ)); 3354 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ)); 3355 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ)); 3356 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ)); 3357 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ)); 3358 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ)); 3359 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ)); 3360 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOOLocal_C",MatSetPreallocationCOOLocal_IS)); 3361 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_IS)); 3362 PetscCall(PetscObjectChangeTypeName((PetscObject)A,MATIS)); 3363 PetscFunctionReturn(0); 3364 } 3365