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