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