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