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