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