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