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 PetscErrorCode ierr; 615 616 PetscFunctionBegin; 617 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) { 618 ISLocalToGlobalMapping rl2g,cl2g; 619 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 620 ierr = MatSetSizes(C,A->cmap->n,A->rmap->n,A->cmap->N,A->rmap->N);CHKERRQ(ierr); 621 ierr = MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 622 ierr = MatSetType(C,MATIS);CHKERRQ(ierr); 623 ierr = MatGetLocalToGlobalMapping(A,&rl2g,&cl2g);CHKERRQ(ierr); 624 ierr = MatSetLocalToGlobalMapping(C,cl2g,rl2g);CHKERRQ(ierr); 625 } else { 626 C = *B; 627 } 628 629 /* perform local transposition */ 630 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr); 631 ierr = MatTranspose(lA,MAT_INITIAL_MATRIX,&lC);CHKERRQ(ierr); 632 ierr = MatISSetLocalMat(C,lC);CHKERRQ(ierr); 633 ierr = MatDestroy(&lC);CHKERRQ(ierr); 634 635 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) { 636 *B = C; 637 } else { 638 ierr = MatHeaderMerge(A,&C);CHKERRQ(ierr); 639 } 640 PetscFunctionReturn(0); 641 } 642 643 #undef __FUNCT__ 644 #define __FUNCT__ "MatDiagonalSet_IS" 645 PetscErrorCode MatDiagonalSet_IS(Mat A,Vec D,InsertMode insmode) 646 { 647 Mat_IS *is = (Mat_IS*)A->data; 648 PetscErrorCode ierr; 649 650 PetscFunctionBegin; 651 if (!D) { /* this code branch is used by MatShift_IS */ 652 ierr = VecSet(is->y,1.);CHKERRQ(ierr); 653 } else { 654 ierr = VecScatterBegin(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 655 ierr = VecScatterEnd(is->rctx,D,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 656 } 657 ierr = VecPointwiseDivide(is->y,is->y,is->counter);CHKERRQ(ierr); 658 ierr = MatDiagonalSet(is->A,is->y,insmode);CHKERRQ(ierr); 659 PetscFunctionReturn(0); 660 } 661 662 #undef __FUNCT__ 663 #define __FUNCT__ "MatShift_IS" 664 PetscErrorCode MatShift_IS(Mat A,PetscScalar a) 665 { 666 PetscErrorCode ierr; 667 668 PetscFunctionBegin; 669 ierr = MatDiagonalSet_IS(A,NULL,ADD_VALUES);CHKERRQ(ierr); 670 PetscFunctionReturn(0); 671 } 672 673 #undef __FUNCT__ 674 #define __FUNCT__ "MatSetValuesLocal_SubMat_IS" 675 static PetscErrorCode MatSetValuesLocal_SubMat_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 676 { 677 PetscErrorCode ierr; 678 Mat_IS *is = (Mat_IS*)A->data; 679 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 680 681 PetscFunctionBegin; 682 #if defined(PETSC_USE_DEBUG) 683 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); 684 #endif 685 ierr = ISLocalToGlobalMappingApply(A->rmap->mapping,m,rows,rows_l);CHKERRQ(ierr); 686 ierr = ISLocalToGlobalMappingApply(A->cmap->mapping,n,cols,cols_l);CHKERRQ(ierr); 687 ierr = MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 688 PetscFunctionReturn(0); 689 } 690 691 #undef __FUNCT__ 692 #define __FUNCT__ "MatSetValuesBlockedLocal_SubMat_IS" 693 static PetscErrorCode MatSetValuesBlockedLocal_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 block indices must be <= %D: they are %D %D",MATIS_MAX_ENTRIES_INSERTION,m,n); 702 #endif 703 ierr = ISLocalToGlobalMappingApplyBlock(A->rmap->mapping,m,rows,rows_l);CHKERRQ(ierr); 704 ierr = ISLocalToGlobalMappingApplyBlock(A->cmap->mapping,n,cols,cols_l);CHKERRQ(ierr); 705 ierr = MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 706 PetscFunctionReturn(0); 707 } 708 709 #undef __FUNCT__ 710 #define __FUNCT__ "PetscLayoutMapLocal_Private" 711 static PetscErrorCode PetscLayoutMapLocal_Private(PetscLayout map,PetscInt N,const PetscInt idxs[], PetscInt *on,PetscInt **oidxs,PetscInt **ogidxs) 712 { 713 PetscInt *owners = map->range; 714 PetscInt n = map->n; 715 PetscSF sf; 716 PetscInt *lidxs,*work = NULL; 717 PetscSFNode *ridxs; 718 PetscMPIInt rank; 719 PetscInt r, p = 0, len = 0; 720 PetscErrorCode ierr; 721 722 PetscFunctionBegin; 723 if (on) *on = 0; /* squelch -Wmaybe-uninitialized */ 724 /* Create SF where leaves are input idxs and roots are owned idxs (code adapted from MatZeroRowsMapLocal_Private) */ 725 ierr = MPI_Comm_rank(map->comm,&rank);CHKERRQ(ierr); 726 ierr = PetscMalloc1(n,&lidxs);CHKERRQ(ierr); 727 for (r = 0; r < n; ++r) lidxs[r] = -1; 728 ierr = PetscMalloc1(N,&ridxs);CHKERRQ(ierr); 729 for (r = 0; r < N; ++r) { 730 const PetscInt idx = idxs[r]; 731 if (idx < 0 || map->N <= idx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index %D out of range [0,%D)",idx,map->N); 732 if (idx < owners[p] || owners[p+1] <= idx) { /* short-circuit the search if the last p owns this idx too */ 733 ierr = PetscLayoutFindOwner(map,idx,&p);CHKERRQ(ierr); 734 } 735 ridxs[r].rank = p; 736 ridxs[r].index = idxs[r] - owners[p]; 737 } 738 ierr = PetscSFCreate(map->comm,&sf);CHKERRQ(ierr); 739 ierr = PetscSFSetGraph(sf,n,N,NULL,PETSC_OWN_POINTER,ridxs,PETSC_OWN_POINTER);CHKERRQ(ierr); 740 ierr = PetscSFReduceBegin(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);CHKERRQ(ierr); 741 ierr = PetscSFReduceEnd(sf,MPIU_INT,(PetscInt*)idxs,lidxs,MPI_LOR);CHKERRQ(ierr); 742 if (ogidxs) { /* communicate global idxs */ 743 PetscInt cum = 0,start,*work2; 744 745 ierr = PetscMalloc1(n,&work);CHKERRQ(ierr); 746 ierr = PetscCalloc1(N,&work2);CHKERRQ(ierr); 747 for (r = 0; r < N; ++r) if (idxs[r] >=0) cum++; 748 ierr = MPI_Scan(&cum,&start,1,MPIU_INT,MPI_SUM,map->comm);CHKERRQ(ierr); 749 start -= cum; 750 cum = 0; 751 for (r = 0; r < N; ++r) if (idxs[r] >=0) work2[r] = start+cum++; 752 ierr = PetscSFReduceBegin(sf,MPIU_INT,work2,work,MPIU_REPLACE);CHKERRQ(ierr); 753 ierr = PetscSFReduceEnd(sf,MPIU_INT,work2,work,MPIU_REPLACE);CHKERRQ(ierr); 754 ierr = PetscFree(work2);CHKERRQ(ierr); 755 } 756 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 757 /* Compress and put in indices */ 758 for (r = 0; r < n; ++r) 759 if (lidxs[r] >= 0) { 760 if (work) work[len] = work[r]; 761 lidxs[len++] = r; 762 } 763 if (on) *on = len; 764 if (oidxs) *oidxs = lidxs; 765 if (ogidxs) *ogidxs = work; 766 PetscFunctionReturn(0); 767 } 768 769 #undef __FUNCT__ 770 #define __FUNCT__ "MatGetSubMatrix_IS" 771 static PetscErrorCode MatGetSubMatrix_IS(Mat mat,IS irow,IS icol,MatReuse scall,Mat *newmat) 772 { 773 Mat locmat,newlocmat; 774 Mat_IS *newmatis; 775 #if defined(PETSC_USE_DEBUG) 776 Vec rtest,ltest; 777 const PetscScalar *array; 778 #endif 779 const PetscInt *idxs; 780 PetscInt i,m,n; 781 PetscErrorCode ierr; 782 783 PetscFunctionBegin; 784 if (scall == MAT_REUSE_MATRIX) { 785 PetscBool ismatis; 786 787 ierr = PetscObjectTypeCompare((PetscObject)*newmat,MATIS,&ismatis);CHKERRQ(ierr); 788 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Not of MATIS type"); 789 newmatis = (Mat_IS*)(*newmat)->data; 790 if (!newmatis->getsub_ris) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local row IS"); 791 if (!newmatis->getsub_cis) SETERRQ(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_ARG_WRONG,"Cannot reuse matrix! Misses local col IS"); 792 } 793 /* irow and icol may not have duplicate entries */ 794 #if defined(PETSC_USE_DEBUG) 795 ierr = MatCreateVecs(mat,<est,&rtest);CHKERRQ(ierr); 796 ierr = ISGetLocalSize(irow,&n);CHKERRQ(ierr); 797 ierr = ISGetIndices(irow,&idxs);CHKERRQ(ierr); 798 for (i=0;i<n;i++) { 799 ierr = VecSetValue(rtest,idxs[i],1.0,ADD_VALUES);CHKERRQ(ierr); 800 } 801 ierr = VecAssemblyBegin(rtest);CHKERRQ(ierr); 802 ierr = VecAssemblyEnd(rtest);CHKERRQ(ierr); 803 ierr = VecGetLocalSize(rtest,&n);CHKERRQ(ierr); 804 ierr = VecGetOwnershipRange(rtest,&m,NULL);CHKERRQ(ierr); 805 ierr = VecGetArrayRead(rtest,&array);CHKERRQ(ierr); 806 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])); 807 ierr = VecRestoreArrayRead(rtest,&array);CHKERRQ(ierr); 808 ierr = ISRestoreIndices(irow,&idxs);CHKERRQ(ierr); 809 ierr = ISGetLocalSize(icol,&n);CHKERRQ(ierr); 810 ierr = ISGetIndices(icol,&idxs);CHKERRQ(ierr); 811 for (i=0;i<n;i++) { 812 ierr = VecSetValue(ltest,idxs[i],1.0,ADD_VALUES);CHKERRQ(ierr); 813 } 814 ierr = VecAssemblyBegin(ltest);CHKERRQ(ierr); 815 ierr = VecAssemblyEnd(ltest);CHKERRQ(ierr); 816 ierr = VecGetLocalSize(ltest,&n);CHKERRQ(ierr); 817 ierr = VecGetOwnershipRange(ltest,&m,NULL);CHKERRQ(ierr); 818 ierr = VecGetArrayRead(ltest,&array);CHKERRQ(ierr); 819 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])); 820 ierr = VecRestoreArrayRead(ltest,&array);CHKERRQ(ierr); 821 ierr = ISRestoreIndices(icol,&idxs);CHKERRQ(ierr); 822 ierr = VecDestroy(&rtest);CHKERRQ(ierr); 823 ierr = VecDestroy(<est);CHKERRQ(ierr); 824 #endif 825 if (scall == MAT_INITIAL_MATRIX) { 826 Mat_IS *matis = (Mat_IS*)mat->data; 827 ISLocalToGlobalMapping rl2g; 828 IS is; 829 PetscInt *lidxs,*lgidxs,*newgidxs; 830 PetscInt ll,newloc; 831 MPI_Comm comm; 832 833 ierr = PetscObjectGetComm((PetscObject)mat,&comm);CHKERRQ(ierr); 834 ierr = ISGetLocalSize(irow,&m);CHKERRQ(ierr); 835 ierr = ISGetLocalSize(icol,&n);CHKERRQ(ierr); 836 ierr = MatCreate(comm,newmat);CHKERRQ(ierr); 837 ierr = MatSetType(*newmat,MATIS);CHKERRQ(ierr); 838 ierr = MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 839 /* communicate irow to their owners in the layout */ 840 ierr = ISGetIndices(irow,&idxs);CHKERRQ(ierr); 841 ierr = PetscLayoutMapLocal_Private(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs);CHKERRQ(ierr); 842 ierr = ISRestoreIndices(irow,&idxs);CHKERRQ(ierr); 843 ierr = MatISSetUpSF(mat);CHKERRQ(ierr); 844 ierr = PetscMemzero(matis->sf_rootdata,matis->sf->nroots*sizeof(PetscInt));CHKERRQ(ierr); 845 for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1; 846 ierr = PetscFree(lidxs);CHKERRQ(ierr); 847 ierr = PetscFree(lgidxs);CHKERRQ(ierr); 848 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 849 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 850 for (i=0,newloc=0;i<matis->sf->nleaves;i++) if (matis->sf_leafdata[i]) newloc++; 851 ierr = PetscMalloc1(newloc,&newgidxs);CHKERRQ(ierr); 852 ierr = PetscMalloc1(newloc,&lidxs);CHKERRQ(ierr); 853 for (i=0,newloc=0;i<matis->sf->nleaves;i++) 854 if (matis->sf_leafdata[i]) { 855 lidxs[newloc] = i; 856 newgidxs[newloc++] = matis->sf_leafdata[i]-1; 857 } 858 ierr = ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 859 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 860 ierr = ISDestroy(&is);CHKERRQ(ierr); 861 /* local is to extract local submatrix */ 862 newmatis = (Mat_IS*)(*newmat)->data; 863 ierr = ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris);CHKERRQ(ierr); 864 if (mat->congruentlayouts == -1) { /* first time we compare rows and cols layouts */ 865 PetscBool cong; 866 ierr = PetscLayoutCompare(mat->rmap,mat->cmap,&cong);CHKERRQ(ierr); 867 if (cong) mat->congruentlayouts = 1; 868 else mat->congruentlayouts = 0; 869 } 870 if (mat->congruentlayouts && irow == icol) { 871 ierr = MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g);CHKERRQ(ierr); 872 ierr = PetscObjectReference((PetscObject)newmatis->getsub_ris);CHKERRQ(ierr); 873 newmatis->getsub_cis = newmatis->getsub_ris; 874 } else { 875 ISLocalToGlobalMapping cl2g; 876 877 /* communicate icol to their owners in the layout */ 878 ierr = ISGetIndices(icol,&idxs);CHKERRQ(ierr); 879 ierr = PetscLayoutMapLocal_Private(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs);CHKERRQ(ierr); 880 ierr = ISRestoreIndices(icol,&idxs);CHKERRQ(ierr); 881 ierr = PetscMemzero(matis->csf_rootdata,matis->csf->nroots*sizeof(PetscInt));CHKERRQ(ierr); 882 for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1; 883 ierr = PetscFree(lidxs);CHKERRQ(ierr); 884 ierr = PetscFree(lgidxs);CHKERRQ(ierr); 885 ierr = PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);CHKERRQ(ierr); 886 ierr = PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);CHKERRQ(ierr); 887 for (i=0,newloc=0;i<matis->csf->nleaves;i++) if (matis->csf_leafdata[i]) newloc++; 888 ierr = PetscMalloc1(newloc,&newgidxs);CHKERRQ(ierr); 889 ierr = PetscMalloc1(newloc,&lidxs);CHKERRQ(ierr); 890 for (i=0,newloc=0;i<matis->csf->nleaves;i++) 891 if (matis->csf_leafdata[i]) { 892 lidxs[newloc] = i; 893 newgidxs[newloc++] = matis->csf_leafdata[i]-1; 894 } 895 ierr = ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 896 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 897 ierr = ISDestroy(&is);CHKERRQ(ierr); 898 /* local is to extract local submatrix */ 899 ierr = ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis);CHKERRQ(ierr); 900 ierr = MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);CHKERRQ(ierr); 901 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 902 } 903 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 904 } else { 905 ierr = MatISGetLocalMat(*newmat,&newlocmat);CHKERRQ(ierr); 906 } 907 ierr = MatISGetLocalMat(mat,&locmat);CHKERRQ(ierr); 908 newmatis = (Mat_IS*)(*newmat)->data; 909 ierr = MatGetSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat);CHKERRQ(ierr); 910 if (scall == MAT_INITIAL_MATRIX) { 911 ierr = MatISSetLocalMat(*newmat,newlocmat);CHKERRQ(ierr); 912 ierr = MatDestroy(&newlocmat);CHKERRQ(ierr); 913 } 914 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 915 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 916 PetscFunctionReturn(0); 917 } 918 919 #undef __FUNCT__ 920 #define __FUNCT__ "MatCopy_IS" 921 static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str) 922 { 923 Mat_IS *a = (Mat_IS*)A->data,*b; 924 PetscBool ismatis; 925 PetscErrorCode ierr; 926 927 PetscFunctionBegin; 928 ierr = PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);CHKERRQ(ierr); 929 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Need to be implemented"); 930 b = (Mat_IS*)B->data; 931 ierr = MatCopy(a->A,b->A,str);CHKERRQ(ierr); 932 PetscFunctionReturn(0); 933 } 934 935 #undef __FUNCT__ 936 #define __FUNCT__ "MatMissingDiagonal_IS" 937 static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool *missing,PetscInt *d) 938 { 939 Vec v; 940 const PetscScalar *array; 941 PetscInt i,n; 942 PetscErrorCode ierr; 943 944 PetscFunctionBegin; 945 *missing = PETSC_FALSE; 946 ierr = MatCreateVecs(A,NULL,&v);CHKERRQ(ierr); 947 ierr = MatGetDiagonal(A,v);CHKERRQ(ierr); 948 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 949 ierr = VecGetArrayRead(v,&array);CHKERRQ(ierr); 950 for (i=0;i<n;i++) if (array[i] == 0.) break; 951 ierr = VecRestoreArrayRead(v,&array);CHKERRQ(ierr); 952 ierr = VecDestroy(&v);CHKERRQ(ierr); 953 if (i != n) *missing = PETSC_TRUE; 954 if (d) { 955 *d = -1; 956 if (*missing) { 957 PetscInt rstart; 958 ierr = MatGetOwnershipRange(A,&rstart,NULL);CHKERRQ(ierr); 959 *d = i+rstart; 960 } 961 } 962 PetscFunctionReturn(0); 963 } 964 965 #undef __FUNCT__ 966 #define __FUNCT__ "MatISSetUpSF_IS" 967 static PetscErrorCode MatISSetUpSF_IS(Mat B) 968 { 969 Mat_IS *matis = (Mat_IS*)(B->data); 970 const PetscInt *gidxs; 971 PetscInt nleaves; 972 PetscErrorCode ierr; 973 974 PetscFunctionBegin; 975 if (matis->sf) PetscFunctionReturn(0); 976 ierr = PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);CHKERRQ(ierr); 977 ierr = ISLocalToGlobalMappingGetIndices(B->rmap->mapping,&gidxs);CHKERRQ(ierr); 978 ierr = ISLocalToGlobalMappingGetSize(B->rmap->mapping,&nleaves);CHKERRQ(ierr); 979 ierr = PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);CHKERRQ(ierr); 980 ierr = ISLocalToGlobalMappingRestoreIndices(B->rmap->mapping,&gidxs);CHKERRQ(ierr); 981 ierr = PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata);CHKERRQ(ierr); 982 if (B->rmap->mapping != B->cmap->mapping) { /* setup SF for columns */ 983 ierr = ISLocalToGlobalMappingGetSize(B->cmap->mapping,&nleaves);CHKERRQ(ierr); 984 ierr = PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);CHKERRQ(ierr); 985 ierr = ISLocalToGlobalMappingGetIndices(B->cmap->mapping,&gidxs);CHKERRQ(ierr); 986 ierr = PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);CHKERRQ(ierr); 987 ierr = ISLocalToGlobalMappingRestoreIndices(B->cmap->mapping,&gidxs);CHKERRQ(ierr); 988 ierr = PetscMalloc2(matis->csf->nroots,&matis->csf_rootdata,matis->csf->nleaves,&matis->csf_leafdata);CHKERRQ(ierr); 989 } else { 990 matis->csf = matis->sf; 991 matis->csf_leafdata = matis->sf_leafdata; 992 matis->csf_rootdata = matis->sf_rootdata; 993 } 994 PetscFunctionReturn(0); 995 } 996 997 #undef __FUNCT__ 998 #define __FUNCT__ "MatISSetPreallocation" 999 /*@ 1000 MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix. 1001 1002 Collective on MPI_Comm 1003 1004 Input Parameters: 1005 + B - the matrix 1006 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1007 (same value is used for all local rows) 1008 . d_nnz - array containing the number of nonzeros in the various rows of the 1009 DIAGONAL portion of the local submatrix (possibly different for each row) 1010 or NULL, if d_nz is used to specify the nonzero structure. 1011 The size of this array is equal to the number of local rows, i.e 'm'. 1012 For matrices that will be factored, you must leave room for (and set) 1013 the diagonal entry even if it is zero. 1014 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1015 submatrix (same value is used for all local rows). 1016 - o_nnz - array containing the number of nonzeros in the various rows of the 1017 OFF-DIAGONAL portion of the local submatrix (possibly different for 1018 each row) or NULL, if o_nz is used to specify the nonzero 1019 structure. The size of this array is equal to the number 1020 of local rows, i.e 'm'. 1021 1022 If the *_nnz parameter is given then the *_nz parameter is ignored 1023 1024 Level: intermediate 1025 1026 Notes: This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition 1027 from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local 1028 matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects. 1029 1030 .keywords: matrix 1031 1032 .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS 1033 @*/ 1034 PetscErrorCode MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1035 { 1036 PetscErrorCode ierr; 1037 1038 PetscFunctionBegin; 1039 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1040 PetscValidType(B,1); 1041 ierr = PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));CHKERRQ(ierr); 1042 PetscFunctionReturn(0); 1043 } 1044 1045 #undef __FUNCT__ 1046 #define __FUNCT__ "MatISSetPreallocation_IS" 1047 static PetscErrorCode MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1048 { 1049 Mat_IS *matis = (Mat_IS*)(B->data); 1050 PetscInt bs,i,nlocalcols; 1051 PetscErrorCode ierr; 1052 1053 PetscFunctionBegin; 1054 if (!matis->A) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"You should first call MatSetLocalToGlobalMapping"); 1055 ierr = MatISSetUpSF(B);CHKERRQ(ierr); 1056 1057 if (!d_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nz; 1058 else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nnz[i]; 1059 1060 if (!o_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nz; 1061 else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nnz[i]; 1062 1063 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1064 ierr = MatGetSize(matis->A,NULL,&nlocalcols);CHKERRQ(ierr); 1065 ierr = MatGetBlockSize(matis->A,&bs);CHKERRQ(ierr); 1066 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1067 1068 for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols); 1069 ierr = MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);CHKERRQ(ierr); 1070 1071 for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs; 1072 ierr = MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);CHKERRQ(ierr); 1073 1074 for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i]-i; 1075 ierr = MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);CHKERRQ(ierr); 1076 PetscFunctionReturn(0); 1077 } 1078 1079 #undef __FUNCT__ 1080 #define __FUNCT__ "MatISSetMPIXAIJPreallocation_Private" 1081 PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce) 1082 { 1083 Mat_IS *matis = (Mat_IS*)(A->data); 1084 PetscInt *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership; 1085 const PetscInt *global_indices_r,*global_indices_c; 1086 PetscInt i,j,bs,rows,cols; 1087 PetscInt lrows,lcols; 1088 PetscInt local_rows,local_cols; 1089 PetscMPIInt nsubdomains; 1090 PetscBool isdense,issbaij; 1091 PetscErrorCode ierr; 1092 1093 PetscFunctionBegin; 1094 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&nsubdomains);CHKERRQ(ierr); 1095 ierr = MatGetSize(A,&rows,&cols);CHKERRQ(ierr); 1096 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 1097 ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr); 1098 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);CHKERRQ(ierr); 1099 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);CHKERRQ(ierr); 1100 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);CHKERRQ(ierr); 1101 if (A->rmap->mapping != A->cmap->mapping) { 1102 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&global_indices_c);CHKERRQ(ierr); 1103 } else { 1104 global_indices_c = global_indices_r; 1105 } 1106 1107 if (issbaij) { 1108 ierr = MatGetRowUpperTriangular(matis->A);CHKERRQ(ierr); 1109 } 1110 /* 1111 An SF reduce is needed to sum up properly on shared rows. 1112 Note that generally preallocation is not exact, since it overestimates nonzeros 1113 */ 1114 ierr = MatISSetUpSF(A);CHKERRQ(ierr); 1115 ierr = MatGetLocalSize(A,&lrows,&lcols);CHKERRQ(ierr); 1116 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);CHKERRQ(ierr); 1117 /* All processes need to compute entire row ownership */ 1118 ierr = PetscMalloc1(rows,&row_ownership);CHKERRQ(ierr); 1119 ierr = MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);CHKERRQ(ierr); 1120 for (i=0;i<nsubdomains;i++) { 1121 for (j=mat_ranges[i];j<mat_ranges[i+1];j++) { 1122 row_ownership[j] = i; 1123 } 1124 } 1125 ierr = MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);CHKERRQ(ierr); 1126 1127 /* 1128 my_dnz and my_onz contains exact contribution to preallocation from each local mat 1129 then, they will be summed up properly. This way, preallocation is always sufficient 1130 */ 1131 ierr = PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);CHKERRQ(ierr); 1132 /* preallocation as a MATAIJ */ 1133 if (isdense) { /* special case for dense local matrices */ 1134 for (i=0;i<local_rows;i++) { 1135 PetscInt index_row = global_indices_r[i]; 1136 for (j=i;j<local_rows;j++) { 1137 PetscInt owner = row_ownership[index_row]; 1138 PetscInt index_col = global_indices_c[j]; 1139 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 1140 my_dnz[i] += 1; 1141 } else { /* offdiag block */ 1142 my_onz[i] += 1; 1143 } 1144 /* same as before, interchanging rows and cols */ 1145 if (i != j) { 1146 owner = row_ownership[index_col]; 1147 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) { 1148 my_dnz[j] += 1; 1149 } else { 1150 my_onz[j] += 1; 1151 } 1152 } 1153 } 1154 } 1155 } else if (matis->A->ops->getrowij) { 1156 const PetscInt *ii,*jj,*jptr; 1157 PetscBool done; 1158 ierr = MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);CHKERRQ(ierr); 1159 if (!done) SETERRQ1(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ called from %s\n",__FUNCT__); 1160 jptr = jj; 1161 for (i=0;i<local_rows;i++) { 1162 PetscInt index_row = global_indices_r[i]; 1163 for (j=0;j<ii[i+1]-ii[i];j++,jptr++) { 1164 PetscInt owner = row_ownership[index_row]; 1165 PetscInt index_col = global_indices_c[*jptr]; 1166 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 1167 my_dnz[i] += 1; 1168 } else { /* offdiag block */ 1169 my_onz[i] += 1; 1170 } 1171 /* same as before, interchanging rows and cols */ 1172 if (issbaij && index_col != index_row) { 1173 owner = row_ownership[index_col]; 1174 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) { 1175 my_dnz[*jptr] += 1; 1176 } else { 1177 my_onz[*jptr] += 1; 1178 } 1179 } 1180 } 1181 } 1182 ierr = MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);CHKERRQ(ierr); 1183 if (!done) SETERRQ1(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called from %s\n",__FUNCT__); 1184 } else { /* loop over rows and use MatGetRow */ 1185 for (i=0;i<local_rows;i++) { 1186 const PetscInt *cols; 1187 PetscInt ncols,index_row = global_indices_r[i]; 1188 ierr = MatGetRow(matis->A,i,&ncols,&cols,NULL);CHKERRQ(ierr); 1189 for (j=0;j<ncols;j++) { 1190 PetscInt owner = row_ownership[index_row]; 1191 PetscInt index_col = global_indices_c[cols[j]]; 1192 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 1193 my_dnz[i] += 1; 1194 } else { /* offdiag block */ 1195 my_onz[i] += 1; 1196 } 1197 /* same as before, interchanging rows and cols */ 1198 if (issbaij && index_col != index_row) { 1199 owner = row_ownership[index_col]; 1200 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) { 1201 my_dnz[cols[j]] += 1; 1202 } else { 1203 my_onz[cols[j]] += 1; 1204 } 1205 } 1206 } 1207 ierr = MatRestoreRow(matis->A,i,&ncols,&cols,NULL);CHKERRQ(ierr); 1208 } 1209 } 1210 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);CHKERRQ(ierr); 1211 if (global_indices_c != global_indices_r) { 1212 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&global_indices_c);CHKERRQ(ierr); 1213 } 1214 ierr = PetscFree(row_ownership);CHKERRQ(ierr); 1215 1216 /* Reduce my_dnz and my_onz */ 1217 if (maxreduce) { 1218 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);CHKERRQ(ierr); 1219 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);CHKERRQ(ierr); 1220 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);CHKERRQ(ierr); 1221 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);CHKERRQ(ierr); 1222 } else { 1223 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);CHKERRQ(ierr); 1224 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);CHKERRQ(ierr); 1225 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);CHKERRQ(ierr); 1226 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);CHKERRQ(ierr); 1227 } 1228 ierr = PetscFree2(my_dnz,my_onz);CHKERRQ(ierr); 1229 1230 /* Resize preallocation if overestimated */ 1231 for (i=0;i<lrows;i++) { 1232 dnz[i] = PetscMin(dnz[i],lcols); 1233 onz[i] = PetscMin(onz[i],cols-lcols); 1234 } 1235 1236 /* Set preallocation */ 1237 ierr = MatMPIAIJSetPreallocation(B,0,dnz,0,onz);CHKERRQ(ierr); 1238 for (i=0;i<lrows/bs;i++) { 1239 dnz[i] = dnz[i*bs]/bs; 1240 onz[i] = onz[i*bs]/bs; 1241 } 1242 ierr = MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 1243 ierr = MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 1244 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1245 if (issbaij) { 1246 ierr = MatRestoreRowUpperTriangular(matis->A);CHKERRQ(ierr); 1247 } 1248 PetscFunctionReturn(0); 1249 } 1250 1251 #undef __FUNCT__ 1252 #define __FUNCT__ "MatISGetMPIXAIJ_IS" 1253 static PetscErrorCode MatISGetMPIXAIJ_IS(Mat mat, MatReuse reuse, Mat *M) 1254 { 1255 Mat_IS *matis = (Mat_IS*)(mat->data); 1256 Mat local_mat; 1257 /* info on mat */ 1258 PetscInt bs,rows,cols,lrows,lcols; 1259 PetscInt local_rows,local_cols; 1260 PetscBool isseqdense,isseqsbaij,isseqaij,isseqbaij; 1261 #if defined (PETSC_USE_DEBUG) 1262 PetscBool lb[4],bb[4]; 1263 #endif 1264 PetscMPIInt nsubdomains; 1265 /* values insertion */ 1266 PetscScalar *array; 1267 /* work */ 1268 PetscErrorCode ierr; 1269 1270 PetscFunctionBegin; 1271 /* get info from mat */ 1272 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&nsubdomains);CHKERRQ(ierr); 1273 if (nsubdomains == 1) { 1274 Mat B; 1275 IS rows,cols; 1276 IS irows,icols; 1277 const PetscInt *ridxs,*cidxs; 1278 1279 ierr = ISLocalToGlobalMappingGetIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 1280 ierr = ISLocalToGlobalMappingGetIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 1281 ierr = ISCreateGeneral(PETSC_COMM_SELF,mat->rmap->n,ridxs,PETSC_USE_POINTER,&rows);CHKERRQ(ierr); 1282 ierr = ISCreateGeneral(PETSC_COMM_SELF,mat->cmap->n,cidxs,PETSC_USE_POINTER,&cols);CHKERRQ(ierr); 1283 ierr = ISLocalToGlobalMappingRestoreIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 1284 ierr = ISLocalToGlobalMappingRestoreIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 1285 ierr = ISSetPermutation(rows);CHKERRQ(ierr); 1286 ierr = ISSetPermutation(cols);CHKERRQ(ierr); 1287 ierr = ISInvertPermutation(rows,mat->rmap->n,&irows);CHKERRQ(ierr); 1288 ierr = ISInvertPermutation(cols,mat->cmap->n,&icols);CHKERRQ(ierr); 1289 ierr = ISDestroy(&cols);CHKERRQ(ierr); 1290 ierr = ISDestroy(&rows);CHKERRQ(ierr); 1291 ierr = MatConvert(matis->A,MATAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 1292 ierr = MatGetSubMatrix(B,irows,icols,reuse,M);CHKERRQ(ierr); 1293 ierr = MatDestroy(&B);CHKERRQ(ierr); 1294 ierr = ISDestroy(&icols);CHKERRQ(ierr); 1295 ierr = ISDestroy(&irows);CHKERRQ(ierr); 1296 PetscFunctionReturn(0); 1297 } 1298 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 1299 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 1300 ierr = MatGetLocalSize(mat,&lrows,&lcols);CHKERRQ(ierr); 1301 ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr); 1302 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 1303 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 1304 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);CHKERRQ(ierr); 1305 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);CHKERRQ(ierr); 1306 if (!isseqdense && !isseqaij && !isseqbaij && !isseqsbaij) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name); 1307 #if defined (PETSC_USE_DEBUG) 1308 lb[0] = isseqdense; 1309 lb[1] = isseqaij; 1310 lb[2] = isseqbaij; 1311 lb[3] = isseqsbaij; 1312 ierr = MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1313 if (!bb[0] && !bb[1] && !bb[2] && !bb[3]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type"); 1314 #endif 1315 1316 if (reuse == MAT_INITIAL_MATRIX) { 1317 ierr = MatCreate(PetscObjectComm((PetscObject)mat),M);CHKERRQ(ierr); 1318 ierr = MatSetSizes(*M,lrows,lcols,rows,cols);CHKERRQ(ierr); 1319 ierr = MatSetBlockSize(*M,bs);CHKERRQ(ierr); 1320 if (!isseqsbaij) { 1321 ierr = MatSetType(*M,MATAIJ);CHKERRQ(ierr); 1322 } else { 1323 ierr = MatSetType(*M,MATSBAIJ);CHKERRQ(ierr); 1324 } 1325 ierr = MatISSetMPIXAIJPreallocation_Private(mat,*M,PETSC_FALSE);CHKERRQ(ierr); 1326 } else { 1327 PetscInt mbs,mrows,mcols,mlrows,mlcols; 1328 /* some checks */ 1329 ierr = MatGetBlockSize(*M,&mbs);CHKERRQ(ierr); 1330 ierr = MatGetSize(*M,&mrows,&mcols);CHKERRQ(ierr); 1331 ierr = MatGetLocalSize(*M,&mlrows,&mlcols);CHKERRQ(ierr); 1332 if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows); 1333 if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols); 1334 if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows); 1335 if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols); 1336 if (mbs != bs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs); 1337 ierr = MatZeroEntries(*M);CHKERRQ(ierr); 1338 } 1339 1340 if (isseqsbaij) { 1341 ierr = MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);CHKERRQ(ierr); 1342 } else { 1343 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1344 local_mat = matis->A; 1345 } 1346 1347 /* Set values */ 1348 ierr = MatSetLocalToGlobalMapping(*M,mat->rmap->mapping,mat->cmap->mapping);CHKERRQ(ierr); 1349 if (isseqdense) { /* special case for dense local matrices */ 1350 PetscInt i,*dummy_rows,*dummy_cols; 1351 1352 if (local_rows != local_cols) { 1353 ierr = PetscMalloc2(local_rows,&dummy_rows,local_cols,&dummy_cols);CHKERRQ(ierr); 1354 for (i=0;i<local_rows;i++) dummy_rows[i] = i; 1355 for (i=0;i<local_cols;i++) dummy_cols[i] = i; 1356 } else { 1357 ierr = PetscMalloc1(local_rows,&dummy_rows);CHKERRQ(ierr); 1358 for (i=0;i<local_rows;i++) dummy_rows[i] = i; 1359 dummy_cols = dummy_rows; 1360 } 1361 ierr = MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 1362 ierr = MatDenseGetArray(local_mat,&array);CHKERRQ(ierr); 1363 ierr = MatSetValuesLocal(*M,local_rows,dummy_rows,local_cols,dummy_cols,array,ADD_VALUES);CHKERRQ(ierr); 1364 ierr = MatDenseRestoreArray(local_mat,&array);CHKERRQ(ierr); 1365 if (dummy_rows != dummy_cols) { 1366 ierr = PetscFree2(dummy_rows,dummy_cols);CHKERRQ(ierr); 1367 } else { 1368 ierr = PetscFree(dummy_rows);CHKERRQ(ierr); 1369 } 1370 } else if (isseqaij) { 1371 PetscInt i,nvtxs,*xadj,*adjncy; 1372 PetscBool done; 1373 1374 ierr = MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1375 if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ called from %s\n",__FUNCT__); 1376 ierr = MatSeqAIJGetArray(local_mat,&array);CHKERRQ(ierr); 1377 for (i=0;i<nvtxs;i++) { 1378 ierr = MatSetValuesLocal(*M,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);CHKERRQ(ierr); 1379 } 1380 ierr = MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1381 if (!done) SETERRQ1(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called from %s\n",__FUNCT__); 1382 ierr = MatSeqAIJRestoreArray(local_mat,&array);CHKERRQ(ierr); 1383 } else { /* very basic values insertion for all other matrix types */ 1384 PetscInt i; 1385 1386 for (i=0;i<local_rows;i++) { 1387 PetscInt j; 1388 const PetscInt *local_indices_cols; 1389 1390 ierr = MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1391 ierr = MatSetValuesLocal(*M,1,&i,j,local_indices_cols,array,ADD_VALUES);CHKERRQ(ierr); 1392 ierr = MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1393 } 1394 } 1395 ierr = MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1396 ierr = MatDestroy(&local_mat);CHKERRQ(ierr); 1397 ierr = MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1398 if (isseqdense) { 1399 ierr = MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 1400 } 1401 PetscFunctionReturn(0); 1402 } 1403 1404 #undef __FUNCT__ 1405 #define __FUNCT__ "MatISGetMPIXAIJ" 1406 /*@ 1407 MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format 1408 1409 Input Parameter: 1410 . mat - the matrix (should be of type MATIS) 1411 . reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1412 1413 Output Parameter: 1414 . newmat - the matrix in AIJ format 1415 1416 Level: developer 1417 1418 Notes: mat and *newmat cannot be the same object when MAT_REUSE_MATRIX is requested. 1419 1420 .seealso: MATIS 1421 @*/ 1422 PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat) 1423 { 1424 PetscErrorCode ierr; 1425 1426 PetscFunctionBegin; 1427 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1428 PetscValidLogicalCollectiveEnum(mat,reuse,2); 1429 PetscValidPointer(newmat,3); 1430 if (reuse != MAT_INITIAL_MATRIX) { 1431 PetscValidHeaderSpecific(*newmat,MAT_CLASSID,3); 1432 PetscCheckSameComm(mat,1,*newmat,3); 1433 if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix"); 1434 } 1435 ierr = PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatReuse,Mat*),(mat,reuse,newmat));CHKERRQ(ierr); 1436 PetscFunctionReturn(0); 1437 } 1438 1439 #undef __FUNCT__ 1440 #define __FUNCT__ "MatDuplicate_IS" 1441 PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat) 1442 { 1443 PetscErrorCode ierr; 1444 Mat_IS *matis = (Mat_IS*)(mat->data); 1445 PetscInt bs,m,n,M,N; 1446 Mat B,localmat; 1447 1448 PetscFunctionBegin; 1449 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 1450 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1451 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 1452 ierr = MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,mat->rmap->mapping,mat->cmap->mapping,&B);CHKERRQ(ierr); 1453 ierr = MatDuplicate(matis->A,op,&localmat);CHKERRQ(ierr); 1454 ierr = MatISSetLocalMat(B,localmat);CHKERRQ(ierr); 1455 ierr = MatDestroy(&localmat);CHKERRQ(ierr); 1456 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1457 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1458 *newmat = B; 1459 PetscFunctionReturn(0); 1460 } 1461 1462 #undef __FUNCT__ 1463 #define __FUNCT__ "MatIsHermitian_IS" 1464 static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg) 1465 { 1466 PetscErrorCode ierr; 1467 Mat_IS *matis = (Mat_IS*)A->data; 1468 PetscBool local_sym; 1469 1470 PetscFunctionBegin; 1471 ierr = MatIsHermitian(matis->A,tol,&local_sym);CHKERRQ(ierr); 1472 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1473 PetscFunctionReturn(0); 1474 } 1475 1476 #undef __FUNCT__ 1477 #define __FUNCT__ "MatIsSymmetric_IS" 1478 static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg) 1479 { 1480 PetscErrorCode ierr; 1481 Mat_IS *matis = (Mat_IS*)A->data; 1482 PetscBool local_sym; 1483 1484 PetscFunctionBegin; 1485 ierr = MatIsSymmetric(matis->A,tol,&local_sym);CHKERRQ(ierr); 1486 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1487 PetscFunctionReturn(0); 1488 } 1489 1490 #undef __FUNCT__ 1491 #define __FUNCT__ "MatIsStructurallySymmetric_IS" 1492 static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool *flg) 1493 { 1494 PetscErrorCode ierr; 1495 Mat_IS *matis = (Mat_IS*)A->data; 1496 PetscBool local_sym; 1497 1498 PetscFunctionBegin; 1499 if (A->rmap->mapping != A->cmap->mapping) { 1500 *flg = PETSC_FALSE; 1501 PetscFunctionReturn(0); 1502 } 1503 ierr = MatIsStructurallySymmetric(matis->A,&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__ "MatDestroy_IS" 1510 static PetscErrorCode MatDestroy_IS(Mat A) 1511 { 1512 PetscErrorCode ierr; 1513 Mat_IS *b = (Mat_IS*)A->data; 1514 1515 PetscFunctionBegin; 1516 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 1517 ierr = VecScatterDestroy(&b->cctx);CHKERRQ(ierr); 1518 ierr = VecScatterDestroy(&b->rctx);CHKERRQ(ierr); 1519 ierr = VecDestroy(&b->x);CHKERRQ(ierr); 1520 ierr = VecDestroy(&b->y);CHKERRQ(ierr); 1521 ierr = VecDestroy(&b->counter);CHKERRQ(ierr); 1522 ierr = ISDestroy(&b->getsub_ris);CHKERRQ(ierr); 1523 ierr = ISDestroy(&b->getsub_cis);CHKERRQ(ierr); 1524 if (b->sf != b->csf) { 1525 ierr = PetscSFDestroy(&b->csf);CHKERRQ(ierr); 1526 ierr = PetscFree2(b->csf_rootdata,b->csf_leafdata);CHKERRQ(ierr); 1527 } 1528 ierr = PetscSFDestroy(&b->sf);CHKERRQ(ierr); 1529 ierr = PetscFree2(b->sf_rootdata,b->sf_leafdata);CHKERRQ(ierr); 1530 ierr = PetscFree(A->data);CHKERRQ(ierr); 1531 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 1532 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);CHKERRQ(ierr); 1533 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);CHKERRQ(ierr); 1534 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);CHKERRQ(ierr); 1535 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);CHKERRQ(ierr); 1536 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",NULL);CHKERRQ(ierr); 1537 PetscFunctionReturn(0); 1538 } 1539 1540 #undef __FUNCT__ 1541 #define __FUNCT__ "MatMult_IS" 1542 static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y) 1543 { 1544 PetscErrorCode ierr; 1545 Mat_IS *is = (Mat_IS*)A->data; 1546 PetscScalar zero = 0.0; 1547 1548 PetscFunctionBegin; 1549 /* scatter the global vector x into the local work vector */ 1550 ierr = VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1551 ierr = VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1552 1553 /* multiply the local matrix */ 1554 ierr = MatMult(is->A,is->x,is->y);CHKERRQ(ierr); 1555 1556 /* scatter product back into global memory */ 1557 ierr = VecSet(y,zero);CHKERRQ(ierr); 1558 ierr = VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1559 ierr = VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1560 PetscFunctionReturn(0); 1561 } 1562 1563 #undef __FUNCT__ 1564 #define __FUNCT__ "MatMultAdd_IS" 1565 static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 1566 { 1567 Vec temp_vec; 1568 PetscErrorCode ierr; 1569 1570 PetscFunctionBegin; /* v3 = v2 + A * v1.*/ 1571 if (v3 != v2) { 1572 ierr = MatMult(A,v1,v3);CHKERRQ(ierr); 1573 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 1574 } else { 1575 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 1576 ierr = MatMult(A,v1,temp_vec);CHKERRQ(ierr); 1577 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 1578 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 1579 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1580 } 1581 PetscFunctionReturn(0); 1582 } 1583 1584 #undef __FUNCT__ 1585 #define __FUNCT__ "MatMultTranspose_IS" 1586 static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x) 1587 { 1588 Mat_IS *is = (Mat_IS*)A->data; 1589 PetscErrorCode ierr; 1590 1591 PetscFunctionBegin; 1592 /* scatter the global vector x into the local work vector */ 1593 ierr = VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1594 ierr = VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1595 1596 /* multiply the local matrix */ 1597 ierr = MatMultTranspose(is->A,is->y,is->x);CHKERRQ(ierr); 1598 1599 /* scatter product back into global vector */ 1600 ierr = VecSet(x,0);CHKERRQ(ierr); 1601 ierr = VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1602 ierr = VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1603 PetscFunctionReturn(0); 1604 } 1605 1606 #undef __FUNCT__ 1607 #define __FUNCT__ "MatMultTransposeAdd_IS" 1608 static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 1609 { 1610 Vec temp_vec; 1611 PetscErrorCode ierr; 1612 1613 PetscFunctionBegin; /* v3 = v2 + A' * v1.*/ 1614 if (v3 != v2) { 1615 ierr = MatMultTranspose(A,v1,v3);CHKERRQ(ierr); 1616 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 1617 } else { 1618 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 1619 ierr = MatMultTranspose(A,v1,temp_vec);CHKERRQ(ierr); 1620 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 1621 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 1622 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1623 } 1624 PetscFunctionReturn(0); 1625 } 1626 1627 #undef __FUNCT__ 1628 #define __FUNCT__ "MatView_IS" 1629 static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer) 1630 { 1631 Mat_IS *a = (Mat_IS*)A->data; 1632 PetscErrorCode ierr; 1633 PetscViewer sviewer; 1634 PetscBool isascii,view = PETSC_TRUE; 1635 1636 PetscFunctionBegin; 1637 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 1638 if (isascii) { 1639 PetscViewerFormat format; 1640 1641 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1642 if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE; 1643 } 1644 if (!view) PetscFunctionReturn(0); 1645 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 1646 ierr = MatView(a->A,sviewer);CHKERRQ(ierr); 1647 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 1648 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1649 PetscFunctionReturn(0); 1650 } 1651 1652 #undef __FUNCT__ 1653 #define __FUNCT__ "MatSetLocalToGlobalMapping_IS" 1654 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping) 1655 { 1656 PetscErrorCode ierr; 1657 PetscInt nr,rbs,nc,cbs; 1658 Mat_IS *is = (Mat_IS*)A->data; 1659 IS from,to; 1660 Vec cglobal,rglobal; 1661 1662 PetscFunctionBegin; 1663 PetscCheckSameComm(A,1,rmapping,2); 1664 PetscCheckSameComm(A,1,cmapping,3); 1665 /* Destroy any previous data */ 1666 ierr = VecDestroy(&is->x);CHKERRQ(ierr); 1667 ierr = VecDestroy(&is->y);CHKERRQ(ierr); 1668 ierr = VecDestroy(&is->counter);CHKERRQ(ierr); 1669 ierr = VecScatterDestroy(&is->rctx);CHKERRQ(ierr); 1670 ierr = VecScatterDestroy(&is->cctx);CHKERRQ(ierr); 1671 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 1672 ierr = PetscSFDestroy(&is->sf);CHKERRQ(ierr); 1673 ierr = PetscFree2(is->sf_rootdata,is->sf_leafdata);CHKERRQ(ierr); 1674 1675 /* Setup Layout and set local to global maps */ 1676 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1677 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1678 ierr = PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);CHKERRQ(ierr); 1679 ierr = PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);CHKERRQ(ierr); 1680 1681 /* Create the local matrix A */ 1682 ierr = ISLocalToGlobalMappingGetSize(rmapping,&nr);CHKERRQ(ierr); 1683 ierr = ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);CHKERRQ(ierr); 1684 ierr = ISLocalToGlobalMappingGetSize(cmapping,&nc);CHKERRQ(ierr); 1685 ierr = ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);CHKERRQ(ierr); 1686 ierr = MatCreate(PETSC_COMM_SELF,&is->A);CHKERRQ(ierr); 1687 ierr = MatSetType(is->A,MATAIJ);CHKERRQ(ierr); 1688 ierr = MatSetSizes(is->A,nr,nc,nr,nc);CHKERRQ(ierr); 1689 ierr = MatSetBlockSizes(is->A,rbs,cbs);CHKERRQ(ierr); 1690 ierr = MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);CHKERRQ(ierr); 1691 ierr = MatAppendOptionsPrefix(is->A,"is_");CHKERRQ(ierr); 1692 ierr = MatSetFromOptions(is->A);CHKERRQ(ierr); 1693 1694 if (!is->islocalref) { /* setup scatters and local vectors for MatMult */ 1695 /* Create the local work vectors */ 1696 ierr = MatCreateVecs(is->A,&is->x,&is->y);CHKERRQ(ierr); 1697 1698 /* setup the global to local scatters */ 1699 ierr = MatCreateVecs(A,&cglobal,&rglobal);CHKERRQ(ierr); 1700 ierr = ISCreateStride(PETSC_COMM_SELF,nr,0,1,&to);CHKERRQ(ierr); 1701 ierr = ISLocalToGlobalMappingApplyIS(rmapping,to,&from);CHKERRQ(ierr); 1702 ierr = VecScatterCreate(rglobal,from,is->y,to,&is->rctx);CHKERRQ(ierr); 1703 if (rmapping != cmapping) { 1704 ierr = ISDestroy(&to);CHKERRQ(ierr); 1705 ierr = ISDestroy(&from);CHKERRQ(ierr); 1706 ierr = ISCreateStride(PETSC_COMM_SELF,nc,0,1,&to);CHKERRQ(ierr); 1707 ierr = ISLocalToGlobalMappingApplyIS(cmapping,to,&from);CHKERRQ(ierr); 1708 ierr = VecScatterCreate(cglobal,from,is->x,to,&is->cctx);CHKERRQ(ierr); 1709 } else { 1710 ierr = PetscObjectReference((PetscObject)is->rctx);CHKERRQ(ierr); 1711 is->cctx = is->rctx; 1712 } 1713 1714 /* interface counter vector (local) */ 1715 ierr = VecDuplicate(is->y,&is->counter);CHKERRQ(ierr); 1716 ierr = VecSet(is->y,1.);CHKERRQ(ierr); 1717 ierr = VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1718 ierr = VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1719 ierr = VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1720 ierr = VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1721 1722 /* free workspace */ 1723 ierr = VecDestroy(&rglobal);CHKERRQ(ierr); 1724 ierr = VecDestroy(&cglobal);CHKERRQ(ierr); 1725 ierr = ISDestroy(&to);CHKERRQ(ierr); 1726 ierr = ISDestroy(&from);CHKERRQ(ierr); 1727 } 1728 ierr = MatSetUp(A);CHKERRQ(ierr); 1729 PetscFunctionReturn(0); 1730 } 1731 1732 #undef __FUNCT__ 1733 #define __FUNCT__ "MatSetValues_IS" 1734 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 1735 { 1736 Mat_IS *is = (Mat_IS*)mat->data; 1737 PetscErrorCode ierr; 1738 #if defined(PETSC_USE_DEBUG) 1739 PetscInt i,zm,zn; 1740 #endif 1741 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1742 1743 PetscFunctionBegin; 1744 #if defined(PETSC_USE_DEBUG) 1745 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); 1746 /* count negative indices */ 1747 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 1748 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 1749 #endif 1750 ierr = ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 1751 ierr = ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 1752 #if defined(PETSC_USE_DEBUG) 1753 /* count negative indices (should be the same as before) */ 1754 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 1755 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 1756 if (zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS"); 1757 if (zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS"); 1758 #endif 1759 ierr = MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 1760 PetscFunctionReturn(0); 1761 } 1762 1763 #undef __FUNCT__ 1764 #define __FUNCT__ "MatSetValuesBlocked_IS" 1765 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 1766 { 1767 Mat_IS *is = (Mat_IS*)mat->data; 1768 PetscErrorCode ierr; 1769 #if defined(PETSC_USE_DEBUG) 1770 PetscInt i,zm,zn; 1771 #endif 1772 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1773 1774 PetscFunctionBegin; 1775 #if defined(PETSC_USE_DEBUG) 1776 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); 1777 /* count negative indices */ 1778 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 1779 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 1780 #endif 1781 ierr = ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 1782 ierr = ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 1783 #if defined(PETSC_USE_DEBUG) 1784 /* count negative indices (should be the same as before) */ 1785 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 1786 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 1787 if (zm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the row indices can not be mapped! Maybe you should not use MATIS"); 1788 if (zn) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Some of the column indices can not be mapped! Maybe you should not use MATIS"); 1789 #endif 1790 ierr = MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 1791 PetscFunctionReturn(0); 1792 } 1793 1794 #undef __FUNCT__ 1795 #define __FUNCT__ "MatSetValuesLocal_IS" 1796 static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 1797 { 1798 PetscErrorCode ierr; 1799 Mat_IS *is = (Mat_IS*)A->data; 1800 1801 PetscFunctionBegin; 1802 ierr = MatSetValues(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 1803 PetscFunctionReturn(0); 1804 } 1805 1806 #undef __FUNCT__ 1807 #define __FUNCT__ "MatSetValuesBlockedLocal_IS" 1808 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 1809 { 1810 PetscErrorCode ierr; 1811 Mat_IS *is = (Mat_IS*)A->data; 1812 1813 PetscFunctionBegin; 1814 ierr = MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 1815 PetscFunctionReturn(0); 1816 } 1817 1818 #undef __FUNCT__ 1819 #define __FUNCT__ "MatISZeroRowsColumnsLocal_Private" 1820 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns) 1821 { 1822 Mat_IS *is = (Mat_IS*)A->data; 1823 PetscErrorCode ierr; 1824 1825 PetscFunctionBegin; 1826 if (!n) { 1827 is->pure_neumann = PETSC_TRUE; 1828 } else { 1829 PetscInt i; 1830 is->pure_neumann = PETSC_FALSE; 1831 1832 if (columns) { 1833 ierr = MatZeroRowsColumns(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 1834 } else { 1835 ierr = MatZeroRows(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 1836 } 1837 if (diag != 0.) { 1838 const PetscScalar *array; 1839 ierr = VecGetArrayRead(is->counter,&array);CHKERRQ(ierr); 1840 for (i=0; i<n; i++) { 1841 ierr = MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);CHKERRQ(ierr); 1842 } 1843 ierr = VecRestoreArrayRead(is->counter,&array);CHKERRQ(ierr); 1844 } 1845 ierr = MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1846 ierr = MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1847 } 1848 PetscFunctionReturn(0); 1849 } 1850 1851 #undef __FUNCT__ 1852 #define __FUNCT__ "MatZeroRowsColumns_Private_IS" 1853 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns) 1854 { 1855 Mat_IS *matis = (Mat_IS*)A->data; 1856 PetscInt nr,nl,len,i; 1857 PetscInt *lrows; 1858 PetscErrorCode ierr; 1859 1860 PetscFunctionBegin; 1861 #if defined(PETSC_USE_DEBUG) 1862 if (columns || diag != 0. || (x && b)) { 1863 PetscBool cong; 1864 ierr = PetscLayoutCompare(A->rmap,A->cmap,&cong);CHKERRQ(ierr); 1865 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"); 1866 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"); 1867 if (!cong && x && b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"A->rmap and A->cmap need to be congruent"); 1868 } 1869 #endif 1870 /* get locally owned rows */ 1871 ierr = PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);CHKERRQ(ierr); 1872 /* fix right hand side if needed */ 1873 if (x && b) { 1874 const PetscScalar *xx; 1875 PetscScalar *bb; 1876 1877 ierr = VecGetArrayRead(x, &xx);CHKERRQ(ierr); 1878 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 1879 for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]]; 1880 ierr = VecRestoreArrayRead(x, &xx);CHKERRQ(ierr); 1881 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 1882 } 1883 /* get rows associated to the local matrices */ 1884 ierr = MatISSetUpSF(A);CHKERRQ(ierr); 1885 ierr = MatGetSize(matis->A,&nl,NULL);CHKERRQ(ierr); 1886 ierr = PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 1887 ierr = PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1888 for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1; 1889 ierr = PetscFree(lrows);CHKERRQ(ierr); 1890 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1891 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1892 ierr = PetscMalloc1(nl,&lrows);CHKERRQ(ierr); 1893 for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i; 1894 ierr = MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);CHKERRQ(ierr); 1895 ierr = PetscFree(lrows);CHKERRQ(ierr); 1896 PetscFunctionReturn(0); 1897 } 1898 1899 #undef __FUNCT__ 1900 #define __FUNCT__ "MatZeroRows_IS" 1901 static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1902 { 1903 PetscErrorCode ierr; 1904 1905 PetscFunctionBegin; 1906 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);CHKERRQ(ierr); 1907 PetscFunctionReturn(0); 1908 } 1909 1910 #undef __FUNCT__ 1911 #define __FUNCT__ "MatZeroRowsColumns_IS" 1912 static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1913 { 1914 PetscErrorCode ierr; 1915 1916 PetscFunctionBegin; 1917 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);CHKERRQ(ierr); 1918 PetscFunctionReturn(0); 1919 } 1920 1921 #undef __FUNCT__ 1922 #define __FUNCT__ "MatAssemblyBegin_IS" 1923 static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type) 1924 { 1925 Mat_IS *is = (Mat_IS*)A->data; 1926 PetscErrorCode ierr; 1927 1928 PetscFunctionBegin; 1929 ierr = MatAssemblyBegin(is->A,type);CHKERRQ(ierr); 1930 PetscFunctionReturn(0); 1931 } 1932 1933 #undef __FUNCT__ 1934 #define __FUNCT__ "MatAssemblyEnd_IS" 1935 static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type) 1936 { 1937 Mat_IS *is = (Mat_IS*)A->data; 1938 PetscErrorCode ierr; 1939 1940 PetscFunctionBegin; 1941 ierr = MatAssemblyEnd(is->A,type);CHKERRQ(ierr); 1942 PetscFunctionReturn(0); 1943 } 1944 1945 #undef __FUNCT__ 1946 #define __FUNCT__ "MatISGetLocalMat_IS" 1947 static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local) 1948 { 1949 Mat_IS *is = (Mat_IS*)mat->data; 1950 1951 PetscFunctionBegin; 1952 *local = is->A; 1953 PetscFunctionReturn(0); 1954 } 1955 1956 #undef __FUNCT__ 1957 #define __FUNCT__ "MatISGetLocalMat" 1958 /*@ 1959 MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix. 1960 1961 Input Parameter: 1962 . mat - the matrix 1963 1964 Output Parameter: 1965 . local - the local matrix 1966 1967 Level: advanced 1968 1969 Notes: 1970 This can be called if you have precomputed the nonzero structure of the 1971 matrix and want to provide it to the inner matrix object to improve the performance 1972 of the MatSetValues() operation. 1973 1974 .seealso: MATIS 1975 @*/ 1976 PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local) 1977 { 1978 PetscErrorCode ierr; 1979 1980 PetscFunctionBegin; 1981 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1982 PetscValidPointer(local,2); 1983 ierr = PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 1984 PetscFunctionReturn(0); 1985 } 1986 1987 #undef __FUNCT__ 1988 #define __FUNCT__ "MatISSetLocalMat_IS" 1989 static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local) 1990 { 1991 Mat_IS *is = (Mat_IS*)mat->data; 1992 PetscInt nrows,ncols,orows,ocols; 1993 PetscErrorCode ierr; 1994 1995 PetscFunctionBegin; 1996 if (is->A) { 1997 ierr = MatGetSize(is->A,&orows,&ocols);CHKERRQ(ierr); 1998 ierr = MatGetSize(local,&nrows,&ncols);CHKERRQ(ierr); 1999 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); 2000 } 2001 ierr = PetscObjectReference((PetscObject)local);CHKERRQ(ierr); 2002 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 2003 is->A = local; 2004 PetscFunctionReturn(0); 2005 } 2006 2007 #undef __FUNCT__ 2008 #define __FUNCT__ "MatISSetLocalMat" 2009 /*@ 2010 MatISSetLocalMat - Replace the local matrix stored inside a MATIS object. 2011 2012 Input Parameter: 2013 . mat - the matrix 2014 . local - the local matrix 2015 2016 Output Parameter: 2017 2018 Level: advanced 2019 2020 Notes: 2021 This can be called if you have precomputed the local matrix and 2022 want to provide it to the matrix object MATIS. 2023 2024 .seealso: MATIS 2025 @*/ 2026 PetscErrorCode MatISSetLocalMat(Mat mat,Mat local) 2027 { 2028 PetscErrorCode ierr; 2029 2030 PetscFunctionBegin; 2031 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2032 PetscValidHeaderSpecific(local,MAT_CLASSID,2); 2033 ierr = PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));CHKERRQ(ierr); 2034 PetscFunctionReturn(0); 2035 } 2036 2037 #undef __FUNCT__ 2038 #define __FUNCT__ "MatZeroEntries_IS" 2039 static PetscErrorCode MatZeroEntries_IS(Mat A) 2040 { 2041 Mat_IS *a = (Mat_IS*)A->data; 2042 PetscErrorCode ierr; 2043 2044 PetscFunctionBegin; 2045 ierr = MatZeroEntries(a->A);CHKERRQ(ierr); 2046 PetscFunctionReturn(0); 2047 } 2048 2049 #undef __FUNCT__ 2050 #define __FUNCT__ "MatScale_IS" 2051 static PetscErrorCode MatScale_IS(Mat A,PetscScalar a) 2052 { 2053 Mat_IS *is = (Mat_IS*)A->data; 2054 PetscErrorCode ierr; 2055 2056 PetscFunctionBegin; 2057 ierr = MatScale(is->A,a);CHKERRQ(ierr); 2058 PetscFunctionReturn(0); 2059 } 2060 2061 #undef __FUNCT__ 2062 #define __FUNCT__ "MatGetDiagonal_IS" 2063 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v) 2064 { 2065 Mat_IS *is = (Mat_IS*)A->data; 2066 PetscErrorCode ierr; 2067 2068 PetscFunctionBegin; 2069 /* get diagonal of the local matrix */ 2070 ierr = MatGetDiagonal(is->A,is->y);CHKERRQ(ierr); 2071 2072 /* scatter diagonal back into global vector */ 2073 ierr = VecSet(v,0);CHKERRQ(ierr); 2074 ierr = VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2075 ierr = VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2076 PetscFunctionReturn(0); 2077 } 2078 2079 #undef __FUNCT__ 2080 #define __FUNCT__ "MatSetOption_IS" 2081 static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg) 2082 { 2083 Mat_IS *a = (Mat_IS*)A->data; 2084 PetscErrorCode ierr; 2085 2086 PetscFunctionBegin; 2087 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 2088 PetscFunctionReturn(0); 2089 } 2090 2091 #undef __FUNCT__ 2092 #define __FUNCT__ "MatAXPY_IS" 2093 static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str) 2094 { 2095 Mat_IS *y = (Mat_IS*)Y->data; 2096 Mat_IS *x; 2097 #if defined(PETSC_USE_DEBUG) 2098 PetscBool ismatis; 2099 #endif 2100 PetscErrorCode ierr; 2101 2102 PetscFunctionBegin; 2103 #if defined(PETSC_USE_DEBUG) 2104 ierr = PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);CHKERRQ(ierr); 2105 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS"); 2106 #endif 2107 x = (Mat_IS*)X->data; 2108 ierr = MatAXPY(y->A,a,x->A,str);CHKERRQ(ierr); 2109 PetscFunctionReturn(0); 2110 } 2111 2112 #undef __FUNCT__ 2113 #define __FUNCT__ "MatGetLocalSubMatrix_IS" 2114 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat) 2115 { 2116 Mat lA; 2117 Mat_IS *matis; 2118 ISLocalToGlobalMapping rl2g,cl2g; 2119 IS is; 2120 const PetscInt *rg,*rl; 2121 PetscInt nrg; 2122 PetscInt N,M,nrl,i,*idxs; 2123 PetscErrorCode ierr; 2124 2125 PetscFunctionBegin; 2126 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 2127 ierr = ISGetLocalSize(row,&nrl);CHKERRQ(ierr); 2128 ierr = ISGetIndices(row,&rl);CHKERRQ(ierr); 2129 ierr = ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);CHKERRQ(ierr); 2130 #if defined(PETSC_USE_DEBUG) 2131 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); 2132 #endif 2133 ierr = PetscMalloc1(nrg,&idxs);CHKERRQ(ierr); 2134 /* map from [0,nrl) to row */ 2135 for (i=0;i<nrl;i++) idxs[i] = rl[i]; 2136 #if defined(PETSC_USE_DEBUG) 2137 for (i=nrl;i<nrg;i++) idxs[i] = nrg; 2138 #else 2139 for (i=nrl;i<nrg;i++) idxs[i] = -1; 2140 #endif 2141 ierr = ISRestoreIndices(row,&rl);CHKERRQ(ierr); 2142 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 2143 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 2144 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 2145 ierr = ISDestroy(&is);CHKERRQ(ierr); 2146 /* compute new l2g map for columns */ 2147 if (col != row || A->rmap->mapping != A->cmap->mapping) { 2148 const PetscInt *cg,*cl; 2149 PetscInt ncg; 2150 PetscInt ncl; 2151 2152 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 2153 ierr = ISGetLocalSize(col,&ncl);CHKERRQ(ierr); 2154 ierr = ISGetIndices(col,&cl);CHKERRQ(ierr); 2155 ierr = ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);CHKERRQ(ierr); 2156 #if defined(PETSC_USE_DEBUG) 2157 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); 2158 #endif 2159 ierr = PetscMalloc1(ncg,&idxs);CHKERRQ(ierr); 2160 /* map from [0,ncl) to col */ 2161 for (i=0;i<ncl;i++) idxs[i] = cl[i]; 2162 #if defined(PETSC_USE_DEBUG) 2163 for (i=ncl;i<ncg;i++) idxs[i] = ncg; 2164 #else 2165 for (i=ncl;i<ncg;i++) idxs[i] = -1; 2166 #endif 2167 ierr = ISRestoreIndices(col,&cl);CHKERRQ(ierr); 2168 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 2169 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 2170 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 2171 ierr = ISDestroy(&is);CHKERRQ(ierr); 2172 } else { 2173 ierr = PetscObjectReference((PetscObject)rl2g);CHKERRQ(ierr); 2174 cl2g = rl2g; 2175 } 2176 /* create the MATIS submatrix */ 2177 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 2178 ierr = MatCreate(PetscObjectComm((PetscObject)A),submat);CHKERRQ(ierr); 2179 ierr = MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2180 ierr = MatSetType(*submat,MATIS);CHKERRQ(ierr); 2181 matis = (Mat_IS*)((*submat)->data); 2182 matis->islocalref = PETSC_TRUE; 2183 ierr = MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);CHKERRQ(ierr); 2184 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr); 2185 ierr = MatISSetLocalMat(*submat,lA);CHKERRQ(ierr); 2186 ierr = MatSetUp(*submat);CHKERRQ(ierr); 2187 ierr = MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2188 ierr = MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2189 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 2190 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 2191 /* remove unsupported ops */ 2192 ierr = PetscMemzero((*submat)->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2193 (*submat)->ops->destroy = MatDestroy_IS; 2194 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS; 2195 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS; 2196 (*submat)->ops->assemblybegin = MatAssemblyBegin_IS; 2197 (*submat)->ops->assemblyend = MatAssemblyEnd_IS; 2198 PetscFunctionReturn(0); 2199 } 2200 2201 #undef __FUNCT__ 2202 #define __FUNCT__ "MatCreateIS" 2203 /*@ 2204 MatCreateIS - Creates a "process" unassembled matrix, assembled on each 2205 process but not across processes. 2206 2207 Input Parameters: 2208 + comm - MPI communicator that will share the matrix 2209 . bs - block size of the matrix 2210 . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products 2211 . rmap - local to global map for rows 2212 - cmap - local to global map for cols 2213 2214 Output Parameter: 2215 . A - the resulting matrix 2216 2217 Level: advanced 2218 2219 Notes: See MATIS for more details. 2220 m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors 2221 used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices. 2222 If either rmap or cmap are NULL, then the matrix is assumed to be square. 2223 2224 .seealso: MATIS, MatSetLocalToGlobalMapping() 2225 @*/ 2226 PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A) 2227 { 2228 PetscErrorCode ierr; 2229 2230 PetscFunctionBegin; 2231 if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings"); 2232 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2233 if (bs > 0) { 2234 ierr = MatSetBlockSize(*A,bs);CHKERRQ(ierr); 2235 } 2236 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 2237 ierr = MatSetType(*A,MATIS);CHKERRQ(ierr); 2238 if (rmap && cmap) { 2239 ierr = MatSetLocalToGlobalMapping(*A,rmap,cmap);CHKERRQ(ierr); 2240 } else if (!rmap) { 2241 ierr = MatSetLocalToGlobalMapping(*A,cmap,cmap);CHKERRQ(ierr); 2242 } else { 2243 ierr = MatSetLocalToGlobalMapping(*A,rmap,rmap);CHKERRQ(ierr); 2244 } 2245 PetscFunctionReturn(0); 2246 } 2247 2248 /*MC 2249 MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP). 2250 This stores the matrices in globally unassembled form. Each processor 2251 assembles only its local Neumann problem and the parallel matrix vector 2252 product is handled "implicitly". 2253 2254 Operations Provided: 2255 + MatMult() 2256 . MatMultAdd() 2257 . MatMultTranspose() 2258 . MatMultTransposeAdd() 2259 . MatZeroEntries() 2260 . MatSetOption() 2261 . MatZeroRows() 2262 . MatSetValues() 2263 . MatSetValuesBlocked() 2264 . MatSetValuesLocal() 2265 . MatSetValuesBlockedLocal() 2266 . MatScale() 2267 . MatGetDiagonal() 2268 . MatMissingDiagonal() 2269 . MatDuplicate() 2270 . MatCopy() 2271 . MatAXPY() 2272 . MatGetSubMatrix() 2273 . MatGetLocalSubMatrix() 2274 . MatTranspose() 2275 - MatSetLocalToGlobalMapping() 2276 2277 Options Database Keys: 2278 . -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions() 2279 2280 Notes: Options prefix for the inner matrix are given by -is_mat_xxx 2281 2282 You must call MatSetLocalToGlobalMapping() before using this matrix type. 2283 2284 You can do matrix preallocation on the local matrix after you obtain it with 2285 MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation() 2286 2287 Level: advanced 2288 2289 .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP 2290 2291 M*/ 2292 2293 #undef __FUNCT__ 2294 #define __FUNCT__ "MatCreate_IS" 2295 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A) 2296 { 2297 PetscErrorCode ierr; 2298 Mat_IS *b; 2299 2300 PetscFunctionBegin; 2301 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 2302 A->data = (void*)b; 2303 2304 /* matrix ops */ 2305 ierr = PetscMemzero(A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2306 A->ops->mult = MatMult_IS; 2307 A->ops->multadd = MatMultAdd_IS; 2308 A->ops->multtranspose = MatMultTranspose_IS; 2309 A->ops->multtransposeadd = MatMultTransposeAdd_IS; 2310 A->ops->destroy = MatDestroy_IS; 2311 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS; 2312 A->ops->setvalues = MatSetValues_IS; 2313 A->ops->setvaluesblocked = MatSetValuesBlocked_IS; 2314 A->ops->setvalueslocal = MatSetValuesLocal_IS; 2315 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS; 2316 A->ops->zerorows = MatZeroRows_IS; 2317 A->ops->zerorowscolumns = MatZeroRowsColumns_IS; 2318 A->ops->assemblybegin = MatAssemblyBegin_IS; 2319 A->ops->assemblyend = MatAssemblyEnd_IS; 2320 A->ops->view = MatView_IS; 2321 A->ops->zeroentries = MatZeroEntries_IS; 2322 A->ops->scale = MatScale_IS; 2323 A->ops->getdiagonal = MatGetDiagonal_IS; 2324 A->ops->setoption = MatSetOption_IS; 2325 A->ops->ishermitian = MatIsHermitian_IS; 2326 A->ops->issymmetric = MatIsSymmetric_IS; 2327 A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS; 2328 A->ops->duplicate = MatDuplicate_IS; 2329 A->ops->missingdiagonal = MatMissingDiagonal_IS; 2330 A->ops->copy = MatCopy_IS; 2331 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 2332 A->ops->getsubmatrix = MatGetSubMatrix_IS; 2333 A->ops->axpy = MatAXPY_IS; 2334 A->ops->diagonalset = MatDiagonalSet_IS; 2335 A->ops->shift = MatShift_IS; 2336 A->ops->transpose = MatTranspose_IS; 2337 A->ops->getinfo = MatGetInfo_IS; 2338 A->ops->diagonalscale = MatDiagonalScale_IS; 2339 2340 /* special MATIS functions */ 2341 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);CHKERRQ(ierr); 2342 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);CHKERRQ(ierr); 2343 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatISGetMPIXAIJ_IS);CHKERRQ(ierr); 2344 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);CHKERRQ(ierr); 2345 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",MatISSetUpSF_IS);CHKERRQ(ierr); 2346 ierr = PetscObjectChangeTypeName((PetscObject)A,MATIS);CHKERRQ(ierr); 2347 PetscFunctionReturn(0); 2348 } 2349