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