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