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