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 PetscBool cong; 1046 MPI_Comm comm; 1047 1048 ierr = PetscObjectGetComm((PetscObject)mat,&comm);CHKERRQ(ierr); 1049 ierr = MatGetBlockSizes(mat,&arbs,&acbs);CHKERRQ(ierr); 1050 ierr = ISGetBlockSize(irow,&irbs);CHKERRQ(ierr); 1051 ierr = ISGetBlockSize(icol,&icbs);CHKERRQ(ierr); 1052 rbs = arbs == irbs ? irbs : 1; 1053 cbs = acbs == icbs ? icbs : 1; 1054 ierr = ISGetLocalSize(irow,&m);CHKERRQ(ierr); 1055 ierr = ISGetLocalSize(icol,&n);CHKERRQ(ierr); 1056 ierr = MatCreate(comm,newmat);CHKERRQ(ierr); 1057 ierr = MatSetType(*newmat,MATIS);CHKERRQ(ierr); 1058 ierr = MatSetSizes(*newmat,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 1059 ierr = MatSetBlockSizes(*newmat,rbs,cbs);CHKERRQ(ierr); 1060 /* communicate irow to their owners in the layout */ 1061 ierr = ISGetIndices(irow,&idxs);CHKERRQ(ierr); 1062 ierr = PetscLayoutMapLocal_Private(mat->rmap,m,idxs,&ll,&lidxs,&lgidxs);CHKERRQ(ierr); 1063 ierr = ISRestoreIndices(irow,&idxs);CHKERRQ(ierr); 1064 ierr = MatISSetUpSF(mat);CHKERRQ(ierr); 1065 ierr = PetscMemzero(matis->sf_rootdata,matis->sf->nroots*sizeof(PetscInt));CHKERRQ(ierr); 1066 for (i=0;i<ll;i++) matis->sf_rootdata[lidxs[i]] = lgidxs[i]+1; 1067 ierr = PetscFree(lidxs);CHKERRQ(ierr); 1068 ierr = PetscFree(lgidxs);CHKERRQ(ierr); 1069 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1070 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1071 for (i=0,newloc=0;i<matis->sf->nleaves;i++) if (matis->sf_leafdata[i]) newloc++; 1072 ierr = PetscMalloc1(newloc,&newgidxs);CHKERRQ(ierr); 1073 ierr = PetscMalloc1(newloc,&lidxs);CHKERRQ(ierr); 1074 for (i=0,newloc=0;i<matis->sf->nleaves;i++) 1075 if (matis->sf_leafdata[i]) { 1076 lidxs[newloc] = i; 1077 newgidxs[newloc++] = matis->sf_leafdata[i]-1; 1078 } 1079 ierr = ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 1080 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 1081 ierr = ISLocalToGlobalMappingSetBlockSize(rl2g,rbs);CHKERRQ(ierr); 1082 ierr = ISDestroy(&is);CHKERRQ(ierr); 1083 /* local is to extract local submatrix */ 1084 newmatis = (Mat_IS*)(*newmat)->data; 1085 ierr = ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_ris);CHKERRQ(ierr); 1086 ierr = MatHasCongruentLayouts(mat,&cong);CHKERRQ(ierr); 1087 if (cong && irow == icol && matis->csf == matis->sf) { 1088 ierr = MatSetLocalToGlobalMapping(*newmat,rl2g,rl2g);CHKERRQ(ierr); 1089 ierr = PetscObjectReference((PetscObject)newmatis->getsub_ris);CHKERRQ(ierr); 1090 newmatis->getsub_cis = newmatis->getsub_ris; 1091 } else { 1092 ISLocalToGlobalMapping cl2g; 1093 1094 /* communicate icol to their owners in the layout */ 1095 ierr = ISGetIndices(icol,&idxs);CHKERRQ(ierr); 1096 ierr = PetscLayoutMapLocal_Private(mat->cmap,n,idxs,&ll,&lidxs,&lgidxs);CHKERRQ(ierr); 1097 ierr = ISRestoreIndices(icol,&idxs);CHKERRQ(ierr); 1098 ierr = PetscMemzero(matis->csf_rootdata,matis->csf->nroots*sizeof(PetscInt));CHKERRQ(ierr); 1099 for (i=0;i<ll;i++) matis->csf_rootdata[lidxs[i]] = lgidxs[i]+1; 1100 ierr = PetscFree(lidxs);CHKERRQ(ierr); 1101 ierr = PetscFree(lgidxs);CHKERRQ(ierr); 1102 ierr = PetscSFBcastBegin(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);CHKERRQ(ierr); 1103 ierr = PetscSFBcastEnd(matis->csf,MPIU_INT,matis->csf_rootdata,matis->csf_leafdata);CHKERRQ(ierr); 1104 for (i=0,newloc=0;i<matis->csf->nleaves;i++) if (matis->csf_leafdata[i]) newloc++; 1105 ierr = PetscMalloc1(newloc,&newgidxs);CHKERRQ(ierr); 1106 ierr = PetscMalloc1(newloc,&lidxs);CHKERRQ(ierr); 1107 for (i=0,newloc=0;i<matis->csf->nleaves;i++) 1108 if (matis->csf_leafdata[i]) { 1109 lidxs[newloc] = i; 1110 newgidxs[newloc++] = matis->csf_leafdata[i]-1; 1111 } 1112 ierr = ISCreateGeneral(comm,newloc,newgidxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 1113 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 1114 ierr = ISLocalToGlobalMappingSetBlockSize(cl2g,cbs);CHKERRQ(ierr); 1115 ierr = ISDestroy(&is);CHKERRQ(ierr); 1116 /* local is to extract local submatrix */ 1117 ierr = ISCreateGeneral(comm,newloc,lidxs,PETSC_OWN_POINTER,&newmatis->getsub_cis);CHKERRQ(ierr); 1118 ierr = MatSetLocalToGlobalMapping(*newmat,rl2g,cl2g);CHKERRQ(ierr); 1119 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 1120 } 1121 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 1122 } else { 1123 ierr = MatISGetLocalMat(*newmat,&newlocmat);CHKERRQ(ierr); 1124 } 1125 ierr = MatISGetLocalMat(mat,&locmat);CHKERRQ(ierr); 1126 newmatis = (Mat_IS*)(*newmat)->data; 1127 ierr = MatCreateSubMatrix(locmat,newmatis->getsub_ris,newmatis->getsub_cis,scall,&newlocmat);CHKERRQ(ierr); 1128 if (scall == MAT_INITIAL_MATRIX) { 1129 ierr = MatISSetLocalMat(*newmat,newlocmat);CHKERRQ(ierr); 1130 ierr = MatDestroy(&newlocmat);CHKERRQ(ierr); 1131 } 1132 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1133 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1134 PetscFunctionReturn(0); 1135 } 1136 1137 static PetscErrorCode MatCopy_IS(Mat A,Mat B,MatStructure str) 1138 { 1139 Mat_IS *a = (Mat_IS*)A->data,*b; 1140 PetscBool ismatis; 1141 PetscErrorCode ierr; 1142 1143 PetscFunctionBegin; 1144 ierr = PetscObjectTypeCompare((PetscObject)B,MATIS,&ismatis);CHKERRQ(ierr); 1145 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Need to be implemented"); 1146 b = (Mat_IS*)B->data; 1147 ierr = MatCopy(a->A,b->A,str);CHKERRQ(ierr); 1148 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 1149 PetscFunctionReturn(0); 1150 } 1151 1152 static PetscErrorCode MatMissingDiagonal_IS(Mat A,PetscBool *missing,PetscInt *d) 1153 { 1154 Vec v; 1155 const PetscScalar *array; 1156 PetscInt i,n; 1157 PetscErrorCode ierr; 1158 1159 PetscFunctionBegin; 1160 *missing = PETSC_FALSE; 1161 ierr = MatCreateVecs(A,NULL,&v);CHKERRQ(ierr); 1162 ierr = MatGetDiagonal(A,v);CHKERRQ(ierr); 1163 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 1164 ierr = VecGetArrayRead(v,&array);CHKERRQ(ierr); 1165 for (i=0;i<n;i++) if (array[i] == 0.) break; 1166 ierr = VecRestoreArrayRead(v,&array);CHKERRQ(ierr); 1167 ierr = VecDestroy(&v);CHKERRQ(ierr); 1168 if (i != n) *missing = PETSC_TRUE; 1169 if (d) { 1170 *d = -1; 1171 if (*missing) { 1172 PetscInt rstart; 1173 ierr = MatGetOwnershipRange(A,&rstart,NULL);CHKERRQ(ierr); 1174 *d = i+rstart; 1175 } 1176 } 1177 PetscFunctionReturn(0); 1178 } 1179 1180 static PetscErrorCode MatISSetUpSF_IS(Mat B) 1181 { 1182 Mat_IS *matis = (Mat_IS*)(B->data); 1183 const PetscInt *gidxs; 1184 PetscInt nleaves; 1185 PetscErrorCode ierr; 1186 1187 PetscFunctionBegin; 1188 if (matis->sf) PetscFunctionReturn(0); 1189 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1190 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1191 ierr = PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->sf);CHKERRQ(ierr); 1192 ierr = ISLocalToGlobalMappingGetIndices(B->rmap->mapping,&gidxs);CHKERRQ(ierr); 1193 ierr = ISLocalToGlobalMappingGetSize(B->rmap->mapping,&nleaves);CHKERRQ(ierr); 1194 ierr = PetscSFSetGraphLayout(matis->sf,B->rmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);CHKERRQ(ierr); 1195 ierr = ISLocalToGlobalMappingRestoreIndices(B->rmap->mapping,&gidxs);CHKERRQ(ierr); 1196 ierr = PetscMalloc2(matis->sf->nroots,&matis->sf_rootdata,matis->sf->nleaves,&matis->sf_leafdata);CHKERRQ(ierr); 1197 if (B->rmap->mapping != B->cmap->mapping) { /* setup SF for columns */ 1198 ierr = ISLocalToGlobalMappingGetSize(B->cmap->mapping,&nleaves);CHKERRQ(ierr); 1199 ierr = PetscSFCreate(PetscObjectComm((PetscObject)B),&matis->csf);CHKERRQ(ierr); 1200 ierr = ISLocalToGlobalMappingGetIndices(B->cmap->mapping,&gidxs);CHKERRQ(ierr); 1201 ierr = PetscSFSetGraphLayout(matis->csf,B->cmap,nleaves,NULL,PETSC_OWN_POINTER,gidxs);CHKERRQ(ierr); 1202 ierr = ISLocalToGlobalMappingRestoreIndices(B->cmap->mapping,&gidxs);CHKERRQ(ierr); 1203 ierr = PetscMalloc2(matis->csf->nroots,&matis->csf_rootdata,matis->csf->nleaves,&matis->csf_leafdata);CHKERRQ(ierr); 1204 } else { 1205 matis->csf = matis->sf; 1206 matis->csf_leafdata = matis->sf_leafdata; 1207 matis->csf_rootdata = matis->sf_rootdata; 1208 } 1209 PetscFunctionReturn(0); 1210 } 1211 1212 /*@ 1213 MatISStoreL2L - Store local-to-local operators during the Galerkin process of MatPtAP. 1214 1215 Collective on MPI_Comm 1216 1217 Input Parameters: 1218 + A - the matrix 1219 - store - the boolean flag 1220 1221 Level: advanced 1222 1223 Notes: 1224 1225 .keywords: matrix 1226 1227 .seealso: MatCreate(), MatCreateIS(), MatISSetPreallocation(), MatPtAP() 1228 @*/ 1229 PetscErrorCode MatISStoreL2L(Mat A, PetscBool store) 1230 { 1231 PetscErrorCode ierr; 1232 1233 PetscFunctionBegin; 1234 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1235 PetscValidType(A,1); 1236 PetscValidLogicalCollectiveBool(A,store,2); 1237 ierr = PetscTryMethod(A,"MatISStoreL2L_C",(Mat,PetscBool),(A,store));CHKERRQ(ierr); 1238 PetscFunctionReturn(0); 1239 } 1240 1241 static PetscErrorCode MatISStoreL2L_IS(Mat A, PetscBool store) 1242 { 1243 Mat_IS *matis = (Mat_IS*)(A->data); 1244 PetscErrorCode ierr; 1245 1246 PetscFunctionBegin; 1247 matis->storel2l = store; 1248 if (!store) { 1249 ierr = PetscObjectCompose((PetscObject)(A),"_MatIS_PtAP_l2l",NULL);CHKERRQ(ierr); 1250 } 1251 PetscFunctionReturn(0); 1252 } 1253 1254 /*@ 1255 MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix. 1256 1257 Collective on MPI_Comm 1258 1259 Input Parameters: 1260 + B - the matrix 1261 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1262 (same value is used for all local rows) 1263 . d_nnz - array containing the number of nonzeros in the various rows of the 1264 DIAGONAL portion of the local submatrix (possibly different for each row) 1265 or NULL, if d_nz is used to specify the nonzero structure. 1266 The size of this array is equal to the number of local rows, i.e 'm'. 1267 For matrices that will be factored, you must leave room for (and set) 1268 the diagonal entry even if it is zero. 1269 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1270 submatrix (same value is used for all local rows). 1271 - o_nnz - array containing the number of nonzeros in the various rows of the 1272 OFF-DIAGONAL portion of the local submatrix (possibly different for 1273 each row) or NULL, if o_nz is used to specify the nonzero 1274 structure. The size of this array is equal to the number 1275 of local rows, i.e 'm'. 1276 1277 If the *_nnz parameter is given then the *_nz parameter is ignored 1278 1279 Level: intermediate 1280 1281 Notes: 1282 This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition 1283 from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local 1284 matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects. 1285 1286 .keywords: matrix 1287 1288 .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS 1289 @*/ 1290 PetscErrorCode MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1291 { 1292 PetscErrorCode ierr; 1293 1294 PetscFunctionBegin; 1295 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1296 PetscValidType(B,1); 1297 ierr = PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));CHKERRQ(ierr); 1298 PetscFunctionReturn(0); 1299 } 1300 1301 /* this is used by DMDA */ 1302 PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1303 { 1304 Mat_IS *matis = (Mat_IS*)(B->data); 1305 PetscInt bs,i,nlocalcols; 1306 PetscErrorCode ierr; 1307 1308 PetscFunctionBegin; 1309 if (!matis->A) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"You should first call MatSetLocalToGlobalMapping"); 1310 ierr = MatISSetUpSF(B);CHKERRQ(ierr); 1311 1312 if (!d_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nz; 1313 else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nnz[i]; 1314 1315 if (!o_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nz; 1316 else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nnz[i]; 1317 1318 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1319 ierr = MatGetSize(matis->A,NULL,&nlocalcols);CHKERRQ(ierr); 1320 ierr = MatGetBlockSize(matis->A,&bs);CHKERRQ(ierr); 1321 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1322 1323 for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols); 1324 ierr = MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);CHKERRQ(ierr); 1325 #if defined(PETSC_HAVE_HYPRE) 1326 ierr = MatHYPRESetPreallocation(matis->A,0,matis->sf_leafdata,0,NULL);CHKERRQ(ierr); 1327 #endif 1328 1329 for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs; 1330 ierr = MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);CHKERRQ(ierr); 1331 1332 for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i]-i; 1333 ierr = MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);CHKERRQ(ierr); 1334 1335 /* for other matrix types */ 1336 ierr = MatSetUp(matis->A);CHKERRQ(ierr); 1337 PetscFunctionReturn(0); 1338 } 1339 1340 PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce) 1341 { 1342 Mat_IS *matis = (Mat_IS*)(A->data); 1343 PetscInt *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership; 1344 const PetscInt *global_indices_r,*global_indices_c; 1345 PetscInt i,j,bs,rows,cols; 1346 PetscInt lrows,lcols; 1347 PetscInt local_rows,local_cols; 1348 PetscMPIInt nsubdomains; 1349 PetscBool isdense,issbaij; 1350 PetscErrorCode ierr; 1351 1352 PetscFunctionBegin; 1353 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&nsubdomains);CHKERRQ(ierr); 1354 ierr = MatGetSize(A,&rows,&cols);CHKERRQ(ierr); 1355 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 1356 ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr); 1357 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);CHKERRQ(ierr); 1358 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);CHKERRQ(ierr); 1359 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);CHKERRQ(ierr); 1360 if (A->rmap->mapping != A->cmap->mapping) { 1361 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&global_indices_c);CHKERRQ(ierr); 1362 } else { 1363 global_indices_c = global_indices_r; 1364 } 1365 1366 if (issbaij) { 1367 ierr = MatGetRowUpperTriangular(matis->A);CHKERRQ(ierr); 1368 } 1369 /* 1370 An SF reduce is needed to sum up properly on shared rows. 1371 Note that generally preallocation is not exact, since it overestimates nonzeros 1372 */ 1373 ierr = MatISSetUpSF(A);CHKERRQ(ierr); 1374 ierr = MatGetLocalSize(A,&lrows,&lcols);CHKERRQ(ierr); 1375 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);CHKERRQ(ierr); 1376 /* All processes need to compute entire row ownership */ 1377 ierr = PetscMalloc1(rows,&row_ownership);CHKERRQ(ierr); 1378 ierr = MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);CHKERRQ(ierr); 1379 for (i=0;i<nsubdomains;i++) { 1380 for (j=mat_ranges[i];j<mat_ranges[i+1];j++) { 1381 row_ownership[j] = i; 1382 } 1383 } 1384 ierr = MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);CHKERRQ(ierr); 1385 1386 /* 1387 my_dnz and my_onz contains exact contribution to preallocation from each local mat 1388 then, they will be summed up properly. This way, preallocation is always sufficient 1389 */ 1390 ierr = PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);CHKERRQ(ierr); 1391 /* preallocation as a MATAIJ */ 1392 if (isdense) { /* special case for dense local matrices */ 1393 for (i=0;i<local_rows;i++) { 1394 PetscInt owner = row_ownership[global_indices_r[i]]; 1395 for (j=0;j<local_cols;j++) { 1396 PetscInt index_col = global_indices_c[j]; 1397 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 1398 my_dnz[i] += 1; 1399 } else { /* offdiag block */ 1400 my_onz[i] += 1; 1401 } 1402 } 1403 } 1404 } else if (matis->A->ops->getrowij) { 1405 const PetscInt *ii,*jj,*jptr; 1406 PetscBool done; 1407 ierr = MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);CHKERRQ(ierr); 1408 if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ"); 1409 jptr = jj; 1410 for (i=0;i<local_rows;i++) { 1411 PetscInt index_row = global_indices_r[i]; 1412 for (j=0;j<ii[i+1]-ii[i];j++,jptr++) { 1413 PetscInt owner = row_ownership[index_row]; 1414 PetscInt index_col = global_indices_c[*jptr]; 1415 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 1416 my_dnz[i] += 1; 1417 } else { /* offdiag block */ 1418 my_onz[i] += 1; 1419 } 1420 /* same as before, interchanging rows and cols */ 1421 if (issbaij && index_col != index_row) { 1422 owner = row_ownership[index_col]; 1423 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) { 1424 my_dnz[*jptr] += 1; 1425 } else { 1426 my_onz[*jptr] += 1; 1427 } 1428 } 1429 } 1430 } 1431 ierr = MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);CHKERRQ(ierr); 1432 if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ"); 1433 } else { /* loop over rows and use MatGetRow */ 1434 for (i=0;i<local_rows;i++) { 1435 const PetscInt *cols; 1436 PetscInt ncols,index_row = global_indices_r[i]; 1437 ierr = MatGetRow(matis->A,i,&ncols,&cols,NULL);CHKERRQ(ierr); 1438 for (j=0;j<ncols;j++) { 1439 PetscInt owner = row_ownership[index_row]; 1440 PetscInt index_col = global_indices_c[cols[j]]; 1441 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 1442 my_dnz[i] += 1; 1443 } else { /* offdiag block */ 1444 my_onz[i] += 1; 1445 } 1446 /* same as before, interchanging rows and cols */ 1447 if (issbaij && index_col != index_row) { 1448 owner = row_ownership[index_col]; 1449 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) { 1450 my_dnz[cols[j]] += 1; 1451 } else { 1452 my_onz[cols[j]] += 1; 1453 } 1454 } 1455 } 1456 ierr = MatRestoreRow(matis->A,i,&ncols,&cols,NULL);CHKERRQ(ierr); 1457 } 1458 } 1459 if (global_indices_c != global_indices_r) { 1460 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&global_indices_c);CHKERRQ(ierr); 1461 } 1462 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);CHKERRQ(ierr); 1463 ierr = PetscFree(row_ownership);CHKERRQ(ierr); 1464 1465 /* Reduce my_dnz and my_onz */ 1466 if (maxreduce) { 1467 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);CHKERRQ(ierr); 1468 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);CHKERRQ(ierr); 1469 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);CHKERRQ(ierr); 1470 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);CHKERRQ(ierr); 1471 } else { 1472 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);CHKERRQ(ierr); 1473 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);CHKERRQ(ierr); 1474 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);CHKERRQ(ierr); 1475 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);CHKERRQ(ierr); 1476 } 1477 ierr = PetscFree2(my_dnz,my_onz);CHKERRQ(ierr); 1478 1479 /* Resize preallocation if overestimated */ 1480 for (i=0;i<lrows;i++) { 1481 dnz[i] = PetscMin(dnz[i],lcols); 1482 onz[i] = PetscMin(onz[i],cols-lcols); 1483 } 1484 1485 /* Set preallocation */ 1486 ierr = MatSeqAIJSetPreallocation(B,0,dnz);CHKERRQ(ierr); 1487 ierr = MatMPIAIJSetPreallocation(B,0,dnz,0,onz);CHKERRQ(ierr); 1488 for (i=0;i<lrows/bs;i++) { 1489 dnz[i] = dnz[i*bs]/bs; 1490 onz[i] = onz[i*bs]/bs; 1491 } 1492 ierr = MatSeqBAIJSetPreallocation(B,bs,0,dnz);CHKERRQ(ierr); 1493 ierr = MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 1494 ierr = MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 1495 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1496 if (issbaij) { 1497 ierr = MatRestoreRowUpperTriangular(matis->A);CHKERRQ(ierr); 1498 } 1499 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1500 PetscFunctionReturn(0); 1501 } 1502 1503 static PetscErrorCode MatISGetMPIXAIJ_IS(Mat mat, MatReuse reuse, Mat *M) 1504 { 1505 Mat_IS *matis = (Mat_IS*)(mat->data); 1506 Mat local_mat; 1507 /* info on mat */ 1508 PetscInt bs,rows,cols,lrows,lcols; 1509 PetscInt local_rows,local_cols; 1510 PetscBool isseqdense,isseqsbaij,isseqaij,isseqbaij; 1511 #if defined (PETSC_USE_DEBUG) 1512 PetscBool lb[4],bb[4]; 1513 #endif 1514 PetscMPIInt nsubdomains; 1515 /* values insertion */ 1516 PetscScalar *array; 1517 /* work */ 1518 PetscErrorCode ierr; 1519 1520 PetscFunctionBegin; 1521 /* get info from mat */ 1522 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&nsubdomains);CHKERRQ(ierr); 1523 if (nsubdomains == 1) { 1524 Mat B; 1525 IS rows,cols; 1526 IS irows,icols; 1527 const PetscInt *ridxs,*cidxs; 1528 1529 ierr = ISLocalToGlobalMappingGetIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 1530 ierr = ISLocalToGlobalMappingGetIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 1531 ierr = ISCreateGeneral(PETSC_COMM_SELF,mat->rmap->n,ridxs,PETSC_USE_POINTER,&rows);CHKERRQ(ierr); 1532 ierr = ISCreateGeneral(PETSC_COMM_SELF,mat->cmap->n,cidxs,PETSC_USE_POINTER,&cols);CHKERRQ(ierr); 1533 ierr = ISSetPermutation(rows);CHKERRQ(ierr); 1534 ierr = ISSetPermutation(cols);CHKERRQ(ierr); 1535 ierr = ISInvertPermutation(rows,mat->rmap->n,&irows);CHKERRQ(ierr); 1536 ierr = ISInvertPermutation(cols,mat->cmap->n,&icols);CHKERRQ(ierr); 1537 ierr = ISLocalToGlobalMappingRestoreIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 1538 ierr = ISLocalToGlobalMappingRestoreIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 1539 ierr = ISDestroy(&cols);CHKERRQ(ierr); 1540 ierr = ISDestroy(&rows);CHKERRQ(ierr); 1541 ierr = MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 1542 ierr = MatCreateSubMatrix(B,irows,icols,reuse,M);CHKERRQ(ierr); 1543 ierr = MatDestroy(&B);CHKERRQ(ierr); 1544 ierr = ISDestroy(&icols);CHKERRQ(ierr); 1545 ierr = ISDestroy(&irows);CHKERRQ(ierr); 1546 PetscFunctionReturn(0); 1547 } 1548 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 1549 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 1550 ierr = MatGetLocalSize(mat,&lrows,&lcols);CHKERRQ(ierr); 1551 ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr); 1552 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 1553 ierr = PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 1554 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);CHKERRQ(ierr); 1555 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);CHKERRQ(ierr); 1556 if (!isseqdense && !isseqaij && !isseqbaij && !isseqsbaij) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name); 1557 #if defined (PETSC_USE_DEBUG) 1558 lb[0] = isseqdense; 1559 lb[1] = isseqaij; 1560 lb[2] = isseqbaij; 1561 lb[3] = isseqsbaij; 1562 ierr = MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1563 if (!bb[0] && !bb[1] && !bb[2] && !bb[3]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type"); 1564 #endif 1565 1566 if (reuse == MAT_INITIAL_MATRIX) { 1567 ierr = MatCreate(PetscObjectComm((PetscObject)mat),M);CHKERRQ(ierr); 1568 ierr = MatSetSizes(*M,lrows,lcols,rows,cols);CHKERRQ(ierr); 1569 if (!isseqsbaij) { 1570 ierr = MatSetType(*M,MATAIJ);CHKERRQ(ierr); 1571 } else { 1572 ierr = MatSetType(*M,MATSBAIJ);CHKERRQ(ierr); 1573 } 1574 ierr = MatSetBlockSize(*M,bs);CHKERRQ(ierr); 1575 ierr = MatISSetMPIXAIJPreallocation_Private(mat,*M,PETSC_FALSE);CHKERRQ(ierr); 1576 } else { 1577 PetscInt mbs,mrows,mcols,mlrows,mlcols; 1578 /* some checks */ 1579 ierr = MatGetBlockSize(*M,&mbs);CHKERRQ(ierr); 1580 ierr = MatGetSize(*M,&mrows,&mcols);CHKERRQ(ierr); 1581 ierr = MatGetLocalSize(*M,&mlrows,&mlcols);CHKERRQ(ierr); 1582 if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows); 1583 if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols); 1584 if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows); 1585 if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols); 1586 if (mbs != bs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs); 1587 ierr = MatZeroEntries(*M);CHKERRQ(ierr); 1588 } 1589 1590 if (isseqsbaij) { 1591 ierr = MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);CHKERRQ(ierr); 1592 } else { 1593 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1594 local_mat = matis->A; 1595 } 1596 1597 /* Set values */ 1598 ierr = MatSetLocalToGlobalMapping(*M,mat->rmap->mapping,mat->cmap->mapping);CHKERRQ(ierr); 1599 if (isseqdense) { /* special case for dense local matrices */ 1600 PetscInt i,*dummy; 1601 1602 ierr = PetscMalloc1(PetscMax(local_rows,local_cols),&dummy);CHKERRQ(ierr); 1603 for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i; 1604 ierr = MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 1605 ierr = MatDenseGetArray(local_mat,&array);CHKERRQ(ierr); 1606 ierr = MatSetValuesLocal(*M,local_rows,dummy,local_cols,dummy,array,ADD_VALUES);CHKERRQ(ierr); 1607 ierr = MatDenseRestoreArray(local_mat,&array);CHKERRQ(ierr); 1608 ierr = PetscFree(dummy);CHKERRQ(ierr); 1609 } else if (isseqaij) { 1610 PetscInt i,nvtxs,*xadj,*adjncy; 1611 PetscBool done; 1612 1613 ierr = MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1614 if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ"); 1615 ierr = MatSeqAIJGetArray(local_mat,&array);CHKERRQ(ierr); 1616 for (i=0;i<nvtxs;i++) { 1617 ierr = MatSetValuesLocal(*M,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);CHKERRQ(ierr); 1618 } 1619 ierr = MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1620 if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ"); 1621 ierr = MatSeqAIJRestoreArray(local_mat,&array);CHKERRQ(ierr); 1622 } else { /* very basic values insertion for all other matrix types */ 1623 PetscInt i; 1624 1625 for (i=0;i<local_rows;i++) { 1626 PetscInt j; 1627 const PetscInt *local_indices_cols; 1628 1629 ierr = MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1630 ierr = MatSetValuesLocal(*M,1,&i,j,local_indices_cols,array,ADD_VALUES);CHKERRQ(ierr); 1631 ierr = MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1632 } 1633 } 1634 ierr = MatAssemblyBegin(*M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1635 ierr = MatDestroy(&local_mat);CHKERRQ(ierr); 1636 ierr = MatAssemblyEnd(*M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1637 if (isseqdense) { 1638 ierr = MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 1639 } 1640 PetscFunctionReturn(0); 1641 } 1642 1643 /*@ 1644 MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format 1645 1646 Input Parameter: 1647 . mat - the matrix (should be of type MATIS) 1648 . reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1649 1650 Output Parameter: 1651 . newmat - the matrix in AIJ format 1652 1653 Level: developer 1654 1655 Notes: 1656 mat and *newmat cannot be the same object when MAT_REUSE_MATRIX is requested. 1657 1658 .seealso: MATIS 1659 @*/ 1660 PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat) 1661 { 1662 PetscErrorCode ierr; 1663 1664 PetscFunctionBegin; 1665 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1666 PetscValidLogicalCollectiveEnum(mat,reuse,2); 1667 PetscValidPointer(newmat,3); 1668 if (reuse != MAT_INITIAL_MATRIX) { 1669 PetscValidHeaderSpecific(*newmat,MAT_CLASSID,3); 1670 PetscCheckSameComm(mat,1,*newmat,3); 1671 if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix"); 1672 } 1673 ierr = PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatReuse,Mat*),(mat,reuse,newmat));CHKERRQ(ierr); 1674 PetscFunctionReturn(0); 1675 } 1676 1677 PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat) 1678 { 1679 PetscErrorCode ierr; 1680 Mat_IS *matis = (Mat_IS*)(mat->data); 1681 PetscInt bs,m,n,M,N; 1682 Mat B,localmat; 1683 1684 PetscFunctionBegin; 1685 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 1686 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1687 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 1688 ierr = MatCreateIS(PetscObjectComm((PetscObject)mat),bs,m,n,M,N,mat->rmap->mapping,mat->cmap->mapping,&B);CHKERRQ(ierr); 1689 ierr = MatDuplicate(matis->A,op,&localmat);CHKERRQ(ierr); 1690 ierr = MatISSetLocalMat(B,localmat);CHKERRQ(ierr); 1691 ierr = MatDestroy(&localmat);CHKERRQ(ierr); 1692 if (matis->sf) { 1693 Mat_IS *bmatis = (Mat_IS*)(B->data); 1694 1695 ierr = PetscObjectReference((PetscObject)matis->sf);CHKERRQ(ierr); 1696 bmatis->sf = matis->sf; 1697 ierr = PetscMalloc2(matis->sf->nroots,&bmatis->sf_rootdata,matis->sf->nleaves,&bmatis->sf_leafdata);CHKERRQ(ierr); 1698 if (matis->sf != matis->csf) { 1699 ierr = PetscObjectReference((PetscObject)matis->csf);CHKERRQ(ierr); 1700 bmatis->csf = matis->csf; 1701 ierr = PetscMalloc2(matis->csf->nroots,&bmatis->csf_rootdata,matis->csf->nleaves,&bmatis->csf_leafdata);CHKERRQ(ierr); 1702 } else { 1703 bmatis->csf = bmatis->sf; 1704 bmatis->csf_leafdata = bmatis->sf_leafdata; 1705 bmatis->csf_rootdata = bmatis->sf_rootdata; 1706 } 1707 } 1708 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1709 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1710 *newmat = B; 1711 PetscFunctionReturn(0); 1712 } 1713 1714 static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg) 1715 { 1716 PetscErrorCode ierr; 1717 Mat_IS *matis = (Mat_IS*)A->data; 1718 PetscBool local_sym; 1719 1720 PetscFunctionBegin; 1721 ierr = MatIsHermitian(matis->A,tol,&local_sym);CHKERRQ(ierr); 1722 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1723 PetscFunctionReturn(0); 1724 } 1725 1726 static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg) 1727 { 1728 PetscErrorCode ierr; 1729 Mat_IS *matis = (Mat_IS*)A->data; 1730 PetscBool local_sym; 1731 1732 PetscFunctionBegin; 1733 ierr = MatIsSymmetric(matis->A,tol,&local_sym);CHKERRQ(ierr); 1734 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1735 PetscFunctionReturn(0); 1736 } 1737 1738 static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool *flg) 1739 { 1740 PetscErrorCode ierr; 1741 Mat_IS *matis = (Mat_IS*)A->data; 1742 PetscBool local_sym; 1743 1744 PetscFunctionBegin; 1745 if (A->rmap->mapping != A->cmap->mapping) { 1746 *flg = PETSC_FALSE; 1747 PetscFunctionReturn(0); 1748 } 1749 ierr = MatIsStructurallySymmetric(matis->A,&local_sym);CHKERRQ(ierr); 1750 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1751 PetscFunctionReturn(0); 1752 } 1753 1754 static PetscErrorCode MatDestroy_IS(Mat A) 1755 { 1756 PetscErrorCode ierr; 1757 Mat_IS *b = (Mat_IS*)A->data; 1758 1759 PetscFunctionBegin; 1760 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 1761 ierr = VecScatterDestroy(&b->cctx);CHKERRQ(ierr); 1762 ierr = VecScatterDestroy(&b->rctx);CHKERRQ(ierr); 1763 ierr = VecDestroy(&b->x);CHKERRQ(ierr); 1764 ierr = VecDestroy(&b->y);CHKERRQ(ierr); 1765 ierr = VecDestroy(&b->counter);CHKERRQ(ierr); 1766 ierr = ISDestroy(&b->getsub_ris);CHKERRQ(ierr); 1767 ierr = ISDestroy(&b->getsub_cis);CHKERRQ(ierr); 1768 if (b->sf != b->csf) { 1769 ierr = PetscSFDestroy(&b->csf);CHKERRQ(ierr); 1770 ierr = PetscFree2(b->csf_rootdata,b->csf_leafdata);CHKERRQ(ierr); 1771 } 1772 ierr = PetscSFDestroy(&b->sf);CHKERRQ(ierr); 1773 ierr = PetscFree2(b->sf_rootdata,b->sf_leafdata);CHKERRQ(ierr); 1774 ierr = PetscFree(A->data);CHKERRQ(ierr); 1775 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 1776 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);CHKERRQ(ierr); 1777 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);CHKERRQ(ierr); 1778 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);CHKERRQ(ierr); 1779 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);CHKERRQ(ierr); 1780 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",NULL);CHKERRQ(ierr); 1781 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL);CHKERRQ(ierr); 1782 PetscFunctionReturn(0); 1783 } 1784 1785 static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y) 1786 { 1787 PetscErrorCode ierr; 1788 Mat_IS *is = (Mat_IS*)A->data; 1789 PetscScalar zero = 0.0; 1790 1791 PetscFunctionBegin; 1792 /* scatter the global vector x into the local work vector */ 1793 ierr = VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1794 ierr = VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1795 1796 /* multiply the local matrix */ 1797 ierr = MatMult(is->A,is->x,is->y);CHKERRQ(ierr); 1798 1799 /* scatter product back into global memory */ 1800 ierr = VecSet(y,zero);CHKERRQ(ierr); 1801 ierr = VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1802 ierr = VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1803 PetscFunctionReturn(0); 1804 } 1805 1806 static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 1807 { 1808 Vec temp_vec; 1809 PetscErrorCode ierr; 1810 1811 PetscFunctionBegin; /* v3 = v2 + A * v1.*/ 1812 if (v3 != v2) { 1813 ierr = MatMult(A,v1,v3);CHKERRQ(ierr); 1814 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 1815 } else { 1816 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 1817 ierr = MatMult(A,v1,temp_vec);CHKERRQ(ierr); 1818 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 1819 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 1820 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1821 } 1822 PetscFunctionReturn(0); 1823 } 1824 1825 static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x) 1826 { 1827 Mat_IS *is = (Mat_IS*)A->data; 1828 PetscErrorCode ierr; 1829 1830 PetscFunctionBegin; 1831 /* scatter the global vector x into the local work vector */ 1832 ierr = VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1833 ierr = VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1834 1835 /* multiply the local matrix */ 1836 ierr = MatMultTranspose(is->A,is->y,is->x);CHKERRQ(ierr); 1837 1838 /* scatter product back into global vector */ 1839 ierr = VecSet(x,0);CHKERRQ(ierr); 1840 ierr = VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1841 ierr = VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1842 PetscFunctionReturn(0); 1843 } 1844 1845 static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 1846 { 1847 Vec temp_vec; 1848 PetscErrorCode ierr; 1849 1850 PetscFunctionBegin; /* v3 = v2 + A' * v1.*/ 1851 if (v3 != v2) { 1852 ierr = MatMultTranspose(A,v1,v3);CHKERRQ(ierr); 1853 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 1854 } else { 1855 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 1856 ierr = MatMultTranspose(A,v1,temp_vec);CHKERRQ(ierr); 1857 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 1858 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 1859 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1860 } 1861 PetscFunctionReturn(0); 1862 } 1863 1864 static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer) 1865 { 1866 Mat_IS *a = (Mat_IS*)A->data; 1867 PetscErrorCode ierr; 1868 PetscViewer sviewer; 1869 PetscBool isascii,view = PETSC_TRUE; 1870 1871 PetscFunctionBegin; 1872 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 1873 if (isascii) { 1874 PetscViewerFormat format; 1875 1876 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1877 if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE; 1878 } 1879 if (!view) PetscFunctionReturn(0); 1880 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 1881 ierr = MatView(a->A,sviewer);CHKERRQ(ierr); 1882 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 1883 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1884 PetscFunctionReturn(0); 1885 } 1886 1887 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping) 1888 { 1889 PetscErrorCode ierr; 1890 PetscInt nr,rbs,nc,cbs; 1891 Mat_IS *is = (Mat_IS*)A->data; 1892 IS from,to; 1893 Vec cglobal,rglobal; 1894 1895 PetscFunctionBegin; 1896 PetscCheckSameComm(A,1,rmapping,2); 1897 PetscCheckSameComm(A,1,cmapping,3); 1898 /* Destroy any previous data */ 1899 ierr = VecDestroy(&is->x);CHKERRQ(ierr); 1900 ierr = VecDestroy(&is->y);CHKERRQ(ierr); 1901 ierr = VecDestroy(&is->counter);CHKERRQ(ierr); 1902 ierr = VecScatterDestroy(&is->rctx);CHKERRQ(ierr); 1903 ierr = VecScatterDestroy(&is->cctx);CHKERRQ(ierr); 1904 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 1905 if (is->csf != is->sf) { 1906 ierr = PetscSFDestroy(&is->csf);CHKERRQ(ierr); 1907 ierr = PetscFree2(is->csf_rootdata,is->csf_leafdata);CHKERRQ(ierr); 1908 } 1909 ierr = PetscSFDestroy(&is->sf);CHKERRQ(ierr); 1910 ierr = PetscFree2(is->sf_rootdata,is->sf_leafdata);CHKERRQ(ierr); 1911 1912 /* Setup Layout and set local to global maps */ 1913 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 1914 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 1915 ierr = ISLocalToGlobalMappingGetSize(rmapping,&nr);CHKERRQ(ierr); 1916 ierr = ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);CHKERRQ(ierr); 1917 ierr = ISLocalToGlobalMappingGetSize(cmapping,&nc);CHKERRQ(ierr); 1918 ierr = ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);CHKERRQ(ierr); 1919 /* check if the two mappings are actually the same for square matrices (DOLFIN passes 2 different objects) */ 1920 if (rmapping != cmapping && A->rmap->N == A->cmap->N) { 1921 PetscBool same,gsame; 1922 1923 same = PETSC_FALSE; 1924 if (nr == nc && cbs == rbs) { 1925 const PetscInt *idxs1,*idxs2; 1926 1927 ierr = ISLocalToGlobalMappingGetBlockIndices(rmapping,&idxs1);CHKERRQ(ierr); 1928 ierr = ISLocalToGlobalMappingGetBlockIndices(cmapping,&idxs2);CHKERRQ(ierr); 1929 ierr = PetscMemcmp(idxs1,idxs2,(nr/rbs)*sizeof(PetscInt),&same);CHKERRQ(ierr); 1930 ierr = ISLocalToGlobalMappingRestoreBlockIndices(rmapping,&idxs1);CHKERRQ(ierr); 1931 ierr = ISLocalToGlobalMappingRestoreBlockIndices(cmapping,&idxs2);CHKERRQ(ierr); 1932 } 1933 ierr = MPIU_Allreduce(&same,&gsame,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 1934 if (gsame) cmapping = rmapping; 1935 } 1936 ierr = PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);CHKERRQ(ierr); 1937 ierr = PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);CHKERRQ(ierr); 1938 1939 /* Create the local matrix A */ 1940 ierr = MatCreate(PETSC_COMM_SELF,&is->A);CHKERRQ(ierr); 1941 ierr = MatSetType(is->A,MATAIJ);CHKERRQ(ierr); 1942 ierr = MatSetSizes(is->A,nr,nc,nr,nc);CHKERRQ(ierr); 1943 ierr = MatSetBlockSizes(is->A,rbs,cbs);CHKERRQ(ierr); 1944 ierr = MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);CHKERRQ(ierr); 1945 ierr = MatAppendOptionsPrefix(is->A,"is_");CHKERRQ(ierr); 1946 ierr = MatSetFromOptions(is->A);CHKERRQ(ierr); 1947 ierr = PetscLayoutSetUp(is->A->rmap);CHKERRQ(ierr); 1948 ierr = PetscLayoutSetUp(is->A->cmap);CHKERRQ(ierr); 1949 1950 if (!is->islocalref) { /* setup scatters and local vectors for MatMult */ 1951 /* Create the local work vectors */ 1952 ierr = MatCreateVecs(is->A,&is->x,&is->y);CHKERRQ(ierr); 1953 1954 /* setup the global to local scatters */ 1955 ierr = MatCreateVecs(A,&cglobal,&rglobal);CHKERRQ(ierr); 1956 ierr = ISCreateStride(PETSC_COMM_SELF,nr,0,1,&to);CHKERRQ(ierr); 1957 ierr = ISLocalToGlobalMappingApplyIS(rmapping,to,&from);CHKERRQ(ierr); 1958 ierr = VecScatterCreate(rglobal,from,is->y,to,&is->rctx);CHKERRQ(ierr); 1959 if (rmapping != cmapping) { 1960 ierr = ISDestroy(&to);CHKERRQ(ierr); 1961 ierr = ISDestroy(&from);CHKERRQ(ierr); 1962 ierr = ISCreateStride(PETSC_COMM_SELF,nc,0,1,&to);CHKERRQ(ierr); 1963 ierr = ISLocalToGlobalMappingApplyIS(cmapping,to,&from);CHKERRQ(ierr); 1964 ierr = VecScatterCreate(cglobal,from,is->x,to,&is->cctx);CHKERRQ(ierr); 1965 } else { 1966 ierr = PetscObjectReference((PetscObject)is->rctx);CHKERRQ(ierr); 1967 is->cctx = is->rctx; 1968 } 1969 1970 /* interface counter vector (local) */ 1971 ierr = VecDuplicate(is->y,&is->counter);CHKERRQ(ierr); 1972 ierr = VecSet(is->y,1.);CHKERRQ(ierr); 1973 ierr = VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1974 ierr = VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1975 ierr = VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1976 ierr = VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1977 1978 /* free workspace */ 1979 ierr = VecDestroy(&rglobal);CHKERRQ(ierr); 1980 ierr = VecDestroy(&cglobal);CHKERRQ(ierr); 1981 ierr = ISDestroy(&to);CHKERRQ(ierr); 1982 ierr = ISDestroy(&from);CHKERRQ(ierr); 1983 } 1984 ierr = MatSetUp(A);CHKERRQ(ierr); 1985 PetscFunctionReturn(0); 1986 } 1987 1988 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 1989 { 1990 Mat_IS *is = (Mat_IS*)mat->data; 1991 PetscErrorCode ierr; 1992 #if defined(PETSC_USE_DEBUG) 1993 PetscInt i,zm,zn; 1994 #endif 1995 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 1996 1997 PetscFunctionBegin; 1998 #if defined(PETSC_USE_DEBUG) 1999 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); 2000 /* count negative indices */ 2001 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 2002 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 2003 #endif 2004 ierr = ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 2005 ierr = ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 2006 #if defined(PETSC_USE_DEBUG) 2007 /* count negative indices (should be the same as before) */ 2008 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 2009 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 2010 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"); 2011 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"); 2012 #endif 2013 ierr = MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 2014 PetscFunctionReturn(0); 2015 } 2016 2017 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2018 { 2019 Mat_IS *is = (Mat_IS*)mat->data; 2020 PetscErrorCode ierr; 2021 #if defined(PETSC_USE_DEBUG) 2022 PetscInt i,zm,zn; 2023 #endif 2024 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2025 2026 PetscFunctionBegin; 2027 #if defined(PETSC_USE_DEBUG) 2028 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); 2029 /* count negative indices */ 2030 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 2031 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 2032 #endif 2033 ierr = ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 2034 ierr = ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 2035 #if defined(PETSC_USE_DEBUG) 2036 /* count negative indices (should be the same as before) */ 2037 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 2038 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 2039 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"); 2040 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"); 2041 #endif 2042 ierr = MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 2043 PetscFunctionReturn(0); 2044 } 2045 2046 static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2047 { 2048 PetscErrorCode ierr; 2049 Mat_IS *is = (Mat_IS*)A->data; 2050 2051 PetscFunctionBegin; 2052 if (is->A->rmap->mapping) { 2053 ierr = MatSetValuesLocal(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2054 } else { 2055 ierr = MatSetValues(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2056 } 2057 PetscFunctionReturn(0); 2058 } 2059 2060 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2061 { 2062 PetscErrorCode ierr; 2063 Mat_IS *is = (Mat_IS*)A->data; 2064 2065 PetscFunctionBegin; 2066 if (is->A->rmap->mapping) { 2067 #if defined(PETSC_USE_DEBUG) 2068 PetscInt ibs,bs; 2069 2070 ierr = ISLocalToGlobalMappingGetBlockSize(is->A->rmap->mapping,&ibs);CHKERRQ(ierr); 2071 ierr = MatGetBlockSize(is->A,&bs);CHKERRQ(ierr); 2072 if (ibs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Different block sizes! local mat %D, local l2g map %D",bs,ibs); 2073 #endif 2074 ierr = MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2075 } else { 2076 ierr = MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2077 } 2078 PetscFunctionReturn(0); 2079 } 2080 2081 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns) 2082 { 2083 Mat_IS *is = (Mat_IS*)A->data; 2084 PetscErrorCode ierr; 2085 2086 PetscFunctionBegin; 2087 if (!n) { 2088 is->pure_neumann = PETSC_TRUE; 2089 } else { 2090 PetscInt i; 2091 is->pure_neumann = PETSC_FALSE; 2092 2093 if (columns) { 2094 ierr = MatZeroRowsColumns(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 2095 } else { 2096 ierr = MatZeroRows(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 2097 } 2098 if (diag != 0.) { 2099 const PetscScalar *array; 2100 ierr = VecGetArrayRead(is->counter,&array);CHKERRQ(ierr); 2101 for (i=0; i<n; i++) { 2102 ierr = MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);CHKERRQ(ierr); 2103 } 2104 ierr = VecRestoreArrayRead(is->counter,&array);CHKERRQ(ierr); 2105 } 2106 ierr = MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2107 ierr = MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2108 } 2109 PetscFunctionReturn(0); 2110 } 2111 2112 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns) 2113 { 2114 Mat_IS *matis = (Mat_IS*)A->data; 2115 PetscInt nr,nl,len,i; 2116 PetscInt *lrows; 2117 PetscErrorCode ierr; 2118 2119 PetscFunctionBegin; 2120 #if defined(PETSC_USE_DEBUG) 2121 if (columns || diag != 0. || (x && b)) { 2122 PetscBool cong; 2123 2124 ierr = PetscLayoutCompare(A->rmap,A->cmap,&cong);CHKERRQ(ierr); 2125 cong = (PetscBool)(cong && matis->sf == matis->csf); 2126 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"); 2127 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"); 2128 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"); 2129 } 2130 #endif 2131 /* get locally owned rows */ 2132 ierr = PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);CHKERRQ(ierr); 2133 /* fix right hand side if needed */ 2134 if (x && b) { 2135 const PetscScalar *xx; 2136 PetscScalar *bb; 2137 2138 ierr = VecGetArrayRead(x, &xx);CHKERRQ(ierr); 2139 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 2140 for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]]; 2141 ierr = VecRestoreArrayRead(x, &xx);CHKERRQ(ierr); 2142 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 2143 } 2144 /* get rows associated to the local matrices */ 2145 ierr = MatISSetUpSF(A);CHKERRQ(ierr); 2146 ierr = MatGetSize(matis->A,&nl,NULL);CHKERRQ(ierr); 2147 ierr = PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 2148 ierr = PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2149 for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1; 2150 ierr = PetscFree(lrows);CHKERRQ(ierr); 2151 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 2152 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 2153 ierr = PetscMalloc1(nl,&lrows);CHKERRQ(ierr); 2154 for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i; 2155 ierr = MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);CHKERRQ(ierr); 2156 ierr = PetscFree(lrows);CHKERRQ(ierr); 2157 PetscFunctionReturn(0); 2158 } 2159 2160 static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2161 { 2162 PetscErrorCode ierr; 2163 2164 PetscFunctionBegin; 2165 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);CHKERRQ(ierr); 2166 PetscFunctionReturn(0); 2167 } 2168 2169 static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2170 { 2171 PetscErrorCode ierr; 2172 2173 PetscFunctionBegin; 2174 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);CHKERRQ(ierr); 2175 PetscFunctionReturn(0); 2176 } 2177 2178 static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type) 2179 { 2180 Mat_IS *is = (Mat_IS*)A->data; 2181 PetscErrorCode ierr; 2182 2183 PetscFunctionBegin; 2184 ierr = MatAssemblyBegin(is->A,type);CHKERRQ(ierr); 2185 PetscFunctionReturn(0); 2186 } 2187 2188 static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type) 2189 { 2190 Mat_IS *is = (Mat_IS*)A->data; 2191 PetscErrorCode ierr; 2192 2193 PetscFunctionBegin; 2194 ierr = MatAssemblyEnd(is->A,type);CHKERRQ(ierr); 2195 /* fix for local empty rows/cols */ 2196 if (is->locempty && type == MAT_FINAL_ASSEMBLY) { 2197 Mat newlA; 2198 ISLocalToGlobalMapping l2g; 2199 IS tis; 2200 const PetscScalar *v; 2201 PetscInt i,n,cf,*loce,*locf; 2202 PetscBool sym; 2203 2204 ierr = MatGetRowMaxAbs(is->A,is->y,NULL);CHKERRQ(ierr); 2205 ierr = MatIsSymmetric(is->A,PETSC_SMALL,&sym);CHKERRQ(ierr); 2206 if (!sym) SETERRQ(PetscObjectComm((PetscObject)is->A),PETSC_ERR_SUP,"Not yet implemented for unsymmetric case"); 2207 ierr = VecGetLocalSize(is->y,&n);CHKERRQ(ierr); 2208 ierr = PetscMalloc1(n,&loce);CHKERRQ(ierr); 2209 ierr = PetscMalloc1(n,&locf);CHKERRQ(ierr); 2210 ierr = VecGetArrayRead(is->y,&v);CHKERRQ(ierr); 2211 for (i=0,cf=0;i<n;i++) { 2212 if (v[i] == 0.0) { 2213 loce[i] = -1; 2214 } else { 2215 loce[i] = cf; 2216 locf[cf++] = i; 2217 } 2218 } 2219 ierr = VecRestoreArrayRead(is->y,&v);CHKERRQ(ierr); 2220 /* extract valid submatrix */ 2221 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is->A),cf,locf,PETSC_USE_POINTER,&tis);CHKERRQ(ierr); 2222 ierr = MatCreateSubMatrix(is->A,tis,tis,MAT_INITIAL_MATRIX,&newlA);CHKERRQ(ierr); 2223 ierr = ISDestroy(&tis);CHKERRQ(ierr); 2224 /* attach local l2g map for successive calls of MatSetValues */ 2225 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is->A),1,n,loce,PETSC_OWN_POINTER,&l2g);CHKERRQ(ierr); 2226 ierr = MatSetLocalToGlobalMapping(newlA,l2g,l2g);CHKERRQ(ierr); 2227 ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr); 2228 /* attach new global l2g map */ 2229 ierr = ISLocalToGlobalMappingApply(A->rmap->mapping,cf,locf,locf);CHKERRQ(ierr); 2230 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,cf,locf,PETSC_OWN_POINTER,&l2g);CHKERRQ(ierr); 2231 ierr = MatSetLocalToGlobalMapping(A,l2g,l2g);CHKERRQ(ierr); 2232 ierr = MatISSetLocalMat(A,newlA);CHKERRQ(ierr); 2233 ierr = MatDestroy(&newlA);CHKERRQ(ierr); 2234 ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr); 2235 } 2236 is->locempty = PETSC_FALSE; 2237 PetscFunctionReturn(0); 2238 } 2239 2240 static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local) 2241 { 2242 Mat_IS *is = (Mat_IS*)mat->data; 2243 2244 PetscFunctionBegin; 2245 *local = is->A; 2246 PetscFunctionReturn(0); 2247 } 2248 2249 static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local) 2250 { 2251 PetscFunctionBegin; 2252 *local = NULL; 2253 PetscFunctionReturn(0); 2254 } 2255 2256 /*@ 2257 MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix. 2258 2259 Input Parameter: 2260 . mat - the matrix 2261 2262 Output Parameter: 2263 . local - the local matrix 2264 2265 Level: advanced 2266 2267 Notes: 2268 This can be called if you have precomputed the nonzero structure of the 2269 matrix and want to provide it to the inner matrix object to improve the performance 2270 of the MatSetValues() operation. 2271 2272 Call MatISRestoreLocalMat() when finished with the local matrix. 2273 2274 .seealso: MATIS 2275 @*/ 2276 PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local) 2277 { 2278 PetscErrorCode ierr; 2279 2280 PetscFunctionBegin; 2281 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2282 PetscValidPointer(local,2); 2283 ierr = PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 2284 PetscFunctionReturn(0); 2285 } 2286 2287 /*@ 2288 MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat() 2289 2290 Input Parameter: 2291 . mat - the matrix 2292 2293 Output Parameter: 2294 . local - the local matrix 2295 2296 Level: advanced 2297 2298 .seealso: MATIS 2299 @*/ 2300 PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local) 2301 { 2302 PetscErrorCode ierr; 2303 2304 PetscFunctionBegin; 2305 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2306 PetscValidPointer(local,2); 2307 ierr = PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 2308 PetscFunctionReturn(0); 2309 } 2310 2311 static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local) 2312 { 2313 Mat_IS *is = (Mat_IS*)mat->data; 2314 PetscInt nrows,ncols,orows,ocols; 2315 PetscErrorCode ierr; 2316 2317 PetscFunctionBegin; 2318 if (is->A) { 2319 ierr = MatGetSize(is->A,&orows,&ocols);CHKERRQ(ierr); 2320 ierr = MatGetSize(local,&nrows,&ncols);CHKERRQ(ierr); 2321 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); 2322 } 2323 ierr = PetscObjectReference((PetscObject)local);CHKERRQ(ierr); 2324 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 2325 is->A = local; 2326 PetscFunctionReturn(0); 2327 } 2328 2329 /*@ 2330 MatISSetLocalMat - Replace the local matrix stored inside a MATIS object. 2331 2332 Input Parameter: 2333 . mat - the matrix 2334 . local - the local matrix 2335 2336 Output Parameter: 2337 2338 Level: advanced 2339 2340 Notes: 2341 This can be called if you have precomputed the local matrix and 2342 want to provide it to the matrix object MATIS. 2343 2344 .seealso: MATIS 2345 @*/ 2346 PetscErrorCode MatISSetLocalMat(Mat mat,Mat local) 2347 { 2348 PetscErrorCode ierr; 2349 2350 PetscFunctionBegin; 2351 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2352 PetscValidHeaderSpecific(local,MAT_CLASSID,2); 2353 ierr = PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));CHKERRQ(ierr); 2354 PetscFunctionReturn(0); 2355 } 2356 2357 static PetscErrorCode MatZeroEntries_IS(Mat A) 2358 { 2359 Mat_IS *a = (Mat_IS*)A->data; 2360 PetscErrorCode ierr; 2361 2362 PetscFunctionBegin; 2363 ierr = MatZeroEntries(a->A);CHKERRQ(ierr); 2364 PetscFunctionReturn(0); 2365 } 2366 2367 static PetscErrorCode MatScale_IS(Mat A,PetscScalar a) 2368 { 2369 Mat_IS *is = (Mat_IS*)A->data; 2370 PetscErrorCode ierr; 2371 2372 PetscFunctionBegin; 2373 ierr = MatScale(is->A,a);CHKERRQ(ierr); 2374 PetscFunctionReturn(0); 2375 } 2376 2377 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v) 2378 { 2379 Mat_IS *is = (Mat_IS*)A->data; 2380 PetscErrorCode ierr; 2381 2382 PetscFunctionBegin; 2383 /* get diagonal of the local matrix */ 2384 ierr = MatGetDiagonal(is->A,is->y);CHKERRQ(ierr); 2385 2386 /* scatter diagonal back into global vector */ 2387 ierr = VecSet(v,0);CHKERRQ(ierr); 2388 ierr = VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2389 ierr = VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2390 PetscFunctionReturn(0); 2391 } 2392 2393 static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg) 2394 { 2395 Mat_IS *a = (Mat_IS*)A->data; 2396 PetscErrorCode ierr; 2397 2398 PetscFunctionBegin; 2399 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 2400 PetscFunctionReturn(0); 2401 } 2402 2403 static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str) 2404 { 2405 Mat_IS *y = (Mat_IS*)Y->data; 2406 Mat_IS *x; 2407 #if defined(PETSC_USE_DEBUG) 2408 PetscBool ismatis; 2409 #endif 2410 PetscErrorCode ierr; 2411 2412 PetscFunctionBegin; 2413 #if defined(PETSC_USE_DEBUG) 2414 ierr = PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);CHKERRQ(ierr); 2415 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS"); 2416 #endif 2417 x = (Mat_IS*)X->data; 2418 ierr = MatAXPY(y->A,a,x->A,str);CHKERRQ(ierr); 2419 PetscFunctionReturn(0); 2420 } 2421 2422 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat) 2423 { 2424 Mat lA; 2425 Mat_IS *matis; 2426 ISLocalToGlobalMapping rl2g,cl2g; 2427 IS is; 2428 const PetscInt *rg,*rl; 2429 PetscInt nrg; 2430 PetscInt N,M,nrl,i,*idxs; 2431 PetscErrorCode ierr; 2432 2433 PetscFunctionBegin; 2434 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 2435 ierr = ISGetLocalSize(row,&nrl);CHKERRQ(ierr); 2436 ierr = ISGetIndices(row,&rl);CHKERRQ(ierr); 2437 ierr = ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);CHKERRQ(ierr); 2438 #if defined(PETSC_USE_DEBUG) 2439 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); 2440 #endif 2441 ierr = PetscMalloc1(nrg,&idxs);CHKERRQ(ierr); 2442 /* map from [0,nrl) to row */ 2443 for (i=0;i<nrl;i++) idxs[i] = rl[i]; 2444 #if defined(PETSC_USE_DEBUG) 2445 for (i=nrl;i<nrg;i++) idxs[i] = nrg; 2446 #else 2447 for (i=nrl;i<nrg;i++) idxs[i] = -1; 2448 #endif 2449 ierr = ISRestoreIndices(row,&rl);CHKERRQ(ierr); 2450 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 2451 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 2452 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 2453 ierr = ISDestroy(&is);CHKERRQ(ierr); 2454 /* compute new l2g map for columns */ 2455 if (col != row || A->rmap->mapping != A->cmap->mapping) { 2456 const PetscInt *cg,*cl; 2457 PetscInt ncg; 2458 PetscInt ncl; 2459 2460 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 2461 ierr = ISGetLocalSize(col,&ncl);CHKERRQ(ierr); 2462 ierr = ISGetIndices(col,&cl);CHKERRQ(ierr); 2463 ierr = ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);CHKERRQ(ierr); 2464 #if defined(PETSC_USE_DEBUG) 2465 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); 2466 #endif 2467 ierr = PetscMalloc1(ncg,&idxs);CHKERRQ(ierr); 2468 /* map from [0,ncl) to col */ 2469 for (i=0;i<ncl;i++) idxs[i] = cl[i]; 2470 #if defined(PETSC_USE_DEBUG) 2471 for (i=ncl;i<ncg;i++) idxs[i] = ncg; 2472 #else 2473 for (i=ncl;i<ncg;i++) idxs[i] = -1; 2474 #endif 2475 ierr = ISRestoreIndices(col,&cl);CHKERRQ(ierr); 2476 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 2477 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 2478 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 2479 ierr = ISDestroy(&is);CHKERRQ(ierr); 2480 } else { 2481 ierr = PetscObjectReference((PetscObject)rl2g);CHKERRQ(ierr); 2482 cl2g = rl2g; 2483 } 2484 /* create the MATIS submatrix */ 2485 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 2486 ierr = MatCreate(PetscObjectComm((PetscObject)A),submat);CHKERRQ(ierr); 2487 ierr = MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2488 ierr = MatSetType(*submat,MATIS);CHKERRQ(ierr); 2489 matis = (Mat_IS*)((*submat)->data); 2490 matis->islocalref = PETSC_TRUE; 2491 ierr = MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);CHKERRQ(ierr); 2492 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr); 2493 ierr = MatISSetLocalMat(*submat,lA);CHKERRQ(ierr); 2494 ierr = MatSetUp(*submat);CHKERRQ(ierr); 2495 ierr = MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2496 ierr = MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2497 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 2498 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 2499 /* remove unsupported ops */ 2500 ierr = PetscMemzero((*submat)->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2501 (*submat)->ops->destroy = MatDestroy_IS; 2502 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS; 2503 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS; 2504 (*submat)->ops->assemblybegin = MatAssemblyBegin_IS; 2505 (*submat)->ops->assemblyend = MatAssemblyEnd_IS; 2506 PetscFunctionReturn(0); 2507 } 2508 2509 static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A) 2510 { 2511 Mat_IS *a = (Mat_IS*)A->data; 2512 PetscErrorCode ierr; 2513 2514 PetscFunctionBegin; 2515 ierr = PetscOptionsHead(PetscOptionsObject,"MATIS options");CHKERRQ(ierr); 2516 ierr = PetscObjectOptionsBegin((PetscObject)A); 2517 ierr = PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns",NULL,a->locempty,&a->locempty,NULL);CHKERRQ(ierr); 2518 ierr = PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL);CHKERRQ(ierr); 2519 ierr = PetscOptionsEnd();CHKERRQ(ierr); 2520 PetscFunctionReturn(0); 2521 } 2522 2523 /*@ 2524 MatCreateIS - Creates a "process" unassembled matrix, assembled on each 2525 process but not across processes. 2526 2527 Input Parameters: 2528 + comm - MPI communicator that will share the matrix 2529 . bs - block size of the matrix 2530 . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products 2531 . rmap - local to global map for rows 2532 - cmap - local to global map for cols 2533 2534 Output Parameter: 2535 . A - the resulting matrix 2536 2537 Level: advanced 2538 2539 Notes: 2540 See MATIS for more details. 2541 m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors 2542 used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices. 2543 If either rmap or cmap are NULL, then the matrix is assumed to be square. 2544 2545 .seealso: MATIS, MatSetLocalToGlobalMapping() 2546 @*/ 2547 PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A) 2548 { 2549 PetscErrorCode ierr; 2550 2551 PetscFunctionBegin; 2552 if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings"); 2553 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2554 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 2555 if (bs > 0) { 2556 ierr = MatSetBlockSize(*A,bs);CHKERRQ(ierr); 2557 } 2558 ierr = MatSetType(*A,MATIS);CHKERRQ(ierr); 2559 if (rmap && cmap) { 2560 ierr = MatSetLocalToGlobalMapping(*A,rmap,cmap);CHKERRQ(ierr); 2561 } else if (!rmap) { 2562 ierr = MatSetLocalToGlobalMapping(*A,cmap,cmap);CHKERRQ(ierr); 2563 } else { 2564 ierr = MatSetLocalToGlobalMapping(*A,rmap,rmap);CHKERRQ(ierr); 2565 } 2566 PetscFunctionReturn(0); 2567 } 2568 2569 /*MC 2570 MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP). 2571 This stores the matrices in globally unassembled form. Each processor 2572 assembles only its local Neumann problem and the parallel matrix vector 2573 product is handled "implicitly". 2574 2575 Operations Provided: 2576 + MatMult() 2577 . MatMultAdd() 2578 . MatMultTranspose() 2579 . MatMultTransposeAdd() 2580 . MatZeroEntries() 2581 . MatSetOption() 2582 . MatZeroRows() 2583 . MatSetValues() 2584 . MatSetValuesBlocked() 2585 . MatSetValuesLocal() 2586 . MatSetValuesBlockedLocal() 2587 . MatScale() 2588 . MatGetDiagonal() 2589 . MatMissingDiagonal() 2590 . MatDuplicate() 2591 . MatCopy() 2592 . MatAXPY() 2593 . MatCreateSubMatrix() 2594 . MatGetLocalSubMatrix() 2595 . MatTranspose() 2596 . MatPtAP() (with P of AIJ type) 2597 - MatSetLocalToGlobalMapping() 2598 2599 Options Database Keys: 2600 + -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions() 2601 . -matis_fixempty - Fixes local matrices in case of empty local rows/columns. 2602 - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of MatPtAP(). 2603 2604 Notes: 2605 Options prefix for the inner matrix are given by -is_mat_xxx 2606 2607 You must call MatSetLocalToGlobalMapping() before using this matrix type. 2608 2609 You can do matrix preallocation on the local matrix after you obtain it with 2610 MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation() 2611 2612 Level: advanced 2613 2614 .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP 2615 2616 M*/ 2617 2618 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A) 2619 { 2620 PetscErrorCode ierr; 2621 Mat_IS *b; 2622 2623 PetscFunctionBegin; 2624 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 2625 A->data = (void*)b; 2626 2627 /* matrix ops */ 2628 ierr = PetscMemzero(A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2629 A->ops->mult = MatMult_IS; 2630 A->ops->multadd = MatMultAdd_IS; 2631 A->ops->multtranspose = MatMultTranspose_IS; 2632 A->ops->multtransposeadd = MatMultTransposeAdd_IS; 2633 A->ops->destroy = MatDestroy_IS; 2634 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS; 2635 A->ops->setvalues = MatSetValues_IS; 2636 A->ops->setvaluesblocked = MatSetValuesBlocked_IS; 2637 A->ops->setvalueslocal = MatSetValuesLocal_IS; 2638 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS; 2639 A->ops->zerorows = MatZeroRows_IS; 2640 A->ops->zerorowscolumns = MatZeroRowsColumns_IS; 2641 A->ops->assemblybegin = MatAssemblyBegin_IS; 2642 A->ops->assemblyend = MatAssemblyEnd_IS; 2643 A->ops->view = MatView_IS; 2644 A->ops->zeroentries = MatZeroEntries_IS; 2645 A->ops->scale = MatScale_IS; 2646 A->ops->getdiagonal = MatGetDiagonal_IS; 2647 A->ops->setoption = MatSetOption_IS; 2648 A->ops->ishermitian = MatIsHermitian_IS; 2649 A->ops->issymmetric = MatIsSymmetric_IS; 2650 A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS; 2651 A->ops->duplicate = MatDuplicate_IS; 2652 A->ops->missingdiagonal = MatMissingDiagonal_IS; 2653 A->ops->copy = MatCopy_IS; 2654 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 2655 A->ops->createsubmatrix = MatCreateSubMatrix_IS; 2656 A->ops->axpy = MatAXPY_IS; 2657 A->ops->diagonalset = MatDiagonalSet_IS; 2658 A->ops->shift = MatShift_IS; 2659 A->ops->transpose = MatTranspose_IS; 2660 A->ops->getinfo = MatGetInfo_IS; 2661 A->ops->diagonalscale = MatDiagonalScale_IS; 2662 A->ops->setfromoptions = MatSetFromOptions_IS; 2663 2664 /* special MATIS functions */ 2665 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);CHKERRQ(ierr); 2666 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);CHKERRQ(ierr); 2667 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);CHKERRQ(ierr); 2668 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatISGetMPIXAIJ_IS);CHKERRQ(ierr); 2669 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);CHKERRQ(ierr); 2670 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",MatISSetUpSF_IS);CHKERRQ(ierr); 2671 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS);CHKERRQ(ierr); 2672 ierr = PetscObjectChangeTypeName((PetscObject)A,MATIS);CHKERRQ(ierr); 2673 PetscFunctionReturn(0); 2674 } 2675