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 /* Resize preallocation if overestimated */ 1798 for (i=0;i<lrows;i++) { 1799 dnz[i] = PetscMin(dnz[i],lcols); 1800 onz[i] = PetscMin(onz[i],cols-lcols); 1801 } 1802 1803 /* Set preallocation */ 1804 ierr = MatSeqAIJSetPreallocation(B,0,dnz);CHKERRQ(ierr); 1805 ierr = MatMPIAIJSetPreallocation(B,0,dnz,0,onz);CHKERRQ(ierr); 1806 for (i=0;i<lrows/bs;i++) { 1807 dnz[i] = dnz[i*bs]/bs; 1808 onz[i] = onz[i*bs]/bs; 1809 } 1810 ierr = MatSeqBAIJSetPreallocation(B,bs,0,dnz);CHKERRQ(ierr); 1811 ierr = MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 1812 ierr = MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 1813 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1814 if (issbaij) { 1815 ierr = MatRestoreRowUpperTriangular(matis->A);CHKERRQ(ierr); 1816 } 1817 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1818 PetscFunctionReturn(0); 1819 } 1820 1821 PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M) 1822 { 1823 Mat_IS *matis = (Mat_IS*)(mat->data); 1824 Mat local_mat,MT; 1825 /* info on mat */ 1826 PetscInt bs,rows,cols,lrows,lcols; 1827 PetscInt local_rows,local_cols; 1828 PetscBool isseqdense,isseqsbaij,isseqaij,isseqbaij; 1829 #if defined (PETSC_USE_DEBUG) 1830 PetscBool lb[4],bb[4]; 1831 #endif 1832 PetscMPIInt size; 1833 /* values insertion */ 1834 PetscScalar *array; 1835 /* work */ 1836 PetscErrorCode ierr; 1837 1838 PetscFunctionBegin; 1839 /* get info from mat */ 1840 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);CHKERRQ(ierr); 1841 if (size == 1 && mat->rmap->N == matis->A->rmap->N && mat->cmap->N == matis->A->cmap->N) { 1842 Mat B; 1843 IS rows,cols; 1844 IS irows,icols; 1845 const PetscInt *ridxs,*cidxs; 1846 PetscInt rbs,cbs; 1847 1848 ierr = ISLocalToGlobalMappingGetBlockIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 1849 ierr = ISLocalToGlobalMappingGetBlockIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 1850 ierr = ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);CHKERRQ(ierr); 1851 ierr = ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);CHKERRQ(ierr); 1852 ierr = ISCreateBlock(PETSC_COMM_SELF,rbs,mat->rmap->n,ridxs,PETSC_USE_POINTER,&rows);CHKERRQ(ierr); 1853 ierr = ISCreateBlock(PETSC_COMM_SELF,cbs,mat->cmap->n,cidxs,PETSC_USE_POINTER,&cols);CHKERRQ(ierr); 1854 ierr = ISSetPermutation(rows);CHKERRQ(ierr); 1855 ierr = ISSetPermutation(cols);CHKERRQ(ierr); 1856 ierr = ISInvertPermutation(rows,mat->rmap->n,&irows);CHKERRQ(ierr); 1857 ierr = ISInvertPermutation(cols,mat->cmap->n,&icols);CHKERRQ(ierr); 1858 ierr = ISLocalToGlobalMappingRestoreBlockIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 1859 ierr = ISLocalToGlobalMappingRestoreBlockIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 1860 ierr = ISDestroy(&cols);CHKERRQ(ierr); 1861 ierr = ISDestroy(&rows);CHKERRQ(ierr); 1862 ierr = MatConvert(matis->A,(rbs == cbs && rbs > 1) ? MATSEQBAIJ : MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 1863 if (reuse != MAT_INPLACE_MATRIX) { 1864 ierr = MatCreateSubMatrix(B,irows,icols,reuse,M);CHKERRQ(ierr); 1865 } else { 1866 Mat C; 1867 1868 ierr = MatCreateSubMatrix(B,irows,icols,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 1869 ierr = MatHeaderReplace(mat,&C);CHKERRQ(ierr); 1870 } 1871 ierr = MatDestroy(&B);CHKERRQ(ierr); 1872 ierr = ISDestroy(&icols);CHKERRQ(ierr); 1873 ierr = ISDestroy(&irows);CHKERRQ(ierr); 1874 PetscFunctionReturn(0); 1875 } 1876 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 1877 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 1878 ierr = MatGetLocalSize(mat,&lrows,&lcols);CHKERRQ(ierr); 1879 ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr); 1880 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 1881 ierr = PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 1882 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);CHKERRQ(ierr); 1883 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);CHKERRQ(ierr); 1884 if (!isseqdense && !isseqaij && !isseqbaij && !isseqsbaij) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name); 1885 #if defined (PETSC_USE_DEBUG) 1886 lb[0] = isseqdense; 1887 lb[1] = isseqaij; 1888 lb[2] = isseqbaij; 1889 lb[3] = isseqsbaij; 1890 ierr = MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1891 if (!bb[0] && !bb[1] && !bb[2] && !bb[3]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type"); 1892 #endif 1893 1894 if (reuse != MAT_REUSE_MATRIX) { 1895 PetscBool isaij; 1896 1897 ierr = MatCreate(PetscObjectComm((PetscObject)mat),&MT);CHKERRQ(ierr); 1898 ierr = MatSetSizes(MT,lrows,lcols,rows,cols);CHKERRQ(ierr); 1899 ierr = PetscStrcmp(mtype,MATAIJ,&isaij);CHKERRQ(ierr); 1900 if (isaij) { 1901 if (!isseqsbaij) { 1902 ierr = MatSetType(MT,MATAIJ);CHKERRQ(ierr); 1903 } else { 1904 ierr = MatSetType(MT,MATSBAIJ);CHKERRQ(ierr); 1905 } 1906 } else { 1907 ierr = MatSetType(MT,mtype);CHKERRQ(ierr); 1908 } 1909 ierr = MatSetBlockSize(MT,bs);CHKERRQ(ierr); 1910 ierr = MatISSetMPIXAIJPreallocation_Private(mat,MT,PETSC_FALSE);CHKERRQ(ierr); 1911 } else { 1912 PetscInt mbs,mrows,mcols,mlrows,mlcols; 1913 1914 /* some checks */ 1915 MT = *M; 1916 ierr = MatGetBlockSize(MT,&mbs);CHKERRQ(ierr); 1917 ierr = MatGetSize(MT,&mrows,&mcols);CHKERRQ(ierr); 1918 ierr = MatGetLocalSize(MT,&mlrows,&mlcols);CHKERRQ(ierr); 1919 if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows); 1920 if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols); 1921 if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows); 1922 if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols); 1923 if (mbs != bs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs); 1924 ierr = MatZeroEntries(MT);CHKERRQ(ierr); 1925 } 1926 1927 if (isseqsbaij) { 1928 ierr = MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);CHKERRQ(ierr); 1929 } else { 1930 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1931 local_mat = matis->A; 1932 } 1933 1934 /* Set values */ 1935 ierr = MatSetLocalToGlobalMapping(MT,mat->rmap->mapping,mat->cmap->mapping);CHKERRQ(ierr); 1936 if (isseqdense) { /* special case for dense local matrices */ 1937 PetscInt i,*dummy; 1938 1939 ierr = PetscMalloc1(PetscMax(local_rows,local_cols),&dummy);CHKERRQ(ierr); 1940 for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i; 1941 ierr = MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 1942 ierr = MatDenseGetArray(local_mat,&array);CHKERRQ(ierr); 1943 ierr = MatSetValuesLocal(MT,local_rows,dummy,local_cols,dummy,array,ADD_VALUES);CHKERRQ(ierr); 1944 ierr = MatDenseRestoreArray(local_mat,&array);CHKERRQ(ierr); 1945 ierr = PetscFree(dummy);CHKERRQ(ierr); 1946 } else if (isseqaij) { 1947 PetscInt i,nvtxs,*xadj,*adjncy; 1948 PetscBool done; 1949 1950 ierr = MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1951 if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ"); 1952 ierr = MatSeqAIJGetArray(local_mat,&array);CHKERRQ(ierr); 1953 for (i=0;i<nvtxs;i++) { 1954 ierr = MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);CHKERRQ(ierr); 1955 } 1956 ierr = MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1957 if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ"); 1958 ierr = MatSeqAIJRestoreArray(local_mat,&array);CHKERRQ(ierr); 1959 } else { /* very basic values insertion for all other matrix types */ 1960 PetscInt i; 1961 1962 for (i=0;i<local_rows;i++) { 1963 PetscInt j; 1964 const PetscInt *local_indices_cols; 1965 1966 ierr = MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1967 ierr = MatSetValuesLocal(MT,1,&i,j,local_indices_cols,array,ADD_VALUES);CHKERRQ(ierr); 1968 ierr = MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1969 } 1970 } 1971 ierr = MatAssemblyBegin(MT,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1972 ierr = MatDestroy(&local_mat);CHKERRQ(ierr); 1973 ierr = MatAssemblyEnd(MT,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1974 if (isseqdense) { 1975 ierr = MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 1976 } 1977 if (reuse == MAT_INPLACE_MATRIX) { 1978 ierr = MatHeaderReplace(mat,&MT);CHKERRQ(ierr); 1979 } else if (reuse == MAT_INITIAL_MATRIX) { 1980 *M = MT; 1981 } 1982 PetscFunctionReturn(0); 1983 } 1984 1985 /*@ 1986 MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format 1987 1988 Input Parameter: 1989 . mat - the matrix (should be of type MATIS) 1990 . reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1991 1992 Output Parameter: 1993 . newmat - the matrix in AIJ format 1994 1995 Level: developer 1996 1997 Notes: 1998 This function has been deprecated and it will be removed in future releases. Update your code to use the MatConvert() interface. 1999 2000 .seealso: MATIS, MatConvert() 2001 @*/ 2002 PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat) 2003 { 2004 PetscErrorCode ierr; 2005 2006 PetscFunctionBegin; 2007 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2008 PetscValidLogicalCollectiveEnum(mat,reuse,2); 2009 PetscValidPointer(newmat,3); 2010 if (reuse == MAT_REUSE_MATRIX) { 2011 PetscValidHeaderSpecific(*newmat,MAT_CLASSID,3); 2012 PetscCheckSameComm(mat,1,*newmat,3); 2013 if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix"); 2014 } 2015 ierr = PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatType,MatReuse,Mat*),(mat,MATAIJ,reuse,newmat));CHKERRQ(ierr); 2016 PetscFunctionReturn(0); 2017 } 2018 2019 PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat) 2020 { 2021 PetscErrorCode ierr; 2022 Mat_IS *matis = (Mat_IS*)(mat->data); 2023 PetscInt rbs,cbs,m,n,M,N; 2024 Mat B,localmat; 2025 2026 PetscFunctionBegin; 2027 ierr = ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);CHKERRQ(ierr); 2028 ierr = ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);CHKERRQ(ierr); 2029 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 2030 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 2031 ierr = MatCreateIS(PetscObjectComm((PetscObject)mat),rbs == cbs ? rbs : 1,m,n,M,N,mat->rmap->mapping,mat->cmap->mapping,&B);CHKERRQ(ierr); 2032 ierr = MatDuplicate(matis->A,op,&localmat);CHKERRQ(ierr); 2033 ierr = MatISSetLocalMat(B,localmat);CHKERRQ(ierr); 2034 ierr = MatDestroy(&localmat);CHKERRQ(ierr); 2035 if (matis->sf) { 2036 Mat_IS *bmatis = (Mat_IS*)(B->data); 2037 2038 ierr = PetscObjectReference((PetscObject)matis->sf);CHKERRQ(ierr); 2039 bmatis->sf = matis->sf; 2040 ierr = PetscMalloc2(matis->sf->nroots,&bmatis->sf_rootdata,matis->sf->nleaves,&bmatis->sf_leafdata);CHKERRQ(ierr); 2041 if (matis->sf != matis->csf) { 2042 ierr = PetscObjectReference((PetscObject)matis->csf);CHKERRQ(ierr); 2043 bmatis->csf = matis->csf; 2044 ierr = PetscMalloc2(matis->csf->nroots,&bmatis->csf_rootdata,matis->csf->nleaves,&bmatis->csf_leafdata);CHKERRQ(ierr); 2045 } else { 2046 bmatis->csf = bmatis->sf; 2047 bmatis->csf_leafdata = bmatis->sf_leafdata; 2048 bmatis->csf_rootdata = bmatis->sf_rootdata; 2049 } 2050 } 2051 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2052 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2053 *newmat = B; 2054 PetscFunctionReturn(0); 2055 } 2056 2057 static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg) 2058 { 2059 PetscErrorCode ierr; 2060 Mat_IS *matis = (Mat_IS*)A->data; 2061 PetscBool local_sym; 2062 2063 PetscFunctionBegin; 2064 ierr = MatIsHermitian(matis->A,tol,&local_sym);CHKERRQ(ierr); 2065 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2066 PetscFunctionReturn(0); 2067 } 2068 2069 static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg) 2070 { 2071 PetscErrorCode ierr; 2072 Mat_IS *matis = (Mat_IS*)A->data; 2073 PetscBool local_sym; 2074 2075 PetscFunctionBegin; 2076 ierr = MatIsSymmetric(matis->A,tol,&local_sym);CHKERRQ(ierr); 2077 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2078 PetscFunctionReturn(0); 2079 } 2080 2081 static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool *flg) 2082 { 2083 PetscErrorCode ierr; 2084 Mat_IS *matis = (Mat_IS*)A->data; 2085 PetscBool local_sym; 2086 2087 PetscFunctionBegin; 2088 if (A->rmap->mapping != A->cmap->mapping) { 2089 *flg = PETSC_FALSE; 2090 PetscFunctionReturn(0); 2091 } 2092 ierr = MatIsStructurallySymmetric(matis->A,&local_sym);CHKERRQ(ierr); 2093 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2094 PetscFunctionReturn(0); 2095 } 2096 2097 static PetscErrorCode MatDestroy_IS(Mat A) 2098 { 2099 PetscErrorCode ierr; 2100 Mat_IS *b = (Mat_IS*)A->data; 2101 2102 PetscFunctionBegin; 2103 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 2104 ierr = VecScatterDestroy(&b->cctx);CHKERRQ(ierr); 2105 ierr = VecScatterDestroy(&b->rctx);CHKERRQ(ierr); 2106 ierr = VecDestroy(&b->x);CHKERRQ(ierr); 2107 ierr = VecDestroy(&b->y);CHKERRQ(ierr); 2108 ierr = VecDestroy(&b->counter);CHKERRQ(ierr); 2109 ierr = ISDestroy(&b->getsub_ris);CHKERRQ(ierr); 2110 ierr = ISDestroy(&b->getsub_cis);CHKERRQ(ierr); 2111 if (b->sf != b->csf) { 2112 ierr = PetscSFDestroy(&b->csf);CHKERRQ(ierr); 2113 ierr = PetscFree2(b->csf_rootdata,b->csf_leafdata);CHKERRQ(ierr); 2114 } else b->csf = NULL; 2115 ierr = PetscSFDestroy(&b->sf);CHKERRQ(ierr); 2116 ierr = PetscFree2(b->sf_rootdata,b->sf_leafdata);CHKERRQ(ierr); 2117 ierr = PetscFree(A->data);CHKERRQ(ierr); 2118 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 2119 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);CHKERRQ(ierr); 2120 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);CHKERRQ(ierr); 2121 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);CHKERRQ(ierr); 2122 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);CHKERRQ(ierr); 2123 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",NULL);CHKERRQ(ierr); 2124 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL);CHKERRQ(ierr); 2125 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",NULL);CHKERRQ(ierr); 2126 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL);CHKERRQ(ierr); 2127 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL);CHKERRQ(ierr); 2128 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL);CHKERRQ(ierr); 2129 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL);CHKERRQ(ierr); 2130 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL);CHKERRQ(ierr); 2131 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL);CHKERRQ(ierr); 2132 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL);CHKERRQ(ierr); 2133 PetscFunctionReturn(0); 2134 } 2135 2136 static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y) 2137 { 2138 PetscErrorCode ierr; 2139 Mat_IS *is = (Mat_IS*)A->data; 2140 PetscScalar zero = 0.0; 2141 2142 PetscFunctionBegin; 2143 /* scatter the global vector x into the local work vector */ 2144 ierr = VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2145 ierr = VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2146 2147 /* multiply the local matrix */ 2148 ierr = MatMult(is->A,is->x,is->y);CHKERRQ(ierr); 2149 2150 /* scatter product back into global memory */ 2151 ierr = VecSet(y,zero);CHKERRQ(ierr); 2152 ierr = VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2153 ierr = VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2154 PetscFunctionReturn(0); 2155 } 2156 2157 static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 2158 { 2159 Vec temp_vec; 2160 PetscErrorCode ierr; 2161 2162 PetscFunctionBegin; /* v3 = v2 + A * v1.*/ 2163 if (v3 != v2) { 2164 ierr = MatMult(A,v1,v3);CHKERRQ(ierr); 2165 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 2166 } else { 2167 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 2168 ierr = MatMult(A,v1,temp_vec);CHKERRQ(ierr); 2169 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 2170 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 2171 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 2172 } 2173 PetscFunctionReturn(0); 2174 } 2175 2176 static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x) 2177 { 2178 Mat_IS *is = (Mat_IS*)A->data; 2179 PetscErrorCode ierr; 2180 2181 PetscFunctionBegin; 2182 /* scatter the global vector x into the local work vector */ 2183 ierr = VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2184 ierr = VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2185 2186 /* multiply the local matrix */ 2187 ierr = MatMultTranspose(is->A,is->y,is->x);CHKERRQ(ierr); 2188 2189 /* scatter product back into global vector */ 2190 ierr = VecSet(x,0);CHKERRQ(ierr); 2191 ierr = VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2192 ierr = VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2193 PetscFunctionReturn(0); 2194 } 2195 2196 static PetscErrorCode MatMultTransposeAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 2197 { 2198 Vec temp_vec; 2199 PetscErrorCode ierr; 2200 2201 PetscFunctionBegin; /* v3 = v2 + A' * v1.*/ 2202 if (v3 != v2) { 2203 ierr = MatMultTranspose(A,v1,v3);CHKERRQ(ierr); 2204 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 2205 } else { 2206 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 2207 ierr = MatMultTranspose(A,v1,temp_vec);CHKERRQ(ierr); 2208 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 2209 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 2210 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 2211 } 2212 PetscFunctionReturn(0); 2213 } 2214 2215 static PetscErrorCode MatView_IS(Mat A,PetscViewer viewer) 2216 { 2217 Mat_IS *a = (Mat_IS*)A->data; 2218 PetscErrorCode ierr; 2219 PetscViewer sviewer; 2220 PetscBool isascii,view = PETSC_TRUE; 2221 2222 PetscFunctionBegin; 2223 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 2224 if (isascii) { 2225 PetscViewerFormat format; 2226 2227 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 2228 if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE; 2229 } 2230 if (!view) PetscFunctionReturn(0); 2231 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 2232 ierr = MatView(a->A,sviewer);CHKERRQ(ierr); 2233 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 2234 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2235 PetscFunctionReturn(0); 2236 } 2237 2238 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping) 2239 { 2240 PetscErrorCode ierr; 2241 PetscInt nr,rbs,nc,cbs; 2242 Mat_IS *is = (Mat_IS*)A->data; 2243 IS from,to; 2244 Vec cglobal,rglobal; 2245 2246 PetscFunctionBegin; 2247 PetscCheckSameComm(A,1,rmapping,2); 2248 PetscCheckSameComm(A,1,cmapping,3); 2249 /* Destroy any previous data */ 2250 ierr = VecDestroy(&is->x);CHKERRQ(ierr); 2251 ierr = VecDestroy(&is->y);CHKERRQ(ierr); 2252 ierr = VecDestroy(&is->counter);CHKERRQ(ierr); 2253 ierr = VecScatterDestroy(&is->rctx);CHKERRQ(ierr); 2254 ierr = VecScatterDestroy(&is->cctx);CHKERRQ(ierr); 2255 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 2256 if (is->csf != is->sf) { 2257 ierr = PetscSFDestroy(&is->csf);CHKERRQ(ierr); 2258 ierr = PetscFree2(is->csf_rootdata,is->csf_leafdata);CHKERRQ(ierr); 2259 } else is->csf = NULL; 2260 ierr = PetscSFDestroy(&is->sf);CHKERRQ(ierr); 2261 ierr = PetscFree2(is->sf_rootdata,is->sf_leafdata);CHKERRQ(ierr); 2262 2263 /* Setup Layout and set local to global maps */ 2264 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 2265 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 2266 ierr = ISLocalToGlobalMappingGetSize(rmapping,&nr);CHKERRQ(ierr); 2267 ierr = ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);CHKERRQ(ierr); 2268 ierr = ISLocalToGlobalMappingGetSize(cmapping,&nc);CHKERRQ(ierr); 2269 ierr = ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);CHKERRQ(ierr); 2270 /* check if the two mappings are actually the same for square matrices (DOLFIN passes 2 different objects) */ 2271 if (rmapping != cmapping && A->rmap->N == A->cmap->N) { 2272 PetscBool same,gsame; 2273 2274 same = PETSC_FALSE; 2275 if (nr == nc && cbs == rbs) { 2276 const PetscInt *idxs1,*idxs2; 2277 2278 ierr = ISLocalToGlobalMappingGetBlockIndices(rmapping,&idxs1);CHKERRQ(ierr); 2279 ierr = ISLocalToGlobalMappingGetBlockIndices(cmapping,&idxs2);CHKERRQ(ierr); 2280 ierr = PetscMemcmp(idxs1,idxs2,(nr/rbs)*sizeof(PetscInt),&same);CHKERRQ(ierr); 2281 ierr = ISLocalToGlobalMappingRestoreBlockIndices(rmapping,&idxs1);CHKERRQ(ierr); 2282 ierr = ISLocalToGlobalMappingRestoreBlockIndices(cmapping,&idxs2);CHKERRQ(ierr); 2283 } 2284 ierr = MPIU_Allreduce(&same,&gsame,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2285 if (gsame) cmapping = rmapping; 2286 } 2287 ierr = PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);CHKERRQ(ierr); 2288 ierr = PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);CHKERRQ(ierr); 2289 2290 /* Create the local matrix A */ 2291 ierr = MatCreate(PETSC_COMM_SELF,&is->A);CHKERRQ(ierr); 2292 ierr = MatSetType(is->A,MATAIJ);CHKERRQ(ierr); 2293 ierr = MatSetSizes(is->A,nr,nc,nr,nc);CHKERRQ(ierr); 2294 ierr = MatSetBlockSizes(is->A,rbs,cbs);CHKERRQ(ierr); 2295 ierr = MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);CHKERRQ(ierr); 2296 ierr = MatAppendOptionsPrefix(is->A,"is_");CHKERRQ(ierr); 2297 ierr = MatSetFromOptions(is->A);CHKERRQ(ierr); 2298 ierr = PetscLayoutSetUp(is->A->rmap);CHKERRQ(ierr); 2299 ierr = PetscLayoutSetUp(is->A->cmap);CHKERRQ(ierr); 2300 2301 if (!is->islocalref) { /* setup scatters and local vectors for MatMult */ 2302 /* Create the local work vectors */ 2303 ierr = MatCreateVecs(is->A,&is->x,&is->y);CHKERRQ(ierr); 2304 2305 /* setup the global to local scatters */ 2306 ierr = MatCreateVecs(A,&cglobal,&rglobal);CHKERRQ(ierr); 2307 ierr = ISCreateStride(PETSC_COMM_SELF,nr,0,1,&to);CHKERRQ(ierr); 2308 ierr = ISLocalToGlobalMappingApplyIS(rmapping,to,&from);CHKERRQ(ierr); 2309 ierr = VecScatterCreate(rglobal,from,is->y,to,&is->rctx);CHKERRQ(ierr); 2310 if (rmapping != cmapping) { 2311 ierr = ISDestroy(&to);CHKERRQ(ierr); 2312 ierr = ISDestroy(&from);CHKERRQ(ierr); 2313 ierr = ISCreateStride(PETSC_COMM_SELF,nc,0,1,&to);CHKERRQ(ierr); 2314 ierr = ISLocalToGlobalMappingApplyIS(cmapping,to,&from);CHKERRQ(ierr); 2315 ierr = VecScatterCreate(cglobal,from,is->x,to,&is->cctx);CHKERRQ(ierr); 2316 } else { 2317 ierr = PetscObjectReference((PetscObject)is->rctx);CHKERRQ(ierr); 2318 is->cctx = is->rctx; 2319 } 2320 2321 /* interface counter vector (local) */ 2322 ierr = VecDuplicate(is->y,&is->counter);CHKERRQ(ierr); 2323 ierr = VecSet(is->y,1.);CHKERRQ(ierr); 2324 ierr = VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2325 ierr = VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2326 ierr = VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2327 ierr = VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2328 2329 /* free workspace */ 2330 ierr = VecDestroy(&rglobal);CHKERRQ(ierr); 2331 ierr = VecDestroy(&cglobal);CHKERRQ(ierr); 2332 ierr = ISDestroy(&to);CHKERRQ(ierr); 2333 ierr = ISDestroy(&from);CHKERRQ(ierr); 2334 } 2335 ierr = MatSetUp(A);CHKERRQ(ierr); 2336 PetscFunctionReturn(0); 2337 } 2338 2339 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2340 { 2341 Mat_IS *is = (Mat_IS*)mat->data; 2342 PetscErrorCode ierr; 2343 #if defined(PETSC_USE_DEBUG) 2344 PetscInt i,zm,zn; 2345 #endif 2346 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2347 2348 PetscFunctionBegin; 2349 #if defined(PETSC_USE_DEBUG) 2350 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); 2351 /* count negative indices */ 2352 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 2353 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 2354 #endif 2355 ierr = ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 2356 ierr = ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 2357 #if defined(PETSC_USE_DEBUG) 2358 /* count negative indices (should be the same as before) */ 2359 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 2360 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 2361 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"); 2362 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"); 2363 #endif 2364 ierr = MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 2365 PetscFunctionReturn(0); 2366 } 2367 2368 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2369 { 2370 Mat_IS *is = (Mat_IS*)mat->data; 2371 PetscErrorCode ierr; 2372 #if defined(PETSC_USE_DEBUG) 2373 PetscInt i,zm,zn; 2374 #endif 2375 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2376 2377 PetscFunctionBegin; 2378 #if defined(PETSC_USE_DEBUG) 2379 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); 2380 /* count negative indices */ 2381 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 2382 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 2383 #endif 2384 ierr = ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 2385 ierr = ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 2386 #if defined(PETSC_USE_DEBUG) 2387 /* count negative indices (should be the same as before) */ 2388 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 2389 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 2390 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"); 2391 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"); 2392 #endif 2393 ierr = MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 2394 PetscFunctionReturn(0); 2395 } 2396 2397 static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2398 { 2399 PetscErrorCode ierr; 2400 Mat_IS *is = (Mat_IS*)A->data; 2401 2402 PetscFunctionBegin; 2403 if (is->A->rmap->mapping) { 2404 ierr = MatSetValuesLocal(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2405 } else { 2406 ierr = MatSetValues(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2407 } 2408 PetscFunctionReturn(0); 2409 } 2410 2411 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2412 { 2413 PetscErrorCode ierr; 2414 Mat_IS *is = (Mat_IS*)A->data; 2415 2416 PetscFunctionBegin; 2417 if (is->A->rmap->mapping) { 2418 #if defined(PETSC_USE_DEBUG) 2419 PetscInt ibs,bs; 2420 2421 ierr = ISLocalToGlobalMappingGetBlockSize(is->A->rmap->mapping,&ibs);CHKERRQ(ierr); 2422 ierr = MatGetBlockSize(is->A,&bs);CHKERRQ(ierr); 2423 if (ibs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Different block sizes! local mat %D, local l2g map %D",bs,ibs); 2424 #endif 2425 ierr = MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2426 } else { 2427 ierr = MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2428 } 2429 PetscFunctionReturn(0); 2430 } 2431 2432 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns) 2433 { 2434 Mat_IS *is = (Mat_IS*)A->data; 2435 PetscErrorCode ierr; 2436 2437 PetscFunctionBegin; 2438 if (!n) { 2439 is->pure_neumann = PETSC_TRUE; 2440 } else { 2441 PetscInt i; 2442 is->pure_neumann = PETSC_FALSE; 2443 2444 if (columns) { 2445 ierr = MatZeroRowsColumns(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 2446 } else { 2447 ierr = MatZeroRows(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 2448 } 2449 if (diag != 0.) { 2450 const PetscScalar *array; 2451 ierr = VecGetArrayRead(is->counter,&array);CHKERRQ(ierr); 2452 for (i=0; i<n; i++) { 2453 ierr = MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);CHKERRQ(ierr); 2454 } 2455 ierr = VecRestoreArrayRead(is->counter,&array);CHKERRQ(ierr); 2456 } 2457 ierr = MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2458 ierr = MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2459 } 2460 PetscFunctionReturn(0); 2461 } 2462 2463 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns) 2464 { 2465 Mat_IS *matis = (Mat_IS*)A->data; 2466 PetscInt nr,nl,len,i; 2467 PetscInt *lrows; 2468 PetscErrorCode ierr; 2469 2470 PetscFunctionBegin; 2471 #if defined(PETSC_USE_DEBUG) 2472 if (columns || diag != 0. || (x && b)) { 2473 PetscBool cong; 2474 2475 ierr = PetscLayoutCompare(A->rmap,A->cmap,&cong);CHKERRQ(ierr); 2476 cong = (PetscBool)(cong && matis->sf == matis->csf); 2477 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"); 2478 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"); 2479 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"); 2480 } 2481 #endif 2482 /* get locally owned rows */ 2483 ierr = PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);CHKERRQ(ierr); 2484 /* fix right hand side if needed */ 2485 if (x && b) { 2486 const PetscScalar *xx; 2487 PetscScalar *bb; 2488 2489 ierr = VecGetArrayRead(x, &xx);CHKERRQ(ierr); 2490 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 2491 for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]]; 2492 ierr = VecRestoreArrayRead(x, &xx);CHKERRQ(ierr); 2493 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 2494 } 2495 /* get rows associated to the local matrices */ 2496 ierr = MatISSetUpSF(A);CHKERRQ(ierr); 2497 ierr = MatGetSize(matis->A,&nl,NULL);CHKERRQ(ierr); 2498 ierr = PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 2499 ierr = PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2500 for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1; 2501 ierr = PetscFree(lrows);CHKERRQ(ierr); 2502 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 2503 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 2504 ierr = PetscMalloc1(nl,&lrows);CHKERRQ(ierr); 2505 for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i; 2506 ierr = MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);CHKERRQ(ierr); 2507 ierr = PetscFree(lrows);CHKERRQ(ierr); 2508 PetscFunctionReturn(0); 2509 } 2510 2511 static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2512 { 2513 PetscErrorCode ierr; 2514 2515 PetscFunctionBegin; 2516 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);CHKERRQ(ierr); 2517 PetscFunctionReturn(0); 2518 } 2519 2520 static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2521 { 2522 PetscErrorCode ierr; 2523 2524 PetscFunctionBegin; 2525 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);CHKERRQ(ierr); 2526 PetscFunctionReturn(0); 2527 } 2528 2529 static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type) 2530 { 2531 Mat_IS *is = (Mat_IS*)A->data; 2532 PetscErrorCode ierr; 2533 2534 PetscFunctionBegin; 2535 ierr = MatAssemblyBegin(is->A,type);CHKERRQ(ierr); 2536 PetscFunctionReturn(0); 2537 } 2538 2539 static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type) 2540 { 2541 Mat_IS *is = (Mat_IS*)A->data; 2542 PetscErrorCode ierr; 2543 2544 PetscFunctionBegin; 2545 ierr = MatAssemblyEnd(is->A,type);CHKERRQ(ierr); 2546 /* fix for local empty rows/cols */ 2547 if (is->locempty && type == MAT_FINAL_ASSEMBLY) { 2548 Mat newlA; 2549 ISLocalToGlobalMapping rl2g,cl2g; 2550 IS nzr,nzc; 2551 PetscInt nr,nc,nnzr,nnzc; 2552 PetscBool lnewl2g,newl2g; 2553 2554 ierr = MatGetSize(is->A,&nr,&nc);CHKERRQ(ierr); 2555 ierr = MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_FALSE,PETSC_SMALL,&nzr);CHKERRQ(ierr); 2556 if (!nzr) { 2557 ierr = ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&nzr);CHKERRQ(ierr); 2558 } 2559 ierr = MatFindNonzeroRowsOrCols_Basic(is->A,PETSC_TRUE,PETSC_SMALL,&nzc);CHKERRQ(ierr); 2560 if (!nzc) { 2561 ierr = ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&nzc);CHKERRQ(ierr); 2562 } 2563 ierr = ISGetSize(nzr,&nnzr);CHKERRQ(ierr); 2564 ierr = ISGetSize(nzc,&nnzc);CHKERRQ(ierr); 2565 if (nnzr != nr || nnzc != nc) { 2566 ISLocalToGlobalMapping l2g; 2567 IS is1,is2; 2568 2569 /* need new global l2g map */ 2570 lnewl2g = PETSC_TRUE; 2571 ierr = MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2572 2573 /* extract valid submatrix */ 2574 ierr = MatCreateSubMatrix(is->A,nzr,nzc,MAT_INITIAL_MATRIX,&newlA);CHKERRQ(ierr); 2575 2576 /* attach local l2g maps for successive calls of MatSetValues on the local matrix */ 2577 ierr = ISLocalToGlobalMappingCreateIS(nzr,&l2g);CHKERRQ(ierr); 2578 ierr = ISCreateStride(PetscObjectComm((PetscObject)is->A),nr,0,1,&is1);CHKERRQ(ierr); 2579 ierr = ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);CHKERRQ(ierr); 2580 ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr); 2581 ierr = ISLocalToGlobalMappingCreateIS(is2,&rl2g);CHKERRQ(ierr); 2582 ierr = ISDestroy(&is1);CHKERRQ(ierr); 2583 ierr = ISDestroy(&is2);CHKERRQ(ierr); 2584 ierr = ISLocalToGlobalMappingCreateIS(nzc,&l2g);CHKERRQ(ierr); 2585 ierr = ISCreateStride(PetscObjectComm((PetscObject)is->A),nc,0,1,&is1);CHKERRQ(ierr); 2586 ierr = ISGlobalToLocalMappingApplyIS(l2g,IS_GTOLM_MASK,is1,&is2);CHKERRQ(ierr); 2587 ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr); 2588 ierr = ISLocalToGlobalMappingCreateIS(is2,&cl2g);CHKERRQ(ierr); 2589 ierr = ISDestroy(&is1);CHKERRQ(ierr); 2590 ierr = ISDestroy(&is2);CHKERRQ(ierr); 2591 ierr = MatSetLocalToGlobalMapping(newlA,rl2g,cl2g);CHKERRQ(ierr); 2592 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 2593 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 2594 } else { /* local matrix fully populated */ 2595 lnewl2g = PETSC_FALSE; 2596 ierr = MPI_Allreduce(&lnewl2g,&newl2g,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2597 ierr = PetscObjectReference((PetscObject)is->A);CHKERRQ(ierr); 2598 newlA = is->A; 2599 } 2600 /* attach new global l2g map if needed */ 2601 if (newl2g) { 2602 IS gnzr,gnzc; 2603 const PetscInt *grid,*gcid; 2604 2605 ierr = ISLocalToGlobalMappingApplyIS(A->rmap->mapping,nzr,&gnzr);CHKERRQ(ierr); 2606 ierr = ISLocalToGlobalMappingApplyIS(A->cmap->mapping,nzc,&gnzc);CHKERRQ(ierr); 2607 ierr = ISGetIndices(gnzr,&grid);CHKERRQ(ierr); 2608 ierr = ISGetIndices(gnzc,&gcid);CHKERRQ(ierr); 2609 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzr,grid,PETSC_COPY_VALUES,&rl2g);CHKERRQ(ierr); 2610 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,nnzc,gcid,PETSC_COPY_VALUES,&cl2g);CHKERRQ(ierr); 2611 ierr = ISRestoreIndices(gnzr,&grid);CHKERRQ(ierr); 2612 ierr = ISRestoreIndices(gnzc,&gcid);CHKERRQ(ierr); 2613 ierr = ISDestroy(&gnzr);CHKERRQ(ierr); 2614 ierr = ISDestroy(&gnzc);CHKERRQ(ierr); 2615 ierr = MatSetLocalToGlobalMapping(A,rl2g,cl2g);CHKERRQ(ierr); 2616 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 2617 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 2618 } 2619 ierr = MatISSetLocalMat(A,newlA);CHKERRQ(ierr); 2620 ierr = MatDestroy(&newlA);CHKERRQ(ierr); 2621 ierr = ISDestroy(&nzr);CHKERRQ(ierr); 2622 ierr = ISDestroy(&nzc);CHKERRQ(ierr); 2623 is->locempty = PETSC_FALSE; 2624 } 2625 PetscFunctionReturn(0); 2626 } 2627 2628 static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local) 2629 { 2630 Mat_IS *is = (Mat_IS*)mat->data; 2631 2632 PetscFunctionBegin; 2633 *local = is->A; 2634 PetscFunctionReturn(0); 2635 } 2636 2637 static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local) 2638 { 2639 PetscFunctionBegin; 2640 *local = NULL; 2641 PetscFunctionReturn(0); 2642 } 2643 2644 /*@ 2645 MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix. 2646 2647 Input Parameter: 2648 . mat - the matrix 2649 2650 Output Parameter: 2651 . local - the local matrix 2652 2653 Level: advanced 2654 2655 Notes: 2656 This can be called if you have precomputed the nonzero structure of the 2657 matrix and want to provide it to the inner matrix object to improve the performance 2658 of the MatSetValues() operation. 2659 2660 Call MatISRestoreLocalMat() when finished with the local matrix. 2661 2662 .seealso: MATIS 2663 @*/ 2664 PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local) 2665 { 2666 PetscErrorCode ierr; 2667 2668 PetscFunctionBegin; 2669 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2670 PetscValidPointer(local,2); 2671 ierr = PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 2672 PetscFunctionReturn(0); 2673 } 2674 2675 /*@ 2676 MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat() 2677 2678 Input Parameter: 2679 . mat - the matrix 2680 2681 Output Parameter: 2682 . local - the local matrix 2683 2684 Level: advanced 2685 2686 .seealso: MATIS 2687 @*/ 2688 PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local) 2689 { 2690 PetscErrorCode ierr; 2691 2692 PetscFunctionBegin; 2693 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2694 PetscValidPointer(local,2); 2695 ierr = PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 2696 PetscFunctionReturn(0); 2697 } 2698 2699 static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local) 2700 { 2701 Mat_IS *is = (Mat_IS*)mat->data; 2702 PetscInt nrows,ncols,orows,ocols; 2703 PetscErrorCode ierr; 2704 2705 PetscFunctionBegin; 2706 if (is->A) { 2707 ierr = MatGetSize(is->A,&orows,&ocols);CHKERRQ(ierr); 2708 ierr = MatGetSize(local,&nrows,&ncols);CHKERRQ(ierr); 2709 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); 2710 } 2711 ierr = PetscObjectReference((PetscObject)local);CHKERRQ(ierr); 2712 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 2713 is->A = local; 2714 PetscFunctionReturn(0); 2715 } 2716 2717 /*@ 2718 MatISSetLocalMat - Replace the local matrix stored inside a MATIS object. 2719 2720 Input Parameter: 2721 . mat - the matrix 2722 . local - the local matrix 2723 2724 Output Parameter: 2725 2726 Level: advanced 2727 2728 Notes: 2729 This can be called if you have precomputed the local matrix and 2730 want to provide it to the matrix object MATIS. 2731 2732 .seealso: MATIS 2733 @*/ 2734 PetscErrorCode MatISSetLocalMat(Mat mat,Mat local) 2735 { 2736 PetscErrorCode ierr; 2737 2738 PetscFunctionBegin; 2739 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2740 PetscValidHeaderSpecific(local,MAT_CLASSID,2); 2741 ierr = PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));CHKERRQ(ierr); 2742 PetscFunctionReturn(0); 2743 } 2744 2745 static PetscErrorCode MatZeroEntries_IS(Mat A) 2746 { 2747 Mat_IS *a = (Mat_IS*)A->data; 2748 PetscErrorCode ierr; 2749 2750 PetscFunctionBegin; 2751 ierr = MatZeroEntries(a->A);CHKERRQ(ierr); 2752 PetscFunctionReturn(0); 2753 } 2754 2755 static PetscErrorCode MatScale_IS(Mat A,PetscScalar a) 2756 { 2757 Mat_IS *is = (Mat_IS*)A->data; 2758 PetscErrorCode ierr; 2759 2760 PetscFunctionBegin; 2761 ierr = MatScale(is->A,a);CHKERRQ(ierr); 2762 PetscFunctionReturn(0); 2763 } 2764 2765 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v) 2766 { 2767 Mat_IS *is = (Mat_IS*)A->data; 2768 PetscErrorCode ierr; 2769 2770 PetscFunctionBegin; 2771 /* get diagonal of the local matrix */ 2772 ierr = MatGetDiagonal(is->A,is->y);CHKERRQ(ierr); 2773 2774 /* scatter diagonal back into global vector */ 2775 ierr = VecSet(v,0);CHKERRQ(ierr); 2776 ierr = VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2777 ierr = VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2778 PetscFunctionReturn(0); 2779 } 2780 2781 static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg) 2782 { 2783 Mat_IS *a = (Mat_IS*)A->data; 2784 PetscErrorCode ierr; 2785 2786 PetscFunctionBegin; 2787 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 2788 PetscFunctionReturn(0); 2789 } 2790 2791 static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str) 2792 { 2793 Mat_IS *y = (Mat_IS*)Y->data; 2794 Mat_IS *x; 2795 #if defined(PETSC_USE_DEBUG) 2796 PetscBool ismatis; 2797 #endif 2798 PetscErrorCode ierr; 2799 2800 PetscFunctionBegin; 2801 #if defined(PETSC_USE_DEBUG) 2802 ierr = PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);CHKERRQ(ierr); 2803 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS"); 2804 #endif 2805 x = (Mat_IS*)X->data; 2806 ierr = MatAXPY(y->A,a,x->A,str);CHKERRQ(ierr); 2807 PetscFunctionReturn(0); 2808 } 2809 2810 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat) 2811 { 2812 Mat lA; 2813 Mat_IS *matis; 2814 ISLocalToGlobalMapping rl2g,cl2g; 2815 IS is; 2816 const PetscInt *rg,*rl; 2817 PetscInt nrg; 2818 PetscInt N,M,nrl,i,*idxs; 2819 PetscErrorCode ierr; 2820 2821 PetscFunctionBegin; 2822 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 2823 ierr = ISGetLocalSize(row,&nrl);CHKERRQ(ierr); 2824 ierr = ISGetIndices(row,&rl);CHKERRQ(ierr); 2825 ierr = ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);CHKERRQ(ierr); 2826 #if defined(PETSC_USE_DEBUG) 2827 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); 2828 #endif 2829 ierr = PetscMalloc1(nrg,&idxs);CHKERRQ(ierr); 2830 /* map from [0,nrl) to row */ 2831 for (i=0;i<nrl;i++) idxs[i] = rl[i]; 2832 #if defined(PETSC_USE_DEBUG) 2833 for (i=nrl;i<nrg;i++) idxs[i] = nrg; 2834 #else 2835 for (i=nrl;i<nrg;i++) idxs[i] = -1; 2836 #endif 2837 ierr = ISRestoreIndices(row,&rl);CHKERRQ(ierr); 2838 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 2839 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 2840 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 2841 ierr = ISDestroy(&is);CHKERRQ(ierr); 2842 /* compute new l2g map for columns */ 2843 if (col != row || A->rmap->mapping != A->cmap->mapping) { 2844 const PetscInt *cg,*cl; 2845 PetscInt ncg; 2846 PetscInt ncl; 2847 2848 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 2849 ierr = ISGetLocalSize(col,&ncl);CHKERRQ(ierr); 2850 ierr = ISGetIndices(col,&cl);CHKERRQ(ierr); 2851 ierr = ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);CHKERRQ(ierr); 2852 #if defined(PETSC_USE_DEBUG) 2853 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); 2854 #endif 2855 ierr = PetscMalloc1(ncg,&idxs);CHKERRQ(ierr); 2856 /* map from [0,ncl) to col */ 2857 for (i=0;i<ncl;i++) idxs[i] = cl[i]; 2858 #if defined(PETSC_USE_DEBUG) 2859 for (i=ncl;i<ncg;i++) idxs[i] = ncg; 2860 #else 2861 for (i=ncl;i<ncg;i++) idxs[i] = -1; 2862 #endif 2863 ierr = ISRestoreIndices(col,&cl);CHKERRQ(ierr); 2864 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 2865 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 2866 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 2867 ierr = ISDestroy(&is);CHKERRQ(ierr); 2868 } else { 2869 ierr = PetscObjectReference((PetscObject)rl2g);CHKERRQ(ierr); 2870 cl2g = rl2g; 2871 } 2872 /* create the MATIS submatrix */ 2873 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 2874 ierr = MatCreate(PetscObjectComm((PetscObject)A),submat);CHKERRQ(ierr); 2875 ierr = MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2876 ierr = MatSetType(*submat,MATIS);CHKERRQ(ierr); 2877 matis = (Mat_IS*)((*submat)->data); 2878 matis->islocalref = PETSC_TRUE; 2879 ierr = MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);CHKERRQ(ierr); 2880 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr); 2881 ierr = MatISSetLocalMat(*submat,lA);CHKERRQ(ierr); 2882 ierr = MatSetUp(*submat);CHKERRQ(ierr); 2883 ierr = MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2884 ierr = MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2885 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 2886 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 2887 /* remove unsupported ops */ 2888 ierr = PetscMemzero((*submat)->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2889 (*submat)->ops->destroy = MatDestroy_IS; 2890 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS; 2891 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS; 2892 (*submat)->ops->assemblybegin = MatAssemblyBegin_IS; 2893 (*submat)->ops->assemblyend = MatAssemblyEnd_IS; 2894 PetscFunctionReturn(0); 2895 } 2896 2897 static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A) 2898 { 2899 Mat_IS *a = (Mat_IS*)A->data; 2900 PetscErrorCode ierr; 2901 2902 PetscFunctionBegin; 2903 ierr = PetscOptionsHead(PetscOptionsObject,"MATIS options");CHKERRQ(ierr); 2904 ierr = PetscObjectOptionsBegin((PetscObject)A); 2905 ierr = PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns","MatISFixLocalEmpty",a->locempty,&a->locempty,NULL);CHKERRQ(ierr); 2906 ierr = PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL);CHKERRQ(ierr); 2907 ierr = PetscOptionsEnd();CHKERRQ(ierr); 2908 PetscFunctionReturn(0); 2909 } 2910 2911 /*@ 2912 MatCreateIS - Creates a "process" unassembled matrix, assembled on each 2913 process but not across processes. 2914 2915 Input Parameters: 2916 + comm - MPI communicator that will share the matrix 2917 . bs - block size of the matrix 2918 . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products 2919 . rmap - local to global map for rows 2920 - cmap - local to global map for cols 2921 2922 Output Parameter: 2923 . A - the resulting matrix 2924 2925 Level: advanced 2926 2927 Notes: 2928 See MATIS for more details. 2929 m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors 2930 used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices. 2931 If either rmap or cmap are NULL, then the matrix is assumed to be square. 2932 2933 .seealso: MATIS, MatSetLocalToGlobalMapping() 2934 @*/ 2935 PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A) 2936 { 2937 PetscErrorCode ierr; 2938 2939 PetscFunctionBegin; 2940 if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings"); 2941 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2942 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 2943 if (bs > 0) { 2944 ierr = MatSetBlockSize(*A,bs);CHKERRQ(ierr); 2945 } 2946 ierr = MatSetType(*A,MATIS);CHKERRQ(ierr); 2947 if (rmap && cmap) { 2948 ierr = MatSetLocalToGlobalMapping(*A,rmap,cmap);CHKERRQ(ierr); 2949 } else if (!rmap) { 2950 ierr = MatSetLocalToGlobalMapping(*A,cmap,cmap);CHKERRQ(ierr); 2951 } else { 2952 ierr = MatSetLocalToGlobalMapping(*A,rmap,rmap);CHKERRQ(ierr); 2953 } 2954 PetscFunctionReturn(0); 2955 } 2956 2957 /*MC 2958 MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP). 2959 This stores the matrices in globally unassembled form. Each processor 2960 assembles only its local Neumann problem and the parallel matrix vector 2961 product is handled "implicitly". 2962 2963 Operations Provided: 2964 + MatMult() 2965 . MatMultAdd() 2966 . MatMultTranspose() 2967 . MatMultTransposeAdd() 2968 . MatZeroEntries() 2969 . MatSetOption() 2970 . MatZeroRows() 2971 . MatSetValues() 2972 . MatSetValuesBlocked() 2973 . MatSetValuesLocal() 2974 . MatSetValuesBlockedLocal() 2975 . MatScale() 2976 . MatGetDiagonal() 2977 . MatMissingDiagonal() 2978 . MatDuplicate() 2979 . MatCopy() 2980 . MatAXPY() 2981 . MatCreateSubMatrix() 2982 . MatGetLocalSubMatrix() 2983 . MatTranspose() 2984 . MatPtAP() (with P of AIJ type) 2985 - MatSetLocalToGlobalMapping() 2986 2987 Options Database Keys: 2988 + -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions() 2989 . -matis_fixempty - Fixes local matrices in case of empty local rows/columns. 2990 - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of MatPtAP(). 2991 2992 Notes: 2993 Options prefix for the inner matrix are given by -is_mat_xxx 2994 2995 You must call MatSetLocalToGlobalMapping() before using this matrix type. 2996 2997 You can do matrix preallocation on the local matrix after you obtain it with 2998 MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation() 2999 3000 Level: advanced 3001 3002 .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP 3003 3004 M*/ 3005 3006 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A) 3007 { 3008 PetscErrorCode ierr; 3009 Mat_IS *b; 3010 3011 PetscFunctionBegin; 3012 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 3013 A->data = (void*)b; 3014 3015 /* matrix ops */ 3016 ierr = PetscMemzero(A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 3017 A->ops->mult = MatMult_IS; 3018 A->ops->multadd = MatMultAdd_IS; 3019 A->ops->multtranspose = MatMultTranspose_IS; 3020 A->ops->multtransposeadd = MatMultTransposeAdd_IS; 3021 A->ops->destroy = MatDestroy_IS; 3022 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS; 3023 A->ops->setvalues = MatSetValues_IS; 3024 A->ops->setvaluesblocked = MatSetValuesBlocked_IS; 3025 A->ops->setvalueslocal = MatSetValuesLocal_IS; 3026 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS; 3027 A->ops->zerorows = MatZeroRows_IS; 3028 A->ops->zerorowscolumns = MatZeroRowsColumns_IS; 3029 A->ops->assemblybegin = MatAssemblyBegin_IS; 3030 A->ops->assemblyend = MatAssemblyEnd_IS; 3031 A->ops->view = MatView_IS; 3032 A->ops->zeroentries = MatZeroEntries_IS; 3033 A->ops->scale = MatScale_IS; 3034 A->ops->getdiagonal = MatGetDiagonal_IS; 3035 A->ops->setoption = MatSetOption_IS; 3036 A->ops->ishermitian = MatIsHermitian_IS; 3037 A->ops->issymmetric = MatIsSymmetric_IS; 3038 A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS; 3039 A->ops->duplicate = MatDuplicate_IS; 3040 A->ops->missingdiagonal = MatMissingDiagonal_IS; 3041 A->ops->copy = MatCopy_IS; 3042 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 3043 A->ops->createsubmatrix = MatCreateSubMatrix_IS; 3044 A->ops->axpy = MatAXPY_IS; 3045 A->ops->diagonalset = MatDiagonalSet_IS; 3046 A->ops->shift = MatShift_IS; 3047 A->ops->transpose = MatTranspose_IS; 3048 A->ops->getinfo = MatGetInfo_IS; 3049 A->ops->diagonalscale = MatDiagonalScale_IS; 3050 A->ops->setfromoptions = MatSetFromOptions_IS; 3051 3052 /* special MATIS functions */ 3053 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);CHKERRQ(ierr); 3054 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);CHKERRQ(ierr); 3055 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);CHKERRQ(ierr); 3056 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3057 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);CHKERRQ(ierr); 3058 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",MatISSetUpSF_IS);CHKERRQ(ierr); 3059 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS);CHKERRQ(ierr); 3060 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISFixLocalEmpty_C",MatISFixLocalEmpty_IS);CHKERRQ(ierr); 3061 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3062 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3063 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3064 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3065 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3066 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3067 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 3068 ierr = PetscObjectChangeTypeName((PetscObject)A,MATIS);CHKERRQ(ierr); 3069 PetscFunctionReturn(0); 3070 } 3071