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