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