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 PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat) 803 { 804 Mat **nest,*snest,**rnest,lA,B; 805 IS *iscol,*isrow,*islrow,*islcol; 806 ISLocalToGlobalMapping rl2g,cl2g; 807 MPI_Comm comm; 808 PetscInt *lr,*lc,*l2gidxs; 809 PetscInt i,j,nr,nc,rbs,cbs; 810 PetscBool convert,lreuse,*istrans; 811 PetscErrorCode ierr; 812 813 PetscFunctionBegin; 814 ierr = MatNestGetSubMats(A,&nr,&nc,&nest);CHKERRQ(ierr); 815 lreuse = PETSC_FALSE; 816 rnest = NULL; 817 if (reuse == MAT_REUSE_MATRIX) { 818 PetscBool ismatis,isnest; 819 820 ierr = PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);CHKERRQ(ierr); 821 if (!ismatis) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type); 822 ierr = MatISGetLocalMat(*newmat,&lA);CHKERRQ(ierr); 823 ierr = PetscObjectTypeCompare((PetscObject)lA,MATNEST,&isnest);CHKERRQ(ierr); 824 if (isnest) { 825 ierr = MatNestGetSubMats(lA,&i,&j,&rnest);CHKERRQ(ierr); 826 lreuse = (PetscBool)(i == nr && j == nc); 827 if (!lreuse) rnest = NULL; 828 } 829 } 830 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 831 ierr = PetscCalloc2(nr,&lr,nc,&lc);CHKERRQ(ierr); 832 ierr = PetscCalloc6(nr,&isrow,nc,&iscol, 833 nr,&islrow,nc,&islcol, 834 nr*nc,&snest,nr*nc,&istrans);CHKERRQ(ierr); 835 ierr = MatNestGetISs(A,isrow,iscol);CHKERRQ(ierr); 836 for (i=0;i<nr;i++) { 837 for (j=0;j<nc;j++) { 838 PetscBool ismatis; 839 PetscInt l1,l2,lb1,lb2,ij=i*nc+j; 840 841 /* Null matrix pointers are allowed in MATNEST */ 842 if (!nest[i][j]) continue; 843 844 /* Nested matrices should be of type MATIS */ 845 ierr = PetscObjectTypeCompare((PetscObject)nest[i][j],MATTRANSPOSEMAT,&istrans[ij]);CHKERRQ(ierr); 846 if (istrans[ij]) { 847 Mat T,lT; 848 ierr = MatTransposeGetMat(nest[i][j],&T);CHKERRQ(ierr); 849 ierr = PetscObjectTypeCompare((PetscObject)T,MATIS,&ismatis);CHKERRQ(ierr); 850 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); 851 ierr = MatISGetLocalMat(T,&lT);CHKERRQ(ierr); 852 ierr = MatCreateTranspose(lT,&snest[ij]);CHKERRQ(ierr); 853 } else { 854 ierr = PetscObjectTypeCompare((PetscObject)nest[i][j],MATIS,&ismatis);CHKERRQ(ierr); 855 if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) is not of type MATIS",i,j); 856 ierr = MatISGetLocalMat(nest[i][j],&snest[ij]);CHKERRQ(ierr); 857 } 858 859 /* Check compatibility of local sizes */ 860 ierr = MatGetSize(snest[ij],&l1,&l2);CHKERRQ(ierr); 861 ierr = MatGetBlockSizes(snest[ij],&lb1,&lb2);CHKERRQ(ierr); 862 if (!l1 || !l2) continue; 863 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); 864 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); 865 lr[i] = l1; 866 lc[j] = l2; 867 868 /* check compatibilty for local matrix reusage */ 869 if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE; 870 } 871 } 872 873 #if defined (PETSC_USE_DEBUG) 874 /* Check compatibility of l2g maps for rows */ 875 for (i=0;i<nr;i++) { 876 rl2g = NULL; 877 for (j=0;j<nc;j++) { 878 PetscInt n1,n2; 879 880 if (!nest[i][j]) continue; 881 if (istrans[i*nc+j]) { 882 Mat T; 883 884 ierr = MatTransposeGetMat(nest[i][j],&T);CHKERRQ(ierr); 885 ierr = MatGetLocalToGlobalMapping(T,NULL,&cl2g);CHKERRQ(ierr); 886 } else { 887 ierr = MatGetLocalToGlobalMapping(nest[i][j],&cl2g,NULL);CHKERRQ(ierr); 888 } 889 ierr = ISLocalToGlobalMappingGetSize(cl2g,&n1);CHKERRQ(ierr); 890 if (!n1) continue; 891 if (!rl2g) { 892 rl2g = cl2g; 893 } else { 894 const PetscInt *idxs1,*idxs2; 895 PetscBool same; 896 897 ierr = ISLocalToGlobalMappingGetSize(rl2g,&n2);CHKERRQ(ierr); 898 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); 899 ierr = ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);CHKERRQ(ierr); 900 ierr = ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);CHKERRQ(ierr); 901 ierr = PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);CHKERRQ(ierr); 902 ierr = ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);CHKERRQ(ierr); 903 ierr = ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);CHKERRQ(ierr); 904 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); 905 } 906 } 907 } 908 /* Check compatibility of l2g maps for columns */ 909 for (i=0;i<nc;i++) { 910 rl2g = NULL; 911 for (j=0;j<nr;j++) { 912 PetscInt n1,n2; 913 914 if (!nest[j][i]) continue; 915 if (istrans[j*nc+i]) { 916 Mat T; 917 918 ierr = MatTransposeGetMat(nest[j][i],&T);CHKERRQ(ierr); 919 ierr = MatGetLocalToGlobalMapping(T,&cl2g,NULL);CHKERRQ(ierr); 920 } else { 921 ierr = MatGetLocalToGlobalMapping(nest[j][i],NULL,&cl2g);CHKERRQ(ierr); 922 } 923 ierr = ISLocalToGlobalMappingGetSize(cl2g,&n1);CHKERRQ(ierr); 924 if (!n1) continue; 925 if (!rl2g) { 926 rl2g = cl2g; 927 } else { 928 const PetscInt *idxs1,*idxs2; 929 PetscBool same; 930 931 ierr = ISLocalToGlobalMappingGetSize(rl2g,&n2);CHKERRQ(ierr); 932 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); 933 ierr = ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);CHKERRQ(ierr); 934 ierr = ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);CHKERRQ(ierr); 935 ierr = PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);CHKERRQ(ierr); 936 ierr = ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);CHKERRQ(ierr); 937 ierr = ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);CHKERRQ(ierr); 938 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); 939 } 940 } 941 } 942 #endif 943 944 B = NULL; 945 if (reuse != MAT_REUSE_MATRIX) { 946 PetscInt stl; 947 948 /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */ 949 for (i=0,stl=0;i<nr;i++) stl += lr[i]; 950 ierr = PetscMalloc1(stl,&l2gidxs);CHKERRQ(ierr); 951 for (i=0,stl=0;i<nr;i++) { 952 Mat usedmat; 953 Mat_IS *matis; 954 const PetscInt *idxs; 955 956 /* local IS for local NEST */ 957 ierr = ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i]);CHKERRQ(ierr); 958 959 /* l2gmap */ 960 j = 0; 961 usedmat = nest[i][j]; 962 while (!usedmat && j < nc-1) usedmat = nest[i][++j]; 963 if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid row mat"); 964 965 if (istrans[i*nc+j]) { 966 Mat T; 967 ierr = MatTransposeGetMat(usedmat,&T);CHKERRQ(ierr); 968 usedmat = T; 969 } 970 matis = (Mat_IS*)(usedmat->data); 971 ierr = ISGetIndices(isrow[i],&idxs);CHKERRQ(ierr); 972 if (istrans[i*nc+j]) { 973 ierr = PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl);CHKERRQ(ierr); 974 ierr = PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl);CHKERRQ(ierr); 975 } else { 976 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl);CHKERRQ(ierr); 977 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl);CHKERRQ(ierr); 978 } 979 ierr = ISRestoreIndices(isrow[i],&idxs);CHKERRQ(ierr); 980 stl += lr[i]; 981 } 982 ierr = ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&rl2g);CHKERRQ(ierr); 983 984 /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */ 985 for (i=0,stl=0;i<nc;i++) stl += lc[i]; 986 ierr = PetscMalloc1(stl,&l2gidxs);CHKERRQ(ierr); 987 for (i=0,stl=0;i<nc;i++) { 988 Mat usedmat; 989 Mat_IS *matis; 990 const PetscInt *idxs; 991 992 /* local IS for local NEST */ 993 ierr = ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i]);CHKERRQ(ierr); 994 995 /* l2gmap */ 996 j = 0; 997 usedmat = nest[j][i]; 998 while (!usedmat && j < nr-1) usedmat = nest[++j][i]; 999 if (!usedmat) SETERRQ(comm,PETSC_ERR_SUP,"Cannot find valid column mat"); 1000 if (istrans[j*nc+i]) { 1001 Mat T; 1002 ierr = MatTransposeGetMat(usedmat,&T);CHKERRQ(ierr); 1003 usedmat = T; 1004 } 1005 matis = (Mat_IS*)(usedmat->data); 1006 ierr = ISGetIndices(iscol[i],&idxs);CHKERRQ(ierr); 1007 if (istrans[j*nc+i]) { 1008 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl);CHKERRQ(ierr); 1009 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl);CHKERRQ(ierr); 1010 } else { 1011 ierr = PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl);CHKERRQ(ierr); 1012 ierr = PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl);CHKERRQ(ierr); 1013 } 1014 ierr = ISRestoreIndices(iscol[i],&idxs);CHKERRQ(ierr); 1015 stl += lc[i]; 1016 } 1017 ierr = ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&cl2g);CHKERRQ(ierr); 1018 1019 /* Create MATIS */ 1020 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 1021 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1022 ierr = MatGetBlockSizes(A,&rbs,&cbs);CHKERRQ(ierr); 1023 ierr = MatSetBlockSizes(B,rbs,cbs);CHKERRQ(ierr); 1024 ierr = MatSetType(B,MATIS);CHKERRQ(ierr); 1025 ierr = MatISSetLocalMatType(B,MATNEST);CHKERRQ(ierr); 1026 { /* hack : avoid setup of scatters */ 1027 Mat_IS *matis = (Mat_IS*)(B->data); 1028 matis->islocalref = PETSC_TRUE; 1029 } 1030 ierr = MatSetLocalToGlobalMapping(B,rl2g,cl2g);CHKERRQ(ierr); 1031 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 1032 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 1033 ierr = MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);CHKERRQ(ierr); 1034 ierr = MatNestSetVecType(lA,VECNEST);CHKERRQ(ierr); 1035 for (i=0;i<nr*nc;i++) { 1036 if (istrans[i]) { 1037 ierr = MatDestroy(&snest[i]);CHKERRQ(ierr); 1038 } 1039 } 1040 ierr = MatISSetLocalMat(B,lA);CHKERRQ(ierr); 1041 ierr = MatDestroy(&lA);CHKERRQ(ierr); 1042 { /* hack : setup of scatters done here */ 1043 Mat_IS *matis = (Mat_IS*)(B->data); 1044 1045 matis->islocalref = PETSC_FALSE; 1046 ierr = MatISSetUpScatters_Private(B);CHKERRQ(ierr); 1047 } 1048 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1049 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1050 if (reuse == MAT_INPLACE_MATRIX) { 1051 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 1052 } else { 1053 *newmat = B; 1054 } 1055 } else { 1056 if (lreuse) { 1057 ierr = MatISGetLocalMat(*newmat,&lA);CHKERRQ(ierr); 1058 for (i=0;i<nr;i++) { 1059 for (j=0;j<nc;j++) { 1060 if (snest[i*nc+j]) { 1061 ierr = MatNestSetSubMat(lA,i,j,snest[i*nc+j]);CHKERRQ(ierr); 1062 if (istrans[i*nc+j]) { 1063 ierr = MatDestroy(&snest[i*nc+j]);CHKERRQ(ierr); 1064 } 1065 } 1066 } 1067 } 1068 } else { 1069 PetscInt stl; 1070 for (i=0,stl=0;i<nr;i++) { 1071 ierr = ISCreateStride(PETSC_COMM_SELF,lr[i],stl,1,&islrow[i]);CHKERRQ(ierr); 1072 stl += lr[i]; 1073 } 1074 for (i=0,stl=0;i<nc;i++) { 1075 ierr = ISCreateStride(PETSC_COMM_SELF,lc[i],stl,1,&islcol[i]);CHKERRQ(ierr); 1076 stl += lc[i]; 1077 } 1078 ierr = MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);CHKERRQ(ierr); 1079 for (i=0;i<nr*nc;i++) { 1080 if (istrans[i]) { 1081 ierr = MatDestroy(&snest[i]);CHKERRQ(ierr); 1082 } 1083 } 1084 ierr = MatISSetLocalMat(*newmat,lA);CHKERRQ(ierr); 1085 ierr = MatDestroy(&lA);CHKERRQ(ierr); 1086 } 1087 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1088 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1089 } 1090 1091 /* Create local matrix in MATNEST format */ 1092 convert = PETSC_FALSE; 1093 ierr = PetscOptionsGetBool(NULL,((PetscObject)A)->prefix,"-matis_convert_local_nest",&convert,NULL);CHKERRQ(ierr); 1094 if (convert) { 1095 Mat M; 1096 MatISLocalFields lf; 1097 PetscContainer c; 1098 1099 ierr = MatISGetLocalMat(*newmat,&lA);CHKERRQ(ierr); 1100 ierr = MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M);CHKERRQ(ierr); 1101 ierr = MatISSetLocalMat(*newmat,M);CHKERRQ(ierr); 1102 ierr = MatDestroy(&M);CHKERRQ(ierr); 1103 1104 /* attach local fields to the matrix */ 1105 ierr = PetscNew(&lf);CHKERRQ(ierr); 1106 ierr = PetscCalloc2(nr,&lf->rf,nc,&lf->cf);CHKERRQ(ierr); 1107 for (i=0;i<nr;i++) { 1108 PetscInt n,st; 1109 1110 ierr = ISGetLocalSize(islrow[i],&n);CHKERRQ(ierr); 1111 ierr = ISStrideGetInfo(islrow[i],&st,NULL);CHKERRQ(ierr); 1112 ierr = ISCreateStride(comm,n,st,1,&lf->rf[i]);CHKERRQ(ierr); 1113 } 1114 for (i=0;i<nc;i++) { 1115 PetscInt n,st; 1116 1117 ierr = ISGetLocalSize(islcol[i],&n);CHKERRQ(ierr); 1118 ierr = ISStrideGetInfo(islcol[i],&st,NULL);CHKERRQ(ierr); 1119 ierr = ISCreateStride(comm,n,st,1,&lf->cf[i]);CHKERRQ(ierr); 1120 } 1121 lf->nr = nr; 1122 lf->nc = nc; 1123 ierr = PetscContainerCreate(PetscObjectComm((PetscObject)(*newmat)),&c);CHKERRQ(ierr); 1124 ierr = PetscContainerSetPointer(c,lf);CHKERRQ(ierr); 1125 ierr = PetscContainerSetUserDestroy(c,MatISContainerDestroyFields_Private);CHKERRQ(ierr); 1126 ierr = PetscObjectCompose((PetscObject)(*newmat),"_convert_nest_lfields",(PetscObject)c);CHKERRQ(ierr); 1127 ierr = PetscContainerDestroy(&c);CHKERRQ(ierr); 1128 } 1129 1130 /* Free workspace */ 1131 for (i=0;i<nr;i++) { 1132 ierr = ISDestroy(&islrow[i]);CHKERRQ(ierr); 1133 } 1134 for (i=0;i<nc;i++) { 1135 ierr = ISDestroy(&islcol[i]);CHKERRQ(ierr); 1136 } 1137 ierr = PetscFree6(isrow,iscol,islrow,islcol,snest,istrans);CHKERRQ(ierr); 1138 ierr = PetscFree2(lr,lc);CHKERRQ(ierr); 1139 PetscFunctionReturn(0); 1140 } 1141 1142 static PetscErrorCode MatDiagonalScale_IS(Mat A, Vec l, Vec r) 1143 { 1144 Mat_IS *matis = (Mat_IS*)A->data; 1145 Vec ll,rr; 1146 const PetscScalar *Y,*X; 1147 PetscScalar *x,*y; 1148 PetscErrorCode ierr; 1149 1150 PetscFunctionBegin; 1151 if (l) { 1152 ll = matis->y; 1153 ierr = VecGetArrayRead(l,&Y);CHKERRQ(ierr); 1154 ierr = VecGetArray(ll,&y);CHKERRQ(ierr); 1155 ierr = PetscSFBcastBegin(matis->sf,MPIU_SCALAR,Y,y);CHKERRQ(ierr); 1156 } else { 1157 ll = NULL; 1158 } 1159 if (r) { 1160 rr = matis->x; 1161 ierr = VecGetArrayRead(r,&X);CHKERRQ(ierr); 1162 ierr = VecGetArray(rr,&x);CHKERRQ(ierr); 1163 ierr = PetscSFBcastBegin(matis->csf,MPIU_SCALAR,X,x);CHKERRQ(ierr); 1164 } else { 1165 rr = NULL; 1166 } 1167 if (ll) { 1168 ierr = PetscSFBcastEnd(matis->sf,MPIU_SCALAR,Y,y);CHKERRQ(ierr); 1169 ierr = VecRestoreArrayRead(l,&Y);CHKERRQ(ierr); 1170 ierr = VecRestoreArray(ll,&y);CHKERRQ(ierr); 1171 } 1172 if (rr) { 1173 ierr = PetscSFBcastEnd(matis->csf,MPIU_SCALAR,X,x);CHKERRQ(ierr); 1174 ierr = VecRestoreArrayRead(r,&X);CHKERRQ(ierr); 1175 ierr = VecRestoreArray(rr,&x);CHKERRQ(ierr); 1176 } 1177 ierr = MatDiagonalScale(matis->A,ll,rr);CHKERRQ(ierr); 1178 PetscFunctionReturn(0); 1179 } 1180 1181 static PetscErrorCode MatGetInfo_IS(Mat A,MatInfoType flag,MatInfo *ginfo) 1182 { 1183 Mat_IS *matis = (Mat_IS*)A->data; 1184 MatInfo info; 1185 PetscReal isend[6],irecv[6]; 1186 PetscInt bs; 1187 PetscErrorCode ierr; 1188 1189 PetscFunctionBegin; 1190 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 1191 if (matis->A->ops->getinfo) { 1192 ierr = MatGetInfo(matis->A,MAT_LOCAL,&info);CHKERRQ(ierr); 1193 isend[0] = info.nz_used; 1194 isend[1] = info.nz_allocated; 1195 isend[2] = info.nz_unneeded; 1196 isend[3] = info.memory; 1197 isend[4] = info.mallocs; 1198 } else { 1199 isend[0] = 0.; 1200 isend[1] = 0.; 1201 isend[2] = 0.; 1202 isend[3] = 0.; 1203 isend[4] = 0.; 1204 } 1205 isend[5] = matis->A->num_ass; 1206 if (flag == MAT_LOCAL) { 1207 ginfo->nz_used = isend[0]; 1208 ginfo->nz_allocated = isend[1]; 1209 ginfo->nz_unneeded = isend[2]; 1210 ginfo->memory = isend[3]; 1211 ginfo->mallocs = isend[4]; 1212 ginfo->assemblies = isend[5]; 1213 } else if (flag == MAT_GLOBAL_MAX) { 1214 ierr = MPIU_Allreduce(isend,irecv,6,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1215 1216 ginfo->nz_used = irecv[0]; 1217 ginfo->nz_allocated = irecv[1]; 1218 ginfo->nz_unneeded = irecv[2]; 1219 ginfo->memory = irecv[3]; 1220 ginfo->mallocs = irecv[4]; 1221 ginfo->assemblies = irecv[5]; 1222 } else if (flag == MAT_GLOBAL_SUM) { 1223 ierr = MPIU_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1224 1225 ginfo->nz_used = irecv[0]; 1226 ginfo->nz_allocated = irecv[1]; 1227 ginfo->nz_unneeded = irecv[2]; 1228 ginfo->memory = irecv[3]; 1229 ginfo->mallocs = irecv[4]; 1230 ginfo->assemblies = A->num_ass; 1231 } 1232 ginfo->block_size = bs; 1233 ginfo->fill_ratio_given = 0; 1234 ginfo->fill_ratio_needed = 0; 1235 ginfo->factor_mallocs = 0; 1236 PetscFunctionReturn(0); 1237 } 1238 1239 PetscErrorCode MatTranspose_IS(Mat A,MatReuse reuse,Mat *B) 1240 { 1241 Mat C,lC,lA; 1242 PetscErrorCode ierr; 1243 1244 PetscFunctionBegin; 1245 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) { 1246 ISLocalToGlobalMapping rl2g,cl2g; 1247 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 1248 ierr = MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N);CHKERRQ(ierr); 1249 ierr = MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 1250 ierr = MatSetType(C,MATIS);CHKERRQ(ierr); 1251 ierr = MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);CHKERRQ(ierr); 1252 ierr = MatSetLocalToGlobalMapping(C,cl2g,rl2g);CHKERRQ(ierr); 1253 } else { 1254 C = *B; 1255 } 1256 1257 /* perform local transposition */ 1258 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr); 1259 ierr = MatTranspose(lA,MAT_INITIAL_MATRIX,&lC);CHKERRQ(ierr); 1260 ierr = MatISSetLocalMat(C,lC);CHKERRQ(ierr); 1261 ierr = MatDestroy(&lC);CHKERRQ(ierr); 1262 1263 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) { 1264 *B = C; 1265 } else { 1266 ierr = MatHeaderMerge(A,&C);CHKERRQ(ierr); 1267 } 1268 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1269 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1270 PetscFunctionReturn(0); 1271 } 1272 1273 PetscErrorCode MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode) 1274 { 1275 Mat_IS *is = (Mat_IS*)A->data; 1276 PetscErrorCode ierr; 1277 1278 PetscFunctionBegin; 1279 if (D) { /* MatShift_IS pass D = NULL */ 1280 ierr = VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1281 ierr = VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1282 } 1283 ierr = VecPointwiseDivide(is->y,is->y,is->counter);CHKERRQ(ierr); 1284 ierr = MatDiagonalSet(is->A,is->y,insmode);CHKERRQ(ierr); 1285 PetscFunctionReturn(0); 1286 } 1287 1288 PetscErrorCode MatShift_IS(Mat A,PetscScalar a) 1289 { 1290 Mat_IS *is = (Mat_IS*)A->data; 1291 PetscErrorCode ierr; 1292 1293 PetscFunctionBegin; 1294 ierr = VecSet(is->y,a);CHKERRQ(ierr); 1295 ierr = MatDiagonalSet_IS(A,NULL,ADD_VALUES);CHKERRQ(ierr); 1296 PetscFunctionReturn(0); 1297 } 1298 1299 static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 1300 { 1301 PetscErrorCode ierr; 1302 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1303 1304 PetscFunctionBegin; 1305 #if defined(PETSC_USE_DEBUG) 1306 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); 1307 #endif 1308 ierr = ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l);CHKERRQ(ierr); 1309 ierr = ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l);CHKERRQ(ierr); 1310 ierr = MatSetValuesLocal_IS(A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 1311 PetscFunctionReturn(0); 1312 } 1313 1314 static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 1315 { 1316 PetscErrorCode ierr; 1317 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1318 1319 PetscFunctionBegin; 1320 #if defined(PETSC_USE_DEBUG) 1321 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); 1322 #endif 1323 ierr = ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l);CHKERRQ(ierr); 1324 ierr = ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l);CHKERRQ(ierr); 1325 ierr = MatSetValuesBlockedLocal_IS(A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 1326 PetscFunctionReturn(0); 1327 } 1328 1329 static PetscErrorCode PetscLayoutMapLocal_Private(PetscLayout map,PetscInt N,const PetscInt idxs[], PetscInt *on,PetscInt **oidxs,PetscInt **ogidxs) 1330 { 1331 PetscInt *owners = map->range; 1332 PetscInt n = map->n; 1333 PetscSF sf; 1334 PetscInt *lidxs,*work = NULL; 1335 PetscSFNode *ridxs; 1336 PetscMPIInt rank; 1337 PetscInt r, p = 0, len = 0; 1338 PetscErrorCode ierr; 1339 1340 PetscFunctionBegin; 1341 if (on) *on = 0; /* squelch -Wmaybe-uninitialized */ 1342 /* Create SF where leaves are input idxs and roots are owned idxs (code adapted from MatZeroRowsMapLocal_Private) */ 1343 ierr = MPI_Comm_rank(map->comm,&rank);CHKERRQ(ierr); 1344 ierr = PetscMalloc1(n,&lidxs);CHKERRQ(ierr); 1345 for (r = 0; r < n; ++r) lidxs[r] = -1; 1346 ierr = PetscMalloc1(N,&ridxs);CHKERRQ(ierr); 1347 for (r = 0; r < N; ++r) { 1348 const PetscInt idx = idxs[r]; 1349 if (idx < 0 || map->N <= idx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D out of range [0,%D)",idx,map->N); 1350 if (idx < owners[p] || owners[p+1] <= idx) { /* short-circuit the search if the last p owns this idx too */ 1351 ierr = PetscLayoutFindOwner(map,idx,&p);CHKERRQ(ierr); 1352 } 1353 ridxs[r].rank = p; 1354 ridxs[r].index = idxs[r] - owners[p]; 1355 } 1356 ierr = PetscSFCreate(map->comm,&sf);CHKERRQ(ierr); 1357 ierr = PetscSFSetGraph(sf,n,N,NULL,PETSC_OWN_POINTER,ridxs,PETSC_OWN_POINTER);CHKERRQ(ierr); 1358 ierr = PetscSFReduceBegin(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);CHKERRQ(ierr); 1359 ierr = PetscSFReduceEnd(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);CHKERRQ(ierr); 1360 if (ogidxs) { /* communicate global idxs */ 1361 PetscInt cum = 0,start,*work2; 1362 1363 ierr = PetscMalloc1(n,&work);CHKERRQ(ierr); 1364 ierr = PetscCalloc1(N,&work2);CHKERRQ(ierr); 1365 for (r = 0; r < N; ++r) if (idxs[r] >=0) cum++; 1366 ierr = MPI_Scan(&cum,&start,1,MPIU_INT,MPI_SUM,map->comm);CHKERRQ(ierr); 1367 start -= cum; 1368 cum = 0; 1369 for (r = 0; r < N; ++r) if (idxs[r] >=0) work2[r] = start+cum++; 1370 ierr = PetscSFReduceBegin(sf,MPIU_INT,work2,work,MPIU_REPLACE);CHKERRQ(ierr); 1371 ierr = PetscSFReduceEnd(sf,MPIU_INT,work2,work,MPIU_REPLACE);CHKERRQ(ierr); 1372 ierr = PetscFree(work2);CHKERRQ(ierr); 1373 } 1374 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 1375 /* Compress and put in indices */ 1376 for (r = 0; r < n; ++r) 1377 if (lidxs[r] >= 0) { 1378 if (work) work[len] = work[r]; 1379 lidxs[len++] = r; 1380 } 1381 if (on) *on = len; 1382 if (oidxs) *oidxs = lidxs; 1383 if (ogidxs) *ogidxs = work; 1384 PetscFunctionReturn(0); 1385 } 1386 1387 static PetscErrorCode MatCreateSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat) 1388 { 1389 Mat locmat,newlocmat; 1390 Mat_IS *newmatis; 1391 #if defined(PETSC_USE_DEBUG) 1392 Vec rtest,ltest; 1393 const PetscScalar *array; 1394 #endif 1395 const PetscInt *idxs; 1396 PetscInt i,m,n; 1397 PetscErrorCode ierr; 1398 1399 PetscFunctionBegin; 1400 if (scall == MAT_REUSE_MATRIX) { 1401 PetscBool ismatis; 1402 1403 ierr = PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);CHKERRQ(ierr); 1404 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Not of MATIS type"); 1405 newmatis = (Mat_IS*)(*newmat)->data; 1406 if (!newmatis->getsub_ris) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local row IS"); 1407 if (!newmatis->getsub_cis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local col IS"); 1408 } 1409 /* irow and icol may not have duplicate entries */ 1410 #if defined(PETSC_USE_DEBUG) 1411 ierr = MatCreateVecs(mat,<est,&rtest);CHKERRQ(ierr); 1412 ierr = ISGetLocalSize(irow,&n);CHKERRQ(ierr); 1413 ierr = ISGetIndices(irow,&idxs);CHKERRQ(ierr); 1414 for (i=0;i<n;i++) { 1415 ierr = VecSetValue(rtest,idxs[i],1.0,ADD_VALUES);CHKERRQ(ierr); 1416 } 1417 ierr = VecAssemblyBegin(rtest);CHKERRQ(ierr); 1418 ierr = VecAssemblyEnd(rtest);CHKERRQ(ierr); 1419 ierr = VecGetLocalSize(rtest,&n);CHKERRQ(ierr); 1420 ierr = VecGetOwnershipRange(rtest,&m,NULL);CHKERRQ(ierr); 1421 ierr = VecGetArrayRead(rtest,&array);CHKERRQ(ierr); 1422 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])); 1423 ierr = VecRestoreArrayRead(rtest,&array);CHKERRQ(ierr); 1424 ierr = ISRestoreIndices(irow,&idxs);CHKERRQ(ierr); 1425 ierr = ISGetLocalSize(icol,&n);CHKERRQ(ierr); 1426 ierr = ISGetIndices(icol,&idxs);CHKERRQ(ierr); 1427 for (i=0;i<n;i++) { 1428 ierr = VecSetValue(ltest,idxs[i],1.0,ADD_VALUES);CHKERRQ(ierr); 1429 } 1430 ierr = VecAssemblyBegin(ltest);CHKERRQ(ierr); 1431 ierr = VecAssemblyEnd(ltest);CHKERRQ(ierr); 1432 ierr = VecGetLocalSize(ltest,&n);CHKERRQ(ierr); 1433 ierr = VecGetOwnershipRange(ltest,&m,NULL);CHKERRQ(ierr); 1434 ierr = VecGetArrayRead(ltest,&array);CHKERRQ(ierr); 1435 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])); 1436 ierr = VecRestoreArrayRead(ltest,&array);CHKERRQ(ierr); 1437 ierr = ISRestoreIndices(icol,&idxs);CHKERRQ(ierr); 1438 ierr = VecDestroy(&rtest);CHKERRQ(ierr); 1439 ierr = VecDestroy(<est);CHKERRQ(ierr); 1440 #endif 1441 if (scall == MAT_INITIAL_MATRIX) { 1442 Mat_IS *matis = (Mat_IS*)mat->data; 1443 ISLocalToGlobalMapping rl2g; 1444 IS is; 1445 PetscInt *lidxs,*lgidxs,*newgidxs; 1446 PetscInt ll,newloc,irbs,icbs,arbs,acbs,rbs,cbs; 1447 PetscBool cong; 1448 MPI_Comm comm; 1449 1450 ierr = PetscObjectGetComm((PetscObject)mat,&comm);CHKERRQ(ierr); 1451 ierr = MatGetBlockSizes(mat,&arbs,&acbs);CHKERRQ(ierr); 1452 ierr = ISGetBlockSize(irow,&irbs);CHKERRQ(ierr); 1453 ierr = ISGetBlockSize(icol,&icbs);CHKERRQ(ierr); 1454 rbs = arbs == irbs ? irbs : 1; 1455 cbs = acbs == icbs ? icbs : 1; 1456 ierr = ISGetLocalSize(irow,&m);CHKERRQ(ierr); 1457 ierr = ISGetLocalSize(icol,&n);CHKERRQ(ierr); 1458 ierr = MatCreate(comm,newmat);CHKERRQ(ierr); 1459 ierr = MatSetType(*newmat,MATIS);CHKERRQ(ierr); 1460 ierr = MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 1461 ierr = MatSetBlockSizes(*newmat,rbs,cbs);CHKERRQ(ierr); 1462 /* communicate irow to their owners in the layout */ 1463 ierr = ISGetIndices(irow,&idxs);CHKERRQ(ierr); 1464 ierr = PetscLayoutMapLocal_Private(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs);CHKERRQ(ierr); 1465 ierr = ISRestoreIndices(irow,&idxs);CHKERRQ(ierr); 1466 ierr = PetscMemzero(matis->sf_rootdata,matis->sf->nroots*sizeof(PetscInt));CHKERRQ(ierr); 1467 for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1; 1468 ierr = PetscFree(lidxs);CHKERRQ(ierr); 1469 ierr = PetscFree(lgidxs);CHKERRQ(ierr); 1470 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1471 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1472 for (i=0,newloc=0;i<matis->sf->nleaves;i++) if (matis->sf_leafdata[i]) newloc++; 1473 ierr = PetscMalloc1(newloc,&newgidxs);CHKERRQ(ierr); 1474 ierr = PetscMalloc1(newloc,&lidxs);CHKERRQ(ierr); 1475 for (i=0,newloc=0;i<matis->sf->nleaves;i++) 1476 if (matis->sf_leafdata[i]) { 1477 lidxs[newloc] = i; 1478 newgidxs[newloc++] = matis->sf_leafdata[i]-1; 1479 } 1480 ierr = ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 1481 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 1482 ierr = ISLocalToGlobalMappingSetBlockSize(rl2g,rbs);CHKERRQ(ierr); 1483 ierr = ISDestroy(&is);CHKERRQ(ierr); 1484 /* local is to extract local submatrix */ 1485 newmatis = (Mat_IS*)(*newmat)->data; 1486 ierr = ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris);CHKERRQ(ierr); 1487 ierr = MatHasCongruentLayouts(mat,&cong);CHKERRQ(ierr); 1488 if (cong && irow == icol && matis->csf == matis->sf) { 1489 ierr = MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g);CHKERRQ(ierr); 1490 ierr = PetscObjectReference((PetscObject)newmatis->getsub_ris);CHKERRQ(ierr); 1491 newmatis->getsub_cis = newmatis->getsub_ris; 1492 } else { 1493 ISLocalToGlobalMapping cl2g; 1494 1495 /* communicate icol to their owners in the layout */ 1496 ierr = ISGetIndices(icol,&idxs);CHKERRQ(ierr); 1497 ierr = PetscLayoutMapLocal_Private(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs);CHKERRQ(ierr); 1498 ierr = ISRestoreIndices(icol,&idxs);CHKERRQ(ierr); 1499 ierr = PetscMemzero(matis->csf_rootdata,matis->csf->nroots*sizeof(PetscInt));CHKERRQ(ierr); 1500 for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1; 1501 ierr = PetscFree(lidxs);CHKERRQ(ierr); 1502 ierr = PetscFree(lgidxs);CHKERRQ(ierr); 1503 ierr = PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);CHKERRQ(ierr); 1504 ierr = PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);CHKERRQ(ierr); 1505 for (i=0,newloc=0;i<matis->csf->nleaves;i++) if (matis->csf_leafdata[i]) newloc++; 1506 ierr = PetscMalloc1(newloc,&newgidxs);CHKERRQ(ierr); 1507 ierr = PetscMalloc1(newloc,&lidxs);CHKERRQ(ierr); 1508 for (i=0,newloc=0;i<matis->csf->nleaves;i++) 1509 if (matis->csf_leafdata[i]) { 1510 lidxs[newloc] = i; 1511 newgidxs[newloc++] = matis->csf_leafdata[i]-1; 1512 } 1513 ierr = ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 1514 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 1515 ierr = ISLocalToGlobalMappingSetBlockSize(cl2g,cbs);CHKERRQ(ierr); 1516 ierr = ISDestroy(&is);CHKERRQ(ierr); 1517 /* local is to extract local submatrix */ 1518 ierr = ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis);CHKERRQ(ierr); 1519 ierr = MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);CHKERRQ(ierr); 1520 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 1521 } 1522 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 1523 } else { 1524 ierr = MatISGetLocalMat(*newmat,&newlocmat);CHKERRQ(ierr); 1525 } 1526 ierr = MatISGetLocalMat(mat,&locmat);CHKERRQ(ierr); 1527 newmatis = (Mat_IS*)(*newmat)->data; 1528 ierr = MatCreateSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat);CHKERRQ(ierr); 1529 if (scall == MAT_INITIAL_MATRIX) { 1530 ierr = MatISSetLocalMat(*newmat,newlocmat);CHKERRQ(ierr); 1531 ierr = MatDestroy(&newlocmat);CHKERRQ(ierr); 1532 } 1533 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1534 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1535 PetscFunctionReturn(0); 1536 } 1537 1538 static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str) 1539 { 1540 Mat_IS *a = (Mat_IS*)A->data,*b; 1541 PetscBool ismatis; 1542 PetscErrorCode ierr; 1543 1544 PetscFunctionBegin; 1545 ierr = PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);CHKERRQ(ierr); 1546 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Need to be implemented"); 1547 b = (Mat_IS*)B->data; 1548 ierr = MatCopy(a->A,b->A,str);CHKERRQ(ierr); 1549 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 1550 PetscFunctionReturn(0); 1551 } 1552 1553 static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool *missing,PetscInt *d) 1554 { 1555 Vec v; 1556 const PetscScalar *array; 1557 PetscInt i,n; 1558 PetscErrorCode ierr; 1559 1560 PetscFunctionBegin; 1561 *missing = PETSC_FALSE; 1562 ierr = MatCreateVecs(A,NULL,&v);CHKERRQ(ierr); 1563 ierr = MatGetDiagonal(A,v);CHKERRQ(ierr); 1564 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 1565 ierr = VecGetArrayRead(v,&array);CHKERRQ(ierr); 1566 for (i=0;i<n;i++) if (array[i] == 0.) break; 1567 ierr = VecRestoreArrayRead(v,&array);CHKERRQ(ierr); 1568 ierr = VecDestroy(&v);CHKERRQ(ierr); 1569 if (i != n) *missing = PETSC_TRUE; 1570 if (d) { 1571 *d = -1; 1572 if (*missing) { 1573 PetscInt rstart; 1574 ierr = MatGetOwnershipRange(A,&rstart,NULL);CHKERRQ(ierr); 1575 *d = i+rstart; 1576 } 1577 } 1578 PetscFunctionReturn(0); 1579 } 1580 1581 static PetscErrorCode MatISSetUpSF_IS(Mat B) 1582 { 1583 Mat_IS *matis = (Mat_IS*)(B->data); 1584 const PetscInt *gidxs; 1585 PetscInt nleaves; 1586 PetscErrorCode ierr; 1587 1588 PetscFunctionBegin; 1589 if (matis->sf) PetscFunctionReturn(0); 1590 ierr = PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);CHKERRQ(ierr); 1591 ierr = ISLocalToGlobalMappingGetIndices(B->rmap->mapping,&gidxs);CHKERRQ(ierr); 1592 ierr = ISLocalToGlobalMappingGetSize(B->rmap->mapping,&nleaves);CHKERRQ(ierr); 1593 ierr = PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);CHKERRQ(ierr); 1594 ierr = ISLocalToGlobalMappingRestoreIndices(B->rmap->mapping,&gidxs);CHKERRQ(ierr); 1595 ierr = PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata);CHKERRQ(ierr); 1596 if (B->rmap->mapping != B->cmap->mapping) { /* setup SF for columns */ 1597 ierr = ISLocalToGlobalMappingGetSize(B->cmap->mapping,&nleaves);CHKERRQ(ierr); 1598 ierr = PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);CHKERRQ(ierr); 1599 ierr = ISLocalToGlobalMappingGetIndices(B->cmap->mapping,&gidxs);CHKERRQ(ierr); 1600 ierr = PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);CHKERRQ(ierr); 1601 ierr = ISLocalToGlobalMappingRestoreIndices(B->cmap->mapping,&gidxs);CHKERRQ(ierr); 1602 ierr = PetscMalloc2(matis->csf->nroots,&matis->csf_rootdata,matis->csf->nleaves,&matis->csf_leafdata);CHKERRQ(ierr); 1603 } else { 1604 matis->csf = matis->sf; 1605 matis->csf_leafdata = matis->sf_leafdata; 1606 matis->csf_rootdata = matis->sf_rootdata; 1607 } 1608 PetscFunctionReturn(0); 1609 } 1610 1611 /*@ 1612 MatISStoreL2L - Store local-to-local operators during the Galerkin process of MatPtAP. 1613 1614 Collective on MPI_Comm 1615 1616 Input Parameters: 1617 + A - the matrix 1618 - store - the boolean flag 1619 1620 Level: advanced 1621 1622 Notes: 1623 1624 .keywords: matrix 1625 1626 .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatPtAP() 1627 @*/ 1628 PetscErrorCode MatISStoreL2L(Mat A, PetscBool store) 1629 { 1630 PetscErrorCode ierr; 1631 1632 PetscFunctionBegin; 1633 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1634 PetscValidType(A,1); 1635 PetscValidLogicalCollectiveBool(A,store,2); 1636 ierr = PetscTryMethod(A,"MatISStoreL2L_C",(Mat,PetscBool),(A,store));CHKERRQ(ierr); 1637 PetscFunctionReturn(0); 1638 } 1639 1640 static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store) 1641 { 1642 Mat_IS *matis = (Mat_IS*)(A->data); 1643 PetscErrorCode ierr; 1644 1645 PetscFunctionBegin; 1646 matis->storel2l = store; 1647 if (!store) { 1648 ierr = PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",NULL);CHKERRQ(ierr); 1649 } 1650 PetscFunctionReturn(0); 1651 } 1652 1653 /*@ 1654 MatISFixLocalEmpty - Compress out zero local rows from the local matrices 1655 1656 Collective on MPI_Comm 1657 1658 Input Parameters: 1659 + A - the matrix 1660 - fix - the boolean flag 1661 1662 Level: advanced 1663 1664 Notes: When fix is true, new local matrices and l2g maps are generated during the final assembly process. 1665 1666 .keywords: matrix 1667 1668 .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatAssemblyEnd(), MAT_FINAL_ASSEMBLY 1669 @*/ 1670 PetscErrorCode MatISFixLocalEmpty(Mat A, PetscBool fix) 1671 { 1672 PetscErrorCode ierr; 1673 1674 PetscFunctionBegin; 1675 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1676 PetscValidType(A,1); 1677 PetscValidLogicalCollectiveBool(A,fix,2); 1678 ierr = PetscTryMethod(A,"MatISFixLocalEmpty_C",(Mat,PetscBool),(A,fix));CHKERRQ(ierr); 1679 PetscFunctionReturn(0); 1680 } 1681 1682 static PetscErrorCode MatISFixLocalEmpty_IS(Mat A, PetscBool fix) 1683 { 1684 Mat_IS *matis = (Mat_IS*)(A->data); 1685 1686 PetscFunctionBegin; 1687 matis->locempty = fix; 1688 PetscFunctionReturn(0); 1689 } 1690 1691 /*@ 1692 MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix. 1693 1694 Collective on MPI_Comm 1695 1696 Input Parameters: 1697 + B - the matrix 1698 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1699 (same value is used for all local rows) 1700 . d_nnz - array containing the number of nonzeros in the various rows of the 1701 DIAGONAL portion of the local submatrix (possibly different for each row) 1702 or NULL, if d_nz is used to specify the nonzero structure. 1703 The size of this array is equal to the number of local rows, i.e 'm'. 1704 For matrices that will be factored, you must leave room for (and set) 1705 the diagonal entry even if it is zero. 1706 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1707 submatrix (same value is used for all local rows). 1708 - o_nnz - array containing the number of nonzeros in the various rows of the 1709 OFF-DIAGONAL portion of the local submatrix (possibly different for 1710 each row) or NULL, if o_nz is used to specify the nonzero 1711 structure. The size of this array is equal to the number 1712 of local rows, i.e 'm'. 1713 1714 If the *_nnz parameter is given then the *_nz parameter is ignored 1715 1716 Level: intermediate 1717 1718 Notes: 1719 This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition 1720 from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local 1721 matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects. 1722 1723 .keywords: matrix 1724 1725 .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS 1726 @*/ 1727 PetscErrorCode MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1728 { 1729 PetscErrorCode ierr; 1730 1731 PetscFunctionBegin; 1732 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1733 PetscValidType(B,1); 1734 ierr = PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));CHKERRQ(ierr); 1735 PetscFunctionReturn(0); 1736 } 1737 1738 /* this is used by DMDA */ 1739 PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1740 { 1741 Mat_IS *matis = (Mat_IS*)(B->data); 1742 PetscInt bs,i,nlocalcols; 1743 PetscErrorCode ierr; 1744 1745 PetscFunctionBegin; 1746 if (!matis->A) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"You should first call MatSetLocalToGlobalMapping"); 1747 1748 if (!d_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nz; 1749 else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nnz[i]; 1750 1751 if (!o_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nz; 1752 else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nnz[i]; 1753 1754 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1755 ierr = MatGetSize(matis->A,NULL,&nlocalcols);CHKERRQ(ierr); 1756 ierr = MatGetBlockSize(matis->A,&bs);CHKERRQ(ierr); 1757 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1758 1759 for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols); 1760 ierr = MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);CHKERRQ(ierr); 1761 #if defined(PETSC_HAVE_HYPRE) 1762 ierr = MatHYPRESetPreallocation(matis->A,0,matis->sf_leafdata,0,NULL);CHKERRQ(ierr); 1763 #endif 1764 1765 for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs; 1766 ierr = MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);CHKERRQ(ierr); 1767 1768 nlocalcols /= bs; 1769 for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols - i); 1770 ierr = MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);CHKERRQ(ierr); 1771 1772 /* for other matrix types */ 1773 ierr = MatSetUp(matis->A);CHKERRQ(ierr); 1774 PetscFunctionReturn(0); 1775 } 1776 1777 PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce) 1778 { 1779 Mat_IS *matis = (Mat_IS*)(A->data); 1780 PetscInt *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership; 1781 const PetscInt *global_indices_r,*global_indices_c; 1782 PetscInt i,j,bs,rows,cols; 1783 PetscInt lrows,lcols; 1784 PetscInt local_rows,local_cols; 1785 PetscMPIInt size; 1786 PetscBool isdense,issbaij; 1787 PetscErrorCode ierr; 1788 1789 PetscFunctionBegin; 1790 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr); 1791 ierr = MatGetSize(A,&rows,&cols);CHKERRQ(ierr); 1792 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 1793 ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr); 1794 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);CHKERRQ(ierr); 1795 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);CHKERRQ(ierr); 1796 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);CHKERRQ(ierr); 1797 if (A->rmap->mapping != A->cmap->mapping) { 1798 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&global_indices_c);CHKERRQ(ierr); 1799 } else { 1800 global_indices_c = global_indices_r; 1801 } 1802 1803 if (issbaij) { 1804 ierr = MatGetRowUpperTriangular(matis->A);CHKERRQ(ierr); 1805 } 1806 /* 1807 An SF reduce is needed to sum up properly on shared rows. 1808 Note that generally preallocation is not exact, since it overestimates nonzeros 1809 */ 1810 ierr = MatGetLocalSize(A,&lrows,&lcols);CHKERRQ(ierr); 1811 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);CHKERRQ(ierr); 1812 /* All processes need to compute entire row ownership */ 1813 ierr = PetscMalloc1(rows,&row_ownership);CHKERRQ(ierr); 1814 ierr = MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);CHKERRQ(ierr); 1815 for (i=0;i<size;i++) { 1816 for (j=mat_ranges[i];j<mat_ranges[i+1];j++) { 1817 row_ownership[j] = i; 1818 } 1819 } 1820 ierr = MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);CHKERRQ(ierr); 1821 1822 /* 1823 my_dnz and my_onz contains exact contribution to preallocation from each local mat 1824 then, they will be summed up properly. This way, preallocation is always sufficient 1825 */ 1826 ierr = PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);CHKERRQ(ierr); 1827 /* preallocation as a MATAIJ */ 1828 if (isdense) { /* special case for dense local matrices */ 1829 for (i=0;i<local_rows;i++) { 1830 PetscInt owner = row_ownership[global_indices_r[i]]; 1831 for (j=0;j<local_cols;j++) { 1832 PetscInt index_col = global_indices_c[j]; 1833 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 1834 my_dnz[i] += 1; 1835 } else { /* offdiag block */ 1836 my_onz[i] += 1; 1837 } 1838 } 1839 } 1840 } else if (matis->A->ops->getrowij) { 1841 const PetscInt *ii,*jj,*jptr; 1842 PetscBool done; 1843 ierr = MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);CHKERRQ(ierr); 1844 if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ"); 1845 jptr = jj; 1846 for (i=0;i<local_rows;i++) { 1847 PetscInt index_row = global_indices_r[i]; 1848 for (j=0;j<ii[i+1]-ii[i];j++,jptr++) { 1849 PetscInt owner = row_ownership[index_row]; 1850 PetscInt index_col = global_indices_c[*jptr]; 1851 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 1852 my_dnz[i] += 1; 1853 } else { /* offdiag block */ 1854 my_onz[i] += 1; 1855 } 1856 /* same as before, interchanging rows and cols */ 1857 if (issbaij && index_col != index_row) { 1858 owner = row_ownership[index_col]; 1859 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) { 1860 my_dnz[*jptr] += 1; 1861 } else { 1862 my_onz[*jptr] += 1; 1863 } 1864 } 1865 } 1866 } 1867 ierr = MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);CHKERRQ(ierr); 1868 if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ"); 1869 } else { /* loop over rows and use MatGetRow */ 1870 for (i=0;i<local_rows;i++) { 1871 const PetscInt *cols; 1872 PetscInt ncols,index_row = global_indices_r[i]; 1873 ierr = MatGetRow(matis->A,i,&ncols,&cols,NULL);CHKERRQ(ierr); 1874 for (j=0;j<ncols;j++) { 1875 PetscInt owner = row_ownership[index_row]; 1876 PetscInt index_col = global_indices_c[cols[j]]; 1877 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 1878 my_dnz[i] += 1; 1879 } else { /* offdiag block */ 1880 my_onz[i] += 1; 1881 } 1882 /* same as before, interchanging rows and cols */ 1883 if (issbaij && index_col != index_row) { 1884 owner = row_ownership[index_col]; 1885 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) { 1886 my_dnz[cols[j]] += 1; 1887 } else { 1888 my_onz[cols[j]] += 1; 1889 } 1890 } 1891 } 1892 ierr = MatRestoreRow(matis->A,i,&ncols,&cols,NULL);CHKERRQ(ierr); 1893 } 1894 } 1895 if (global_indices_c != global_indices_r) { 1896 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&global_indices_c);CHKERRQ(ierr); 1897 } 1898 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);CHKERRQ(ierr); 1899 ierr = PetscFree(row_ownership);CHKERRQ(ierr); 1900 1901 /* Reduce my_dnz and my_onz */ 1902 if (maxreduce) { 1903 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);CHKERRQ(ierr); 1904 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);CHKERRQ(ierr); 1905 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);CHKERRQ(ierr); 1906 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);CHKERRQ(ierr); 1907 } else { 1908 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);CHKERRQ(ierr); 1909 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);CHKERRQ(ierr); 1910 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);CHKERRQ(ierr); 1911 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);CHKERRQ(ierr); 1912 } 1913 ierr = PetscFree2(my_dnz,my_onz);CHKERRQ(ierr); 1914 1915 /* Resize preallocation if overestimated */ 1916 for (i=0;i<lrows;i++) { 1917 dnz[i] = PetscMin(dnz[i],lcols); 1918 onz[i] = PetscMin(onz[i],cols-lcols); 1919 } 1920 1921 /* Set preallocation */ 1922 ierr = MatSeqAIJSetPreallocation(B,0,dnz);CHKERRQ(ierr); 1923 ierr = MatMPIAIJSetPreallocation(B,0,dnz,0,onz);CHKERRQ(ierr); 1924 for (i=0;i<lrows;i+=bs) { 1925 PetscInt b, d = dnz[i],o = onz[i]; 1926 1927 for (b=1;b<bs;b++) { 1928 d = PetscMax(d,dnz[i+b]); 1929 o = PetscMax(o,onz[i+b]); 1930 } 1931 dnz[i/bs] = PetscMin(d/bs + d%bs,lcols/bs); 1932 onz[i/bs] = PetscMin(o/bs + o%bs,(cols-lcols)/bs); 1933 } 1934 ierr = MatSeqBAIJSetPreallocation(B,bs,0,dnz);CHKERRQ(ierr); 1935 ierr = MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 1936 ierr = MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 1937 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1938 if (issbaij) { 1939 ierr = MatRestoreRowUpperTriangular(matis->A);CHKERRQ(ierr); 1940 } 1941 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1942 PetscFunctionReturn(0); 1943 } 1944 1945 PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M) 1946 { 1947 Mat_IS *matis = (Mat_IS*)(mat->data); 1948 Mat local_mat,MT; 1949 PetscInt rbs,cbs,rows,cols,lrows,lcols; 1950 PetscInt local_rows,local_cols; 1951 PetscBool isseqdense,isseqsbaij,isseqaij,isseqbaij; 1952 #if defined (PETSC_USE_DEBUG) 1953 PetscBool lb[4],bb[4]; 1954 #endif 1955 PetscMPIInt size; 1956 PetscScalar *array; 1957 PetscErrorCode ierr; 1958 1959 PetscFunctionBegin; 1960 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);CHKERRQ(ierr); 1961 if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) { 1962 Mat B; 1963 IS irows = NULL,icols = NULL; 1964 PetscInt rbs,cbs; 1965 1966 ierr = ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);CHKERRQ(ierr); 1967 ierr = ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);CHKERRQ(ierr); 1968 if (reuse != MAT_REUSE_MATRIX) { /* check if l2g maps are one-to-one */ 1969 IS rows,cols; 1970 const PetscInt *ridxs,*cidxs; 1971 PetscInt i,nw,*work; 1972 1973 ierr = ISLocalToGlobalMappingGetBlockIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 1974 ierr = ISLocalToGlobalMappingGetSize(mat->rmap->mapping,&nw);CHKERRQ(ierr); 1975 nw = nw/rbs; 1976 ierr = PetscCalloc1(nw,&work);CHKERRQ(ierr); 1977 for (i=0;i<nw;i++) work[ridxs[i]] += 1; 1978 for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break; 1979 if (i == nw) { 1980 ierr = ISCreateBlock(PETSC_COMM_SELF,rbs,nw,ridxs,PETSC_USE_POINTER,&rows);CHKERRQ(ierr); 1981 ierr = ISSetPermutation(rows);CHKERRQ(ierr); 1982 ierr = ISInvertPermutation(rows,PETSC_DECIDE,&irows);CHKERRQ(ierr); 1983 ierr = ISDestroy(&rows);CHKERRQ(ierr); 1984 } 1985 ierr = ISLocalToGlobalMappingRestoreBlockIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 1986 ierr = PetscFree(work);CHKERRQ(ierr); 1987 if (irows && mat->rmap->mapping != mat->cmap->mapping) { 1988 ierr = ISLocalToGlobalMappingGetBlockIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 1989 ierr = ISLocalToGlobalMappingGetSize(mat->cmap->mapping,&nw);CHKERRQ(ierr); 1990 nw = nw/cbs; 1991 ierr = PetscCalloc1(nw,&work);CHKERRQ(ierr); 1992 for (i=0;i<nw;i++) work[cidxs[i]] += 1; 1993 for (i=0;i<nw;i++) if (!work[i] || work[i] > 1) break; 1994 if (i == nw) { 1995 ierr = ISCreateBlock(PETSC_COMM_SELF,cbs,nw,cidxs,PETSC_USE_POINTER,&cols);CHKERRQ(ierr); 1996 ierr = ISSetPermutation(cols);CHKERRQ(ierr); 1997 ierr = ISInvertPermutation(cols,PETSC_DECIDE,&icols);CHKERRQ(ierr); 1998 ierr = ISDestroy(&cols);CHKERRQ(ierr); 1999 } 2000 ierr = ISLocalToGlobalMappingRestoreBlockIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 2001 ierr = PetscFree(work);CHKERRQ(ierr); 2002 } else if (irows) { 2003 ierr = PetscObjectReference((PetscObject)irows);CHKERRQ(ierr); 2004 icols = irows; 2005 } 2006 } else { 2007 ierr = PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject*)&irows);CHKERRQ(ierr); 2008 ierr = PetscObjectQuery((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject*)&icols);CHKERRQ(ierr); 2009 if (irows) { ierr = PetscObjectReference((PetscObject)irows);CHKERRQ(ierr); } 2010 if (icols) { ierr = PetscObjectReference((PetscObject)icols);CHKERRQ(ierr); } 2011 } 2012 if (!irows || !icols) { 2013 ierr = ISDestroy(&icols);CHKERRQ(ierr); 2014 ierr = ISDestroy(&irows);CHKERRQ(ierr); 2015 goto general_assembly; 2016 } 2017 ierr = MatConvert(matis->A,(rbs == cbs && rbs > 1) ? MATSEQBAIJ : MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 2018 if (reuse != MAT_INPLACE_MATRIX) { 2019 ierr = MatCreateSubMatrix(B,irows,icols,reuse,M);CHKERRQ(ierr); 2020 ierr = PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_irows",(PetscObject)irows);CHKERRQ(ierr); 2021 ierr = PetscObjectCompose((PetscObject)(*M),"_MatIS_IS_XAIJ_icols",(PetscObject)icols);CHKERRQ(ierr); 2022 } else { 2023 Mat C; 2024 2025 ierr = MatCreateSubMatrix(B,irows,icols,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 2026 ierr = MatHeaderReplace(mat,&C);CHKERRQ(ierr); 2027 } 2028 ierr = MatDestroy(&B);CHKERRQ(ierr); 2029 ierr = ISDestroy(&icols);CHKERRQ(ierr); 2030 ierr = ISDestroy(&irows);CHKERRQ(ierr); 2031 PetscFunctionReturn(0); 2032 } 2033 general_assembly: 2034 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 2035 ierr = ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);CHKERRQ(ierr); 2036 ierr = ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);CHKERRQ(ierr); 2037 ierr = MatGetLocalSize(mat,&lrows,&lcols);CHKERRQ(ierr); 2038 ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr); 2039 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 2040 ierr = PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 2041 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);CHKERRQ(ierr); 2042 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);CHKERRQ(ierr); 2043 if (!isseqdense && !isseqaij && !isseqbaij && !isseqsbaij) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name); 2044 #if defined (PETSC_USE_DEBUG) 2045 lb[0] = isseqdense; 2046 lb[1] = isseqaij; 2047 lb[2] = isseqbaij; 2048 lb[3] = isseqsbaij; 2049 ierr = MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 2050 if (!bb[0] && !bb[1] && !bb[2] && !bb[3]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type"); 2051 #endif 2052 2053 if (reuse != MAT_REUSE_MATRIX) { 2054 ierr = MatCreate(PetscObjectComm((PetscObject)mat),&MT);CHKERRQ(ierr); 2055 ierr = MatSetSizes(MT,lrows,lcols,rows,cols);CHKERRQ(ierr); 2056 ierr = MatSetType(MT,mtype);CHKERRQ(ierr); 2057 ierr = MatSetBlockSizes(MT,rbs,cbs);CHKERRQ(ierr); 2058 ierr = MatISSetMPIXAIJPreallocation_Private(mat,MT,PETSC_FALSE);CHKERRQ(ierr); 2059 } else { 2060 PetscInt mrbs,mcbs,mrows,mcols,mlrows,mlcols; 2061 2062 /* some checks */ 2063 MT = *M; 2064 ierr = MatGetBlockSizes(MT,&mrbs,&mcbs);CHKERRQ(ierr); 2065 ierr = MatGetSize(MT,&mrows,&mcols);CHKERRQ(ierr); 2066 ierr = MatGetLocalSize(MT,&mlrows,&mlcols);CHKERRQ(ierr); 2067 if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows); 2068 if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols); 2069 if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows); 2070 if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols); 2071 if (mrbs != rbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong row block size (%d != %d)",rbs,mrbs); 2072 if (mcbs != cbs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong col block size (%d != %d)",cbs,mcbs); 2073 ierr = MatZeroEntries(MT);CHKERRQ(ierr); 2074 } 2075 2076 if (isseqsbaij || isseqbaij) { 2077 ierr = MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&local_mat);CHKERRQ(ierr); 2078 isseqaij = PETSC_TRUE; 2079 } else { 2080 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 2081 local_mat = matis->A; 2082 } 2083 2084 /* Set values */ 2085 ierr = MatSetLocalToGlobalMapping(MT,mat->rmap->mapping,mat->cmap->mapping);CHKERRQ(ierr); 2086 if (isseqdense) { /* special case for dense local matrices */ 2087 PetscInt i,*dummy; 2088 2089 ierr = PetscMalloc1(PetscMax(local_rows,local_cols),&dummy);CHKERRQ(ierr); 2090 for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i; 2091 ierr = MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 2092 ierr = MatDenseGetArray(local_mat,&array);CHKERRQ(ierr); 2093 ierr = MatSetValuesLocal(MT,local_rows,dummy,local_cols,dummy,array,ADD_VALUES);CHKERRQ(ierr); 2094 ierr = MatDenseRestoreArray(local_mat,&array);CHKERRQ(ierr); 2095 ierr = PetscFree(dummy);CHKERRQ(ierr); 2096 } else if (isseqaij) { 2097 const PetscInt *blocks; 2098 PetscInt i,nvtxs,*xadj,*adjncy, nb; 2099 PetscBool done; 2100 2101 ierr = MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 2102 if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ"); 2103 ierr = MatSeqAIJGetArray(local_mat,&array);CHKERRQ(ierr); 2104 ierr = MatGetVariableBlockSizes(local_mat,&nb,&blocks);CHKERRQ(ierr); 2105 if (nb) { /* speed up assembly for special blocked matrices (used by BDDC) */ 2106 PetscInt sum; 2107 2108 for (i=0,sum=0;i<nb;i++) sum += blocks[i]; 2109 if (sum == nvtxs) { 2110 PetscInt r; 2111 2112 for (i=0,r=0;i<nb;i++) { 2113 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]); 2114 ierr = MatSetValuesLocal(MT,blocks[i],adjncy+xadj[r],blocks[i],adjncy+xadj[r],array+xadj[r],ADD_VALUES);CHKERRQ(ierr); 2115 r += blocks[i]; 2116 } 2117 } else { 2118 for (i=0;i<nvtxs;i++) { 2119 ierr = MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);CHKERRQ(ierr); 2120 } 2121 } 2122 } else { 2123 for (i=0;i<nvtxs;i++) { 2124 ierr = MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);CHKERRQ(ierr); 2125 } 2126 } 2127 ierr = MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 2128 if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ"); 2129 ierr = MatSeqAIJRestoreArray(local_mat,&array);CHKERRQ(ierr); 2130 } else { /* very basic values insertion for all other matrix types */ 2131 PetscInt i; 2132 2133 for (i=0;i<local_rows;i++) { 2134 PetscInt j; 2135 const PetscInt *local_indices_cols; 2136 2137 ierr = MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 2138 ierr = MatSetValuesLocal(MT,1,&i,j,local_indices_cols,array,ADD_VALUES);CHKERRQ(ierr); 2139 ierr = MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 2140 } 2141 } 2142 ierr = MatAssemblyBegin(MT,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2143 ierr = MatDestroy(&local_mat);CHKERRQ(ierr); 2144 ierr = MatAssemblyEnd(MT,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2145 if (isseqdense) { 2146 ierr = MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 2147 } 2148 if (reuse == MAT_INPLACE_MATRIX) { 2149 ierr = MatHeaderReplace(mat,&MT);CHKERRQ(ierr); 2150 } else if (reuse == MAT_INITIAL_MATRIX) { 2151 *M = MT; 2152 } 2153 PetscFunctionReturn(0); 2154 } 2155 2156 /*@ 2157 MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format 2158 2159 Input Parameter: 2160 . mat - the matrix (should be of type MATIS) 2161 . reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 2162 2163 Output Parameter: 2164 . newmat - the matrix in AIJ format 2165 2166 Level: developer 2167 2168 Notes: 2169 This function has been deprecated and it will be removed in future releases. Update your code to use the MatConvert() interface. 2170 2171 .seealso: MATIS, MatConvert() 2172 @*/ 2173 PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat) 2174 { 2175 PetscErrorCode ierr; 2176 2177 PetscFunctionBegin; 2178 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2179 PetscValidLogicalCollectiveEnum(mat,reuse,2); 2180 PetscValidPointer(newmat,3); 2181 if (reuse == MAT_REUSE_MATRIX) { 2182 PetscValidHeaderSpecific(*newmat,MAT_CLASSID,3); 2183 PetscCheckSameComm(mat,1,*newmat,3); 2184 if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix"); 2185 } 2186 ierr = PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatType,MatReuse,Mat*),(mat,MATAIJ,reuse,newmat));CHKERRQ(ierr); 2187 PetscFunctionReturn(0); 2188 } 2189 2190 PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat) 2191 { 2192 PetscErrorCode ierr; 2193 Mat_IS *matis = (Mat_IS*)(mat->data); 2194 PetscInt rbs,cbs,m,n,M,N; 2195 Mat B,localmat; 2196 2197 PetscFunctionBegin; 2198 ierr = ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);CHKERRQ(ierr); 2199 ierr = ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);CHKERRQ(ierr); 2200 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 2201 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 2202 ierr = MatCreate(PetscObjectComm((PetscObject)mat),&B);CHKERRQ(ierr); 2203 ierr = MatSetSizes(B,m,n,M,N);CHKERRQ(ierr); 2204 ierr = MatSetBlockSize(B,rbs == cbs ? rbs : 1);CHKERRQ(ierr); 2205 ierr = MatSetType(B,MATIS);CHKERRQ(ierr); 2206 ierr = MatISSetLocalMatType(B,matis->lmattype);CHKERRQ(ierr); 2207 ierr = MatSetLocalToGlobalMapping(B,mat->rmap->mapping,mat->cmap->mapping);CHKERRQ(ierr); 2208 ierr = MatDuplicate(matis->A,op,&localmat);CHKERRQ(ierr); 2209 ierr = MatISSetLocalMat(B,localmat);CHKERRQ(ierr); 2210 ierr = MatDestroy(&localmat);CHKERRQ(ierr); 2211 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2212 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2213 *newmat = B; 2214 PetscFunctionReturn(0); 2215 } 2216 2217 static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg) 2218 { 2219 PetscErrorCode ierr; 2220 Mat_IS *matis = (Mat_IS*)A->data; 2221 PetscBool local_sym; 2222 2223 PetscFunctionBegin; 2224 ierr = MatIsHermitian(matis->A,tol,&local_sym);CHKERRQ(ierr); 2225 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2226 PetscFunctionReturn(0); 2227 } 2228 2229 static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg) 2230 { 2231 PetscErrorCode ierr; 2232 Mat_IS *matis = (Mat_IS*)A->data; 2233 PetscBool local_sym; 2234 2235 PetscFunctionBegin; 2236 ierr = MatIsSymmetric(matis->A,tol,&local_sym);CHKERRQ(ierr); 2237 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2238 PetscFunctionReturn(0); 2239 } 2240 2241 static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool *flg) 2242 { 2243 PetscErrorCode ierr; 2244 Mat_IS *matis = (Mat_IS*)A->data; 2245 PetscBool local_sym; 2246 2247 PetscFunctionBegin; 2248 if (A->rmap->mapping != A->cmap->mapping) { 2249 *flg = PETSC_FALSE; 2250 PetscFunctionReturn(0); 2251 } 2252 ierr = MatIsStructurallySymmetric(matis->A,&local_sym);CHKERRQ(ierr); 2253 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2254 PetscFunctionReturn(0); 2255 } 2256 2257 static PetscErrorCode MatDestroy_IS(Mat A) 2258 { 2259 PetscErrorCode ierr; 2260 Mat_IS *b = (Mat_IS*)A->data; 2261 2262 PetscFunctionBegin; 2263 ierr = PetscFree(b->bdiag);CHKERRQ(ierr); 2264 ierr = PetscFree(b->lmattype);CHKERRQ(ierr); 2265 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 2266 ierr = VecScatterDestroy(&b->cctx);CHKERRQ(ierr); 2267 ierr = VecScatterDestroy(&b->rctx);CHKERRQ(ierr); 2268 ierr = VecDestroy(&b->x);CHKERRQ(ierr); 2269 ierr = VecDestroy(&b->y);CHKERRQ(ierr); 2270 ierr = VecDestroy(&b->counter);CHKERRQ(ierr); 2271 ierr = ISDestroy(&b->getsub_ris);CHKERRQ(ierr); 2272 ierr = ISDestroy(&b->getsub_cis);CHKERRQ(ierr); 2273 if (b->sf != b->csf) { 2274 ierr = PetscSFDestroy(&b->csf);CHKERRQ(ierr); 2275 ierr = PetscFree2(b->csf_rootdata,b->csf_leafdata);CHKERRQ(ierr); 2276 } else b->csf = NULL; 2277 ierr = PetscSFDestroy(&b->sf);CHKERRQ(ierr); 2278 ierr = PetscFree2(b->sf_rootdata,b->sf_leafdata);CHKERRQ(ierr); 2279 ierr = PetscFree(A->data);CHKERRQ(ierr); 2280 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 2281 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",NULL);CHKERRQ(ierr); 2282 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);CHKERRQ(ierr); 2283 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);CHKERRQ(ierr); 2284 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);CHKERRQ(ierr); 2285 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);CHKERRQ(ierr); 2286 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL);CHKERRQ(ierr); 2287 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",NULL);CHKERRQ(ierr); 2288 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL);CHKERRQ(ierr); 2289 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL);CHKERRQ(ierr); 2290 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL);CHKERRQ(ierr); 2291 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL);CHKERRQ(ierr); 2292 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL);CHKERRQ(ierr); 2293 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL);CHKERRQ(ierr); 2294 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL);CHKERRQ(ierr); 2295 PetscFunctionReturn(0); 2296 } 2297 2298 static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y) 2299 { 2300 PetscErrorCode ierr; 2301 Mat_IS *is = (Mat_IS*)A->data; 2302 PetscScalar zero = 0.0; 2303 2304 PetscFunctionBegin; 2305 /* scatter the global vector x into the local work vector */ 2306 ierr = VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2307 ierr = VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2308 2309 /* multiply the local matrix */ 2310 ierr = MatMult(is->A,is->x,is->y);CHKERRQ(ierr); 2311 2312 /* scatter product back into global memory */ 2313 ierr = VecSet(y,zero);CHKERRQ(ierr); 2314 ierr = VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2315 ierr = VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2316 PetscFunctionReturn(0); 2317 } 2318 2319 static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 2320 { 2321 Vec temp_vec; 2322 PetscErrorCode ierr; 2323 2324 PetscFunctionBegin; /* v3 = v2 + A * v1.*/ 2325 if (v3 != v2) { 2326 ierr = MatMult(A,v1,v3);CHKERRQ(ierr); 2327 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 2328 } else { 2329 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 2330 ierr = MatMult(A,v1,temp_vec);CHKERRQ(ierr); 2331 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 2332 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 2333 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 2334 } 2335 PetscFunctionReturn(0); 2336 } 2337 2338 static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x) 2339 { 2340 Mat_IS *is = (Mat_IS*)A->data; 2341 PetscErrorCode ierr; 2342 2343 PetscFunctionBegin; 2344 /* scatter the global vector x into the local work vector */ 2345 ierr = VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2346 ierr = VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2347 2348 /* multiply the local matrix */ 2349 ierr = MatMultTranspose(is->A,is->y,is->x);CHKERRQ(ierr); 2350 2351 /* scatter product back into global vector */ 2352 ierr = VecSet(x,0);CHKERRQ(ierr); 2353 ierr = VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2354 ierr = VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2355 PetscFunctionReturn(0); 2356 } 2357 2358 static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 2359 { 2360 Vec temp_vec; 2361 PetscErrorCode ierr; 2362 2363 PetscFunctionBegin; /* v3 = v2 + A' * v1.*/ 2364 if (v3 != v2) { 2365 ierr = MatMultTranspose(A,v1,v3);CHKERRQ(ierr); 2366 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 2367 } else { 2368 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 2369 ierr = MatMultTranspose(A,v1,temp_vec);CHKERRQ(ierr); 2370 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 2371 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 2372 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 2373 } 2374 PetscFunctionReturn(0); 2375 } 2376 2377 static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer) 2378 { 2379 Mat_IS *a = (Mat_IS*)A->data; 2380 PetscErrorCode ierr; 2381 PetscViewer sviewer; 2382 PetscBool isascii,view = PETSC_TRUE; 2383 2384 PetscFunctionBegin; 2385 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 2386 if (isascii) { 2387 PetscViewerFormat format; 2388 2389 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 2390 if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE; 2391 } 2392 if (!view) PetscFunctionReturn(0); 2393 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 2394 ierr = MatView(a->A,sviewer);CHKERRQ(ierr); 2395 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 2396 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2397 PetscFunctionReturn(0); 2398 } 2399 2400 static PetscErrorCode MatInvertBlockDiagonal_IS(Mat mat,const PetscScalar **values) 2401 { 2402 Mat_IS *is = (Mat_IS*)mat->data; 2403 MPI_Datatype nodeType; 2404 const PetscScalar *lv; 2405 PetscInt bs; 2406 PetscErrorCode ierr; 2407 2408 PetscFunctionBegin; 2409 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 2410 ierr = MatSetBlockSize(is->A,bs);CHKERRQ(ierr); 2411 ierr = MatInvertBlockDiagonal(is->A,&lv);CHKERRQ(ierr); 2412 if (!is->bdiag) { 2413 ierr = PetscMalloc1(bs*mat->rmap->n,&is->bdiag);CHKERRQ(ierr); 2414 } 2415 ierr = MPI_Type_contiguous(bs,MPIU_SCALAR,&nodeType);CHKERRQ(ierr); 2416 ierr = MPI_Type_commit(&nodeType);CHKERRQ(ierr); 2417 ierr = PetscSFReduceBegin(is->sf,nodeType,lv,is->bdiag,MPIU_REPLACE);CHKERRQ(ierr); 2418 ierr = PetscSFReduceEnd(is->sf,nodeType,lv,is->bdiag,MPIU_REPLACE);CHKERRQ(ierr); 2419 ierr = MPI_Type_free(&nodeType);CHKERRQ(ierr); 2420 if (values) *values = is->bdiag; 2421 PetscFunctionReturn(0); 2422 } 2423 2424 static PetscErrorCode MatISSetUpScatters_Private(Mat A) 2425 { 2426 Vec cglobal,rglobal; 2427 IS from; 2428 Mat_IS *is = (Mat_IS*)A->data; 2429 PetscScalar sum; 2430 const PetscInt *garray; 2431 PetscInt nr,rbs,nc,cbs; 2432 PetscBool iscuda; 2433 PetscErrorCode ierr; 2434 2435 PetscFunctionBegin; 2436 ierr = ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nr);CHKERRQ(ierr); 2437 ierr = ISLocalToGlobalMappingGetBlockSize(A->rmap->mapping,&rbs);CHKERRQ(ierr); 2438 ierr = ISLocalToGlobalMappingGetSize(A->cmap->mapping,&nc);CHKERRQ(ierr); 2439 ierr = ISLocalToGlobalMappingGetBlockSize(A->cmap->mapping,&cbs);CHKERRQ(ierr); 2440 ierr = VecDestroy(&is->x);CHKERRQ(ierr); 2441 ierr = VecDestroy(&is->y);CHKERRQ(ierr); 2442 ierr = VecDestroy(&is->counter);CHKERRQ(ierr); 2443 ierr = VecScatterDestroy(&is->rctx);CHKERRQ(ierr); 2444 ierr = VecScatterDestroy(&is->cctx);CHKERRQ(ierr); 2445 ierr = MatCreateVecs(is->A,&is->x,&is->y);CHKERRQ(ierr); 2446 ierr = PetscObjectTypeCompare((PetscObject)is->y,VECSEQCUDA,&iscuda);CHKERRQ(ierr); 2447 if (iscuda) { 2448 ierr = PetscFree(A->defaultvectype);CHKERRQ(ierr); 2449 ierr = PetscStrallocpy(VECCUDA,&A->defaultvectype);CHKERRQ(ierr); 2450 } 2451 ierr = MatCreateVecs(A,&cglobal,&rglobal);CHKERRQ(ierr); 2452 ierr = ISLocalToGlobalMappingGetBlockIndices(A->rmap->mapping,&garray);CHKERRQ(ierr); 2453 ierr = ISCreateBlock(PetscObjectComm((PetscObject)A),rbs,nr/rbs,garray,PETSC_USE_POINTER,&from);CHKERRQ(ierr); 2454 ierr = VecScatterCreateWithData(rglobal,from,is->y,NULL,&is->rctx);CHKERRQ(ierr); 2455 ierr = ISLocalToGlobalMappingRestoreBlockIndices(A->rmap->mapping,&garray);CHKERRQ(ierr); 2456 ierr = ISDestroy(&from);CHKERRQ(ierr); 2457 if (A->rmap->mapping != A->cmap->mapping) { 2458 ierr = ISLocalToGlobalMappingGetBlockIndices(A->cmap->mapping,&garray);CHKERRQ(ierr); 2459 ierr = ISCreateBlock(PetscObjectComm((PetscObject)A),cbs,nc/cbs,garray,PETSC_USE_POINTER,&from);CHKERRQ(ierr); 2460 ierr = VecScatterCreateWithData(cglobal,from,is->x,NULL,&is->cctx);CHKERRQ(ierr); 2461 ierr = ISLocalToGlobalMappingRestoreBlockIndices(A->cmap->mapping,&garray);CHKERRQ(ierr); 2462 ierr = ISDestroy(&from);CHKERRQ(ierr); 2463 } else { 2464 ierr = PetscObjectReference((PetscObject)is->rctx);CHKERRQ(ierr); 2465 is->cctx = is->rctx; 2466 } 2467 ierr = VecDestroy(&cglobal);CHKERRQ(ierr); 2468 2469 /* interface counter vector (local) */ 2470 ierr = VecDuplicate(is->y,&is->counter);CHKERRQ(ierr); 2471 ierr = VecSet(is->y,1.);CHKERRQ(ierr); 2472 ierr = VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2473 ierr = VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2474 ierr = VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2475 ierr = VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2476 2477 /* special functions for block-diagonal matrices */ 2478 ierr = VecSum(rglobal,&sum);CHKERRQ(ierr); 2479 if ((PetscInt)(PetscRealPart(sum)) == A->rmap->N && A->rmap->N == A->cmap->N && A->rmap->mapping == A->cmap->mapping) { 2480 A->ops->invertblockdiagonal = MatInvertBlockDiagonal_IS; 2481 } else { 2482 A->ops->invertblockdiagonal = NULL; 2483 } 2484 ierr = VecDestroy(&rglobal);CHKERRQ(ierr); 2485 2486 /* setup SF for general purpose shared indices based communications */ 2487 ierr = MatISSetUpSF_IS(A);CHKERRQ(ierr); 2488 PetscFunctionReturn(0); 2489 } 2490 2491 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping) 2492 { 2493 PetscErrorCode ierr; 2494 PetscInt nr,rbs,nc,cbs; 2495 Mat_IS *is = (Mat_IS*)A->data; 2496 2497 PetscFunctionBegin; 2498 PetscCheckSameComm(A,1,rmapping,2); 2499 PetscCheckSameComm(A,1,cmapping,3); 2500 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 2501 if (is->csf != is->sf) { 2502 ierr = PetscSFDestroy(&is->csf);CHKERRQ(ierr); 2503 ierr = PetscFree2(is->csf_rootdata,is->csf_leafdata);CHKERRQ(ierr); 2504 } else is->csf = NULL; 2505 ierr = PetscSFDestroy(&is->sf);CHKERRQ(ierr); 2506 ierr = PetscFree2(is->sf_rootdata,is->sf_leafdata);CHKERRQ(ierr); 2507 ierr = PetscFree(is->bdiag);CHKERRQ(ierr); 2508 2509 /* Setup Layout and set local to global maps */ 2510 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 2511 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 2512 ierr = ISLocalToGlobalMappingGetSize(rmapping,&nr);CHKERRQ(ierr); 2513 ierr = ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);CHKERRQ(ierr); 2514 ierr = ISLocalToGlobalMappingGetSize(cmapping,&nc);CHKERRQ(ierr); 2515 ierr = ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);CHKERRQ(ierr); 2516 /* check if the two mappings are actually the same for square matrices (DOLFIN passes 2 different objects) */ 2517 if (rmapping != cmapping && A->rmap->N == A->cmap->N) { 2518 PetscBool same,gsame; 2519 2520 same = PETSC_FALSE; 2521 if (nr == nc && cbs == rbs) { 2522 const PetscInt *idxs1,*idxs2; 2523 2524 ierr = ISLocalToGlobalMappingGetBlockIndices(rmapping,&idxs1);CHKERRQ(ierr); 2525 ierr = ISLocalToGlobalMappingGetBlockIndices(cmapping,&idxs2);CHKERRQ(ierr); 2526 ierr = PetscMemcmp(idxs1,idxs2,(nr/rbs)*sizeof(PetscInt),&same);CHKERRQ(ierr); 2527 ierr = ISLocalToGlobalMappingRestoreBlockIndices(rmapping,&idxs1);CHKERRQ(ierr); 2528 ierr = ISLocalToGlobalMappingRestoreBlockIndices(cmapping,&idxs2);CHKERRQ(ierr); 2529 } 2530 ierr = MPIU_Allreduce(&same,&gsame,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2531 if (gsame) cmapping = rmapping; 2532 } 2533 ierr = PetscLayoutSetBlockSize(A->rmap,rbs);CHKERRQ(ierr); 2534 ierr = PetscLayoutSetBlockSize(A->cmap,cbs);CHKERRQ(ierr); 2535 ierr = PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);CHKERRQ(ierr); 2536 ierr = PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);CHKERRQ(ierr); 2537 2538 /* Create the local matrix A */ 2539 ierr = MatCreate(PETSC_COMM_SELF,&is->A);CHKERRQ(ierr); 2540 ierr = MatSetType(is->A,is->lmattype);CHKERRQ(ierr); 2541 ierr = MatSetSizes(is->A,nr,nc,nr,nc);CHKERRQ(ierr); 2542 ierr = MatSetBlockSizes(is->A,rbs,cbs);CHKERRQ(ierr); 2543 ierr = MatSetOptionsPrefix(is->A,"is_");CHKERRQ(ierr); 2544 ierr = MatAppendOptionsPrefix(is->A,((PetscObject)A)->prefix);CHKERRQ(ierr); 2545 ierr = PetscLayoutSetUp(is->A->rmap);CHKERRQ(ierr); 2546 ierr = PetscLayoutSetUp(is->A->cmap);CHKERRQ(ierr); 2547 2548 if (!is->islocalref) { /* setup scatters and local vectors for MatMult */ 2549 ierr = MatISSetUpScatters_Private(A);CHKERRQ(ierr); 2550 } 2551 ierr = MatSetUp(A);CHKERRQ(ierr); 2552 PetscFunctionReturn(0); 2553 } 2554 2555 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2556 { 2557 Mat_IS *is = (Mat_IS*)mat->data; 2558 PetscErrorCode ierr; 2559 #if defined(PETSC_USE_DEBUG) 2560 PetscInt i,zm,zn; 2561 #endif 2562 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2563 2564 PetscFunctionBegin; 2565 #if defined(PETSC_USE_DEBUG) 2566 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); 2567 /* count negative indices */ 2568 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 2569 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 2570 #endif 2571 ierr = ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 2572 ierr = ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 2573 #if defined(PETSC_USE_DEBUG) 2574 /* count negative indices (should be the same as before) */ 2575 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 2576 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 2577 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"); 2578 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"); 2579 #endif 2580 ierr = MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 2581 PetscFunctionReturn(0); 2582 } 2583 2584 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2585 { 2586 Mat_IS *is = (Mat_IS*)mat->data; 2587 PetscErrorCode ierr; 2588 #if defined(PETSC_USE_DEBUG) 2589 PetscInt i,zm,zn; 2590 #endif 2591 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2592 2593 PetscFunctionBegin; 2594 #if defined(PETSC_USE_DEBUG) 2595 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); 2596 /* count negative indices */ 2597 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 2598 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 2599 #endif 2600 ierr = ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 2601 ierr = ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 2602 #if defined(PETSC_USE_DEBUG) 2603 /* count negative indices (should be the same as before) */ 2604 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 2605 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 2606 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"); 2607 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"); 2608 #endif 2609 ierr = MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 2610 PetscFunctionReturn(0); 2611 } 2612 2613 static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2614 { 2615 PetscErrorCode ierr; 2616 Mat_IS *is = (Mat_IS*)A->data; 2617 2618 PetscFunctionBegin; 2619 if (is->A->rmap->mapping) { 2620 ierr = MatSetValuesLocal(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2621 } else { 2622 ierr = MatSetValues(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2623 } 2624 PetscFunctionReturn(0); 2625 } 2626 2627 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2628 { 2629 PetscErrorCode ierr; 2630 Mat_IS *is = (Mat_IS*)A->data; 2631 2632 PetscFunctionBegin; 2633 if (is->A->rmap->mapping) { 2634 #if defined(PETSC_USE_DEBUG) 2635 PetscInt ibs,bs; 2636 2637 ierr = ISLocalToGlobalMappingGetBlockSize(is->A->rmap->mapping,&ibs);CHKERRQ(ierr); 2638 ierr = MatGetBlockSize(is->A,&bs);CHKERRQ(ierr); 2639 if (ibs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Different block sizes! local mat %D, local l2g map %D",bs,ibs); 2640 #endif 2641 ierr = MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2642 } else { 2643 ierr = MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2644 } 2645 PetscFunctionReturn(0); 2646 } 2647 2648 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns) 2649 { 2650 Mat_IS *is = (Mat_IS*)A->data; 2651 PetscErrorCode ierr; 2652 2653 PetscFunctionBegin; 2654 if (!n) { 2655 is->pure_neumann = PETSC_TRUE; 2656 } else { 2657 PetscInt i; 2658 is->pure_neumann = PETSC_FALSE; 2659 2660 if (columns) { 2661 ierr = MatZeroRowsColumns(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 2662 } else { 2663 ierr = MatZeroRows(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 2664 } 2665 if (diag != 0.) { 2666 const PetscScalar *array; 2667 ierr = VecGetArrayRead(is->counter,&array);CHKERRQ(ierr); 2668 for (i=0; i<n; i++) { 2669 ierr = MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);CHKERRQ(ierr); 2670 } 2671 ierr = VecRestoreArrayRead(is->counter,&array);CHKERRQ(ierr); 2672 } 2673 ierr = MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2674 ierr = MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2675 } 2676 PetscFunctionReturn(0); 2677 } 2678 2679 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns) 2680 { 2681 Mat_IS *matis = (Mat_IS*)A->data; 2682 PetscInt nr,nl,len,i; 2683 PetscInt *lrows; 2684 PetscErrorCode ierr; 2685 2686 PetscFunctionBegin; 2687 #if defined(PETSC_USE_DEBUG) 2688 if (columns || diag != 0. || (x && b)) { 2689 PetscBool cong; 2690 2691 ierr = PetscLayoutCompare(A->rmap,A->cmap,&cong);CHKERRQ(ierr); 2692 cong = (PetscBool)(cong && matis->sf == matis->csf); 2693 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"); 2694 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"); 2695 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"); 2696 } 2697 #endif 2698 /* get locally owned rows */ 2699 ierr = PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);CHKERRQ(ierr); 2700 /* fix right hand side if needed */ 2701 if (x && b) { 2702 const PetscScalar *xx; 2703 PetscScalar *bb; 2704 2705 ierr = VecGetArrayRead(x, &xx);CHKERRQ(ierr); 2706 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 2707 for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]]; 2708 ierr = VecRestoreArrayRead(x, &xx);CHKERRQ(ierr); 2709 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 2710 } 2711 /* get rows associated to the local matrices */ 2712 ierr = MatGetSize(matis->A,&nl,NULL);CHKERRQ(ierr); 2713 ierr = PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 2714 ierr = PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2715 for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1; 2716 ierr = PetscFree(lrows);CHKERRQ(ierr); 2717 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 2718 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 2719 ierr = PetscMalloc1(nl,&lrows);CHKERRQ(ierr); 2720 for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i; 2721 ierr = MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);CHKERRQ(ierr); 2722 ierr = PetscFree(lrows);CHKERRQ(ierr); 2723 PetscFunctionReturn(0); 2724 } 2725 2726 static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2727 { 2728 PetscErrorCode ierr; 2729 2730 PetscFunctionBegin; 2731 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);CHKERRQ(ierr); 2732 PetscFunctionReturn(0); 2733 } 2734 2735 static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2736 { 2737 PetscErrorCode ierr; 2738 2739 PetscFunctionBegin; 2740 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);CHKERRQ(ierr); 2741 PetscFunctionReturn(0); 2742 } 2743 2744 static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type) 2745 { 2746 Mat_IS *is = (Mat_IS*)A->data; 2747 PetscErrorCode ierr; 2748 2749 PetscFunctionBegin; 2750 ierr = MatAssemblyBegin(is->A,type);CHKERRQ(ierr); 2751 PetscFunctionReturn(0); 2752 } 2753 2754 static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type) 2755 { 2756 Mat_IS *is = (Mat_IS*)A->data; 2757 PetscErrorCode ierr; 2758 2759 PetscFunctionBegin; 2760 ierr = MatAssemblyEnd(is->A,type);CHKERRQ(ierr); 2761 /* fix for local empty rows/cols */ 2762 if (is->locempty && type == MAT_FINAL_ASSEMBLY) { 2763 Mat newlA; 2764 ISLocalToGlobalMapping rl2g,cl2g; 2765 IS nzr,nzc; 2766 PetscInt nr,nc,nnzr,nnzc; 2767 PetscBool lnewl2g,newl2g; 2768 2769 ierr = MatGetSize(is->A,&nr,&nc);CHKERRQ(ierr); 2770 ierr = MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_FALSE,PETSC_SMALL,&nzr);CHKERRQ(ierr); 2771 if (!nzr) { 2772 ierr = ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&nzr);CHKERRQ(ierr); 2773 } 2774 ierr = MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_TRUE,PETSC_SMALL,&nzc);CHKERRQ(ierr); 2775 if (!nzc) { 2776 ierr = ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&nzc);CHKERRQ(ierr); 2777 } 2778 ierr = ISGetSize(nzr,&nnzr);CHKERRQ(ierr); 2779 ierr = ISGetSize(nzc,&nnzc);CHKERRQ(ierr); 2780 if (nnzr != nr || nnzc != nc) { 2781 ISLocalToGlobalMapping l2g; 2782 IS is1,is2; 2783 2784 /* need new global l2g map */ 2785 lnewl2g = PETSC_TRUE; 2786 ierr = MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2787 2788 /* extract valid submatrix */ 2789 ierr = MatCreateSubMatrix(is->A,nzr,nzc,MAT_INITIAL_MATRIX,&newlA);CHKERRQ(ierr); 2790 2791 /* attach local l2g maps for successive calls of MatSetValues on the local matrix */ 2792 ierr = ISLocalToGlobalMappingCreateIS(nzr,&l2g);CHKERRQ(ierr); 2793 ierr = ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&is1);CHKERRQ(ierr); 2794 ierr = ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);CHKERRQ(ierr); 2795 ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr); 2796 if (is->A->rmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */ 2797 const PetscInt *idxs1,*idxs2; 2798 PetscInt j,i,nl,*tidxs; 2799 2800 ierr = ISLocalToGlobalMappingGetSize(is->A->rmap->mapping,&nl);CHKERRQ(ierr); 2801 ierr = ISLocalToGlobalMappingGetIndices(is->A->rmap->mapping,&idxs1);CHKERRQ(ierr); 2802 ierr = PetscMalloc1(nl,&tidxs);CHKERRQ(ierr); 2803 ierr = ISGetIndices(is2,&idxs2);CHKERRQ(ierr); 2804 for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++]; 2805 if (j != nr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nr); 2806 ierr = ISLocalToGlobalMappingRestoreIndices(is->A->rmap->mapping,&idxs1);CHKERRQ(ierr); 2807 ierr = ISRestoreIndices(is2,&idxs2);CHKERRQ(ierr); 2808 ierr = ISDestroy(&is2);CHKERRQ(ierr); 2809 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is->A->rmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);CHKERRQ(ierr); 2810 } 2811 ierr = ISLocalToGlobalMappingCreateIS(is2,&rl2g);CHKERRQ(ierr); 2812 ierr = ISDestroy(&is1);CHKERRQ(ierr); 2813 ierr = ISDestroy(&is2);CHKERRQ(ierr); 2814 2815 ierr = ISLocalToGlobalMappingCreateIS(nzc,&l2g);CHKERRQ(ierr); 2816 ierr = ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&is1);CHKERRQ(ierr); 2817 ierr = ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);CHKERRQ(ierr); 2818 ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr); 2819 if (is->A->cmap->mapping) { /* the matrix has a local-to-local map already attached (probably DMDA generated) */ 2820 const PetscInt *idxs1,*idxs2; 2821 PetscInt j,i,nl,*tidxs; 2822 2823 ierr = ISLocalToGlobalMappingGetSize(is->A->cmap->mapping,&nl);CHKERRQ(ierr); 2824 ierr = ISLocalToGlobalMappingGetIndices(is->A->cmap->mapping,&idxs1);CHKERRQ(ierr); 2825 ierr = PetscMalloc1(nl,&tidxs);CHKERRQ(ierr); 2826 ierr = ISGetIndices(is2,&idxs2);CHKERRQ(ierr); 2827 for (i=0,j=0;i<nl;i++) tidxs[i] = idxs1[i] < 0 ? -1 : idxs2[j++]; 2828 if (j != nc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected count %D != %D",j,nc); 2829 ierr = ISLocalToGlobalMappingRestoreIndices(is->A->cmap->mapping,&idxs1);CHKERRQ(ierr); 2830 ierr = ISRestoreIndices(is2,&idxs2);CHKERRQ(ierr); 2831 ierr = ISDestroy(&is2);CHKERRQ(ierr); 2832 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is->A->cmap->mapping),nl,tidxs,PETSC_OWN_POINTER,&is2);CHKERRQ(ierr); 2833 } 2834 ierr = ISLocalToGlobalMappingCreateIS(is2,&cl2g);CHKERRQ(ierr); 2835 ierr = ISDestroy(&is1);CHKERRQ(ierr); 2836 ierr = ISDestroy(&is2);CHKERRQ(ierr); 2837 2838 ierr = MatSetLocalToGlobalMapping(newlA,rl2g,cl2g);CHKERRQ(ierr); 2839 2840 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 2841 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 2842 } else { /* local matrix fully populated */ 2843 lnewl2g = PETSC_FALSE; 2844 ierr = MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2845 ierr = PetscObjectReference((PetscObject)is->A);CHKERRQ(ierr); 2846 newlA = is->A; 2847 } 2848 /* attach new global l2g map if needed */ 2849 if (newl2g) { 2850 IS gnzr,gnzc; 2851 const PetscInt *grid,*gcid; 2852 2853 ierr = ISLocalToGlobalMappingApplyIS(A->rmap->mapping,nzr,&gnzr);CHKERRQ(ierr); 2854 ierr = ISLocalToGlobalMappingApplyIS(A->cmap->mapping,nzc,&gnzc);CHKERRQ(ierr); 2855 ierr = ISGetIndices(gnzr,&grid);CHKERRQ(ierr); 2856 ierr = ISGetIndices(gnzc,&gcid);CHKERRQ(ierr); 2857 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzr,grid,PETSC_COPY_VALUES,&rl2g);CHKERRQ(ierr); 2858 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzc,gcid,PETSC_COPY_VALUES,&cl2g);CHKERRQ(ierr); 2859 ierr = ISRestoreIndices(gnzr,&grid);CHKERRQ(ierr); 2860 ierr = ISRestoreIndices(gnzc,&gcid);CHKERRQ(ierr); 2861 ierr = ISDestroy(&gnzr);CHKERRQ(ierr); 2862 ierr = ISDestroy(&gnzc);CHKERRQ(ierr); 2863 ierr = MatSetLocalToGlobalMapping(A,rl2g,cl2g);CHKERRQ(ierr); 2864 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 2865 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 2866 } 2867 ierr = MatISSetLocalMat(A,newlA);CHKERRQ(ierr); 2868 ierr = MatDestroy(&newlA);CHKERRQ(ierr); 2869 ierr = ISDestroy(&nzr);CHKERRQ(ierr); 2870 ierr = ISDestroy(&nzc);CHKERRQ(ierr); 2871 is->locempty = PETSC_FALSE; 2872 } 2873 PetscFunctionReturn(0); 2874 } 2875 2876 static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local) 2877 { 2878 Mat_IS *is = (Mat_IS*)mat->data; 2879 2880 PetscFunctionBegin; 2881 *local = is->A; 2882 PetscFunctionReturn(0); 2883 } 2884 2885 static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local) 2886 { 2887 PetscFunctionBegin; 2888 *local = NULL; 2889 PetscFunctionReturn(0); 2890 } 2891 2892 /*@ 2893 MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix. 2894 2895 Input Parameter: 2896 . mat - the matrix 2897 2898 Output Parameter: 2899 . local - the local matrix 2900 2901 Level: advanced 2902 2903 Notes: 2904 This can be called if you have precomputed the nonzero structure of the 2905 matrix and want to provide it to the inner matrix object to improve the performance 2906 of the MatSetValues() operation. 2907 2908 Call MatISRestoreLocalMat() when finished with the local matrix. 2909 2910 .seealso: MATIS 2911 @*/ 2912 PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local) 2913 { 2914 PetscErrorCode ierr; 2915 2916 PetscFunctionBegin; 2917 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2918 PetscValidPointer(local,2); 2919 ierr = PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 2920 PetscFunctionReturn(0); 2921 } 2922 2923 /*@ 2924 MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat() 2925 2926 Input Parameter: 2927 . mat - the matrix 2928 2929 Output Parameter: 2930 . local - the local matrix 2931 2932 Level: advanced 2933 2934 .seealso: MATIS 2935 @*/ 2936 PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local) 2937 { 2938 PetscErrorCode ierr; 2939 2940 PetscFunctionBegin; 2941 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2942 PetscValidPointer(local,2); 2943 ierr = PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 2944 PetscFunctionReturn(0); 2945 } 2946 2947 static PetscErrorCode MatISSetLocalMatType_IS(Mat mat,MatType mtype) 2948 { 2949 Mat_IS *is = (Mat_IS*)mat->data; 2950 PetscErrorCode ierr; 2951 2952 PetscFunctionBegin; 2953 if (is->A) { 2954 ierr = MatSetType(is->A,mtype);CHKERRQ(ierr); 2955 } 2956 ierr = PetscFree(is->lmattype);CHKERRQ(ierr); 2957 ierr = PetscStrallocpy(mtype,&is->lmattype);CHKERRQ(ierr); 2958 PetscFunctionReturn(0); 2959 } 2960 2961 /*@ 2962 MatISSetLocalMatType - Specifies the type of local matrix 2963 2964 Input Parameter: 2965 . mat - the matrix 2966 . mtype - the local matrix type 2967 2968 Output Parameter: 2969 2970 Level: advanced 2971 2972 .seealso: MATIS, MatSetType(), MatType 2973 @*/ 2974 PetscErrorCode MatISSetLocalMatType(Mat mat,MatType mtype) 2975 { 2976 PetscErrorCode ierr; 2977 2978 PetscFunctionBegin; 2979 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2980 ierr = PetscUseMethod(mat,"MatISSetLocalMatType_C",(Mat,MatType),(mat,mtype));CHKERRQ(ierr); 2981 PetscFunctionReturn(0); 2982 } 2983 2984 static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local) 2985 { 2986 Mat_IS *is = (Mat_IS*)mat->data; 2987 PetscInt nrows,ncols,orows,ocols; 2988 PetscErrorCode ierr; 2989 MatType mtype,otype; 2990 PetscBool sametype = PETSC_TRUE; 2991 2992 PetscFunctionBegin; 2993 if (is->A) { 2994 ierr = MatGetSize(is->A,&orows,&ocols);CHKERRQ(ierr); 2995 ierr = MatGetSize(local,&nrows,&ncols);CHKERRQ(ierr); 2996 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); 2997 ierr = MatGetType(local,&mtype);CHKERRQ(ierr); 2998 ierr = MatGetType(is->A,&otype);CHKERRQ(ierr); 2999 ierr = PetscStrcmp(mtype,otype,&sametype);CHKERRQ(ierr); 3000 } 3001 ierr = PetscObjectReference((PetscObject)local);CHKERRQ(ierr); 3002 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 3003 is->A = local; 3004 ierr = MatGetType(is->A,&mtype);CHKERRQ(ierr); 3005 ierr = MatISSetLocalMatType(mat,mtype);CHKERRQ(ierr); 3006 if (!sametype && !is->islocalref) { 3007 ierr = MatISSetUpScatters_Private(mat);CHKERRQ(ierr); 3008 } 3009 PetscFunctionReturn(0); 3010 } 3011 3012 /*@ 3013 MatISSetLocalMat - Replace the local matrix stored inside a MATIS object. 3014 3015 Collective on Mat 3016 3017 Input Parameter: 3018 . mat - the matrix 3019 . local - the local matrix 3020 3021 Output Parameter: 3022 3023 Level: advanced 3024 3025 Notes: 3026 This can be called if you have precomputed the local matrix and 3027 want to provide it to the matrix object MATIS. 3028 3029 .seealso: MATIS 3030 @*/ 3031 PetscErrorCode MatISSetLocalMat(Mat mat,Mat local) 3032 { 3033 PetscErrorCode ierr; 3034 3035 PetscFunctionBegin; 3036 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3037 PetscValidHeaderSpecific(local,MAT_CLASSID,2); 3038 ierr = PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));CHKERRQ(ierr); 3039 PetscFunctionReturn(0); 3040 } 3041 3042 static PetscErrorCode MatZeroEntries_IS(Mat A) 3043 { 3044 Mat_IS *a = (Mat_IS*)A->data; 3045 PetscErrorCode ierr; 3046 3047 PetscFunctionBegin; 3048 ierr = MatZeroEntries(a->A);CHKERRQ(ierr); 3049 PetscFunctionReturn(0); 3050 } 3051 3052 static PetscErrorCode MatScale_IS(Mat A,PetscScalar a) 3053 { 3054 Mat_IS *is = (Mat_IS*)A->data; 3055 PetscErrorCode ierr; 3056 3057 PetscFunctionBegin; 3058 ierr = MatScale(is->A,a);CHKERRQ(ierr); 3059 PetscFunctionReturn(0); 3060 } 3061 3062 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v) 3063 { 3064 Mat_IS *is = (Mat_IS*)A->data; 3065 PetscErrorCode ierr; 3066 3067 PetscFunctionBegin; 3068 /* get diagonal of the local matrix */ 3069 ierr = MatGetDiagonal(is->A,is->y);CHKERRQ(ierr); 3070 3071 /* scatter diagonal back into global vector */ 3072 ierr = VecSet(v,0);CHKERRQ(ierr); 3073 ierr = VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 3074 ierr = VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 3075 PetscFunctionReturn(0); 3076 } 3077 3078 static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg) 3079 { 3080 Mat_IS *a = (Mat_IS*)A->data; 3081 PetscErrorCode ierr; 3082 3083 PetscFunctionBegin; 3084 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 3085 PetscFunctionReturn(0); 3086 } 3087 3088 static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str) 3089 { 3090 Mat_IS *y = (Mat_IS*)Y->data; 3091 Mat_IS *x; 3092 #if defined(PETSC_USE_DEBUG) 3093 PetscBool ismatis; 3094 #endif 3095 PetscErrorCode ierr; 3096 3097 PetscFunctionBegin; 3098 #if defined(PETSC_USE_DEBUG) 3099 ierr = PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);CHKERRQ(ierr); 3100 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS"); 3101 #endif 3102 x = (Mat_IS*)X->data; 3103 ierr = MatAXPY(y->A,a,x->A,str);CHKERRQ(ierr); 3104 PetscFunctionReturn(0); 3105 } 3106 3107 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat) 3108 { 3109 Mat lA; 3110 Mat_IS *matis; 3111 ISLocalToGlobalMapping rl2g,cl2g; 3112 IS is; 3113 const PetscInt *rg,*rl; 3114 PetscInt nrg; 3115 PetscInt N,M,nrl,i,*idxs; 3116 PetscErrorCode ierr; 3117 3118 PetscFunctionBegin; 3119 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 3120 ierr = ISGetLocalSize(row,&nrl);CHKERRQ(ierr); 3121 ierr = ISGetIndices(row,&rl);CHKERRQ(ierr); 3122 ierr = ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);CHKERRQ(ierr); 3123 #if defined(PETSC_USE_DEBUG) 3124 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); 3125 #endif 3126 ierr = PetscMalloc1(nrg,&idxs);CHKERRQ(ierr); 3127 /* map from [0,nrl) to row */ 3128 for (i=0;i<nrl;i++) idxs[i] = rl[i]; 3129 for (i=nrl;i<nrg;i++) idxs[i] = -1; 3130 ierr = ISRestoreIndices(row,&rl);CHKERRQ(ierr); 3131 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 3132 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 3133 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 3134 ierr = ISDestroy(&is);CHKERRQ(ierr); 3135 /* compute new l2g map for columns */ 3136 if (col != row || A->rmap->mapping != A->cmap->mapping) { 3137 const PetscInt *cg,*cl; 3138 PetscInt ncg; 3139 PetscInt ncl; 3140 3141 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 3142 ierr = ISGetLocalSize(col,&ncl);CHKERRQ(ierr); 3143 ierr = ISGetIndices(col,&cl);CHKERRQ(ierr); 3144 ierr = ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);CHKERRQ(ierr); 3145 #if defined(PETSC_USE_DEBUG) 3146 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); 3147 #endif 3148 ierr = PetscMalloc1(ncg,&idxs);CHKERRQ(ierr); 3149 /* map from [0,ncl) to col */ 3150 for (i=0;i<ncl;i++) idxs[i] = cl[i]; 3151 for (i=ncl;i<ncg;i++) idxs[i] = -1; 3152 ierr = ISRestoreIndices(col,&cl);CHKERRQ(ierr); 3153 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 3154 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 3155 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 3156 ierr = ISDestroy(&is);CHKERRQ(ierr); 3157 } else { 3158 ierr = PetscObjectReference((PetscObject)rl2g);CHKERRQ(ierr); 3159 cl2g = rl2g; 3160 } 3161 /* create the MATIS submatrix */ 3162 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 3163 ierr = MatCreate(PetscObjectComm((PetscObject)A),submat);CHKERRQ(ierr); 3164 ierr = MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3165 ierr = MatSetType(*submat,MATIS);CHKERRQ(ierr); 3166 matis = (Mat_IS*)((*submat)->data); 3167 matis->islocalref = PETSC_TRUE; 3168 ierr = MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);CHKERRQ(ierr); 3169 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr); 3170 ierr = MatISSetLocalMat(*submat,lA);CHKERRQ(ierr); 3171 ierr = MatSetUp(*submat);CHKERRQ(ierr); 3172 ierr = MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3173 ierr = MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3174 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 3175 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 3176 /* remove unsupported ops */ 3177 ierr = PetscMemzero((*submat)->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 3178 (*submat)->ops->destroy = MatDestroy_IS; 3179 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS; 3180 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS; 3181 (*submat)->ops->assemblybegin = MatAssemblyBegin_IS; 3182 (*submat)->ops->assemblyend = MatAssemblyEnd_IS; 3183 PetscFunctionReturn(0); 3184 } 3185 3186 static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A) 3187 { 3188 Mat_IS *a = (Mat_IS*)A->data; 3189 char type[256]; 3190 PetscBool flg; 3191 PetscErrorCode ierr; 3192 3193 PetscFunctionBegin; 3194 ierr = PetscOptionsHead(PetscOptionsObject,"MATIS options");CHKERRQ(ierr); 3195 ierr = PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns","MatISFixLocalEmpty",a->locempty,&a->locempty,NULL);CHKERRQ(ierr); 3196 ierr = PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL);CHKERRQ(ierr); 3197 ierr = PetscOptionsFList("-matis_localmat_type","Matrix type","MatISSetLocalMatType",MatList,a->lmattype,type,256,&flg);CHKERRQ(ierr); 3198 if (flg) { 3199 ierr = MatISSetLocalMatType(A,type);CHKERRQ(ierr); 3200 } 3201 if (a->A) { 3202 ierr = MatSetFromOptions(a->A);CHKERRQ(ierr); 3203 } 3204 ierr = PetscOptionsTail();CHKERRQ(ierr); 3205 PetscFunctionReturn(0); 3206 } 3207 3208 /*@ 3209 MatCreateIS - Creates a "process" unassembled matrix, assembled on each 3210 process but not across processes. 3211 3212 Input Parameters: 3213 + comm - MPI communicator that will share the matrix 3214 . bs - block size of the matrix 3215 . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products 3216 . rmap - local to global map for rows 3217 - cmap - local to global map for cols 3218 3219 Output Parameter: 3220 . A - the resulting matrix 3221 3222 Level: advanced 3223 3224 Notes: 3225 See MATIS for more details. 3226 m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors 3227 used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices. 3228 If either rmap or cmap are NULL, then the matrix is assumed to be square. 3229 3230 .seealso: MATIS, MatSetLocalToGlobalMapping() 3231 @*/ 3232 PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A) 3233 { 3234 PetscErrorCode ierr; 3235 3236 PetscFunctionBegin; 3237 if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings"); 3238 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3239 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 3240 if (bs > 0) { 3241 ierr = MatSetBlockSize(*A,bs);CHKERRQ(ierr); 3242 } 3243 ierr = MatSetType(*A,MATIS);CHKERRQ(ierr); 3244 if (rmap && cmap) { 3245 ierr = MatSetLocalToGlobalMapping(*A,rmap,cmap);CHKERRQ(ierr); 3246 } else if (!rmap) { 3247 ierr = MatSetLocalToGlobalMapping(*A,cmap,cmap);CHKERRQ(ierr); 3248 } else { 3249 ierr = MatSetLocalToGlobalMapping(*A,rmap,rmap);CHKERRQ(ierr); 3250 } 3251 PetscFunctionReturn(0); 3252 } 3253 3254 /*MC 3255 MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP). 3256 This stores the matrices in globally unassembled form. Each processor assembles only its local Neumann problem and the parallel matrix vector 3257 product is handled "implicitly". 3258 3259 Options Database Keys: 3260 + -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions() 3261 . -matis_fixempty - Fixes local matrices in case of empty local rows/columns. 3262 - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of MatPtAP(). 3263 3264 Notes: 3265 Options prefix for the inner matrix are given by -is_mat_xxx 3266 3267 You must call MatSetLocalToGlobalMapping() before using this matrix type. 3268 3269 You can do matrix preallocation on the local matrix after you obtain it with 3270 MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation() 3271 3272 Level: advanced 3273 3274 .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP 3275 3276 M*/ 3277 3278 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A) 3279 { 3280 PetscErrorCode ierr; 3281 Mat_IS *b; 3282 3283 PetscFunctionBegin; 3284 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 3285 ierr = PetscStrallocpy(MATAIJ,&b->lmattype);CHKERRQ(ierr); 3286 A->data = (void*)b; 3287 3288 /* matrix ops */ 3289 ierr = PetscMemzero(A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 3290 A->ops->mult = MatMult_IS; 3291 A->ops->multadd = MatMultAdd_IS; 3292 A->ops->multtranspose = MatMultTranspose_IS; 3293 A->ops->multtransposeadd = MatMultTransposeAdd_IS; 3294 A->ops->destroy = MatDestroy_IS; 3295 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS; 3296 A->ops->setvalues = MatSetValues_IS; 3297 A->ops->setvaluesblocked = MatSetValuesBlocked_IS; 3298 A->ops->setvalueslocal = MatSetValuesLocal_IS; 3299 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS; 3300 A->ops->zerorows = MatZeroRows_IS; 3301 A->ops->zerorowscolumns = MatZeroRowsColumns_IS; 3302 A->ops->assemblybegin = MatAssemblyBegin_IS; 3303 A->ops->assemblyend = MatAssemblyEnd_IS; 3304 A->ops->view = MatView_IS; 3305 A->ops->zeroentries = MatZeroEntries_IS; 3306 A->ops->scale = MatScale_IS; 3307 A->ops->getdiagonal = MatGetDiagonal_IS; 3308 A->ops->setoption = MatSetOption_IS; 3309 A->ops->ishermitian = MatIsHermitian_IS; 3310 A->ops->issymmetric = MatIsSymmetric_IS; 3311 A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS; 3312 A->ops->duplicate = MatDuplicate_IS; 3313 A->ops->missingdiagonal = MatMissingDiagonal_IS; 3314 A->ops->copy = MatCopy_IS; 3315 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 3316 A->ops->createsubmatrix = MatCreateSubMatrix_IS; 3317 A->ops->axpy = MatAXPY_IS; 3318 A->ops->diagonalset = MatDiagonalSet_IS; 3319 A->ops->shift = MatShift_IS; 3320 A->ops->transpose = MatTranspose_IS; 3321 A->ops->getinfo = MatGetInfo_IS; 3322 A->ops->diagonalscale = MatDiagonalScale_IS; 3323 A->ops->setfromoptions = MatSetFromOptions_IS; 3324 3325 /* special MATIS functions */ 3326 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMatType_C",MatISSetLocalMatType_IS);CHKERRQ(ierr); 3327 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);CHKERRQ(ierr); 3328 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);CHKERRQ(ierr); 3329 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);CHKERRQ(ierr); 3330 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3331 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);CHKERRQ(ierr); 3332 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS);CHKERRQ(ierr); 3333 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",MatISFixLocalEmpty_IS);CHKERRQ(ierr); 3334 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3335 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3336 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3337 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3338 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3339 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3340 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3341 ierr = PetscObjectChangeTypeName((PetscObject)A,MATIS);CHKERRQ(ierr); 3342 PetscFunctionReturn(0); 3343 } 3344