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