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