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