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 <petsc/private/matisimpl.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 if (A->rmap->n) { 597 PetscInt dc,oc,stc,*aux; 598 599 PetscCall(MatGetLocalSize(Ad,NULL,&dc)); 600 PetscCall(MatGetLocalSize(Ao,NULL,&oc)); 601 PetscCheck(!oc || garray,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"garray not present"); 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(MatDestroy(&b->dA)); 2200 PetscCall(MatDestroy(&b->assembledA)); 2201 PetscCall(PetscFree(A->data)); 2202 PetscCall(PetscObjectChangeTypeName((PetscObject)A,NULL)); 2203 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",NULL)); 2204 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL)); 2205 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL)); 2206 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL)); 2207 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL)); 2208 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL)); 2209 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",NULL)); 2210 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalToGlobalMapping_C",NULL)); 2211 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL)); 2212 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL)); 2213 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL)); 2214 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL)); 2215 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL)); 2216 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL)); 2217 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL)); 2218 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOOLocal_C",NULL)); 2219 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",NULL)); 2220 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",NULL)); 2221 PetscFunctionReturn(0); 2222 } 2223 2224 static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y) 2225 { 2226 Mat_IS *is = (Mat_IS*)A->data; 2227 PetscScalar zero = 0.0; 2228 2229 PetscFunctionBegin; 2230 /* scatter the global vector x into the local work vector */ 2231 PetscCall(VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD)); 2232 PetscCall(VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD)); 2233 2234 /* multiply the local matrix */ 2235 PetscCall(MatMult(is->A,is->x,is->y)); 2236 2237 /* scatter product back into global memory */ 2238 PetscCall(VecSet(y,zero)); 2239 PetscCall(VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE)); 2240 PetscCall(VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE)); 2241 PetscFunctionReturn(0); 2242 } 2243 2244 static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 2245 { 2246 Vec temp_vec; 2247 2248 PetscFunctionBegin; /* v3 = v2 + A * v1.*/ 2249 if (v3 != v2) { 2250 PetscCall(MatMult(A,v1,v3)); 2251 PetscCall(VecAXPY(v3,1.0,v2)); 2252 } else { 2253 PetscCall(VecDuplicate(v2,&temp_vec)); 2254 PetscCall(MatMult(A,v1,temp_vec)); 2255 PetscCall(VecAXPY(temp_vec,1.0,v2)); 2256 PetscCall(VecCopy(temp_vec,v3)); 2257 PetscCall(VecDestroy(&temp_vec)); 2258 } 2259 PetscFunctionReturn(0); 2260 } 2261 2262 static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x) 2263 { 2264 Mat_IS *is = (Mat_IS*)A->data; 2265 2266 PetscFunctionBegin; 2267 /* scatter the global vector x into the local work vector */ 2268 PetscCall(VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD)); 2269 PetscCall(VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD)); 2270 2271 /* multiply the local matrix */ 2272 PetscCall(MatMultTranspose(is->A,is->y,is->x)); 2273 2274 /* scatter product back into global vector */ 2275 PetscCall(VecSet(x,0)); 2276 PetscCall(VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE)); 2277 PetscCall(VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE)); 2278 PetscFunctionReturn(0); 2279 } 2280 2281 static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 2282 { 2283 Vec temp_vec; 2284 2285 PetscFunctionBegin; /* v3 = v2 + A' * v1.*/ 2286 if (v3 != v2) { 2287 PetscCall(MatMultTranspose(A,v1,v3)); 2288 PetscCall(VecAXPY(v3,1.0,v2)); 2289 } else { 2290 PetscCall(VecDuplicate(v2,&temp_vec)); 2291 PetscCall(MatMultTranspose(A,v1,temp_vec)); 2292 PetscCall(VecAXPY(temp_vec,1.0,v2)); 2293 PetscCall(VecCopy(temp_vec,v3)); 2294 PetscCall(VecDestroy(&temp_vec)); 2295 } 2296 PetscFunctionReturn(0); 2297 } 2298 2299 static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer) 2300 { 2301 Mat_IS *a = (Mat_IS*)A->data; 2302 PetscViewer sviewer; 2303 PetscBool isascii,view = PETSC_TRUE; 2304 2305 PetscFunctionBegin; 2306 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); 2307 if (isascii) { 2308 PetscViewerFormat format; 2309 2310 PetscCall(PetscViewerGetFormat(viewer,&format)); 2311 if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE; 2312 } 2313 if (!view) PetscFunctionReturn(0); 2314 PetscCall(PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); 2315 PetscCall(MatView(a->A,sviewer)); 2316 PetscCall(PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); 2317 PetscCall(PetscViewerFlush(viewer)); 2318 PetscFunctionReturn(0); 2319 } 2320 2321 static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat,const PetscScalar **values) 2322 { 2323 Mat_IS *is = (Mat_IS*)mat->data; 2324 MPI_Datatype nodeType; 2325 const PetscScalar *lv; 2326 PetscInt bs; 2327 2328 PetscFunctionBegin; 2329 PetscCall(MatGetBlockSize(mat,&bs)); 2330 PetscCall(MatSetBlockSize(is->A,bs)); 2331 PetscCall(MatInvertBlockDiagonal(is->A,&lv)); 2332 if (!is->bdiag) { 2333 PetscCall(PetscMalloc1(bs*mat->rmap->n,&is->bdiag)); 2334 } 2335 PetscCallMPI(MPI_Type_contiguous(bs,MPIU_SCALAR,&nodeType)); 2336 PetscCallMPI(MPI_Type_commit(&nodeType)); 2337 PetscCall(PetscSFReduceBegin(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE)); 2338 PetscCall(PetscSFReduceEnd(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE)); 2339 PetscCallMPI(MPI_Type_free(&nodeType)); 2340 if (values) *values = is->bdiag; 2341 PetscFunctionReturn(0); 2342 } 2343 2344 static PetscErrorCode MatISSetUpScatters_Private(Mat A) 2345 { 2346 Vec cglobal,rglobal; 2347 IS from; 2348 Mat_IS *is = (Mat_IS*)A->data; 2349 PetscScalar sum; 2350 const PetscInt *garray; 2351 PetscInt nr,rbs,nc,cbs; 2352 VecType rtype; 2353 2354 PetscFunctionBegin; 2355 PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping,&nr)); 2356 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping,&rbs)); 2357 PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping,&nc)); 2358 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping,&cbs)); 2359 PetscCall(VecDestroy(&is->x)); 2360 PetscCall(VecDestroy(&is->y)); 2361 PetscCall(VecDestroy(&is->counter)); 2362 PetscCall(VecScatterDestroy(&is->rctx)); 2363 PetscCall(VecScatterDestroy(&is->cctx)); 2364 PetscCall(MatCreateVecs(is->A,&is->x,&is->y)); 2365 PetscCall(VecBindToCPU(is->y,PETSC_TRUE)); 2366 PetscCall(VecGetRootType_Private(is->y,&rtype)); 2367 PetscCall(PetscFree(A->defaultvectype)); 2368 PetscCall(PetscStrallocpy(rtype,&A->defaultvectype)); 2369 PetscCall(MatCreateVecs(A,&cglobal,&rglobal)); 2370 PetscCall(VecBindToCPU(rglobal,PETSC_TRUE)); 2371 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping,&garray)); 2372 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A),rbs,nr/rbs,garray,PETSC_USE_POINTER,&from)); 2373 PetscCall(VecScatterCreate(rglobal,from,is->y,NULL,&is->rctx)); 2374 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping,&garray)); 2375 PetscCall(ISDestroy(&from)); 2376 if (is->rmapping != is->cmapping) { 2377 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping,&garray)); 2378 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A),cbs,nc/cbs,garray,PETSC_USE_POINTER,&from)); 2379 PetscCall(VecScatterCreate(cglobal,from,is->x,NULL,&is->cctx)); 2380 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping,&garray)); 2381 PetscCall(ISDestroy(&from)); 2382 } else { 2383 PetscCall(PetscObjectReference((PetscObject)is->rctx)); 2384 is->cctx = is->rctx; 2385 } 2386 PetscCall(VecDestroy(&cglobal)); 2387 2388 /* interface counter vector (local) */ 2389 PetscCall(VecDuplicate(is->y,&is->counter)); 2390 PetscCall(VecBindToCPU(is->counter,PETSC_TRUE)); 2391 PetscCall(VecSet(is->y,1.)); 2392 PetscCall(VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE)); 2393 PetscCall(VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE)); 2394 PetscCall(VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD)); 2395 PetscCall(VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD)); 2396 PetscCall(VecBindToCPU(is->y,PETSC_FALSE)); 2397 PetscCall(VecBindToCPU(is->counter,PETSC_FALSE)); 2398 2399 /* special functions for block-diagonal matrices */ 2400 PetscCall(VecSum(rglobal,&sum)); 2401 A->ops->invertblockdiagonal = NULL; 2402 if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && is->rmapping == is->cmapping) A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS; 2403 PetscCall(VecDestroy(&rglobal)); 2404 2405 /* setup SF for general purpose shared indices based communications */ 2406 PetscCall(MatISSetUpSF_IS(A)); 2407 PetscFunctionReturn(0); 2408 } 2409 2410 static PetscErrorCode MatISFilterL2GMap(Mat A, ISLocalToGlobalMapping map, ISLocalToGlobalMapping *nmap, ISLocalToGlobalMapping *lmap) 2411 { 2412 IS is; 2413 ISLocalToGlobalMappingType l2gtype; 2414 const PetscInt *idxs; 2415 PetscHSetI ht; 2416 PetscInt *nidxs; 2417 PetscInt i,n,bs,c; 2418 PetscBool flg[] = {PETSC_FALSE,PETSC_FALSE}; 2419 2420 PetscFunctionBegin; 2421 PetscCall(ISLocalToGlobalMappingGetSize(map,&n)); 2422 PetscCall(ISLocalToGlobalMappingGetBlockSize(map,&bs)); 2423 PetscCall(ISLocalToGlobalMappingGetBlockIndices(map,&idxs)); 2424 PetscCall(PetscHSetICreate(&ht)); 2425 PetscCall(PetscMalloc1(n/bs,&nidxs)); 2426 for (i=0,c=0;i<n/bs;i++) { 2427 PetscBool missing; 2428 if (idxs[i] < 0) { flg[0] = PETSC_TRUE; continue; } 2429 PetscCall(PetscHSetIQueryAdd(ht,idxs[i],&missing)); 2430 if (!missing) flg[1] = PETSC_TRUE; 2431 else nidxs[c++] = idxs[i]; 2432 } 2433 PetscCall(PetscHSetIDestroy(&ht)); 2434 PetscCall(MPIU_Allreduce(MPI_IN_PLACE,flg,2,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A))); 2435 if (!flg[0] && !flg[1]) { /* Entries are all non negative and unique */ 2436 *nmap = NULL; 2437 *lmap = NULL; 2438 PetscCall(PetscFree(nidxs)); 2439 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map,&idxs)); 2440 PetscFunctionReturn(0); 2441 } 2442 2443 /* New l2g map without negative or repeated indices */ 2444 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)A),bs,c,nidxs,PETSC_USE_POINTER,&is)); 2445 PetscCall(ISLocalToGlobalMappingCreateIS(is,nmap)); 2446 PetscCall(ISDestroy(&is)); 2447 PetscCall(ISLocalToGlobalMappingGetType(map,&l2gtype)); 2448 PetscCall(ISLocalToGlobalMappingSetType(*nmap,l2gtype)); 2449 2450 /* New local l2g map for repeated indices */ 2451 PetscCall(ISGlobalToLocalMappingApplyBlock(*nmap,IS_GTOLM_MASK,n/bs,idxs,NULL,nidxs)); 2452 PetscCall(ISCreateBlock(PETSC_COMM_SELF,bs,n/bs,nidxs,PETSC_USE_POINTER,&is)); 2453 PetscCall(ISLocalToGlobalMappingCreateIS(is,lmap)); 2454 PetscCall(ISDestroy(&is)); 2455 2456 PetscCall(PetscFree(nidxs)); 2457 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(map,&idxs)); 2458 PetscFunctionReturn(0); 2459 } 2460 2461 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping) 2462 { 2463 Mat_IS *is = (Mat_IS*)A->data; 2464 ISLocalToGlobalMapping localrmapping = NULL, localcmapping = NULL; 2465 PetscBool cong, freem[] = { PETSC_FALSE, PETSC_FALSE }; 2466 PetscInt nr,rbs,nc,cbs; 2467 2468 PetscFunctionBegin; 2469 if (rmapping) PetscCheckSameComm(A,1,rmapping,2); 2470 if (cmapping) PetscCheckSameComm(A,1,cmapping,3); 2471 2472 PetscCall(ISLocalToGlobalMappingDestroy(&is->rmapping)); 2473 PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping)); 2474 PetscCall(PetscLayoutSetUp(A->rmap)); 2475 PetscCall(PetscLayoutSetUp(A->cmap)); 2476 PetscCall(MatHasCongruentLayouts(A,&cong)); 2477 2478 /* If NULL, local space matches global space */ 2479 if (!rmapping) { 2480 IS is; 2481 2482 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),A->rmap->N,0,1,&is)); 2483 PetscCall(ISLocalToGlobalMappingCreateIS(is,&rmapping)); 2484 if (A->rmap->bs > 0) PetscCall(ISLocalToGlobalMappingSetBlockSize(rmapping,A->rmap->bs)); 2485 PetscCall(ISDestroy(&is)); 2486 freem[0] = PETSC_TRUE; 2487 if (!cmapping && cong && A->rmap->bs == A->cmap->bs) cmapping = rmapping; 2488 } else if (!is->islocalref) { /* check if the l2g map has negative or repeated entries */ 2489 PetscCall(MatISFilterL2GMap(A,rmapping,&is->rmapping,&localrmapping)); 2490 if (rmapping == cmapping) { 2491 PetscCall(PetscObjectReference((PetscObject)is->rmapping)); 2492 is->cmapping = is->rmapping; 2493 PetscCall(PetscObjectReference((PetscObject)localrmapping)); 2494 localcmapping = localrmapping; 2495 } 2496 } 2497 if (!cmapping) { 2498 IS is; 2499 2500 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),A->cmap->N,0,1,&is)); 2501 PetscCall(ISLocalToGlobalMappingCreateIS(is,&cmapping)); 2502 if (A->cmap->bs > 0) PetscCall(ISLocalToGlobalMappingSetBlockSize(cmapping,A->cmap->bs)); 2503 PetscCall(ISDestroy(&is)); 2504 freem[1] = PETSC_TRUE; 2505 } else if (cmapping != rmapping && !is->islocalref) { /* check if the l2g map has negative or repeated entries */ 2506 PetscCall(MatISFilterL2GMap(A,cmapping,&is->cmapping,&localcmapping)); 2507 } 2508 if (!is->rmapping) { 2509 PetscCall(PetscObjectReference((PetscObject)rmapping)); 2510 is->rmapping = rmapping; 2511 } 2512 if (!is->cmapping) { 2513 PetscCall(PetscObjectReference((PetscObject)cmapping)); 2514 is->cmapping = cmapping; 2515 } 2516 2517 /* Clean up */ 2518 PetscCall(MatDestroy(&is->A)); 2519 if (is->csf != is->sf) { 2520 PetscCall(PetscSFDestroy(&is->csf)); 2521 PetscCall(PetscFree2(is->csf_rootdata,is->csf_leafdata)); 2522 } else is->csf = NULL; 2523 PetscCall(PetscSFDestroy(&is->sf)); 2524 PetscCall(PetscFree2(is->sf_rootdata,is->sf_leafdata)); 2525 PetscCall(PetscFree(is->bdiag)); 2526 2527 /* check if the two mappings are actually the same for square matrices since MATIS has some optimization for this case 2528 (DOLFIN passes 2 different objects) */ 2529 PetscCall(ISLocalToGlobalMappingGetSize(is->rmapping,&nr)); 2530 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->rmapping,&rbs)); 2531 PetscCall(ISLocalToGlobalMappingGetSize(is->cmapping,&nc)); 2532 PetscCall(ISLocalToGlobalMappingGetBlockSize(is->cmapping,&cbs)); 2533 if (is->rmapping != is->cmapping && cong) { 2534 PetscBool same = PETSC_FALSE; 2535 if (nr == nc && cbs == rbs) { 2536 const PetscInt *idxs1,*idxs2; 2537 2538 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->rmapping,&idxs1)); 2539 PetscCall(ISLocalToGlobalMappingGetBlockIndices(is->cmapping,&idxs2)); 2540 PetscCall(PetscArraycmp(idxs1,idxs2,nr/rbs,&same)); 2541 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping,&idxs1)); 2542 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping,&idxs2)); 2543 } 2544 PetscCall(MPIU_Allreduce(MPI_IN_PLACE,&same,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A))); 2545 if (same) { 2546 PetscCall(ISLocalToGlobalMappingDestroy(&is->cmapping)); 2547 PetscCall(PetscObjectReference((PetscObject)is->rmapping)); 2548 is->cmapping = is->rmapping; 2549 } 2550 } 2551 PetscCall(PetscLayoutSetBlockSize(A->rmap,rbs)); 2552 PetscCall(PetscLayoutSetBlockSize(A->cmap,cbs)); 2553 /* Pass the user defined maps to the layout */ 2554 PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping)); 2555 PetscCall(PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping)); 2556 if (freem[0]) PetscCall(ISLocalToGlobalMappingDestroy(&rmapping)); 2557 if (freem[1]) PetscCall(ISLocalToGlobalMappingDestroy(&cmapping)); 2558 2559 /* Create the local matrix A */ 2560 PetscCall(MatCreate(PETSC_COMM_SELF,&is->A)); 2561 PetscCall(MatSetType(is->A,is->lmattype)); 2562 PetscCall(MatSetSizes(is->A,nr,nc,nr,nc)); 2563 PetscCall(MatSetBlockSizes(is->A,rbs,cbs)); 2564 PetscCall(MatSetOptionsPrefix(is->A,"is_")); 2565 PetscCall(MatAppendOptionsPrefix(is->A,((PetscObject)A)->prefix)); 2566 PetscCall(PetscLayoutSetUp(is->A->rmap)); 2567 PetscCall(PetscLayoutSetUp(is->A->cmap)); 2568 PetscCall(MatSetLocalToGlobalMapping(is->A,localrmapping,localcmapping)); 2569 PetscCall(ISLocalToGlobalMappingDestroy(&localrmapping)); 2570 PetscCall(ISLocalToGlobalMappingDestroy(&localcmapping)); 2571 2572 /* setup scatters and local vectors for MatMult */ 2573 if (!is->islocalref) PetscCall(MatISSetUpScatters_Private(A)); 2574 A->preallocated = PETSC_TRUE; 2575 PetscFunctionReturn(0); 2576 } 2577 2578 static PetscErrorCode MatSetUp_IS(Mat A) 2579 { 2580 ISLocalToGlobalMapping rmap, cmap; 2581 2582 PetscFunctionBegin; 2583 PetscCall(MatGetLocalToGlobalMapping(A,&rmap,&cmap)); 2584 if (!rmap && !cmap) { 2585 PetscCall(MatSetLocalToGlobalMapping(A,NULL,NULL)); 2586 } 2587 PetscFunctionReturn(0); 2588 } 2589 2590 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2591 { 2592 Mat_IS *is = (Mat_IS*)mat->data; 2593 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2594 2595 PetscFunctionBegin; 2596 PetscCall(ISGlobalToLocalMappingApply(is->rmapping,IS_GTOLM_MASK,m,rows,&m,rows_l)); 2597 if (m != n || rows != cols || is->cmapping != is->rmapping) { 2598 PetscCall(ISGlobalToLocalMappingApply(is->cmapping,IS_GTOLM_MASK,n,cols,&n,cols_l)); 2599 PetscCall(MatSetValues(is->A,m,rows_l,n,cols_l,values,addv)); 2600 } else { 2601 PetscCall(MatSetValues(is->A,m,rows_l,m,rows_l,values,addv)); 2602 } 2603 PetscFunctionReturn(0); 2604 } 2605 2606 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2607 { 2608 Mat_IS *is = (Mat_IS*)mat->data; 2609 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2610 2611 PetscFunctionBegin; 2612 PetscCall(ISGlobalToLocalMappingApplyBlock(is->rmapping,IS_GTOLM_MASK,m,rows,&m,rows_l)); 2613 if (m != n || rows != cols || is->cmapping != is->rmapping) { 2614 PetscCall(ISGlobalToLocalMappingApply(is->cmapping,IS_GTOLM_MASK,n,cols,&n,cols_l)); 2615 PetscCall(MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv)); 2616 } else { 2617 PetscCall(MatSetValuesBlocked(is->A,m,rows_l,n,rows_l,values,addv)); 2618 } 2619 PetscFunctionReturn(0); 2620 } 2621 2622 static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2623 { 2624 Mat_IS *is = (Mat_IS*)A->data; 2625 2626 PetscFunctionBegin; 2627 if (is->A->rmap->mapping || is->A->cmap->mapping) { 2628 PetscCall(MatSetValuesLocal(is->A,m,rows,n,cols,values,addv)); 2629 } else { 2630 PetscCall(MatSetValues(is->A,m,rows,n,cols,values,addv)); 2631 } 2632 PetscFunctionReturn(0); 2633 } 2634 2635 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2636 { 2637 Mat_IS *is = (Mat_IS*)A->data; 2638 2639 PetscFunctionBegin; 2640 if (is->A->rmap->mapping || is->A->cmap->mapping) { 2641 PetscCall(MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv)); 2642 } else { 2643 PetscCall(MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv)); 2644 } 2645 PetscFunctionReturn(0); 2646 } 2647 2648 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns) 2649 { 2650 Mat_IS *is = (Mat_IS*)A->data; 2651 2652 PetscFunctionBegin; 2653 if (!n) { 2654 is->pure_neumann = PETSC_TRUE; 2655 } else { 2656 PetscInt i; 2657 is->pure_neumann = PETSC_FALSE; 2658 2659 if (columns) { 2660 PetscCall(MatZeroRowsColumns(is->A,n,rows,diag,NULL,NULL)); 2661 } else { 2662 PetscCall(MatZeroRows(is->A,n,rows,diag,NULL,NULL)); 2663 } 2664 if (diag != 0.) { 2665 const PetscScalar *array; 2666 PetscCall(VecGetArrayRead(is->counter,&array)); 2667 for (i=0; i<n; i++) { 2668 PetscCall(MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES)); 2669 } 2670 PetscCall(VecRestoreArrayRead(is->counter,&array)); 2671 } 2672 PetscCall(MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY)); 2673 PetscCall(MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY)); 2674 } 2675 PetscFunctionReturn(0); 2676 } 2677 2678 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns) 2679 { 2680 Mat_IS *matis = (Mat_IS*)A->data; 2681 PetscInt nr,nl,len,i; 2682 PetscInt *lrows; 2683 2684 PetscFunctionBegin; 2685 if (PetscUnlikelyDebug(columns || diag != 0. || (x && b))) { 2686 PetscBool cong; 2687 2688 PetscCall(PetscLayoutCompare(A->rmap,A->cmap,&cong)); 2689 cong = (PetscBool)(cong && matis->sf == matis->csf); 2690 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"); 2691 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"); 2692 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"); 2693 } 2694 /* get locally owned rows */ 2695 PetscCall(PetscLayoutMapLocal(A->rmap,n,rows,&len,&lrows,NULL)); 2696 /* fix right hand side if needed */ 2697 if (x && b) { 2698 const PetscScalar *xx; 2699 PetscScalar *bb; 2700 2701 PetscCall(VecGetArrayRead(x, &xx)); 2702 PetscCall(VecGetArray(b, &bb)); 2703 for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]]; 2704 PetscCall(VecRestoreArrayRead(x, &xx)); 2705 PetscCall(VecRestoreArray(b, &bb)); 2706 } 2707 /* get rows associated to the local matrices */ 2708 PetscCall(MatGetSize(matis->A,&nl,NULL)); 2709 PetscCall(PetscArrayzero(matis->sf_leafdata,nl)); 2710 PetscCall(PetscArrayzero(matis->sf_rootdata,A->rmap->n)); 2711 for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1; 2712 PetscCall(PetscFree(lrows)); 2713 PetscCall(PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE)); 2714 PetscCall(PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE)); 2715 PetscCall(PetscMalloc1(nl,&lrows)); 2716 for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i; 2717 PetscCall(MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns)); 2718 PetscCall(PetscFree(lrows)); 2719 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 2720 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 2721 PetscFunctionReturn(0); 2722 } 2723 2724 static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2725 { 2726 PetscFunctionBegin; 2727 PetscCall(MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE)); 2728 PetscFunctionReturn(0); 2729 } 2730 2731 static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2732 { 2733 PetscFunctionBegin; 2734 PetscCall(MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE)); 2735 PetscFunctionReturn(0); 2736 } 2737 2738 static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type) 2739 { 2740 Mat_IS *is = (Mat_IS*)A->data; 2741 2742 PetscFunctionBegin; 2743 PetscCall(MatAssemblyBegin(is->A,type)); 2744 PetscFunctionReturn(0); 2745 } 2746 2747 static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type) 2748 { 2749 Mat_IS *is = (Mat_IS*)A->data; 2750 PetscBool lnnz; 2751 2752 PetscFunctionBegin; 2753 PetscCall(MatAssemblyEnd(is->A,type)); 2754 /* fix for local empty rows/cols */ 2755 if (is->locempty && type == MAT_FINAL_ASSEMBLY) { 2756 Mat newlA; 2757 ISLocalToGlobalMapping rl2g,cl2g; 2758 IS nzr,nzc; 2759 PetscInt nr,nc,nnzr,nnzc; 2760 PetscBool lnewl2g,newl2g; 2761 2762 PetscCall(MatGetSize(is->A,&nr,&nc)); 2763 PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_FALSE,PETSC_SMALL,&nzr)); 2764 if (!nzr) { 2765 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&nzr)); 2766 } 2767 PetscCall(MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_TRUE,PETSC_SMALL,&nzc)); 2768 if (!nzc) { 2769 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&nzc)); 2770 } 2771 PetscCall(ISGetSize(nzr,&nnzr)); 2772 PetscCall(ISGetSize(nzc,&nnzc)); 2773 if (nnzr != nr || nnzc != nc) { /* need new global l2g map */ 2774 lnewl2g = PETSC_TRUE; 2775 PetscCallMPI(MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A))); 2776 2777 /* extract valid submatrix */ 2778 PetscCall(MatCreateSubMatrix(is->A,nzr,nzc,MAT_INITIAL_MATRIX,&newlA)); 2779 } else { /* local matrix fully populated */ 2780 lnewl2g = PETSC_FALSE; 2781 PetscCallMPI(MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A))); 2782 PetscCall(PetscObjectReference((PetscObject)is->A)); 2783 newlA = is->A; 2784 } 2785 2786 /* attach new global l2g map if needed */ 2787 if (newl2g) { 2788 IS zr,zc; 2789 const PetscInt *ridxs,*cidxs,*zridxs,*zcidxs; 2790 PetscInt *nidxs,i; 2791 2792 PetscCall(ISComplement(nzr,0,nr,&zr)); 2793 PetscCall(ISComplement(nzc,0,nc,&zc)); 2794 PetscCall(PetscMalloc1(PetscMax(nr,nc),&nidxs)); 2795 PetscCall(ISLocalToGlobalMappingGetIndices(is->rmapping,&ridxs)); 2796 PetscCall(ISLocalToGlobalMappingGetIndices(is->cmapping,&cidxs)); 2797 PetscCall(ISGetIndices(zr,&zridxs)); 2798 PetscCall(ISGetIndices(zc,&zcidxs)); 2799 PetscCall(ISGetLocalSize(zr,&nnzr)); 2800 PetscCall(ISGetLocalSize(zc,&nnzc)); 2801 2802 PetscCall(PetscArraycpy(nidxs,ridxs,nr)); 2803 for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1; 2804 PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nr,nidxs,PETSC_COPY_VALUES,&rl2g)); 2805 PetscCall(PetscArraycpy(nidxs,cidxs,nc)); 2806 for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1; 2807 PetscCall(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nc,nidxs,PETSC_COPY_VALUES,&cl2g)); 2808 2809 PetscCall(ISRestoreIndices(zr,&zridxs)); 2810 PetscCall(ISRestoreIndices(zc,&zcidxs)); 2811 PetscCall(ISLocalToGlobalMappingRestoreIndices(is->rmapping,&ridxs)); 2812 PetscCall(ISLocalToGlobalMappingRestoreIndices(is->cmapping,&cidxs)); 2813 PetscCall(ISDestroy(&nzr)); 2814 PetscCall(ISDestroy(&nzc)); 2815 PetscCall(ISDestroy(&zr)); 2816 PetscCall(ISDestroy(&zc)); 2817 PetscCall(PetscFree(nidxs)); 2818 PetscCall(MatSetLocalToGlobalMapping(A,rl2g,cl2g)); 2819 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 2820 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 2821 } 2822 PetscCall(MatISSetLocalMat(A,newlA)); 2823 PetscCall(MatDestroy(&newlA)); 2824 PetscCall(ISDestroy(&nzr)); 2825 PetscCall(ISDestroy(&nzc)); 2826 is->locempty = PETSC_FALSE; 2827 } 2828 lnnz = (PetscBool)(is->A->nonzerostate == is->lnnzstate); 2829 is->lnnzstate = is->A->nonzerostate; 2830 PetscCall(MPIU_Allreduce(MPI_IN_PLACE,&lnnz,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A))); 2831 if (lnnz) A->nonzerostate++; 2832 PetscFunctionReturn(0); 2833 } 2834 2835 static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local) 2836 { 2837 Mat_IS *is = (Mat_IS*)mat->data; 2838 2839 PetscFunctionBegin; 2840 *local = is->A; 2841 PetscFunctionReturn(0); 2842 } 2843 2844 static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local) 2845 { 2846 PetscFunctionBegin; 2847 *local = NULL; 2848 PetscFunctionReturn(0); 2849 } 2850 2851 /*@ 2852 MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix. 2853 2854 Input Parameter: 2855 . mat - the matrix 2856 2857 Output Parameter: 2858 . local - the local matrix 2859 2860 Level: advanced 2861 2862 Notes: 2863 This can be called if you have precomputed the nonzero structure of the 2864 matrix and want to provide it to the inner matrix object to improve the performance 2865 of the MatSetValues() operation. 2866 2867 Call MatISRestoreLocalMat() when finished with the local matrix. 2868 2869 .seealso: `MATIS` 2870 @*/ 2871 PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local) 2872 { 2873 PetscFunctionBegin; 2874 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2875 PetscValidPointer(local,2); 2876 PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local)); 2877 PetscFunctionReturn(0); 2878 } 2879 2880 /*@ 2881 MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat() 2882 2883 Input Parameter: 2884 . mat - the matrix 2885 2886 Output Parameter: 2887 . local - the local matrix 2888 2889 Level: advanced 2890 2891 .seealso: `MATIS` 2892 @*/ 2893 PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local) 2894 { 2895 PetscFunctionBegin; 2896 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2897 PetscValidPointer(local,2); 2898 PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local)); 2899 PetscFunctionReturn(0); 2900 } 2901 2902 static PetscErrorCode MatISSetLocalMatType_IS(Mat mat,MatType mtype) 2903 { 2904 Mat_IS *is = (Mat_IS*)mat->data; 2905 2906 PetscFunctionBegin; 2907 if (is->A) { 2908 PetscCall(MatSetType(is->A,mtype)); 2909 } 2910 PetscCall(PetscFree(is->lmattype)); 2911 PetscCall(PetscStrallocpy(mtype,&is->lmattype)); 2912 PetscFunctionReturn(0); 2913 } 2914 2915 /*@ 2916 MatISSetLocalMatType - Specifies the type of local matrix 2917 2918 Input Parameters: 2919 + mat - the matrix 2920 - mtype - the local matrix type 2921 2922 Output Parameter: 2923 2924 Level: advanced 2925 2926 .seealso: `MATIS`, `MatSetType()`, `MatType` 2927 @*/ 2928 PetscErrorCode MatISSetLocalMatType(Mat mat,MatType mtype) 2929 { 2930 PetscFunctionBegin; 2931 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2932 PetscUseMethod(mat,"MatISSetLocalMatType_C",(Mat,MatType),(mat,mtype)); 2933 PetscFunctionReturn(0); 2934 } 2935 2936 static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local) 2937 { 2938 Mat_IS *is = (Mat_IS*)mat->data; 2939 PetscInt nrows,ncols,orows,ocols; 2940 MatType mtype,otype; 2941 PetscBool sametype = PETSC_TRUE; 2942 2943 PetscFunctionBegin; 2944 if (is->A && !is->islocalref) { 2945 PetscCall(MatGetSize(is->A,&orows,&ocols)); 2946 PetscCall(MatGetSize(local,&nrows,&ncols)); 2947 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); 2948 PetscCall(MatGetType(local,&mtype)); 2949 PetscCall(MatGetType(is->A,&otype)); 2950 PetscCall(PetscStrcmp(mtype,otype,&sametype)); 2951 } 2952 PetscCall(PetscObjectReference((PetscObject)local)); 2953 PetscCall(MatDestroy(&is->A)); 2954 is->A = local; 2955 PetscCall(MatGetType(is->A,&mtype)); 2956 PetscCall(MatISSetLocalMatType(mat,mtype)); 2957 if (!sametype && !is->islocalref) { 2958 PetscCall(MatISSetUpScatters_Private(mat)); 2959 } 2960 PetscFunctionReturn(0); 2961 } 2962 2963 /*@ 2964 MatISSetLocalMat - Replace the local matrix stored inside a MATIS object. 2965 2966 Collective on Mat 2967 2968 Input Parameters: 2969 + mat - the matrix 2970 - local - the local matrix 2971 2972 Output Parameter: 2973 2974 Level: advanced 2975 2976 Notes: 2977 This can be called if you have precomputed the local matrix and 2978 want to provide it to the matrix object MATIS. 2979 2980 .seealso: `MATIS` 2981 @*/ 2982 PetscErrorCode MatISSetLocalMat(Mat mat,Mat local) 2983 { 2984 PetscFunctionBegin; 2985 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2986 PetscValidHeaderSpecific(local,MAT_CLASSID,2); 2987 PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local)); 2988 PetscFunctionReturn(0); 2989 } 2990 2991 static PetscErrorCode MatZeroEntries_IS(Mat A) 2992 { 2993 Mat_IS *a = (Mat_IS*)A->data; 2994 2995 PetscFunctionBegin; 2996 PetscCall(MatZeroEntries(a->A)); 2997 PetscFunctionReturn(0); 2998 } 2999 3000 static PetscErrorCode MatScale_IS(Mat A,PetscScalar a) 3001 { 3002 Mat_IS *is = (Mat_IS*)A->data; 3003 3004 PetscFunctionBegin; 3005 PetscCall(MatScale(is->A,a)); 3006 PetscFunctionReturn(0); 3007 } 3008 3009 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v) 3010 { 3011 Mat_IS *is = (Mat_IS*)A->data; 3012 3013 PetscFunctionBegin; 3014 /* get diagonal of the local matrix */ 3015 PetscCall(MatGetDiagonal(is->A,is->y)); 3016 3017 /* scatter diagonal back into global vector */ 3018 PetscCall(VecSet(v,0)); 3019 PetscCall(VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE)); 3020 PetscCall(VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE)); 3021 PetscFunctionReturn(0); 3022 } 3023 3024 static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg) 3025 { 3026 Mat_IS *a = (Mat_IS*)A->data; 3027 3028 PetscFunctionBegin; 3029 PetscCall(MatSetOption(a->A,op,flg)); 3030 PetscFunctionReturn(0); 3031 } 3032 3033 static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str) 3034 { 3035 Mat_IS *y = (Mat_IS*)Y->data; 3036 Mat_IS *x; 3037 3038 PetscFunctionBegin; 3039 if (PetscDefined(USE_DEBUG)) { 3040 PetscBool ismatis; 3041 PetscCall(PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis)); 3042 PetscCheck(ismatis,PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS"); 3043 } 3044 x = (Mat_IS*)X->data; 3045 PetscCall(MatAXPY(y->A,a,x->A,str)); 3046 PetscFunctionReturn(0); 3047 } 3048 3049 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat) 3050 { 3051 Mat lA; 3052 Mat_IS *matis = (Mat_IS*)(A->data); 3053 ISLocalToGlobalMapping rl2g,cl2g; 3054 IS is; 3055 const PetscInt *rg,*rl; 3056 PetscInt nrg; 3057 PetscInt N,M,nrl,i,*idxs; 3058 3059 PetscFunctionBegin; 3060 PetscCall(ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg)); 3061 PetscCall(ISGetLocalSize(row,&nrl)); 3062 PetscCall(ISGetIndices(row,&rl)); 3063 PetscCall(ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg)); 3064 if (PetscDefined(USE_DEBUG)) { 3065 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); 3066 } 3067 PetscCall(PetscMalloc1(nrg,&idxs)); 3068 /* map from [0,nrl) to row */ 3069 for (i=0;i<nrl;i++) idxs[i] = rl[i]; 3070 for (i=nrl;i<nrg;i++) idxs[i] = -1; 3071 PetscCall(ISRestoreIndices(row,&rl)); 3072 PetscCall(ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg)); 3073 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is)); 3074 PetscCall(ISLocalToGlobalMappingCreateIS(is,&rl2g)); 3075 PetscCall(ISDestroy(&is)); 3076 /* compute new l2g map for columns */ 3077 if (col != row || matis->rmapping != matis->cmapping || matis->A->rmap->mapping != matis->A->cmap->mapping) { 3078 const PetscInt *cg,*cl; 3079 PetscInt ncg; 3080 PetscInt ncl; 3081 3082 PetscCall(ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg)); 3083 PetscCall(ISGetLocalSize(col,&ncl)); 3084 PetscCall(ISGetIndices(col,&cl)); 3085 PetscCall(ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg)); 3086 if (PetscDefined(USE_DEBUG)) { 3087 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); 3088 } 3089 PetscCall(PetscMalloc1(ncg,&idxs)); 3090 /* map from [0,ncl) to col */ 3091 for (i=0;i<ncl;i++) idxs[i] = cl[i]; 3092 for (i=ncl;i<ncg;i++) idxs[i] = -1; 3093 PetscCall(ISRestoreIndices(col,&cl)); 3094 PetscCall(ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg)); 3095 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is)); 3096 PetscCall(ISLocalToGlobalMappingCreateIS(is,&cl2g)); 3097 PetscCall(ISDestroy(&is)); 3098 } else { 3099 PetscCall(PetscObjectReference((PetscObject)rl2g)); 3100 cl2g = rl2g; 3101 } 3102 /* create the MATIS submatrix */ 3103 PetscCall(MatGetSize(A,&M,&N)); 3104 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),submat)); 3105 PetscCall(MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N)); 3106 PetscCall(MatSetType(*submat,MATIS)); 3107 matis = (Mat_IS*)((*submat)->data); 3108 matis->islocalref = PETSC_TRUE; 3109 PetscCall(MatSetLocalToGlobalMapping(*submat,rl2g,cl2g)); 3110 PetscCall(MatISGetLocalMat(A,&lA)); 3111 PetscCall(MatISSetLocalMat(*submat,lA)); 3112 PetscCall(MatSetUp(*submat)); 3113 PetscCall(MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY)); 3114 PetscCall(MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY)); 3115 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g)); 3116 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g)); 3117 3118 /* remove unsupported ops */ 3119 PetscCall(PetscMemzero((*submat)->ops,sizeof(struct _MatOps))); 3120 (*submat)->ops->destroy = MatDestroy_IS; 3121 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS; 3122 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS; 3123 (*submat)->ops->assemblybegin = MatAssemblyBegin_IS; 3124 (*submat)->ops->assemblyend = MatAssemblyEnd_IS; 3125 PetscFunctionReturn(0); 3126 } 3127 3128 static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A) 3129 { 3130 Mat_IS *a = (Mat_IS*)A->data; 3131 char type[256]; 3132 PetscBool flg; 3133 3134 PetscFunctionBegin; 3135 PetscOptionsHeadBegin(PetscOptionsObject,"MATIS options"); 3136 PetscCall(PetscOptionsBool("-matis_keepassembled","Store an assembled version if needed","MatISKeepAssembled",a->keepassembled,&a->keepassembled,NULL)); 3137 PetscCall(PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns","MatISFixLocalEmpty",a->locempty,&a->locempty,NULL)); 3138 PetscCall(PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL)); 3139 PetscCall(PetscOptionsFList("-matis_localmat_type","Matrix type","MatISSetLocalMatType",MatList,a->lmattype,type,256,&flg)); 3140 if (flg) { 3141 PetscCall(MatISSetLocalMatType(A,type)); 3142 } 3143 if (a->A) { 3144 PetscCall(MatSetFromOptions(a->A)); 3145 } 3146 PetscOptionsHeadEnd(); 3147 PetscFunctionReturn(0); 3148 } 3149 3150 /*@ 3151 MatCreateIS - Creates a "process" unassembled matrix, assembled on each 3152 process but not across processes. 3153 3154 Input Parameters: 3155 + comm - MPI communicator that will share the matrix 3156 . bs - block size of the matrix 3157 . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products 3158 . rmap - local to global map for rows 3159 - cmap - local to global map for cols 3160 3161 Output Parameter: 3162 . A - the resulting matrix 3163 3164 Level: advanced 3165 3166 Notes: 3167 See MATIS for more details. 3168 m and n are NOT related to the size of the map; they represent the size of the local parts of the distributed vectors 3169 used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices. 3170 If rmap (cmap) is NULL, then the local row (column) spaces matches the global space. 3171 3172 .seealso: `MATIS`, `MatSetLocalToGlobalMapping()` 3173 @*/ 3174 PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A) 3175 { 3176 PetscFunctionBegin; 3177 PetscCall(MatCreate(comm,A)); 3178 PetscCall(MatSetSizes(*A,m,n,M,N)); 3179 if (bs > 0) { 3180 PetscCall(MatSetBlockSize(*A,bs)); 3181 } 3182 PetscCall(MatSetType(*A,MATIS)); 3183 PetscCall(MatSetLocalToGlobalMapping(*A,rmap,cmap)); 3184 PetscFunctionReturn(0); 3185 } 3186 3187 static PetscErrorCode MatHasOperation_IS(Mat A, MatOperation op, PetscBool *has) 3188 { 3189 Mat_IS *a = (Mat_IS*)A->data; 3190 static MatOperation tobefiltered[] = { MATOP_MULT_ADD, MATOP_MULT_TRANSPOSE_ADD, MATOP_GET_DIAGONAL_BLOCK, MATOP_INCREASE_OVERLAP }; 3191 3192 PetscFunctionBegin; 3193 *has = PETSC_FALSE; 3194 if (!((void**)A->ops)[op]) PetscFunctionReturn(0); 3195 *has = PETSC_TRUE; 3196 for (PetscInt i = 0; PETSC_STATIC_ARRAY_LENGTH(tobefiltered); i++) if (op == tobefiltered[i]) PetscFunctionReturn(0); 3197 PetscCall(MatHasOperation(a->A,op,has)); 3198 PetscFunctionReturn(0); 3199 } 3200 3201 static PetscErrorCode MatSetValuesCOO_IS(Mat A,const PetscScalar v[],InsertMode imode) 3202 { 3203 Mat_IS *a = (Mat_IS*)A->data; 3204 3205 PetscFunctionBegin; 3206 PetscCall(MatSetValuesCOO(a->A,v,imode)); 3207 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 3208 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 3209 PetscFunctionReturn(0); 3210 } 3211 3212 static PetscErrorCode MatSetPreallocationCOOLocal_IS(Mat A,PetscCount ncoo,PetscInt coo_i[],PetscInt coo_j[]) 3213 { 3214 Mat_IS *a = (Mat_IS*)A->data; 3215 3216 PetscFunctionBegin; 3217 PetscCheck(a->A,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Need to provide l2g map first via MatSetLocalToGlobalMapping"); 3218 if (a->A->rmap->mapping || a->A->cmap->mapping) { 3219 PetscCall(MatSetPreallocationCOOLocal(a->A,ncoo,coo_i,coo_j)); 3220 } else { 3221 PetscCall(MatSetPreallocationCOO(a->A,ncoo,coo_i,coo_j)); 3222 } 3223 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",MatSetValuesCOO_IS)); 3224 A->preallocated = PETSC_TRUE; 3225 PetscFunctionReturn(0); 3226 } 3227 3228 static PetscErrorCode MatSetPreallocationCOO_IS(Mat A,PetscCount ncoo,const PetscInt coo_i[],const PetscInt coo_j[]) 3229 { 3230 Mat_IS *a = (Mat_IS*)A->data; 3231 PetscInt *coo_il, *coo_jl, incoo; 3232 3233 PetscFunctionBegin; 3234 PetscCheck(a->A,PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Need to provide l2g map first via MatSetLocalToGlobalMapping"); 3235 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); 3236 PetscCall(PetscMalloc2(ncoo,&coo_il,ncoo,&coo_jl)); 3237 PetscCall(ISGlobalToLocalMappingApply(a->rmapping,IS_GTOLM_MASK,ncoo,coo_i,&incoo,coo_il)); 3238 PetscCall(ISGlobalToLocalMappingApply(a->cmapping,IS_GTOLM_MASK,ncoo,coo_j,&incoo,coo_jl)); 3239 PetscCall(MatSetPreallocationCOO(a->A,ncoo,coo_il,coo_jl)); 3240 PetscCall(PetscFree2(coo_il,coo_jl)); 3241 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",MatSetValuesCOO_IS)); 3242 A->preallocated = PETSC_TRUE; 3243 PetscFunctionReturn(0); 3244 } 3245 3246 static PetscErrorCode MatISGetAssembled_Private(Mat A, Mat *tA) 3247 { 3248 Mat_IS *a = (Mat_IS*)A->data; 3249 PetscObjectState Astate, aAstate = PETSC_MIN_INT; 3250 PetscObjectState Annzstate, aAnnzstate = PETSC_MIN_INT; 3251 3252 PetscFunctionBegin; 3253 PetscCall(PetscObjectStateGet((PetscObject)A,&Astate)); 3254 Annzstate = A->nonzerostate; 3255 if (a->assembledA) { 3256 PetscCall(PetscObjectStateGet((PetscObject)a->assembledA,&aAstate)); 3257 aAnnzstate = a->assembledA->nonzerostate; 3258 } 3259 if (aAnnzstate != Annzstate) PetscCall(MatDestroy(&a->assembledA)); 3260 if (Astate != aAstate || !a->assembledA) { 3261 MatType aAtype; 3262 PetscMPIInt size; 3263 PetscInt rbs, cbs, bs; 3264 3265 /* the assembled form is used as temporary storage for parallel operations 3266 like createsubmatrices and the like, do not waste device memory */ 3267 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A),&size)); 3268 PetscCall(ISLocalToGlobalMappingGetBlockSize(a->cmapping,&cbs)); 3269 PetscCall(ISLocalToGlobalMappingGetBlockSize(a->rmapping,&rbs)); 3270 bs = rbs == cbs ? rbs : 1; 3271 if (a->assembledA) PetscCall(MatGetType(a->assembledA,&aAtype)); 3272 else if (size > 1) aAtype = bs > 1 ? MATMPIBAIJ : MATMPIAIJ; 3273 else aAtype = bs > 1 ? MATSEQBAIJ : MATSEQAIJ; 3274 3275 PetscCall(MatConvert(A,aAtype,a->assembledA ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX,&a->assembledA)); 3276 PetscCall(PetscObjectStateSet((PetscObject)a->assembledA,Astate)); 3277 a->assembledA->nonzerostate = Annzstate; 3278 } 3279 PetscCall(PetscObjectReference((PetscObject)a->assembledA)); 3280 *tA = a->assembledA; 3281 if (!a->keepassembled) PetscCall(MatDestroy(&a->assembledA)); 3282 PetscFunctionReturn(0); 3283 } 3284 3285 static PetscErrorCode MatISRestoreAssembled_Private(Mat A, Mat *tA) 3286 { 3287 PetscFunctionBegin; 3288 PetscCall(MatDestroy(tA)); 3289 *tA = NULL; 3290 PetscFunctionReturn(0); 3291 } 3292 3293 static PetscErrorCode MatGetDiagonalBlock_IS(Mat A, Mat *dA) 3294 { 3295 Mat_IS *a = (Mat_IS*)A->data; 3296 PetscObjectState Astate, dAstate = PETSC_MIN_INT; 3297 3298 PetscFunctionBegin; 3299 PetscCall(PetscObjectStateGet((PetscObject)A,&Astate)); 3300 if (a->dA) PetscCall(PetscObjectStateGet((PetscObject)a->dA,&dAstate)); 3301 if (Astate != dAstate) { 3302 Mat tA; 3303 MatType ltype; 3304 3305 PetscCall(MatDestroy(&a->dA)); 3306 PetscCall(MatISGetAssembled_Private(A,&tA)); 3307 PetscCall(MatGetDiagonalBlock(tA,&a->dA)); 3308 PetscCall(MatPropagateSymmetryOptions(tA,a->dA)); 3309 PetscCall(MatGetType(a->A,<ype)); 3310 PetscCall(MatConvert(a->dA,ltype,MAT_INPLACE_MATRIX,&a->dA)); 3311 PetscCall(PetscObjectReference((PetscObject)a->dA)); 3312 PetscCall(MatISRestoreAssembled_Private(A,&tA)); 3313 PetscCall(PetscObjectStateSet((PetscObject)a->dA,Astate)); 3314 } 3315 *dA = a->dA; 3316 PetscFunctionReturn(0); 3317 } 3318 3319 static PetscErrorCode MatCreateSubMatrices_IS(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse reuse,Mat *submat[]) 3320 { 3321 Mat tA; 3322 3323 PetscFunctionBegin; 3324 PetscCall(MatISGetAssembled_Private(A,&tA)); 3325 PetscCall(MatCreateSubMatrices(tA,n,irow,icol,reuse,submat)); 3326 /* MatCreateSubMatrices_MPIAIJ is a mess at the moment */ 3327 #if 0 3328 { 3329 Mat_IS *a = (Mat_IS*)A->data; 3330 MatType ltype; 3331 VecType vtype; 3332 char *flg; 3333 3334 PetscCall(MatGetType(a->A,<ype)); 3335 PetscCall(MatGetVecType(a->A,&vtype)); 3336 PetscCall(PetscStrstr(vtype,"cuda",&flg)); 3337 if (!flg) PetscCall(PetscStrstr(vtype,"hip",&flg)); 3338 if (!flg) PetscCall(PetscStrstr(vtype,"kokkos",&flg)); 3339 if (flg) { 3340 for (PetscInt i = 0; i < n; i++) { 3341 Mat sA = (*submat)[i]; 3342 3343 PetscCall(MatConvert(sA,ltype,MAT_INPLACE_MATRIX,&sA)); 3344 (*submat)[i] = sA; 3345 } 3346 } 3347 } 3348 #endif 3349 PetscCall(MatISRestoreAssembled_Private(A,&tA)); 3350 PetscFunctionReturn(0); 3351 } 3352 3353 static PetscErrorCode MatIncreaseOverlap_IS(Mat A, PetscInt n, IS is[], PetscInt ov) 3354 { 3355 Mat tA; 3356 3357 PetscFunctionBegin; 3358 PetscCall(MatISGetAssembled_Private(A,&tA)); 3359 PetscCall(MatIncreaseOverlap(tA,n,is,ov)); 3360 PetscCall(MatISRestoreAssembled_Private(A,&tA)); 3361 PetscFunctionReturn(0); 3362 } 3363 3364 /*@ 3365 MatISGetLocalToGlobalMapping - Gets the local-to-global numbering of the MATIS object 3366 3367 Not Collective 3368 3369 Input Parameter: 3370 . A - the matrix 3371 3372 Output Parameters: 3373 + rmapping - row mapping 3374 - cmapping - column mapping 3375 3376 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. 3377 3378 Level: advanced 3379 3380 .seealso: `MatSetLocalToGlobalMapping()` 3381 @*/ 3382 PetscErrorCode MatISGetLocalToGlobalMapping(Mat A,ISLocalToGlobalMapping *rmapping,ISLocalToGlobalMapping *cmapping) 3383 { 3384 PetscFunctionBegin; 3385 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3386 PetscValidType(A,1); 3387 if (rmapping) PetscValidPointer(rmapping,2); 3388 if (cmapping) PetscValidPointer(cmapping,3); 3389 PetscUseMethod(A,"MatISGetLocalToGlobalMapping_C",(Mat,ISLocalToGlobalMapping*,ISLocalToGlobalMapping*),(A,rmapping,cmapping)); 3390 PetscFunctionReturn(0); 3391 } 3392 3393 static PetscErrorCode MatISGetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping *r, ISLocalToGlobalMapping *c) 3394 { 3395 Mat_IS *a = (Mat_IS*)A->data; 3396 3397 PetscFunctionBegin; 3398 if (r) *r = a->rmapping; 3399 if (c) *c = a->cmapping; 3400 PetscFunctionReturn(0); 3401 } 3402 3403 /*MC 3404 MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP). 3405 This stores the matrices in globally unassembled form. Each processor assembles only its local Neumann problem and the parallel matrix vector 3406 product is handled "implicitly". 3407 3408 Options Database Keys: 3409 + -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions() 3410 . -matis_fixempty - Fixes local matrices in case of empty local rows/columns. 3411 - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of MatPtAP(). 3412 3413 Notes: 3414 Options prefix for the inner matrix are given by -is_mat_xxx 3415 3416 You must call MatSetLocalToGlobalMapping() before using this matrix type. 3417 3418 You can do matrix preallocation on the local matrix after you obtain it with 3419 MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation() 3420 3421 Level: advanced 3422 3423 .seealso: `Mat`, `MatISGetLocalMat()`, `MatSetLocalToGlobalMapping()`, `MatISSetPreallocation()`, `MatCreateIS()`, `PCBDDC`, `KSPFETIDP` 3424 3425 M*/ 3426 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A) 3427 { 3428 Mat_IS *a; 3429 3430 PetscFunctionBegin; 3431 PetscCall(PetscNewLog(A,&a)); 3432 PetscCall(PetscStrallocpy(MATAIJ,&a->lmattype)); 3433 A->data = (void*)a; 3434 3435 /* matrix ops */ 3436 PetscCall(PetscMemzero(A->ops,sizeof(struct _MatOps))); 3437 A->ops->mult = MatMult_IS; 3438 A->ops->multadd = MatMultAdd_IS; 3439 A->ops->multtranspose = MatMultTranspose_IS; 3440 A->ops->multtransposeadd = MatMultTransposeAdd_IS; 3441 A->ops->destroy = MatDestroy_IS; 3442 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS; 3443 A->ops->setvalues = MatSetValues_IS; 3444 A->ops->setvaluesblocked = MatSetValuesBlocked_IS; 3445 A->ops->setvalueslocal = MatSetValuesLocal_IS; 3446 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS; 3447 A->ops->zerorows = MatZeroRows_IS; 3448 A->ops->zerorowscolumns = MatZeroRowsColumns_IS; 3449 A->ops->assemblybegin = MatAssemblyBegin_IS; 3450 A->ops->assemblyend = MatAssemblyEnd_IS; 3451 A->ops->view = MatView_IS; 3452 A->ops->zeroentries = MatZeroEntries_IS; 3453 A->ops->scale = MatScale_IS; 3454 A->ops->getdiagonal = MatGetDiagonal_IS; 3455 A->ops->setoption = MatSetOption_IS; 3456 A->ops->ishermitian = MatIsHermitian_IS; 3457 A->ops->issymmetric = MatIsSymmetric_IS; 3458 A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS; 3459 A->ops->duplicate = MatDuplicate_IS; 3460 A->ops->missingdiagonal = MatMissingDiagonal_IS; 3461 A->ops->copy = MatCopy_IS; 3462 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 3463 A->ops->createsubmatrix = MatCreateSubMatrix_IS; 3464 A->ops->axpy = MatAXPY_IS; 3465 A->ops->diagonalset = MatDiagonalSet_IS; 3466 A->ops->shift = MatShift_IS; 3467 A->ops->transpose = MatTranspose_IS; 3468 A->ops->getinfo = MatGetInfo_IS; 3469 A->ops->diagonalscale = MatDiagonalScale_IS; 3470 A->ops->setfromoptions = MatSetFromOptions_IS; 3471 A->ops->setup = MatSetUp_IS; 3472 A->ops->hasoperation = MatHasOperation_IS; 3473 A->ops->getdiagonalblock = MatGetDiagonalBlock_IS; 3474 A->ops->createsubmatrices = MatCreateSubMatrices_IS; 3475 A->ops->increaseoverlap = MatIncreaseOverlap_IS; 3476 3477 /* special MATIS functions */ 3478 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",MatISSetLocalMatType_IS)); 3479 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS)); 3480 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS)); 3481 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS)); 3482 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ)); 3483 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS)); 3484 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS)); 3485 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",MatISFixLocalEmpty_IS)); 3486 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalToGlobalMapping_C",MatISGetLocalToGlobalMapping_IS)); 3487 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ)); 3488 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ)); 3489 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ)); 3490 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ)); 3491 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ)); 3492 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ)); 3493 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ)); 3494 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOOLocal_C",MatSetPreallocationCOOLocal_IS)); 3495 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_IS)); 3496 PetscCall(PetscObjectChangeTypeName((PetscObject)A,MATIS)); 3497 PetscFunctionReturn(0); 3498 } 3499