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