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