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->lmattype);CHKERRQ(ierr); 2316 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 2317 ierr = VecScatterDestroy(&b->cctx);CHKERRQ(ierr); 2318 ierr = VecScatterDestroy(&b->rctx);CHKERRQ(ierr); 2319 ierr = VecDestroy(&b->x);CHKERRQ(ierr); 2320 ierr = VecDestroy(&b->y);CHKERRQ(ierr); 2321 ierr = VecDestroy(&b->counter);CHKERRQ(ierr); 2322 ierr = ISDestroy(&b->getsub_ris);CHKERRQ(ierr); 2323 ierr = ISDestroy(&b->getsub_cis);CHKERRQ(ierr); 2324 if (b->sf != b->csf) { 2325 ierr = PetscSFDestroy(&b->csf);CHKERRQ(ierr); 2326 ierr = PetscFree2(b->csf_rootdata,b->csf_leafdata);CHKERRQ(ierr); 2327 } else b->csf = NULL; 2328 ierr = PetscSFDestroy(&b->sf);CHKERRQ(ierr); 2329 ierr = PetscFree2(b->sf_rootdata,b->sf_leafdata);CHKERRQ(ierr); 2330 ierr = PetscFree(A->data);CHKERRQ(ierr); 2331 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 2332 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",NULL);CHKERRQ(ierr); 2333 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);CHKERRQ(ierr); 2334 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);CHKERRQ(ierr); 2335 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);CHKERRQ(ierr); 2336 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);CHKERRQ(ierr); 2337 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",NULL);CHKERRQ(ierr); 2338 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL);CHKERRQ(ierr); 2339 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",NULL);CHKERRQ(ierr); 2340 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL);CHKERRQ(ierr); 2341 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL);CHKERRQ(ierr); 2342 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL);CHKERRQ(ierr); 2343 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL);CHKERRQ(ierr); 2344 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL);CHKERRQ(ierr); 2345 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL);CHKERRQ(ierr); 2346 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL);CHKERRQ(ierr); 2347 PetscFunctionReturn(0); 2348 } 2349 2350 static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y) 2351 { 2352 PetscErrorCode ierr; 2353 Mat_IS *is = (Mat_IS*)A->data; 2354 PetscScalar zero = 0.0; 2355 2356 PetscFunctionBegin; 2357 /* scatter the global vector x into the local work vector */ 2358 ierr = VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2359 ierr = VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2360 2361 /* multiply the local matrix */ 2362 ierr = MatMult(is->A,is->x,is->y);CHKERRQ(ierr); 2363 2364 /* scatter product back into global memory */ 2365 ierr = VecSet(y,zero);CHKERRQ(ierr); 2366 ierr = VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2367 ierr = VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2368 PetscFunctionReturn(0); 2369 } 2370 2371 static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 2372 { 2373 Vec temp_vec; 2374 PetscErrorCode ierr; 2375 2376 PetscFunctionBegin; /* v3 = v2 + A * v1.*/ 2377 if (v3 != v2) { 2378 ierr = MatMult(A,v1,v3);CHKERRQ(ierr); 2379 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 2380 } else { 2381 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 2382 ierr = MatMult(A,v1,temp_vec);CHKERRQ(ierr); 2383 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 2384 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 2385 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 2386 } 2387 PetscFunctionReturn(0); 2388 } 2389 2390 static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x) 2391 { 2392 Mat_IS *is = (Mat_IS*)A->data; 2393 PetscErrorCode ierr; 2394 2395 PetscFunctionBegin; 2396 /* scatter the global vector x into the local work vector */ 2397 ierr = VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2398 ierr = VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2399 2400 /* multiply the local matrix */ 2401 ierr = MatMultTranspose(is->A,is->y,is->x);CHKERRQ(ierr); 2402 2403 /* scatter product back into global vector */ 2404 ierr = VecSet(x,0);CHKERRQ(ierr); 2405 ierr = VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2406 ierr = VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2407 PetscFunctionReturn(0); 2408 } 2409 2410 static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 2411 { 2412 Vec temp_vec; 2413 PetscErrorCode ierr; 2414 2415 PetscFunctionBegin; /* v3 = v2 + A' * v1.*/ 2416 if (v3 != v2) { 2417 ierr = MatMultTranspose(A,v1,v3);CHKERRQ(ierr); 2418 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 2419 } else { 2420 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 2421 ierr = MatMultTranspose(A,v1,temp_vec);CHKERRQ(ierr); 2422 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 2423 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 2424 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 2425 } 2426 PetscFunctionReturn(0); 2427 } 2428 2429 static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer) 2430 { 2431 Mat_IS *a = (Mat_IS*)A->data; 2432 PetscErrorCode ierr; 2433 PetscViewer sviewer; 2434 PetscBool isascii,view = PETSC_TRUE; 2435 2436 PetscFunctionBegin; 2437 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 2438 if (isascii) { 2439 PetscViewerFormat format; 2440 2441 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 2442 if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE; 2443 } 2444 if (!view) PetscFunctionReturn(0); 2445 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 2446 ierr = MatView(a->A,sviewer);CHKERRQ(ierr); 2447 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 2448 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2449 PetscFunctionReturn(0); 2450 } 2451 2452 static PetscErrorCode MatISSetUpScatters_Private(Mat A) 2453 { 2454 Vec cglobal,rglobal; 2455 IS from; 2456 Mat_IS *is = (Mat_IS*)A->data; 2457 const PetscInt *garray; 2458 PetscInt nr,rbs,nc,cbs; 2459 PetscBool iscuda; 2460 PetscErrorCode ierr; 2461 2462 PetscFunctionBegin; 2463 ierr = ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nr);CHKERRQ(ierr); 2464 ierr = ISLocalToGlobalMappingGetBlockSize(A->rmap->mapping,&rbs);CHKERRQ(ierr); 2465 ierr = ISLocalToGlobalMappingGetSize(A->cmap->mapping,&nc);CHKERRQ(ierr); 2466 ierr = ISLocalToGlobalMappingGetBlockSize(A->cmap->mapping,&cbs);CHKERRQ(ierr); 2467 ierr = VecDestroy(&is->x);CHKERRQ(ierr); 2468 ierr = VecDestroy(&is->y);CHKERRQ(ierr); 2469 ierr = VecDestroy(&is->counter);CHKERRQ(ierr); 2470 ierr = VecScatterDestroy(&is->rctx);CHKERRQ(ierr); 2471 ierr = VecScatterDestroy(&is->cctx);CHKERRQ(ierr); 2472 ierr = MatCreateVecs(is->A,&is->x,&is->y);CHKERRQ(ierr); 2473 ierr = PetscObjectTypeCompare((PetscObject)is->y,VECSEQCUDA,&iscuda);CHKERRQ(ierr); 2474 if (iscuda) { 2475 ierr = PetscFree(A->defaultvectype);CHKERRQ(ierr); 2476 ierr = PetscStrallocpy(VECCUDA,&A->defaultvectype);CHKERRQ(ierr); 2477 } 2478 ierr = MatCreateVecs(A,&cglobal,&rglobal);CHKERRQ(ierr); 2479 ierr = ISLocalToGlobalMappingGetBlockIndices(A->rmap->mapping,&garray);CHKERRQ(ierr); 2480 ierr = ISCreateBlock(PetscObjectComm((PetscObject)A),rbs,nr/rbs,garray,PETSC_USE_POINTER,&from);CHKERRQ(ierr); 2481 ierr = VecScatterCreateWithData(rglobal,from,is->y,NULL,&is->rctx);CHKERRQ(ierr); 2482 ierr = ISLocalToGlobalMappingRestoreBlockIndices(A->rmap->mapping,&garray);CHKERRQ(ierr); 2483 ierr = ISDestroy(&from);CHKERRQ(ierr); 2484 if (A->rmap->mapping != A->cmap->mapping) { 2485 ierr = ISLocalToGlobalMappingGetBlockIndices(A->cmap->mapping,&garray);CHKERRQ(ierr); 2486 ierr = ISCreateBlock(PetscObjectComm((PetscObject)A),cbs,nc/cbs,garray,PETSC_USE_POINTER,&from);CHKERRQ(ierr); 2487 ierr = VecScatterCreateWithData(cglobal,from,is->x,NULL,&is->cctx);CHKERRQ(ierr); 2488 ierr = ISLocalToGlobalMappingRestoreBlockIndices(A->cmap->mapping,&garray);CHKERRQ(ierr); 2489 ierr = ISDestroy(&from);CHKERRQ(ierr); 2490 } else { 2491 ierr = PetscObjectReference((PetscObject)is->rctx);CHKERRQ(ierr); 2492 is->cctx = is->rctx; 2493 } 2494 /* interface counter vector (local) */ 2495 ierr = VecDuplicate(is->y,&is->counter);CHKERRQ(ierr); 2496 ierr = VecSet(is->y,1.);CHKERRQ(ierr); 2497 ierr = VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2498 ierr = VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2499 ierr = VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2500 ierr = VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2501 ierr = VecDestroy(&rglobal);CHKERRQ(ierr); 2502 ierr = VecDestroy(&cglobal);CHKERRQ(ierr); 2503 ierr = ISDestroy(&from);CHKERRQ(ierr); 2504 PetscFunctionReturn(0); 2505 } 2506 2507 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping) 2508 { 2509 PetscErrorCode ierr; 2510 PetscInt nr,rbs,nc,cbs; 2511 Mat_IS *is = (Mat_IS*)A->data; 2512 2513 PetscFunctionBegin; 2514 PetscCheckSameComm(A,1,rmapping,2); 2515 PetscCheckSameComm(A,1,cmapping,3); 2516 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 2517 if (is->csf != is->sf) { 2518 ierr = PetscSFDestroy(&is->csf);CHKERRQ(ierr); 2519 ierr = PetscFree2(is->csf_rootdata,is->csf_leafdata);CHKERRQ(ierr); 2520 } else is->csf = NULL; 2521 ierr = PetscSFDestroy(&is->sf);CHKERRQ(ierr); 2522 ierr = PetscFree2(is->sf_rootdata,is->sf_leafdata);CHKERRQ(ierr); 2523 2524 /* Setup Layout and set local to global maps */ 2525 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 2526 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 2527 ierr = ISLocalToGlobalMappingGetSize(rmapping,&nr);CHKERRQ(ierr); 2528 ierr = ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);CHKERRQ(ierr); 2529 ierr = ISLocalToGlobalMappingGetSize(cmapping,&nc);CHKERRQ(ierr); 2530 ierr = ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);CHKERRQ(ierr); 2531 /* check if the two mappings are actually the same for square matrices (DOLFIN passes 2 different objects) */ 2532 if (rmapping != cmapping && A->rmap->N == A->cmap->N) { 2533 PetscBool same,gsame; 2534 2535 same = PETSC_FALSE; 2536 if (nr == nc && cbs == rbs) { 2537 const PetscInt *idxs1,*idxs2; 2538 2539 ierr = ISLocalToGlobalMappingGetBlockIndices(rmapping,&idxs1);CHKERRQ(ierr); 2540 ierr = ISLocalToGlobalMappingGetBlockIndices(cmapping,&idxs2);CHKERRQ(ierr); 2541 ierr = PetscMemcmp(idxs1,idxs2,(nr/rbs)*sizeof(PetscInt),&same);CHKERRQ(ierr); 2542 ierr = ISLocalToGlobalMappingRestoreBlockIndices(rmapping,&idxs1);CHKERRQ(ierr); 2543 ierr = ISLocalToGlobalMappingRestoreBlockIndices(cmapping,&idxs2);CHKERRQ(ierr); 2544 } 2545 ierr = MPIU_Allreduce(&same,&gsame,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2546 if (gsame) cmapping = rmapping; 2547 } 2548 ierr = PetscLayoutSetBlockSize(A->rmap,rbs);CHKERRQ(ierr); 2549 ierr = PetscLayoutSetBlockSize(A->cmap,cbs);CHKERRQ(ierr); 2550 ierr = PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);CHKERRQ(ierr); 2551 ierr = PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);CHKERRQ(ierr); 2552 2553 /* Create the local matrix A */ 2554 ierr = MatCreate(PETSC_COMM_SELF,&is->A);CHKERRQ(ierr); 2555 ierr = MatSetType(is->A,is->lmattype);CHKERRQ(ierr); 2556 ierr = MatSetSizes(is->A,nr,nc,nr,nc);CHKERRQ(ierr); 2557 ierr = MatSetBlockSizes(is->A,rbs,cbs);CHKERRQ(ierr); 2558 ierr = MatSetOptionsPrefix(is->A,"is_");CHKERRQ(ierr); 2559 ierr = MatAppendOptionsPrefix(is->A,((PetscObject)A)->prefix);CHKERRQ(ierr); 2560 ierr = PetscLayoutSetUp(is->A->rmap);CHKERRQ(ierr); 2561 ierr = PetscLayoutSetUp(is->A->cmap);CHKERRQ(ierr); 2562 2563 if (!is->islocalref) { /* setup scatters and local vectors for MatMult */ 2564 ierr = MatISSetUpScatters_Private(A);CHKERRQ(ierr); 2565 } 2566 ierr = MatSetUp(A);CHKERRQ(ierr); 2567 PetscFunctionReturn(0); 2568 } 2569 2570 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2571 { 2572 Mat_IS *is = (Mat_IS*)mat->data; 2573 PetscErrorCode ierr; 2574 #if defined(PETSC_USE_DEBUG) 2575 PetscInt i,zm,zn; 2576 #endif 2577 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2578 2579 PetscFunctionBegin; 2580 #if defined(PETSC_USE_DEBUG) 2581 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); 2582 /* count negative indices */ 2583 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 2584 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 2585 #endif 2586 ierr = ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 2587 ierr = ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 2588 #if defined(PETSC_USE_DEBUG) 2589 /* count negative indices (should be the same as before) */ 2590 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 2591 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 2592 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"); 2593 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"); 2594 #endif 2595 ierr = MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 2596 PetscFunctionReturn(0); 2597 } 2598 2599 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2600 { 2601 Mat_IS *is = (Mat_IS*)mat->data; 2602 PetscErrorCode ierr; 2603 #if defined(PETSC_USE_DEBUG) 2604 PetscInt i,zm,zn; 2605 #endif 2606 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2607 2608 PetscFunctionBegin; 2609 #if defined(PETSC_USE_DEBUG) 2610 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); 2611 /* count negative indices */ 2612 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 2613 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 2614 #endif 2615 ierr = ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 2616 ierr = ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 2617 #if defined(PETSC_USE_DEBUG) 2618 /* count negative indices (should be the same as before) */ 2619 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 2620 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 2621 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"); 2622 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"); 2623 #endif 2624 ierr = MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 2625 PetscFunctionReturn(0); 2626 } 2627 2628 static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2629 { 2630 PetscErrorCode ierr; 2631 Mat_IS *is = (Mat_IS*)A->data; 2632 2633 PetscFunctionBegin; 2634 if (is->A->rmap->mapping) { 2635 ierr = MatSetValuesLocal(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2636 } else { 2637 ierr = MatSetValues(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2638 } 2639 PetscFunctionReturn(0); 2640 } 2641 2642 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2643 { 2644 PetscErrorCode ierr; 2645 Mat_IS *is = (Mat_IS*)A->data; 2646 2647 PetscFunctionBegin; 2648 if (is->A->rmap->mapping) { 2649 #if defined(PETSC_USE_DEBUG) 2650 PetscInt ibs,bs; 2651 2652 ierr = ISLocalToGlobalMappingGetBlockSize(is->A->rmap->mapping,&ibs);CHKERRQ(ierr); 2653 ierr = MatGetBlockSize(is->A,&bs);CHKERRQ(ierr); 2654 if (ibs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Different block sizes! local mat %D, local l2g map %D",bs,ibs); 2655 #endif 2656 ierr = MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2657 } else { 2658 ierr = MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2659 } 2660 PetscFunctionReturn(0); 2661 } 2662 2663 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns) 2664 { 2665 Mat_IS *is = (Mat_IS*)A->data; 2666 PetscErrorCode ierr; 2667 2668 PetscFunctionBegin; 2669 if (!n) { 2670 is->pure_neumann = PETSC_TRUE; 2671 } else { 2672 PetscInt i; 2673 is->pure_neumann = PETSC_FALSE; 2674 2675 if (columns) { 2676 ierr = MatZeroRowsColumns(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 2677 } else { 2678 ierr = MatZeroRows(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 2679 } 2680 if (diag != 0.) { 2681 const PetscScalar *array; 2682 ierr = VecGetArrayRead(is->counter,&array);CHKERRQ(ierr); 2683 for (i=0; i<n; i++) { 2684 ierr = MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);CHKERRQ(ierr); 2685 } 2686 ierr = VecRestoreArrayRead(is->counter,&array);CHKERRQ(ierr); 2687 } 2688 ierr = MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2689 ierr = MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2690 } 2691 PetscFunctionReturn(0); 2692 } 2693 2694 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns) 2695 { 2696 Mat_IS *matis = (Mat_IS*)A->data; 2697 PetscInt nr,nl,len,i; 2698 PetscInt *lrows; 2699 PetscErrorCode ierr; 2700 2701 PetscFunctionBegin; 2702 #if defined(PETSC_USE_DEBUG) 2703 if (columns || diag != 0. || (x && b)) { 2704 PetscBool cong; 2705 2706 ierr = PetscLayoutCompare(A->rmap,A->cmap,&cong);CHKERRQ(ierr); 2707 cong = (PetscBool)(cong && matis->sf == matis->csf); 2708 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"); 2709 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"); 2710 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"); 2711 } 2712 #endif 2713 /* get locally owned rows */ 2714 ierr = PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);CHKERRQ(ierr); 2715 /* fix right hand side if needed */ 2716 if (x && b) { 2717 const PetscScalar *xx; 2718 PetscScalar *bb; 2719 2720 ierr = VecGetArrayRead(x, &xx);CHKERRQ(ierr); 2721 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 2722 for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]]; 2723 ierr = VecRestoreArrayRead(x, &xx);CHKERRQ(ierr); 2724 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 2725 } 2726 /* get rows associated to the local matrices */ 2727 ierr = MatISSetUpSF(A);CHKERRQ(ierr); 2728 ierr = MatGetSize(matis->A,&nl,NULL);CHKERRQ(ierr); 2729 ierr = PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 2730 ierr = PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2731 for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1; 2732 ierr = PetscFree(lrows);CHKERRQ(ierr); 2733 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 2734 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 2735 ierr = PetscMalloc1(nl,&lrows);CHKERRQ(ierr); 2736 for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i; 2737 ierr = MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);CHKERRQ(ierr); 2738 ierr = PetscFree(lrows);CHKERRQ(ierr); 2739 PetscFunctionReturn(0); 2740 } 2741 2742 static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2743 { 2744 PetscErrorCode ierr; 2745 2746 PetscFunctionBegin; 2747 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);CHKERRQ(ierr); 2748 PetscFunctionReturn(0); 2749 } 2750 2751 static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2752 { 2753 PetscErrorCode ierr; 2754 2755 PetscFunctionBegin; 2756 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);CHKERRQ(ierr); 2757 PetscFunctionReturn(0); 2758 } 2759 2760 static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type) 2761 { 2762 Mat_IS *is = (Mat_IS*)A->data; 2763 PetscErrorCode ierr; 2764 2765 PetscFunctionBegin; 2766 ierr = MatAssemblyBegin(is->A,type);CHKERRQ(ierr); 2767 PetscFunctionReturn(0); 2768 } 2769 2770 static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type) 2771 { 2772 Mat_IS *is = (Mat_IS*)A->data; 2773 PetscErrorCode ierr; 2774 2775 PetscFunctionBegin; 2776 ierr = MatAssemblyEnd(is->A,type);CHKERRQ(ierr); 2777 /* fix for local empty rows/cols */ 2778 if (is->locempty && type == MAT_FINAL_ASSEMBLY) { 2779 Mat newlA; 2780 ISLocalToGlobalMapping rl2g,cl2g; 2781 IS nzr,nzc; 2782 PetscInt nr,nc,nnzr,nnzc; 2783 PetscBool lnewl2g,newl2g; 2784 2785 ierr = MatGetSize(is->A,&nr,&nc);CHKERRQ(ierr); 2786 ierr = MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_FALSE,PETSC_SMALL,&nzr);CHKERRQ(ierr); 2787 if (!nzr) { 2788 ierr = ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&nzr);CHKERRQ(ierr); 2789 } 2790 ierr = MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_TRUE,PETSC_SMALL,&nzc);CHKERRQ(ierr); 2791 if (!nzc) { 2792 ierr = ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&nzc);CHKERRQ(ierr); 2793 } 2794 ierr = ISGetSize(nzr,&nnzr);CHKERRQ(ierr); 2795 ierr = ISGetSize(nzc,&nnzc);CHKERRQ(ierr); 2796 if (nnzr != nr || nnzc != nc) { 2797 ISLocalToGlobalMapping l2g; 2798 IS is1,is2; 2799 2800 /* need new global l2g map */ 2801 lnewl2g = PETSC_TRUE; 2802 ierr = MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2803 2804 /* extract valid submatrix */ 2805 ierr = MatCreateSubMatrix(is->A,nzr,nzc,MAT_INITIAL_MATRIX,&newlA);CHKERRQ(ierr); 2806 2807 /* attach local l2g maps for successive calls of MatSetValues on the local matrix */ 2808 ierr = ISLocalToGlobalMappingCreateIS(nzr,&l2g);CHKERRQ(ierr); 2809 ierr = ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&is1);CHKERRQ(ierr); 2810 ierr = ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);CHKERRQ(ierr); 2811 ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr); 2812 if (is->A->rmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */ 2813 const PetscInt *idxs1,*idxs2; 2814 PetscInt j,i,nl,*tidxs; 2815 2816 ierr = ISLocalToGlobalMappingGetSize(is->A->rmap->mapping,&nl);CHKERRQ(ierr); 2817 ierr = ISLocalToGlobalMappingGetIndices(is->A->rmap->mapping,&idxs1);CHKERRQ(ierr); 2818 ierr = PetscMalloc1(nl,&tidxs);CHKERRQ(ierr); 2819 ierr = ISGetIndices(is2,&idxs2);CHKERRQ(ierr); 2820 for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++]; 2821 if (j != nr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nr); 2822 ierr = ISLocalToGlobalMappingRestoreIndices(is->A->rmap->mapping,&idxs1);CHKERRQ(ierr); 2823 ierr = ISRestoreIndices(is2,&idxs2);CHKERRQ(ierr); 2824 ierr = ISDestroy(&is2);CHKERRQ(ierr); 2825 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is->A->rmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);CHKERRQ(ierr); 2826 } 2827 ierr = ISLocalToGlobalMappingCreateIS(is2,&rl2g);CHKERRQ(ierr); 2828 ierr = ISDestroy(&is1);CHKERRQ(ierr); 2829 ierr = ISDestroy(&is2);CHKERRQ(ierr); 2830 2831 ierr = ISLocalToGlobalMappingCreateIS(nzc,&l2g);CHKERRQ(ierr); 2832 ierr = ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&is1);CHKERRQ(ierr); 2833 ierr = ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);CHKERRQ(ierr); 2834 ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr); 2835 if (is->A->cmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */ 2836 const PetscInt *idxs1,*idxs2; 2837 PetscInt j,i,nl,*tidxs; 2838 2839 ierr = ISLocalToGlobalMappingGetSize(is->A->cmap->mapping,&nl);CHKERRQ(ierr); 2840 ierr = ISLocalToGlobalMappingGetIndices(is->A->cmap->mapping,&idxs1);CHKERRQ(ierr); 2841 ierr = PetscMalloc1(nl,&tidxs);CHKERRQ(ierr); 2842 ierr = ISGetIndices(is2,&idxs2);CHKERRQ(ierr); 2843 for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++]; 2844 if (j != nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nc); 2845 ierr = ISLocalToGlobalMappingRestoreIndices(is->A->cmap->mapping,&idxs1);CHKERRQ(ierr); 2846 ierr = ISRestoreIndices(is2,&idxs2);CHKERRQ(ierr); 2847 ierr = ISDestroy(&is2);CHKERRQ(ierr); 2848 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is->A->cmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);CHKERRQ(ierr); 2849 } 2850 ierr = ISLocalToGlobalMappingCreateIS(is2,&cl2g);CHKERRQ(ierr); 2851 ierr = ISDestroy(&is1);CHKERRQ(ierr); 2852 ierr = ISDestroy(&is2);CHKERRQ(ierr); 2853 2854 ierr = MatSetLocalToGlobalMapping(newlA,rl2g,cl2g);CHKERRQ(ierr); 2855 2856 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 2857 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 2858 } else { /* local matrix fully populated */ 2859 lnewl2g = PETSC_FALSE; 2860 ierr = MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2861 ierr = PetscObjectReference((PetscObject)is->A);CHKERRQ(ierr); 2862 newlA = is->A; 2863 } 2864 /* attach new global l2g map if needed */ 2865 if (newl2g) { 2866 IS gnzr,gnzc; 2867 const PetscInt *grid,*gcid; 2868 2869 ierr = ISLocalToGlobalMappingApplyIS(A->rmap->mapping,nzr,&gnzr);CHKERRQ(ierr); 2870 ierr = ISLocalToGlobalMappingApplyIS(A->cmap->mapping,nzc,&gnzc);CHKERRQ(ierr); 2871 ierr = ISGetIndices(gnzr,&grid);CHKERRQ(ierr); 2872 ierr = ISGetIndices(gnzc,&gcid);CHKERRQ(ierr); 2873 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzr,grid,PETSC_COPY_VALUES,&rl2g);CHKERRQ(ierr); 2874 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzc,gcid,PETSC_COPY_VALUES,&cl2g);CHKERRQ(ierr); 2875 ierr = ISRestoreIndices(gnzr,&grid);CHKERRQ(ierr); 2876 ierr = ISRestoreIndices(gnzc,&gcid);CHKERRQ(ierr); 2877 ierr = ISDestroy(&gnzr);CHKERRQ(ierr); 2878 ierr = ISDestroy(&gnzc);CHKERRQ(ierr); 2879 ierr = MatSetLocalToGlobalMapping(A,rl2g,cl2g);CHKERRQ(ierr); 2880 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 2881 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 2882 } 2883 ierr = MatISSetLocalMat(A,newlA);CHKERRQ(ierr); 2884 ierr = MatDestroy(&newlA);CHKERRQ(ierr); 2885 ierr = ISDestroy(&nzr);CHKERRQ(ierr); 2886 ierr = ISDestroy(&nzc);CHKERRQ(ierr); 2887 is->locempty = PETSC_FALSE; 2888 } 2889 PetscFunctionReturn(0); 2890 } 2891 2892 static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local) 2893 { 2894 Mat_IS *is = (Mat_IS*)mat->data; 2895 2896 PetscFunctionBegin; 2897 *local = is->A; 2898 PetscFunctionReturn(0); 2899 } 2900 2901 static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local) 2902 { 2903 PetscFunctionBegin; 2904 *local = NULL; 2905 PetscFunctionReturn(0); 2906 } 2907 2908 /*@ 2909 MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix. 2910 2911 Input Parameter: 2912 . mat - the matrix 2913 2914 Output Parameter: 2915 . local - the local matrix 2916 2917 Level: advanced 2918 2919 Notes: 2920 This can be called if you have precomputed the nonzero structure of the 2921 matrix and want to provide it to the inner matrix object to improve the performance 2922 of the MatSetValues() operation. 2923 2924 Call MatISRestoreLocalMat() when finished with the local matrix. 2925 2926 .seealso: MATIS 2927 @*/ 2928 PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local) 2929 { 2930 PetscErrorCode ierr; 2931 2932 PetscFunctionBegin; 2933 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2934 PetscValidPointer(local,2); 2935 ierr = PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 2936 PetscFunctionReturn(0); 2937 } 2938 2939 /*@ 2940 MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat() 2941 2942 Input Parameter: 2943 . mat - the matrix 2944 2945 Output Parameter: 2946 . local - the local matrix 2947 2948 Level: advanced 2949 2950 .seealso: MATIS 2951 @*/ 2952 PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local) 2953 { 2954 PetscErrorCode ierr; 2955 2956 PetscFunctionBegin; 2957 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2958 PetscValidPointer(local,2); 2959 ierr = PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 2960 PetscFunctionReturn(0); 2961 } 2962 2963 static PetscErrorCode MatISSetLocalMatType_IS(Mat mat,MatType mtype) 2964 { 2965 Mat_IS *is = (Mat_IS*)mat->data; 2966 PetscErrorCode ierr; 2967 2968 PetscFunctionBegin; 2969 if (is->A) { 2970 ierr = MatSetType(is->A,mtype);CHKERRQ(ierr); 2971 } 2972 ierr = PetscFree(is->lmattype);CHKERRQ(ierr); 2973 ierr = PetscStrallocpy(mtype,&is->lmattype);CHKERRQ(ierr); 2974 PetscFunctionReturn(0); 2975 } 2976 2977 /*@ 2978 MatISSetLocalMatType - Specifies the type of local matrix 2979 2980 Input Parameter: 2981 . mat - the matrix 2982 . mtype - the local matrix type 2983 2984 Output Parameter: 2985 2986 Level: advanced 2987 2988 .seealso: MATIS, MatSetType(), MatType 2989 @*/ 2990 PetscErrorCode MatISSetLocalMatType(Mat mat,MatType mtype) 2991 { 2992 PetscErrorCode ierr; 2993 2994 PetscFunctionBegin; 2995 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2996 ierr = PetscUseMethod(mat,"MatISSetLocalMatType_C",(Mat,MatType),(mat,mtype));CHKERRQ(ierr); 2997 PetscFunctionReturn(0); 2998 } 2999 3000 static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local) 3001 { 3002 Mat_IS *is = (Mat_IS*)mat->data; 3003 PetscInt nrows,ncols,orows,ocols; 3004 PetscErrorCode ierr; 3005 MatType mtype,otype; 3006 PetscBool sametype = PETSC_TRUE; 3007 3008 PetscFunctionBegin; 3009 if (is->A) { 3010 ierr = MatGetSize(is->A,&orows,&ocols);CHKERRQ(ierr); 3011 ierr = MatGetSize(local,&nrows,&ncols);CHKERRQ(ierr); 3012 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); 3013 ierr = MatGetType(local,&mtype);CHKERRQ(ierr); 3014 ierr = MatGetType(is->A,&otype);CHKERRQ(ierr); 3015 ierr = PetscStrcmp(mtype,otype,&sametype);CHKERRQ(ierr); 3016 } 3017 ierr = PetscObjectReference((PetscObject)local);CHKERRQ(ierr); 3018 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 3019 is->A = local; 3020 ierr = MatGetType(is->A,&mtype);CHKERRQ(ierr); 3021 ierr = MatISSetLocalMatType(mat,mtype);CHKERRQ(ierr); 3022 if (!sametype && !is->islocalref) { 3023 ierr = MatISSetUpScatters_Private(mat);CHKERRQ(ierr); 3024 } 3025 PetscFunctionReturn(0); 3026 } 3027 3028 /*@ 3029 MatISSetLocalMat - Replace the local matrix stored inside a MATIS object. 3030 3031 Collective on Mat 3032 3033 Input Parameter: 3034 . mat - the matrix 3035 . local - the local matrix 3036 3037 Output Parameter: 3038 3039 Level: advanced 3040 3041 Notes: 3042 This can be called if you have precomputed the local matrix and 3043 want to provide it to the matrix object MATIS. 3044 3045 .seealso: MATIS 3046 @*/ 3047 PetscErrorCode MatISSetLocalMat(Mat mat,Mat local) 3048 { 3049 PetscErrorCode ierr; 3050 3051 PetscFunctionBegin; 3052 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3053 PetscValidHeaderSpecific(local,MAT_CLASSID,2); 3054 ierr = PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));CHKERRQ(ierr); 3055 PetscFunctionReturn(0); 3056 } 3057 3058 static PetscErrorCode MatZeroEntries_IS(Mat A) 3059 { 3060 Mat_IS *a = (Mat_IS*)A->data; 3061 PetscErrorCode ierr; 3062 3063 PetscFunctionBegin; 3064 ierr = MatZeroEntries(a->A);CHKERRQ(ierr); 3065 PetscFunctionReturn(0); 3066 } 3067 3068 static PetscErrorCode MatScale_IS(Mat A,PetscScalar a) 3069 { 3070 Mat_IS *is = (Mat_IS*)A->data; 3071 PetscErrorCode ierr; 3072 3073 PetscFunctionBegin; 3074 ierr = MatScale(is->A,a);CHKERRQ(ierr); 3075 PetscFunctionReturn(0); 3076 } 3077 3078 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v) 3079 { 3080 Mat_IS *is = (Mat_IS*)A->data; 3081 PetscErrorCode ierr; 3082 3083 PetscFunctionBegin; 3084 /* get diagonal of the local matrix */ 3085 ierr = MatGetDiagonal(is->A,is->y);CHKERRQ(ierr); 3086 3087 /* scatter diagonal back into global vector */ 3088 ierr = VecSet(v,0);CHKERRQ(ierr); 3089 ierr = VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 3090 ierr = VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 3091 PetscFunctionReturn(0); 3092 } 3093 3094 static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg) 3095 { 3096 Mat_IS *a = (Mat_IS*)A->data; 3097 PetscErrorCode ierr; 3098 3099 PetscFunctionBegin; 3100 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 3101 PetscFunctionReturn(0); 3102 } 3103 3104 static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str) 3105 { 3106 Mat_IS *y = (Mat_IS*)Y->data; 3107 Mat_IS *x; 3108 #if defined(PETSC_USE_DEBUG) 3109 PetscBool ismatis; 3110 #endif 3111 PetscErrorCode ierr; 3112 3113 PetscFunctionBegin; 3114 #if defined(PETSC_USE_DEBUG) 3115 ierr = PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);CHKERRQ(ierr); 3116 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS"); 3117 #endif 3118 x = (Mat_IS*)X->data; 3119 ierr = MatAXPY(y->A,a,x->A,str);CHKERRQ(ierr); 3120 PetscFunctionReturn(0); 3121 } 3122 3123 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat) 3124 { 3125 Mat lA; 3126 Mat_IS *matis; 3127 ISLocalToGlobalMapping rl2g,cl2g; 3128 IS is; 3129 const PetscInt *rg,*rl; 3130 PetscInt nrg; 3131 PetscInt N,M,nrl,i,*idxs; 3132 PetscErrorCode ierr; 3133 3134 PetscFunctionBegin; 3135 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 3136 ierr = ISGetLocalSize(row,&nrl);CHKERRQ(ierr); 3137 ierr = ISGetIndices(row,&rl);CHKERRQ(ierr); 3138 ierr = ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);CHKERRQ(ierr); 3139 #if defined(PETSC_USE_DEBUG) 3140 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); 3141 #endif 3142 ierr = PetscMalloc1(nrg,&idxs);CHKERRQ(ierr); 3143 /* map from [0,nrl) to row */ 3144 for (i=0;i<nrl;i++) idxs[i] = rl[i]; 3145 #if defined(PETSC_USE_DEBUG) 3146 for (i=nrl;i<nrg;i++) idxs[i] = nrg; 3147 #else 3148 for (i=nrl;i<nrg;i++) idxs[i] = -1; 3149 #endif 3150 ierr = ISRestoreIndices(row,&rl);CHKERRQ(ierr); 3151 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 3152 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 3153 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 3154 ierr = ISDestroy(&is);CHKERRQ(ierr); 3155 /* compute new l2g map for columns */ 3156 if (col != row || A->rmap->mapping != A->cmap->mapping) { 3157 const PetscInt *cg,*cl; 3158 PetscInt ncg; 3159 PetscInt ncl; 3160 3161 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 3162 ierr = ISGetLocalSize(col,&ncl);CHKERRQ(ierr); 3163 ierr = ISGetIndices(col,&cl);CHKERRQ(ierr); 3164 ierr = ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);CHKERRQ(ierr); 3165 #if defined(PETSC_USE_DEBUG) 3166 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); 3167 #endif 3168 ierr = PetscMalloc1(ncg,&idxs);CHKERRQ(ierr); 3169 /* map from [0,ncl) to col */ 3170 for (i=0;i<ncl;i++) idxs[i] = cl[i]; 3171 #if defined(PETSC_USE_DEBUG) 3172 for (i=ncl;i<ncg;i++) idxs[i] = ncg; 3173 #else 3174 for (i=ncl;i<ncg;i++) idxs[i] = -1; 3175 #endif 3176 ierr = ISRestoreIndices(col,&cl);CHKERRQ(ierr); 3177 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 3178 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 3179 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 3180 ierr = ISDestroy(&is);CHKERRQ(ierr); 3181 } else { 3182 ierr = PetscObjectReference((PetscObject)rl2g);CHKERRQ(ierr); 3183 cl2g = rl2g; 3184 } 3185 /* create the MATIS submatrix */ 3186 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 3187 ierr = MatCreate(PetscObjectComm((PetscObject)A),submat);CHKERRQ(ierr); 3188 ierr = MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3189 ierr = MatSetType(*submat,MATIS);CHKERRQ(ierr); 3190 matis = (Mat_IS*)((*submat)->data); 3191 matis->islocalref = PETSC_TRUE; 3192 ierr = MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);CHKERRQ(ierr); 3193 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr); 3194 ierr = MatISSetLocalMat(*submat,lA);CHKERRQ(ierr); 3195 ierr = MatSetUp(*submat);CHKERRQ(ierr); 3196 ierr = MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3197 ierr = MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3198 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 3199 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 3200 /* remove unsupported ops */ 3201 ierr = PetscMemzero((*submat)->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 3202 (*submat)->ops->destroy = MatDestroy_IS; 3203 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS; 3204 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS; 3205 (*submat)->ops->assemblybegin = MatAssemblyBegin_IS; 3206 (*submat)->ops->assemblyend = MatAssemblyEnd_IS; 3207 PetscFunctionReturn(0); 3208 } 3209 3210 static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A) 3211 { 3212 Mat_IS *a = (Mat_IS*)A->data; 3213 char type[256]; 3214 PetscBool flg; 3215 PetscErrorCode ierr; 3216 3217 PetscFunctionBegin; 3218 ierr = PetscOptionsHead(PetscOptionsObject,"MATIS options");CHKERRQ(ierr); 3219 ierr = PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns","MatISFixLocalEmpty",a->locempty,&a->locempty,NULL);CHKERRQ(ierr); 3220 ierr = PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL);CHKERRQ(ierr); 3221 ierr = PetscOptionsFList("-matis_localmat_type","Matrix type","MatISSetLocalMatType",MatList,a->lmattype,type,256,&flg);CHKERRQ(ierr); 3222 if (flg) { 3223 ierr = MatISSetLocalMatType(A,type);CHKERRQ(ierr); 3224 } 3225 if (a->A) { 3226 ierr = MatSetFromOptions(a->A);CHKERRQ(ierr); 3227 } 3228 ierr = PetscOptionsTail();CHKERRQ(ierr); 3229 PetscFunctionReturn(0); 3230 } 3231 3232 /*@ 3233 MatCreateIS - Creates a "process" unassembled matrix, assembled on each 3234 process but not across processes. 3235 3236 Input Parameters: 3237 + comm - MPI communicator that will share the matrix 3238 . bs - block size of the matrix 3239 . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products 3240 . rmap - local to global map for rows 3241 - cmap - local to global map for cols 3242 3243 Output Parameter: 3244 . A - the resulting matrix 3245 3246 Level: advanced 3247 3248 Notes: 3249 See MATIS for more details. 3250 m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors 3251 used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices. 3252 If either rmap or cmap are NULL, then the matrix is assumed to be square. 3253 3254 .seealso: MATIS, MatSetLocalToGlobalMapping() 3255 @*/ 3256 PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A) 3257 { 3258 PetscErrorCode ierr; 3259 3260 PetscFunctionBegin; 3261 if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings"); 3262 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3263 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 3264 if (bs > 0) { 3265 ierr = MatSetBlockSize(*A,bs);CHKERRQ(ierr); 3266 } 3267 ierr = MatSetType(*A,MATIS);CHKERRQ(ierr); 3268 if (rmap && cmap) { 3269 ierr = MatSetLocalToGlobalMapping(*A,rmap,cmap);CHKERRQ(ierr); 3270 } else if (!rmap) { 3271 ierr = MatSetLocalToGlobalMapping(*A,cmap,cmap);CHKERRQ(ierr); 3272 } else { 3273 ierr = MatSetLocalToGlobalMapping(*A,rmap,rmap);CHKERRQ(ierr); 3274 } 3275 PetscFunctionReturn(0); 3276 } 3277 3278 /*MC 3279 MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP). 3280 This stores the matrices in globally unassembled form. Each processor 3281 assembles only its local Neumann problem and the parallel matrix vector 3282 product is handled "implicitly". 3283 3284 Operations Provided: 3285 + MatMult() 3286 . MatMultAdd() 3287 . MatMultTranspose() 3288 . MatMultTransposeAdd() 3289 . MatZeroEntries() 3290 . MatSetOption() 3291 . MatZeroRows() 3292 . MatSetValues() 3293 . MatSetValuesBlocked() 3294 . MatSetValuesLocal() 3295 . MatSetValuesBlockedLocal() 3296 . MatScale() 3297 . MatGetDiagonal() 3298 . MatMissingDiagonal() 3299 . MatDuplicate() 3300 . MatCopy() 3301 . MatAXPY() 3302 . MatCreateSubMatrix() 3303 . MatGetLocalSubMatrix() 3304 . MatTranspose() 3305 . MatPtAP() (with P of AIJ type) 3306 - MatSetLocalToGlobalMapping() 3307 3308 Options Database Keys: 3309 + -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions() 3310 . -matis_fixempty - Fixes local matrices in case of empty local rows/columns. 3311 - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of MatPtAP(). 3312 3313 Notes: 3314 Options prefix for the inner matrix are given by -is_mat_xxx 3315 3316 You must call MatSetLocalToGlobalMapping() before using this matrix type. 3317 3318 You can do matrix preallocation on the local matrix after you obtain it with 3319 MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation() 3320 3321 Level: advanced 3322 3323 .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP 3324 3325 M*/ 3326 3327 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A) 3328 { 3329 PetscErrorCode ierr; 3330 Mat_IS *b; 3331 3332 PetscFunctionBegin; 3333 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 3334 ierr = PetscStrallocpy(MATAIJ,&b->lmattype);CHKERRQ(ierr); 3335 A->data = (void*)b; 3336 3337 /* matrix ops */ 3338 ierr = PetscMemzero(A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 3339 A->ops->mult = MatMult_IS; 3340 A->ops->multadd = MatMultAdd_IS; 3341 A->ops->multtranspose = MatMultTranspose_IS; 3342 A->ops->multtransposeadd = MatMultTransposeAdd_IS; 3343 A->ops->destroy = MatDestroy_IS; 3344 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS; 3345 A->ops->setvalues = MatSetValues_IS; 3346 A->ops->setvaluesblocked = MatSetValuesBlocked_IS; 3347 A->ops->setvalueslocal = MatSetValuesLocal_IS; 3348 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS; 3349 A->ops->zerorows = MatZeroRows_IS; 3350 A->ops->zerorowscolumns = MatZeroRowsColumns_IS; 3351 A->ops->assemblybegin = MatAssemblyBegin_IS; 3352 A->ops->assemblyend = MatAssemblyEnd_IS; 3353 A->ops->view = MatView_IS; 3354 A->ops->zeroentries = MatZeroEntries_IS; 3355 A->ops->scale = MatScale_IS; 3356 A->ops->getdiagonal = MatGetDiagonal_IS; 3357 A->ops->setoption = MatSetOption_IS; 3358 A->ops->ishermitian = MatIsHermitian_IS; 3359 A->ops->issymmetric = MatIsSymmetric_IS; 3360 A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS; 3361 A->ops->duplicate = MatDuplicate_IS; 3362 A->ops->missingdiagonal = MatMissingDiagonal_IS; 3363 A->ops->copy = MatCopy_IS; 3364 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 3365 A->ops->createsubmatrix = MatCreateSubMatrix_IS; 3366 A->ops->axpy = MatAXPY_IS; 3367 A->ops->diagonalset = MatDiagonalSet_IS; 3368 A->ops->shift = MatShift_IS; 3369 A->ops->transpose = MatTranspose_IS; 3370 A->ops->getinfo = MatGetInfo_IS; 3371 A->ops->diagonalscale = MatDiagonalScale_IS; 3372 A->ops->setfromoptions = MatSetFromOptions_IS; 3373 3374 /* special MATIS functions */ 3375 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",MatISSetLocalMatType_IS);CHKERRQ(ierr); 3376 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);CHKERRQ(ierr); 3377 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);CHKERRQ(ierr); 3378 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);CHKERRQ(ierr); 3379 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3380 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);CHKERRQ(ierr); 3381 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",MatISSetUpSF_IS);CHKERRQ(ierr); 3382 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS);CHKERRQ(ierr); 3383 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",MatISFixLocalEmpty_IS);CHKERRQ(ierr); 3384 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3385 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3386 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3387 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3388 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3389 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3390 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3391 ierr = PetscObjectChangeTypeName((PetscObject)A,MATIS);CHKERRQ(ierr); 3392 PetscFunctionReturn(0); 3393 } 3394