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