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