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