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