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 IS irows,icols; 1133 const PetscInt *ridxs,*cidxs; 1134 1135 ierr = ISLocalToGlobalMappingGetIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 1136 ierr = ISLocalToGlobalMappingGetIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 1137 ierr = ISCreateGeneral(PETSC_COMM_SELF,mat->rmap->n,ridxs,PETSC_USE_POINTER,&rows);CHKERRQ(ierr); 1138 ierr = ISCreateGeneral(PETSC_COMM_SELF,mat->cmap->n,cidxs,PETSC_USE_POINTER,&cols);CHKERRQ(ierr); 1139 ierr = ISLocalToGlobalMappingRestoreIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 1140 ierr = ISLocalToGlobalMappingRestoreIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 1141 ierr = ISSetPermutation(rows);CHKERRQ(ierr); 1142 ierr = ISSetPermutation(cols);CHKERRQ(ierr); 1143 ierr = ISInvertPermutation(rows,mat->rmap->n,&irows);CHKERRQ(ierr); 1144 ierr = ISInvertPermutation(cols,mat->cmap->n,&icols);CHKERRQ(ierr); 1145 ierr = ISDestroy(&cols);CHKERRQ(ierr); 1146 ierr = ISDestroy(&rows);CHKERRQ(ierr); 1147 ierr = MatConvert(matis->A,MATAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 1148 ierr = MatGetSubMatrix(B,irows,icols,reuse,M);CHKERRQ(ierr); 1149 ierr = MatDestroy(&B);CHKERRQ(ierr); 1150 ierr = ISDestroy(&icols);CHKERRQ(ierr); 1151 ierr = ISDestroy(&irows);CHKERRQ(ierr); 1152 PetscFunctionReturn(0); 1153 } 1154 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 1155 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 1156 ierr = MatGetLocalSize(mat,&lrows,&lcols);CHKERRQ(ierr); 1157 ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr); 1158 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 1159 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 1160 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);CHKERRQ(ierr); 1161 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);CHKERRQ(ierr); 1162 if (!isseqdense && !isseqaij && !isseqbaij && !isseqsbaij) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name); 1163 #if defined (PETSC_USE_DEBUG) 1164 lb[0] = isseqdense; 1165 lb[1] = isseqaij; 1166 lb[2] = isseqbaij; 1167 lb[3] = isseqsbaij; 1168 ierr = MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1169 if (!bb[0] && !bb[1] && !bb[2] && !bb[3]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type"); 1170 #endif 1171 1172 if (reuse == MAT_INITIAL_MATRIX) { 1173 ierr = MatCreate(PetscObjectComm((PetscObject)mat),M);CHKERRQ(ierr); 1174 ierr = MatSetSizes(*M,lrows,lcols,rows,cols);CHKERRQ(ierr); 1175 ierr = MatSetBlockSize(*M,bs);CHKERRQ(ierr); 1176 if (!isseqsbaij) { 1177 ierr = MatSetType(*M,MATAIJ);CHKERRQ(ierr); 1178 } else { 1179 ierr = MatSetType(*M,MATSBAIJ);CHKERRQ(ierr); 1180 } 1181 ierr = MatISSetMPIXAIJPreallocation_Private(mat,*M,PETSC_FALSE);CHKERRQ(ierr); 1182 } else { 1183 PetscInt mbs,mrows,mcols,mlrows,mlcols; 1184 /* some checks */ 1185 ierr = MatGetBlockSize(*M,&mbs);CHKERRQ(ierr); 1186 ierr = MatGetSize(*M,&mrows,&mcols);CHKERRQ(ierr); 1187 ierr = MatGetLocalSize(*M,&mlrows,&mlcols);CHKERRQ(ierr); 1188 if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows); 1189 if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols); 1190 if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows); 1191 if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols); 1192 if (mbs != bs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs); 1193 ierr = MatZeroEntries(*M);CHKERRQ(ierr); 1194 } 1195 1196 if (isseqsbaij) { 1197 ierr = MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);CHKERRQ(ierr); 1198 } else { 1199 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1200 local_mat = matis->A; 1201 } 1202 1203 /* Set values */ 1204 ierr = MatSetLocalToGlobalMapping(*M,mat->rmap->mapping,mat->cmap->mapping);CHKERRQ(ierr); 1205 if (isseqdense) { /* special case for dense local matrices */ 1206 PetscInt i,*dummy_rows,*dummy_cols; 1207 1208 if (local_rows != local_cols) { 1209 ierr = PetscMalloc2(local_rows,&dummy_rows,local_cols,&dummy_cols);CHKERRQ(ierr); 1210 for (i=0;i<local_rows;i++) dummy_rows[i] = i; 1211 for (i=0;i<local_cols;i++) dummy_cols[i] = i; 1212 } else { 1213 ierr = PetscMalloc1(local_rows,&dummy_rows);CHKERRQ(ierr); 1214 for (i=0;i<local_rows;i++) dummy_rows[i] = i; 1215 dummy_cols = dummy_rows; 1216 } 1217 ierr = MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 1218 ierr = MatDenseGetArray(local_mat,&array);CHKERRQ(ierr); 1219 ierr = MatSetValuesLocal(*M,local_rows,dummy_rows,local_cols,dummy_cols,array,ADD_VALUES);CHKERRQ(ierr); 1220 ierr = MatDenseRestoreArray(local_mat,&array);CHKERRQ(ierr); 1221 if (dummy_rows != dummy_cols) { 1222 ierr = PetscFree2(dummy_rows,dummy_cols);CHKERRQ(ierr); 1223 } else { 1224 ierr = PetscFree(dummy_rows);CHKERRQ(ierr); 1225 } 1226 } else if (isseqaij) { 1227 PetscInt i,nvtxs,*xadj,*adjncy; 1228 PetscBool done; 1229 1230 ierr = MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1231 if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ called from %s\n",__FUNCT__); 1232 ierr = MatSeqAIJGetArray(local_mat,&array);CHKERRQ(ierr); 1233 for (i=0;i<nvtxs;i++) { 1234 ierr = MatSetValuesLocal(*M,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);CHKERRQ(ierr); 1235 } 1236 ierr = MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1237 if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called from %s\n",__FUNCT__); 1238 ierr = MatSeqAIJRestoreArray(local_mat,&array);CHKERRQ(ierr); 1239 } else { /* very basic values insertion for all other matrix types */ 1240 PetscInt i; 1241 1242 for (i=0;i<local_rows;i++) { 1243 PetscInt j; 1244 const PetscInt *local_indices_cols; 1245 1246 ierr = MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1247 ierr = MatSetValuesLocal(*M,1,&i,j,local_indices_cols,array,ADD_VALUES);CHKERRQ(ierr); 1248 ierr = MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1249 } 1250 } 1251 ierr = MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1252 ierr = MatDestroy(&local_mat);CHKERRQ(ierr); 1253 ierr = MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1254 if (isseqdense) { 1255 ierr = MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 1256 } 1257 PetscFunctionReturn(0); 1258 } 1259 1260 #undef __FUNCT__ 1261 #define __FUNCT__ "MatISGetMPIXAIJ" 1262 /*@ 1263 MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format 1264 1265 Input Parameter: 1266 . mat - the matrix (should be of type MATIS) 1267 . reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1268 1269 Output Parameter: 1270 . newmat - the matrix in AIJ format 1271 1272 Level: developer 1273 1274 Notes: mat and *newmat cannot be the same object when MAT_REUSE_MATRIX is requested. 1275 1276 .seealso: MATIS 1277 @*/ 1278 PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat) 1279 { 1280 PetscErrorCode ierr; 1281 1282 PetscFunctionBegin; 1283 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1284 PetscValidLogicalCollectiveEnum(mat,reuse,2); 1285 PetscValidPointer(newmat,3); 1286 if (reuse != MAT_INITIAL_MATRIX) { 1287 PetscValidHeaderSpecific(*newmat,MAT_CLASSID,3); 1288 PetscCheckSameComm(mat,1,*newmat,3); 1289 if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix"); 1290 } 1291 ierr = PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatReuse,Mat*),(mat,reuse,newmat));CHKERRQ(ierr); 1292 PetscFunctionReturn(0); 1293 } 1294 1295 #undef __FUNCT__ 1296 #define __FUNCT__ "MatDuplicate_IS" 1297 PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat) 1298 { 1299 PetscErrorCode ierr; 1300 Mat_IS *matis = (Mat_IS*)(mat->data); 1301 PetscInt bs,m,n,M,N; 1302 Mat B,localmat; 1303 1304 PetscFunctionBegin; 1305 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 1306 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1307 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 1308 ierr = MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,mat->rmap->mapping,mat->cmap->mapping,&B);CHKERRQ(ierr); 1309 ierr = MatDuplicate(matis->A,op,&localmat);CHKERRQ(ierr); 1310 ierr = MatISSetLocalMat(B,localmat);CHKERRQ(ierr); 1311 ierr = MatDestroy(&localmat);CHKERRQ(ierr); 1312 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1313 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1314 *newmat = B; 1315 PetscFunctionReturn(0); 1316 } 1317 1318 #undef __FUNCT__ 1319 #define __FUNCT__ "MatIsHermitian_IS" 1320 static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg) 1321 { 1322 PetscErrorCode ierr; 1323 Mat_IS *matis = (Mat_IS*)A->data; 1324 PetscBool local_sym; 1325 1326 PetscFunctionBegin; 1327 ierr = MatIsHermitian(matis->A,tol,&local_sym);CHKERRQ(ierr); 1328 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1329 PetscFunctionReturn(0); 1330 } 1331 1332 #undef __FUNCT__ 1333 #define __FUNCT__ "MatIsSymmetric_IS" 1334 static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg) 1335 { 1336 PetscErrorCode ierr; 1337 Mat_IS *matis = (Mat_IS*)A->data; 1338 PetscBool local_sym; 1339 1340 PetscFunctionBegin; 1341 ierr = MatIsSymmetric(matis->A,tol,&local_sym);CHKERRQ(ierr); 1342 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1343 PetscFunctionReturn(0); 1344 } 1345 1346 #undef __FUNCT__ 1347 #define __FUNCT__ "MatIsStructurallySymmetric_IS" 1348 static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool *flg) 1349 { 1350 PetscErrorCode ierr; 1351 Mat_IS *matis = (Mat_IS*)A->data; 1352 PetscBool local_sym; 1353 1354 PetscFunctionBegin; 1355 if (A->rmap->mapping != A->cmap->mapping) { 1356 *flg = PETSC_FALSE; 1357 PetscFunctionReturn(0); 1358 } 1359 ierr = MatIsStructurallySymmetric(matis->A,&local_sym);CHKERRQ(ierr); 1360 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1361 PetscFunctionReturn(0); 1362 } 1363 1364 #undef __FUNCT__ 1365 #define __FUNCT__ "MatDestroy_IS" 1366 static PetscErrorCode MatDestroy_IS(Mat A) 1367 { 1368 PetscErrorCode ierr; 1369 Mat_IS *b = (Mat_IS*)A->data; 1370 1371 PetscFunctionBegin; 1372 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 1373 ierr = VecScatterDestroy(&b->cctx);CHKERRQ(ierr); 1374 ierr = VecScatterDestroy(&b->rctx);CHKERRQ(ierr); 1375 ierr = VecDestroy(&b->x);CHKERRQ(ierr); 1376 ierr = VecDestroy(&b->y);CHKERRQ(ierr); 1377 ierr = VecDestroy(&b->counter);CHKERRQ(ierr); 1378 ierr = ISDestroy(&b->getsub_ris);CHKERRQ(ierr); 1379 ierr = ISDestroy(&b->getsub_cis);CHKERRQ(ierr); 1380 if (b->sf != b->csf) { 1381 ierr = PetscSFDestroy(&b->csf);CHKERRQ(ierr); 1382 ierr = PetscFree2(b->csf_rootdata,b->csf_leafdata);CHKERRQ(ierr); 1383 } 1384 ierr = PetscSFDestroy(&b->sf);CHKERRQ(ierr); 1385 ierr = PetscFree2(b->sf_rootdata,b->sf_leafdata);CHKERRQ(ierr); 1386 ierr = PetscFree(A->data);CHKERRQ(ierr); 1387 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 1388 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);CHKERRQ(ierr); 1389 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);CHKERRQ(ierr); 1390 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);CHKERRQ(ierr); 1391 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);CHKERRQ(ierr); 1392 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",NULL);CHKERRQ(ierr); 1393 PetscFunctionReturn(0); 1394 } 1395 1396 #undef __FUNCT__ 1397 #define __FUNCT__ "MatMult_IS" 1398 static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y) 1399 { 1400 PetscErrorCode ierr; 1401 Mat_IS *is = (Mat_IS*)A->data; 1402 PetscScalar zero = 0.0; 1403 1404 PetscFunctionBegin; 1405 /* scatter the global vector x into the local work vector */ 1406 ierr = VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1407 ierr = VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1408 1409 /* multiply the local matrix */ 1410 ierr = MatMult(is->A,is->x,is->y);CHKERRQ(ierr); 1411 1412 /* scatter product back into global memory */ 1413 ierr = VecSet(y,zero);CHKERRQ(ierr); 1414 ierr = VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1415 ierr = VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1416 PetscFunctionReturn(0); 1417 } 1418 1419 #undef __FUNCT__ 1420 #define __FUNCT__ "MatMultAdd_IS" 1421 static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 1422 { 1423 Vec temp_vec; 1424 PetscErrorCode ierr; 1425 1426 PetscFunctionBegin; /* v3 = v2 + A * v1.*/ 1427 if (v3 != v2) { 1428 ierr = MatMult(A,v1,v3);CHKERRQ(ierr); 1429 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 1430 } else { 1431 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 1432 ierr = MatMult(A,v1,temp_vec);CHKERRQ(ierr); 1433 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 1434 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 1435 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1436 } 1437 PetscFunctionReturn(0); 1438 } 1439 1440 #undef __FUNCT__ 1441 #define __FUNCT__ "MatMultTranspose_IS" 1442 static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x) 1443 { 1444 Mat_IS *is = (Mat_IS*)A->data; 1445 PetscErrorCode ierr; 1446 1447 PetscFunctionBegin; 1448 /* scatter the global vector x into the local work vector */ 1449 ierr = VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1450 ierr = VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1451 1452 /* multiply the local matrix */ 1453 ierr = MatMultTranspose(is->A,is->y,is->x);CHKERRQ(ierr); 1454 1455 /* scatter product back into global vector */ 1456 ierr = VecSet(x,0);CHKERRQ(ierr); 1457 ierr = VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1458 ierr = VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1459 PetscFunctionReturn(0); 1460 } 1461 1462 #undef __FUNCT__ 1463 #define __FUNCT__ "MatMultTransposeAdd_IS" 1464 static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 1465 { 1466 Vec temp_vec; 1467 PetscErrorCode ierr; 1468 1469 PetscFunctionBegin; /* v3 = v2 + A' * v1.*/ 1470 if (v3 != v2) { 1471 ierr = MatMultTranspose(A,v1,v3);CHKERRQ(ierr); 1472 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 1473 } else { 1474 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 1475 ierr = MatMultTranspose(A,v1,temp_vec);CHKERRQ(ierr); 1476 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 1477 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 1478 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1479 } 1480 PetscFunctionReturn(0); 1481 } 1482 1483 #undef __FUNCT__ 1484 #define __FUNCT__ "MatView_IS" 1485 static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer) 1486 { 1487 Mat_IS *a = (Mat_IS*)A->data; 1488 PetscErrorCode ierr; 1489 PetscViewer sviewer; 1490 PetscBool isascii,view = PETSC_TRUE; 1491 1492 PetscFunctionBegin; 1493 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 1494 if (isascii) { 1495 PetscViewerFormat format; 1496 1497 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1498 if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE; 1499 } 1500 if (!view) PetscFunctionReturn(0); 1501 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 1502 ierr = MatView(a->A,sviewer);CHKERRQ(ierr); 1503 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 1504 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1505 PetscFunctionReturn(0); 1506 } 1507 1508 #undef __FUNCT__ 1509 #define __FUNCT__ "MatSetLocalToGlobalMapping_IS" 1510 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping) 1511 { 1512 PetscErrorCode ierr; 1513 PetscInt nr,rbs,nc,cbs; 1514 Mat_IS *is = (Mat_IS*)A->data; 1515 IS from,to; 1516 Vec cglobal,rglobal; 1517 1518 PetscFunctionBegin; 1519 PetscCheckSameComm(A,1,rmapping,2); 1520 PetscCheckSameComm(A,1,cmapping,3); 1521 /* Destroy any previous data */ 1522 ierr = VecDestroy(&is->x);CHKERRQ(ierr); 1523 ierr = VecDestroy(&is->y);CHKERRQ(ierr); 1524 ierr = VecDestroy(&is->counter);CHKERRQ(ierr); 1525 ierr = VecScatterDestroy(&is->rctx);CHKERRQ(ierr); 1526 ierr = VecScatterDestroy(&is->cctx);CHKERRQ(ierr); 1527 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 1528 ierr = PetscSFDestroy(&is->sf);CHKERRQ(ierr); 1529 ierr = PetscFree2(is->sf_rootdata,is->sf_leafdata);CHKERRQ(ierr); 1530 1531 /* Setup Layout and set local to global maps */ 1532 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1533 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1534 ierr = PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);CHKERRQ(ierr); 1535 ierr = PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);CHKERRQ(ierr); 1536 1537 /* Create the local matrix A */ 1538 ierr = ISLocalToGlobalMappingGetSize(rmapping,&nr);CHKERRQ(ierr); 1539 ierr = ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);CHKERRQ(ierr); 1540 ierr = ISLocalToGlobalMappingGetSize(cmapping,&nc);CHKERRQ(ierr); 1541 ierr = ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);CHKERRQ(ierr); 1542 ierr = MatCreate(PETSC_COMM_SELF,&is->A);CHKERRQ(ierr); 1543 ierr = MatSetType(is->A,MATAIJ);CHKERRQ(ierr); 1544 ierr = MatSetSizes(is->A,nr,nc,nr,nc);CHKERRQ(ierr); 1545 ierr = MatSetBlockSizes(is->A,rbs,cbs);CHKERRQ(ierr); 1546 ierr = MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);CHKERRQ(ierr); 1547 ierr = MatAppendOptionsPrefix(is->A,"is_");CHKERRQ(ierr); 1548 ierr = MatSetFromOptions(is->A);CHKERRQ(ierr); 1549 1550 if (!is->islocalref) { /* setup scatters and local vectors for MatMult */ 1551 /* Create the local work vectors */ 1552 ierr = MatCreateVecs(is->A,&is->x,&is->y);CHKERRQ(ierr); 1553 1554 /* setup the global to local scatters */ 1555 ierr = MatCreateVecs(A,&cglobal,&rglobal);CHKERRQ(ierr); 1556 ierr = ISCreateStride(PETSC_COMM_SELF,nr,0,1,&to);CHKERRQ(ierr); 1557 ierr = ISLocalToGlobalMappingApplyIS(rmapping,to,&from);CHKERRQ(ierr); 1558 ierr = VecScatterCreate(rglobal,from,is->y,to,&is->rctx);CHKERRQ(ierr); 1559 if (rmapping != cmapping) { 1560 ierr = ISDestroy(&to);CHKERRQ(ierr); 1561 ierr = ISDestroy(&from);CHKERRQ(ierr); 1562 ierr = ISCreateStride(PETSC_COMM_SELF,nc,0,1,&to);CHKERRQ(ierr); 1563 ierr = ISLocalToGlobalMappingApplyIS(cmapping,to,&from);CHKERRQ(ierr); 1564 ierr = VecScatterCreate(cglobal,from,is->x,to,&is->cctx);CHKERRQ(ierr); 1565 } else { 1566 ierr = PetscObjectReference((PetscObject)is->rctx);CHKERRQ(ierr); 1567 is->cctx = is->rctx; 1568 } 1569 1570 /* interface counter vector (local) */ 1571 ierr = VecDuplicate(is->y,&is->counter);CHKERRQ(ierr); 1572 ierr = VecSet(is->y,1.);CHKERRQ(ierr); 1573 ierr = VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1574 ierr = VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1575 ierr = VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1576 ierr = VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1577 1578 /* free workspace */ 1579 ierr = VecDestroy(&rglobal);CHKERRQ(ierr); 1580 ierr = VecDestroy(&cglobal);CHKERRQ(ierr); 1581 ierr = ISDestroy(&to);CHKERRQ(ierr); 1582 ierr = ISDestroy(&from);CHKERRQ(ierr); 1583 } 1584 ierr = MatSetUp(A);CHKERRQ(ierr); 1585 PetscFunctionReturn(0); 1586 } 1587 1588 #undef __FUNCT__ 1589 #define __FUNCT__ "MatSetValues_IS" 1590 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 1591 { 1592 Mat_IS *is = (Mat_IS*)mat->data; 1593 PetscErrorCode ierr; 1594 #if defined(PETSC_USE_DEBUG) 1595 PetscInt i,zm,zn; 1596 #endif 1597 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1598 1599 PetscFunctionBegin; 1600 #if defined(PETSC_USE_DEBUG) 1601 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); 1602 /* count negative indices */ 1603 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 1604 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 1605 #endif 1606 ierr = ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 1607 ierr = ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 1608 #if defined(PETSC_USE_DEBUG) 1609 /* count negative indices (should be the same as before) */ 1610 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 1611 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 1612 if (zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS"); 1613 if (zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS"); 1614 #endif 1615 ierr = MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 1616 PetscFunctionReturn(0); 1617 } 1618 1619 #undef __FUNCT__ 1620 #define __FUNCT__ "MatSetValuesBlocked_IS" 1621 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 1622 { 1623 Mat_IS *is = (Mat_IS*)mat->data; 1624 PetscErrorCode ierr; 1625 #if defined(PETSC_USE_DEBUG) 1626 PetscInt i,zm,zn; 1627 #endif 1628 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1629 1630 PetscFunctionBegin; 1631 #if defined(PETSC_USE_DEBUG) 1632 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); 1633 /* count negative indices */ 1634 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 1635 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 1636 #endif 1637 ierr = ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 1638 ierr = ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 1639 #if defined(PETSC_USE_DEBUG) 1640 /* count negative indices (should be the same as before) */ 1641 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 1642 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 1643 if (zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS"); 1644 if (zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS"); 1645 #endif 1646 ierr = MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 1647 PetscFunctionReturn(0); 1648 } 1649 1650 #undef __FUNCT__ 1651 #define __FUNCT__ "MatSetValuesLocal_IS" 1652 static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 1653 { 1654 PetscErrorCode ierr; 1655 Mat_IS *is = (Mat_IS*)A->data; 1656 1657 PetscFunctionBegin; 1658 ierr = MatSetValues(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 1659 PetscFunctionReturn(0); 1660 } 1661 1662 #undef __FUNCT__ 1663 #define __FUNCT__ "MatSetValuesBlockedLocal_IS" 1664 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 1665 { 1666 PetscErrorCode ierr; 1667 Mat_IS *is = (Mat_IS*)A->data; 1668 1669 PetscFunctionBegin; 1670 ierr = MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 1671 PetscFunctionReturn(0); 1672 } 1673 1674 #undef __FUNCT__ 1675 #define __FUNCT__ "MatISZeroRowsColumnsLocal_Private" 1676 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns) 1677 { 1678 Mat_IS *is = (Mat_IS*)A->data; 1679 PetscErrorCode ierr; 1680 1681 PetscFunctionBegin; 1682 if (!n) { 1683 is->pure_neumann = PETSC_TRUE; 1684 } else { 1685 PetscInt i; 1686 is->pure_neumann = PETSC_FALSE; 1687 1688 if (columns) { 1689 ierr = MatZeroRowsColumns(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 1690 } else { 1691 ierr = MatZeroRows(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 1692 } 1693 if (diag != 0.) { 1694 const PetscScalar *array; 1695 ierr = VecGetArrayRead(is->counter,&array);CHKERRQ(ierr); 1696 for (i=0; i<n; i++) { 1697 ierr = MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);CHKERRQ(ierr); 1698 } 1699 ierr = VecRestoreArrayRead(is->counter,&array);CHKERRQ(ierr); 1700 } 1701 ierr = MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1702 ierr = MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1703 } 1704 PetscFunctionReturn(0); 1705 } 1706 1707 #undef __FUNCT__ 1708 #define __FUNCT__ "MatZeroRowsColumns_Private_IS" 1709 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns) 1710 { 1711 Mat_IS *matis = (Mat_IS*)A->data; 1712 PetscInt nr,nl,len,i; 1713 PetscInt *lrows; 1714 PetscErrorCode ierr; 1715 1716 PetscFunctionBegin; 1717 #if defined(PETSC_USE_DEBUG) 1718 if (columns || diag != 0. || (x && b)) { 1719 PetscBool cong; 1720 ierr = PetscLayoutCompare(A->rmap,A->cmap,&cong);CHKERRQ(ierr); 1721 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"); 1722 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"); 1723 if (!cong && x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"A->rmap and A->cmap need to be congruent"); 1724 } 1725 #endif 1726 /* get locally owned rows */ 1727 ierr = PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);CHKERRQ(ierr); 1728 /* fix right hand side if needed */ 1729 if (x && b) { 1730 const PetscScalar *xx; 1731 PetscScalar *bb; 1732 1733 ierr = VecGetArrayRead(x, &xx);CHKERRQ(ierr); 1734 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 1735 for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]]; 1736 ierr = VecRestoreArrayRead(x, &xx);CHKERRQ(ierr); 1737 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 1738 } 1739 /* get rows associated to the local matrices */ 1740 ierr = MatISSetUpSF(A);CHKERRQ(ierr); 1741 ierr = MatGetSize(matis->A,&nl,NULL);CHKERRQ(ierr); 1742 ierr = PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 1743 ierr = PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1744 for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1; 1745 ierr = PetscFree(lrows);CHKERRQ(ierr); 1746 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1747 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1748 ierr = PetscMalloc1(nl,&lrows);CHKERRQ(ierr); 1749 for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i; 1750 ierr = MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);CHKERRQ(ierr); 1751 ierr = PetscFree(lrows);CHKERRQ(ierr); 1752 PetscFunctionReturn(0); 1753 } 1754 1755 #undef __FUNCT__ 1756 #define __FUNCT__ "MatZeroRows_IS" 1757 static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1758 { 1759 PetscErrorCode ierr; 1760 1761 PetscFunctionBegin; 1762 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);CHKERRQ(ierr); 1763 PetscFunctionReturn(0); 1764 } 1765 1766 #undef __FUNCT__ 1767 #define __FUNCT__ "MatZeroRowsColumns_IS" 1768 static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1769 { 1770 PetscErrorCode ierr; 1771 1772 PetscFunctionBegin; 1773 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);CHKERRQ(ierr); 1774 PetscFunctionReturn(0); 1775 } 1776 1777 #undef __FUNCT__ 1778 #define __FUNCT__ "MatAssemblyBegin_IS" 1779 static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type) 1780 { 1781 Mat_IS *is = (Mat_IS*)A->data; 1782 PetscErrorCode ierr; 1783 1784 PetscFunctionBegin; 1785 ierr = MatAssemblyBegin(is->A,type);CHKERRQ(ierr); 1786 PetscFunctionReturn(0); 1787 } 1788 1789 #undef __FUNCT__ 1790 #define __FUNCT__ "MatAssemblyEnd_IS" 1791 static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type) 1792 { 1793 Mat_IS *is = (Mat_IS*)A->data; 1794 PetscErrorCode ierr; 1795 1796 PetscFunctionBegin; 1797 ierr = MatAssemblyEnd(is->A,type);CHKERRQ(ierr); 1798 PetscFunctionReturn(0); 1799 } 1800 1801 #undef __FUNCT__ 1802 #define __FUNCT__ "MatISGetLocalMat_IS" 1803 static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local) 1804 { 1805 Mat_IS *is = (Mat_IS*)mat->data; 1806 1807 PetscFunctionBegin; 1808 *local = is->A; 1809 PetscFunctionReturn(0); 1810 } 1811 1812 #undef __FUNCT__ 1813 #define __FUNCT__ "MatISGetLocalMat" 1814 /*@ 1815 MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix. 1816 1817 Input Parameter: 1818 . mat - the matrix 1819 1820 Output Parameter: 1821 . local - the local matrix 1822 1823 Level: advanced 1824 1825 Notes: 1826 This can be called if you have precomputed the nonzero structure of the 1827 matrix and want to provide it to the inner matrix object to improve the performance 1828 of the MatSetValues() operation. 1829 1830 .seealso: MATIS 1831 @*/ 1832 PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local) 1833 { 1834 PetscErrorCode ierr; 1835 1836 PetscFunctionBegin; 1837 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1838 PetscValidPointer(local,2); 1839 ierr = PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 1840 PetscFunctionReturn(0); 1841 } 1842 1843 #undef __FUNCT__ 1844 #define __FUNCT__ "MatISSetLocalMat_IS" 1845 static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local) 1846 { 1847 Mat_IS *is = (Mat_IS*)mat->data; 1848 PetscInt nrows,ncols,orows,ocols; 1849 PetscErrorCode ierr; 1850 1851 PetscFunctionBegin; 1852 if (is->A) { 1853 ierr = MatGetSize(is->A,&orows,&ocols);CHKERRQ(ierr); 1854 ierr = MatGetSize(local,&nrows,&ncols);CHKERRQ(ierr); 1855 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); 1856 } 1857 ierr = PetscObjectReference((PetscObject)local);CHKERRQ(ierr); 1858 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 1859 is->A = local; 1860 PetscFunctionReturn(0); 1861 } 1862 1863 #undef __FUNCT__ 1864 #define __FUNCT__ "MatISSetLocalMat" 1865 /*@ 1866 MatISSetLocalMat - Replace the local matrix stored inside a MATIS object. 1867 1868 Input Parameter: 1869 . mat - the matrix 1870 . local - the local matrix 1871 1872 Output Parameter: 1873 1874 Level: advanced 1875 1876 Notes: 1877 This can be called if you have precomputed the local matrix and 1878 want to provide it to the matrix object MATIS. 1879 1880 .seealso: MATIS 1881 @*/ 1882 PetscErrorCode MatISSetLocalMat(Mat mat,Mat local) 1883 { 1884 PetscErrorCode ierr; 1885 1886 PetscFunctionBegin; 1887 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1888 PetscValidHeaderSpecific(local,MAT_CLASSID,2); 1889 ierr = PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));CHKERRQ(ierr); 1890 PetscFunctionReturn(0); 1891 } 1892 1893 #undef __FUNCT__ 1894 #define __FUNCT__ "MatZeroEntries_IS" 1895 static PetscErrorCode MatZeroEntries_IS(Mat A) 1896 { 1897 Mat_IS *a = (Mat_IS*)A->data; 1898 PetscErrorCode ierr; 1899 1900 PetscFunctionBegin; 1901 ierr = MatZeroEntries(a->A);CHKERRQ(ierr); 1902 PetscFunctionReturn(0); 1903 } 1904 1905 #undef __FUNCT__ 1906 #define __FUNCT__ "MatScale_IS" 1907 static PetscErrorCode MatScale_IS(Mat A,PetscScalar a) 1908 { 1909 Mat_IS *is = (Mat_IS*)A->data; 1910 PetscErrorCode ierr; 1911 1912 PetscFunctionBegin; 1913 ierr = MatScale(is->A,a);CHKERRQ(ierr); 1914 PetscFunctionReturn(0); 1915 } 1916 1917 #undef __FUNCT__ 1918 #define __FUNCT__ "MatGetDiagonal_IS" 1919 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v) 1920 { 1921 Mat_IS *is = (Mat_IS*)A->data; 1922 PetscErrorCode ierr; 1923 1924 PetscFunctionBegin; 1925 /* get diagonal of the local matrix */ 1926 ierr = MatGetDiagonal(is->A,is->y);CHKERRQ(ierr); 1927 1928 /* scatter diagonal back into global vector */ 1929 ierr = VecSet(v,0);CHKERRQ(ierr); 1930 ierr = VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1931 ierr = VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1932 PetscFunctionReturn(0); 1933 } 1934 1935 #undef __FUNCT__ 1936 #define __FUNCT__ "MatSetOption_IS" 1937 static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg) 1938 { 1939 Mat_IS *a = (Mat_IS*)A->data; 1940 PetscErrorCode ierr; 1941 1942 PetscFunctionBegin; 1943 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 1944 PetscFunctionReturn(0); 1945 } 1946 1947 #undef __FUNCT__ 1948 #define __FUNCT__ "MatAXPY_IS" 1949 static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str) 1950 { 1951 Mat_IS *y = (Mat_IS*)Y->data; 1952 Mat_IS *x; 1953 #if defined(PETSC_USE_DEBUG) 1954 PetscBool ismatis; 1955 #endif 1956 PetscErrorCode ierr; 1957 1958 PetscFunctionBegin; 1959 #if defined(PETSC_USE_DEBUG) 1960 ierr = PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);CHKERRQ(ierr); 1961 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS"); 1962 #endif 1963 x = (Mat_IS*)X->data; 1964 ierr = MatAXPY(y->A,a,x->A,str);CHKERRQ(ierr); 1965 PetscFunctionReturn(0); 1966 } 1967 1968 #undef __FUNCT__ 1969 #define __FUNCT__ "MatGetLocalSubMatrix_IS" 1970 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat) 1971 { 1972 Mat lA; 1973 Mat_IS *matis; 1974 ISLocalToGlobalMapping rl2g,cl2g; 1975 IS is; 1976 const PetscInt *rg,*rl; 1977 PetscInt nrg; 1978 PetscInt N,M,nrl,i,*idxs; 1979 PetscErrorCode ierr; 1980 1981 PetscFunctionBegin; 1982 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 1983 ierr = ISGetLocalSize(row,&nrl);CHKERRQ(ierr); 1984 ierr = ISGetIndices(row,&rl);CHKERRQ(ierr); 1985 ierr = ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);CHKERRQ(ierr); 1986 #if defined(PETSC_USE_DEBUG) 1987 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); 1988 #endif 1989 ierr = PetscMalloc1(nrg,&idxs);CHKERRQ(ierr); 1990 /* map from [0,nrl) to row */ 1991 for (i=0;i<nrl;i++) idxs[i] = rl[i]; 1992 #if defined(PETSC_USE_DEBUG) 1993 for (i=nrl;i<nrg;i++) idxs[i] = nrg; 1994 #else 1995 for (i=nrl;i<nrg;i++) idxs[i] = -1; 1996 #endif 1997 ierr = ISRestoreIndices(row,&rl);CHKERRQ(ierr); 1998 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 1999 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 2000 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 2001 ierr = ISDestroy(&is);CHKERRQ(ierr); 2002 /* compute new l2g map for columns */ 2003 if (col != row || A->rmap->mapping != A->cmap->mapping) { 2004 const PetscInt *cg,*cl; 2005 PetscInt ncg; 2006 PetscInt ncl; 2007 2008 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 2009 ierr = ISGetLocalSize(col,&ncl);CHKERRQ(ierr); 2010 ierr = ISGetIndices(col,&cl);CHKERRQ(ierr); 2011 ierr = ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);CHKERRQ(ierr); 2012 #if defined(PETSC_USE_DEBUG) 2013 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); 2014 #endif 2015 ierr = PetscMalloc1(ncg,&idxs);CHKERRQ(ierr); 2016 /* map from [0,ncl) to col */ 2017 for (i=0;i<ncl;i++) idxs[i] = cl[i]; 2018 #if defined(PETSC_USE_DEBUG) 2019 for (i=ncl;i<ncg;i++) idxs[i] = ncg; 2020 #else 2021 for (i=ncl;i<ncg;i++) idxs[i] = -1; 2022 #endif 2023 ierr = ISRestoreIndices(col,&cl);CHKERRQ(ierr); 2024 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 2025 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 2026 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 2027 ierr = ISDestroy(&is);CHKERRQ(ierr); 2028 } else { 2029 ierr = PetscObjectReference((PetscObject)rl2g);CHKERRQ(ierr); 2030 cl2g = rl2g; 2031 } 2032 /* create the MATIS submatrix */ 2033 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 2034 ierr = MatCreate(PetscObjectComm((PetscObject)A),submat);CHKERRQ(ierr); 2035 ierr = MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2036 ierr = MatSetType(*submat,MATIS);CHKERRQ(ierr); 2037 matis = (Mat_IS*)((*submat)->data); 2038 matis->islocalref = PETSC_TRUE; 2039 ierr = MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);CHKERRQ(ierr); 2040 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr); 2041 ierr = MatISSetLocalMat(*submat,lA);CHKERRQ(ierr); 2042 ierr = MatSetUp(*submat);CHKERRQ(ierr); 2043 ierr = MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2044 ierr = MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2045 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 2046 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 2047 /* remove unsupported ops */ 2048 ierr = PetscMemzero((*submat)->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2049 (*submat)->ops->destroy = MatDestroy_IS; 2050 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS; 2051 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS; 2052 (*submat)->ops->assemblybegin = MatAssemblyBegin_IS; 2053 (*submat)->ops->assemblyend = MatAssemblyEnd_IS; 2054 PetscFunctionReturn(0); 2055 } 2056 2057 #undef __FUNCT__ 2058 #define __FUNCT__ "MatCreateIS" 2059 /*@ 2060 MatCreateIS - Creates a "process" unassembled matrix, assembled on each 2061 process but not across processes. 2062 2063 Input Parameters: 2064 + comm - MPI communicator that will share the matrix 2065 . bs - block size of the matrix 2066 . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products 2067 . rmap - local to global map for rows 2068 - cmap - local to global map for cols 2069 2070 Output Parameter: 2071 . A - the resulting matrix 2072 2073 Level: advanced 2074 2075 Notes: See MATIS for more details. 2076 m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors 2077 used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices. 2078 If either rmap or cmap are NULL, then the matrix is assumed to be square. 2079 2080 .seealso: MATIS, MatSetLocalToGlobalMapping() 2081 @*/ 2082 PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A) 2083 { 2084 PetscErrorCode ierr; 2085 2086 PetscFunctionBegin; 2087 if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings"); 2088 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2089 if (bs > 0) { 2090 ierr = MatSetBlockSize(*A,bs);CHKERRQ(ierr); 2091 } 2092 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 2093 ierr = MatSetType(*A,MATIS);CHKERRQ(ierr); 2094 if (rmap && cmap) { 2095 ierr = MatSetLocalToGlobalMapping(*A,rmap,cmap);CHKERRQ(ierr); 2096 } else if (!rmap) { 2097 ierr = MatSetLocalToGlobalMapping(*A,cmap,cmap);CHKERRQ(ierr); 2098 } else { 2099 ierr = MatSetLocalToGlobalMapping(*A,rmap,rmap);CHKERRQ(ierr); 2100 } 2101 PetscFunctionReturn(0); 2102 } 2103 2104 /*MC 2105 MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP). 2106 This stores the matrices in globally unassembled form. Each processor 2107 assembles only its local Neumann problem and the parallel matrix vector 2108 product is handled "implicitly". 2109 2110 Operations Provided: 2111 + MatMult() 2112 . MatMultAdd() 2113 . MatMultTranspose() 2114 . MatMultTransposeAdd() 2115 . MatZeroEntries() 2116 . MatSetOption() 2117 . MatZeroRows() 2118 . MatSetValues() 2119 . MatSetValuesBlocked() 2120 . MatSetValuesLocal() 2121 . MatSetValuesBlockedLocal() 2122 . MatScale() 2123 . MatGetDiagonal() 2124 . MatMissingDiagonal() 2125 . MatDuplicate() 2126 . MatCopy() 2127 . MatAXPY() 2128 . MatGetSubMatrix() 2129 . MatGetLocalSubMatrix() 2130 . MatTranspose() 2131 - MatSetLocalToGlobalMapping() 2132 2133 Options Database Keys: 2134 . -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions() 2135 2136 Notes: Options prefix for the inner matrix are given by -is_mat_xxx 2137 2138 You must call MatSetLocalToGlobalMapping() before using this matrix type. 2139 2140 You can do matrix preallocation on the local matrix after you obtain it with 2141 MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation() 2142 2143 Level: advanced 2144 2145 .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP 2146 2147 M*/ 2148 2149 #undef __FUNCT__ 2150 #define __FUNCT__ "MatCreate_IS" 2151 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A) 2152 { 2153 PetscErrorCode ierr; 2154 Mat_IS *b; 2155 2156 PetscFunctionBegin; 2157 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 2158 A->data = (void*)b; 2159 2160 /* matrix ops */ 2161 ierr = PetscMemzero(A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2162 A->ops->mult = MatMult_IS; 2163 A->ops->multadd = MatMultAdd_IS; 2164 A->ops->multtranspose = MatMultTranspose_IS; 2165 A->ops->multtransposeadd = MatMultTransposeAdd_IS; 2166 A->ops->destroy = MatDestroy_IS; 2167 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS; 2168 A->ops->setvalues = MatSetValues_IS; 2169 A->ops->setvaluesblocked = MatSetValuesBlocked_IS; 2170 A->ops->setvalueslocal = MatSetValuesLocal_IS; 2171 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS; 2172 A->ops->zerorows = MatZeroRows_IS; 2173 A->ops->zerorowscolumns = MatZeroRowsColumns_IS; 2174 A->ops->assemblybegin = MatAssemblyBegin_IS; 2175 A->ops->assemblyend = MatAssemblyEnd_IS; 2176 A->ops->view = MatView_IS; 2177 A->ops->zeroentries = MatZeroEntries_IS; 2178 A->ops->scale = MatScale_IS; 2179 A->ops->getdiagonal = MatGetDiagonal_IS; 2180 A->ops->setoption = MatSetOption_IS; 2181 A->ops->ishermitian = MatIsHermitian_IS; 2182 A->ops->issymmetric = MatIsSymmetric_IS; 2183 A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS; 2184 A->ops->duplicate = MatDuplicate_IS; 2185 A->ops->missingdiagonal = MatMissingDiagonal_IS; 2186 A->ops->copy = MatCopy_IS; 2187 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 2188 A->ops->getsubmatrix = MatGetSubMatrix_IS; 2189 A->ops->axpy = MatAXPY_IS; 2190 A->ops->diagonalset = MatDiagonalSet_IS; 2191 A->ops->shift = MatShift_IS; 2192 A->ops->transpose = MatTranspose_IS; 2193 2194 /* special MATIS functions */ 2195 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);CHKERRQ(ierr); 2196 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);CHKERRQ(ierr); 2197 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatISGetMPIXAIJ_IS);CHKERRQ(ierr); 2198 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);CHKERRQ(ierr); 2199 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",MatISSetUpSF_IS);CHKERRQ(ierr); 2200 ierr = PetscObjectChangeTypeName((PetscObject)A,MATIS);CHKERRQ(ierr); 2201 PetscFunctionReturn(0); 2202 } 2203