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