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