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