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