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