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