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