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 CHKERRQ(MatDestroySubMatrices(ptap->ris1 ? 2 : 1,&ptap->lP)); 26 CHKERRQ(ISDestroy(&ptap->cis0)); 27 CHKERRQ(ISDestroy(&ptap->cis1)); 28 CHKERRQ(ISDestroy(&ptap->ris0)); 29 CHKERRQ(ISDestroy(&ptap->ris1)); 30 CHKERRQ(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 CHKERRQ(PetscObjectQuery((PetscObject)C,"_MatIS_PtAP",(PetscObject*)&c)); 46 PetscCheck(c,PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Missing PtAP information"); 47 CHKERRQ(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 CHKERRQ(MatCreateSubMatrices(P,n,ris,cis,reuse,&ptap->lP)); 55 56 CHKERRQ(MatISGetLocalMat(A,&lA)); 57 CHKERRQ(MatISGetLocalMat(C,&lC)); 58 if (ptap->ris1) { /* unsymmetric A mapping */ 59 Mat lPt; 60 61 CHKERRQ(MatTranspose(ptap->lP[1],MAT_INITIAL_MATRIX,&lPt)); 62 CHKERRQ(MatMatMatMult(lPt,lA,ptap->lP[0],reuse,ptap->fill,&lC)); 63 if (matis->storel2l) { 64 CHKERRQ(PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",(PetscObject)lPt)); 65 } 66 CHKERRQ(MatDestroy(&lPt)); 67 } else { 68 CHKERRQ(MatPtAP(lA,ptap->lP[0],reuse,ptap->fill,&lC)); 69 if (matis->storel2l) { 70 CHKERRQ(PetscObjectCompose((PetscObject)C,"_MatIS_PtAP_l2l",(PetscObject)ptap->lP[0])); 71 } 72 } 73 if (reuse == MAT_INITIAL_MATRIX) { 74 CHKERRQ(MatISSetLocalMat(C,lC)); 75 CHKERRQ(MatDestroy(&lC)); 76 } 77 CHKERRQ(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 78 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject)PT,&comm)); 96 bs = 1; 97 CHKERRQ(PetscObjectBaseTypeCompare((PetscObject)PT,MATMPIAIJ,&ismpiaij)); 98 CHKERRQ(PetscObjectBaseTypeCompare((PetscObject)PT,MATMPIBAIJ,&ismpibaij)); 99 CHKERRQ(PetscObjectBaseTypeCompare((PetscObject)PT,MATSEQAIJ,&isseqaij)); 100 CHKERRQ(PetscObjectTypeCompare((PetscObject)PT,MATSEQBAIJ,&isseqbaij)); 101 if (isseqaij || isseqbaij) { 102 Pd = PT; 103 Po = NULL; 104 garray = NULL; 105 } else if (ismpiaij) { 106 CHKERRQ(MatMPIAIJGetSeqAIJ(PT,&Pd,&Po,&garray)); 107 } else if (ismpibaij) { 108 CHKERRQ(MatMPIBAIJGetSeqBAIJ(PT,&Pd,&Po,&garray)); 109 CHKERRQ(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 CHKERRQ(MatFindNonzeroRowsOrCols_Basic(Po,PETSC_TRUE,PETSC_SMALL,&zo)); 118 } 119 CHKERRQ(MatFindNonzeroRowsOrCols_Basic(Pd,PETSC_TRUE,PETSC_SMALL,&zd)); 120 121 CHKERRQ(MatGetLocalSize(PT,NULL,&dc)); 122 CHKERRQ(MatGetOwnershipRangeColumn(PT,&stc,NULL)); 123 if (Po) CHKERRQ(MatGetLocalSize(Po,NULL,&oc)); 124 else oc = 0; 125 CHKERRQ(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 CHKERRQ(ISSetBlockSize(zd,bs)); 132 CHKERRQ(ISGetLocalSize(zd,&nz)); 133 CHKERRQ(ISGetIndices(zd,&idxs)); 134 ctd = nz/bs; 135 for (i=0; i<ctd; i++) aux[i] = (idxs[bs*i]+stc)/bs; 136 CHKERRQ(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 CHKERRQ(ISSetBlockSize(zo,bs)); 147 CHKERRQ(ISGetLocalSize(zo,&nz)); 148 CHKERRQ(ISGetIndices(zo,&idxs)); 149 cto = nz/bs; 150 for (i=0; i<cto; i++) aux[i+ctd] = garray[idxs[bs*i]/bs]; 151 CHKERRQ(ISRestoreIndices(zo,&idxs)); 152 } else { 153 cto = oc/bs; 154 for (i=0; i<cto; i++) aux[i+ctd] = garray[i]; 155 } 156 CHKERRQ(ISCreateBlock(comm,bs,ctd+cto,aux,PETSC_OWN_POINTER,cis)); 157 CHKERRQ(ISDestroy(&zd)); 158 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject)A,&comm)); 175 CHKERRQ(MatSetType(C,MATIS)); 176 CHKERRQ(MatISGetLocalMat(A,&lA)); 177 CHKERRQ(MatGetType(lA,&lmtype)); 178 CHKERRQ(MatISSetLocalMatType(C,lmtype)); 179 CHKERRQ(MatGetSize(P,NULL,&N)); 180 CHKERRQ(MatGetLocalSize(P,NULL,&dc)); 181 CHKERRQ(MatSetSizes(C,dc,dc,N,N)); 182 /* Not sure about this 183 CHKERRQ(MatGetBlockSizes(P,NULL,&ibs)); 184 CHKERRQ(MatSetBlockSize(*C,ibs)); 185 */ 186 187 CHKERRQ(PetscNew(&ptap)); 188 CHKERRQ(PetscContainerCreate(PETSC_COMM_SELF,&c)); 189 CHKERRQ(PetscContainerSetPointer(c,ptap)); 190 CHKERRQ(PetscContainerSetUserDestroy(c,MatISContainerDestroyPtAP_Private)); 191 CHKERRQ(PetscObjectCompose((PetscObject)C,"_MatIS_PtAP",(PetscObject)c)); 192 CHKERRQ(PetscContainerDestroy(&c)); 193 ptap->fill = fill; 194 195 CHKERRQ(MatISGetLocalToGlobalMapping(A,&rl2g,&cl2g)); 196 197 CHKERRQ(ISLocalToGlobalMappingGetBlockSize(cl2g,&ibs)); 198 CHKERRQ(ISLocalToGlobalMappingGetSize(cl2g,&N)); 199 CHKERRQ(ISLocalToGlobalMappingGetBlockIndices(cl2g,&garray)); 200 CHKERRQ(ISCreateBlock(comm,ibs,N/ibs,garray,PETSC_COPY_VALUES,&ptap->ris0)); 201 CHKERRQ(ISLocalToGlobalMappingRestoreBlockIndices(cl2g,&garray)); 202 203 CHKERRQ(MatCreateSubMatrix(P,ptap->ris0,NULL,MAT_INITIAL_MATRIX,&PT)); 204 CHKERRQ(MatGetNonzeroColumnsLocal_Private(PT,&ptap->cis0)); 205 CHKERRQ(ISLocalToGlobalMappingCreateIS(ptap->cis0,&Ccl2g)); 206 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingGetSize(rl2g,&N1)); 214 CHKERRQ(ISLocalToGlobalMappingGetBlockSize(rl2g,&ibs1)); 215 CHKERRQ(ISLocalToGlobalMappingGetBlockIndices(rl2g,&garray)); 216 CHKERRQ(ISCreateBlock(comm,ibs,N/ibs,garray,PETSC_COPY_VALUES,&ptap->ris1)); 217 CHKERRQ(ISLocalToGlobalMappingRestoreBlockIndices(rl2g,&garray)); 218 if (ibs1 == ibs && N1 == N) { /* check if the l2gmaps are the same */ 219 const PetscInt *i1,*i2; 220 221 CHKERRQ(ISBlockGetIndices(ptap->ris0,&i1)); 222 CHKERRQ(ISBlockGetIndices(ptap->ris1,&i2)); 223 CHKERRQ(PetscArraycmp(i1,i2,N,&lsame)); 224 } 225 CHKERRMPI(MPIU_Allreduce(&lsame,&same,1,MPIU_BOOL,MPI_LAND,comm)); 226 if (same) { 227 CHKERRQ(ISDestroy(&ptap->ris1)); 228 } else { 229 CHKERRQ(MatCreateSubMatrix(P,ptap->ris1,NULL,MAT_INITIAL_MATRIX,&PT)); 230 CHKERRQ(MatGetNonzeroColumnsLocal_Private(PT,&ptap->cis1)); 231 CHKERRQ(ISLocalToGlobalMappingCreateIS(ptap->cis1,&Crl2g)); 232 CHKERRQ(MatDestroy(&PT)); 233 } 234 } 235 /* Not sure about this 236 if (!Crl2g) { 237 CHKERRQ(MatGetBlockSize(C,&ibs)); 238 CHKERRQ(ISLocalToGlobalMappingSetBlockSize(Ccl2g,ibs)); 239 } 240 */ 241 CHKERRQ(MatSetLocalToGlobalMapping(C,Crl2g ? Crl2g : Ccl2g,Ccl2g)); 242 CHKERRQ(ISLocalToGlobalMappingDestroy(&Crl2g)); 243 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISDestroy(&lf->rf[i])); 289 } 290 for (i=0;i<lf->nc;i++) { 291 CHKERRQ(ISDestroy(&lf->cf[i])); 292 } 293 CHKERRQ(PetscFree2(lf->rf,lf->cf)); 294 CHKERRQ(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 CHKERRQ(MatGetBlockSize(A,&bs)); 309 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)A),A->rmap->n/bs,0,1,&is)); 310 if (bs > 1) { 311 IS is2; 312 PetscInt i,*aux; 313 314 CHKERRQ(ISGetLocalSize(is,&i)); 315 CHKERRQ(ISGetIndices(is,(const PetscInt**)&aux)); 316 CHKERRQ(ISCreateBlock(PetscObjectComm((PetscObject)A),bs,i,aux,PETSC_COPY_VALUES,&is2)); 317 CHKERRQ(ISRestoreIndices(is,(const PetscInt**)&aux)); 318 CHKERRQ(ISDestroy(&is)); 319 is = is2; 320 } 321 CHKERRQ(ISSetIdentity(is)); 322 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&rl2g)); 323 CHKERRQ(ISDestroy(&is)); 324 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)A),A->cmap->n/bs,0,1,&is)); 325 if (bs > 1) { 326 IS is2; 327 PetscInt i,*aux; 328 329 CHKERRQ(ISGetLocalSize(is,&i)); 330 CHKERRQ(ISGetIndices(is,(const PetscInt**)&aux)); 331 CHKERRQ(ISCreateBlock(PetscObjectComm((PetscObject)A),bs,i,aux,PETSC_COPY_VALUES,&is2)); 332 CHKERRQ(ISRestoreIndices(is,(const PetscInt**)&aux)); 333 CHKERRQ(ISDestroy(&is)); 334 is = is2; 335 } 336 CHKERRQ(ISSetIdentity(is)); 337 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&cl2g)); 338 CHKERRQ(ISDestroy(&is)); 339 CHKERRQ(MatCreateIS(PetscObjectComm((PetscObject)A),bs,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N,rl2g,cl2g,&B)); 340 CHKERRQ(ISLocalToGlobalMappingDestroy(&rl2g)); 341 CHKERRQ(ISLocalToGlobalMappingDestroy(&cl2g)); 342 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&lB)); 343 if (reuse == MAT_INITIAL_MATRIX) *newmat = B; 344 } else { 345 B = *newmat; 346 CHKERRQ(PetscObjectReference((PetscObject)A)); 347 lB = A; 348 } 349 CHKERRQ(MatISSetLocalMat(B,lB)); 350 CHKERRQ(MatDestroy(&lB)); 351 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 352 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 353 if (reuse == MAT_INPLACE_MATRIX) { 354 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingGetNodeInfo(matis->rmapping,&n,&ecount,&eneighs)); 372 PetscCheckFalse(m != n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected %" PetscInt_FMT " != %" PetscInt_FMT,m,n); 373 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingRestoreNodeInfo(matis->rmapping,&n,&ecount,&eneighs)); 392 CHKERRQ(MatSeqAIJRestoreArray(matis->A,&aa)); 393 CHKERRQ(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");CHKERRQ(ierr); 416 CHKERRQ(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();CHKERRQ(ierr); 418 if (mode == MAT_IS_DISASSEMBLE_L2G_MAT) { 419 CHKERRQ(MatGetLocalToGlobalMapping(A,l2g,NULL)); 420 PetscFunctionReturn(0); 421 } 422 CHKERRQ(PetscObjectGetComm((PetscObject)A,&comm)); 423 CHKERRQ(PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ ,&ismpiaij)); 424 CHKERRQ(PetscObjectBaseTypeCompare((PetscObject)A,MATMPIBAIJ,&ismpibaij)); 425 CHKERRQ(MatGetBlockSize(A,&bs)); 426 switch (mode) { 427 case MAT_IS_DISASSEMBLE_L2G_ND: 428 CHKERRQ(MatPartitioningCreate(comm,&part)); 429 CHKERRQ(MatPartitioningSetAdjacency(part,A)); 430 CHKERRQ(PetscObjectSetOptionsPrefix((PetscObject)part,((PetscObject)A)->prefix)); 431 CHKERRQ(MatPartitioningSetFromOptions(part)); 432 CHKERRQ(MatPartitioningApplyND(part,&ndmap)); 433 CHKERRQ(MatPartitioningDestroy(&part)); 434 CHKERRQ(ISBuildTwoSided(ndmap,NULL,&ndsub)); 435 CHKERRQ(MatMPIAIJSetUseScalableIncreaseOverlap(A,PETSC_TRUE)); 436 CHKERRQ(MatIncreaseOverlap(A,1,&ndsub,1)); 437 CHKERRQ(ISLocalToGlobalMappingCreateIS(ndsub,l2g)); 438 439 /* it may happen that a separator node is not properly shared */ 440 CHKERRQ(ISLocalToGlobalMappingGetNodeInfo(*l2g,&nl,&ncount,NULL)); 441 CHKERRQ(PetscSFCreate(comm,&sf)); 442 CHKERRQ(ISLocalToGlobalMappingGetIndices(*l2g,&garray)); 443 CHKERRQ(PetscSFSetGraphLayout(sf,A->rmap,nl,NULL,PETSC_OWN_POINTER,garray)); 444 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(*l2g,&garray)); 445 CHKERRQ(PetscCalloc1(A->rmap->n,&ndmapc)); 446 CHKERRQ(PetscSFReduceBegin(sf,MPIU_INT,ncount,ndmapc,MPI_REPLACE)); 447 CHKERRQ(PetscSFReduceEnd(sf,MPIU_INT,ncount,ndmapc,MPI_REPLACE)); 448 CHKERRQ(ISLocalToGlobalMappingRestoreNodeInfo(*l2g,NULL,&ncount,NULL)); 449 CHKERRQ(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 CHKERRMPI(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);CHKERRQ(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 CHKERRQ(PetscMalloc1(nl,&lndmapi)); 470 CHKERRQ(PetscSFBcastBegin(sf,MPIU_INT,dnz,lndmapi,MPI_REPLACE)); 471 472 /* compute adjacency of isolated separators node */ 473 CHKERRQ(PetscMalloc1(gcnt,&workis)); 474 for (i = 0, cnt = 0; i < A->rmap->n; i++) { 475 if (ndmapi[i] < 0 && ndmapc[i] < 2) { 476 CHKERRQ(ISCreateStride(comm,1,i+A->rmap->rstart,1,&workis[cnt++])); 477 } 478 } 479 for (i = cnt; i < gcnt; i++) { 480 CHKERRQ(ISCreateStride(comm,0,0,1,&workis[i])); 481 } 482 for (i = 0; i < gcnt; i++) { 483 CHKERRQ(PetscObjectSetName((PetscObject)workis[i],"ISOLATED")); 484 CHKERRQ(ISViewFromOptions(workis[i],NULL,"-view_isolated_separators")); 485 } 486 487 /* no communications since all the ISes correspond to locally owned rows */ 488 CHKERRQ(MatIncreaseOverlap(A,gcnt,workis,1)); 489 490 /* end communicate global id of separators */ 491 CHKERRQ(PetscSFBcastEnd(sf,MPIU_INT,dnz,lndmapi,MPI_REPLACE)); 492 493 /* communicate new layers : create a matrix and transpose it */ 494 CHKERRQ(PetscArrayzero(dnz,A->rmap->n)); 495 CHKERRQ(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 CHKERRQ(ISGetLocalSize(workis[j],&s)); 502 CHKERRQ(ISGetIndices(workis[j],&idxs)); 503 CHKERRQ(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 CHKERRQ(PetscObjectSetName((PetscObject)workis[i],"EXTENDED")); 511 CHKERRQ(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 CHKERRQ(PetscMalloc1(j,&vals)); 516 for (i = 0; i < j; i++) vals[i] = 1.0; 517 518 CHKERRQ(MatCreate(comm,&A2)); 519 CHKERRQ(MatSetType(A2,MATMPIAIJ)); 520 CHKERRQ(MatSetSizes(A2,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 521 CHKERRQ(MatMPIAIJSetPreallocation(A2,0,dnz,0,onz)); 522 CHKERRQ(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 CHKERRQ(ISGetIndices(workis[j],&idxs)); 529 CHKERRQ(MatSetValues(A2,1,&row,s,idxs,vals,INSERT_VALUES)); 530 CHKERRQ(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 CHKERRQ(PetscFree(vals)); 536 CHKERRQ(MatAssemblyBegin(A2,MAT_FINAL_ASSEMBLY)); 537 CHKERRQ(MatAssemblyEnd(A2,MAT_FINAL_ASSEMBLY)); 538 CHKERRQ(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 CHKERRQ(ISCreateGeneral(comm,j,lndmapi,PETSC_USE_POINTER,&is)); 544 CHKERRQ(MatMPIAIJGetLocalMatCondensed(A2,MAT_INITIAL_MATRIX,&is,NULL,&A3)); 545 CHKERRQ(ISDestroy(&is)); 546 CHKERRQ(MatDestroy(&A2)); 547 548 /* extend local to global map to include connected isolated separators */ 549 CHKERRQ(PetscObjectQuery((PetscObject)A3,"_petsc_GetLocalMatCondensed_iscol",(PetscObject*)&is)); 550 PetscCheck(is,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing column map"); 551 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&ll2g)); 552 CHKERRQ(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 CHKERRQ(ISCreateGeneral(PETSC_COMM_SELF,ii[i],jj,PETSC_COPY_VALUES,&is)); 555 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingApplyIS(ll2g,is,&is2)); 558 CHKERRQ(ISDestroy(&is)); 559 CHKERRQ(ISLocalToGlobalMappingDestroy(&ll2g)); 560 561 /* add new nodes to the local-to-global map */ 562 CHKERRQ(ISLocalToGlobalMappingDestroy(l2g)); 563 CHKERRQ(ISExpand(ndsub,is2,&is)); 564 CHKERRQ(ISDestroy(&is2)); 565 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,l2g)); 566 CHKERRQ(ISDestroy(&is)); 567 568 CHKERRQ(MatDestroy(&A3)); 569 CHKERRQ(PetscFree(lndmapi)); 570 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 571 for (i = 0; i < gcnt; i++) { 572 CHKERRQ(ISDestroy(&workis[i])); 573 } 574 CHKERRQ(PetscFree(workis)); 575 } 576 CHKERRQ(ISRestoreIndices(ndmap,&ndmapi)); 577 CHKERRQ(PetscSFDestroy(&sf)); 578 CHKERRQ(PetscFree(ndmapc)); 579 CHKERRQ(ISDestroy(&ndmap)); 580 CHKERRQ(ISDestroy(&ndsub)); 581 CHKERRQ(ISLocalToGlobalMappingSetBlockSize(*l2g,bs)); 582 CHKERRQ(ISLocalToGlobalMappingViewFromOptions(*l2g,NULL,"-matis_nd_l2g_view")); 583 break; 584 case MAT_IS_DISASSEMBLE_L2G_NATURAL: 585 if (ismpiaij) { 586 CHKERRQ(MatMPIAIJGetSeqAIJ(A,&Ad,&Ao,&garray)); 587 } else if (ismpibaij) { 588 CHKERRQ(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 CHKERRQ(MatGetLocalSize(A,NULL,&dc)); 595 CHKERRQ(MatGetLocalSize(Ao,NULL,&oc)); 596 CHKERRQ(MatGetOwnershipRangeColumn(A,&stc,NULL)); 597 CHKERRQ(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] = garray[i]; 600 CHKERRQ(ISCreateBlock(comm,bs,(dc+oc)/bs,aux,PETSC_OWN_POINTER,&is)); 601 } else { 602 CHKERRQ(ISCreateBlock(comm,1,0,NULL,PETSC_OWN_POINTER,&is)); 603 } 604 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,l2g)); 605 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject)A,&comm)); 632 CHKERRMPI(MPI_Comm_size(comm,&size)); 633 if (size == 1) { 634 CHKERRQ(MatConvert_SeqXAIJ_IS(A,type,reuse,newmat)); 635 PetscFunctionReturn(0); 636 } 637 if (reuse != MAT_REUSE_MATRIX && A->cmap->N == A->rmap->N) { 638 CHKERRQ(MatMPIXAIJComputeLocalToGlobalMapping_Private(A,&rl2g)); 639 CHKERRQ(MatCreate(comm,&B)); 640 CHKERRQ(MatSetType(B,MATIS)); 641 CHKERRQ(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 642 CHKERRQ(MatSetLocalToGlobalMapping(B,rl2g,rl2g)); 643 CHKERRQ(MatGetBlockSize(A,&bs)); 644 CHKERRQ(MatSetBlockSize(B,bs)); 645 CHKERRQ(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 CHKERRQ(MatISGetLocalToGlobalMapping(B,&rl2g,&cl2g)); 657 CHKERRQ(ISLocalToGlobalMappingGetBlockIndices(rl2g,&ridx)); 658 CHKERRQ(ISLocalToGlobalMappingGetBlockIndices(cl2g,&cidx)); 659 CHKERRQ(ISLocalToGlobalMappingGetSize(rl2g,&nr)); 660 CHKERRQ(ISLocalToGlobalMappingGetSize(cl2g,&nc)); 661 CHKERRQ(ISLocalToGlobalMappingGetBlockSize(rl2g,&rbs)); 662 CHKERRQ(ISLocalToGlobalMappingGetBlockSize(cl2g,&cbs)); 663 CHKERRQ(ISCreateBlock(comm,rbs,nr/rbs,ridx,PETSC_USE_POINTER,&rows)); 664 if (rl2g != cl2g) { 665 CHKERRQ(ISCreateBlock(comm,cbs,nc/cbs,cidx,PETSC_USE_POINTER,&cols)); 666 } else { 667 CHKERRQ(PetscObjectReference((PetscObject)rows)); 668 cols = rows; 669 } 670 CHKERRQ(MatISGetLocalMat(B,&lA)); 671 CHKERRQ(MatCreateSubMatrices(A,1,&rows,&cols,MAT_INITIAL_MATRIX,&newlA)); 672 CHKERRQ(MatConvert(newlA[0],MATSEQAIJ,MAT_INPLACE_MATRIX,&newlA[0])); 673 CHKERRQ(ISLocalToGlobalMappingRestoreBlockIndices(rl2g,&ridx)); 674 CHKERRQ(ISLocalToGlobalMappingRestoreBlockIndices(cl2g,&cidx)); 675 CHKERRQ(ISDestroy(&rows)); 676 CHKERRQ(ISDestroy(&cols)); 677 if (!lA->preallocated) { /* first time */ 678 CHKERRQ(MatDuplicate(newlA[0],MAT_COPY_VALUES,&lA)); 679 CHKERRQ(MatISSetLocalMat(B,lA)); 680 CHKERRQ(PetscObjectDereference((PetscObject)lA)); 681 } 682 CHKERRQ(MatCopy(newlA[0],lA,SAME_NONZERO_PATTERN)); 683 CHKERRQ(MatDestroySubMatrices(1,&newlA)); 684 CHKERRQ(MatISScaleDisassembling_Private(B)); 685 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 686 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 687 if (was_inplace) CHKERRQ(MatHeaderReplace(A,&B)); 688 else *newmat = B; 689 PetscFunctionReturn(0); 690 } 691 /* rectangular case, just compress out the column space */ 692 CHKERRQ(PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ ,&ismpiaij)); 693 CHKERRQ(PetscObjectBaseTypeCompare((PetscObject)A,MATMPIBAIJ,&ismpibaij)); 694 if (ismpiaij) { 695 bs = 1; 696 CHKERRQ(MatMPIAIJGetSeqAIJ(A,&Ad,&Ao,&garray)); 697 } else if (ismpibaij) { 698 CHKERRQ(MatGetBlockSize(A,&bs)); 699 CHKERRQ(MatMPIBAIJGetSeqBAIJ(A,&Ad,&Ao,&garray)); 700 CHKERRQ(MatConvert(Ad,MATSEQAIJ,MAT_INITIAL_MATRIX,&Ad)); 701 CHKERRQ(MatConvert(Ao,MATSEQAIJ,MAT_INITIAL_MATRIX,&Ao)); 702 } else SETERRQ(comm,PETSC_ERR_SUP,"Type %s",((PetscObject)A)->type_name); 703 CHKERRQ(MatSeqAIJGetArray(Ad,&dd)); 704 CHKERRQ(MatSeqAIJGetArray(Ao,&od)); 705 PetscCheck(garray,comm,PETSC_ERR_ARG_WRONGSTATE,"garray not present"); 706 707 /* access relevant information from MPIAIJ */ 708 CHKERRQ(MatGetOwnershipRange(A,&str,NULL)); 709 CHKERRQ(MatGetOwnershipRangeColumn(A,&stc,NULL)); 710 CHKERRQ(MatGetLocalSize(A,&dr,&dc)); 711 CHKERRQ(MatGetLocalSize(Ao,NULL,&oc)); 712 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISCreateStride(comm,dr/bs,str/bs,1,&is)); 722 if (bs > 1) { 723 IS is2; 724 725 CHKERRQ(ISGetLocalSize(is,&i)); 726 CHKERRQ(ISGetIndices(is,(const PetscInt**)&aux)); 727 CHKERRQ(ISCreateBlock(comm,bs,i,aux,PETSC_COPY_VALUES,&is2)); 728 CHKERRQ(ISRestoreIndices(is,(const PetscInt**)&aux)); 729 CHKERRQ(ISDestroy(&is)); 730 is = is2; 731 } 732 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&rl2g)); 733 CHKERRQ(ISDestroy(&is)); 734 if (dr) { 735 CHKERRQ(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 CHKERRQ(ISCreateBlock(comm,bs,(dc+oc)/bs,aux,PETSC_OWN_POINTER,&is)); 739 lc = dc+oc; 740 } else { 741 CHKERRQ(ISCreateBlock(comm,bs,0,NULL,PETSC_OWN_POINTER,&is)); 742 lc = 0; 743 } 744 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&cl2g)); 745 CHKERRQ(ISDestroy(&is)); 746 747 /* create MATIS object */ 748 CHKERRQ(MatCreate(comm,&B)); 749 CHKERRQ(MatSetSizes(B,dr,dc,PETSC_DECIDE,PETSC_DECIDE)); 750 CHKERRQ(MatSetType(B,MATIS)); 751 CHKERRQ(MatSetBlockSize(B,bs)); 752 CHKERRQ(MatSetLocalToGlobalMapping(B,rl2g,cl2g)); 753 CHKERRQ(ISLocalToGlobalMappingDestroy(&rl2g)); 754 CHKERRQ(ISLocalToGlobalMappingDestroy(&cl2g)); 755 756 /* merge local matrices */ 757 CHKERRQ(PetscMalloc1(nnz+dr+1,&aux)); 758 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSeqAIJRestoreArray(Ad,&dd)); 776 CHKERRQ(MatSeqAIJRestoreArray(Ao,&od)); 777 778 ii = aux; 779 jj = aux+dr+1; 780 aa = data; 781 CHKERRQ(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 CHKERRQ(PetscContainerCreate(PETSC_COMM_SELF,&c)); 790 CHKERRQ(PetscContainerSetPointer(c,ptrs[i])); 791 CHKERRQ(PetscContainerSetUserDestroy(c,PetscContainerUserDestroyDefault)); 792 CHKERRQ(PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c)); 793 CHKERRQ(PetscContainerDestroy(&c)); 794 } 795 if (ismpibaij) { /* destroy converted local matrices */ 796 CHKERRQ(MatDestroy(&Ad)); 797 CHKERRQ(MatDestroy(&Ao)); 798 } 799 800 /* finalize matrix */ 801 CHKERRQ(MatISSetLocalMat(B,lA)); 802 CHKERRQ(MatDestroy(&lA)); 803 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 804 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 805 if (reuse == MAT_INPLACE_MATRIX) { 806 CHKERRQ(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 CHKERRQ(MatNestGetSubMats(A,&nr,&nc,&nest)); 823 lreuse = PETSC_FALSE; 824 rnest = NULL; 825 if (reuse == MAT_REUSE_MATRIX) { 826 PetscBool ismatis,isnest; 827 828 CHKERRQ(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 CHKERRQ(MatISGetLocalMat(*newmat,&lA)); 831 CHKERRQ(PetscObjectTypeCompare((PetscObject)lA,MATNEST,&isnest)); 832 if (isnest) { 833 CHKERRQ(MatNestGetSubMats(lA,&i,&j,&rnest)); 834 lreuse = (PetscBool)(i == nr && j == nc); 835 if (!lreuse) rnest = NULL; 836 } 837 } 838 CHKERRQ(PetscObjectGetComm((PetscObject)A,&comm)); 839 CHKERRQ(PetscCalloc2(nr,&lr,nc,&lc)); 840 CHKERRQ(PetscCalloc6(nr,&isrow,nc,&iscol,nr,&islrow,nc,&islcol,nr*nc,&snest,nr*nc,&istrans)); 841 CHKERRQ(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 CHKERRQ(PetscObjectTypeCompare((PetscObject)nest[i][j],MATTRANSPOSEMAT,&istrans[ij])); 852 if (istrans[ij]) { 853 Mat T,lT; 854 CHKERRQ(MatTransposeGetMat(nest[i][j],&T)); 855 CHKERRQ(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 CHKERRQ(MatISGetLocalMat(T,&lT)); 858 CHKERRQ(MatCreateTranspose(lT,&snest[ij])); 859 } else { 860 CHKERRQ(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 CHKERRQ(MatISGetLocalMat(nest[i][j],&snest[ij])); 863 } 864 865 /* Check compatibility of local sizes */ 866 CHKERRQ(MatGetSize(snest[ij],&l1,&l2)); 867 CHKERRQ(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 CHKERRQ(MatTransposeGetMat(nest[i][j],&T)); 891 CHKERRQ(MatISGetLocalToGlobalMapping(T,NULL,&cl2g)); 892 } else { 893 CHKERRQ(MatISGetLocalToGlobalMapping(nest[i][j],&cl2g,NULL)); 894 } 895 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingGetIndices(cl2g,&idxs1)); 906 CHKERRQ(ISLocalToGlobalMappingGetIndices(rl2g,&idxs2)); 907 CHKERRQ(PetscArraycmp(idxs1,idxs2,n1,&same)); 908 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1)); 909 CHKERRQ(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 CHKERRQ(MatTransposeGetMat(nest[j][i],&T)); 925 CHKERRQ(MatISGetLocalToGlobalMapping(T,&cl2g,NULL)); 926 } else { 927 CHKERRQ(MatISGetLocalToGlobalMapping(nest[j][i],NULL,&cl2g)); 928 } 929 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingGetIndices(cl2g,&idxs1)); 940 CHKERRQ(ISLocalToGlobalMappingGetIndices(rl2g,&idxs2)); 941 CHKERRQ(PetscArraycmp(idxs1,idxs2,n1,&same)); 942 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1)); 943 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatTransposeGetMat(usedmat,&T)); 974 usedmat = T; 975 } 976 matis = (Mat_IS*)(usedmat->data); 977 CHKERRQ(ISGetIndices(isrow[i],&idxs)); 978 if (istrans[i*nc+j]) { 979 CHKERRQ(PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 980 CHKERRQ(PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 981 } else { 982 CHKERRQ(PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 983 CHKERRQ(PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 984 } 985 CHKERRQ(ISRestoreIndices(isrow[i],&idxs)); 986 stl += lr[i]; 987 } 988 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatTransposeGetMat(usedmat,&T)); 1009 usedmat = T; 1010 } 1011 matis = (Mat_IS*)(usedmat->data); 1012 CHKERRQ(ISGetIndices(iscol[i],&idxs)); 1013 if (istrans[j*nc+i]) { 1014 CHKERRQ(PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 1015 CHKERRQ(PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 1016 } else { 1017 CHKERRQ(PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 1018 CHKERRQ(PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl,MPI_REPLACE)); 1019 } 1020 CHKERRQ(ISRestoreIndices(iscol[i],&idxs)); 1021 stl += lc[i]; 1022 } 1023 CHKERRQ(ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&cl2g)); 1024 1025 /* Create MATIS */ 1026 CHKERRQ(MatCreate(comm,&B)); 1027 CHKERRQ(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 1028 CHKERRQ(MatGetBlockSizes(A,&rbs,&cbs)); 1029 CHKERRQ(MatSetBlockSizes(B,rbs,cbs)); 1030 CHKERRQ(MatSetType(B,MATIS)); 1031 CHKERRQ(MatISSetLocalMatType(B,MATNEST)); 1032 { /* hack : avoid setup of scatters */ 1033 Mat_IS *matis = (Mat_IS*)(B->data); 1034 matis->islocalref = PETSC_TRUE; 1035 } 1036 CHKERRQ(MatSetLocalToGlobalMapping(B,rl2g,cl2g)); 1037 CHKERRQ(ISLocalToGlobalMappingDestroy(&rl2g)); 1038 CHKERRQ(ISLocalToGlobalMappingDestroy(&cl2g)); 1039 CHKERRQ(MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA)); 1040 CHKERRQ(MatNestSetVecType(lA,VECNEST)); 1041 for (i=0;i<nr*nc;i++) { 1042 if (istrans[i]) { 1043 CHKERRQ(MatDestroy(&snest[i])); 1044 } 1045 } 1046 CHKERRQ(MatISSetLocalMat(B,lA)); 1047 CHKERRQ(MatDestroy(&lA)); 1048 { /* hack : setup of scatters done here */ 1049 Mat_IS *matis = (Mat_IS*)(B->data); 1050 1051 matis->islocalref = PETSC_FALSE; 1052 CHKERRQ(MatISSetUpScatters_Private(B)); 1053 } 1054 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1055 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1056 if (reuse == MAT_INPLACE_MATRIX) { 1057 CHKERRQ(MatHeaderReplace(A,&B)); 1058 } else { 1059 *newmat = B; 1060 } 1061 } else { 1062 if (lreuse) { 1063 CHKERRQ(MatISGetLocalMat(*newmat,&lA)); 1064 for (i=0;i<nr;i++) { 1065 for (j=0;j<nc;j++) { 1066 if (snest[i*nc+j]) { 1067 CHKERRQ(MatNestSetSubMat(lA,i,j,snest[i*nc+j])); 1068 if (istrans[i*nc+j]) { 1069 CHKERRQ(MatDestroy(&snest[i*nc+j])); 1070 } 1071 } 1072 } 1073 } 1074 } else { 1075 PetscInt stl; 1076 for (i=0,stl=0;i<nr;i++) { 1077 CHKERRQ(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 CHKERRQ(ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i])); 1082 stl += lc[i]; 1083 } 1084 CHKERRQ(MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA)); 1085 for (i=0;i<nr*nc;i++) { 1086 if (istrans[i]) { 1087 CHKERRQ(MatDestroy(&snest[i])); 1088 } 1089 } 1090 CHKERRQ(MatISSetLocalMat(*newmat,lA)); 1091 CHKERRQ(MatDestroy(&lA)); 1092 } 1093 CHKERRQ(MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY)); 1094 CHKERRQ(MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY)); 1095 } 1096 1097 /* Create local matrix in MATNEST format */ 1098 convert = PETSC_FALSE; 1099 CHKERRQ(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 CHKERRQ(MatISGetLocalMat(*newmat,&lA)); 1106 CHKERRQ(MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M)); 1107 CHKERRQ(MatISSetLocalMat(*newmat,M)); 1108 CHKERRQ(MatDestroy(&M)); 1109 1110 /* attach local fields to the matrix */ 1111 CHKERRQ(PetscNew(&lf)); 1112 CHKERRQ(PetscMalloc2(nr,&lf->rf,nc,&lf->cf)); 1113 for (i=0;i<nr;i++) { 1114 PetscInt n,st; 1115 1116 CHKERRQ(ISGetLocalSize(islrow[i],&n)); 1117 CHKERRQ(ISStrideGetInfo(islrow[i],&st,NULL)); 1118 CHKERRQ(ISCreateStride(comm,n,st,1,&lf->rf[i])); 1119 } 1120 for (i=0;i<nc;i++) { 1121 PetscInt n,st; 1122 1123 CHKERRQ(ISGetLocalSize(islcol[i],&n)); 1124 CHKERRQ(ISStrideGetInfo(islcol[i],&st,NULL)); 1125 CHKERRQ(ISCreateStride(comm,n,st,1,&lf->cf[i])); 1126 } 1127 lf->nr = nr; 1128 lf->nc = nc; 1129 CHKERRQ(PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)),&c)); 1130 CHKERRQ(PetscContainerSetPointer(c,lf)); 1131 CHKERRQ(PetscContainerSetUserDestroy(c,MatISContainerDestroyFields_Private)); 1132 CHKERRQ(PetscObjectCompose((PetscObject)(*newmat),"_convert_nest_lfields",(PetscObject)c)); 1133 CHKERRQ(PetscContainerDestroy(&c)); 1134 } 1135 1136 /* Free workspace */ 1137 for (i=0;i<nr;i++) { 1138 CHKERRQ(ISDestroy(&islrow[i])); 1139 } 1140 for (i=0;i<nc;i++) { 1141 CHKERRQ(ISDestroy(&islcol[i])); 1142 } 1143 CHKERRQ(PetscFree6(isrow,iscol,islrow,islcol,snest,istrans)); 1144 CHKERRQ(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 CHKERRQ(VecGetArrayRead(l,&Y)); 1159 CHKERRQ(VecGetArray(ll,&y)); 1160 CHKERRQ(PetscSFBcastBegin(matis->sf,MPIU_SCALAR,Y,y,MPI_REPLACE)); 1161 } else { 1162 ll = NULL; 1163 } 1164 if (r) { 1165 rr = matis->x; 1166 CHKERRQ(VecGetArrayRead(r,&X)); 1167 CHKERRQ(VecGetArray(rr,&x)); 1168 CHKERRQ(PetscSFBcastBegin(matis->csf,MPIU_SCALAR,X,x,MPI_REPLACE)); 1169 } else { 1170 rr = NULL; 1171 } 1172 if (ll) { 1173 CHKERRQ(PetscSFBcastEnd(matis->sf,MPIU_SCALAR,Y,y,MPI_REPLACE)); 1174 CHKERRQ(VecRestoreArrayRead(l,&Y)); 1175 CHKERRQ(VecRestoreArray(ll,&y)); 1176 } 1177 if (rr) { 1178 CHKERRQ(PetscSFBcastEnd(matis->csf,MPIU_SCALAR,X,x,MPI_REPLACE)); 1179 CHKERRQ(VecRestoreArrayRead(r,&X)); 1180 CHKERRQ(VecRestoreArray(rr,&x)); 1181 } 1182 CHKERRQ(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 CHKERRQ(MatGetBlockSize(A,&bs)); 1195 if (matis->A->ops->getinfo) { 1196 CHKERRQ(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 CHKERRMPI(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 CHKERRMPI(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 CHKERRQ(MatCreate(PetscObjectComm((PetscObject)A),&C)); 1251 CHKERRQ(MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N)); 1252 CHKERRQ(MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs))); 1253 CHKERRQ(MatSetType(C,MATIS)); 1254 CHKERRQ(MatGetLocalToGlobalMapping(A,&rl2g,&cl2g)); 1255 CHKERRQ(MatSetLocalToGlobalMapping(C,cl2g,rl2g)); 1256 } else C = *B; 1257 1258 /* perform local transposition */ 1259 CHKERRQ(MatISGetLocalMat(A,&lA)); 1260 CHKERRQ(MatTranspose(lA,MAT_INITIAL_MATRIX,&lC)); 1261 CHKERRQ(MatSetLocalToGlobalMapping(lC,lA->cmap->mapping,lA->rmap->mapping)); 1262 CHKERRQ(MatISSetLocalMat(C,lC)); 1263 CHKERRQ(MatDestroy(&lC)); 1264 1265 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) { 1266 *B = C; 1267 } else { 1268 CHKERRQ(MatHeaderMerge(A,&C)); 1269 } 1270 CHKERRQ(MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY)); 1271 CHKERRQ(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 CHKERRQ(VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD)); 1282 CHKERRQ(VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD)); 1283 } 1284 CHKERRQ(VecPointwiseDivide(is->y,is->y,is->counter)); 1285 CHKERRQ(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 CHKERRQ(VecSet(is->y,a)); 1295 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l)); 1306 CHKERRQ(ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l)); 1307 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l)); 1318 CHKERRQ(ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l)); 1319 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatCreateVecs(mat,<est,&rtest)); 1346 CHKERRQ(ISGetLocalSize(irow,&n)); 1347 CHKERRQ(ISGetIndices(irow,&idxs)); 1348 for (i=0;i<n;i++) { 1349 CHKERRQ(VecSetValue(rtest,idxs[i],1.0,ADD_VALUES)); 1350 } 1351 CHKERRQ(VecAssemblyBegin(rtest)); 1352 CHKERRQ(VecAssemblyEnd(rtest)); 1353 CHKERRQ(VecGetLocalSize(rtest,&n)); 1354 CHKERRQ(VecGetOwnershipRange(rtest,&m,NULL)); 1355 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(rtest,&array)); 1358 CHKERRQ(ISRestoreIndices(irow,&idxs)); 1359 CHKERRQ(ISGetLocalSize(icol,&n)); 1360 CHKERRQ(ISGetIndices(icol,&idxs)); 1361 for (i=0;i<n;i++) { 1362 CHKERRQ(VecSetValue(ltest,idxs[i],1.0,ADD_VALUES)); 1363 } 1364 CHKERRQ(VecAssemblyBegin(ltest)); 1365 CHKERRQ(VecAssemblyEnd(ltest)); 1366 CHKERRQ(VecGetLocalSize(ltest,&n)); 1367 CHKERRQ(VecGetOwnershipRange(ltest,&m,NULL)); 1368 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(ltest,&array)); 1371 CHKERRQ(ISRestoreIndices(icol,&idxs)); 1372 CHKERRQ(VecDestroy(&rtest)); 1373 CHKERRQ(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 CHKERRQ(PetscObjectGetComm((PetscObject)mat,&comm)); 1385 CHKERRQ(MatGetBlockSizes(mat,&arbs,&acbs)); 1386 CHKERRQ(ISGetBlockSize(irow,&irbs)); 1387 CHKERRQ(ISGetBlockSize(icol,&icbs)); 1388 rbs = arbs == irbs ? irbs : 1; 1389 cbs = acbs == icbs ? icbs : 1; 1390 CHKERRQ(ISGetLocalSize(irow,&m)); 1391 CHKERRQ(ISGetLocalSize(icol,&n)); 1392 CHKERRQ(MatCreate(comm,newmat)); 1393 CHKERRQ(MatSetType(*newmat,MATIS)); 1394 CHKERRQ(MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE)); 1395 CHKERRQ(MatSetBlockSizes(*newmat,rbs,cbs)); 1396 /* communicate irow to their owners in the layout */ 1397 CHKERRQ(ISGetIndices(irow,&idxs)); 1398 CHKERRQ(PetscLayoutMapLocal(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs)); 1399 CHKERRQ(ISRestoreIndices(irow,&idxs)); 1400 CHKERRQ(PetscArrayzero(matis->sf_rootdata,matis->sf->nroots)); 1401 for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1; 1402 CHKERRQ(PetscFree(lidxs)); 1403 CHKERRQ(PetscFree(lgidxs)); 1404 CHKERRQ(PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE)); 1405 CHKERRQ(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 CHKERRQ(PetscMalloc1(newloc,&newgidxs)); 1408 CHKERRQ(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 CHKERRQ(ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is)); 1415 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&rl2g)); 1416 CHKERRQ(ISLocalToGlobalMappingSetBlockSize(rl2g,rbs)); 1417 CHKERRQ(ISDestroy(&is)); 1418 /* local is to extract local submatrix */ 1419 newmatis = (Mat_IS*)(*newmat)->data; 1420 CHKERRQ(ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris)); 1421 CHKERRQ(MatHasCongruentLayouts(mat,&cong)); 1422 if (cong && irow == icol && matis->csf == matis->sf) { 1423 CHKERRQ(MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g)); 1424 CHKERRQ(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 CHKERRQ(ISGetIndices(icol,&idxs)); 1431 CHKERRQ(PetscLayoutMapLocal(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs)); 1432 CHKERRQ(ISRestoreIndices(icol,&idxs)); 1433 CHKERRQ(PetscArrayzero(matis->csf_rootdata,matis->csf->nroots)); 1434 for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1; 1435 CHKERRQ(PetscFree(lidxs)); 1436 CHKERRQ(PetscFree(lgidxs)); 1437 CHKERRQ(PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata,MPI_REPLACE)); 1438 CHKERRQ(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 CHKERRQ(PetscMalloc1(newloc,&newgidxs)); 1441 CHKERRQ(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 CHKERRQ(ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is)); 1448 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&cl2g)); 1449 CHKERRQ(ISLocalToGlobalMappingSetBlockSize(cl2g,cbs)); 1450 CHKERRQ(ISDestroy(&is)); 1451 /* local is to extract local submatrix */ 1452 CHKERRQ(ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis)); 1453 CHKERRQ(MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g)); 1454 CHKERRQ(ISLocalToGlobalMappingDestroy(&cl2g)); 1455 } 1456 CHKERRQ(ISLocalToGlobalMappingDestroy(&rl2g)); 1457 } else { 1458 CHKERRQ(MatISGetLocalMat(*newmat,&newlocmat)); 1459 } 1460 CHKERRQ(MatISGetLocalMat(mat,&locmat)); 1461 newmatis = (Mat_IS*)(*newmat)->data; 1462 CHKERRQ(MatCreateSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat)); 1463 if (scall == MAT_INITIAL_MATRIX) { 1464 CHKERRQ(MatISSetLocalMat(*newmat,newlocmat)); 1465 CHKERRQ(MatDestroy(&newlocmat)); 1466 } 1467 CHKERRQ(MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY)); 1468 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatCopy(a->A,b->A,str)); 1482 CHKERRQ(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 CHKERRQ(MatCreateVecs(A,NULL,&v)); 1495 CHKERRQ(MatGetDiagonal(A,v)); 1496 CHKERRQ(VecGetLocalSize(v,&n)); 1497 CHKERRQ(VecGetArrayRead(v,&array)); 1498 for (i=0;i<n;i++) if (array[i] == 0.) break; 1499 CHKERRQ(VecRestoreArrayRead(v,&array)); 1500 CHKERRQ(VecDestroy(&v)); 1501 if (i != n) *missing = PETSC_TRUE; 1502 if (d) { 1503 *d = -1; 1504 if (*missing) { 1505 PetscInt rstart; 1506 CHKERRQ(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 CHKERRQ(PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf)); 1522 CHKERRQ(ISLocalToGlobalMappingGetIndices(matis->rmapping,&gidxs)); 1523 CHKERRQ(ISLocalToGlobalMappingGetSize(matis->rmapping,&nleaves)); 1524 CHKERRQ(PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs)); 1525 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(matis->rmapping,&gidxs)); 1526 CHKERRQ(PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata)); 1527 if (matis->rmapping != matis->cmapping) { /* setup SF for columns */ 1528 CHKERRQ(ISLocalToGlobalMappingGetSize(matis->cmapping,&nleaves)); 1529 CHKERRQ(PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf)); 1530 CHKERRQ(ISLocalToGlobalMappingGetIndices(matis->cmapping,&gidxs)); 1531 CHKERRQ(PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs)); 1532 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(matis->cmapping,&gidxs)); 1533 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE)); 1671 CHKERRQ(MatGetSize(matis->A,NULL,&nlocalcols)); 1672 CHKERRQ(MatGetBlockSize(matis->A,&bs)); 1673 CHKERRQ(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 CHKERRQ(MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata)); 1677 #if defined(PETSC_HAVE_HYPRE) 1678 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata)); 1694 1695 /* for other matrix types */ 1696 CHKERRQ(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 CHKERRMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A),&size)); 1714 CHKERRQ(MatGetSize(A,&rows,&cols)); 1715 CHKERRQ(MatGetBlockSize(A,&bs)); 1716 CHKERRQ(MatGetSize(matis->A,&local_rows,&local_cols)); 1717 CHKERRQ(PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense)); 1718 CHKERRQ(PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij)); 1719 CHKERRQ(ISLocalToGlobalMappingGetIndices(matis->rmapping,&global_indices_r)); 1720 if (matis->rmapping != matis->cmapping) { 1721 CHKERRQ(ISLocalToGlobalMappingGetIndices(matis->cmapping,&global_indices_c)); 1722 } else global_indices_c = global_indices_r; 1723 1724 if (issbaij) CHKERRQ(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 CHKERRQ(MatGetLocalSize(A,&lrows,&lcols)); 1730 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);CHKERRQ(ierr); 1731 /* All processes need to compute entire row ownership */ 1732 CHKERRQ(PetscMalloc1(rows,&row_ownership)); 1733 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatRestoreRow(matis->A,i,&ncols,&cols,NULL)); 1810 } 1811 } 1812 if (global_indices_c != global_indices_r) { 1813 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(matis->cmapping,&global_indices_c)); 1814 } 1815 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(matis->rmapping,&global_indices_r)); 1816 CHKERRQ(PetscFree(row_ownership)); 1817 1818 /* Reduce my_dnz and my_onz */ 1819 if (maxreduce) { 1820 CHKERRQ(PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX)); 1821 CHKERRQ(PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX)); 1822 CHKERRQ(PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX)); 1823 CHKERRQ(PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX)); 1824 } else { 1825 CHKERRQ(PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM)); 1826 CHKERRQ(PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM)); 1827 CHKERRQ(PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM)); 1828 CHKERRQ(PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM)); 1829 } 1830 CHKERRQ(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 CHKERRQ(MatSetBlockSizesFromMats(B,A,A)); 1840 CHKERRQ(MatSeqAIJSetPreallocation(B,0,dnz)); 1841 CHKERRQ(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 CHKERRQ(MatSeqBAIJSetPreallocation(B,bs,0,dnz)); 1853 CHKERRQ(MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz)); 1854 CHKERRQ(MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz)); 1855 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1856 if (issbaij) CHKERRQ(MatRestoreRowUpperTriangular(matis->A)); 1857 CHKERRQ(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 CHKERRMPI(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 CHKERRQ(ISLocalToGlobalMappingGetBlockSize(matis->rmapping,&rbs)); 1879 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingGetBlockIndices(matis->rmapping,&ridxs)); 1886 CHKERRQ(ISLocalToGlobalMappingGetSize(matis->rmapping,&nw)); 1887 nw = nw/rbs; 1888 CHKERRQ(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 CHKERRQ(ISCreateBlock(PETSC_COMM_SELF,rbs,nw,ridxs,PETSC_USE_POINTER,&rows)); 1893 CHKERRQ(ISSetPermutation(rows)); 1894 CHKERRQ(ISInvertPermutation(rows,PETSC_DECIDE,&irows)); 1895 CHKERRQ(ISDestroy(&rows)); 1896 } 1897 CHKERRQ(ISLocalToGlobalMappingRestoreBlockIndices(matis->rmapping,&ridxs)); 1898 CHKERRQ(PetscFree(work)); 1899 if (irows && matis->rmapping != matis->cmapping) { 1900 CHKERRQ(ISLocalToGlobalMappingGetBlockIndices(matis->cmapping,&cidxs)); 1901 CHKERRQ(ISLocalToGlobalMappingGetSize(matis->cmapping,&nw)); 1902 nw = nw/cbs; 1903 CHKERRQ(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 CHKERRQ(ISCreateBlock(PETSC_COMM_SELF,cbs,nw,cidxs,PETSC_USE_POINTER,&cols)); 1908 CHKERRQ(ISSetPermutation(cols)); 1909 CHKERRQ(ISInvertPermutation(cols,PETSC_DECIDE,&icols)); 1910 CHKERRQ(ISDestroy(&cols)); 1911 } 1912 CHKERRQ(ISLocalToGlobalMappingRestoreBlockIndices(matis->cmapping,&cidxs)); 1913 CHKERRQ(PetscFree(work)); 1914 } else if (irows) { 1915 CHKERRQ(PetscObjectReference((PetscObject)irows)); 1916 icols = irows; 1917 } 1918 } else { 1919 CHKERRQ(PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject*)&irows)); 1920 CHKERRQ(PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject*)&icols)); 1921 if (irows) CHKERRQ(PetscObjectReference((PetscObject)irows)); 1922 if (icols) CHKERRQ(PetscObjectReference((PetscObject)icols)); 1923 } 1924 if (!irows || !icols) { 1925 CHKERRQ(ISDestroy(&icols)); 1926 CHKERRQ(ISDestroy(&irows)); 1927 goto general_assembly; 1928 } 1929 CHKERRQ(MatConvert(matis->A,mtype,MAT_INITIAL_MATRIX,&B)); 1930 if (reuse != MAT_INPLACE_MATRIX) { 1931 CHKERRQ(MatCreateSubMatrix(B,irows,icols,reuse,M)); 1932 CHKERRQ(PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject)irows)); 1933 CHKERRQ(PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject)icols)); 1934 } else { 1935 Mat C; 1936 1937 CHKERRQ(MatCreateSubMatrix(B,irows,icols,MAT_INITIAL_MATRIX,&C)); 1938 CHKERRQ(MatHeaderReplace(mat,&C)); 1939 } 1940 CHKERRQ(MatDestroy(&B)); 1941 CHKERRQ(ISDestroy(&icols)); 1942 CHKERRQ(ISDestroy(&irows)); 1943 PetscFunctionReturn(0); 1944 } 1945 general_assembly: 1946 CHKERRQ(MatGetSize(mat,&rows,&cols)); 1947 CHKERRQ(ISLocalToGlobalMappingGetBlockSize(matis->rmapping,&rbs)); 1948 CHKERRQ(ISLocalToGlobalMappingGetBlockSize(matis->cmapping,&cbs)); 1949 CHKERRQ(MatGetLocalSize(mat,&lrows,&lcols)); 1950 CHKERRQ(MatGetSize(matis->A,&local_rows,&local_cols)); 1951 CHKERRQ(PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense)); 1952 CHKERRQ(PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij)); 1953 CHKERRQ(PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij)); 1954 CHKERRQ(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 CHKERRMPI(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 CHKERRQ(MatCreate(PetscObjectComm((PetscObject)mat),&MT)); 1969 CHKERRQ(MatSetSizes(MT,lrows,lcols,rows,cols)); 1970 CHKERRQ(MatSetType(MT,mtype)); 1971 CHKERRQ(MatSetBlockSizes(MT,rbs,cbs)); 1972 CHKERRQ(MatISSetMPIXAIJPreallocation_Private(mat,MT,PETSC_FALSE)); 1973 } else { 1974 PetscInt mrbs,mcbs,mrows,mcols,mlrows,mlcols; 1975 1976 /* some checks */ 1977 MT = *M; 1978 CHKERRQ(MatGetBlockSizes(MT,&mrbs,&mcbs)); 1979 CHKERRQ(MatGetSize(MT,&mrows,&mcols)); 1980 CHKERRQ(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 CHKERRQ(MatZeroEntries(MT)); 1988 } 1989 1990 if (isseqsbaij || isseqbaij) { 1991 CHKERRQ(MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&local_mat)); 1992 isseqaij = PETSC_TRUE; 1993 } else { 1994 CHKERRQ(PetscObjectReference((PetscObject)matis->A)); 1995 local_mat = matis->A; 1996 } 1997 1998 /* Set values */ 1999 CHKERRQ(MatSetLocalToGlobalMapping(MT,matis->rmapping,matis->cmapping)); 2000 if (isseqdense) { /* special case for dense local matrices */ 2001 PetscInt i,*dummy; 2002 2003 CHKERRQ(PetscMalloc1(PetscMax(local_rows,local_cols),&dummy)); 2004 for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i; 2005 CHKERRQ(MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_FALSE)); 2006 CHKERRQ(MatDenseGetArrayRead(local_mat,&array)); 2007 CHKERRQ(MatSetValuesLocal(MT,local_rows,dummy,local_cols,dummy,array,ADD_VALUES)); 2008 CHKERRQ(MatDenseRestoreArrayRead(local_mat,&array)); 2009 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSeqAIJGetArray(local_mat,&sarray)); 2019 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],sarray+xadj[i],ADD_VALUES)); 2040 } 2041 } 2042 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatGetRow(local_mat,i,&j,&local_indices_cols,&array)); 2053 CHKERRQ(MatSetValuesLocal(MT,1,&i,j,local_indices_cols,array,ADD_VALUES)); 2054 CHKERRQ(MatRestoreRow(local_mat,i,&j,&local_indices_cols,&array)); 2055 } 2056 } 2057 CHKERRQ(MatAssemblyBegin(MT,MAT_FINAL_ASSEMBLY)); 2058 CHKERRQ(MatDestroy(&local_mat)); 2059 CHKERRQ(MatAssemblyEnd(MT,MAT_FINAL_ASSEMBLY)); 2060 if (isseqdense) { 2061 CHKERRQ(MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_TRUE)); 2062 } 2063 if (reuse == MAT_INPLACE_MATRIX) { 2064 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs)); 2111 CHKERRQ(ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs)); 2112 CHKERRQ(MatGetSize(mat,&M,&N)); 2113 CHKERRQ(MatGetLocalSize(mat,&m,&n)); 2114 CHKERRQ(MatCreate(PetscObjectComm((PetscObject)mat),&B)); 2115 CHKERRQ(MatSetSizes(B,m,n,M,N)); 2116 CHKERRQ(MatSetBlockSize(B,rbs == cbs ? rbs : 1)); 2117 CHKERRQ(MatSetType(B,MATIS)); 2118 CHKERRQ(MatISSetLocalMatType(B,matis->lmattype)); 2119 CHKERRQ(MatSetLocalToGlobalMapping(B,mat->rmap->mapping,mat->cmap->mapping)); 2120 CHKERRQ(MatDuplicate(matis->A,op,&localmat)); 2121 CHKERRQ(MatSetLocalToGlobalMapping(localmat,matis->A->rmap->mapping,matis->A->cmap->mapping)); 2122 CHKERRQ(MatISSetLocalMat(B,localmat)); 2123 CHKERRQ(MatDestroy(&localmat)); 2124 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 2125 CHKERRQ(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 CHKERRQ(MatIsHermitian(matis->A,tol,&local_sym)); 2137 CHKERRMPI(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 CHKERRQ(MatIsSymmetric(matis->A,tol,&local_sym)); 2152 CHKERRMPI(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 CHKERRQ(MatIsStructurallySymmetric(matis->A,&local_sym)); 2167 CHKERRMPI(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 CHKERRQ(PetscFree(b->bdiag)); 2177 CHKERRQ(PetscFree(b->lmattype)); 2178 CHKERRQ(MatDestroy(&b->A)); 2179 CHKERRQ(VecScatterDestroy(&b->cctx)); 2180 CHKERRQ(VecScatterDestroy(&b->rctx)); 2181 CHKERRQ(VecDestroy(&b->x)); 2182 CHKERRQ(VecDestroy(&b->y)); 2183 CHKERRQ(VecDestroy(&b->counter)); 2184 CHKERRQ(ISDestroy(&b->getsub_ris)); 2185 CHKERRQ(ISDestroy(&b->getsub_cis)); 2186 if (b->sf != b->csf) { 2187 CHKERRQ(PetscSFDestroy(&b->csf)); 2188 CHKERRQ(PetscFree2(b->csf_rootdata,b->csf_leafdata)); 2189 } else b->csf = NULL; 2190 CHKERRQ(PetscSFDestroy(&b->sf)); 2191 CHKERRQ(PetscFree2(b->sf_rootdata,b->sf_leafdata)); 2192 CHKERRQ(ISLocalToGlobalMappingDestroy(&b->rmapping)); 2193 CHKERRQ(ISLocalToGlobalMappingDestroy(&b->cmapping)); 2194 CHKERRQ(PetscFree(A->data)); 2195 CHKERRQ(PetscObjectChangeTypeName((PetscObject)A,NULL)); 2196 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",NULL)); 2197 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL)); 2198 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL)); 2199 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL)); 2200 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL)); 2201 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL)); 2202 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",NULL)); 2203 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalToGlobalMapping_C",NULL)); 2204 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL)); 2205 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL)); 2206 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL)); 2207 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL)); 2208 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL)); 2209 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL)); 2210 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL)); 2211 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOOLocal_C",NULL)); 2212 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",NULL)); 2213 CHKERRQ(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 CHKERRQ(VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD)); 2225 CHKERRQ(VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD)); 2226 2227 /* multiply the local matrix */ 2228 CHKERRQ(MatMult(is->A,is->x,is->y)); 2229 2230 /* scatter product back into global memory */ 2231 CHKERRQ(VecSet(y,zero)); 2232 CHKERRQ(VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE)); 2233 CHKERRQ(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 CHKERRQ(MatMult(A,v1,v3)); 2244 CHKERRQ(VecAXPY(v3,1.0,v2)); 2245 } else { 2246 CHKERRQ(VecDuplicate(v2,&temp_vec)); 2247 CHKERRQ(MatMult(A,v1,temp_vec)); 2248 CHKERRQ(VecAXPY(temp_vec,1.0,v2)); 2249 CHKERRQ(VecCopy(temp_vec,v3)); 2250 CHKERRQ(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 CHKERRQ(VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD)); 2262 CHKERRQ(VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD)); 2263 2264 /* multiply the local matrix */ 2265 CHKERRQ(MatMultTranspose(is->A,is->y,is->x)); 2266 2267 /* scatter product back into global vector */ 2268 CHKERRQ(VecSet(x,0)); 2269 CHKERRQ(VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE)); 2270 CHKERRQ(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 CHKERRQ(MatMultTranspose(A,v1,v3)); 2281 CHKERRQ(VecAXPY(v3,1.0,v2)); 2282 } else { 2283 CHKERRQ(VecDuplicate(v2,&temp_vec)); 2284 CHKERRQ(MatMultTranspose(A,v1,temp_vec)); 2285 CHKERRQ(VecAXPY(temp_vec,1.0,v2)); 2286 CHKERRQ(VecCopy(temp_vec,v3)); 2287 CHKERRQ(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 CHKERRQ(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); 2300 if (isascii) { 2301 PetscViewerFormat format; 2302 2303 CHKERRQ(PetscViewerGetFormat(viewer,&format)); 2304 if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE; 2305 } 2306 if (!view) PetscFunctionReturn(0); 2307 CHKERRQ(PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); 2308 CHKERRQ(MatView(a->A,sviewer)); 2309 CHKERRQ(PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); 2310 CHKERRQ(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 CHKERRQ(MatGetBlockSize(mat,&bs)); 2323 CHKERRQ(MatSetBlockSize(is->A,bs)); 2324 CHKERRQ(MatInvertBlockDiagonal(is->A,&lv)); 2325 if (!is->bdiag) { 2326 CHKERRQ(PetscMalloc1(bs*mat->rmap->n,&is->bdiag)); 2327 } 2328 CHKERRMPI(MPI_Type_contiguous(bs,MPIU_SCALAR,&nodeType)); 2329 CHKERRMPI(MPI_Type_commit(&nodeType)); 2330 CHKERRQ(PetscSFReduceBegin(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE)); 2331 CHKERRQ(PetscSFReduceEnd(is->sf,nodeType,lv,is->bdiag,MPI_REPLACE)); 2332 CHKERRMPI(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 CHKERRQ(ISLocalToGlobalMappingGetSize(is->rmapping,&nr)); 2349 CHKERRQ(ISLocalToGlobalMappingGetBlockSize(is->rmapping,&rbs)); 2350 CHKERRQ(ISLocalToGlobalMappingGetSize(is->cmapping,&nc)); 2351 CHKERRQ(ISLocalToGlobalMappingGetBlockSize(is->cmapping,&cbs)); 2352 CHKERRQ(VecDestroy(&is->x)); 2353 CHKERRQ(VecDestroy(&is->y)); 2354 CHKERRQ(VecDestroy(&is->counter)); 2355 CHKERRQ(VecScatterDestroy(&is->rctx)); 2356 CHKERRQ(VecScatterDestroy(&is->cctx)); 2357 CHKERRQ(MatCreateVecs(is->A,&is->x,&is->y)); 2358 CHKERRQ(VecBindToCPU(is->y,PETSC_TRUE)); 2359 CHKERRQ(VecGetRootType_Private(is->y,&rtype)); 2360 CHKERRQ(PetscFree(A->defaultvectype)); 2361 CHKERRQ(PetscStrallocpy(rtype,&A->defaultvectype)); 2362 CHKERRQ(MatCreateVecs(A,&cglobal,&rglobal)); 2363 CHKERRQ(VecBindToCPU(rglobal,PETSC_TRUE)); 2364 CHKERRQ(ISLocalToGlobalMappingGetBlockIndices(is->rmapping,&garray)); 2365 CHKERRQ(ISCreateBlock(PetscObjectComm((PetscObject)A),rbs,nr/rbs,garray,PETSC_USE_POINTER,&from)); 2366 CHKERRQ(VecScatterCreate(rglobal,from,is->y,NULL,&is->rctx)); 2367 CHKERRQ(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping,&garray)); 2368 CHKERRQ(ISDestroy(&from)); 2369 if (is->rmapping != is->cmapping) { 2370 CHKERRQ(ISLocalToGlobalMappingGetBlockIndices(is->cmapping,&garray)); 2371 CHKERRQ(ISCreateBlock(PetscObjectComm((PetscObject)A),cbs,nc/cbs,garray,PETSC_USE_POINTER,&from)); 2372 CHKERRQ(VecScatterCreate(cglobal,from,is->x,NULL,&is->cctx)); 2373 CHKERRQ(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping,&garray)); 2374 CHKERRQ(ISDestroy(&from)); 2375 } else { 2376 CHKERRQ(PetscObjectReference((PetscObject)is->rctx)); 2377 is->cctx = is->rctx; 2378 } 2379 CHKERRQ(VecDestroy(&cglobal)); 2380 2381 /* interface counter vector (local) */ 2382 CHKERRQ(VecDuplicate(is->y,&is->counter)); 2383 CHKERRQ(VecBindToCPU(is->counter,PETSC_TRUE)); 2384 CHKERRQ(VecSet(is->y,1.)); 2385 CHKERRQ(VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE)); 2386 CHKERRQ(VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE)); 2387 CHKERRQ(VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD)); 2388 CHKERRQ(VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD)); 2389 CHKERRQ(VecBindToCPU(is->y,PETSC_FALSE)); 2390 CHKERRQ(VecBindToCPU(is->counter,PETSC_FALSE)); 2391 2392 /* special functions for block-diagonal matrices */ 2393 CHKERRQ(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 CHKERRQ(VecDestroy(&rglobal)); 2397 2398 /* setup SF for general purpose shared indices based communications */ 2399 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingGetSize(map,&n)); 2415 CHKERRQ(ISLocalToGlobalMappingGetBlockSize(map,&bs)); 2416 CHKERRQ(ISLocalToGlobalMappingGetBlockIndices(map,&idxs)); 2417 CHKERRQ(PetscHSetICreate(&ht)); 2418 CHKERRQ(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 CHKERRQ(PetscHSetIQueryAdd(ht,idxs[i],&missing)); 2423 if (!missing) flg[1] = PETSC_TRUE; 2424 else nidxs[c++] = idxs[i]; 2425 } 2426 CHKERRQ(PetscHSetIDestroy(&ht)); 2427 CHKERRMPI(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 CHKERRQ(PetscFree(nidxs)); 2432 CHKERRQ(ISLocalToGlobalMappingRestoreBlockIndices(map,&idxs)); 2433 PetscFunctionReturn(0); 2434 } 2435 2436 /* New l2g map without negative or repeated indices */ 2437 CHKERRQ(ISCreateBlock(PetscObjectComm((PetscObject)A),bs,c,nidxs,PETSC_USE_POINTER,&is)); 2438 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,nmap)); 2439 CHKERRQ(ISDestroy(&is)); 2440 CHKERRQ(ISLocalToGlobalMappingGetType(map,&l2gtype)); 2441 CHKERRQ(ISLocalToGlobalMappingSetType(*nmap,l2gtype)); 2442 2443 /* New local l2g map for repeated indices */ 2444 CHKERRQ(ISGlobalToLocalMappingApplyBlock(*nmap,IS_GTOLM_MASK,n/bs,idxs,NULL,nidxs)); 2445 CHKERRQ(ISCreateBlock(PETSC_COMM_SELF,bs,n/bs,nidxs,PETSC_USE_POINTER,&is)); 2446 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,lmap)); 2447 CHKERRQ(ISDestroy(&is)); 2448 2449 CHKERRQ(PetscFree(nidxs)); 2450 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingDestroy(&is->rmapping)); 2466 CHKERRQ(ISLocalToGlobalMappingDestroy(&is->cmapping)); 2467 CHKERRQ(PetscLayoutSetUp(A->rmap)); 2468 CHKERRQ(PetscLayoutSetUp(A->cmap)); 2469 CHKERRQ(MatHasCongruentLayouts(A,&cong)); 2470 2471 /* If NULL, local space matches global space */ 2472 if (!rmapping) { 2473 IS is; 2474 2475 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)A),A->rmap->N,0,1,&is)); 2476 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&rmapping)); 2477 if (A->rmap->bs > 0) CHKERRQ(ISLocalToGlobalMappingSetBlockSize(rmapping,A->rmap->bs)); 2478 CHKERRQ(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 CHKERRQ(MatISFilterL2GMap(A,rmapping,&is->rmapping,&localrmapping)); 2483 if (rmapping == cmapping) { 2484 CHKERRQ(PetscObjectReference((PetscObject)is->rmapping)); 2485 is->cmapping = is->rmapping; 2486 CHKERRQ(PetscObjectReference((PetscObject)localrmapping)); 2487 localcmapping = localrmapping; 2488 } 2489 } 2490 if (!cmapping) { 2491 IS is; 2492 2493 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)A),A->cmap->N,0,1,&is)); 2494 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&cmapping)); 2495 if (A->cmap->bs > 0) CHKERRQ(ISLocalToGlobalMappingSetBlockSize(cmapping,A->cmap->bs)); 2496 CHKERRQ(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 CHKERRQ(MatISFilterL2GMap(A,cmapping,&is->cmapping,&localcmapping)); 2500 } 2501 if (!is->rmapping) { 2502 CHKERRQ(PetscObjectReference((PetscObject)rmapping)); 2503 is->rmapping = rmapping; 2504 } 2505 if (!is->cmapping) { 2506 CHKERRQ(PetscObjectReference((PetscObject)cmapping)); 2507 is->cmapping = cmapping; 2508 } 2509 2510 /* Clean up */ 2511 CHKERRQ(MatDestroy(&is->A)); 2512 if (is->csf != is->sf) { 2513 CHKERRQ(PetscSFDestroy(&is->csf)); 2514 CHKERRQ(PetscFree2(is->csf_rootdata,is->csf_leafdata)); 2515 } else is->csf = NULL; 2516 CHKERRQ(PetscSFDestroy(&is->sf)); 2517 CHKERRQ(PetscFree2(is->sf_rootdata,is->sf_leafdata)); 2518 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingGetSize(is->rmapping,&nr)); 2523 CHKERRQ(ISLocalToGlobalMappingGetBlockSize(is->rmapping,&rbs)); 2524 CHKERRQ(ISLocalToGlobalMappingGetSize(is->cmapping,&nc)); 2525 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingGetBlockIndices(is->rmapping,&idxs1)); 2532 CHKERRQ(ISLocalToGlobalMappingGetBlockIndices(is->cmapping,&idxs2)); 2533 CHKERRQ(PetscArraycmp(idxs1,idxs2,nr/rbs,&same)); 2534 CHKERRQ(ISLocalToGlobalMappingRestoreBlockIndices(is->rmapping,&idxs1)); 2535 CHKERRQ(ISLocalToGlobalMappingRestoreBlockIndices(is->cmapping,&idxs2)); 2536 } 2537 CHKERRMPI(MPIU_Allreduce(MPI_IN_PLACE,&same,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A))); 2538 if (same) { 2539 CHKERRQ(ISLocalToGlobalMappingDestroy(&is->cmapping)); 2540 CHKERRQ(PetscObjectReference((PetscObject)is->rmapping)); 2541 is->cmapping = is->rmapping; 2542 } 2543 } 2544 CHKERRQ(PetscLayoutSetBlockSize(A->rmap,rbs)); 2545 CHKERRQ(PetscLayoutSetBlockSize(A->cmap,cbs)); 2546 /* Pass the user defined maps to the layout */ 2547 CHKERRQ(PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping)); 2548 CHKERRQ(PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping)); 2549 if (freem[0]) CHKERRQ(ISLocalToGlobalMappingDestroy(&rmapping)); 2550 if (freem[1]) CHKERRQ(ISLocalToGlobalMappingDestroy(&cmapping)); 2551 2552 /* Create the local matrix A */ 2553 CHKERRQ(MatCreate(PETSC_COMM_SELF,&is->A)); 2554 CHKERRQ(MatSetType(is->A,is->lmattype)); 2555 CHKERRQ(MatSetSizes(is->A,nr,nc,nr,nc)); 2556 CHKERRQ(MatSetBlockSizes(is->A,rbs,cbs)); 2557 CHKERRQ(MatSetOptionsPrefix(is->A,"is_")); 2558 CHKERRQ(MatAppendOptionsPrefix(is->A,((PetscObject)A)->prefix)); 2559 CHKERRQ(PetscLayoutSetUp(is->A->rmap)); 2560 CHKERRQ(PetscLayoutSetUp(is->A->cmap)); 2561 CHKERRQ(MatSetLocalToGlobalMapping(is->A,localrmapping,localcmapping)); 2562 CHKERRQ(ISLocalToGlobalMappingDestroy(&localrmapping)); 2563 CHKERRQ(ISLocalToGlobalMappingDestroy(&localcmapping)); 2564 2565 /* setup scatters and local vectors for MatMult */ 2566 if (!is->islocalref) CHKERRQ(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 CHKERRQ(MatGetLocalToGlobalMapping(A,&rmap,&cmap)); 2577 if (!rmap && !cmap) { 2578 CHKERRQ(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 CHKERRQ(ISGlobalToLocalMappingApply(is->rmapping,IS_GTOLM_MASK,m,rows,&m,rows_l)); 2590 if (m != n || rows != cols || is->cmapping != is->rmapping) { 2591 CHKERRQ(ISGlobalToLocalMappingApply(is->cmapping,IS_GTOLM_MASK,n,cols,&n,cols_l)); 2592 CHKERRQ(MatSetValues(is->A,m,rows_l,n,cols_l,values,addv)); 2593 } else { 2594 CHKERRQ(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 CHKERRQ(ISGlobalToLocalMappingApplyBlock(is->rmapping,IS_GTOLM_MASK,m,rows,&m,rows_l)); 2606 if (m != n || rows != cols || is->cmapping != is->rmapping) { 2607 CHKERRQ(ISGlobalToLocalMappingApply(is->cmapping,IS_GTOLM_MASK,n,cols,&n,cols_l)); 2608 CHKERRQ(MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv)); 2609 } else { 2610 CHKERRQ(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 CHKERRQ(MatSetValuesLocal(is->A,m,rows,n,cols,values,addv)); 2622 } else { 2623 CHKERRQ(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 CHKERRQ(MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv)); 2635 } else { 2636 CHKERRQ(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 CHKERRQ(MatZeroRowsColumns(is->A,n,rows,diag,NULL,NULL)); 2654 } else { 2655 CHKERRQ(MatZeroRows(is->A,n,rows,diag,NULL,NULL)); 2656 } 2657 if (diag != 0.) { 2658 const PetscScalar *array; 2659 CHKERRQ(VecGetArrayRead(is->counter,&array)); 2660 for (i=0; i<n; i++) { 2661 CHKERRQ(MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES)); 2662 } 2663 CHKERRQ(VecRestoreArrayRead(is->counter,&array)); 2664 } 2665 CHKERRQ(MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY)); 2666 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecGetArrayRead(x, &xx)); 2695 CHKERRQ(VecGetArray(b, &bb)); 2696 for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]]; 2697 CHKERRQ(VecRestoreArrayRead(x, &xx)); 2698 CHKERRQ(VecRestoreArray(b, &bb)); 2699 } 2700 /* get rows associated to the local matrices */ 2701 CHKERRQ(MatGetSize(matis->A,&nl,NULL)); 2702 CHKERRQ(PetscArrayzero(matis->sf_leafdata,nl)); 2703 CHKERRQ(PetscArrayzero(matis->sf_rootdata,A->rmap->n)); 2704 for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1; 2705 CHKERRQ(PetscFree(lrows)); 2706 CHKERRQ(PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE)); 2707 CHKERRQ(PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata,MPI_REPLACE)); 2708 CHKERRQ(PetscMalloc1(nl,&lrows)); 2709 for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i; 2710 CHKERRQ(MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns)); 2711 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatGetSize(is->A,&nr,&nc)); 2753 CHKERRQ(MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_FALSE,PETSC_SMALL,&nzr)); 2754 if (!nzr) { 2755 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&nzr)); 2756 } 2757 CHKERRQ(MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_TRUE,PETSC_SMALL,&nzc)); 2758 if (!nzc) { 2759 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&nzc)); 2760 } 2761 CHKERRQ(ISGetSize(nzr,&nnzr)); 2762 CHKERRQ(ISGetSize(nzc,&nnzc)); 2763 if (nnzr != nr || nnzc != nc) { /* need new global l2g map */ 2764 lnewl2g = PETSC_TRUE; 2765 CHKERRMPI(MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A))); 2766 2767 /* extract valid submatrix */ 2768 CHKERRQ(MatCreateSubMatrix(is->A,nzr,nzc,MAT_INITIAL_MATRIX,&newlA)); 2769 } else { /* local matrix fully populated */ 2770 lnewl2g = PETSC_FALSE; 2771 CHKERRMPI(MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A))); 2772 CHKERRQ(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 CHKERRQ(ISComplement(nzr,0,nr,&zr)); 2783 CHKERRQ(ISComplement(nzc,0,nc,&zc)); 2784 CHKERRQ(PetscMalloc1(PetscMax(nr,nc),&nidxs)); 2785 CHKERRQ(ISLocalToGlobalMappingGetIndices(is->rmapping,&ridxs)); 2786 CHKERRQ(ISLocalToGlobalMappingGetIndices(is->cmapping,&cidxs)); 2787 CHKERRQ(ISGetIndices(zr,&zridxs)); 2788 CHKERRQ(ISGetIndices(zc,&zcidxs)); 2789 CHKERRQ(ISGetLocalSize(zr,&nnzr)); 2790 CHKERRQ(ISGetLocalSize(zc,&nnzc)); 2791 2792 CHKERRQ(PetscArraycpy(nidxs,ridxs,nr)); 2793 for (i = 0; i < nnzr; i++) nidxs[zridxs[i]] = -1; 2794 CHKERRQ(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nr,nidxs,PETSC_COPY_VALUES,&rl2g)); 2795 CHKERRQ(PetscArraycpy(nidxs,cidxs,nc)); 2796 for (i = 0; i < nnzc; i++) nidxs[zcidxs[i]] = -1; 2797 CHKERRQ(ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nc,nidxs,PETSC_COPY_VALUES,&cl2g)); 2798 2799 CHKERRQ(ISRestoreIndices(zr,&zridxs)); 2800 CHKERRQ(ISRestoreIndices(zc,&zcidxs)); 2801 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(is->rmapping,&ridxs)); 2802 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(is->cmapping,&cidxs)); 2803 CHKERRQ(ISDestroy(&nzr)); 2804 CHKERRQ(ISDestroy(&nzc)); 2805 CHKERRQ(ISDestroy(&zr)); 2806 CHKERRQ(ISDestroy(&zc)); 2807 CHKERRQ(PetscFree(nidxs)); 2808 CHKERRQ(MatSetLocalToGlobalMapping(A,rl2g,cl2g)); 2809 CHKERRQ(ISLocalToGlobalMappingDestroy(&rl2g)); 2810 CHKERRQ(ISLocalToGlobalMappingDestroy(&cl2g)); 2811 } 2812 CHKERRQ(MatISSetLocalMat(A,newlA)); 2813 CHKERRQ(MatDestroy(&newlA)); 2814 CHKERRQ(ISDestroy(&nzr)); 2815 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSetType(is->A,mtype)); 2895 } 2896 CHKERRQ(PetscFree(is->lmattype)); 2897 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatGetSize(is->A,&orows,&ocols)); 2932 CHKERRQ(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 CHKERRQ(MatGetType(local,&mtype)); 2935 CHKERRQ(MatGetType(is->A,&otype)); 2936 CHKERRQ(PetscStrcmp(mtype,otype,&sametype)); 2937 } 2938 CHKERRQ(PetscObjectReference((PetscObject)local)); 2939 CHKERRQ(MatDestroy(&is->A)); 2940 is->A = local; 2941 CHKERRQ(MatGetType(is->A,&mtype)); 2942 CHKERRQ(MatISSetLocalMatType(mat,mtype)); 2943 if (!sametype && !is->islocalref) { 2944 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatGetDiagonal(is->A,is->y)); 3002 3003 /* scatter diagonal back into global vector */ 3004 CHKERRQ(VecSet(v,0)); 3005 CHKERRQ(VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE)); 3006 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg)); 3047 CHKERRQ(ISGetLocalSize(row,&nrl)); 3048 CHKERRQ(ISGetIndices(row,&rl)); 3049 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISRestoreIndices(row,&rl)); 3058 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg)); 3059 CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is)); 3060 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&rl2g)); 3061 CHKERRQ(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 CHKERRQ(ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg)); 3069 CHKERRQ(ISGetLocalSize(col,&ncl)); 3070 CHKERRQ(ISGetIndices(col,&cl)); 3071 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISRestoreIndices(col,&cl)); 3080 CHKERRQ(ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg)); 3081 CHKERRQ(ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is)); 3082 CHKERRQ(ISLocalToGlobalMappingCreateIS(is,&cl2g)); 3083 CHKERRQ(ISDestroy(&is)); 3084 } else { 3085 CHKERRQ(PetscObjectReference((PetscObject)rl2g)); 3086 cl2g = rl2g; 3087 } 3088 /* create the MATIS submatrix */ 3089 CHKERRQ(MatGetSize(A,&M,&N)); 3090 CHKERRQ(MatCreate(PetscObjectComm((PetscObject)A),submat)); 3091 CHKERRQ(MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N)); 3092 CHKERRQ(MatSetType(*submat,MATIS)); 3093 matis = (Mat_IS*)((*submat)->data); 3094 matis->islocalref = PETSC_TRUE; 3095 CHKERRQ(MatSetLocalToGlobalMapping(*submat,rl2g,cl2g)); 3096 CHKERRQ(MatISGetLocalMat(A,&lA)); 3097 CHKERRQ(MatISSetLocalMat(*submat,lA)); 3098 CHKERRQ(MatSetUp(*submat)); 3099 CHKERRQ(MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY)); 3100 CHKERRQ(MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY)); 3101 CHKERRQ(ISLocalToGlobalMappingDestroy(&rl2g)); 3102 CHKERRQ(ISLocalToGlobalMappingDestroy(&cl2g)); 3103 3104 /* remove unsupported ops */ 3105 CHKERRQ(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 CHKERRQ(PetscOptionsHead(PetscOptionsObject,"MATIS options")); 3122 CHKERRQ(PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns","MatISFixLocalEmpty",a->locempty,&a->locempty,NULL)); 3123 CHKERRQ(PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL)); 3124 CHKERRQ(PetscOptionsFList("-matis_localmat_type","Matrix type","MatISSetLocalMatType",MatList,a->lmattype,type,256,&flg)); 3125 if (flg) { 3126 CHKERRQ(MatISSetLocalMatType(A,type)); 3127 } 3128 if (a->A) { 3129 CHKERRQ(MatSetFromOptions(a->A)); 3130 } 3131 CHKERRQ(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 CHKERRQ(MatCreate(comm,A)); 3163 CHKERRQ(MatSetSizes(*A,m,n,M,N)); 3164 if (bs > 0) { 3165 CHKERRQ(MatSetBlockSize(*A,bs)); 3166 } 3167 CHKERRQ(MatSetType(*A,MATIS)); 3168 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSetValuesCOO(a->A,v,imode)); 3189 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 3190 CHKERRQ(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 CHKERRQ(MatSetPreallocationCOOLocal(a->A,ncoo,coo_i,coo_j)); 3202 } else { 3203 CHKERRQ(MatSetPreallocationCOO(a->A,ncoo,coo_i,coo_j)); 3204 } 3205 CHKERRQ(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 CHKERRQ(PetscMalloc2(ncoo,&coo_il,ncoo,&coo_jl)); 3219 CHKERRQ(ISGlobalToLocalMappingApply(a->rmapping,IS_GTOLM_MASK,ncoo,coo_i,&incoo,coo_il)); 3220 CHKERRQ(ISGlobalToLocalMappingApply(a->cmapping,IS_GTOLM_MASK,ncoo,coo_j,&incoo,coo_jl)); 3221 CHKERRQ(MatSetPreallocationCOO(a->A,ncoo,coo_il,coo_jl)); 3222 CHKERRQ(PetscFree2(coo_il,coo_jl)); 3223 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscNewLog(A,&a)); 3296 CHKERRQ(PetscStrallocpy(MATAIJ,&a->lmattype)); 3297 A->data = (void*)a; 3298 3299 /* matrix ops */ 3300 CHKERRQ(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 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",MatISSetLocalMatType_IS)); 3340 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS)); 3341 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS)); 3342 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS)); 3343 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ)); 3344 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS)); 3345 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS)); 3346 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",MatISFixLocalEmpty_IS)); 3347 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalToGlobalMapping_C",MatISGetLocalToGlobalMapping_IS)); 3348 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ)); 3349 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ)); 3350 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ)); 3351 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ)); 3352 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ)); 3353 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ)); 3354 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ)); 3355 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOOLocal_C",MatSetPreallocationCOOLocal_IS)); 3356 CHKERRQ(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_IS)); 3357 CHKERRQ(PetscObjectChangeTypeName((PetscObject)A,MATIS)); 3358 PetscFunctionReturn(0); 3359 } 3360