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 12 #define MATIS_MAX_ENTRIES_INSERTION 2048 13 14 static PetscErrorCode MatISComputeSF_Private(Mat); 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "MatConvert_Nest_IS" 18 PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat A,MatType type,MatReuse reuse,Mat *newmat) 19 { 20 Mat **nest,*snest,**rnest,lA,B; 21 IS *iscol,*isrow,*islrow,*islcol; 22 ISLocalToGlobalMapping rl2g,cl2g; 23 MPI_Comm comm; 24 PetscInt *lr,*lc,*lrb,*lcb,*l2gidxs; 25 PetscInt sti,stl,i,j,nr,nc,rbs,cbs; 26 PetscBool convert,lreuse; 27 PetscErrorCode ierr; 28 29 PetscFunctionBeginUser; 30 ierr = MatNestGetSubMats(A,&nr,&nc,&nest);CHKERRQ(ierr); 31 lreuse = PETSC_FALSE; 32 rnest = NULL; 33 if (reuse == MAT_REUSE_MATRIX) { 34 PetscBool ismatis,isnest; 35 36 ierr = PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);CHKERRQ(ierr); 37 if (!ismatis) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type); 38 ierr = MatISGetLocalMat(*newmat,&lA);CHKERRQ(ierr); 39 ierr = PetscObjectTypeCompare((PetscObject)lA,MATNEST,&isnest);CHKERRQ(ierr); 40 if (isnest) { 41 ierr = MatNestGetSubMats(lA,&i,&j,&rnest);CHKERRQ(ierr); 42 lreuse = (PetscBool)(i == nr && j == nc); 43 if (!lreuse) rnest = NULL; 44 } 45 } 46 ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 47 ierr = PetscCalloc4(nr,&lr,nc,&lc,nr,&lrb,nc,&lcb);CHKERRQ(ierr); 48 ierr = PetscCalloc5(nr,&isrow,nc,&iscol, 49 nr,&islrow,nc,&islcol, 50 nr*nc,&snest);CHKERRQ(ierr); 51 ierr = MatNestGetISs(A,isrow,iscol);CHKERRQ(ierr); 52 for (i=0;i<nr;i++) { 53 for (j=0;j<nc;j++) { 54 PetscBool ismatis; 55 PetscInt l1,l2,ij=i*nc+j; 56 57 /* Null matrix pointers are allowed in MATNEST */ 58 if (!nest[i][j]) continue; 59 60 /* Nested matrices should be of type MATIS */ 61 ierr = PetscObjectTypeCompare((PetscObject)nest[i][j],MATIS,&ismatis);CHKERRQ(ierr); 62 if (!ismatis) SETERRQ2(comm,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) is not of type MATIS",i,j); 63 64 /* Check compatibility of local sizes */ 65 ierr = MatISGetLocalMat(nest[i][j],&snest[ij]);CHKERRQ(ierr); 66 ierr = MatGetSize(snest[ij],&l1,&l2);CHKERRQ(ierr); 67 if (!l1 || !l2) continue; 68 if (lr[i] && l1 != lr[i]) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid local size %D != %D",i,j,lr[i],l1); 69 if (lc[j] && l2 != lc[j]) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid local size %D != %D",i,j,lc[j],l2); 70 lr[i] = l1; 71 lc[j] = l2; 72 ierr = MatGetBlockSizes(snest[ij],&l1,&l2);CHKERRQ(ierr); 73 if (lrb[i] && l1 != lrb[i]) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid block size %D != %D",i,j,lrb[i],l1); 74 if (lcb[j] && l2 != lcb[j]) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid block size %D != %D",i,j,lcb[j],l2); 75 lrb[i] = l1; 76 lcb[j] = l2; 77 78 /* check compatibilty for local matrix reusage */ 79 if (rnest && !rnest[i][j] != !snest[ij]) lreuse = PETSC_FALSE; 80 } 81 } 82 83 #if defined (PETSC_USE_DEBUG) 84 /* Check compatibility of l2g maps for rows */ 85 for (i=0;i<nr;i++) { 86 rl2g = NULL; 87 for (j=0;j<nc;j++) { 88 PetscInt n1,n2; 89 90 if (!nest[i][j]) continue; 91 ierr = MatGetLocalToGlobalMapping(nest[i][j],&cl2g,NULL);CHKERRQ(ierr); 92 ierr = ISLocalToGlobalMappingGetSize(cl2g,&n1);CHKERRQ(ierr); 93 if (!n1) continue; 94 if (!rl2g) { 95 rl2g = cl2g; 96 } else { 97 const PetscInt *idxs1,*idxs2; 98 PetscBool same; 99 100 ierr = ISLocalToGlobalMappingGetSize(rl2g,&n2);CHKERRQ(ierr); 101 if (n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid row l2gmap size %D != %D",i,j,n1,n2); 102 ierr = ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);CHKERRQ(ierr); 103 ierr = ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);CHKERRQ(ierr); 104 ierr = PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);CHKERRQ(ierr); 105 ierr = ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);CHKERRQ(ierr); 106 ierr = ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);CHKERRQ(ierr); 107 if (!same) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid row l2gmap",i,j); 108 } 109 } 110 } 111 /* Check compatibility of l2g maps for columns */ 112 for (i=0;i<nc;i++) { 113 rl2g = NULL; 114 for (j=0;j<nr;j++) { 115 PetscInt n1,n2; 116 117 if (!nest[j][i]) continue; 118 ierr = MatGetLocalToGlobalMapping(nest[j][i],NULL,&cl2g);CHKERRQ(ierr); 119 ierr = ISLocalToGlobalMappingGetSize(cl2g,&n1);CHKERRQ(ierr); 120 if (!n1) continue; 121 if (!rl2g) { 122 rl2g = cl2g; 123 } else { 124 const PetscInt *idxs1,*idxs2; 125 PetscBool same; 126 127 ierr = ISLocalToGlobalMappingGetSize(rl2g,&n2);CHKERRQ(ierr); 128 if (n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid column l2gmap size %D != %D",j,i,n1,n2); 129 ierr = ISLocalToGlobalMappingGetIndices(cl2g,&idxs1);CHKERRQ(ierr); 130 ierr = ISLocalToGlobalMappingGetIndices(rl2g,&idxs2);CHKERRQ(ierr); 131 ierr = PetscMemcmp(idxs1,idxs2,n1*sizeof(PetscInt),&same);CHKERRQ(ierr); 132 ierr = ISLocalToGlobalMappingRestoreIndices(cl2g,&idxs1);CHKERRQ(ierr); 133 ierr = ISLocalToGlobalMappingRestoreIndices(rl2g,&idxs2);CHKERRQ(ierr); 134 if (!same) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot convert from MATNEST to MATIS! Matrix block (%D,%D) has invalid column l2gmap",j,i); 135 } 136 } 137 } 138 #endif 139 140 B = NULL; 141 if (reuse != MAT_REUSE_MATRIX) { 142 /* Create l2g map for the rows of the new matrix and index sets for the local MATNEST */ 143 for (i=0,stl=0;i<nr;i++) stl += lr[i]; 144 ierr = PetscMalloc1(stl,&l2gidxs);CHKERRQ(ierr); 145 for (i=0,sti=0,stl=0;i<nr;i++) { 146 Mat usedmat; 147 Mat_IS *matis; 148 const PetscInt *idxs; 149 PetscInt n = lr[i]/lrb[i]; 150 151 /* local IS for local NEST */ 152 ierr = ISCreateStride(PETSC_COMM_SELF,n,sti,lrb[i],&islrow[i]);CHKERRQ(ierr); 153 sti += n; 154 155 /* l2gmap */ 156 j = 0; 157 usedmat = nest[i][j]; 158 while (!usedmat && j < nc) usedmat = nest[i][j++]; 159 matis = (Mat_IS*)(usedmat->data); 160 if (!matis->sf) { 161 ierr = MatISComputeSF_Private(usedmat);CHKERRQ(ierr); 162 } 163 ierr = ISGetIndices(isrow[i],&idxs);CHKERRQ(ierr); 164 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,idxs,l2gidxs+stl);CHKERRQ(ierr); 165 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,idxs,l2gidxs+stl);CHKERRQ(ierr); 166 ierr = ISRestoreIndices(isrow[i],&idxs);CHKERRQ(ierr); 167 stl += lr[i]; 168 } 169 ierr = ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&rl2g);CHKERRQ(ierr); 170 171 /* Create l2g map for columns of the new matrix and index sets for the local MATNEST */ 172 for (i=0,stl=0;i<nc;i++) stl += lc[i]; 173 ierr = PetscMalloc1(stl,&l2gidxs);CHKERRQ(ierr); 174 for (i=0,sti=0,stl=0;i<nc;i++) { 175 Mat usedmat; 176 Mat_IS *matis; 177 const PetscInt *idxs; 178 PetscInt n = lc[i]/lcb[i]; 179 180 /* local IS for local NEST */ 181 ierr = ISCreateStride(PETSC_COMM_SELF,n,sti,lcb[i],&islcol[i]);CHKERRQ(ierr); 182 sti += n; 183 184 /* l2gmap */ 185 j = 0; 186 usedmat = nest[j][i]; 187 while (!usedmat && j < nr) usedmat = nest[j++][i]; 188 matis = (Mat_IS*)(usedmat->data); 189 if (!matis->sf) { 190 ierr = MatISComputeSF_Private(usedmat);CHKERRQ(ierr); 191 } 192 ierr = ISGetIndices(iscol[i],&idxs);CHKERRQ(ierr); 193 ierr = PetscSFBcastBegin(matis->csf,MPIU_INT,idxs,l2gidxs+stl);CHKERRQ(ierr); 194 ierr = PetscSFBcastEnd(matis->csf,MPIU_INT,idxs,l2gidxs+stl);CHKERRQ(ierr); 195 ierr = ISRestoreIndices(iscol[i],&idxs);CHKERRQ(ierr); 196 stl += lc[i]; 197 } 198 ierr = ISLocalToGlobalMappingCreate(comm,1,stl,l2gidxs,PETSC_OWN_POINTER,&cl2g);CHKERRQ(ierr); 199 200 /* Create MATIS */ 201 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 202 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 203 ierr = MatGetBlockSizes(A,&rbs,&cbs);CHKERRQ(ierr); 204 ierr = MatSetBlockSizes(B,rbs,cbs);CHKERRQ(ierr); 205 ierr = MatSetType(B,MATIS);CHKERRQ(ierr); 206 ierr = MatSetLocalToGlobalMapping(B,rl2g,cl2g);CHKERRQ(ierr); 207 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 208 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 209 ierr = MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);CHKERRQ(ierr); 210 ierr = MatISSetLocalMat(B,lA);CHKERRQ(ierr); 211 ierr = MatDestroy(&lA);CHKERRQ(ierr); 212 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 213 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 214 if (reuse == MAT_INPLACE_MATRIX) { 215 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 216 } else { 217 *newmat = B; 218 } 219 } else { 220 if (lreuse) { 221 ierr = MatISGetLocalMat(*newmat,&lA);CHKERRQ(ierr); 222 for (i=0;i<nr;i++) { 223 for (j=0;j<nc;j++) { 224 if (snest[i*nc+j]) { 225 ierr = MatNestSetSubMat(lA,i,j,snest[i*nc+j]);CHKERRQ(ierr); 226 } 227 } 228 } 229 } else { 230 for (i=0,sti=0;i<nr;i++) { 231 PetscInt n = lr[i]/lrb[i]; 232 233 ierr = ISCreateStride(PETSC_COMM_SELF,n,sti,lrb[i],&islrow[i]);CHKERRQ(ierr); 234 sti += n; 235 } 236 for (i=0,sti=0;i<nc;i++) { 237 PetscInt n = lc[i]/lcb[i]; 238 239 ierr = ISCreateStride(PETSC_COMM_SELF,n,sti,lcb[i],&islcol[i]);CHKERRQ(ierr); 240 sti += n; 241 } 242 ierr = MatCreateNest(PETSC_COMM_SELF,nr,islrow,nc,islcol,snest,&lA);CHKERRQ(ierr); 243 ierr = MatISSetLocalMat(*newmat,lA);CHKERRQ(ierr); 244 ierr = MatDestroy(&lA);CHKERRQ(ierr); 245 } 246 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 247 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 248 } 249 250 /* Free workspace */ 251 for (i=0;i<nr;i++) { 252 ierr = ISDestroy(&islrow[i]);CHKERRQ(ierr); 253 } 254 for (i=0;i<nc;i++) { 255 ierr = ISDestroy(&islcol[i]);CHKERRQ(ierr); 256 } 257 ierr = PetscFree4(lr,lc,lrb,lcb);CHKERRQ(ierr); 258 ierr = PetscFree5(isrow,iscol,islrow,islcol,snest);CHKERRQ(ierr); 259 260 /* Create local matrix in MATNEST format */ 261 convert = PETSC_FALSE; 262 ierr = PetscOptionsGetBool(NULL,((PetscObject)A)->prefix,"-matis_convert_local_nest",&convert,NULL);CHKERRQ(ierr); 263 if (convert) { 264 Mat M; 265 266 ierr = MatISGetLocalMat(*newmat,&lA);CHKERRQ(ierr); 267 ierr = MatConvert(lA,MATAIJ,MAT_INITIAL_MATRIX,&M);CHKERRQ(ierr); 268 ierr = MatISSetLocalMat(*newmat,M);CHKERRQ(ierr); 269 ierr = MatDestroy(&M);CHKERRQ(ierr); 270 } 271 272 PetscFunctionReturn(0); 273 } 274 275 #undef __FUNCT__ 276 #define __FUNCT__ "MatTranspose_IS" 277 PetscErrorCode MatTranspose_IS(Mat A,MatReuse reuse,Mat *B) 278 { 279 Mat C,lC,lA; 280 ISLocalToGlobalMapping rl2g,cl2g; 281 PetscBool notset = PETSC_FALSE,cong = PETSC_TRUE; 282 PetscErrorCode ierr; 283 284 PetscFunctionBegin; 285 if (reuse == MAT_REUSE_MATRIX) { 286 PetscBool rcong,ccong; 287 288 ierr = PetscLayoutCompare((*B)->cmap,A->rmap,&rcong);CHKERRQ(ierr); 289 ierr = PetscLayoutCompare((*B)->rmap,A->cmap,&ccong);CHKERRQ(ierr); 290 cong = (PetscBool)(rcong || ccong); 291 } 292 if (reuse == MAT_INITIAL_MATRIX || *B == A || !cong) { 293 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 294 } else { 295 ierr = PetscObjectTypeCompare((PetscObject)*B,MATIS,¬set);CHKERRQ(ierr); 296 C = *B; 297 } 298 299 if (!notset) { 300 ierr = MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N);CHKERRQ(ierr); 301 ierr = MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 302 ierr = MatSetType(C,MATIS);CHKERRQ(ierr); 303 ierr = MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);CHKERRQ(ierr); 304 ierr = MatSetLocalToGlobalMapping(C,cl2g,rl2g);CHKERRQ(ierr); 305 } 306 307 /* perform local transposition */ 308 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr); 309 ierr = MatTranspose(lA,MAT_INITIAL_MATRIX,&lC);CHKERRQ(ierr); 310 ierr = MatISSetLocalMat(C,lC);CHKERRQ(ierr); 311 ierr = MatDestroy(&lC);CHKERRQ(ierr); 312 313 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 314 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 315 316 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 317 if (!cong) { 318 ierr = MatDestroy(B);CHKERRQ(ierr); 319 } 320 *B = C; 321 } else { 322 ierr = MatHeaderMerge(A,&C);CHKERRQ(ierr); 323 } 324 PetscFunctionReturn(0); 325 } 326 327 #undef __FUNCT__ 328 #define __FUNCT__ "MatDiagonalSet_IS" 329 PetscErrorCode MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode) 330 { 331 Mat_IS *is = (Mat_IS*)A->data; 332 PetscErrorCode ierr; 333 334 PetscFunctionBegin; 335 if (!D) { /* this code branch is used by MatShift_IS */ 336 ierr = VecSet(is->y,1.);CHKERRQ(ierr); 337 } else { 338 ierr = VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 339 ierr = VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 340 } 341 ierr = VecPointwiseDivide(is->y,is->y,is->counter);CHKERRQ(ierr); 342 ierr = MatDiagonalSet(is->A,is->y,insmode);CHKERRQ(ierr); 343 PetscFunctionReturn(0); 344 } 345 346 #undef __FUNCT__ 347 #define __FUNCT__ "MatShift_IS" 348 PetscErrorCode MatShift_IS(Mat A,PetscScalar a) 349 { 350 PetscErrorCode ierr; 351 352 PetscFunctionBegin; 353 ierr = MatDiagonalSet_IS(A,NULL,ADD_VALUES);CHKERRQ(ierr); 354 PetscFunctionReturn(0); 355 } 356 357 #undef __FUNCT__ 358 #define __FUNCT__ "MatSetValuesLocal_SubMat_IS" 359 static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 360 { 361 PetscErrorCode ierr; 362 Mat_IS *is = (Mat_IS*)A->data; 363 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 364 365 PetscFunctionBegin; 366 #if defined(PETSC_USE_DEBUG) 367 if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n); 368 #endif 369 ierr = ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l);CHKERRQ(ierr); 370 ierr = ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l);CHKERRQ(ierr); 371 ierr = MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 372 PetscFunctionReturn(0); 373 } 374 375 #undef __FUNCT__ 376 #define __FUNCT__ "MatSetValuesBlockedLocal_SubMat_IS" 377 static PetscErrorCode MatSetValuesBlockedLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 378 { 379 PetscErrorCode ierr; 380 Mat_IS *is = (Mat_IS*)A->data; 381 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 382 383 PetscFunctionBegin; 384 #if defined(PETSC_USE_DEBUG) 385 if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column block indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n); 386 #endif 387 ierr = ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l);CHKERRQ(ierr); 388 ierr = ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l);CHKERRQ(ierr); 389 ierr = MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 390 PetscFunctionReturn(0); 391 } 392 393 #undef __FUNCT__ 394 #define __FUNCT__ "PetscLayoutMapLocal_Private" 395 static PetscErrorCode PetscLayoutMapLocal_Private(PetscLayout map,PetscInt N,const PetscInt idxs[], PetscInt *on,PetscInt **oidxs,PetscInt **ogidxs) 396 { 397 PetscInt *owners = map->range; 398 PetscInt n = map->n; 399 PetscSF sf; 400 PetscInt *lidxs,*work = NULL; 401 PetscSFNode *ridxs; 402 PetscMPIInt rank; 403 PetscInt r, p = 0, len = 0; 404 PetscErrorCode ierr; 405 406 PetscFunctionBegin; 407 /* Create SF where leaves are input idxs and roots are owned idxs (code adapted from MatZeroRowsMapLocal_Private) */ 408 ierr = MPI_Comm_rank(map->comm,&rank);CHKERRQ(ierr); 409 ierr = PetscMalloc1(n,&lidxs);CHKERRQ(ierr); 410 for (r = 0; r < n; ++r) lidxs[r] = -1; 411 ierr = PetscMalloc1(N,&ridxs);CHKERRQ(ierr); 412 for (r = 0; r < N; ++r) { 413 const PetscInt idx = idxs[r]; 414 if (idx < 0 || map->N <= idx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D out of range [0,%D)",idx,map->N); 415 if (idx < owners[p] || owners[p+1] <= idx) { /* short-circuit the search if the last p owns this idx too */ 416 ierr = PetscLayoutFindOwner(map,idx,&p);CHKERRQ(ierr); 417 } 418 ridxs[r].rank = p; 419 ridxs[r].index = idxs[r] - owners[p]; 420 } 421 ierr = PetscSFCreate(map->comm,&sf);CHKERRQ(ierr); 422 ierr = PetscSFSetGraph(sf,n,N,NULL,PETSC_OWN_POINTER,ridxs,PETSC_OWN_POINTER);CHKERRQ(ierr); 423 ierr = PetscSFReduceBegin(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);CHKERRQ(ierr); 424 ierr = PetscSFReduceEnd(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);CHKERRQ(ierr); 425 if (ogidxs) { /* communicate global idxs */ 426 PetscInt cum = 0,start,*work2; 427 428 ierr = PetscMalloc1(n,&work);CHKERRQ(ierr); 429 ierr = PetscCalloc1(N,&work2);CHKERRQ(ierr); 430 for (r = 0; r < N; ++r) if (idxs[r] >=0) cum++; 431 ierr = MPI_Scan(&cum,&start,1,MPIU_INT,MPI_SUM,map->comm);CHKERRQ(ierr); 432 start -= cum; 433 cum = 0; 434 for (r = 0; r < N; ++r) if (idxs[r] >=0) work2[r] = start+cum++; 435 ierr = PetscSFReduceBegin(sf,MPIU_INT,work2,work,MPIU_REPLACE);CHKERRQ(ierr); 436 ierr = PetscSFReduceEnd(sf,MPIU_INT,work2,work,MPIU_REPLACE);CHKERRQ(ierr); 437 ierr = PetscFree(work2);CHKERRQ(ierr); 438 } 439 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 440 /* Compress and put in indices */ 441 for (r = 0; r < n; ++r) 442 if (lidxs[r] >= 0) { 443 if (work) work[len] = work[r]; 444 lidxs[len++] = r; 445 } 446 if (on) *on = len; 447 if (oidxs) *oidxs = lidxs; 448 if (ogidxs) *ogidxs = work; 449 PetscFunctionReturn(0); 450 } 451 452 #undef __FUNCT__ 453 #define __FUNCT__ "MatGetSubMatrix_IS" 454 static PetscErrorCode MatGetSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat) 455 { 456 Mat locmat,newlocmat; 457 Mat_IS *newmatis; 458 #if defined(PETSC_USE_DEBUG) 459 Vec rtest,ltest; 460 const PetscScalar *array; 461 #endif 462 const PetscInt *idxs; 463 PetscInt i,m,n; 464 PetscErrorCode ierr; 465 466 PetscFunctionBegin; 467 if (scall == MAT_REUSE_MATRIX) { 468 PetscBool ismatis; 469 470 ierr = PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);CHKERRQ(ierr); 471 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Not of MATIS type"); 472 newmatis = (Mat_IS*)(*newmat)->data; 473 if (!newmatis->getsub_ris) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local row IS"); 474 if (!newmatis->getsub_cis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local col IS"); 475 } 476 /* irow and icol may not have duplicate entries */ 477 #if defined(PETSC_USE_DEBUG) 478 ierr = MatCreateVecs(mat,<est,&rtest);CHKERRQ(ierr); 479 ierr = ISGetLocalSize(irow,&n);CHKERRQ(ierr); 480 ierr = ISGetIndices(irow,&idxs);CHKERRQ(ierr); 481 for (i=0;i<n;i++) { 482 ierr = VecSetValue(rtest,idxs[i],1.0,ADD_VALUES);CHKERRQ(ierr); 483 } 484 ierr = VecAssemblyBegin(rtest);CHKERRQ(ierr); 485 ierr = VecAssemblyEnd(rtest);CHKERRQ(ierr); 486 ierr = VecGetLocalSize(rtest,&n);CHKERRQ(ierr); 487 ierr = VecGetOwnershipRange(rtest,&m,NULL);CHKERRQ(ierr); 488 ierr = VecGetArrayRead(rtest,&array);CHKERRQ(ierr); 489 for (i=0;i<n;i++) if (array[i] != 0. && array[i] != 1.) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Index %D counted %D times! Irow may not have duplicate entries",i+m,(PetscInt)PetscRealPart(array[i])); 490 ierr = VecRestoreArrayRead(rtest,&array);CHKERRQ(ierr); 491 ierr = ISRestoreIndices(irow,&idxs);CHKERRQ(ierr); 492 ierr = ISGetLocalSize(icol,&n);CHKERRQ(ierr); 493 ierr = ISGetIndices(icol,&idxs);CHKERRQ(ierr); 494 for (i=0;i<n;i++) { 495 ierr = VecSetValue(ltest,idxs[i],1.0,ADD_VALUES);CHKERRQ(ierr); 496 } 497 ierr = VecAssemblyBegin(ltest);CHKERRQ(ierr); 498 ierr = VecAssemblyEnd(ltest);CHKERRQ(ierr); 499 ierr = VecGetLocalSize(ltest,&n);CHKERRQ(ierr); 500 ierr = VecGetOwnershipRange(ltest,&m,NULL);CHKERRQ(ierr); 501 ierr = VecGetArrayRead(ltest,&array);CHKERRQ(ierr); 502 for (i=0;i<n;i++) if (array[i] != 0. && array[i] != 1.) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Index %D counted %D times! Icol may not have duplicate entries",i+m,(PetscInt)PetscRealPart(array[i])); 503 ierr = VecRestoreArrayRead(ltest,&array);CHKERRQ(ierr); 504 ierr = ISRestoreIndices(icol,&idxs);CHKERRQ(ierr); 505 ierr = VecDestroy(&rtest);CHKERRQ(ierr); 506 ierr = VecDestroy(<est);CHKERRQ(ierr); 507 #endif 508 if (scall == MAT_INITIAL_MATRIX) { 509 Mat_IS *matis = (Mat_IS*)mat->data; 510 ISLocalToGlobalMapping rl2g; 511 IS is; 512 PetscInt *lidxs,*lgidxs,*newgidxs; 513 PetscInt ll,newloc; 514 MPI_Comm comm; 515 516 ierr = PetscObjectGetComm((PetscObject)mat,&comm);CHKERRQ(ierr); 517 ierr = ISGetLocalSize(irow,&m);CHKERRQ(ierr); 518 ierr = ISGetLocalSize(icol,&n);CHKERRQ(ierr); 519 ierr = MatCreate(comm,newmat);CHKERRQ(ierr); 520 ierr = MatSetType(*newmat,MATIS);CHKERRQ(ierr); 521 ierr = MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 522 /* communicate irow to their owners in the layout */ 523 ierr = ISGetIndices(irow,&idxs);CHKERRQ(ierr); 524 ierr = PetscLayoutMapLocal_Private(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs);CHKERRQ(ierr); 525 ierr = ISRestoreIndices(irow,&idxs);CHKERRQ(ierr); 526 if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */ 527 ierr = MatISComputeSF_Private(mat);CHKERRQ(ierr); 528 } 529 ierr = PetscMemzero(matis->sf_rootdata,matis->sf_nroots*sizeof(PetscInt));CHKERRQ(ierr); 530 for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1; 531 ierr = PetscFree(lidxs);CHKERRQ(ierr); 532 ierr = PetscFree(lgidxs);CHKERRQ(ierr); 533 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 534 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 535 for (i=0,newloc=0;i<matis->sf_nleaves;i++) if (matis->sf_leafdata[i]) newloc++; 536 ierr = PetscMalloc1(newloc,&newgidxs);CHKERRQ(ierr); 537 ierr = PetscMalloc1(newloc,&lidxs);CHKERRQ(ierr); 538 for (i=0,newloc=0;i<matis->sf_nleaves;i++) 539 if (matis->sf_leafdata[i]) { 540 lidxs[newloc] = i; 541 newgidxs[newloc++] = matis->sf_leafdata[i]-1; 542 } 543 ierr = ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 544 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 545 ierr = ISDestroy(&is);CHKERRQ(ierr); 546 /* local is to extract local submatrix */ 547 newmatis = (Mat_IS*)(*newmat)->data; 548 ierr = ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris);CHKERRQ(ierr); 549 if (mat->congruentlayouts == -1) { /* first time we compare rows and cols layouts */ 550 PetscBool cong; 551 ierr = PetscLayoutCompare(mat->rmap,mat->cmap,&cong);CHKERRQ(ierr); 552 if (cong) mat->congruentlayouts = 1; 553 else mat->congruentlayouts = 0; 554 } 555 if (mat->congruentlayouts && irow == icol) { 556 ierr = MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g);CHKERRQ(ierr); 557 ierr = PetscObjectReference((PetscObject)newmatis->getsub_ris);CHKERRQ(ierr); 558 newmatis->getsub_cis = newmatis->getsub_ris; 559 } else { 560 ISLocalToGlobalMapping cl2g; 561 562 /* communicate icol to their owners in the layout */ 563 ierr = ISGetIndices(icol,&idxs);CHKERRQ(ierr); 564 ierr = PetscLayoutMapLocal_Private(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs);CHKERRQ(ierr); 565 ierr = ISRestoreIndices(icol,&idxs);CHKERRQ(ierr); 566 ierr = PetscMemzero(matis->csf_rootdata,matis->csf_nroots*sizeof(PetscInt));CHKERRQ(ierr); 567 for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1; 568 ierr = PetscFree(lidxs);CHKERRQ(ierr); 569 ierr = PetscFree(lgidxs);CHKERRQ(ierr); 570 ierr = PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);CHKERRQ(ierr); 571 ierr = PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);CHKERRQ(ierr); 572 for (i=0,newloc=0;i<matis->csf_nleaves;i++) if (matis->csf_leafdata[i]) newloc++; 573 ierr = PetscMalloc1(newloc,&newgidxs);CHKERRQ(ierr); 574 ierr = PetscMalloc1(newloc,&lidxs);CHKERRQ(ierr); 575 for (i=0,newloc=0;i<matis->csf_nleaves;i++) 576 if (matis->csf_leafdata[i]) { 577 lidxs[newloc] = i; 578 newgidxs[newloc++] = matis->csf_leafdata[i]-1; 579 } 580 ierr = ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 581 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 582 ierr = ISDestroy(&is);CHKERRQ(ierr); 583 /* local is to extract local submatrix */ 584 ierr = ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis);CHKERRQ(ierr); 585 ierr = MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);CHKERRQ(ierr); 586 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 587 } 588 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 589 } else { 590 ierr = MatISGetLocalMat(*newmat,&newlocmat);CHKERRQ(ierr); 591 } 592 ierr = MatISGetLocalMat(mat,&locmat);CHKERRQ(ierr); 593 newmatis = (Mat_IS*)(*newmat)->data; 594 ierr = MatGetSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat);CHKERRQ(ierr); 595 if (scall == MAT_INITIAL_MATRIX) { 596 ierr = MatISSetLocalMat(*newmat,newlocmat);CHKERRQ(ierr); 597 ierr = MatDestroy(&newlocmat);CHKERRQ(ierr); 598 } 599 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 600 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 601 PetscFunctionReturn(0); 602 } 603 604 #undef __FUNCT__ 605 #define __FUNCT__ "MatCopy_IS" 606 static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str) 607 { 608 Mat_IS *a = (Mat_IS*)A->data,*b; 609 PetscBool ismatis; 610 PetscErrorCode ierr; 611 612 PetscFunctionBegin; 613 ierr = PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);CHKERRQ(ierr); 614 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Need to be implemented"); 615 b = (Mat_IS*)B->data; 616 ierr = MatCopy(a->A,b->A,str);CHKERRQ(ierr); 617 PetscFunctionReturn(0); 618 } 619 620 #undef __FUNCT__ 621 #define __FUNCT__ "MatMissingDiagonal_IS" 622 static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool *missing,PetscInt *d) 623 { 624 Vec v; 625 const PetscScalar *array; 626 PetscInt i,n; 627 PetscErrorCode ierr; 628 629 PetscFunctionBegin; 630 *missing = PETSC_FALSE; 631 ierr = MatCreateVecs(A,NULL,&v);CHKERRQ(ierr); 632 ierr = MatGetDiagonal(A,v);CHKERRQ(ierr); 633 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 634 ierr = VecGetArrayRead(v,&array);CHKERRQ(ierr); 635 for (i=0;i<n;i++) if (array[i] == 0.) break; 636 ierr = VecRestoreArrayRead(v,&array);CHKERRQ(ierr); 637 ierr = VecDestroy(&v);CHKERRQ(ierr); 638 if (i != n) *missing = PETSC_TRUE; 639 if (d) { 640 *d = -1; 641 if (*missing) { 642 PetscInt rstart; 643 ierr = MatGetOwnershipRange(A,&rstart,NULL);CHKERRQ(ierr); 644 *d = i+rstart; 645 } 646 } 647 PetscFunctionReturn(0); 648 } 649 650 #undef __FUNCT__ 651 #define __FUNCT__ "MatISComputeSF_Private" 652 static PetscErrorCode MatISComputeSF_Private(Mat B) 653 { 654 Mat_IS *matis = (Mat_IS*)(B->data); 655 const PetscInt *gidxs; 656 PetscErrorCode ierr; 657 658 PetscFunctionBegin; 659 ierr = MatGetSize(matis->A,&matis->sf_nleaves,NULL);CHKERRQ(ierr); 660 ierr = MatGetLocalSize(B,&matis->sf_nroots,NULL);CHKERRQ(ierr); 661 ierr = PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);CHKERRQ(ierr); 662 ierr = ISLocalToGlobalMappingGetIndices(B->rmap->mapping,&gidxs);CHKERRQ(ierr); 663 ierr = PetscSFSetGraphLayout(matis->sf,B->rmap,matis->sf_nleaves,NULL,PETSC_OWN_POINTER,gidxs);CHKERRQ(ierr); 664 ierr = ISLocalToGlobalMappingRestoreIndices(B->rmap->mapping,&gidxs);CHKERRQ(ierr); 665 ierr = PetscMalloc2(matis->sf_nroots,&matis->sf_rootdata,matis->sf_nleaves,&matis->sf_leafdata);CHKERRQ(ierr); 666 if (B->rmap->mapping != B->cmap->mapping) { /* setup SF for columns */ 667 ierr = MatGetSize(matis->A,NULL,&matis->csf_nleaves);CHKERRQ(ierr); 668 ierr = MatGetLocalSize(B,NULL,&matis->csf_nroots);CHKERRQ(ierr); 669 ierr = PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);CHKERRQ(ierr); 670 ierr = ISLocalToGlobalMappingGetIndices(B->cmap->mapping,&gidxs);CHKERRQ(ierr); 671 ierr = PetscSFSetGraphLayout(matis->csf,B->cmap,matis->csf_nleaves,NULL,PETSC_OWN_POINTER,gidxs);CHKERRQ(ierr); 672 ierr = ISLocalToGlobalMappingRestoreIndices(B->cmap->mapping,&gidxs);CHKERRQ(ierr); 673 ierr = PetscMalloc2(matis->csf_nroots,&matis->csf_rootdata,matis->csf_nleaves,&matis->csf_leafdata);CHKERRQ(ierr); 674 } else { 675 matis->csf = matis->sf; 676 matis->csf_nleaves = matis->sf_nleaves; 677 matis->csf_nroots = matis->sf_nroots; 678 matis->csf_leafdata = matis->sf_leafdata; 679 matis->csf_rootdata = matis->sf_rootdata; 680 } 681 PetscFunctionReturn(0); 682 } 683 684 #undef __FUNCT__ 685 #define __FUNCT__ "MatISSetPreallocation" 686 /*@ 687 MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix. 688 689 Collective on MPI_Comm 690 691 Input Parameters: 692 + B - the matrix 693 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 694 (same value is used for all local rows) 695 . d_nnz - array containing the number of nonzeros in the various rows of the 696 DIAGONAL portion of the local submatrix (possibly different for each row) 697 or NULL, if d_nz is used to specify the nonzero structure. 698 The size of this array is equal to the number of local rows, i.e 'm'. 699 For matrices that will be factored, you must leave room for (and set) 700 the diagonal entry even if it is zero. 701 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 702 submatrix (same value is used for all local rows). 703 - o_nnz - array containing the number of nonzeros in the various rows of the 704 OFF-DIAGONAL portion of the local submatrix (possibly different for 705 each row) or NULL, if o_nz is used to specify the nonzero 706 structure. The size of this array is equal to the number 707 of local rows, i.e 'm'. 708 709 If the *_nnz parameter is given then the *_nz parameter is ignored 710 711 Level: intermediate 712 713 Notes: This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition 714 from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local 715 matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects. 716 717 .keywords: matrix 718 719 .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS 720 @*/ 721 PetscErrorCode MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 722 { 723 PetscErrorCode ierr; 724 725 PetscFunctionBegin; 726 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 727 PetscValidType(B,1); 728 ierr = PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));CHKERRQ(ierr); 729 PetscFunctionReturn(0); 730 } 731 732 #undef __FUNCT__ 733 #define __FUNCT__ "MatISSetPreallocation_IS" 734 static PetscErrorCode MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 735 { 736 Mat_IS *matis = (Mat_IS*)(B->data); 737 PetscInt bs,i,nlocalcols; 738 PetscErrorCode ierr; 739 740 PetscFunctionBegin; 741 if (!matis->A) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"You should first call MatSetLocalToGlobalMapping"); 742 if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */ 743 ierr = MatISComputeSF_Private(B);CHKERRQ(ierr); 744 } 745 if (!d_nnz) { 746 for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] = d_nz; 747 } else { 748 for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] = d_nnz[i]; 749 } 750 if (!o_nnz) { 751 for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] += o_nz; 752 } else { 753 for (i=0;i<matis->sf_nroots;i++) matis->sf_rootdata[i] += o_nnz[i]; 754 } 755 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 756 ierr = MatGetSize(matis->A,NULL,&nlocalcols);CHKERRQ(ierr); 757 ierr = MatGetBlockSize(matis->A,&bs);CHKERRQ(ierr); 758 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 759 for (i=0;i<matis->sf_nleaves;i++) { 760 matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols); 761 } 762 ierr = MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);CHKERRQ(ierr); 763 for (i=0;i<matis->sf_nleaves/bs;i++) { 764 matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs; 765 } 766 ierr = MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);CHKERRQ(ierr); 767 for (i=0;i<matis->sf_nleaves/bs;i++) { 768 matis->sf_leafdata[i] = matis->sf_leafdata[i]-i; 769 } 770 ierr = MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);CHKERRQ(ierr); 771 PetscFunctionReturn(0); 772 } 773 774 #undef __FUNCT__ 775 #define __FUNCT__ "MatISSetMPIXAIJPreallocation_Private" 776 PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce) 777 { 778 Mat_IS *matis = (Mat_IS*)(A->data); 779 PetscInt *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership; 780 const PetscInt *global_indices_r,*global_indices_c; 781 PetscInt i,j,bs,rows,cols; 782 PetscInt lrows,lcols; 783 PetscInt local_rows,local_cols; 784 PetscMPIInt nsubdomains; 785 PetscBool isdense,issbaij; 786 PetscErrorCode ierr; 787 788 PetscFunctionBegin; 789 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&nsubdomains);CHKERRQ(ierr); 790 ierr = MatGetSize(A,&rows,&cols);CHKERRQ(ierr); 791 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 792 ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr); 793 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);CHKERRQ(ierr); 794 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);CHKERRQ(ierr); 795 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);CHKERRQ(ierr); 796 if (A->rmap->mapping != A->cmap->mapping) { 797 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&global_indices_c);CHKERRQ(ierr); 798 } else { 799 global_indices_c = global_indices_r; 800 } 801 802 if (issbaij) { 803 ierr = MatGetRowUpperTriangular(matis->A);CHKERRQ(ierr); 804 } 805 /* 806 An SF reduce is needed to sum up properly on shared rows. 807 Note that generally preallocation is not exact, since it overestimates nonzeros 808 */ 809 if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */ 810 ierr = MatISComputeSF_Private(A);CHKERRQ(ierr); 811 } 812 ierr = MatGetLocalSize(A,&lrows,&lcols);CHKERRQ(ierr); 813 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);CHKERRQ(ierr); 814 /* All processes need to compute entire row ownership */ 815 ierr = PetscMalloc1(rows,&row_ownership);CHKERRQ(ierr); 816 ierr = MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);CHKERRQ(ierr); 817 for (i=0;i<nsubdomains;i++) { 818 for (j=mat_ranges[i];j<mat_ranges[i+1];j++) { 819 row_ownership[j] = i; 820 } 821 } 822 ierr = MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);CHKERRQ(ierr); 823 824 /* 825 my_dnz and my_onz contains exact contribution to preallocation from each local mat 826 then, they will be summed up properly. This way, preallocation is always sufficient 827 */ 828 ierr = PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);CHKERRQ(ierr); 829 /* preallocation as a MATAIJ */ 830 if (isdense) { /* special case for dense local matrices */ 831 for (i=0;i<local_rows;i++) { 832 PetscInt index_row = global_indices_r[i]; 833 for (j=i;j<local_rows;j++) { 834 PetscInt owner = row_ownership[index_row]; 835 PetscInt index_col = global_indices_c[j]; 836 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 837 my_dnz[i] += 1; 838 } else { /* offdiag block */ 839 my_onz[i] += 1; 840 } 841 /* same as before, interchanging rows and cols */ 842 if (i != j) { 843 owner = row_ownership[index_col]; 844 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) { 845 my_dnz[j] += 1; 846 } else { 847 my_onz[j] += 1; 848 } 849 } 850 } 851 } 852 } else if (matis->A->ops->getrowij) { 853 const PetscInt *ii,*jj,*jptr; 854 PetscBool done; 855 ierr = MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);CHKERRQ(ierr); 856 if (!done) SETERRQ1(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ called from %s\n",__FUNCT__); 857 jptr = jj; 858 for (i=0;i<local_rows;i++) { 859 PetscInt index_row = global_indices_r[i]; 860 for (j=0;j<ii[i+1]-ii[i];j++,jptr++) { 861 PetscInt owner = row_ownership[index_row]; 862 PetscInt index_col = global_indices_c[*jptr]; 863 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 864 my_dnz[i] += 1; 865 } else { /* offdiag block */ 866 my_onz[i] += 1; 867 } 868 /* same as before, interchanging rows and cols */ 869 if (issbaij && index_col != index_row) { 870 owner = row_ownership[index_col]; 871 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) { 872 my_dnz[*jptr] += 1; 873 } else { 874 my_onz[*jptr] += 1; 875 } 876 } 877 } 878 } 879 ierr = MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);CHKERRQ(ierr); 880 if (!done) SETERRQ1(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called from %s\n",__FUNCT__); 881 } else { /* loop over rows and use MatGetRow */ 882 for (i=0;i<local_rows;i++) { 883 const PetscInt *cols; 884 PetscInt ncols,index_row = global_indices_r[i]; 885 ierr = MatGetRow(matis->A,i,&ncols,&cols,NULL);CHKERRQ(ierr); 886 for (j=0;j<ncols;j++) { 887 PetscInt owner = row_ownership[index_row]; 888 PetscInt index_col = global_indices_c[cols[j]]; 889 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 890 my_dnz[i] += 1; 891 } else { /* offdiag block */ 892 my_onz[i] += 1; 893 } 894 /* same as before, interchanging rows and cols */ 895 if (issbaij && index_col != index_row) { 896 owner = row_ownership[index_col]; 897 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) { 898 my_dnz[cols[j]] += 1; 899 } else { 900 my_onz[cols[j]] += 1; 901 } 902 } 903 } 904 ierr = MatRestoreRow(matis->A,i,&ncols,&cols,NULL);CHKERRQ(ierr); 905 } 906 } 907 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);CHKERRQ(ierr); 908 if (global_indices_c != global_indices_r) { 909 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&global_indices_c);CHKERRQ(ierr); 910 } 911 ierr = PetscFree(row_ownership);CHKERRQ(ierr); 912 913 /* Reduce my_dnz and my_onz */ 914 if (maxreduce) { 915 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);CHKERRQ(ierr); 916 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);CHKERRQ(ierr); 917 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);CHKERRQ(ierr); 918 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);CHKERRQ(ierr); 919 } else { 920 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);CHKERRQ(ierr); 921 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);CHKERRQ(ierr); 922 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);CHKERRQ(ierr); 923 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);CHKERRQ(ierr); 924 } 925 ierr = PetscFree2(my_dnz,my_onz);CHKERRQ(ierr); 926 927 /* Resize preallocation if overestimated */ 928 for (i=0;i<lrows;i++) { 929 dnz[i] = PetscMin(dnz[i],lcols); 930 onz[i] = PetscMin(onz[i],cols-lcols); 931 } 932 933 /* Set preallocation */ 934 ierr = MatMPIAIJSetPreallocation(B,0,dnz,0,onz);CHKERRQ(ierr); 935 for (i=0;i<lrows/bs;i++) { 936 dnz[i] = dnz[i*bs]/bs; 937 onz[i] = onz[i*bs]/bs; 938 } 939 ierr = MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 940 ierr = MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 941 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 942 if (issbaij) { 943 ierr = MatRestoreRowUpperTriangular(matis->A);CHKERRQ(ierr); 944 } 945 PetscFunctionReturn(0); 946 } 947 948 #undef __FUNCT__ 949 #define __FUNCT__ "MatISGetMPIXAIJ_IS" 950 static PetscErrorCode MatISGetMPIXAIJ_IS(Mat mat, MatReuse reuse, Mat *M) 951 { 952 Mat_IS *matis = (Mat_IS*)(mat->data); 953 Mat local_mat; 954 /* info on mat */ 955 PetscInt bs,rows,cols,lrows,lcols; 956 PetscInt local_rows,local_cols; 957 PetscBool isdense,issbaij,isseqaij; 958 PetscMPIInt nsubdomains; 959 /* values insertion */ 960 PetscScalar *array; 961 /* work */ 962 PetscErrorCode ierr; 963 964 PetscFunctionBegin; 965 /* get info from mat */ 966 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&nsubdomains);CHKERRQ(ierr); 967 if (nsubdomains == 1) { 968 Mat B; 969 IS rows,cols; 970 const PetscInt *ridxs,*cidxs; 971 972 ierr = ISLocalToGlobalMappingGetIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 973 ierr = ISLocalToGlobalMappingGetIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 974 ierr = ISCreateGeneral(PETSC_COMM_SELF,mat->rmap->n,ridxs,PETSC_USE_POINTER,&rows);CHKERRQ(ierr); 975 ierr = ISCreateGeneral(PETSC_COMM_SELF,mat->cmap->n,cidxs,PETSC_USE_POINTER,&cols);CHKERRQ(ierr); 976 ierr = MatConvert(matis->A,MATAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 977 ierr = MatGetSubMatrix(B,rows,cols,reuse,M);CHKERRQ(ierr); 978 ierr = MatDestroy(&B);CHKERRQ(ierr); 979 ierr = ISDestroy(&cols);CHKERRQ(ierr); 980 ierr = ISDestroy(&rows);CHKERRQ(ierr); 981 ierr = ISLocalToGlobalMappingRestoreIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 982 ierr = ISLocalToGlobalMappingRestoreIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 983 PetscFunctionReturn(0); 984 } 985 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 986 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 987 ierr = MatGetLocalSize(mat,&lrows,&lcols);CHKERRQ(ierr); 988 ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr); 989 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);CHKERRQ(ierr); 990 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 991 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);CHKERRQ(ierr); 992 993 if (reuse == MAT_INITIAL_MATRIX) { 994 MatType new_mat_type; 995 PetscBool issbaij_red; 996 997 /* determining new matrix type */ 998 ierr = MPIU_Allreduce(&issbaij,&issbaij_red,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 999 if (issbaij_red) { 1000 new_mat_type = MATSBAIJ; 1001 } else { 1002 if (bs>1) { 1003 new_mat_type = MATBAIJ; 1004 } else { 1005 new_mat_type = MATAIJ; 1006 } 1007 } 1008 1009 ierr = MatCreate(PetscObjectComm((PetscObject)mat),M);CHKERRQ(ierr); 1010 ierr = MatSetSizes(*M,lrows,lcols,rows,cols);CHKERRQ(ierr); 1011 ierr = MatSetBlockSize(*M,bs);CHKERRQ(ierr); 1012 ierr = MatSetType(*M,new_mat_type);CHKERRQ(ierr); 1013 ierr = MatISSetMPIXAIJPreallocation_Private(mat,*M,PETSC_FALSE);CHKERRQ(ierr); 1014 } else { 1015 PetscInt mbs,mrows,mcols,mlrows,mlcols; 1016 /* some checks */ 1017 ierr = MatGetBlockSize(*M,&mbs);CHKERRQ(ierr); 1018 ierr = MatGetSize(*M,&mrows,&mcols);CHKERRQ(ierr); 1019 ierr = MatGetLocalSize(*M,&mlrows,&mlcols);CHKERRQ(ierr); 1020 if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows); 1021 if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols); 1022 if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows); 1023 if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols); 1024 if (mbs != bs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs); 1025 ierr = MatZeroEntries(*M);CHKERRQ(ierr); 1026 } 1027 1028 if (issbaij) { 1029 ierr = MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);CHKERRQ(ierr); 1030 } else { 1031 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1032 local_mat = matis->A; 1033 } 1034 1035 /* Set values */ 1036 ierr = MatSetLocalToGlobalMapping(*M,mat->rmap->mapping,mat->cmap->mapping);CHKERRQ(ierr); 1037 if (isdense) { /* special case for dense local matrices */ 1038 PetscInt i,*dummy_rows,*dummy_cols; 1039 1040 if (local_rows != local_cols) { 1041 ierr = PetscMalloc2(local_rows,&dummy_rows,local_cols,&dummy_cols);CHKERRQ(ierr); 1042 for (i=0;i<local_rows;i++) dummy_rows[i] = i; 1043 for (i=0;i<local_cols;i++) dummy_cols[i] = i; 1044 } else { 1045 ierr = PetscMalloc1(local_rows,&dummy_rows);CHKERRQ(ierr); 1046 for (i=0;i<local_rows;i++) dummy_rows[i] = i; 1047 dummy_cols = dummy_rows; 1048 } 1049 ierr = MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 1050 ierr = MatDenseGetArray(local_mat,&array);CHKERRQ(ierr); 1051 ierr = MatSetValuesLocal(*M,local_rows,dummy_rows,local_cols,dummy_cols,array,ADD_VALUES);CHKERRQ(ierr); 1052 ierr = MatDenseRestoreArray(local_mat,&array);CHKERRQ(ierr); 1053 if (dummy_rows != dummy_cols) { 1054 ierr = PetscFree2(dummy_rows,dummy_cols);CHKERRQ(ierr); 1055 } else { 1056 ierr = PetscFree(dummy_rows);CHKERRQ(ierr); 1057 } 1058 } else if (isseqaij) { 1059 PetscInt i,nvtxs,*xadj,*adjncy; 1060 PetscBool done; 1061 1062 ierr = MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1063 if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ called from %s\n",__FUNCT__); 1064 ierr = MatSeqAIJGetArray(local_mat,&array);CHKERRQ(ierr); 1065 for (i=0;i<nvtxs;i++) { 1066 ierr = MatSetValuesLocal(*M,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);CHKERRQ(ierr); 1067 } 1068 ierr = MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1069 if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called from %s\n",__FUNCT__); 1070 ierr = MatSeqAIJRestoreArray(local_mat,&array);CHKERRQ(ierr); 1071 } else { /* very basic values insertion for all other matrix types */ 1072 PetscInt i; 1073 1074 for (i=0;i<local_rows;i++) { 1075 PetscInt j; 1076 const PetscInt *local_indices_cols; 1077 1078 ierr = MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1079 ierr = MatSetValuesLocal(*M,1,&i,j,local_indices_cols,array,ADD_VALUES);CHKERRQ(ierr); 1080 ierr = MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1081 } 1082 } 1083 ierr = MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1084 ierr = MatDestroy(&local_mat);CHKERRQ(ierr); 1085 ierr = MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1086 if (isdense) { 1087 ierr = MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 1088 } 1089 PetscFunctionReturn(0); 1090 } 1091 1092 #undef __FUNCT__ 1093 #define __FUNCT__ "MatISGetMPIXAIJ" 1094 /*@ 1095 MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format 1096 1097 Input Parameter: 1098 . mat - the matrix (should be of type MATIS) 1099 . reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1100 1101 Output Parameter: 1102 . newmat - the matrix in AIJ format 1103 1104 Level: developer 1105 1106 Notes: mat and *newmat cannot be the same object when MAT_REUSE_MATRIX is requested. 1107 1108 .seealso: MATIS 1109 @*/ 1110 PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat) 1111 { 1112 PetscErrorCode ierr; 1113 1114 PetscFunctionBegin; 1115 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1116 PetscValidLogicalCollectiveEnum(mat,reuse,2); 1117 PetscValidPointer(newmat,3); 1118 if (reuse != MAT_INITIAL_MATRIX) { 1119 PetscValidHeaderSpecific(*newmat,MAT_CLASSID,3); 1120 PetscCheckSameComm(mat,1,*newmat,3); 1121 if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix"); 1122 } 1123 ierr = PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatReuse,Mat*),(mat,reuse,newmat));CHKERRQ(ierr); 1124 PetscFunctionReturn(0); 1125 } 1126 1127 #undef __FUNCT__ 1128 #define __FUNCT__ "MatDuplicate_IS" 1129 PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat) 1130 { 1131 PetscErrorCode ierr; 1132 Mat_IS *matis = (Mat_IS*)(mat->data); 1133 PetscInt bs,m,n,M,N; 1134 Mat B,localmat; 1135 1136 PetscFunctionBegin; 1137 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 1138 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1139 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 1140 ierr = MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,mat->rmap->mapping,mat->cmap->mapping,&B);CHKERRQ(ierr); 1141 ierr = MatDuplicate(matis->A,op,&localmat);CHKERRQ(ierr); 1142 ierr = MatISSetLocalMat(B,localmat);CHKERRQ(ierr); 1143 ierr = MatDestroy(&localmat);CHKERRQ(ierr); 1144 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1145 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1146 *newmat = B; 1147 PetscFunctionReturn(0); 1148 } 1149 1150 #undef __FUNCT__ 1151 #define __FUNCT__ "MatIsHermitian_IS" 1152 static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg) 1153 { 1154 PetscErrorCode ierr; 1155 Mat_IS *matis = (Mat_IS*)A->data; 1156 PetscBool local_sym; 1157 1158 PetscFunctionBegin; 1159 ierr = MatIsHermitian(matis->A,tol,&local_sym);CHKERRQ(ierr); 1160 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1161 PetscFunctionReturn(0); 1162 } 1163 1164 #undef __FUNCT__ 1165 #define __FUNCT__ "MatIsSymmetric_IS" 1166 static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg) 1167 { 1168 PetscErrorCode ierr; 1169 Mat_IS *matis = (Mat_IS*)A->data; 1170 PetscBool local_sym; 1171 1172 PetscFunctionBegin; 1173 ierr = MatIsSymmetric(matis->A,tol,&local_sym);CHKERRQ(ierr); 1174 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1175 PetscFunctionReturn(0); 1176 } 1177 1178 #undef __FUNCT__ 1179 #define __FUNCT__ "MatDestroy_IS" 1180 static PetscErrorCode MatDestroy_IS(Mat A) 1181 { 1182 PetscErrorCode ierr; 1183 Mat_IS *b = (Mat_IS*)A->data; 1184 1185 PetscFunctionBegin; 1186 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 1187 ierr = VecScatterDestroy(&b->cctx);CHKERRQ(ierr); 1188 ierr = VecScatterDestroy(&b->rctx);CHKERRQ(ierr); 1189 ierr = VecDestroy(&b->x);CHKERRQ(ierr); 1190 ierr = VecDestroy(&b->y);CHKERRQ(ierr); 1191 ierr = VecDestroy(&b->counter);CHKERRQ(ierr); 1192 ierr = ISDestroy(&b->getsub_ris);CHKERRQ(ierr); 1193 ierr = ISDestroy(&b->getsub_cis);CHKERRQ(ierr); 1194 if (b->sf != b->csf) { 1195 ierr = PetscSFDestroy(&b->csf);CHKERRQ(ierr); 1196 ierr = PetscFree2(b->csf_rootdata,b->csf_leafdata);CHKERRQ(ierr); 1197 } 1198 ierr = PetscSFDestroy(&b->sf);CHKERRQ(ierr); 1199 ierr = PetscFree2(b->sf_rootdata,b->sf_leafdata);CHKERRQ(ierr); 1200 ierr = PetscFree(A->data);CHKERRQ(ierr); 1201 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 1202 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);CHKERRQ(ierr); 1203 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);CHKERRQ(ierr); 1204 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);CHKERRQ(ierr); 1205 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);CHKERRQ(ierr); 1206 PetscFunctionReturn(0); 1207 } 1208 1209 #undef __FUNCT__ 1210 #define __FUNCT__ "MatMult_IS" 1211 static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y) 1212 { 1213 PetscErrorCode ierr; 1214 Mat_IS *is = (Mat_IS*)A->data; 1215 PetscScalar zero = 0.0; 1216 1217 PetscFunctionBegin; 1218 /* scatter the global vector x into the local work vector */ 1219 ierr = VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1220 ierr = VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1221 1222 /* multiply the local matrix */ 1223 ierr = MatMult(is->A,is->x,is->y);CHKERRQ(ierr); 1224 1225 /* scatter product back into global memory */ 1226 ierr = VecSet(y,zero);CHKERRQ(ierr); 1227 ierr = VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1228 ierr = VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1229 PetscFunctionReturn(0); 1230 } 1231 1232 #undef __FUNCT__ 1233 #define __FUNCT__ "MatMultAdd_IS" 1234 static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 1235 { 1236 Vec temp_vec; 1237 PetscErrorCode ierr; 1238 1239 PetscFunctionBegin; /* v3 = v2 + A * v1.*/ 1240 if (v3 != v2) { 1241 ierr = MatMult(A,v1,v3);CHKERRQ(ierr); 1242 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 1243 } else { 1244 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 1245 ierr = MatMult(A,v1,temp_vec);CHKERRQ(ierr); 1246 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 1247 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 1248 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1249 } 1250 PetscFunctionReturn(0); 1251 } 1252 1253 #undef __FUNCT__ 1254 #define __FUNCT__ "MatMultTranspose_IS" 1255 static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x) 1256 { 1257 Mat_IS *is = (Mat_IS*)A->data; 1258 PetscErrorCode ierr; 1259 1260 PetscFunctionBegin; 1261 /* scatter the global vector x into the local work vector */ 1262 ierr = VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1263 ierr = VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1264 1265 /* multiply the local matrix */ 1266 ierr = MatMultTranspose(is->A,is->y,is->x);CHKERRQ(ierr); 1267 1268 /* scatter product back into global vector */ 1269 ierr = VecSet(x,0);CHKERRQ(ierr); 1270 ierr = VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1271 ierr = VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1272 PetscFunctionReturn(0); 1273 } 1274 1275 #undef __FUNCT__ 1276 #define __FUNCT__ "MatMultTransposeAdd_IS" 1277 static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 1278 { 1279 Vec temp_vec; 1280 PetscErrorCode ierr; 1281 1282 PetscFunctionBegin; /* v3 = v2 + A' * v1.*/ 1283 if (v3 != v2) { 1284 ierr = MatMultTranspose(A,v1,v3);CHKERRQ(ierr); 1285 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 1286 } else { 1287 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 1288 ierr = MatMultTranspose(A,v1,temp_vec);CHKERRQ(ierr); 1289 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 1290 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 1291 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1292 } 1293 PetscFunctionReturn(0); 1294 } 1295 1296 #undef __FUNCT__ 1297 #define __FUNCT__ "MatView_IS" 1298 static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer) 1299 { 1300 Mat_IS *a = (Mat_IS*)A->data; 1301 PetscErrorCode ierr; 1302 PetscViewer sviewer; 1303 1304 PetscFunctionBegin; 1305 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 1306 ierr = MatView(a->A,sviewer);CHKERRQ(ierr); 1307 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 1308 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1309 PetscFunctionReturn(0); 1310 } 1311 1312 #undef __FUNCT__ 1313 #define __FUNCT__ "MatSetLocalToGlobalMapping_IS" 1314 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping) 1315 { 1316 PetscErrorCode ierr; 1317 PetscInt nr,rbs,nc,cbs; 1318 Mat_IS *is = (Mat_IS*)A->data; 1319 IS from,to; 1320 Vec cglobal,rglobal; 1321 1322 PetscFunctionBegin; 1323 PetscCheckSameComm(A,1,rmapping,2); 1324 PetscCheckSameComm(A,1,cmapping,3); 1325 /* Destroy any previous data */ 1326 ierr = VecDestroy(&is->x);CHKERRQ(ierr); 1327 ierr = VecDestroy(&is->y);CHKERRQ(ierr); 1328 ierr = VecDestroy(&is->counter);CHKERRQ(ierr); 1329 ierr = VecScatterDestroy(&is->rctx);CHKERRQ(ierr); 1330 ierr = VecScatterDestroy(&is->cctx);CHKERRQ(ierr); 1331 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 1332 ierr = PetscSFDestroy(&is->sf);CHKERRQ(ierr); 1333 ierr = PetscFree2(is->sf_rootdata,is->sf_leafdata);CHKERRQ(ierr); 1334 1335 /* Setup Layout and set local to global maps */ 1336 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1337 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1338 ierr = PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);CHKERRQ(ierr); 1339 ierr = PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);CHKERRQ(ierr); 1340 1341 /* Create the local matrix A */ 1342 ierr = ISLocalToGlobalMappingGetSize(rmapping,&nr);CHKERRQ(ierr); 1343 ierr = ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);CHKERRQ(ierr); 1344 ierr = ISLocalToGlobalMappingGetSize(cmapping,&nc);CHKERRQ(ierr); 1345 ierr = ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);CHKERRQ(ierr); 1346 ierr = MatCreate(PETSC_COMM_SELF,&is->A);CHKERRQ(ierr); 1347 ierr = MatSetType(is->A,MATAIJ);CHKERRQ(ierr); 1348 ierr = MatSetSizes(is->A,nr,nc,nr,nc);CHKERRQ(ierr); 1349 ierr = MatSetBlockSizes(is->A,rbs,cbs);CHKERRQ(ierr); 1350 ierr = MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);CHKERRQ(ierr); 1351 ierr = MatAppendOptionsPrefix(is->A,"is_");CHKERRQ(ierr); 1352 ierr = MatSetFromOptions(is->A);CHKERRQ(ierr); 1353 1354 if (!is->islocalref) { /* setup scatters and local vectors for MatMult */ 1355 /* Create the local work vectors */ 1356 ierr = MatCreateVecs(is->A,&is->x,&is->y);CHKERRQ(ierr); 1357 1358 /* setup the global to local scatters */ 1359 ierr = MatCreateVecs(A,&cglobal,&rglobal);CHKERRQ(ierr); 1360 ierr = ISCreateStride(PETSC_COMM_SELF,nr,0,1,&to);CHKERRQ(ierr); 1361 ierr = ISLocalToGlobalMappingApplyIS(rmapping,to,&from);CHKERRQ(ierr); 1362 ierr = VecScatterCreate(rglobal,from,is->y,to,&is->rctx);CHKERRQ(ierr); 1363 if (rmapping != cmapping) { 1364 ierr = ISDestroy(&to);CHKERRQ(ierr); 1365 ierr = ISDestroy(&from);CHKERRQ(ierr); 1366 ierr = ISCreateStride(PETSC_COMM_SELF,nc,0,1,&to);CHKERRQ(ierr); 1367 ierr = ISLocalToGlobalMappingApplyIS(cmapping,to,&from);CHKERRQ(ierr); 1368 ierr = VecScatterCreate(cglobal,from,is->x,to,&is->cctx);CHKERRQ(ierr); 1369 } else { 1370 ierr = PetscObjectReference((PetscObject)is->rctx);CHKERRQ(ierr); 1371 is->cctx = is->rctx; 1372 } 1373 1374 /* interface counter vector (local) */ 1375 ierr = VecDuplicate(is->y,&is->counter);CHKERRQ(ierr); 1376 ierr = VecSet(is->y,1.);CHKERRQ(ierr); 1377 ierr = VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1378 ierr = VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1379 ierr = VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1380 ierr = VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1381 1382 /* free workspace */ 1383 ierr = VecDestroy(&rglobal);CHKERRQ(ierr); 1384 ierr = VecDestroy(&cglobal);CHKERRQ(ierr); 1385 ierr = ISDestroy(&to);CHKERRQ(ierr); 1386 ierr = ISDestroy(&from);CHKERRQ(ierr); 1387 } 1388 ierr = MatSetUp(A);CHKERRQ(ierr); 1389 PetscFunctionReturn(0); 1390 } 1391 1392 #undef __FUNCT__ 1393 #define __FUNCT__ "MatSetValues_IS" 1394 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 1395 { 1396 Mat_IS *is = (Mat_IS*)mat->data; 1397 PetscErrorCode ierr; 1398 #if defined(PETSC_USE_DEBUG) 1399 PetscInt i,zm,zn; 1400 #endif 1401 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1402 1403 PetscFunctionBegin; 1404 #if defined(PETSC_USE_DEBUG) 1405 if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n); 1406 /* count negative indices */ 1407 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 1408 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 1409 #endif 1410 ierr = ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 1411 ierr = ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 1412 #if defined(PETSC_USE_DEBUG) 1413 /* count negative indices (should be the same as before) */ 1414 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 1415 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 1416 if (zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS"); 1417 if (zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS"); 1418 #endif 1419 ierr = MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 1420 PetscFunctionReturn(0); 1421 } 1422 1423 #undef __FUNCT__ 1424 #define __FUNCT__ "MatSetValuesBlocked_IS" 1425 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 1426 { 1427 Mat_IS *is = (Mat_IS*)mat->data; 1428 PetscErrorCode ierr; 1429 #if defined(PETSC_USE_DEBUG) 1430 PetscInt i,zm,zn; 1431 #endif 1432 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1433 1434 PetscFunctionBegin; 1435 #if defined(PETSC_USE_DEBUG) 1436 if (m > MATIS_MAX_ENTRIES_INSERTION || n > MATIS_MAX_ENTRIES_INSERTION) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of row/column block indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n); 1437 /* count negative indices */ 1438 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 1439 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 1440 #endif 1441 ierr = ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 1442 ierr = ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 1443 #if defined(PETSC_USE_DEBUG) 1444 /* count negative indices (should be the same as before) */ 1445 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 1446 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 1447 if (zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS"); 1448 if (zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS"); 1449 #endif 1450 ierr = MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 1451 PetscFunctionReturn(0); 1452 } 1453 1454 #undef __FUNCT__ 1455 #define __FUNCT__ "MatSetValuesLocal_IS" 1456 static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 1457 { 1458 PetscErrorCode ierr; 1459 Mat_IS *is = (Mat_IS*)A->data; 1460 1461 PetscFunctionBegin; 1462 ierr = MatSetValues(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 1463 PetscFunctionReturn(0); 1464 } 1465 1466 #undef __FUNCT__ 1467 #define __FUNCT__ "MatSetValuesBlockedLocal_IS" 1468 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 1469 { 1470 PetscErrorCode ierr; 1471 Mat_IS *is = (Mat_IS*)A->data; 1472 1473 PetscFunctionBegin; 1474 ierr = MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 1475 PetscFunctionReturn(0); 1476 } 1477 1478 #undef __FUNCT__ 1479 #define __FUNCT__ "MatISZeroRowsColumnsLocal_Private" 1480 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns) 1481 { 1482 Mat_IS *is = (Mat_IS*)A->data; 1483 PetscErrorCode ierr; 1484 1485 PetscFunctionBegin; 1486 if (!n) { 1487 is->pure_neumann = PETSC_TRUE; 1488 } else { 1489 PetscInt i; 1490 is->pure_neumann = PETSC_FALSE; 1491 1492 if (columns) { 1493 ierr = MatZeroRowsColumns(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 1494 } else { 1495 ierr = MatZeroRows(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 1496 } 1497 if (diag != 0.) { 1498 const PetscScalar *array; 1499 ierr = VecGetArrayRead(is->counter,&array);CHKERRQ(ierr); 1500 for (i=0; i<n; i++) { 1501 ierr = MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);CHKERRQ(ierr); 1502 } 1503 ierr = VecRestoreArrayRead(is->counter,&array);CHKERRQ(ierr); 1504 } 1505 ierr = MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1506 ierr = MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1507 } 1508 PetscFunctionReturn(0); 1509 } 1510 1511 #undef __FUNCT__ 1512 #define __FUNCT__ "MatZeroRowsColumns_Private_IS" 1513 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns) 1514 { 1515 Mat_IS *matis = (Mat_IS*)A->data; 1516 PetscInt nr,nl,len,i; 1517 PetscInt *lrows; 1518 PetscErrorCode ierr; 1519 1520 PetscFunctionBegin; 1521 #if defined(PETSC_USE_DEBUG) 1522 if (columns || diag != 0. || (x && b)) { 1523 PetscBool cong; 1524 ierr = PetscLayoutCompare(A->rmap,A->cmap,&cong);CHKERRQ(ierr); 1525 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 for MATIS"); 1526 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 for MATIS"); 1527 if (!cong && x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"A->rmap and A->cmap need to be congruent"); 1528 } 1529 #endif 1530 /* get locally owned rows */ 1531 ierr = PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);CHKERRQ(ierr); 1532 /* fix right hand side if needed */ 1533 if (x && b) { 1534 const PetscScalar *xx; 1535 PetscScalar *bb; 1536 1537 ierr = VecGetArrayRead(x, &xx);CHKERRQ(ierr); 1538 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 1539 for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]]; 1540 ierr = VecRestoreArrayRead(x, &xx);CHKERRQ(ierr); 1541 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 1542 } 1543 /* get rows associated to the local matrices */ 1544 if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */ 1545 ierr = MatISComputeSF_Private(A);CHKERRQ(ierr); 1546 } 1547 ierr = MatGetSize(matis->A,&nl,NULL);CHKERRQ(ierr); 1548 ierr = PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 1549 ierr = PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1550 for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1; 1551 ierr = PetscFree(lrows);CHKERRQ(ierr); 1552 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1553 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1554 ierr = PetscMalloc1(nl,&lrows);CHKERRQ(ierr); 1555 for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i; 1556 ierr = MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);CHKERRQ(ierr); 1557 ierr = PetscFree(lrows);CHKERRQ(ierr); 1558 PetscFunctionReturn(0); 1559 } 1560 1561 #undef __FUNCT__ 1562 #define __FUNCT__ "MatZeroRows_IS" 1563 static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1564 { 1565 PetscErrorCode ierr; 1566 1567 PetscFunctionBegin; 1568 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);CHKERRQ(ierr); 1569 PetscFunctionReturn(0); 1570 } 1571 1572 #undef __FUNCT__ 1573 #define __FUNCT__ "MatZeroRowsColumns_IS" 1574 static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1575 { 1576 PetscErrorCode ierr; 1577 1578 PetscFunctionBegin; 1579 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);CHKERRQ(ierr); 1580 PetscFunctionReturn(0); 1581 } 1582 1583 #undef __FUNCT__ 1584 #define __FUNCT__ "MatAssemblyBegin_IS" 1585 static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type) 1586 { 1587 Mat_IS *is = (Mat_IS*)A->data; 1588 PetscErrorCode ierr; 1589 1590 PetscFunctionBegin; 1591 ierr = MatAssemblyBegin(is->A,type);CHKERRQ(ierr); 1592 PetscFunctionReturn(0); 1593 } 1594 1595 #undef __FUNCT__ 1596 #define __FUNCT__ "MatAssemblyEnd_IS" 1597 static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type) 1598 { 1599 Mat_IS *is = (Mat_IS*)A->data; 1600 PetscErrorCode ierr; 1601 1602 PetscFunctionBegin; 1603 ierr = MatAssemblyEnd(is->A,type);CHKERRQ(ierr); 1604 PetscFunctionReturn(0); 1605 } 1606 1607 #undef __FUNCT__ 1608 #define __FUNCT__ "MatISGetLocalMat_IS" 1609 static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local) 1610 { 1611 Mat_IS *is = (Mat_IS*)mat->data; 1612 1613 PetscFunctionBegin; 1614 *local = is->A; 1615 PetscFunctionReturn(0); 1616 } 1617 1618 #undef __FUNCT__ 1619 #define __FUNCT__ "MatISGetLocalMat" 1620 /*@ 1621 MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix. 1622 1623 Input Parameter: 1624 . mat - the matrix 1625 1626 Output Parameter: 1627 . local - the local matrix 1628 1629 Level: advanced 1630 1631 Notes: 1632 This can be called if you have precomputed the nonzero structure of the 1633 matrix and want to provide it to the inner matrix object to improve the performance 1634 of the MatSetValues() operation. 1635 1636 .seealso: MATIS 1637 @*/ 1638 PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local) 1639 { 1640 PetscErrorCode ierr; 1641 1642 PetscFunctionBegin; 1643 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1644 PetscValidPointer(local,2); 1645 ierr = PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 1646 PetscFunctionReturn(0); 1647 } 1648 1649 #undef __FUNCT__ 1650 #define __FUNCT__ "MatISSetLocalMat_IS" 1651 static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local) 1652 { 1653 Mat_IS *is = (Mat_IS*)mat->data; 1654 PetscInt nrows,ncols,orows,ocols; 1655 PetscErrorCode ierr; 1656 1657 PetscFunctionBegin; 1658 if (is->A) { 1659 ierr = MatGetSize(is->A,&orows,&ocols);CHKERRQ(ierr); 1660 ierr = MatGetSize(local,&nrows,&ncols);CHKERRQ(ierr); 1661 if (orows != nrows || ocols != ncols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local MATIS matrix should be of size %Dx%D (you passed a %Dx%D matrix)",orows,ocols,nrows,ncols); 1662 } 1663 ierr = PetscObjectReference((PetscObject)local);CHKERRQ(ierr); 1664 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 1665 is->A = local; 1666 PetscFunctionReturn(0); 1667 } 1668 1669 #undef __FUNCT__ 1670 #define __FUNCT__ "MatISSetLocalMat" 1671 /*@ 1672 MatISSetLocalMat - Replace the local matrix stored inside a MATIS object. 1673 1674 Input Parameter: 1675 . mat - the matrix 1676 . local - the local matrix 1677 1678 Output Parameter: 1679 1680 Level: advanced 1681 1682 Notes: 1683 This can be called if you have precomputed the local matrix and 1684 want to provide it to the matrix object MATIS. 1685 1686 .seealso: MATIS 1687 @*/ 1688 PetscErrorCode MatISSetLocalMat(Mat mat,Mat local) 1689 { 1690 PetscErrorCode ierr; 1691 1692 PetscFunctionBegin; 1693 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1694 PetscValidHeaderSpecific(local,MAT_CLASSID,2); 1695 ierr = PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));CHKERRQ(ierr); 1696 PetscFunctionReturn(0); 1697 } 1698 1699 #undef __FUNCT__ 1700 #define __FUNCT__ "MatZeroEntries_IS" 1701 static PetscErrorCode MatZeroEntries_IS(Mat A) 1702 { 1703 Mat_IS *a = (Mat_IS*)A->data; 1704 PetscErrorCode ierr; 1705 1706 PetscFunctionBegin; 1707 ierr = MatZeroEntries(a->A);CHKERRQ(ierr); 1708 PetscFunctionReturn(0); 1709 } 1710 1711 #undef __FUNCT__ 1712 #define __FUNCT__ "MatScale_IS" 1713 static PetscErrorCode MatScale_IS(Mat A,PetscScalar a) 1714 { 1715 Mat_IS *is = (Mat_IS*)A->data; 1716 PetscErrorCode ierr; 1717 1718 PetscFunctionBegin; 1719 ierr = MatScale(is->A,a);CHKERRQ(ierr); 1720 PetscFunctionReturn(0); 1721 } 1722 1723 #undef __FUNCT__ 1724 #define __FUNCT__ "MatGetDiagonal_IS" 1725 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v) 1726 { 1727 Mat_IS *is = (Mat_IS*)A->data; 1728 PetscErrorCode ierr; 1729 1730 PetscFunctionBegin; 1731 /* get diagonal of the local matrix */ 1732 ierr = MatGetDiagonal(is->A,is->y);CHKERRQ(ierr); 1733 1734 /* scatter diagonal back into global vector */ 1735 ierr = VecSet(v,0);CHKERRQ(ierr); 1736 ierr = VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1737 ierr = VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1738 PetscFunctionReturn(0); 1739 } 1740 1741 #undef __FUNCT__ 1742 #define __FUNCT__ "MatSetOption_IS" 1743 static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg) 1744 { 1745 Mat_IS *a = (Mat_IS*)A->data; 1746 PetscErrorCode ierr; 1747 1748 PetscFunctionBegin; 1749 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 1750 PetscFunctionReturn(0); 1751 } 1752 1753 #undef __FUNCT__ 1754 #define __FUNCT__ "MatAXPY_IS" 1755 static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str) 1756 { 1757 Mat_IS *y = (Mat_IS*)Y->data; 1758 Mat_IS *x; 1759 #if defined(PETSC_USE_DEBUG) 1760 PetscBool ismatis; 1761 #endif 1762 PetscErrorCode ierr; 1763 1764 PetscFunctionBegin; 1765 #if defined(PETSC_USE_DEBUG) 1766 ierr = PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);CHKERRQ(ierr); 1767 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS"); 1768 #endif 1769 x = (Mat_IS*)X->data; 1770 ierr = MatAXPY(y->A,a,x->A,str);CHKERRQ(ierr); 1771 PetscFunctionReturn(0); 1772 } 1773 1774 #undef __FUNCT__ 1775 #define __FUNCT__ "MatGetLocalSubMatrix_IS" 1776 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat) 1777 { 1778 Mat lA; 1779 Mat_IS *matis; 1780 ISLocalToGlobalMapping rl2g,cl2g; 1781 IS is; 1782 const PetscInt *rg,*rl; 1783 PetscInt nrg; 1784 PetscInt N,M,nrl,i,*idxs; 1785 PetscErrorCode ierr; 1786 1787 PetscFunctionBegin; 1788 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 1789 ierr = ISGetLocalSize(row,&nrl);CHKERRQ(ierr); 1790 ierr = ISGetIndices(row,&rl);CHKERRQ(ierr); 1791 ierr = ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);CHKERRQ(ierr); 1792 #if defined(PETSC_USE_DEBUG) 1793 for (i=0;i<nrl;i++) if (rl[i]>=nrg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row index %D -> %D greater than maximum possible %D",i,rl[i],nrg); 1794 #endif 1795 ierr = PetscMalloc1(nrg,&idxs);CHKERRQ(ierr); 1796 /* map from [0,nrl) to row */ 1797 for (i=0;i<nrl;i++) idxs[i] = rl[i]; 1798 #if defined(PETSC_USE_DEBUG) 1799 for (i=nrl;i<nrg;i++) idxs[i] = nrg; 1800 #else 1801 for (i=nrl;i<nrg;i++) idxs[i] = -1; 1802 #endif 1803 ierr = ISRestoreIndices(row,&rl);CHKERRQ(ierr); 1804 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 1805 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 1806 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 1807 ierr = ISDestroy(&is);CHKERRQ(ierr); 1808 /* compute new l2g map for columns */ 1809 if (col != row || A->rmap->mapping != A->cmap->mapping) { 1810 const PetscInt *cg,*cl; 1811 PetscInt ncg; 1812 PetscInt ncl; 1813 1814 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 1815 ierr = ISGetLocalSize(col,&ncl);CHKERRQ(ierr); 1816 ierr = ISGetIndices(col,&cl);CHKERRQ(ierr); 1817 ierr = ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);CHKERRQ(ierr); 1818 #if defined(PETSC_USE_DEBUG) 1819 for (i=0;i<ncl;i++) if (cl[i]>=ncg) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local column index %D -> %D greater than maximum possible %D",i,cl[i],ncg); 1820 #endif 1821 ierr = PetscMalloc1(ncg,&idxs);CHKERRQ(ierr); 1822 /* map from [0,ncl) to col */ 1823 for (i=0;i<ncl;i++) idxs[i] = cl[i]; 1824 #if defined(PETSC_USE_DEBUG) 1825 for (i=ncl;i<ncg;i++) idxs[i] = ncg; 1826 #else 1827 for (i=ncl;i<ncg;i++) idxs[i] = -1; 1828 #endif 1829 ierr = ISRestoreIndices(col,&cl);CHKERRQ(ierr); 1830 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 1831 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 1832 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 1833 ierr = ISDestroy(&is);CHKERRQ(ierr); 1834 } else { 1835 ierr = PetscObjectReference((PetscObject)rl2g);CHKERRQ(ierr); 1836 cl2g = rl2g; 1837 } 1838 /* create the MATIS submatrix */ 1839 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 1840 ierr = MatCreate(PetscObjectComm((PetscObject)A),submat);CHKERRQ(ierr); 1841 ierr = MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 1842 ierr = MatSetType(*submat,MATIS);CHKERRQ(ierr); 1843 matis = (Mat_IS*)((*submat)->data); 1844 matis->islocalref = PETSC_TRUE; 1845 ierr = MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);CHKERRQ(ierr); 1846 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr); 1847 ierr = MatISSetLocalMat(*submat,lA);CHKERRQ(ierr); 1848 ierr = MatSetUp(*submat);CHKERRQ(ierr); 1849 ierr = MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1850 ierr = MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1851 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 1852 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 1853 /* remove unsupported ops */ 1854 ierr = PetscMemzero((*submat)->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1855 (*submat)->ops->destroy = MatDestroy_IS; 1856 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS; 1857 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS; 1858 (*submat)->ops->assemblybegin = MatAssemblyBegin_IS; 1859 (*submat)->ops->assemblyend = MatAssemblyEnd_IS; 1860 PetscFunctionReturn(0); 1861 } 1862 1863 #undef __FUNCT__ 1864 #define __FUNCT__ "MatCreateIS" 1865 /*@ 1866 MatCreateIS - Creates a "process" unassembled matrix, assembled on each 1867 process but not across processes. 1868 1869 Input Parameters: 1870 + comm - MPI communicator that will share the matrix 1871 . bs - block size of the matrix 1872 . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products 1873 . rmap - local to global map for rows 1874 - cmap - local to global map for cols 1875 1876 Output Parameter: 1877 . A - the resulting matrix 1878 1879 Level: advanced 1880 1881 Notes: See MATIS for more details. 1882 m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors 1883 used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices. 1884 If either rmap or cmap are NULL, then the matrix is assumed to be square. 1885 1886 .seealso: MATIS, MatSetLocalToGlobalMapping() 1887 @*/ 1888 PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A) 1889 { 1890 PetscErrorCode ierr; 1891 1892 PetscFunctionBegin; 1893 if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings"); 1894 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1895 if (bs > 0) { 1896 ierr = MatSetBlockSize(*A,bs);CHKERRQ(ierr); 1897 } 1898 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 1899 ierr = MatSetType(*A,MATIS);CHKERRQ(ierr); 1900 if (rmap && cmap) { 1901 ierr = MatSetLocalToGlobalMapping(*A,rmap,cmap);CHKERRQ(ierr); 1902 } else if (!rmap) { 1903 ierr = MatSetLocalToGlobalMapping(*A,cmap,cmap);CHKERRQ(ierr); 1904 } else { 1905 ierr = MatSetLocalToGlobalMapping(*A,rmap,rmap);CHKERRQ(ierr); 1906 } 1907 PetscFunctionReturn(0); 1908 } 1909 1910 /*MC 1911 MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP). 1912 This stores the matrices in globally unassembled form. Each processor 1913 assembles only its local Neumann problem and the parallel matrix vector 1914 product is handled "implicitly". 1915 1916 Operations Provided: 1917 + MatMult() 1918 . MatMultAdd() 1919 . MatMultTranspose() 1920 . MatMultTransposeAdd() 1921 . MatZeroEntries() 1922 . MatSetOption() 1923 . MatZeroRows() 1924 . MatSetValues() 1925 . MatSetValuesBlocked() 1926 . MatSetValuesLocal() 1927 . MatSetValuesBlockedLocal() 1928 . MatScale() 1929 . MatGetDiagonal() 1930 . MatMissingDiagonal() 1931 . MatDuplicate() 1932 . MatCopy() 1933 . MatAXPY() 1934 . MatGetSubMatrix() 1935 . MatGetLocalSubMatrix() 1936 . MatTranspose() 1937 - MatSetLocalToGlobalMapping() 1938 1939 Options Database Keys: 1940 . -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions() 1941 1942 Notes: Options prefix for the inner matrix are given by -is_mat_xxx 1943 1944 You must call MatSetLocalToGlobalMapping() before using this matrix type. 1945 1946 You can do matrix preallocation on the local matrix after you obtain it with 1947 MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation() 1948 1949 Level: advanced 1950 1951 .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP 1952 1953 M*/ 1954 1955 #undef __FUNCT__ 1956 #define __FUNCT__ "MatCreate_IS" 1957 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A) 1958 { 1959 PetscErrorCode ierr; 1960 Mat_IS *b; 1961 1962 PetscFunctionBegin; 1963 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 1964 A->data = (void*)b; 1965 1966 /* matrix ops */ 1967 ierr = PetscMemzero(A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1968 A->ops->mult = MatMult_IS; 1969 A->ops->multadd = MatMultAdd_IS; 1970 A->ops->multtranspose = MatMultTranspose_IS; 1971 A->ops->multtransposeadd = MatMultTransposeAdd_IS; 1972 A->ops->destroy = MatDestroy_IS; 1973 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS; 1974 A->ops->setvalues = MatSetValues_IS; 1975 A->ops->setvaluesblocked = MatSetValuesBlocked_IS; 1976 A->ops->setvalueslocal = MatSetValuesLocal_IS; 1977 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS; 1978 A->ops->zerorows = MatZeroRows_IS; 1979 A->ops->zerorowscolumns = MatZeroRowsColumns_IS; 1980 A->ops->assemblybegin = MatAssemblyBegin_IS; 1981 A->ops->assemblyend = MatAssemblyEnd_IS; 1982 A->ops->view = MatView_IS; 1983 A->ops->zeroentries = MatZeroEntries_IS; 1984 A->ops->scale = MatScale_IS; 1985 A->ops->getdiagonal = MatGetDiagonal_IS; 1986 A->ops->setoption = MatSetOption_IS; 1987 A->ops->ishermitian = MatIsHermitian_IS; 1988 A->ops->issymmetric = MatIsSymmetric_IS; 1989 A->ops->duplicate = MatDuplicate_IS; 1990 A->ops->missingdiagonal = MatMissingDiagonal_IS; 1991 A->ops->copy = MatCopy_IS; 1992 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 1993 A->ops->getsubmatrix = MatGetSubMatrix_IS; 1994 A->ops->axpy = MatAXPY_IS; 1995 A->ops->diagonalset = MatDiagonalSet_IS; 1996 A->ops->shift = MatShift_IS; 1997 A->ops->transpose = MatTranspose_IS; 1998 1999 /* special MATIS functions */ 2000 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);CHKERRQ(ierr); 2001 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);CHKERRQ(ierr); 2002 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatISGetMPIXAIJ_IS);CHKERRQ(ierr); 2003 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);CHKERRQ(ierr); 2004 ierr = PetscObjectChangeTypeName((PetscObject)A,MATIS);CHKERRQ(ierr); 2005 PetscFunctionReturn(0); 2006 } 2007