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