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