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 MatISSetPreallocation - Preallocates memory for a MATIS parallel matrix. 1536 1537 Collective on MPI_Comm 1538 1539 Input Parameters: 1540 + B - the matrix 1541 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1542 (same value is used for all local rows) 1543 . d_nnz - array containing the number of nonzeros in the various rows of the 1544 DIAGONAL portion of the local submatrix (possibly different for each row) 1545 or NULL, if d_nz is used to specify the nonzero structure. 1546 The size of this array is equal to the number of local rows, i.e 'm'. 1547 For matrices that will be factored, you must leave room for (and set) 1548 the diagonal entry even if it is zero. 1549 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1550 submatrix (same value is used for all local rows). 1551 - o_nnz - array containing the number of nonzeros in the various rows of the 1552 OFF-DIAGONAL portion of the local submatrix (possibly different for 1553 each row) or NULL, if o_nz is used to specify the nonzero 1554 structure. The size of this array is equal to the number 1555 of local rows, i.e 'm'. 1556 1557 If the *_nnz parameter is given then the *_nz parameter is ignored 1558 1559 Level: intermediate 1560 1561 Notes: 1562 This function has the same interface as the MPIAIJ preallocation routine in order to simplify the transition 1563 from the asssembled format to the unassembled one. It overestimates the preallocation of MATIS local 1564 matrices; for exact preallocation, the user should set the preallocation directly on local matrix objects. 1565 1566 .keywords: matrix 1567 1568 .seealso: MatCreate(), MatCreateIS(), MatMPIAIJSetPreallocation(), MatISGetLocalMat(), MATIS 1569 @*/ 1570 PetscErrorCode MatISSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1571 { 1572 PetscErrorCode ierr; 1573 1574 PetscFunctionBegin; 1575 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1576 PetscValidType(B,1); 1577 ierr = PetscTryMethod(B,"MatISSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz));CHKERRQ(ierr); 1578 PetscFunctionReturn(0); 1579 } 1580 1581 /* this is used by DMDA */ 1582 PETSC_EXTERN PetscErrorCode MatISSetPreallocation_IS(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1583 { 1584 Mat_IS *matis = (Mat_IS*)(B->data); 1585 PetscInt bs,i,nlocalcols; 1586 PetscErrorCode ierr; 1587 1588 PetscFunctionBegin; 1589 if (!matis->A) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"You should first call MatSetLocalToGlobalMapping"); 1590 ierr = MatISSetUpSF(B);CHKERRQ(ierr); 1591 1592 if (!d_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nz; 1593 else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] = d_nnz[i]; 1594 1595 if (!o_nnz) for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nz; 1596 else for (i=0;i<matis->sf->nroots;i++) matis->sf_rootdata[i] += o_nnz[i]; 1597 1598 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1599 ierr = MatGetSize(matis->A,NULL,&nlocalcols);CHKERRQ(ierr); 1600 ierr = MatGetBlockSize(matis->A,&bs);CHKERRQ(ierr); 1601 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 1602 1603 for (i=0;i<matis->sf->nleaves;i++) matis->sf_leafdata[i] = PetscMin(matis->sf_leafdata[i],nlocalcols); 1604 ierr = MatSeqAIJSetPreallocation(matis->A,0,matis->sf_leafdata);CHKERRQ(ierr); 1605 #if defined(PETSC_HAVE_HYPRE) 1606 ierr = MatHYPRESetPreallocation(matis->A,0,matis->sf_leafdata,0,NULL);CHKERRQ(ierr); 1607 #endif 1608 1609 for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i*bs]/bs; 1610 ierr = MatSeqBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);CHKERRQ(ierr); 1611 1612 for (i=0;i<matis->sf->nleaves/bs;i++) matis->sf_leafdata[i] = matis->sf_leafdata[i]-i; 1613 ierr = MatSeqSBAIJSetPreallocation(matis->A,bs,0,matis->sf_leafdata);CHKERRQ(ierr); 1614 1615 /* for other matrix types */ 1616 ierr = MatSetUp(matis->A);CHKERRQ(ierr); 1617 PetscFunctionReturn(0); 1618 } 1619 1620 PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce) 1621 { 1622 Mat_IS *matis = (Mat_IS*)(A->data); 1623 PetscInt *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership; 1624 const PetscInt *global_indices_r,*global_indices_c; 1625 PetscInt i,j,bs,rows,cols; 1626 PetscInt lrows,lcols; 1627 PetscInt local_rows,local_cols; 1628 PetscMPIInt nsubdomains; 1629 PetscBool isdense,issbaij; 1630 PetscErrorCode ierr; 1631 1632 PetscFunctionBegin; 1633 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&nsubdomains);CHKERRQ(ierr); 1634 ierr = MatGetSize(A,&rows,&cols);CHKERRQ(ierr); 1635 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 1636 ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr); 1637 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);CHKERRQ(ierr); 1638 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);CHKERRQ(ierr); 1639 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);CHKERRQ(ierr); 1640 if (A->rmap->mapping != A->cmap->mapping) { 1641 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&global_indices_c);CHKERRQ(ierr); 1642 } else { 1643 global_indices_c = global_indices_r; 1644 } 1645 1646 if (issbaij) { 1647 ierr = MatGetRowUpperTriangular(matis->A);CHKERRQ(ierr); 1648 } 1649 /* 1650 An SF reduce is needed to sum up properly on shared rows. 1651 Note that generally preallocation is not exact, since it overestimates nonzeros 1652 */ 1653 ierr = MatISSetUpSF(A);CHKERRQ(ierr); 1654 ierr = MatGetLocalSize(A,&lrows,&lcols);CHKERRQ(ierr); 1655 ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);CHKERRQ(ierr); 1656 /* All processes need to compute entire row ownership */ 1657 ierr = PetscMalloc1(rows,&row_ownership);CHKERRQ(ierr); 1658 ierr = MatGetOwnershipRanges(A,(const PetscInt**)&mat_ranges);CHKERRQ(ierr); 1659 for (i=0;i<nsubdomains;i++) { 1660 for (j=mat_ranges[i];j<mat_ranges[i+1];j++) { 1661 row_ownership[j] = i; 1662 } 1663 } 1664 ierr = MatGetOwnershipRangesColumn(A,(const PetscInt**)&mat_ranges);CHKERRQ(ierr); 1665 1666 /* 1667 my_dnz and my_onz contains exact contribution to preallocation from each local mat 1668 then, they will be summed up properly. This way, preallocation is always sufficient 1669 */ 1670 ierr = PetscCalloc2(local_rows,&my_dnz,local_rows,&my_onz);CHKERRQ(ierr); 1671 /* preallocation as a MATAIJ */ 1672 if (isdense) { /* special case for dense local matrices */ 1673 for (i=0;i<local_rows;i++) { 1674 PetscInt owner = row_ownership[global_indices_r[i]]; 1675 for (j=0;j<local_cols;j++) { 1676 PetscInt index_col = global_indices_c[j]; 1677 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 1678 my_dnz[i] += 1; 1679 } else { /* offdiag block */ 1680 my_onz[i] += 1; 1681 } 1682 } 1683 } 1684 } else if (matis->A->ops->getrowij) { 1685 const PetscInt *ii,*jj,*jptr; 1686 PetscBool done; 1687 ierr = MatGetRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);CHKERRQ(ierr); 1688 if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatGetRowIJ"); 1689 jptr = jj; 1690 for (i=0;i<local_rows;i++) { 1691 PetscInt index_row = global_indices_r[i]; 1692 for (j=0;j<ii[i+1]-ii[i];j++,jptr++) { 1693 PetscInt owner = row_ownership[index_row]; 1694 PetscInt index_col = global_indices_c[*jptr]; 1695 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 1696 my_dnz[i] += 1; 1697 } else { /* offdiag block */ 1698 my_onz[i] += 1; 1699 } 1700 /* same as before, interchanging rows and cols */ 1701 if (issbaij && index_col != index_row) { 1702 owner = row_ownership[index_col]; 1703 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) { 1704 my_dnz[*jptr] += 1; 1705 } else { 1706 my_onz[*jptr] += 1; 1707 } 1708 } 1709 } 1710 } 1711 ierr = MatRestoreRowIJ(matis->A,0,PETSC_FALSE,PETSC_FALSE,&local_rows,&ii,&jj,&done);CHKERRQ(ierr); 1712 if (!done) SETERRQ(PetscObjectComm((PetscObject)(matis->A)),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ"); 1713 } else { /* loop over rows and use MatGetRow */ 1714 for (i=0;i<local_rows;i++) { 1715 const PetscInt *cols; 1716 PetscInt ncols,index_row = global_indices_r[i]; 1717 ierr = MatGetRow(matis->A,i,&ncols,&cols,NULL);CHKERRQ(ierr); 1718 for (j=0;j<ncols;j++) { 1719 PetscInt owner = row_ownership[index_row]; 1720 PetscInt index_col = global_indices_c[cols[j]]; 1721 if (index_col > mat_ranges[owner]-1 && index_col < mat_ranges[owner+1] ) { /* diag block */ 1722 my_dnz[i] += 1; 1723 } else { /* offdiag block */ 1724 my_onz[i] += 1; 1725 } 1726 /* same as before, interchanging rows and cols */ 1727 if (issbaij && index_col != index_row) { 1728 owner = row_ownership[index_col]; 1729 if (index_row > mat_ranges[owner]-1 && index_row < mat_ranges[owner+1] ) { 1730 my_dnz[cols[j]] += 1; 1731 } else { 1732 my_onz[cols[j]] += 1; 1733 } 1734 } 1735 } 1736 ierr = MatRestoreRow(matis->A,i,&ncols,&cols,NULL);CHKERRQ(ierr); 1737 } 1738 } 1739 if (global_indices_c != global_indices_r) { 1740 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&global_indices_c);CHKERRQ(ierr); 1741 } 1742 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&global_indices_r);CHKERRQ(ierr); 1743 ierr = PetscFree(row_ownership);CHKERRQ(ierr); 1744 1745 /* Reduce my_dnz and my_onz */ 1746 if (maxreduce) { 1747 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);CHKERRQ(ierr); 1748 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);CHKERRQ(ierr); 1749 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_MAX);CHKERRQ(ierr); 1750 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_MAX);CHKERRQ(ierr); 1751 } else { 1752 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);CHKERRQ(ierr); 1753 ierr = PetscSFReduceBegin(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);CHKERRQ(ierr); 1754 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_dnz,dnz,MPI_SUM);CHKERRQ(ierr); 1755 ierr = PetscSFReduceEnd(matis->sf,MPIU_INT,my_onz,onz,MPI_SUM);CHKERRQ(ierr); 1756 } 1757 ierr = PetscFree2(my_dnz,my_onz);CHKERRQ(ierr); 1758 1759 /* Resize preallocation if overestimated */ 1760 for (i=0;i<lrows;i++) { 1761 dnz[i] = PetscMin(dnz[i],lcols); 1762 onz[i] = PetscMin(onz[i],cols-lcols); 1763 } 1764 1765 /* Set preallocation */ 1766 ierr = MatSeqAIJSetPreallocation(B,0,dnz);CHKERRQ(ierr); 1767 ierr = MatMPIAIJSetPreallocation(B,0,dnz,0,onz);CHKERRQ(ierr); 1768 for (i=0;i<lrows/bs;i++) { 1769 dnz[i] = dnz[i*bs]/bs; 1770 onz[i] = onz[i*bs]/bs; 1771 } 1772 ierr = MatSeqBAIJSetPreallocation(B,bs,0,dnz);CHKERRQ(ierr); 1773 ierr = MatMPIBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 1774 ierr = MatMPISBAIJSetPreallocation(B,bs,0,dnz,0,onz);CHKERRQ(ierr); 1775 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1776 if (issbaij) { 1777 ierr = MatRestoreRowUpperTriangular(matis->A);CHKERRQ(ierr); 1778 } 1779 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1780 PetscFunctionReturn(0); 1781 } 1782 1783 PETSC_INTERN PetscErrorCode MatConvert_IS_XAIJ(Mat mat, MatType mtype, MatReuse reuse, Mat *M) 1784 { 1785 Mat_IS *matis = (Mat_IS*)(mat->data); 1786 Mat local_mat,MT; 1787 /* info on mat */ 1788 PetscInt bs,rows,cols,lrows,lcols; 1789 PetscInt local_rows,local_cols; 1790 PetscBool isseqdense,isseqsbaij,isseqaij,isseqbaij; 1791 #if defined (PETSC_USE_DEBUG) 1792 PetscBool lb[4],bb[4]; 1793 #endif 1794 PetscMPIInt nsubdomains; 1795 /* values insertion */ 1796 PetscScalar *array; 1797 /* work */ 1798 PetscErrorCode ierr; 1799 1800 PetscFunctionBegin; 1801 /* get info from mat */ 1802 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&nsubdomains);CHKERRQ(ierr); 1803 if (nsubdomains == 1) { 1804 Mat B; 1805 IS rows,cols; 1806 IS irows,icols; 1807 const PetscInt *ridxs,*cidxs; 1808 PetscInt rbs,cbs; 1809 1810 ierr = ISLocalToGlobalMappingGetBlockIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 1811 ierr = ISLocalToGlobalMappingGetBlockIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 1812 ierr = ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);CHKERRQ(ierr); 1813 ierr = ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);CHKERRQ(ierr); 1814 ierr = ISCreateBlock(PETSC_COMM_SELF,rbs,mat->rmap->n,ridxs,PETSC_USE_POINTER,&rows);CHKERRQ(ierr); 1815 ierr = ISCreateBlock(PETSC_COMM_SELF,cbs,mat->cmap->n,cidxs,PETSC_USE_POINTER,&cols);CHKERRQ(ierr); 1816 ierr = ISSetPermutation(rows);CHKERRQ(ierr); 1817 ierr = ISSetPermutation(cols);CHKERRQ(ierr); 1818 ierr = ISInvertPermutation(rows,mat->rmap->n,&irows);CHKERRQ(ierr); 1819 ierr = ISInvertPermutation(cols,mat->cmap->n,&icols);CHKERRQ(ierr); 1820 ierr = ISLocalToGlobalMappingRestoreBlockIndices(mat->rmap->mapping,&ridxs);CHKERRQ(ierr); 1821 ierr = ISLocalToGlobalMappingRestoreBlockIndices(mat->cmap->mapping,&cidxs);CHKERRQ(ierr); 1822 ierr = ISDestroy(&cols);CHKERRQ(ierr); 1823 ierr = ISDestroy(&rows);CHKERRQ(ierr); 1824 ierr = MatConvert(matis->A,(rbs == cbs && rbs > 1) ? MATSEQBAIJ : MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 1825 if (reuse != MAT_INPLACE_MATRIX) { 1826 ierr = MatCreateSubMatrix(B,irows,icols,reuse,M);CHKERRQ(ierr); 1827 } else { 1828 Mat C; 1829 1830 ierr = MatCreateSubMatrix(B,irows,icols,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr); 1831 ierr = MatHeaderReplace(mat,&C);CHKERRQ(ierr); 1832 } 1833 ierr = MatDestroy(&B);CHKERRQ(ierr); 1834 ierr = ISDestroy(&icols);CHKERRQ(ierr); 1835 ierr = ISDestroy(&irows);CHKERRQ(ierr); 1836 PetscFunctionReturn(0); 1837 } 1838 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 1839 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 1840 ierr = MatGetLocalSize(mat,&lrows,&lcols);CHKERRQ(ierr); 1841 ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr); 1842 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 1843 ierr = PetscObjectBaseTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 1844 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQBAIJ,&isseqbaij);CHKERRQ(ierr); 1845 ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&isseqsbaij);CHKERRQ(ierr); 1846 if (!isseqdense && !isseqaij && !isseqbaij && !isseqsbaij) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for matrix type %s",((PetscObject)(matis->A))->type_name); 1847 #if defined (PETSC_USE_DEBUG) 1848 lb[0] = isseqdense; 1849 lb[1] = isseqaij; 1850 lb[2] = isseqbaij; 1851 lb[3] = isseqsbaij; 1852 ierr = MPIU_Allreduce(lb,bb,4,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 1853 if (!bb[0] && !bb[1] && !bb[2] && !bb[3]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Local matrices must have the same type"); 1854 #endif 1855 1856 if (reuse != MAT_REUSE_MATRIX) { 1857 PetscBool isaij; 1858 1859 ierr = MatCreate(PetscObjectComm((PetscObject)mat),&MT);CHKERRQ(ierr); 1860 ierr = MatSetSizes(MT,lrows,lcols,rows,cols);CHKERRQ(ierr); 1861 ierr = PetscStrcmp(mtype,MATAIJ,&isaij);CHKERRQ(ierr); 1862 if (isaij) { 1863 if (!isseqsbaij) { 1864 ierr = MatSetType(MT,MATAIJ);CHKERRQ(ierr); 1865 } else { 1866 ierr = MatSetType(MT,MATSBAIJ);CHKERRQ(ierr); 1867 } 1868 } else { 1869 ierr = MatSetType(MT,mtype);CHKERRQ(ierr); 1870 } 1871 ierr = MatSetBlockSize(MT,bs);CHKERRQ(ierr); 1872 ierr = MatISSetMPIXAIJPreallocation_Private(mat,MT,PETSC_FALSE);CHKERRQ(ierr); 1873 } else { 1874 PetscInt mbs,mrows,mcols,mlrows,mlcols; 1875 1876 /* some checks */ 1877 MT = *M; 1878 ierr = MatGetBlockSize(MT,&mbs);CHKERRQ(ierr); 1879 ierr = MatGetSize(MT,&mrows,&mcols);CHKERRQ(ierr); 1880 ierr = MatGetLocalSize(MT,&mlrows,&mlcols);CHKERRQ(ierr); 1881 if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows); 1882 if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols); 1883 if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows); 1884 if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols); 1885 if (mbs != bs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs); 1886 ierr = MatZeroEntries(MT);CHKERRQ(ierr); 1887 } 1888 1889 if (isseqsbaij) { 1890 ierr = MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);CHKERRQ(ierr); 1891 } else { 1892 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1893 local_mat = matis->A; 1894 } 1895 1896 /* Set values */ 1897 ierr = MatSetLocalToGlobalMapping(MT,mat->rmap->mapping,mat->cmap->mapping);CHKERRQ(ierr); 1898 if (isseqdense) { /* special case for dense local matrices */ 1899 PetscInt i,*dummy; 1900 1901 ierr = PetscMalloc1(PetscMax(local_rows,local_cols),&dummy);CHKERRQ(ierr); 1902 for (i=0;i<PetscMax(local_rows,local_cols);i++) dummy[i] = i; 1903 ierr = MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 1904 ierr = MatDenseGetArray(local_mat,&array);CHKERRQ(ierr); 1905 ierr = MatSetValuesLocal(MT,local_rows,dummy,local_cols,dummy,array,ADD_VALUES);CHKERRQ(ierr); 1906 ierr = MatDenseRestoreArray(local_mat,&array);CHKERRQ(ierr); 1907 ierr = PetscFree(dummy);CHKERRQ(ierr); 1908 } else if (isseqaij) { 1909 PetscInt i,nvtxs,*xadj,*adjncy; 1910 PetscBool done; 1911 1912 ierr = MatGetRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1913 if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatGetRowIJ"); 1914 ierr = MatSeqAIJGetArray(local_mat,&array);CHKERRQ(ierr); 1915 for (i=0;i<nvtxs;i++) { 1916 ierr = MatSetValuesLocal(MT,1,&i,xadj[i+1]-xadj[i],adjncy+xadj[i],array+xadj[i],ADD_VALUES);CHKERRQ(ierr); 1917 } 1918 ierr = MatRestoreRowIJ(local_mat,0,PETSC_FALSE,PETSC_FALSE,&nvtxs,(const PetscInt**)&xadj,(const PetscInt**)&adjncy,&done);CHKERRQ(ierr); 1919 if (!done) SETERRQ(PetscObjectComm((PetscObject)local_mat),PETSC_ERR_PLIB,"Error in MatRestoreRowIJ"); 1920 ierr = MatSeqAIJRestoreArray(local_mat,&array);CHKERRQ(ierr); 1921 } else { /* very basic values insertion for all other matrix types */ 1922 PetscInt i; 1923 1924 for (i=0;i<local_rows;i++) { 1925 PetscInt j; 1926 const PetscInt *local_indices_cols; 1927 1928 ierr = MatGetRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1929 ierr = MatSetValuesLocal(MT,1,&i,j,local_indices_cols,array,ADD_VALUES);CHKERRQ(ierr); 1930 ierr = MatRestoreRow(local_mat,i,&j,&local_indices_cols,(const PetscScalar**)&array);CHKERRQ(ierr); 1931 } 1932 } 1933 ierr = MatAssemblyBegin(MT,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1934 ierr = MatDestroy(&local_mat);CHKERRQ(ierr); 1935 ierr = MatAssemblyEnd(MT,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1936 if (isseqdense) { 1937 ierr = MatSetOption(MT,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 1938 } 1939 if (reuse == MAT_INPLACE_MATRIX) { 1940 ierr = MatHeaderReplace(mat,&MT);CHKERRQ(ierr); 1941 } else if (reuse == MAT_INITIAL_MATRIX) { 1942 *M = MT; 1943 } 1944 PetscFunctionReturn(0); 1945 } 1946 1947 /*@ 1948 MatISGetMPIXAIJ - Converts MATIS matrix into a parallel AIJ format 1949 1950 Input Parameter: 1951 . mat - the matrix (should be of type MATIS) 1952 . reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1953 1954 Output Parameter: 1955 . newmat - the matrix in AIJ format 1956 1957 Level: developer 1958 1959 Notes: 1960 This function has been deprecated and it will be removed in future releases. Update your code to use the MatConvert() interface. 1961 1962 .seealso: MATIS, MatConvert() 1963 @*/ 1964 PetscErrorCode MatISGetMPIXAIJ(Mat mat, MatReuse reuse, Mat *newmat) 1965 { 1966 PetscErrorCode ierr; 1967 1968 PetscFunctionBegin; 1969 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1970 PetscValidLogicalCollectiveEnum(mat,reuse,2); 1971 PetscValidPointer(newmat,3); 1972 if (reuse == MAT_REUSE_MATRIX) { 1973 PetscValidHeaderSpecific(*newmat,MAT_CLASSID,3); 1974 PetscCheckSameComm(mat,1,*newmat,3); 1975 if (mat == *newmat) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse the same matrix"); 1976 } 1977 ierr = PetscUseMethod(mat,"MatISGetMPIXAIJ_C",(Mat,MatType,MatReuse,Mat*),(mat,MATAIJ,reuse,newmat));CHKERRQ(ierr); 1978 PetscFunctionReturn(0); 1979 } 1980 1981 PetscErrorCode MatDuplicate_IS(Mat mat,MatDuplicateOption op,Mat *newmat) 1982 { 1983 PetscErrorCode ierr; 1984 Mat_IS *matis = (Mat_IS*)(mat->data); 1985 PetscInt rbs,cbs,m,n,M,N; 1986 Mat B,localmat; 1987 1988 PetscFunctionBegin; 1989 ierr = ISLocalToGlobalMappingGetBlockSize(mat->rmap->mapping,&rbs);CHKERRQ(ierr); 1990 ierr = ISLocalToGlobalMappingGetBlockSize(mat->cmap->mapping,&cbs);CHKERRQ(ierr); 1991 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1992 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 1993 ierr = MatCreateIS(PetscObjectComm((PetscObject)mat),rbs == cbs ? rbs : 1,m,n,M,N,mat->rmap->mapping,mat->cmap->mapping,&B);CHKERRQ(ierr); 1994 ierr = MatDuplicate(matis->A,op,&localmat);CHKERRQ(ierr); 1995 ierr = MatISSetLocalMat(B,localmat);CHKERRQ(ierr); 1996 ierr = MatDestroy(&localmat);CHKERRQ(ierr); 1997 if (matis->sf) { 1998 Mat_IS *bmatis = (Mat_IS*)(B->data); 1999 2000 ierr = PetscObjectReference((PetscObject)matis->sf);CHKERRQ(ierr); 2001 bmatis->sf = matis->sf; 2002 ierr = PetscMalloc2(matis->sf->nroots,&bmatis->sf_rootdata,matis->sf->nleaves,&bmatis->sf_leafdata);CHKERRQ(ierr); 2003 if (matis->sf != matis->csf) { 2004 ierr = PetscObjectReference((PetscObject)matis->csf);CHKERRQ(ierr); 2005 bmatis->csf = matis->csf; 2006 ierr = PetscMalloc2(matis->csf->nroots,&bmatis->csf_rootdata,matis->csf->nleaves,&bmatis->csf_leafdata);CHKERRQ(ierr); 2007 } else { 2008 bmatis->csf = bmatis->sf; 2009 bmatis->csf_leafdata = bmatis->sf_leafdata; 2010 bmatis->csf_rootdata = bmatis->sf_rootdata; 2011 } 2012 } 2013 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2014 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2015 *newmat = B; 2016 PetscFunctionReturn(0); 2017 } 2018 2019 static PetscErrorCode MatIsHermitian_IS(Mat A,PetscReal tol,PetscBool *flg) 2020 { 2021 PetscErrorCode ierr; 2022 Mat_IS *matis = (Mat_IS*)A->data; 2023 PetscBool local_sym; 2024 2025 PetscFunctionBegin; 2026 ierr = MatIsHermitian(matis->A,tol,&local_sym);CHKERRQ(ierr); 2027 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2028 PetscFunctionReturn(0); 2029 } 2030 2031 static PetscErrorCode MatIsSymmetric_IS(Mat A,PetscReal tol,PetscBool *flg) 2032 { 2033 PetscErrorCode ierr; 2034 Mat_IS *matis = (Mat_IS*)A->data; 2035 PetscBool local_sym; 2036 2037 PetscFunctionBegin; 2038 ierr = MatIsSymmetric(matis->A,tol,&local_sym);CHKERRQ(ierr); 2039 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2040 PetscFunctionReturn(0); 2041 } 2042 2043 static PetscErrorCode MatIsStructurallySymmetric_IS(Mat A,PetscBool *flg) 2044 { 2045 PetscErrorCode ierr; 2046 Mat_IS *matis = (Mat_IS*)A->data; 2047 PetscBool local_sym; 2048 2049 PetscFunctionBegin; 2050 if (A->rmap->mapping != A->cmap->mapping) { 2051 *flg = PETSC_FALSE; 2052 PetscFunctionReturn(0); 2053 } 2054 ierr = MatIsStructurallySymmetric(matis->A,&local_sym);CHKERRQ(ierr); 2055 ierr = MPIU_Allreduce(&local_sym,flg,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2056 PetscFunctionReturn(0); 2057 } 2058 2059 static PetscErrorCode MatDestroy_IS(Mat A) 2060 { 2061 PetscErrorCode ierr; 2062 Mat_IS *b = (Mat_IS*)A->data; 2063 2064 PetscFunctionBegin; 2065 ierr = MatDestroy(&b->A);CHKERRQ(ierr); 2066 ierr = VecScatterDestroy(&b->cctx);CHKERRQ(ierr); 2067 ierr = VecScatterDestroy(&b->rctx);CHKERRQ(ierr); 2068 ierr = VecDestroy(&b->x);CHKERRQ(ierr); 2069 ierr = VecDestroy(&b->y);CHKERRQ(ierr); 2070 ierr = VecDestroy(&b->counter);CHKERRQ(ierr); 2071 ierr = ISDestroy(&b->getsub_ris);CHKERRQ(ierr); 2072 ierr = ISDestroy(&b->getsub_cis);CHKERRQ(ierr); 2073 if (b->sf != b->csf) { 2074 ierr = PetscSFDestroy(&b->csf);CHKERRQ(ierr); 2075 ierr = PetscFree2(b->csf_rootdata,b->csf_leafdata);CHKERRQ(ierr); 2076 } 2077 ierr = PetscSFDestroy(&b->sf);CHKERRQ(ierr); 2078 ierr = PetscFree2(b->sf_rootdata,b->sf_leafdata);CHKERRQ(ierr); 2079 ierr = PetscFree(A->data);CHKERRQ(ierr); 2080 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 2081 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",NULL);CHKERRQ(ierr); 2082 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",NULL);CHKERRQ(ierr); 2083 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",NULL);CHKERRQ(ierr); 2084 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",NULL);CHKERRQ(ierr); 2085 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",NULL);CHKERRQ(ierr); 2086 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",NULL);CHKERRQ(ierr); 2087 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",NULL);CHKERRQ(ierr); 2088 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",NULL);CHKERRQ(ierr); 2089 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",NULL);CHKERRQ(ierr); 2090 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",NULL);CHKERRQ(ierr); 2091 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",NULL);CHKERRQ(ierr); 2092 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",NULL);CHKERRQ(ierr); 2093 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",NULL);CHKERRQ(ierr); 2094 PetscFunctionReturn(0); 2095 } 2096 2097 static PetscErrorCode MatMult_IS(Mat A,Vec x,Vec y) 2098 { 2099 PetscErrorCode ierr; 2100 Mat_IS *is = (Mat_IS*)A->data; 2101 PetscScalar zero = 0.0; 2102 2103 PetscFunctionBegin; 2104 /* scatter the global vector x into the local work vector */ 2105 ierr = VecScatterBegin(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2106 ierr = VecScatterEnd(is->cctx,x,is->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2107 2108 /* multiply the local matrix */ 2109 ierr = MatMult(is->A,is->x,is->y);CHKERRQ(ierr); 2110 2111 /* scatter product back into global memory */ 2112 ierr = VecSet(y,zero);CHKERRQ(ierr); 2113 ierr = VecScatterBegin(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2114 ierr = VecScatterEnd(is->rctx,is->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2115 PetscFunctionReturn(0); 2116 } 2117 2118 static PetscErrorCode MatMultAdd_IS(Mat A,Vec v1,Vec v2,Vec v3) 2119 { 2120 Vec temp_vec; 2121 PetscErrorCode ierr; 2122 2123 PetscFunctionBegin; /* v3 = v2 + A * v1.*/ 2124 if (v3 != v2) { 2125 ierr = MatMult(A,v1,v3);CHKERRQ(ierr); 2126 ierr = VecAXPY(v3,1.0,v2);CHKERRQ(ierr); 2127 } else { 2128 ierr = VecDuplicate(v2,&temp_vec);CHKERRQ(ierr); 2129 ierr = MatMult(A,v1,temp_vec);CHKERRQ(ierr); 2130 ierr = VecAXPY(temp_vec,1.0,v2);CHKERRQ(ierr); 2131 ierr = VecCopy(temp_vec,v3);CHKERRQ(ierr); 2132 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 2133 } 2134 PetscFunctionReturn(0); 2135 } 2136 2137 static PetscErrorCode MatMultTranspose_IS(Mat A,Vec y,Vec x) 2138 { 2139 Mat_IS *is = (Mat_IS*)A->data; 2140 PetscErrorCode ierr; 2141 2142 PetscFunctionBegin; 2143 /* scatter the global vector x into the local work vector */ 2144 ierr = VecScatterBegin(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2145 ierr = VecScatterEnd(is->rctx,y,is->y,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2146 2147 /* multiply the local matrix */ 2148 ierr = MatMultTranspose(is->A,is->y,is->x);CHKERRQ(ierr); 2149 2150 /* scatter product back into global vector */ 2151 ierr = VecSet(x,0);CHKERRQ(ierr); 2152 ierr = VecScatterBegin(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2153 ierr = VecScatterEnd(is->cctx,is->x,x,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2154 PetscFunctionReturn(0); 2155 } 2156 2157 static PetscErrorCode MatMultTransposeAdd_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 = MatMultTranspose(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 = MatMultTranspose(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 MatView_IS(Mat A,PetscViewer viewer) 2177 { 2178 Mat_IS *a = (Mat_IS*)A->data; 2179 PetscErrorCode ierr; 2180 PetscViewer sviewer; 2181 PetscBool isascii,view = PETSC_TRUE; 2182 2183 PetscFunctionBegin; 2184 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 2185 if (isascii) { 2186 PetscViewerFormat format; 2187 2188 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 2189 if (format == PETSC_VIEWER_ASCII_INFO) view = PETSC_FALSE; 2190 } 2191 if (!view) PetscFunctionReturn(0); 2192 ierr = PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 2193 ierr = MatView(a->A,sviewer);CHKERRQ(ierr); 2194 ierr = PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);CHKERRQ(ierr); 2195 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2196 PetscFunctionReturn(0); 2197 } 2198 2199 static PetscErrorCode MatSetLocalToGlobalMapping_IS(Mat A,ISLocalToGlobalMapping rmapping,ISLocalToGlobalMapping cmapping) 2200 { 2201 PetscErrorCode ierr; 2202 PetscInt nr,rbs,nc,cbs; 2203 Mat_IS *is = (Mat_IS*)A->data; 2204 IS from,to; 2205 Vec cglobal,rglobal; 2206 2207 PetscFunctionBegin; 2208 PetscCheckSameComm(A,1,rmapping,2); 2209 PetscCheckSameComm(A,1,cmapping,3); 2210 /* Destroy any previous data */ 2211 ierr = VecDestroy(&is->x);CHKERRQ(ierr); 2212 ierr = VecDestroy(&is->y);CHKERRQ(ierr); 2213 ierr = VecDestroy(&is->counter);CHKERRQ(ierr); 2214 ierr = VecScatterDestroy(&is->rctx);CHKERRQ(ierr); 2215 ierr = VecScatterDestroy(&is->cctx);CHKERRQ(ierr); 2216 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 2217 if (is->csf != is->sf) { 2218 ierr = PetscSFDestroy(&is->csf);CHKERRQ(ierr); 2219 ierr = PetscFree2(is->csf_rootdata,is->csf_leafdata);CHKERRQ(ierr); 2220 } 2221 ierr = PetscSFDestroy(&is->sf);CHKERRQ(ierr); 2222 ierr = PetscFree2(is->sf_rootdata,is->sf_leafdata);CHKERRQ(ierr); 2223 2224 /* Setup Layout and set local to global maps */ 2225 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 2226 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 2227 ierr = ISLocalToGlobalMappingGetSize(rmapping,&nr);CHKERRQ(ierr); 2228 ierr = ISLocalToGlobalMappingGetBlockSize(rmapping,&rbs);CHKERRQ(ierr); 2229 ierr = ISLocalToGlobalMappingGetSize(cmapping,&nc);CHKERRQ(ierr); 2230 ierr = ISLocalToGlobalMappingGetBlockSize(cmapping,&cbs);CHKERRQ(ierr); 2231 /* check if the two mappings are actually the same for square matrices (DOLFIN passes 2 different objects) */ 2232 if (rmapping != cmapping && A->rmap->N == A->cmap->N) { 2233 PetscBool same,gsame; 2234 2235 same = PETSC_FALSE; 2236 if (nr == nc && cbs == rbs) { 2237 const PetscInt *idxs1,*idxs2; 2238 2239 ierr = ISLocalToGlobalMappingGetBlockIndices(rmapping,&idxs1);CHKERRQ(ierr); 2240 ierr = ISLocalToGlobalMappingGetBlockIndices(cmapping,&idxs2);CHKERRQ(ierr); 2241 ierr = PetscMemcmp(idxs1,idxs2,(nr/rbs)*sizeof(PetscInt),&same);CHKERRQ(ierr); 2242 ierr = ISLocalToGlobalMappingRestoreBlockIndices(rmapping,&idxs1);CHKERRQ(ierr); 2243 ierr = ISLocalToGlobalMappingRestoreBlockIndices(cmapping,&idxs2);CHKERRQ(ierr); 2244 } 2245 ierr = MPIU_Allreduce(&same,&gsame,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));CHKERRQ(ierr); 2246 if (gsame) cmapping = rmapping; 2247 } 2248 ierr = PetscLayoutSetISLocalToGlobalMapping(A->rmap,rmapping);CHKERRQ(ierr); 2249 ierr = PetscLayoutSetISLocalToGlobalMapping(A->cmap,cmapping);CHKERRQ(ierr); 2250 2251 /* Create the local matrix A */ 2252 ierr = MatCreate(PETSC_COMM_SELF,&is->A);CHKERRQ(ierr); 2253 ierr = MatSetType(is->A,MATAIJ);CHKERRQ(ierr); 2254 ierr = MatSetSizes(is->A,nr,nc,nr,nc);CHKERRQ(ierr); 2255 ierr = MatSetBlockSizes(is->A,rbs,cbs);CHKERRQ(ierr); 2256 ierr = MatSetOptionsPrefix(is->A,((PetscObject)A)->prefix);CHKERRQ(ierr); 2257 ierr = MatAppendOptionsPrefix(is->A,"is_");CHKERRQ(ierr); 2258 ierr = MatSetFromOptions(is->A);CHKERRQ(ierr); 2259 ierr = PetscLayoutSetUp(is->A->rmap);CHKERRQ(ierr); 2260 ierr = PetscLayoutSetUp(is->A->cmap);CHKERRQ(ierr); 2261 2262 if (!is->islocalref) { /* setup scatters and local vectors for MatMult */ 2263 /* Create the local work vectors */ 2264 ierr = MatCreateVecs(is->A,&is->x,&is->y);CHKERRQ(ierr); 2265 2266 /* setup the global to local scatters */ 2267 ierr = MatCreateVecs(A,&cglobal,&rglobal);CHKERRQ(ierr); 2268 ierr = ISCreateStride(PETSC_COMM_SELF,nr,0,1,&to);CHKERRQ(ierr); 2269 ierr = ISLocalToGlobalMappingApplyIS(rmapping,to,&from);CHKERRQ(ierr); 2270 ierr = VecScatterCreate(rglobal,from,is->y,to,&is->rctx);CHKERRQ(ierr); 2271 if (rmapping != cmapping) { 2272 ierr = ISDestroy(&to);CHKERRQ(ierr); 2273 ierr = ISDestroy(&from);CHKERRQ(ierr); 2274 ierr = ISCreateStride(PETSC_COMM_SELF,nc,0,1,&to);CHKERRQ(ierr); 2275 ierr = ISLocalToGlobalMappingApplyIS(cmapping,to,&from);CHKERRQ(ierr); 2276 ierr = VecScatterCreate(cglobal,from,is->x,to,&is->cctx);CHKERRQ(ierr); 2277 } else { 2278 ierr = PetscObjectReference((PetscObject)is->rctx);CHKERRQ(ierr); 2279 is->cctx = is->rctx; 2280 } 2281 2282 /* interface counter vector (local) */ 2283 ierr = VecDuplicate(is->y,&is->counter);CHKERRQ(ierr); 2284 ierr = VecSet(is->y,1.);CHKERRQ(ierr); 2285 ierr = VecScatterBegin(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2286 ierr = VecScatterEnd(is->rctx,is->y,rglobal,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2287 ierr = VecScatterBegin(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2288 ierr = VecScatterEnd(is->rctx,rglobal,is->counter,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2289 2290 /* free workspace */ 2291 ierr = VecDestroy(&rglobal);CHKERRQ(ierr); 2292 ierr = VecDestroy(&cglobal);CHKERRQ(ierr); 2293 ierr = ISDestroy(&to);CHKERRQ(ierr); 2294 ierr = ISDestroy(&from);CHKERRQ(ierr); 2295 } 2296 ierr = MatSetUp(A);CHKERRQ(ierr); 2297 PetscFunctionReturn(0); 2298 } 2299 2300 static PetscErrorCode MatSetValues_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2301 { 2302 Mat_IS *is = (Mat_IS*)mat->data; 2303 PetscErrorCode ierr; 2304 #if defined(PETSC_USE_DEBUG) 2305 PetscInt i,zm,zn; 2306 #endif 2307 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2308 2309 PetscFunctionBegin; 2310 #if defined(PETSC_USE_DEBUG) 2311 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); 2312 /* count negative indices */ 2313 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 2314 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 2315 #endif 2316 ierr = ISGlobalToLocalMappingApply(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 2317 ierr = ISGlobalToLocalMappingApply(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 2318 #if defined(PETSC_USE_DEBUG) 2319 /* count negative indices (should be the same as before) */ 2320 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 2321 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 2322 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"); 2323 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"); 2324 #endif 2325 ierr = MatSetValues(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 2326 PetscFunctionReturn(0); 2327 } 2328 2329 static PetscErrorCode MatSetValuesBlocked_IS(Mat mat, PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols, const PetscScalar *values, InsertMode addv) 2330 { 2331 Mat_IS *is = (Mat_IS*)mat->data; 2332 PetscErrorCode ierr; 2333 #if defined(PETSC_USE_DEBUG) 2334 PetscInt i,zm,zn; 2335 #endif 2336 PetscInt rows_l[MATIS_MAX_ENTRIES_INSERTION],cols_l[MATIS_MAX_ENTRIES_INSERTION]; 2337 2338 PetscFunctionBegin; 2339 #if defined(PETSC_USE_DEBUG) 2340 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); 2341 /* count negative indices */ 2342 for (i=0,zm=0;i<m;i++) if (rows[i] < 0) zm++; 2343 for (i=0,zn=0;i<n;i++) if (cols[i] < 0) zn++; 2344 #endif 2345 ierr = ISGlobalToLocalMappingApplyBlock(mat->rmap->mapping,IS_GTOLM_MASK,m,rows,&m,rows_l);CHKERRQ(ierr); 2346 ierr = ISGlobalToLocalMappingApplyBlock(mat->cmap->mapping,IS_GTOLM_MASK,n,cols,&n,cols_l);CHKERRQ(ierr); 2347 #if defined(PETSC_USE_DEBUG) 2348 /* count negative indices (should be the same as before) */ 2349 for (i=0;i<m;i++) if (rows_l[i] < 0) zm--; 2350 for (i=0;i<n;i++) if (cols_l[i] < 0) zn--; 2351 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"); 2352 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"); 2353 #endif 2354 ierr = MatSetValuesBlocked(is->A,m,rows_l,n,cols_l,values,addv);CHKERRQ(ierr); 2355 PetscFunctionReturn(0); 2356 } 2357 2358 static PetscErrorCode MatSetValuesLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2359 { 2360 PetscErrorCode ierr; 2361 Mat_IS *is = (Mat_IS*)A->data; 2362 2363 PetscFunctionBegin; 2364 if (is->A->rmap->mapping) { 2365 ierr = MatSetValuesLocal(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2366 } else { 2367 ierr = MatSetValues(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2368 } 2369 PetscFunctionReturn(0); 2370 } 2371 2372 static PetscErrorCode MatSetValuesBlockedLocal_IS(Mat A,PetscInt m,const PetscInt *rows, PetscInt n,const PetscInt *cols,const PetscScalar *values,InsertMode addv) 2373 { 2374 PetscErrorCode ierr; 2375 Mat_IS *is = (Mat_IS*)A->data; 2376 2377 PetscFunctionBegin; 2378 if (is->A->rmap->mapping) { 2379 #if defined(PETSC_USE_DEBUG) 2380 PetscInt ibs,bs; 2381 2382 ierr = ISLocalToGlobalMappingGetBlockSize(is->A->rmap->mapping,&ibs);CHKERRQ(ierr); 2383 ierr = MatGetBlockSize(is->A,&bs);CHKERRQ(ierr); 2384 if (ibs != bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Different block sizes! local mat %D, local l2g map %D",bs,ibs); 2385 #endif 2386 ierr = MatSetValuesBlockedLocal(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2387 } else { 2388 ierr = MatSetValuesBlocked(is->A,m,rows,n,cols,values,addv);CHKERRQ(ierr); 2389 } 2390 PetscFunctionReturn(0); 2391 } 2392 2393 static PetscErrorCode MatISZeroRowsColumnsLocal_Private(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,PetscBool columns) 2394 { 2395 Mat_IS *is = (Mat_IS*)A->data; 2396 PetscErrorCode ierr; 2397 2398 PetscFunctionBegin; 2399 if (!n) { 2400 is->pure_neumann = PETSC_TRUE; 2401 } else { 2402 PetscInt i; 2403 is->pure_neumann = PETSC_FALSE; 2404 2405 if (columns) { 2406 ierr = MatZeroRowsColumns(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 2407 } else { 2408 ierr = MatZeroRows(is->A,n,rows,diag,0,0);CHKERRQ(ierr); 2409 } 2410 if (diag != 0.) { 2411 const PetscScalar *array; 2412 ierr = VecGetArrayRead(is->counter,&array);CHKERRQ(ierr); 2413 for (i=0; i<n; i++) { 2414 ierr = MatSetValue(is->A,rows[i],rows[i],diag/(array[rows[i]]),INSERT_VALUES);CHKERRQ(ierr); 2415 } 2416 ierr = VecRestoreArrayRead(is->counter,&array);CHKERRQ(ierr); 2417 } 2418 ierr = MatAssemblyBegin(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2419 ierr = MatAssemblyEnd(is->A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2420 } 2421 PetscFunctionReturn(0); 2422 } 2423 2424 static PetscErrorCode MatZeroRowsColumns_Private_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b,PetscBool columns) 2425 { 2426 Mat_IS *matis = (Mat_IS*)A->data; 2427 PetscInt nr,nl,len,i; 2428 PetscInt *lrows; 2429 PetscErrorCode ierr; 2430 2431 PetscFunctionBegin; 2432 #if defined(PETSC_USE_DEBUG) 2433 if (columns || diag != 0. || (x && b)) { 2434 PetscBool cong; 2435 2436 ierr = PetscLayoutCompare(A->rmap,A->cmap,&cong);CHKERRQ(ierr); 2437 cong = (PetscBool)(cong && matis->sf == matis->csf); 2438 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"); 2439 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"); 2440 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"); 2441 } 2442 #endif 2443 /* get locally owned rows */ 2444 ierr = PetscLayoutMapLocal_Private(A->rmap,n,rows,&len,&lrows,NULL);CHKERRQ(ierr); 2445 /* fix right hand side if needed */ 2446 if (x && b) { 2447 const PetscScalar *xx; 2448 PetscScalar *bb; 2449 2450 ierr = VecGetArrayRead(x, &xx);CHKERRQ(ierr); 2451 ierr = VecGetArray(b, &bb);CHKERRQ(ierr); 2452 for (i=0;i<len;++i) bb[lrows[i]] = diag*xx[lrows[i]]; 2453 ierr = VecRestoreArrayRead(x, &xx);CHKERRQ(ierr); 2454 ierr = VecRestoreArray(b, &bb);CHKERRQ(ierr); 2455 } 2456 /* get rows associated to the local matrices */ 2457 ierr = MatISSetUpSF(A);CHKERRQ(ierr); 2458 ierr = MatGetSize(matis->A,&nl,NULL);CHKERRQ(ierr); 2459 ierr = PetscMemzero(matis->sf_leafdata,nl*sizeof(PetscInt));CHKERRQ(ierr); 2460 ierr = PetscMemzero(matis->sf_rootdata,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2461 for (i=0;i<len;i++) matis->sf_rootdata[lrows[i]] = 1; 2462 ierr = PetscFree(lrows);CHKERRQ(ierr); 2463 ierr = PetscSFBcastBegin(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 2464 ierr = PetscSFBcastEnd(matis->sf,MPIU_INT,matis->sf_rootdata,matis->sf_leafdata);CHKERRQ(ierr); 2465 ierr = PetscMalloc1(nl,&lrows);CHKERRQ(ierr); 2466 for (i=0,nr=0;i<nl;i++) if (matis->sf_leafdata[i]) lrows[nr++] = i; 2467 ierr = MatISZeroRowsColumnsLocal_Private(A,nr,lrows,diag,columns);CHKERRQ(ierr); 2468 ierr = PetscFree(lrows);CHKERRQ(ierr); 2469 PetscFunctionReturn(0); 2470 } 2471 2472 static PetscErrorCode MatZeroRows_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2473 { 2474 PetscErrorCode ierr; 2475 2476 PetscFunctionBegin; 2477 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_FALSE);CHKERRQ(ierr); 2478 PetscFunctionReturn(0); 2479 } 2480 2481 static PetscErrorCode MatZeroRowsColumns_IS(Mat A,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2482 { 2483 PetscErrorCode ierr; 2484 2485 PetscFunctionBegin; 2486 ierr = MatZeroRowsColumns_Private_IS(A,n,rows,diag,x,b,PETSC_TRUE);CHKERRQ(ierr); 2487 PetscFunctionReturn(0); 2488 } 2489 2490 static PetscErrorCode MatAssemblyBegin_IS(Mat A,MatAssemblyType type) 2491 { 2492 Mat_IS *is = (Mat_IS*)A->data; 2493 PetscErrorCode ierr; 2494 2495 PetscFunctionBegin; 2496 ierr = MatAssemblyBegin(is->A,type);CHKERRQ(ierr); 2497 PetscFunctionReturn(0); 2498 } 2499 2500 static PetscErrorCode MatAssemblyEnd_IS(Mat A,MatAssemblyType type) 2501 { 2502 Mat_IS *is = (Mat_IS*)A->data; 2503 PetscErrorCode ierr; 2504 2505 PetscFunctionBegin; 2506 ierr = MatAssemblyEnd(is->A,type);CHKERRQ(ierr); 2507 /* fix for local empty rows/cols */ 2508 if (is->locempty && type == MAT_FINAL_ASSEMBLY) { 2509 Mat newlA; 2510 ISLocalToGlobalMapping l2g; 2511 IS tis; 2512 const PetscScalar *v; 2513 PetscInt i,n,cf,*loce,*locf; 2514 PetscBool sym; 2515 2516 ierr = MatGetRowMaxAbs(is->A,is->y,NULL);CHKERRQ(ierr); 2517 ierr = MatIsSymmetric(is->A,PETSC_SMALL,&sym);CHKERRQ(ierr); 2518 if (!sym) SETERRQ(PetscObjectComm((PetscObject)is->A),PETSC_ERR_SUP,"Not yet implemented for unsymmetric case"); 2519 ierr = VecGetLocalSize(is->y,&n);CHKERRQ(ierr); 2520 ierr = PetscMalloc1(n,&loce);CHKERRQ(ierr); 2521 ierr = PetscMalloc1(n,&locf);CHKERRQ(ierr); 2522 ierr = VecGetArrayRead(is->y,&v);CHKERRQ(ierr); 2523 for (i=0,cf=0;i<n;i++) { 2524 if (v[i] == 0.0) { 2525 loce[i] = -1; 2526 } else { 2527 loce[i] = cf; 2528 locf[cf++] = i; 2529 } 2530 } 2531 ierr = VecRestoreArrayRead(is->y,&v);CHKERRQ(ierr); 2532 /* extract valid submatrix */ 2533 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)is->A),cf,locf,PETSC_USE_POINTER,&tis);CHKERRQ(ierr); 2534 ierr = MatCreateSubMatrix(is->A,tis,tis,MAT_INITIAL_MATRIX,&newlA);CHKERRQ(ierr); 2535 ierr = ISDestroy(&tis);CHKERRQ(ierr); 2536 /* attach local l2g map for successive calls of MatSetValues */ 2537 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)is->A),1,n,loce,PETSC_OWN_POINTER,&l2g);CHKERRQ(ierr); 2538 ierr = MatSetLocalToGlobalMapping(newlA,l2g,l2g);CHKERRQ(ierr); 2539 ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr); 2540 /* attach new global l2g map */ 2541 ierr = ISLocalToGlobalMappingApply(A->rmap->mapping,cf,locf,locf);CHKERRQ(ierr); 2542 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,cf,locf,PETSC_OWN_POINTER,&l2g);CHKERRQ(ierr); 2543 ierr = MatSetLocalToGlobalMapping(A,l2g,l2g);CHKERRQ(ierr); 2544 ierr = MatISSetLocalMat(A,newlA);CHKERRQ(ierr); 2545 ierr = MatDestroy(&newlA);CHKERRQ(ierr); 2546 ierr = ISLocalToGlobalMappingDestroy(&l2g);CHKERRQ(ierr); 2547 } 2548 is->locempty = PETSC_FALSE; 2549 PetscFunctionReturn(0); 2550 } 2551 2552 static PetscErrorCode MatISGetLocalMat_IS(Mat mat,Mat *local) 2553 { 2554 Mat_IS *is = (Mat_IS*)mat->data; 2555 2556 PetscFunctionBegin; 2557 *local = is->A; 2558 PetscFunctionReturn(0); 2559 } 2560 2561 static PetscErrorCode MatISRestoreLocalMat_IS(Mat mat,Mat *local) 2562 { 2563 PetscFunctionBegin; 2564 *local = NULL; 2565 PetscFunctionReturn(0); 2566 } 2567 2568 /*@ 2569 MatISGetLocalMat - Gets the local matrix stored inside a MATIS matrix. 2570 2571 Input Parameter: 2572 . mat - the matrix 2573 2574 Output Parameter: 2575 . local - the local matrix 2576 2577 Level: advanced 2578 2579 Notes: 2580 This can be called if you have precomputed the nonzero structure of the 2581 matrix and want to provide it to the inner matrix object to improve the performance 2582 of the MatSetValues() operation. 2583 2584 Call MatISRestoreLocalMat() when finished with the local matrix. 2585 2586 .seealso: MATIS 2587 @*/ 2588 PetscErrorCode MatISGetLocalMat(Mat mat,Mat *local) 2589 { 2590 PetscErrorCode ierr; 2591 2592 PetscFunctionBegin; 2593 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2594 PetscValidPointer(local,2); 2595 ierr = PetscUseMethod(mat,"MatISGetLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 2596 PetscFunctionReturn(0); 2597 } 2598 2599 /*@ 2600 MatISRestoreLocalMat - Restores the local matrix obtained with MatISGetLocalMat() 2601 2602 Input Parameter: 2603 . mat - the matrix 2604 2605 Output Parameter: 2606 . local - the local matrix 2607 2608 Level: advanced 2609 2610 .seealso: MATIS 2611 @*/ 2612 PetscErrorCode MatISRestoreLocalMat(Mat mat,Mat *local) 2613 { 2614 PetscErrorCode ierr; 2615 2616 PetscFunctionBegin; 2617 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2618 PetscValidPointer(local,2); 2619 ierr = PetscUseMethod(mat,"MatISRestoreLocalMat_C",(Mat,Mat*),(mat,local));CHKERRQ(ierr); 2620 PetscFunctionReturn(0); 2621 } 2622 2623 static PetscErrorCode MatISSetLocalMat_IS(Mat mat,Mat local) 2624 { 2625 Mat_IS *is = (Mat_IS*)mat->data; 2626 PetscInt nrows,ncols,orows,ocols; 2627 PetscErrorCode ierr; 2628 2629 PetscFunctionBegin; 2630 if (is->A) { 2631 ierr = MatGetSize(is->A,&orows,&ocols);CHKERRQ(ierr); 2632 ierr = MatGetSize(local,&nrows,&ncols);CHKERRQ(ierr); 2633 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); 2634 } 2635 ierr = PetscObjectReference((PetscObject)local);CHKERRQ(ierr); 2636 ierr = MatDestroy(&is->A);CHKERRQ(ierr); 2637 is->A = local; 2638 PetscFunctionReturn(0); 2639 } 2640 2641 /*@ 2642 MatISSetLocalMat - Replace the local matrix stored inside a MATIS object. 2643 2644 Input Parameter: 2645 . mat - the matrix 2646 . local - the local matrix 2647 2648 Output Parameter: 2649 2650 Level: advanced 2651 2652 Notes: 2653 This can be called if you have precomputed the local matrix and 2654 want to provide it to the matrix object MATIS. 2655 2656 .seealso: MATIS 2657 @*/ 2658 PetscErrorCode MatISSetLocalMat(Mat mat,Mat local) 2659 { 2660 PetscErrorCode ierr; 2661 2662 PetscFunctionBegin; 2663 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2664 PetscValidHeaderSpecific(local,MAT_CLASSID,2); 2665 ierr = PetscUseMethod(mat,"MatISSetLocalMat_C",(Mat,Mat),(mat,local));CHKERRQ(ierr); 2666 PetscFunctionReturn(0); 2667 } 2668 2669 static PetscErrorCode MatZeroEntries_IS(Mat A) 2670 { 2671 Mat_IS *a = (Mat_IS*)A->data; 2672 PetscErrorCode ierr; 2673 2674 PetscFunctionBegin; 2675 ierr = MatZeroEntries(a->A);CHKERRQ(ierr); 2676 PetscFunctionReturn(0); 2677 } 2678 2679 static PetscErrorCode MatScale_IS(Mat A,PetscScalar a) 2680 { 2681 Mat_IS *is = (Mat_IS*)A->data; 2682 PetscErrorCode ierr; 2683 2684 PetscFunctionBegin; 2685 ierr = MatScale(is->A,a);CHKERRQ(ierr); 2686 PetscFunctionReturn(0); 2687 } 2688 2689 static PetscErrorCode MatGetDiagonal_IS(Mat A, Vec v) 2690 { 2691 Mat_IS *is = (Mat_IS*)A->data; 2692 PetscErrorCode ierr; 2693 2694 PetscFunctionBegin; 2695 /* get diagonal of the local matrix */ 2696 ierr = MatGetDiagonal(is->A,is->y);CHKERRQ(ierr); 2697 2698 /* scatter diagonal back into global vector */ 2699 ierr = VecSet(v,0);CHKERRQ(ierr); 2700 ierr = VecScatterBegin(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2701 ierr = VecScatterEnd(is->rctx,is->y,v,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2702 PetscFunctionReturn(0); 2703 } 2704 2705 static PetscErrorCode MatSetOption_IS(Mat A,MatOption op,PetscBool flg) 2706 { 2707 Mat_IS *a = (Mat_IS*)A->data; 2708 PetscErrorCode ierr; 2709 2710 PetscFunctionBegin; 2711 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 2712 PetscFunctionReturn(0); 2713 } 2714 2715 static PetscErrorCode MatAXPY_IS(Mat Y,PetscScalar a,Mat X,MatStructure str) 2716 { 2717 Mat_IS *y = (Mat_IS*)Y->data; 2718 Mat_IS *x; 2719 #if defined(PETSC_USE_DEBUG) 2720 PetscBool ismatis; 2721 #endif 2722 PetscErrorCode ierr; 2723 2724 PetscFunctionBegin; 2725 #if defined(PETSC_USE_DEBUG) 2726 ierr = PetscObjectTypeCompare((PetscObject)X,MATIS,&ismatis);CHKERRQ(ierr); 2727 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot call MatAXPY(Y,a,X,str) with X not of type MATIS"); 2728 #endif 2729 x = (Mat_IS*)X->data; 2730 ierr = MatAXPY(y->A,a,x->A,str);CHKERRQ(ierr); 2731 PetscFunctionReturn(0); 2732 } 2733 2734 static PetscErrorCode MatGetLocalSubMatrix_IS(Mat A,IS row,IS col,Mat *submat) 2735 { 2736 Mat lA; 2737 Mat_IS *matis; 2738 ISLocalToGlobalMapping rl2g,cl2g; 2739 IS is; 2740 const PetscInt *rg,*rl; 2741 PetscInt nrg; 2742 PetscInt N,M,nrl,i,*idxs; 2743 PetscErrorCode ierr; 2744 2745 PetscFunctionBegin; 2746 ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 2747 ierr = ISGetLocalSize(row,&nrl);CHKERRQ(ierr); 2748 ierr = ISGetIndices(row,&rl);CHKERRQ(ierr); 2749 ierr = ISLocalToGlobalMappingGetSize(A->rmap->mapping,&nrg);CHKERRQ(ierr); 2750 #if defined(PETSC_USE_DEBUG) 2751 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); 2752 #endif 2753 ierr = PetscMalloc1(nrg,&idxs);CHKERRQ(ierr); 2754 /* map from [0,nrl) to row */ 2755 for (i=0;i<nrl;i++) idxs[i] = rl[i]; 2756 #if defined(PETSC_USE_DEBUG) 2757 for (i=nrl;i<nrg;i++) idxs[i] = nrg; 2758 #else 2759 for (i=nrl;i<nrg;i++) idxs[i] = -1; 2760 #endif 2761 ierr = ISRestoreIndices(row,&rl);CHKERRQ(ierr); 2762 ierr = ISLocalToGlobalMappingRestoreIndices(A->rmap->mapping,&rg);CHKERRQ(ierr); 2763 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),nrg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 2764 ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g);CHKERRQ(ierr); 2765 ierr = ISDestroy(&is);CHKERRQ(ierr); 2766 /* compute new l2g map for columns */ 2767 if (col != row || A->rmap->mapping != A->cmap->mapping) { 2768 const PetscInt *cg,*cl; 2769 PetscInt ncg; 2770 PetscInt ncl; 2771 2772 ierr = ISLocalToGlobalMappingGetIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 2773 ierr = ISGetLocalSize(col,&ncl);CHKERRQ(ierr); 2774 ierr = ISGetIndices(col,&cl);CHKERRQ(ierr); 2775 ierr = ISLocalToGlobalMappingGetSize(A->cmap->mapping,&ncg);CHKERRQ(ierr); 2776 #if defined(PETSC_USE_DEBUG) 2777 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); 2778 #endif 2779 ierr = PetscMalloc1(ncg,&idxs);CHKERRQ(ierr); 2780 /* map from [0,ncl) to col */ 2781 for (i=0;i<ncl;i++) idxs[i] = cl[i]; 2782 #if defined(PETSC_USE_DEBUG) 2783 for (i=ncl;i<ncg;i++) idxs[i] = ncg; 2784 #else 2785 for (i=ncl;i<ncg;i++) idxs[i] = -1; 2786 #endif 2787 ierr = ISRestoreIndices(col,&cl);CHKERRQ(ierr); 2788 ierr = ISLocalToGlobalMappingRestoreIndices(A->cmap->mapping,&cg);CHKERRQ(ierr); 2789 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)A),ncg,idxs,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 2790 ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g);CHKERRQ(ierr); 2791 ierr = ISDestroy(&is);CHKERRQ(ierr); 2792 } else { 2793 ierr = PetscObjectReference((PetscObject)rl2g);CHKERRQ(ierr); 2794 cl2g = rl2g; 2795 } 2796 /* create the MATIS submatrix */ 2797 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 2798 ierr = MatCreate(PetscObjectComm((PetscObject)A),submat);CHKERRQ(ierr); 2799 ierr = MatSetSizes(*submat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2800 ierr = MatSetType(*submat,MATIS);CHKERRQ(ierr); 2801 matis = (Mat_IS*)((*submat)->data); 2802 matis->islocalref = PETSC_TRUE; 2803 ierr = MatSetLocalToGlobalMapping(*submat,rl2g,cl2g);CHKERRQ(ierr); 2804 ierr = MatISGetLocalMat(A,&lA);CHKERRQ(ierr); 2805 ierr = MatISSetLocalMat(*submat,lA);CHKERRQ(ierr); 2806 ierr = MatSetUp(*submat);CHKERRQ(ierr); 2807 ierr = MatAssemblyBegin(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2808 ierr = MatAssemblyEnd(*submat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2809 ierr = ISLocalToGlobalMappingDestroy(&rl2g);CHKERRQ(ierr); 2810 ierr = ISLocalToGlobalMappingDestroy(&cl2g);CHKERRQ(ierr); 2811 /* remove unsupported ops */ 2812 ierr = PetscMemzero((*submat)->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2813 (*submat)->ops->destroy = MatDestroy_IS; 2814 (*submat)->ops->setvalueslocal = MatSetValuesLocal_SubMat_IS; 2815 (*submat)->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_SubMat_IS; 2816 (*submat)->ops->assemblybegin = MatAssemblyBegin_IS; 2817 (*submat)->ops->assemblyend = MatAssemblyEnd_IS; 2818 PetscFunctionReturn(0); 2819 } 2820 2821 static PetscErrorCode MatSetFromOptions_IS(PetscOptionItems *PetscOptionsObject, Mat A) 2822 { 2823 Mat_IS *a = (Mat_IS*)A->data; 2824 PetscErrorCode ierr; 2825 2826 PetscFunctionBegin; 2827 ierr = PetscOptionsHead(PetscOptionsObject,"MATIS options");CHKERRQ(ierr); 2828 ierr = PetscObjectOptionsBegin((PetscObject)A); 2829 ierr = PetscOptionsBool("-matis_fixempty","Fix local matrices in case of empty local rows/columns",NULL,a->locempty,&a->locempty,NULL);CHKERRQ(ierr); 2830 ierr = PetscOptionsBool("-matis_storel2l","Store local-to-local matrices generated from PtAP operations","MatISStoreL2L",a->storel2l,&a->storel2l,NULL);CHKERRQ(ierr); 2831 ierr = PetscOptionsEnd();CHKERRQ(ierr); 2832 PetscFunctionReturn(0); 2833 } 2834 2835 /*@ 2836 MatCreateIS - Creates a "process" unassembled matrix, assembled on each 2837 process but not across processes. 2838 2839 Input Parameters: 2840 + comm - MPI communicator that will share the matrix 2841 . bs - block size of the matrix 2842 . m,n,M,N - local and/or global sizes of the left and right vector used in matrix vector products 2843 . rmap - local to global map for rows 2844 - cmap - local to global map for cols 2845 2846 Output Parameter: 2847 . A - the resulting matrix 2848 2849 Level: advanced 2850 2851 Notes: 2852 See MATIS for more details. 2853 m and n are NOT related to the size of the map; they represent the size of the local parts of the vectors 2854 used in MatMult operations. The sizes of rmap and cmap define the size of the local matrices. 2855 If either rmap or cmap are NULL, then the matrix is assumed to be square. 2856 2857 .seealso: MATIS, MatSetLocalToGlobalMapping() 2858 @*/ 2859 PetscErrorCode MatCreateIS(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt M,PetscInt N,ISLocalToGlobalMapping rmap,ISLocalToGlobalMapping cmap,Mat *A) 2860 { 2861 PetscErrorCode ierr; 2862 2863 PetscFunctionBegin; 2864 if (!rmap && !cmap) SETERRQ(comm,PETSC_ERR_USER,"You need to provide at least one of the mappings"); 2865 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2866 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 2867 if (bs > 0) { 2868 ierr = MatSetBlockSize(*A,bs);CHKERRQ(ierr); 2869 } 2870 ierr = MatSetType(*A,MATIS);CHKERRQ(ierr); 2871 if (rmap && cmap) { 2872 ierr = MatSetLocalToGlobalMapping(*A,rmap,cmap);CHKERRQ(ierr); 2873 } else if (!rmap) { 2874 ierr = MatSetLocalToGlobalMapping(*A,cmap,cmap);CHKERRQ(ierr); 2875 } else { 2876 ierr = MatSetLocalToGlobalMapping(*A,rmap,rmap);CHKERRQ(ierr); 2877 } 2878 PetscFunctionReturn(0); 2879 } 2880 2881 /*MC 2882 MATIS - MATIS = "is" - A matrix type to be used for using the non-overlapping domain decomposition methods (e.g. PCBDDC or KSPFETIDP). 2883 This stores the matrices in globally unassembled form. Each processor 2884 assembles only its local Neumann problem and the parallel matrix vector 2885 product is handled "implicitly". 2886 2887 Operations Provided: 2888 + MatMult() 2889 . MatMultAdd() 2890 . MatMultTranspose() 2891 . MatMultTransposeAdd() 2892 . MatZeroEntries() 2893 . MatSetOption() 2894 . MatZeroRows() 2895 . MatSetValues() 2896 . MatSetValuesBlocked() 2897 . MatSetValuesLocal() 2898 . MatSetValuesBlockedLocal() 2899 . MatScale() 2900 . MatGetDiagonal() 2901 . MatMissingDiagonal() 2902 . MatDuplicate() 2903 . MatCopy() 2904 . MatAXPY() 2905 . MatCreateSubMatrix() 2906 . MatGetLocalSubMatrix() 2907 . MatTranspose() 2908 . MatPtAP() (with P of AIJ type) 2909 - MatSetLocalToGlobalMapping() 2910 2911 Options Database Keys: 2912 + -mat_type is - sets the matrix type to "is" during a call to MatSetFromOptions() 2913 . -matis_fixempty - Fixes local matrices in case of empty local rows/columns. 2914 - -matis_storel2l - stores the local-to-local operators generated by the Galerkin process of MatPtAP(). 2915 2916 Notes: 2917 Options prefix for the inner matrix are given by -is_mat_xxx 2918 2919 You must call MatSetLocalToGlobalMapping() before using this matrix type. 2920 2921 You can do matrix preallocation on the local matrix after you obtain it with 2922 MatISGetLocalMat(); otherwise, you could use MatISSetPreallocation() 2923 2924 Level: advanced 2925 2926 .seealso: Mat, MatISGetLocalMat(), MatSetLocalToGlobalMapping(), MatISSetPreallocation(), MatCreateIS(), PCBDDC, KSPFETIDP 2927 2928 M*/ 2929 2930 PETSC_EXTERN PetscErrorCode MatCreate_IS(Mat A) 2931 { 2932 PetscErrorCode ierr; 2933 Mat_IS *b; 2934 2935 PetscFunctionBegin; 2936 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 2937 A->data = (void*)b; 2938 2939 /* matrix ops */ 2940 ierr = PetscMemzero(A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2941 A->ops->mult = MatMult_IS; 2942 A->ops->multadd = MatMultAdd_IS; 2943 A->ops->multtranspose = MatMultTranspose_IS; 2944 A->ops->multtransposeadd = MatMultTransposeAdd_IS; 2945 A->ops->destroy = MatDestroy_IS; 2946 A->ops->setlocaltoglobalmapping = MatSetLocalToGlobalMapping_IS; 2947 A->ops->setvalues = MatSetValues_IS; 2948 A->ops->setvaluesblocked = MatSetValuesBlocked_IS; 2949 A->ops->setvalueslocal = MatSetValuesLocal_IS; 2950 A->ops->setvaluesblockedlocal = MatSetValuesBlockedLocal_IS; 2951 A->ops->zerorows = MatZeroRows_IS; 2952 A->ops->zerorowscolumns = MatZeroRowsColumns_IS; 2953 A->ops->assemblybegin = MatAssemblyBegin_IS; 2954 A->ops->assemblyend = MatAssemblyEnd_IS; 2955 A->ops->view = MatView_IS; 2956 A->ops->zeroentries = MatZeroEntries_IS; 2957 A->ops->scale = MatScale_IS; 2958 A->ops->getdiagonal = MatGetDiagonal_IS; 2959 A->ops->setoption = MatSetOption_IS; 2960 A->ops->ishermitian = MatIsHermitian_IS; 2961 A->ops->issymmetric = MatIsSymmetric_IS; 2962 A->ops->isstructurallysymmetric = MatIsStructurallySymmetric_IS; 2963 A->ops->duplicate = MatDuplicate_IS; 2964 A->ops->missingdiagonal = MatMissingDiagonal_IS; 2965 A->ops->copy = MatCopy_IS; 2966 A->ops->getlocalsubmatrix = MatGetLocalSubMatrix_IS; 2967 A->ops->createsubmatrix = MatCreateSubMatrix_IS; 2968 A->ops->axpy = MatAXPY_IS; 2969 A->ops->diagonalset = MatDiagonalSet_IS; 2970 A->ops->shift = MatShift_IS; 2971 A->ops->transpose = MatTranspose_IS; 2972 A->ops->getinfo = MatGetInfo_IS; 2973 A->ops->diagonalscale = MatDiagonalScale_IS; 2974 A->ops->setfromoptions = MatSetFromOptions_IS; 2975 2976 /* special MATIS functions */ 2977 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetLocalMat_C",MatISGetLocalMat_IS);CHKERRQ(ierr); 2978 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISRestoreLocalMat_C",MatISRestoreLocalMat_IS);CHKERRQ(ierr); 2979 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetLocalMat_C",MatISSetLocalMat_IS);CHKERRQ(ierr); 2980 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISGetMPIXAIJ_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 2981 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetPreallocation_C",MatISSetPreallocation_IS);CHKERRQ(ierr); 2982 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISSetUpSF_C",MatISSetUpSF_IS);CHKERRQ(ierr); 2983 ierr = PetscObjectComposeFunction((PetscObject)A,"MatISStoreL2L_C",MatISStoreL2L_IS);CHKERRQ(ierr); 2984 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpiaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 2985 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpibaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 2986 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_mpisbaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 2987 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 2988 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqbaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 2989 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_seqsbaij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 2990 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_is_aij_C",MatConvert_IS_XAIJ);CHKERRQ(ierr); 2991 ierr = PetscObjectChangeTypeName((PetscObject)A,MATIS);CHKERRQ(ierr); 2992 PetscFunctionReturn(0); 2993 } 2994