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