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