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