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 /* set preallocation */ 933 ierr = MatMPIAIJSetPreallocation(B,0,dnz,0,onz);CHKERRQ(ierr); 934 for (i=0;i<lrows/bs;i++) { 935 dnz[i] = dnz[i*bs]/bs; 936 onz[i] = onz[i*bs]/bs; 937 } 938 ierr = MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 939 ierr = MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 940 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 941 if (issbaij) { 942 ierr = MatRestoreRowUpperTriangular(matis->A);CHKERRQ(ierr); 943 } 944 PetscFunctionReturn(0); 945 } 946 947 #undef __FUNCT__ 948 #define __FUNCT__ "MatISGetMPIXAIJ_IS" 949 static PetscErrorCode MatISGetMPIXAIJ_IS(Mat mat, MatReuse reuse, Mat *M) 950 { 951 Mat_IS *matis = (Mat_IS*)(mat->data); 952 Mat local_mat; 953 /* info on mat */ 954 PetscInt bs,rows,cols,lrows,lcols; 955 PetscInt local_rows,local_cols; 956 PetscBool isdense,issbaij,isseqaij; 957 PetscMPIInt nsubdomains; 958 /* values insertion */ 959 PetscScalar *array; 960 /* work */ 961 PetscErrorCode ierr; 962 963 PetscFunctionBegin; 964 /* get info from mat */ 965 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&nsubdomains);CHKERRQ(ierr); 966 if (nsubdomains == 1) { 967 if (reuse == MAT_INITIAL_MATRIX) { 968 ierr = MatDuplicate(matis->A,MAT_COPY_VALUES,&(*M));CHKERRQ(ierr); 969 } else { 970 ierr = MatCopy(matis->A,*M,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 971 } 972 PetscFunctionReturn(0); 973 } 974 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 975 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 976 ierr = MatGetLocalSize(mat,&lrows,&lcols);CHKERRQ(ierr); 977 ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr); 978 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);CHKERRQ(ierr); 979 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 980 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);CHKERRQ(ierr); 981 982 if (reuse == MAT_INITIAL_MATRIX) { 983 MatType new_mat_type; 984 PetscBool issbaij_red; 985 986 /* determining new matrix type */ 987 ierr = MPIU_Allreduce(&issbaij,&issbaij_red,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 988 if (issbaij_red) { 989 new_mat_type = MATSBAIJ; 990 } else { 991 if (bs>1) { 992 new_mat_type = MATBAIJ; 993 } else { 994 new_mat_type = MATAIJ; 995 } 996 } 997 998 ierr = MatCreate(PetscObjectComm((PetscObject)mat),M);CHKERRQ(ierr); 999 ierr = MatSetSizes(*M,lrows,lcols,rows,cols);CHKERRQ(ierr); 1000 ierr = MatSetBlockSize(*M,bs);CHKERRQ(ierr); 1001 ierr = MatSetType(*M,new_mat_type);CHKERRQ(ierr); 1002 ierr = MatISSetMPIXAIJPreallocation_Private(mat,*M,PETSC_FALSE);CHKERRQ(ierr); 1003 } else { 1004 PetscInt mbs,mrows,mcols,mlrows,mlcols; 1005 /* some checks */ 1006 ierr = MatGetBlockSize(*M,&mbs);CHKERRQ(ierr); 1007 ierr = MatGetSize(*M,&mrows,&mcols);CHKERRQ(ierr); 1008 ierr = MatGetLocalSize(*M,&mlrows,&mlcols);CHKERRQ(ierr); 1009 if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows); 1010 if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols); 1011 if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows); 1012 if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols); 1013 if (mbs != bs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs); 1014 ierr = MatZeroEntries(*M);CHKERRQ(ierr); 1015 } 1016 1017 if (issbaij) { 1018 ierr = MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);CHKERRQ(ierr); 1019 } else { 1020 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1021 local_mat = matis->A; 1022 } 1023 1024 /* Set values */ 1025 ierr = MatSetLocalToGlobalMapping(*M,mat->rmap->mapping,mat->cmap->mapping);CHKERRQ(ierr); 1026 if (isdense) { /* special case for dense local matrices */ 1027 PetscInt i,*dummy_rows,*dummy_cols; 1028 1029 if (local_rows != local_cols) { 1030 ierr = PetscMalloc2(local_rows,&dummy_rows,local_cols,&dummy_cols);CHKERRQ(ierr); 1031 for (i=0;i<local_rows;i++) dummy_rows[i] = i; 1032 for (i=0;i<local_cols;i++) dummy_cols[i] = i; 1033 } else { 1034 ierr = PetscMalloc1(local_rows,&dummy_rows);CHKERRQ(ierr); 1035 for (i=0;i<local_rows;i++) dummy_rows[i] = i; 1036 dummy_cols = dummy_rows; 1037 } 1038 ierr = MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 1039 ierr = MatDenseGetArray(local_mat,&array);CHKERRQ(ierr); 1040 ierr = MatSetValuesLocal(*M,local_rows,dummy_rows,local_cols,dummy_cols,array,ADD_VALUES);CHKERRQ(ierr); 1041 ierr = MatDenseRestoreArray(local_mat,&array);CHKERRQ(ierr); 1042 if (dummy_rows != dummy_cols) { 1043 ierr = PetscFree2(dummy_rows,dummy_cols);CHKERRQ(ierr); 1044 } else { 1045 ierr = PetscFree(dummy_rows);CHKERRQ(ierr); 1046 } 1047 } else if (isseqaij) { 1048 PetscInt i,nvtxs,*xadj,*adjncy; 1049 PetscBool done; 1050 1051 ierr = MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1052 if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ called from %s\n",__FUNCT__); 1053 ierr = MatSeqAIJGetArray(local_mat,&array);CHKERRQ(ierr); 1054 for (i=0;i<nvtxs;i++) { 1055 ierr = MatSetValuesLocal(*M,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);CHKERRQ(ierr); 1056 } 1057 ierr = MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1058 if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called from %s\n",__FUNCT__); 1059 ierr = MatSeqAIJRestoreArray(local_mat,&array);CHKERRQ(ierr); 1060 } else { /* very basic values insertion for all other matrix types */ 1061 PetscInt i; 1062 1063 for (i=0;i<local_rows;i++) { 1064 PetscInt j; 1065 const PetscInt *local_indices_cols; 1066 1067 ierr = MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1068 ierr = MatSetValuesLocal(*M,1,&i,j,local_indices_cols,array,ADD_VALUES);CHKERRQ(ierr); 1069 ierr = MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1070 } 1071 } 1072 ierr = MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1073 ierr = MatDestroy(&local_mat);CHKERRQ(ierr); 1074 ierr = MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1075 if (isdense) { 1076 ierr = MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 1077 } 1078 PetscFunctionReturn(0); 1079 } 1080 1081 #undef __FUNCT__ 1082 #define __FUNCT__ "MatISGetMPIXAIJ" 1083 /*@ 1084 MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format 1085 1086 Input Parameter: 1087 . mat - the matrix (should be of type MATIS) 1088 . reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1089 1090 Output Parameter: 1091 . newmat - the matrix in AIJ format 1092 1093 Level: developer 1094 1095 Notes: mat and *newmat cannot be the same object when MAT_REUSE_MATRIX is requested. 1096 1097 .seealso: MATIS 1098 @*/ 1099 PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat) 1100 { 1101 PetscErrorCode ierr; 1102 1103 PetscFunctionBegin; 1104 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1105 PetscValidLogicalCollectiveEnum(mat,reuse,2); 1106 PetscValidPointer(newmat,3); 1107 if (reuse != MAT_INITIAL_MATRIX) { 1108 PetscValidHeaderSpecific(*newmat,MAT_CLASSID,3); 1109 PetscCheckSameComm(mat,1,*newmat,3); 1110 if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix"); 1111 } 1112 ierr = PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatReuse,Mat*),(mat,reuse,newmat));CHKERRQ(ierr); 1113 PetscFunctionReturn(0); 1114 } 1115 1116 #undef __FUNCT__ 1117 #define __FUNCT__ "MatDuplicate_IS" 1118 PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat) 1119 { 1120 PetscErrorCode ierr; 1121 Mat_IS *matis = (Mat_IS*)(mat->data); 1122 PetscInt bs,m,n,M,N; 1123 Mat B,localmat; 1124 1125 PetscFunctionBegin; 1126 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 1127 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1128 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 1129 ierr = MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,mat->rmap->mapping,mat->cmap->mapping,&B);CHKERRQ(ierr); 1130 ierr = MatDuplicate(matis->A,op,&localmat);CHKERRQ(ierr); 1131 ierr = MatISSetLocalMat(B,localmat);CHKERRQ(ierr); 1132 ierr = MatDestroy(&localmat);CHKERRQ(ierr); 1133 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1134 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1135 *newmat = B; 1136 PetscFunctionReturn(0); 1137 } 1138 1139 #undef __FUNCT__ 1140 #define __FUNCT__ "MatIsHermitian_IS" 1141 static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg) 1142 { 1143 PetscErrorCode ierr; 1144 Mat_IS *matis = (Mat_IS*)A->data; 1145 PetscBool local_sym; 1146 1147 PetscFunctionBegin; 1148 ierr = MatIsHermitian(matis->A,tol,&local_sym);CHKERRQ(ierr); 1149 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1150 PetscFunctionReturn(0); 1151 } 1152 1153 #undef __FUNCT__ 1154 #define __FUNCT__ "MatIsSymmetric_IS" 1155 static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg) 1156 { 1157 PetscErrorCode ierr; 1158 Mat_IS *matis = (Mat_IS*)A->data; 1159 PetscBool local_sym; 1160 1161 PetscFunctionBegin; 1162 ierr = MatIsSymmetric(matis->A,tol,&local_sym);CHKERRQ(ierr); 1163 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1164 PetscFunctionReturn(0); 1165 } 1166 1167 #undef __FUNCT__ 1168 #define __FUNCT__ "MatDestroy_IS" 1169 static PetscErrorCode MatDestroy_IS(Mat A) 1170 { 1171 PetscErrorCode ierr; 1172 Mat_IS *b = (Mat_IS*)A->data; 1173 1174 PetscFunctionBegin; 1175 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 1176 ierr = VecScatterDestroy(&b->cctx);CHKERRQ(ierr); 1177 ierr = VecScatterDestroy(&b->rctx);CHKERRQ(ierr); 1178 ierr = VecDestroy(&b->x);CHKERRQ(ierr); 1179 ierr = VecDestroy(&b->y);CHKERRQ(ierr); 1180 ierr = VecDestroy(&b->counter);CHKERRQ(ierr); 1181 ierr = ISDestroy(&b->getsub_ris);CHKERRQ(ierr); 1182 ierr = ISDestroy(&b->getsub_cis);CHKERRQ(ierr); 1183 if (b->sf != b->csf) { 1184 ierr = PetscSFDestroy(&b->csf);CHKERRQ(ierr); 1185 ierr = PetscFree2(b->csf_rootdata,b->csf_leafdata);CHKERRQ(ierr); 1186 } 1187 ierr = PetscSFDestroy(&b->sf);CHKERRQ(ierr); 1188 ierr = PetscFree2(b->sf_rootdata,b->sf_leafdata);CHKERRQ(ierr); 1189 ierr = PetscFree(A->data);CHKERRQ(ierr); 1190 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 1191 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);CHKERRQ(ierr); 1192 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);CHKERRQ(ierr); 1193 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);CHKERRQ(ierr); 1194 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);CHKERRQ(ierr); 1195 PetscFunctionReturn(0); 1196 } 1197 1198 #undef __FUNCT__ 1199 #define __FUNCT__ "MatMult_IS" 1200 static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y) 1201 { 1202 PetscErrorCode ierr; 1203 Mat_IS *is = (Mat_IS*)A->data; 1204 PetscScalar zero = 0.0; 1205 1206 PetscFunctionBegin; 1207 /* scatter the global vector x into the local work vector */ 1208 ierr = VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1209 ierr = VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1210 1211 /* multiply the local matrix */ 1212 ierr = MatMult(is->A,is->x,is->y);CHKERRQ(ierr); 1213 1214 /* scatter product back into global memory */ 1215 ierr = VecSet(y,zero);CHKERRQ(ierr); 1216 ierr = VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1217 ierr = VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1218 PetscFunctionReturn(0); 1219 } 1220 1221 #undef __FUNCT__ 1222 #define __FUNCT__ "MatMultAdd_IS" 1223 static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 1224 { 1225 Vec temp_vec; 1226 PetscErrorCode ierr; 1227 1228 PetscFunctionBegin; /* v3 = v2 + A * v1.*/ 1229 if (v3 != v2) { 1230 ierr = MatMult(A,v1,v3);CHKERRQ(ierr); 1231 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 1232 } else { 1233 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 1234 ierr = MatMult(A,v1,temp_vec);CHKERRQ(ierr); 1235 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 1236 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 1237 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1238 } 1239 PetscFunctionReturn(0); 1240 } 1241 1242 #undef __FUNCT__ 1243 #define __FUNCT__ "MatMultTranspose_IS" 1244 static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x) 1245 { 1246 Mat_IS *is = (Mat_IS*)A->data; 1247 PetscErrorCode ierr; 1248 1249 PetscFunctionBegin; 1250 /* scatter the global vector x into the local work vector */ 1251 ierr = VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1252 ierr = VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1253 1254 /* multiply the local matrix */ 1255 ierr = MatMultTranspose(is->A,is->y,is->x);CHKERRQ(ierr); 1256 1257 /* scatter product back into global vector */ 1258 ierr = VecSet(x,0);CHKERRQ(ierr); 1259 ierr = VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1260 ierr = VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1261 PetscFunctionReturn(0); 1262 } 1263 1264 #undef __FUNCT__ 1265 #define __FUNCT__ "MatMultTransposeAdd_IS" 1266 static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 1267 { 1268 Vec temp_vec; 1269 PetscErrorCode ierr; 1270 1271 PetscFunctionBegin; /* v3 = v2 + A' * v1.*/ 1272 if (v3 != v2) { 1273 ierr = MatMultTranspose(A,v1,v3);CHKERRQ(ierr); 1274 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 1275 } else { 1276 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 1277 ierr = MatMultTranspose(A,v1,temp_vec);CHKERRQ(ierr); 1278 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 1279 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 1280 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1281 } 1282 PetscFunctionReturn(0); 1283 } 1284 1285 #undef __FUNCT__ 1286 #define __FUNCT__ "MatView_IS" 1287 static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer) 1288 { 1289 Mat_IS *a = (Mat_IS*)A->data; 1290 PetscErrorCode ierr; 1291 PetscViewer sviewer; 1292 1293 PetscFunctionBegin; 1294 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 1295 ierr = MatView(a->A,sviewer);CHKERRQ(ierr); 1296 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 1297 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1298 PetscFunctionReturn(0); 1299 } 1300 1301 #undef __FUNCT__ 1302 #define __FUNCT__ "MatSetLocalToGlobalMapping_IS" 1303 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping) 1304 { 1305 PetscErrorCode ierr; 1306 PetscInt nr,rbs,nc,cbs; 1307 Mat_IS *is = (Mat_IS*)A->data; 1308 IS from,to; 1309 Vec cglobal,rglobal; 1310 1311 PetscFunctionBegin; 1312 PetscCheckSameComm(A,1,rmapping,2); 1313 PetscCheckSameComm(A,1,cmapping,3); 1314 /* Destroy any previous data */ 1315 ierr = VecDestroy(&is->x);CHKERRQ(ierr); 1316 ierr = VecDestroy(&is->y);CHKERRQ(ierr); 1317 ierr = VecDestroy(&is->counter);CHKERRQ(ierr); 1318 ierr = VecScatterDestroy(&is->rctx);CHKERRQ(ierr); 1319 ierr = VecScatterDestroy(&is->cctx);CHKERRQ(ierr); 1320 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 1321 ierr = PetscSFDestroy(&is->sf);CHKERRQ(ierr); 1322 ierr = PetscFree2(is->sf_rootdata,is->sf_leafdata);CHKERRQ(ierr); 1323 1324 /* Setup Layout and set local to global maps */ 1325 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1326 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1327 ierr = PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);CHKERRQ(ierr); 1328 ierr = PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);CHKERRQ(ierr); 1329 1330 /* Create the local matrix A */ 1331 ierr = ISLocalToGlobalMappingGetSize(rmapping,&nr);CHKERRQ(ierr); 1332 ierr = ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);CHKERRQ(ierr); 1333 ierr = ISLocalToGlobalMappingGetSize(cmapping,&nc);CHKERRQ(ierr); 1334 ierr = ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);CHKERRQ(ierr); 1335 ierr = MatCreate(PETSC_COMM_SELF,&is->A);CHKERRQ(ierr); 1336 ierr = MatSetType(is->A,MATAIJ);CHKERRQ(ierr); 1337 ierr = MatSetSizes(is->A,nr,nc,nr,nc);CHKERRQ(ierr); 1338 ierr = MatSetBlockSizes(is->A,rbs,cbs);CHKERRQ(ierr); 1339 ierr = MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);CHKERRQ(ierr); 1340 ierr = MatAppendOptionsPrefix(is->A,"is_");CHKERRQ(ierr); 1341 ierr = MatSetFromOptions(is->A);CHKERRQ(ierr); 1342 1343 if (!is->islocalref) { /* setup scatters and local vectors for MatMult */ 1344 /* Create the local work vectors */ 1345 ierr = MatCreateVecs(is->A,&is->x,&is->y);CHKERRQ(ierr); 1346 1347 /* setup the global to local scatters */ 1348 ierr = MatCreateVecs(A,&cglobal,&rglobal);CHKERRQ(ierr); 1349 ierr = ISCreateStride(PETSC_COMM_SELF,nr,0,1,&to);CHKERRQ(ierr); 1350 ierr = ISLocalToGlobalMappingApplyIS(rmapping,to,&from);CHKERRQ(ierr); 1351 ierr = VecScatterCreate(rglobal,from,is->y,to,&is->rctx);CHKERRQ(ierr); 1352 if (rmapping != cmapping) { 1353 ierr = ISDestroy(&to);CHKERRQ(ierr); 1354 ierr = ISDestroy(&from);CHKERRQ(ierr); 1355 ierr = ISCreateStride(PETSC_COMM_SELF,nc,0,1,&to);CHKERRQ(ierr); 1356 ierr = ISLocalToGlobalMappingApplyIS(cmapping,to,&from);CHKERRQ(ierr); 1357 ierr = VecScatterCreate(cglobal,from,is->x,to,&is->cctx);CHKERRQ(ierr); 1358 } else { 1359 ierr = PetscObjectReference((PetscObject)is->rctx);CHKERRQ(ierr); 1360 is->cctx = is->rctx; 1361 } 1362 1363 /* interface counter vector (local) */ 1364 ierr = VecDuplicate(is->y,&is->counter);CHKERRQ(ierr); 1365 ierr = VecSet(is->y,1.);CHKERRQ(ierr); 1366 ierr = VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1367 ierr = VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1368 ierr = VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1369 ierr = VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1370 1371 /* free workspace */ 1372 ierr = VecDestroy(&rglobal);CHKERRQ(ierr); 1373 ierr = VecDestroy(&cglobal);CHKERRQ(ierr); 1374 ierr = ISDestroy(&to);CHKERRQ(ierr); 1375 ierr = ISDestroy(&from);CHKERRQ(ierr); 1376 } 1377 ierr = MatSetUp(A);CHKERRQ(ierr); 1378 PetscFunctionReturn(0); 1379 } 1380 1381 #undef __FUNCT__ 1382 #define __FUNCT__ "MatSetValues_IS" 1383 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 1384 { 1385 Mat_IS *is = (Mat_IS*)mat->data; 1386 PetscErrorCode ierr; 1387 #if defined(PETSC_USE_DEBUG) 1388 PetscInt i,zm,zn; 1389 #endif 1390 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1391 1392 PetscFunctionBegin; 1393 #if defined(PETSC_USE_DEBUG) 1394 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); 1395 /* count negative indices */ 1396 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 1397 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 1398 #endif 1399 ierr = ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 1400 ierr = ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 1401 #if defined(PETSC_USE_DEBUG) 1402 /* count negative indices (should be the same as before) */ 1403 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 1404 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 1405 if (zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS"); 1406 if (zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS"); 1407 #endif 1408 ierr = MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 1409 PetscFunctionReturn(0); 1410 } 1411 1412 #undef __FUNCT__ 1413 #define __FUNCT__ "MatSetValuesBlocked_IS" 1414 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 1415 { 1416 Mat_IS *is = (Mat_IS*)mat->data; 1417 PetscErrorCode ierr; 1418 #if defined(PETSC_USE_DEBUG) 1419 PetscInt i,zm,zn; 1420 #endif 1421 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1422 1423 PetscFunctionBegin; 1424 #if defined(PETSC_USE_DEBUG) 1425 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); 1426 /* count negative indices */ 1427 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 1428 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 1429 #endif 1430 ierr = ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 1431 ierr = ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 1432 #if defined(PETSC_USE_DEBUG) 1433 /* count negative indices (should be the same as before) */ 1434 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 1435 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 1436 if (zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS"); 1437 if (zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS"); 1438 #endif 1439 ierr = MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 1440 PetscFunctionReturn(0); 1441 } 1442 1443 #undef __FUNCT__ 1444 #define __FUNCT__ "MatSetValuesLocal_IS" 1445 static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 1446 { 1447 PetscErrorCode ierr; 1448 Mat_IS *is = (Mat_IS*)A->data; 1449 1450 PetscFunctionBegin; 1451 ierr = MatSetValues(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 1452 PetscFunctionReturn(0); 1453 } 1454 1455 #undef __FUNCT__ 1456 #define __FUNCT__ "MatSetValuesBlockedLocal_IS" 1457 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 1458 { 1459 PetscErrorCode ierr; 1460 Mat_IS *is = (Mat_IS*)A->data; 1461 1462 PetscFunctionBegin; 1463 ierr = MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 1464 PetscFunctionReturn(0); 1465 } 1466 1467 #undef __FUNCT__ 1468 #define __FUNCT__ "MatISZeroRowsColumnsLocal_Private" 1469 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns) 1470 { 1471 Mat_IS *is = (Mat_IS*)A->data; 1472 PetscErrorCode ierr; 1473 1474 PetscFunctionBegin; 1475 if (!n) { 1476 is->pure_neumann = PETSC_TRUE; 1477 } else { 1478 PetscInt i; 1479 is->pure_neumann = PETSC_FALSE; 1480 1481 if (columns) { 1482 ierr = MatZeroRowsColumns(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 1483 } else { 1484 ierr = MatZeroRows(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 1485 } 1486 if (diag != 0.) { 1487 const PetscScalar *array; 1488 ierr = VecGetArrayRead(is->counter,&array);CHKERRQ(ierr); 1489 for (i=0; i<n; i++) { 1490 ierr = MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);CHKERRQ(ierr); 1491 } 1492 ierr = VecRestoreArrayRead(is->counter,&array);CHKERRQ(ierr); 1493 } 1494 ierr = MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1495 ierr = MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1496 } 1497 PetscFunctionReturn(0); 1498 } 1499 1500 #undef __FUNCT__ 1501 #define __FUNCT__ "MatZeroRowsColumns_Private_IS" 1502 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns) 1503 { 1504 Mat_IS *matis = (Mat_IS*)A->data; 1505 PetscInt nr,nl,len,i; 1506 PetscInt *lrows; 1507 PetscErrorCode ierr; 1508 1509 PetscFunctionBegin; 1510 #if defined(PETSC_USE_DEBUG) 1511 if (columns || diag != 0. || (x && b)) { 1512 PetscBool cong; 1513 ierr = PetscLayoutCompare(A->rmap,A->cmap,&cong);CHKERRQ(ierr); 1514 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"); 1515 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"); 1516 if (!cong && x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"A->rmap and A->cmap need to be congruent"); 1517 } 1518 #endif 1519 /* get locally owned rows */ 1520 ierr = PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);CHKERRQ(ierr); 1521 /* fix right hand side if needed */ 1522 if (x && b) { 1523 const PetscScalar *xx; 1524 PetscScalar *bb; 1525 1526 ierr = VecGetArrayRead(x, &xx);CHKERRQ(ierr); 1527 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 1528 for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]]; 1529 ierr = VecRestoreArrayRead(x, &xx);CHKERRQ(ierr); 1530 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 1531 } 1532 /* get rows associated to the local matrices */ 1533 if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */ 1534 ierr = MatISComputeSF_Private(A);CHKERRQ(ierr); 1535 } 1536 ierr = MatGetSize(matis->A,&nl,NULL);CHKERRQ(ierr); 1537 ierr = PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 1538 ierr = PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1539 for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1; 1540 ierr = PetscFree(lrows);CHKERRQ(ierr); 1541 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1542 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1543 ierr = PetscMalloc1(nl,&lrows);CHKERRQ(ierr); 1544 for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i; 1545 ierr = MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);CHKERRQ(ierr); 1546 ierr = PetscFree(lrows);CHKERRQ(ierr); 1547 PetscFunctionReturn(0); 1548 } 1549 1550 #undef __FUNCT__ 1551 #define __FUNCT__ "MatZeroRows_IS" 1552 static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1553 { 1554 PetscErrorCode ierr; 1555 1556 PetscFunctionBegin; 1557 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);CHKERRQ(ierr); 1558 PetscFunctionReturn(0); 1559 } 1560 1561 #undef __FUNCT__ 1562 #define __FUNCT__ "MatZeroRowsColumns_IS" 1563 static PetscErrorCode MatZeroRowsColumns_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_TRUE);CHKERRQ(ierr); 1569 PetscFunctionReturn(0); 1570 } 1571 1572 #undef __FUNCT__ 1573 #define __FUNCT__ "MatAssemblyBegin_IS" 1574 static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type) 1575 { 1576 Mat_IS *is = (Mat_IS*)A->data; 1577 PetscErrorCode ierr; 1578 1579 PetscFunctionBegin; 1580 ierr = MatAssemblyBegin(is->A,type);CHKERRQ(ierr); 1581 PetscFunctionReturn(0); 1582 } 1583 1584 #undef __FUNCT__ 1585 #define __FUNCT__ "MatAssemblyEnd_IS" 1586 static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type) 1587 { 1588 Mat_IS *is = (Mat_IS*)A->data; 1589 PetscErrorCode ierr; 1590 1591 PetscFunctionBegin; 1592 ierr = MatAssemblyEnd(is->A,type);CHKERRQ(ierr); 1593 PetscFunctionReturn(0); 1594 } 1595 1596 #undef __FUNCT__ 1597 #define __FUNCT__ "MatISGetLocalMat_IS" 1598 static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local) 1599 { 1600 Mat_IS *is = (Mat_IS*)mat->data; 1601 1602 PetscFunctionBegin; 1603 *local = is->A; 1604 PetscFunctionReturn(0); 1605 } 1606 1607 #undef __FUNCT__ 1608 #define __FUNCT__ "MatISGetLocalMat" 1609 /*@ 1610 MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix. 1611 1612 Input Parameter: 1613 . mat - the matrix 1614 1615 Output Parameter: 1616 . local - the local matrix 1617 1618 Level: advanced 1619 1620 Notes: 1621 This can be called if you have precomputed the nonzero structure of the 1622 matrix and want to provide it to the inner matrix object to improve the performance 1623 of the MatSetValues() operation. 1624 1625 .seealso: MATIS 1626 @*/ 1627 PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local) 1628 { 1629 PetscErrorCode ierr; 1630 1631 PetscFunctionBegin; 1632 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1633 PetscValidPointer(local,2); 1634 ierr = PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 1635 PetscFunctionReturn(0); 1636 } 1637 1638 #undef __FUNCT__ 1639 #define __FUNCT__ "MatISSetLocalMat_IS" 1640 static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local) 1641 { 1642 Mat_IS *is = (Mat_IS*)mat->data; 1643 PetscInt nrows,ncols,orows,ocols; 1644 PetscErrorCode ierr; 1645 1646 PetscFunctionBegin; 1647 if (is->A) { 1648 ierr = MatGetSize(is->A,&orows,&ocols);CHKERRQ(ierr); 1649 ierr = MatGetSize(local,&nrows,&ncols);CHKERRQ(ierr); 1650 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); 1651 } 1652 ierr = PetscObjectReference((PetscObject)local);CHKERRQ(ierr); 1653 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 1654 is->A = local; 1655 PetscFunctionReturn(0); 1656 } 1657 1658 #undef __FUNCT__ 1659 #define __FUNCT__ "MatISSetLocalMat" 1660 /*@ 1661 MatISSetLocalMat - Replace the local matrix stored inside a MATIS object. 1662 1663 Input Parameter: 1664 . mat - the matrix 1665 . local - the local matrix 1666 1667 Output Parameter: 1668 1669 Level: advanced 1670 1671 Notes: 1672 This can be called if you have precomputed the local matrix and 1673 want to provide it to the matrix object MATIS. 1674 1675 .seealso: MATIS 1676 @*/ 1677 PetscErrorCode MatISSetLocalMat(Mat mat,Mat local) 1678 { 1679 PetscErrorCode ierr; 1680 1681 PetscFunctionBegin; 1682 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1683 PetscValidHeaderSpecific(local,MAT_CLASSID,2); 1684 ierr = PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));CHKERRQ(ierr); 1685 PetscFunctionReturn(0); 1686 } 1687 1688 #undef __FUNCT__ 1689 #define __FUNCT__ "MatZeroEntries_IS" 1690 static PetscErrorCode MatZeroEntries_IS(Mat A) 1691 { 1692 Mat_IS *a = (Mat_IS*)A->data; 1693 PetscErrorCode ierr; 1694 1695 PetscFunctionBegin; 1696 ierr = MatZeroEntries(a->A);CHKERRQ(ierr); 1697 PetscFunctionReturn(0); 1698 } 1699 1700 #undef __FUNCT__ 1701 #define __FUNCT__ "MatScale_IS" 1702 static PetscErrorCode MatScale_IS(Mat A,PetscScalar a) 1703 { 1704 Mat_IS *is = (Mat_IS*)A->data; 1705 PetscErrorCode ierr; 1706 1707 PetscFunctionBegin; 1708 ierr = MatScale(is->A,a);CHKERRQ(ierr); 1709 PetscFunctionReturn(0); 1710 } 1711 1712 #undef __FUNCT__ 1713 #define __FUNCT__ "MatGetDiagonal_IS" 1714 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v) 1715 { 1716 Mat_IS *is = (Mat_IS*)A->data; 1717 PetscErrorCode ierr; 1718 1719 PetscFunctionBegin; 1720 /* get diagonal of the local matrix */ 1721 ierr = MatGetDiagonal(is->A,is->y);CHKERRQ(ierr); 1722 1723 /* scatter diagonal back into global vector */ 1724 ierr = VecSet(v,0);CHKERRQ(ierr); 1725 ierr = VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1726 ierr = VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1727 PetscFunctionReturn(0); 1728 } 1729 1730 #undef __FUNCT__ 1731 #define __FUNCT__ "MatSetOption_IS" 1732 static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg) 1733 { 1734 Mat_IS *a = (Mat_IS*)A->data; 1735 PetscErrorCode ierr; 1736 1737 PetscFunctionBegin; 1738 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 1739 PetscFunctionReturn(0); 1740 } 1741 1742 #undef __FUNCT__ 1743 #define __FUNCT__ "MatAXPY_IS" 1744 static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str) 1745 { 1746 Mat_IS *y = (Mat_IS*)Y->data; 1747 Mat_IS *x; 1748 #if defined(PETSC_USE_DEBUG) 1749 PetscBool ismatis; 1750 #endif 1751 PetscErrorCode ierr; 1752 1753 PetscFunctionBegin; 1754 #if defined(PETSC_USE_DEBUG) 1755 ierr = PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);CHKERRQ(ierr); 1756 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS"); 1757 #endif 1758 x = (Mat_IS*)X->data; 1759 ierr = MatAXPY(y->A,a,x->A,str);CHKERRQ(ierr); 1760 PetscFunctionReturn(0); 1761 } 1762 1763 #undef __FUNCT__ 1764 #define __FUNCT__ "MatGetLocalSubMatrix_IS" 1765 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat) 1766 { 1767 Mat lA; 1768 Mat_IS *matis; 1769 ISLocalToGlobalMapping rl2g,cl2g; 1770 IS is; 1771 const PetscInt *rg,*rl; 1772 PetscInt nrg; 1773 PetscInt N,M,nrl,i,*idxs; 1774 PetscErrorCode ierr; 1775 1776 PetscFunctionBegin; 1777 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 1778 ierr = ISGetLocalSize(row,&nrl);CHKERRQ(ierr); 1779 ierr = ISGetIndices(row,&rl);CHKERRQ(ierr); 1780 ierr = ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);CHKERRQ(ierr); 1781 #if defined(PETSC_USE_DEBUG) 1782 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); 1783 #endif 1784 ierr = PetscMalloc1(nrg,&idxs);CHKERRQ(ierr); 1785 /* map from [0,nrl) to row */ 1786 for (i=0;i<nrl;i++) idxs[i] = rl[i]; 1787 #if defined(PETSC_USE_DEBUG) 1788 for (i=nrl;i<nrg;i++) idxs[i] = nrg; 1789 #else 1790 for (i=nrl;i<nrg;i++) idxs[i] = -1; 1791 #endif 1792 ierr = ISRestoreIndices(row,&rl);CHKERRQ(ierr); 1793 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 1794 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 1795 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 1796 ierr = ISDestroy(&is);CHKERRQ(ierr); 1797 /* compute new l2g map for columns */ 1798 if (col != row || A->rmap->mapping != A->cmap->mapping) { 1799 const PetscInt *cg,*cl; 1800 PetscInt ncg; 1801 PetscInt ncl; 1802 1803 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 1804 ierr = ISGetLocalSize(col,&ncl);CHKERRQ(ierr); 1805 ierr = ISGetIndices(col,&cl);CHKERRQ(ierr); 1806 ierr = ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);CHKERRQ(ierr); 1807 #if defined(PETSC_USE_DEBUG) 1808 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); 1809 #endif 1810 ierr = PetscMalloc1(ncg,&idxs);CHKERRQ(ierr); 1811 /* map from [0,ncl) to col */ 1812 for (i=0;i<ncl;i++) idxs[i] = cl[i]; 1813 #if defined(PETSC_USE_DEBUG) 1814 for (i=ncl;i<ncg;i++) idxs[i] = ncg; 1815 #else 1816 for (i=ncl;i<ncg;i++) idxs[i] = -1; 1817 #endif 1818 ierr = ISRestoreIndices(col,&cl);CHKERRQ(ierr); 1819 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 1820 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 1821 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 1822 ierr = ISDestroy(&is);CHKERRQ(ierr); 1823 } else { 1824 ierr = PetscObjectReference((PetscObject)rl2g);CHKERRQ(ierr); 1825 cl2g = rl2g; 1826 } 1827 /* create the MATIS submatrix */ 1828 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 1829 ierr = MatCreate(PetscObjectComm((PetscObject)A),submat);CHKERRQ(ierr); 1830 ierr = MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 1831 ierr = MatSetType(*submat,MATIS);CHKERRQ(ierr); 1832 matis = (Mat_IS*)((*submat)->data); 1833 matis->islocalref = PETSC_TRUE; 1834 ierr = MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);CHKERRQ(ierr); 1835 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr); 1836 ierr = MatISSetLocalMat(*submat,lA);CHKERRQ(ierr); 1837 ierr = MatSetUp(*submat);CHKERRQ(ierr); 1838 ierr = MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1839 ierr = MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1840 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 1841 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 1842 /* remove unsupported ops */ 1843 ierr = PetscMemzero((*submat)->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1844 (*submat)->ops->destroy = MatDestroy_IS; 1845 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS; 1846 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS; 1847 (*submat)->ops->assemblybegin = MatAssemblyBegin_IS; 1848 (*submat)->ops->assemblyend = MatAssemblyEnd_IS; 1849 PetscFunctionReturn(0); 1850 } 1851 1852 #undef __FUNCT__ 1853 #define __FUNCT__ "MatCreateIS" 1854 /*@ 1855 MatCreateIS - Creates a "process" unassembled matrix, assembled on each 1856 process but not across processes. 1857 1858 Input Parameters: 1859 + comm - MPI communicator that will share the matrix 1860 . bs - block size of the matrix 1861 . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products 1862 . rmap - local to global map for rows 1863 - cmap - local to global map for cols 1864 1865 Output Parameter: 1866 . A - the resulting matrix 1867 1868 Level: advanced 1869 1870 Notes: See MATIS for more details. 1871 m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors 1872 used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices. 1873 If either rmap or cmap are NULL, then the matrix is assumed to be square. 1874 1875 .seealso: MATIS, MatSetLocalToGlobalMapping() 1876 @*/ 1877 PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A) 1878 { 1879 PetscErrorCode ierr; 1880 1881 PetscFunctionBegin; 1882 if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings"); 1883 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1884 if (bs > 0) { 1885 ierr = MatSetBlockSize(*A,bs);CHKERRQ(ierr); 1886 } 1887 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 1888 ierr = MatSetType(*A,MATIS);CHKERRQ(ierr); 1889 if (rmap && cmap) { 1890 ierr = MatSetLocalToGlobalMapping(*A,rmap,cmap);CHKERRQ(ierr); 1891 } else if (!rmap) { 1892 ierr = MatSetLocalToGlobalMapping(*A,cmap,cmap);CHKERRQ(ierr); 1893 } else { 1894 ierr = MatSetLocalToGlobalMapping(*A,rmap,rmap);CHKERRQ(ierr); 1895 } 1896 PetscFunctionReturn(0); 1897 } 1898 1899 /*MC 1900 MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP). 1901 This stores the matrices in globally unassembled form. Each processor 1902 assembles only its local Neumann problem and the parallel matrix vector 1903 product is handled "implicitly". 1904 1905 Operations Provided: 1906 + MatMult() 1907 . MatMultAdd() 1908 . MatMultTranspose() 1909 . MatMultTransposeAdd() 1910 . MatZeroEntries() 1911 . MatSetOption() 1912 . MatZeroRows() 1913 . MatSetValues() 1914 . MatSetValuesBlocked() 1915 . MatSetValuesLocal() 1916 . MatSetValuesBlockedLocal() 1917 . MatScale() 1918 . MatGetDiagonal() 1919 . MatMissingDiagonal() 1920 . MatDuplicate() 1921 . MatCopy() 1922 . MatAXPY() 1923 . MatGetSubMatrix() 1924 . MatGetLocalSubMatrix() 1925 . MatTranspose() 1926 - MatSetLocalToGlobalMapping() 1927 1928 Options Database Keys: 1929 . -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions() 1930 1931 Notes: Options prefix for the inner matrix are given by -is_mat_xxx 1932 1933 You must call MatSetLocalToGlobalMapping() before using this matrix type. 1934 1935 You can do matrix preallocation on the local matrix after you obtain it with 1936 MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation() 1937 1938 Level: advanced 1939 1940 .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP 1941 1942 M*/ 1943 1944 #undef __FUNCT__ 1945 #define __FUNCT__ "MatCreate_IS" 1946 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A) 1947 { 1948 PetscErrorCode ierr; 1949 Mat_IS *b; 1950 1951 PetscFunctionBegin; 1952 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 1953 A->data = (void*)b; 1954 1955 /* matrix ops */ 1956 ierr = PetscMemzero(A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1957 A->ops->mult = MatMult_IS; 1958 A->ops->multadd = MatMultAdd_IS; 1959 A->ops->multtranspose = MatMultTranspose_IS; 1960 A->ops->multtransposeadd = MatMultTransposeAdd_IS; 1961 A->ops->destroy = MatDestroy_IS; 1962 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS; 1963 A->ops->setvalues = MatSetValues_IS; 1964 A->ops->setvaluesblocked = MatSetValuesBlocked_IS; 1965 A->ops->setvalueslocal = MatSetValuesLocal_IS; 1966 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS; 1967 A->ops->zerorows = MatZeroRows_IS; 1968 A->ops->zerorowscolumns = MatZeroRowsColumns_IS; 1969 A->ops->assemblybegin = MatAssemblyBegin_IS; 1970 A->ops->assemblyend = MatAssemblyEnd_IS; 1971 A->ops->view = MatView_IS; 1972 A->ops->zeroentries = MatZeroEntries_IS; 1973 A->ops->scale = MatScale_IS; 1974 A->ops->getdiagonal = MatGetDiagonal_IS; 1975 A->ops->setoption = MatSetOption_IS; 1976 A->ops->ishermitian = MatIsHermitian_IS; 1977 A->ops->issymmetric = MatIsSymmetric_IS; 1978 A->ops->duplicate = MatDuplicate_IS; 1979 A->ops->missingdiagonal = MatMissingDiagonal_IS; 1980 A->ops->copy = MatCopy_IS; 1981 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 1982 A->ops->getsubmatrix = MatGetSubMatrix_IS; 1983 A->ops->axpy = MatAXPY_IS; 1984 A->ops->diagonalset = MatDiagonalSet_IS; 1985 A->ops->shift = MatShift_IS; 1986 A->ops->transpose = MatTranspose_IS; 1987 1988 /* special MATIS functions */ 1989 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);CHKERRQ(ierr); 1990 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);CHKERRQ(ierr); 1991 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatISGetMPIXAIJ_IS);CHKERRQ(ierr); 1992 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);CHKERRQ(ierr); 1993 ierr = PetscObjectChangeTypeName((PetscObject)A,MATIS);CHKERRQ(ierr); 1994 PetscFunctionReturn(0); 1995 } 1996