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