1 #include <../src/mat/impls/aij/mpi/mpiaij.h> /*I "petscmat.h" I*/ 2 #include <../src/mat/impls/sell/mpi/mpisell.h> /*I "petscmat.h" I*/ 3 #include <petsc/private/vecimpl.h> 4 #include <petsc/private/isimpl.h> 5 #include <petscblaslapack.h> 6 #include <petscsf.h> 7 8 /*MC 9 MATSELL - MATSELL = "sell" - A matrix type to be used for sparse matrices. 10 11 This matrix type is identical to MATSEQSELL when constructed with a single process communicator, 12 and MATMPISELL otherwise. As a result, for single process communicators, 13 MatSeqSELLSetPreallocation is supported, and similarly MatMPISELLSetPreallocation is supported 14 for communicators controlling multiple processes. It is recommended that you call both of 15 the above preallocation routines for simplicity. 16 17 Options Database Keys: 18 . -mat_type sell - sets the matrix type to "sell" during a call to MatSetFromOptions() 19 20 Level: beginner 21 22 .seealso: MatCreateSELL(), MatCreateSeqSELL(), MATSEQSELL, MATMPISELL 23 M*/ 24 25 PetscErrorCode MatDiagonalSet_MPISELL(Mat Y,Vec D,InsertMode is) 26 { 27 Mat_MPISELL *sell=(Mat_MPISELL*)Y->data; 28 29 PetscFunctionBegin; 30 if (Y->assembled && Y->rmap->rstart == Y->cmap->rstart && Y->rmap->rend == Y->cmap->rend) { 31 PetscCall(MatDiagonalSet(sell->A,D,is)); 32 } else { 33 PetscCall(MatDiagonalSet_Default(Y,D,is)); 34 } 35 PetscFunctionReturn(0); 36 } 37 38 /* 39 Local utility routine that creates a mapping from the global column 40 number to the local number in the off-diagonal part of the local 41 storage of the matrix. When PETSC_USE_CTABLE is used this is scalable at 42 a slightly higher hash table cost; without it it is not scalable (each processor 43 has an order N integer array but is fast to acess. 44 */ 45 PetscErrorCode MatCreateColmap_MPISELL_Private(Mat mat) 46 { 47 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 48 PetscInt n=sell->B->cmap->n,i; 49 50 PetscFunctionBegin; 51 PetscCheck(sell->garray,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MPISELL Matrix was assembled but is missing garray"); 52 #if defined(PETSC_USE_CTABLE) 53 PetscCall(PetscTableCreate(n,mat->cmap->N+1,&sell->colmap)); 54 for (i=0; i<n; i++) { 55 PetscCall(PetscTableAdd(sell->colmap,sell->garray[i]+1,i+1,INSERT_VALUES)); 56 } 57 #else 58 PetscCall(PetscCalloc1(mat->cmap->N+1,&sell->colmap)); 59 PetscCall(PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N+1)*sizeof(PetscInt))); 60 for (i=0; i<n; i++) sell->colmap[sell->garray[i]] = i+1; 61 #endif 62 PetscFunctionReturn(0); 63 } 64 65 #define MatSetValues_SeqSELL_A_Private(row,col,value,addv,orow,ocol) \ 66 { \ 67 if (col <= lastcol1) low1 = 0; \ 68 else high1 = nrow1; \ 69 lastcol1 = col; \ 70 while (high1-low1 > 5) { \ 71 t = (low1+high1)/2; \ 72 if (*(cp1+8*t) > col) high1 = t; \ 73 else low1 = t; \ 74 } \ 75 for (_i=low1; _i<high1; _i++) { \ 76 if (*(cp1+8*_i) > col) break; \ 77 if (*(cp1+8*_i) == col) { \ 78 if (addv == ADD_VALUES) *(vp1+8*_i) += value; \ 79 else *(vp1+8*_i) = value; \ 80 goto a_noinsert; \ 81 } \ 82 } \ 83 if (value == 0.0 && ignorezeroentries) {low1 = 0; high1 = nrow1;goto a_noinsert;} \ 84 if (nonew == 1) {low1 = 0; high1 = nrow1; goto a_noinsert;} \ 85 PetscCheckFalse(nonew == -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", orow, ocol); \ 86 MatSeqXSELLReallocateSELL(A,am,1,nrow1,a->sliidx,row/8,row,col,a->colidx,a->val,cp1,vp1,nonew,MatScalar); \ 87 /* shift up all the later entries in this row */ \ 88 for (ii=nrow1-1; ii>=_i; ii--) { \ 89 *(cp1+8*(ii+1)) = *(cp1+8*ii); \ 90 *(vp1+8*(ii+1)) = *(vp1+8*ii); \ 91 } \ 92 *(cp1+8*_i) = col; \ 93 *(vp1+8*_i) = value; \ 94 a->nz++; nrow1++; A->nonzerostate++; \ 95 a_noinsert: ; \ 96 a->rlen[row] = nrow1; \ 97 } 98 99 #define MatSetValues_SeqSELL_B_Private(row,col,value,addv,orow,ocol) \ 100 { \ 101 if (col <= lastcol2) low2 = 0; \ 102 else high2 = nrow2; \ 103 lastcol2 = col; \ 104 while (high2-low2 > 5) { \ 105 t = (low2+high2)/2; \ 106 if (*(cp2+8*t) > col) high2 = t; \ 107 else low2 = t; \ 108 } \ 109 for (_i=low2; _i<high2; _i++) { \ 110 if (*(cp2+8*_i) > col) break; \ 111 if (*(cp2+8*_i) == col) { \ 112 if (addv == ADD_VALUES) *(vp2+8*_i) += value; \ 113 else *(vp2+8*_i) = value; \ 114 goto b_noinsert; \ 115 } \ 116 } \ 117 if (value == 0.0 && ignorezeroentries) {low2 = 0; high2 = nrow2; goto b_noinsert;} \ 118 if (nonew == 1) {low2 = 0; high2 = nrow2; goto b_noinsert;} \ 119 PetscCheckFalse(nonew == -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", orow, ocol); \ 120 MatSeqXSELLReallocateSELL(B,bm,1,nrow2,b->sliidx,row/8,row,col,b->colidx,b->val,cp2,vp2,nonew,MatScalar); \ 121 /* shift up all the later entries in this row */ \ 122 for (ii=nrow2-1; ii>=_i; ii--) { \ 123 *(cp2+8*(ii+1)) = *(cp2+8*ii); \ 124 *(vp2+8*(ii+1)) = *(vp2+8*ii); \ 125 } \ 126 *(cp2+8*_i) = col; \ 127 *(vp2+8*_i) = value; \ 128 b->nz++; nrow2++; B->nonzerostate++; \ 129 b_noinsert: ; \ 130 b->rlen[row] = nrow2; \ 131 } 132 133 PetscErrorCode MatSetValues_MPISELL(Mat mat,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode addv) 134 { 135 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 136 PetscScalar value; 137 PetscInt i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend,shift1,shift2; 138 PetscInt cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col; 139 PetscBool roworiented=sell->roworiented; 140 141 /* Some Variables required in the macro */ 142 Mat A=sell->A; 143 Mat_SeqSELL *a=(Mat_SeqSELL*)A->data; 144 PetscBool ignorezeroentries=a->ignorezeroentries,found; 145 Mat B=sell->B; 146 Mat_SeqSELL *b=(Mat_SeqSELL*)B->data; 147 PetscInt *cp1,*cp2,ii,_i,nrow1,nrow2,low1,high1,low2,high2,t,lastcol1,lastcol2; 148 MatScalar *vp1,*vp2; 149 150 PetscFunctionBegin; 151 for (i=0; i<m; i++) { 152 if (im[i] < 0) continue; 153 PetscCheck(im[i] < mat->rmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,im[i],mat->rmap->N-1); 154 if (im[i] >= rstart && im[i] < rend) { 155 row = im[i] - rstart; 156 lastcol1 = -1; 157 shift1 = a->sliidx[row>>3]+(row&0x07); /* starting index of the row */ 158 cp1 = a->colidx+shift1; 159 vp1 = a->val+shift1; 160 nrow1 = a->rlen[row]; 161 low1 = 0; 162 high1 = nrow1; 163 lastcol2 = -1; 164 shift2 = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */ 165 cp2 = b->colidx+shift2; 166 vp2 = b->val+shift2; 167 nrow2 = b->rlen[row]; 168 low2 = 0; 169 high2 = nrow2; 170 171 for (j=0; j<n; j++) { 172 if (roworiented) value = v[i*n+j]; 173 else value = v[i+j*m]; 174 if (ignorezeroentries && value == 0.0 && (addv == ADD_VALUES)) continue; 175 if (in[j] >= cstart && in[j] < cend) { 176 col = in[j] - cstart; 177 MatSetValue_SeqSELL_Private(A,row,col,value,addv,im[i],in[j],cp1,vp1,lastcol1,low1,high1); /* set one value */ 178 } else if (in[j] < 0) continue; 179 else PetscCheck(in[j] < mat->cmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[j],mat->cmap->N-1); 180 else { 181 if (mat->was_assembled) { 182 if (!sell->colmap) { 183 PetscCall(MatCreateColmap_MPISELL_Private(mat)); 184 } 185 #if defined(PETSC_USE_CTABLE) 186 PetscCall(PetscTableFind(sell->colmap,in[j]+1,&col)); 187 col--; 188 #else 189 col = sell->colmap[in[j]] - 1; 190 #endif 191 if (col < 0 && !((Mat_SeqSELL*)(sell->B->data))->nonew) { 192 PetscCall(MatDisAssemble_MPISELL(mat)); 193 col = in[j]; 194 /* Reinitialize the variables required by MatSetValues_SeqSELL_B_Private() */ 195 B = sell->B; 196 b = (Mat_SeqSELL*)B->data; 197 shift2 = b->sliidx[row>>3]+(row&0x07); /* starting index of the row */ 198 cp2 = b->colidx+shift2; 199 vp2 = b->val+shift2; 200 nrow2 = b->rlen[row]; 201 low2 = 0; 202 high2 = nrow2; 203 } else PetscCheckFalse(col < 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", im[i], in[j]); 204 } else col = in[j]; 205 MatSetValue_SeqSELL_Private(B,row,col,value,addv,im[i],in[j],cp2,vp2,lastcol2,low2,high2); /* set one value */ 206 } 207 } 208 } else { 209 PetscCheck(!mat->nooffprocentries,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set",im[i]); 210 if (!sell->donotstash) { 211 mat->assembled = PETSC_FALSE; 212 if (roworiented) { 213 PetscCall(MatStashValuesRow_Private(&mat->stash,im[i],n,in,v+i*n,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)))); 214 } else { 215 PetscCall(MatStashValuesCol_Private(&mat->stash,im[i],n,in,v+i,m,(PetscBool)(ignorezeroentries && (addv == ADD_VALUES)))); 216 } 217 } 218 } 219 } 220 PetscFunctionReturn(0); 221 } 222 223 PetscErrorCode MatGetValues_MPISELL(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[]) 224 { 225 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 226 PetscInt i,j,rstart=mat->rmap->rstart,rend=mat->rmap->rend; 227 PetscInt cstart=mat->cmap->rstart,cend=mat->cmap->rend,row,col; 228 229 PetscFunctionBegin; 230 for (i=0; i<m; i++) { 231 if (idxm[i] < 0) continue; /* negative row */ 232 PetscCheck(idxm[i] < mat->rmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,idxm[i],mat->rmap->N-1); 233 if (idxm[i] >= rstart && idxm[i] < rend) { 234 row = idxm[i] - rstart; 235 for (j=0; j<n; j++) { 236 if (idxn[j] < 0) continue; /* negative column */ 237 PetscCheck(idxn[j] < mat->cmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,idxn[j],mat->cmap->N-1); 238 if (idxn[j] >= cstart && idxn[j] < cend) { 239 col = idxn[j] - cstart; 240 PetscCall(MatGetValues(sell->A,1,&row,1,&col,v+i*n+j)); 241 } else { 242 if (!sell->colmap) { 243 PetscCall(MatCreateColmap_MPISELL_Private(mat)); 244 } 245 #if defined(PETSC_USE_CTABLE) 246 PetscCall(PetscTableFind(sell->colmap,idxn[j]+1,&col)); 247 col--; 248 #else 249 col = sell->colmap[idxn[j]] - 1; 250 #endif 251 if ((col < 0) || (sell->garray[col] != idxn[j])) *(v+i*n+j) = 0.0; 252 else { 253 PetscCall(MatGetValues(sell->B,1,&row,1,&col,v+i*n+j)); 254 } 255 } 256 } 257 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only local values currently supported"); 258 } 259 PetscFunctionReturn(0); 260 } 261 262 extern PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat,Vec,Vec); 263 264 PetscErrorCode MatAssemblyBegin_MPISELL(Mat mat,MatAssemblyType mode) 265 { 266 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 267 PetscInt nstash,reallocs; 268 269 PetscFunctionBegin; 270 if (sell->donotstash || mat->nooffprocentries) PetscFunctionReturn(0); 271 272 PetscCall(MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range)); 273 PetscCall(MatStashGetInfo_Private(&mat->stash,&nstash,&reallocs)); 274 PetscCall(PetscInfo(sell->A,"Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n",nstash,reallocs)); 275 PetscFunctionReturn(0); 276 } 277 278 PetscErrorCode MatAssemblyEnd_MPISELL(Mat mat,MatAssemblyType mode) 279 { 280 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 281 PetscMPIInt n; 282 PetscInt i,flg; 283 PetscInt *row,*col; 284 PetscScalar *val; 285 PetscBool other_disassembled; 286 /* do not use 'b = (Mat_SeqSELL*)sell->B->data' as B can be reset in disassembly */ 287 PetscFunctionBegin; 288 if (!sell->donotstash && !mat->nooffprocentries) { 289 while (1) { 290 PetscCall(MatStashScatterGetMesg_Private(&mat->stash,&n,&row,&col,&val,&flg)); 291 if (!flg) break; 292 293 for (i=0; i<n; i++) { /* assemble one by one */ 294 PetscCall(MatSetValues_MPISELL(mat,1,row+i,1,col+i,val+i,mat->insertmode)); 295 } 296 } 297 PetscCall(MatStashScatterEnd_Private(&mat->stash)); 298 } 299 PetscCall(MatAssemblyBegin(sell->A,mode)); 300 PetscCall(MatAssemblyEnd(sell->A,mode)); 301 302 /* 303 determine if any processor has disassembled, if so we must 304 also disassemble ourselfs, in order that we may reassemble. 305 */ 306 /* 307 if nonzero structure of submatrix B cannot change then we know that 308 no processor disassembled thus we can skip this stuff 309 */ 310 if (!((Mat_SeqSELL*)sell->B->data)->nonew) { 311 PetscCall(MPIU_Allreduce(&mat->was_assembled,&other_disassembled,1,MPIU_BOOL,MPI_PROD,PetscObjectComm((PetscObject)mat))); 312 PetscCheckFalse(mat->was_assembled && !other_disassembled,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatDisAssemble not implemented yet"); 313 } 314 if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) { 315 PetscCall(MatSetUpMultiply_MPISELL(mat)); 316 } 317 /* 318 PetscCall(MatSetOption(sell->B,MAT_USE_INODES,PETSC_FALSE)); 319 */ 320 PetscCall(MatAssemblyBegin(sell->B,mode)); 321 PetscCall(MatAssemblyEnd(sell->B,mode)); 322 PetscCall(PetscFree2(sell->rowvalues,sell->rowindices)); 323 sell->rowvalues = NULL; 324 PetscCall(VecDestroy(&sell->diag)); 325 326 /* if no new nonzero locations are allowed in matrix then only set the matrix state the first time through */ 327 if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqSELL*)(sell->A->data))->nonew) { 328 PetscObjectState state = sell->A->nonzerostate + sell->B->nonzerostate; 329 PetscCall(MPIU_Allreduce(&state,&mat->nonzerostate,1,MPIU_INT64,MPI_SUM,PetscObjectComm((PetscObject)mat))); 330 } 331 PetscFunctionReturn(0); 332 } 333 334 PetscErrorCode MatZeroEntries_MPISELL(Mat A) 335 { 336 Mat_MPISELL *l=(Mat_MPISELL*)A->data; 337 338 PetscFunctionBegin; 339 PetscCall(MatZeroEntries(l->A)); 340 PetscCall(MatZeroEntries(l->B)); 341 PetscFunctionReturn(0); 342 } 343 344 PetscErrorCode MatMult_MPISELL(Mat A,Vec xx,Vec yy) 345 { 346 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 347 PetscInt nt; 348 349 PetscFunctionBegin; 350 PetscCall(VecGetLocalSize(xx,&nt)); 351 PetscCheckFalse(nt != A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible partition of A (%" PetscInt_FMT ") and xx (%" PetscInt_FMT ")",A->cmap->n,nt); 352 PetscCall(VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD)); 353 PetscCall((*a->A->ops->mult)(a->A,xx,yy)); 354 PetscCall(VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD)); 355 PetscCall((*a->B->ops->multadd)(a->B,a->lvec,yy,yy)); 356 PetscFunctionReturn(0); 357 } 358 359 PetscErrorCode MatMultDiagonalBlock_MPISELL(Mat A,Vec bb,Vec xx) 360 { 361 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 362 363 PetscFunctionBegin; 364 PetscCall(MatMultDiagonalBlock(a->A,bb,xx)); 365 PetscFunctionReturn(0); 366 } 367 368 PetscErrorCode MatMultAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz) 369 { 370 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 371 372 PetscFunctionBegin; 373 PetscCall(VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD)); 374 PetscCall((*a->A->ops->multadd)(a->A,xx,yy,zz)); 375 PetscCall(VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD)); 376 PetscCall((*a->B->ops->multadd)(a->B,a->lvec,zz,zz)); 377 PetscFunctionReturn(0); 378 } 379 380 PetscErrorCode MatMultTranspose_MPISELL(Mat A,Vec xx,Vec yy) 381 { 382 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 383 384 PetscFunctionBegin; 385 /* do nondiagonal part */ 386 PetscCall((*a->B->ops->multtranspose)(a->B,xx,a->lvec)); 387 /* do local part */ 388 PetscCall((*a->A->ops->multtranspose)(a->A,xx,yy)); 389 /* add partial results together */ 390 PetscCall(VecScatterBegin(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE)); 391 PetscCall(VecScatterEnd(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE)); 392 PetscFunctionReturn(0); 393 } 394 395 PetscErrorCode MatIsTranspose_MPISELL(Mat Amat,Mat Bmat,PetscReal tol,PetscBool *f) 396 { 397 MPI_Comm comm; 398 Mat_MPISELL *Asell=(Mat_MPISELL*)Amat->data,*Bsell; 399 Mat Adia=Asell->A,Bdia,Aoff,Boff,*Aoffs,*Boffs; 400 IS Me,Notme; 401 PetscInt M,N,first,last,*notme,i; 402 PetscMPIInt size; 403 404 PetscFunctionBegin; 405 /* Easy test: symmetric diagonal block */ 406 Bsell = (Mat_MPISELL*)Bmat->data; Bdia = Bsell->A; 407 PetscCall(MatIsTranspose(Adia,Bdia,tol,f)); 408 if (!*f) PetscFunctionReturn(0); 409 PetscCall(PetscObjectGetComm((PetscObject)Amat,&comm)); 410 PetscCallMPI(MPI_Comm_size(comm,&size)); 411 if (size == 1) PetscFunctionReturn(0); 412 413 /* Hard test: off-diagonal block. This takes a MatCreateSubMatrix. */ 414 PetscCall(MatGetSize(Amat,&M,&N)); 415 PetscCall(MatGetOwnershipRange(Amat,&first,&last)); 416 PetscCall(PetscMalloc1(N-last+first,¬me)); 417 for (i=0; i<first; i++) notme[i] = i; 418 for (i=last; i<M; i++) notme[i-last+first] = i; 419 PetscCall(ISCreateGeneral(MPI_COMM_SELF,N-last+first,notme,PETSC_COPY_VALUES,&Notme)); 420 PetscCall(ISCreateStride(MPI_COMM_SELF,last-first,first,1,&Me)); 421 PetscCall(MatCreateSubMatrices(Amat,1,&Me,&Notme,MAT_INITIAL_MATRIX,&Aoffs)); 422 Aoff = Aoffs[0]; 423 PetscCall(MatCreateSubMatrices(Bmat,1,&Notme,&Me,MAT_INITIAL_MATRIX,&Boffs)); 424 Boff = Boffs[0]; 425 PetscCall(MatIsTranspose(Aoff,Boff,tol,f)); 426 PetscCall(MatDestroyMatrices(1,&Aoffs)); 427 PetscCall(MatDestroyMatrices(1,&Boffs)); 428 PetscCall(ISDestroy(&Me)); 429 PetscCall(ISDestroy(&Notme)); 430 PetscCall(PetscFree(notme)); 431 PetscFunctionReturn(0); 432 } 433 434 PetscErrorCode MatMultTransposeAdd_MPISELL(Mat A,Vec xx,Vec yy,Vec zz) 435 { 436 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 437 438 PetscFunctionBegin; 439 /* do nondiagonal part */ 440 PetscCall((*a->B->ops->multtranspose)(a->B,xx,a->lvec)); 441 /* do local part */ 442 PetscCall((*a->A->ops->multtransposeadd)(a->A,xx,yy,zz)); 443 /* add partial results together */ 444 PetscCall(VecScatterBegin(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE)); 445 PetscCall(VecScatterEnd(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE)); 446 PetscFunctionReturn(0); 447 } 448 449 /* 450 This only works correctly for square matrices where the subblock A->A is the 451 diagonal block 452 */ 453 PetscErrorCode MatGetDiagonal_MPISELL(Mat A,Vec v) 454 { 455 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 456 457 PetscFunctionBegin; 458 PetscCheckFalse(A->rmap->N != A->cmap->N,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Supports only square matrix where A->A is diag block"); 459 PetscCheckFalse(A->rmap->rstart != A->cmap->rstart || A->rmap->rend != A->cmap->rend,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"row partition must equal col partition"); 460 PetscCall(MatGetDiagonal(a->A,v)); 461 PetscFunctionReturn(0); 462 } 463 464 PetscErrorCode MatScale_MPISELL(Mat A,PetscScalar aa) 465 { 466 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 467 468 PetscFunctionBegin; 469 PetscCall(MatScale(a->A,aa)); 470 PetscCall(MatScale(a->B,aa)); 471 PetscFunctionReturn(0); 472 } 473 474 PetscErrorCode MatDestroy_MPISELL(Mat mat) 475 { 476 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 477 478 PetscFunctionBegin; 479 #if defined(PETSC_USE_LOG) 480 PetscLogObjectState((PetscObject)mat,"Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT,mat->rmap->N,mat->cmap->N); 481 #endif 482 PetscCall(MatStashDestroy_Private(&mat->stash)); 483 PetscCall(VecDestroy(&sell->diag)); 484 PetscCall(MatDestroy(&sell->A)); 485 PetscCall(MatDestroy(&sell->B)); 486 #if defined(PETSC_USE_CTABLE) 487 PetscCall(PetscTableDestroy(&sell->colmap)); 488 #else 489 PetscCall(PetscFree(sell->colmap)); 490 #endif 491 PetscCall(PetscFree(sell->garray)); 492 PetscCall(VecDestroy(&sell->lvec)); 493 PetscCall(VecScatterDestroy(&sell->Mvctx)); 494 PetscCall(PetscFree2(sell->rowvalues,sell->rowindices)); 495 PetscCall(PetscFree(sell->ld)); 496 PetscCall(PetscFree(mat->data)); 497 498 PetscCall(PetscObjectChangeTypeName((PetscObject)mat,NULL)); 499 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatStoreValues_C",NULL)); 500 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatRetrieveValues_C",NULL)); 501 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatIsTranspose_C",NULL)); 502 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatMPISELLSetPreallocation_C",NULL)); 503 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatConvert_mpisell_mpiaij_C",NULL)); 504 PetscCall(PetscObjectComposeFunction((PetscObject)mat,"MatDiagonalScaleLocal_C",NULL)); 505 PetscFunctionReturn(0); 506 } 507 508 #include <petscdraw.h> 509 PetscErrorCode MatView_MPISELL_ASCIIorDraworSocket(Mat mat,PetscViewer viewer) 510 { 511 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 512 PetscErrorCode ierr; 513 PetscMPIInt rank=sell->rank,size=sell->size; 514 PetscBool isdraw,iascii,isbinary; 515 PetscViewer sviewer; 516 PetscViewerFormat format; 517 518 PetscFunctionBegin; 519 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 520 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 521 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 522 if (iascii) { 523 PetscCall(PetscViewerGetFormat(viewer,&format)); 524 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 525 MatInfo info; 526 PetscInt *inodes; 527 528 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)mat),&rank)); 529 PetscCall(MatGetInfo(mat,MAT_LOCAL,&info)); 530 PetscCall(MatInodeGetInodeSizes(sell->A,NULL,&inodes,NULL)); 531 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 532 if (!inodes) { 533 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", not using I-node routines\n", 534 rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory);PetscCall(ierr); 535 } else { 536 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Local rows %" PetscInt_FMT " nz %" PetscInt_FMT " nz alloced %" PetscInt_FMT " mem %" PetscInt_FMT ", using I-node routines\n", 537 rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory);PetscCall(ierr); 538 } 539 PetscCall(MatGetInfo(sell->A,MAT_LOCAL,&info)); 540 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"[%d] on-diagonal part: nz %" PetscInt_FMT " \n",rank,(PetscInt)info.nz_used)); 541 PetscCall(MatGetInfo(sell->B,MAT_LOCAL,&info)); 542 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer,"[%d] off-diagonal part: nz %" PetscInt_FMT " \n",rank,(PetscInt)info.nz_used)); 543 PetscCall(PetscViewerFlush(viewer)); 544 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 545 PetscCall(PetscViewerASCIIPrintf(viewer,"Information on VecScatter used in matrix-vector product: \n")); 546 PetscCall(VecScatterView(sell->Mvctx,viewer)); 547 PetscFunctionReturn(0); 548 } else if (format == PETSC_VIEWER_ASCII_INFO) { 549 PetscInt inodecount,inodelimit,*inodes; 550 PetscCall(MatInodeGetInodeSizes(sell->A,&inodecount,&inodes,&inodelimit)); 551 if (inodes) { 552 PetscCall(PetscViewerASCIIPrintf(viewer,"using I-node (on process 0) routines: found %" PetscInt_FMT " nodes, limit used is %" PetscInt_FMT "\n",inodecount,inodelimit)); 553 } else { 554 PetscCall(PetscViewerASCIIPrintf(viewer,"not using I-node (on process 0) routines\n")); 555 } 556 PetscFunctionReturn(0); 557 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 558 PetscFunctionReturn(0); 559 } 560 } else if (isbinary) { 561 if (size == 1) { 562 PetscCall(PetscObjectSetName((PetscObject)sell->A,((PetscObject)mat)->name)); 563 PetscCall(MatView(sell->A,viewer)); 564 } else { 565 /* PetscCall(MatView_MPISELL_Binary(mat,viewer)); */ 566 } 567 PetscFunctionReturn(0); 568 } else if (isdraw) { 569 PetscDraw draw; 570 PetscBool isnull; 571 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 572 PetscCall(PetscDrawIsNull(draw,&isnull)); 573 if (isnull) PetscFunctionReturn(0); 574 } 575 576 { 577 /* assemble the entire matrix onto first processor. */ 578 Mat A; 579 Mat_SeqSELL *Aloc; 580 PetscInt M=mat->rmap->N,N=mat->cmap->N,*acolidx,row,col,i,j; 581 MatScalar *aval; 582 PetscBool isnonzero; 583 584 PetscCall(MatCreate(PetscObjectComm((PetscObject)mat),&A)); 585 if (rank == 0) { 586 PetscCall(MatSetSizes(A,M,N,M,N)); 587 } else { 588 PetscCall(MatSetSizes(A,0,0,M,N)); 589 } 590 /* This is just a temporary matrix, so explicitly using MATMPISELL is probably best */ 591 PetscCall(MatSetType(A,MATMPISELL)); 592 PetscCall(MatMPISELLSetPreallocation(A,0,NULL,0,NULL)); 593 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE)); 594 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)A)); 595 596 /* copy over the A part */ 597 Aloc = (Mat_SeqSELL*)sell->A->data; 598 acolidx = Aloc->colidx; aval = Aloc->val; 599 for (i=0; i<Aloc->totalslices; i++) { /* loop over slices */ 600 for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) { 601 isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]); 602 if (isnonzero) { /* check the mask bit */ 603 row = (i<<3)+(j&0x07) + mat->rmap->rstart; /* i<<3 is the starting row of this slice */ 604 col = *acolidx + mat->rmap->rstart; 605 PetscCall(MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES)); 606 } 607 aval++; acolidx++; 608 } 609 } 610 611 /* copy over the B part */ 612 Aloc = (Mat_SeqSELL*)sell->B->data; 613 acolidx = Aloc->colidx; aval = Aloc->val; 614 for (i=0; i<Aloc->totalslices; i++) { 615 for (j=Aloc->sliidx[i]; j<Aloc->sliidx[i+1]; j++) { 616 isnonzero = (PetscBool)((j-Aloc->sliidx[i])/8 < Aloc->rlen[(i<<3)+(j&0x07)]); 617 if (isnonzero) { 618 row = (i<<3)+(j&0x07) + mat->rmap->rstart; 619 col = sell->garray[*acolidx]; 620 PetscCall(MatSetValues(A,1,&row,1,&col,aval,INSERT_VALUES)); 621 } 622 aval++; acolidx++; 623 } 624 } 625 626 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 627 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 628 /* 629 Everyone has to call to draw the matrix since the graphics waits are 630 synchronized across all processors that share the PetscDraw object 631 */ 632 PetscCall(PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); 633 if (rank == 0) { 634 PetscCall(PetscObjectSetName((PetscObject)((Mat_MPISELL*)(A->data))->A,((PetscObject)mat)->name)); 635 PetscCall(MatView_SeqSELL(((Mat_MPISELL*)(A->data))->A,sviewer)); 636 } 637 PetscCall(PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer)); 638 PetscCall(PetscViewerFlush(viewer)); 639 PetscCall(MatDestroy(&A)); 640 } 641 PetscFunctionReturn(0); 642 } 643 644 PetscErrorCode MatView_MPISELL(Mat mat,PetscViewer viewer) 645 { 646 PetscBool iascii,isdraw,issocket,isbinary; 647 648 PetscFunctionBegin; 649 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 650 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 651 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 652 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket)); 653 if (iascii || isdraw || isbinary || issocket) { 654 PetscCall(MatView_MPISELL_ASCIIorDraworSocket(mat,viewer)); 655 } 656 PetscFunctionReturn(0); 657 } 658 659 PetscErrorCode MatGetGhosts_MPISELL(Mat mat,PetscInt *nghosts,const PetscInt *ghosts[]) 660 { 661 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 662 663 PetscFunctionBegin; 664 PetscCall(MatGetSize(sell->B,NULL,nghosts)); 665 if (ghosts) *ghosts = sell->garray; 666 PetscFunctionReturn(0); 667 } 668 669 PetscErrorCode MatGetInfo_MPISELL(Mat matin,MatInfoType flag,MatInfo *info) 670 { 671 Mat_MPISELL *mat=(Mat_MPISELL*)matin->data; 672 Mat A=mat->A,B=mat->B; 673 PetscLogDouble isend[5],irecv[5]; 674 675 PetscFunctionBegin; 676 info->block_size = 1.0; 677 PetscCall(MatGetInfo(A,MAT_LOCAL,info)); 678 679 isend[0] = info->nz_used; isend[1] = info->nz_allocated; isend[2] = info->nz_unneeded; 680 isend[3] = info->memory; isend[4] = info->mallocs; 681 682 PetscCall(MatGetInfo(B,MAT_LOCAL,info)); 683 684 isend[0] += info->nz_used; isend[1] += info->nz_allocated; isend[2] += info->nz_unneeded; 685 isend[3] += info->memory; isend[4] += info->mallocs; 686 if (flag == MAT_LOCAL) { 687 info->nz_used = isend[0]; 688 info->nz_allocated = isend[1]; 689 info->nz_unneeded = isend[2]; 690 info->memory = isend[3]; 691 info->mallocs = isend[4]; 692 } else if (flag == MAT_GLOBAL_MAX) { 693 PetscCall(MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_MAX,PetscObjectComm((PetscObject)matin))); 694 695 info->nz_used = irecv[0]; 696 info->nz_allocated = irecv[1]; 697 info->nz_unneeded = irecv[2]; 698 info->memory = irecv[3]; 699 info->mallocs = irecv[4]; 700 } else if (flag == MAT_GLOBAL_SUM) { 701 PetscCall(MPIU_Allreduce(isend,irecv,5,MPIU_PETSCLOGDOUBLE,MPI_SUM,PetscObjectComm((PetscObject)matin))); 702 703 info->nz_used = irecv[0]; 704 info->nz_allocated = irecv[1]; 705 info->nz_unneeded = irecv[2]; 706 info->memory = irecv[3]; 707 info->mallocs = irecv[4]; 708 } 709 info->fill_ratio_given = 0; /* no parallel LU/ILU/Cholesky */ 710 info->fill_ratio_needed = 0; 711 info->factor_mallocs = 0; 712 PetscFunctionReturn(0); 713 } 714 715 PetscErrorCode MatSetOption_MPISELL(Mat A,MatOption op,PetscBool flg) 716 { 717 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 718 719 PetscFunctionBegin; 720 switch (op) { 721 case MAT_NEW_NONZERO_LOCATIONS: 722 case MAT_NEW_NONZERO_ALLOCATION_ERR: 723 case MAT_UNUSED_NONZERO_LOCATION_ERR: 724 case MAT_KEEP_NONZERO_PATTERN: 725 case MAT_NEW_NONZERO_LOCATION_ERR: 726 case MAT_USE_INODES: 727 case MAT_IGNORE_ZERO_ENTRIES: 728 MatCheckPreallocated(A,1); 729 PetscCall(MatSetOption(a->A,op,flg)); 730 PetscCall(MatSetOption(a->B,op,flg)); 731 break; 732 case MAT_ROW_ORIENTED: 733 MatCheckPreallocated(A,1); 734 a->roworiented = flg; 735 736 PetscCall(MatSetOption(a->A,op,flg)); 737 PetscCall(MatSetOption(a->B,op,flg)); 738 break; 739 case MAT_FORCE_DIAGONAL_ENTRIES: 740 case MAT_SORTED_FULL: 741 PetscCall(PetscInfo(A,"Option %s ignored\n",MatOptions[op])); 742 break; 743 case MAT_IGNORE_OFF_PROC_ENTRIES: 744 a->donotstash = flg; 745 break; 746 case MAT_SPD: 747 A->spd_set = PETSC_TRUE; 748 A->spd = flg; 749 if (flg) { 750 A->symmetric = PETSC_TRUE; 751 A->structurally_symmetric = PETSC_TRUE; 752 A->symmetric_set = PETSC_TRUE; 753 A->structurally_symmetric_set = PETSC_TRUE; 754 } 755 break; 756 case MAT_SYMMETRIC: 757 MatCheckPreallocated(A,1); 758 PetscCall(MatSetOption(a->A,op,flg)); 759 break; 760 case MAT_STRUCTURALLY_SYMMETRIC: 761 MatCheckPreallocated(A,1); 762 PetscCall(MatSetOption(a->A,op,flg)); 763 break; 764 case MAT_HERMITIAN: 765 MatCheckPreallocated(A,1); 766 PetscCall(MatSetOption(a->A,op,flg)); 767 break; 768 case MAT_SYMMETRY_ETERNAL: 769 MatCheckPreallocated(A,1); 770 PetscCall(MatSetOption(a->A,op,flg)); 771 break; 772 default: 773 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 774 } 775 PetscFunctionReturn(0); 776 } 777 778 PetscErrorCode MatDiagonalScale_MPISELL(Mat mat,Vec ll,Vec rr) 779 { 780 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 781 Mat a=sell->A,b=sell->B; 782 PetscInt s1,s2,s3; 783 784 PetscFunctionBegin; 785 PetscCall(MatGetLocalSize(mat,&s2,&s3)); 786 if (rr) { 787 PetscCall(VecGetLocalSize(rr,&s1)); 788 PetscCheckFalse(s1!=s3,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"right vector non-conforming local size"); 789 /* Overlap communication with computation. */ 790 PetscCall(VecScatterBegin(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD)); 791 } 792 if (ll) { 793 PetscCall(VecGetLocalSize(ll,&s1)); 794 PetscCheckFalse(s1!=s2,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"left vector non-conforming local size"); 795 PetscCall((*b->ops->diagonalscale)(b,ll,NULL)); 796 } 797 /* scale the diagonal block */ 798 PetscCall((*a->ops->diagonalscale)(a,ll,rr)); 799 800 if (rr) { 801 /* Do a scatter end and then right scale the off-diagonal block */ 802 PetscCall(VecScatterEnd(sell->Mvctx,rr,sell->lvec,INSERT_VALUES,SCATTER_FORWARD)); 803 PetscCall((*b->ops->diagonalscale)(b,NULL,sell->lvec)); 804 } 805 PetscFunctionReturn(0); 806 } 807 808 PetscErrorCode MatSetUnfactored_MPISELL(Mat A) 809 { 810 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 811 812 PetscFunctionBegin; 813 PetscCall(MatSetUnfactored(a->A)); 814 PetscFunctionReturn(0); 815 } 816 817 PetscErrorCode MatEqual_MPISELL(Mat A,Mat B,PetscBool *flag) 818 { 819 Mat_MPISELL *matB=(Mat_MPISELL*)B->data,*matA=(Mat_MPISELL*)A->data; 820 Mat a,b,c,d; 821 PetscBool flg; 822 823 PetscFunctionBegin; 824 a = matA->A; b = matA->B; 825 c = matB->A; d = matB->B; 826 827 PetscCall(MatEqual(a,c,&flg)); 828 if (flg) { 829 PetscCall(MatEqual(b,d,&flg)); 830 } 831 PetscCall(MPIU_Allreduce(&flg,flag,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A))); 832 PetscFunctionReturn(0); 833 } 834 835 PetscErrorCode MatCopy_MPISELL(Mat A,Mat B,MatStructure str) 836 { 837 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 838 Mat_MPISELL *b=(Mat_MPISELL*)B->data; 839 840 PetscFunctionBegin; 841 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 842 if ((str != SAME_NONZERO_PATTERN) || (A->ops->copy != B->ops->copy)) { 843 /* because of the column compression in the off-processor part of the matrix a->B, 844 the number of columns in a->B and b->B may be different, hence we cannot call 845 the MatCopy() directly on the two parts. If need be, we can provide a more 846 efficient copy than the MatCopy_Basic() by first uncompressing the a->B matrices 847 then copying the submatrices */ 848 PetscCall(MatCopy_Basic(A,B,str)); 849 } else { 850 PetscCall(MatCopy(a->A,b->A,str)); 851 PetscCall(MatCopy(a->B,b->B,str)); 852 } 853 PetscFunctionReturn(0); 854 } 855 856 PetscErrorCode MatSetUp_MPISELL(Mat A) 857 { 858 PetscFunctionBegin; 859 PetscCall(MatMPISELLSetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL)); 860 PetscFunctionReturn(0); 861 } 862 863 extern PetscErrorCode MatConjugate_SeqSELL(Mat); 864 865 PetscErrorCode MatConjugate_MPISELL(Mat mat) 866 { 867 PetscFunctionBegin; 868 if (PetscDefined(USE_COMPLEX)) { 869 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 870 871 PetscCall(MatConjugate_SeqSELL(sell->A)); 872 PetscCall(MatConjugate_SeqSELL(sell->B)); 873 } 874 PetscFunctionReturn(0); 875 } 876 877 PetscErrorCode MatRealPart_MPISELL(Mat A) 878 { 879 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 880 881 PetscFunctionBegin; 882 PetscCall(MatRealPart(a->A)); 883 PetscCall(MatRealPart(a->B)); 884 PetscFunctionReturn(0); 885 } 886 887 PetscErrorCode MatImaginaryPart_MPISELL(Mat A) 888 { 889 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 890 891 PetscFunctionBegin; 892 PetscCall(MatImaginaryPart(a->A)); 893 PetscCall(MatImaginaryPart(a->B)); 894 PetscFunctionReturn(0); 895 } 896 897 PetscErrorCode MatInvertBlockDiagonal_MPISELL(Mat A,const PetscScalar **values) 898 { 899 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 900 901 PetscFunctionBegin; 902 PetscCall(MatInvertBlockDiagonal(a->A,values)); 903 A->factorerrortype = a->A->factorerrortype; 904 PetscFunctionReturn(0); 905 } 906 907 static PetscErrorCode MatSetRandom_MPISELL(Mat x,PetscRandom rctx) 908 { 909 Mat_MPISELL *sell=(Mat_MPISELL*)x->data; 910 911 PetscFunctionBegin; 912 PetscCall(MatSetRandom(sell->A,rctx)); 913 PetscCall(MatSetRandom(sell->B,rctx)); 914 PetscCall(MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY)); 915 PetscCall(MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY)); 916 PetscFunctionReturn(0); 917 } 918 919 PetscErrorCode MatSetFromOptions_MPISELL(PetscOptionItems *PetscOptionsObject,Mat A) 920 { 921 PetscFunctionBegin; 922 PetscCall(PetscOptionsHead(PetscOptionsObject,"MPISELL options")); 923 PetscCall(PetscOptionsTail()); 924 PetscFunctionReturn(0); 925 } 926 927 PetscErrorCode MatShift_MPISELL(Mat Y,PetscScalar a) 928 { 929 Mat_MPISELL *msell=(Mat_MPISELL*)Y->data; 930 Mat_SeqSELL *sell=(Mat_SeqSELL*)msell->A->data; 931 932 PetscFunctionBegin; 933 if (!Y->preallocated) { 934 PetscCall(MatMPISELLSetPreallocation(Y,1,NULL,0,NULL)); 935 } else if (!sell->nz) { 936 PetscInt nonew = sell->nonew; 937 PetscCall(MatSeqSELLSetPreallocation(msell->A,1,NULL)); 938 sell->nonew = nonew; 939 } 940 PetscCall(MatShift_Basic(Y,a)); 941 PetscFunctionReturn(0); 942 } 943 944 PetscErrorCode MatMissingDiagonal_MPISELL(Mat A,PetscBool *missing,PetscInt *d) 945 { 946 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 947 948 PetscFunctionBegin; 949 PetscCheckFalse(A->rmap->n != A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only works for square matrices"); 950 PetscCall(MatMissingDiagonal(a->A,missing,d)); 951 if (d) { 952 PetscInt rstart; 953 PetscCall(MatGetOwnershipRange(A,&rstart,NULL)); 954 *d += rstart; 955 956 } 957 PetscFunctionReturn(0); 958 } 959 960 PetscErrorCode MatGetDiagonalBlock_MPISELL(Mat A,Mat *a) 961 { 962 PetscFunctionBegin; 963 *a = ((Mat_MPISELL*)A->data)->A; 964 PetscFunctionReturn(0); 965 } 966 967 /* -------------------------------------------------------------------*/ 968 static struct _MatOps MatOps_Values = {MatSetValues_MPISELL, 969 NULL, 970 NULL, 971 MatMult_MPISELL, 972 /* 4*/ MatMultAdd_MPISELL, 973 MatMultTranspose_MPISELL, 974 MatMultTransposeAdd_MPISELL, 975 NULL, 976 NULL, 977 NULL, 978 /*10*/ NULL, 979 NULL, 980 NULL, 981 MatSOR_MPISELL, 982 NULL, 983 /*15*/ MatGetInfo_MPISELL, 984 MatEqual_MPISELL, 985 MatGetDiagonal_MPISELL, 986 MatDiagonalScale_MPISELL, 987 NULL, 988 /*20*/ MatAssemblyBegin_MPISELL, 989 MatAssemblyEnd_MPISELL, 990 MatSetOption_MPISELL, 991 MatZeroEntries_MPISELL, 992 /*24*/ NULL, 993 NULL, 994 NULL, 995 NULL, 996 NULL, 997 /*29*/ MatSetUp_MPISELL, 998 NULL, 999 NULL, 1000 MatGetDiagonalBlock_MPISELL, 1001 NULL, 1002 /*34*/ MatDuplicate_MPISELL, 1003 NULL, 1004 NULL, 1005 NULL, 1006 NULL, 1007 /*39*/ NULL, 1008 NULL, 1009 NULL, 1010 MatGetValues_MPISELL, 1011 MatCopy_MPISELL, 1012 /*44*/ NULL, 1013 MatScale_MPISELL, 1014 MatShift_MPISELL, 1015 MatDiagonalSet_MPISELL, 1016 NULL, 1017 /*49*/ MatSetRandom_MPISELL, 1018 NULL, 1019 NULL, 1020 NULL, 1021 NULL, 1022 /*54*/ MatFDColoringCreate_MPIXAIJ, 1023 NULL, 1024 MatSetUnfactored_MPISELL, 1025 NULL, 1026 NULL, 1027 /*59*/ NULL, 1028 MatDestroy_MPISELL, 1029 MatView_MPISELL, 1030 NULL, 1031 NULL, 1032 /*64*/ NULL, 1033 NULL, 1034 NULL, 1035 NULL, 1036 NULL, 1037 /*69*/ NULL, 1038 NULL, 1039 NULL, 1040 NULL, 1041 NULL, 1042 NULL, 1043 /*75*/ MatFDColoringApply_AIJ, /* reuse AIJ function */ 1044 MatSetFromOptions_MPISELL, 1045 NULL, 1046 NULL, 1047 NULL, 1048 /*80*/ NULL, 1049 NULL, 1050 NULL, 1051 /*83*/ NULL, 1052 NULL, 1053 NULL, 1054 NULL, 1055 NULL, 1056 NULL, 1057 /*89*/ NULL, 1058 NULL, 1059 NULL, 1060 NULL, 1061 NULL, 1062 /*94*/ NULL, 1063 NULL, 1064 NULL, 1065 NULL, 1066 NULL, 1067 /*99*/ NULL, 1068 NULL, 1069 NULL, 1070 MatConjugate_MPISELL, 1071 NULL, 1072 /*104*/NULL, 1073 MatRealPart_MPISELL, 1074 MatImaginaryPart_MPISELL, 1075 NULL, 1076 NULL, 1077 /*109*/NULL, 1078 NULL, 1079 NULL, 1080 NULL, 1081 MatMissingDiagonal_MPISELL, 1082 /*114*/NULL, 1083 NULL, 1084 MatGetGhosts_MPISELL, 1085 NULL, 1086 NULL, 1087 /*119*/NULL, 1088 NULL, 1089 NULL, 1090 NULL, 1091 NULL, 1092 /*124*/NULL, 1093 NULL, 1094 MatInvertBlockDiagonal_MPISELL, 1095 NULL, 1096 NULL, 1097 /*129*/NULL, 1098 NULL, 1099 NULL, 1100 NULL, 1101 NULL, 1102 /*134*/NULL, 1103 NULL, 1104 NULL, 1105 NULL, 1106 NULL, 1107 /*139*/NULL, 1108 NULL, 1109 NULL, 1110 MatFDColoringSetUp_MPIXAIJ, 1111 NULL, 1112 /*144*/NULL, 1113 NULL, 1114 NULL, 1115 NULL 1116 }; 1117 1118 /* ----------------------------------------------------------------------------------------*/ 1119 1120 PetscErrorCode MatStoreValues_MPISELL(Mat mat) 1121 { 1122 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1123 1124 PetscFunctionBegin; 1125 PetscCall(MatStoreValues(sell->A)); 1126 PetscCall(MatStoreValues(sell->B)); 1127 PetscFunctionReturn(0); 1128 } 1129 1130 PetscErrorCode MatRetrieveValues_MPISELL(Mat mat) 1131 { 1132 Mat_MPISELL *sell=(Mat_MPISELL*)mat->data; 1133 1134 PetscFunctionBegin; 1135 PetscCall(MatRetrieveValues(sell->A)); 1136 PetscCall(MatRetrieveValues(sell->B)); 1137 PetscFunctionReturn(0); 1138 } 1139 1140 PetscErrorCode MatMPISELLSetPreallocation_MPISELL(Mat B,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[]) 1141 { 1142 Mat_MPISELL *b; 1143 1144 PetscFunctionBegin; 1145 PetscCall(PetscLayoutSetUp(B->rmap)); 1146 PetscCall(PetscLayoutSetUp(B->cmap)); 1147 b = (Mat_MPISELL*)B->data; 1148 1149 if (!B->preallocated) { 1150 /* Explicitly create 2 MATSEQSELL matrices. */ 1151 PetscCall(MatCreate(PETSC_COMM_SELF,&b->A)); 1152 PetscCall(MatSetSizes(b->A,B->rmap->n,B->cmap->n,B->rmap->n,B->cmap->n)); 1153 PetscCall(MatSetBlockSizesFromMats(b->A,B,B)); 1154 PetscCall(MatSetType(b->A,MATSEQSELL)); 1155 PetscCall(PetscLogObjectParent((PetscObject)B,(PetscObject)b->A)); 1156 PetscCall(MatCreate(PETSC_COMM_SELF,&b->B)); 1157 PetscCall(MatSetSizes(b->B,B->rmap->n,B->cmap->N,B->rmap->n,B->cmap->N)); 1158 PetscCall(MatSetBlockSizesFromMats(b->B,B,B)); 1159 PetscCall(MatSetType(b->B,MATSEQSELL)); 1160 PetscCall(PetscLogObjectParent((PetscObject)B,(PetscObject)b->B)); 1161 } 1162 1163 PetscCall(MatSeqSELLSetPreallocation(b->A,d_rlenmax,d_rlen)); 1164 PetscCall(MatSeqSELLSetPreallocation(b->B,o_rlenmax,o_rlen)); 1165 B->preallocated = PETSC_TRUE; 1166 B->was_assembled = PETSC_FALSE; 1167 /* 1168 critical for MatAssemblyEnd to work. 1169 MatAssemblyBegin checks it to set up was_assembled 1170 and MatAssemblyEnd checks was_assembled to determine whether to build garray 1171 */ 1172 B->assembled = PETSC_FALSE; 1173 PetscFunctionReturn(0); 1174 } 1175 1176 PetscErrorCode MatDuplicate_MPISELL(Mat matin,MatDuplicateOption cpvalues,Mat *newmat) 1177 { 1178 Mat mat; 1179 Mat_MPISELL *a,*oldmat=(Mat_MPISELL*)matin->data; 1180 1181 PetscFunctionBegin; 1182 *newmat = NULL; 1183 PetscCall(MatCreate(PetscObjectComm((PetscObject)matin),&mat)); 1184 PetscCall(MatSetSizes(mat,matin->rmap->n,matin->cmap->n,matin->rmap->N,matin->cmap->N)); 1185 PetscCall(MatSetBlockSizesFromMats(mat,matin,matin)); 1186 PetscCall(MatSetType(mat,((PetscObject)matin)->type_name)); 1187 a = (Mat_MPISELL*)mat->data; 1188 1189 mat->factortype = matin->factortype; 1190 mat->assembled = PETSC_TRUE; 1191 mat->insertmode = NOT_SET_VALUES; 1192 mat->preallocated = PETSC_TRUE; 1193 1194 a->size = oldmat->size; 1195 a->rank = oldmat->rank; 1196 a->donotstash = oldmat->donotstash; 1197 a->roworiented = oldmat->roworiented; 1198 a->rowindices = NULL; 1199 a->rowvalues = NULL; 1200 a->getrowactive = PETSC_FALSE; 1201 1202 PetscCall(PetscLayoutReference(matin->rmap,&mat->rmap)); 1203 PetscCall(PetscLayoutReference(matin->cmap,&mat->cmap)); 1204 1205 if (oldmat->colmap) { 1206 #if defined(PETSC_USE_CTABLE) 1207 PetscCall(PetscTableCreateCopy(oldmat->colmap,&a->colmap)); 1208 #else 1209 PetscCall(PetscMalloc1(mat->cmap->N,&a->colmap)); 1210 PetscCall(PetscLogObjectMemory((PetscObject)mat,(mat->cmap->N)*sizeof(PetscInt))); 1211 PetscCall(PetscArraycpy(a->colmap,oldmat->colmap,mat->cmap->N)); 1212 #endif 1213 } else a->colmap = NULL; 1214 if (oldmat->garray) { 1215 PetscInt len; 1216 len = oldmat->B->cmap->n; 1217 PetscCall(PetscMalloc1(len+1,&a->garray)); 1218 PetscCall(PetscLogObjectMemory((PetscObject)mat,len*sizeof(PetscInt))); 1219 if (len) PetscCall(PetscArraycpy(a->garray,oldmat->garray,len)); 1220 } else a->garray = NULL; 1221 1222 PetscCall(VecDuplicate(oldmat->lvec,&a->lvec)); 1223 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->lvec)); 1224 PetscCall(VecScatterCopy(oldmat->Mvctx,&a->Mvctx)); 1225 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->Mvctx)); 1226 PetscCall(MatDuplicate(oldmat->A,cpvalues,&a->A)); 1227 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->A)); 1228 PetscCall(MatDuplicate(oldmat->B,cpvalues,&a->B)); 1229 PetscCall(PetscLogObjectParent((PetscObject)mat,(PetscObject)a->B)); 1230 PetscCall(PetscFunctionListDuplicate(((PetscObject)matin)->qlist,&((PetscObject)mat)->qlist)); 1231 *newmat = mat; 1232 PetscFunctionReturn(0); 1233 } 1234 1235 /*@C 1236 MatMPISELLSetPreallocation - Preallocates memory for a sparse parallel matrix in sell format. 1237 For good matrix assembly performance the user should preallocate the matrix storage by 1238 setting the parameters d_nz (or d_nnz) and o_nz (or o_nnz). 1239 1240 Collective 1241 1242 Input Parameters: 1243 + B - the matrix 1244 . d_nz - number of nonzeros per row in DIAGONAL portion of local submatrix 1245 (same value is used for all local rows) 1246 . d_nnz - array containing the number of nonzeros in the various rows of the 1247 DIAGONAL portion of the local submatrix (possibly different for each row) 1248 or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure. 1249 The size of this array is equal to the number of local rows, i.e 'm'. 1250 For matrices that will be factored, you must leave room for (and set) 1251 the diagonal entry even if it is zero. 1252 . o_nz - number of nonzeros per row in the OFF-DIAGONAL portion of local 1253 submatrix (same value is used for all local rows). 1254 - o_nnz - array containing the number of nonzeros in the various rows of the 1255 OFF-DIAGONAL portion of the local submatrix (possibly different for 1256 each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero 1257 structure. The size of this array is equal to the number 1258 of local rows, i.e 'm'. 1259 1260 If the *_nnz parameter is given then the *_nz parameter is ignored 1261 1262 The stored row and column indices begin with zero. 1263 1264 The parallel matrix is partitioned such that the first m0 rows belong to 1265 process 0, the next m1 rows belong to process 1, the next m2 rows belong 1266 to process 2 etc.. where m0,m1,m2... are the input parameter 'm'. 1267 1268 The DIAGONAL portion of the local submatrix of a processor can be defined 1269 as the submatrix which is obtained by extraction the part corresponding to 1270 the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the 1271 first row that belongs to the processor, r2 is the last row belonging to 1272 the this processor, and c1-c2 is range of indices of the local part of a 1273 vector suitable for applying the matrix to. This is an mxn matrix. In the 1274 common case of a square matrix, the row and column ranges are the same and 1275 the DIAGONAL part is also square. The remaining portion of the local 1276 submatrix (mxN) constitute the OFF-DIAGONAL portion. 1277 1278 If o_nnz, d_nnz are specified, then o_nz, and d_nz are ignored. 1279 1280 You can call MatGetInfo() to get information on how effective the preallocation was; 1281 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1282 You can also run with the option -info and look for messages with the string 1283 malloc in them to see if additional memory allocation was needed. 1284 1285 Example usage: 1286 1287 Consider the following 8x8 matrix with 34 non-zero values, that is 1288 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1289 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1290 as follows: 1291 1292 .vb 1293 1 2 0 | 0 3 0 | 0 4 1294 Proc0 0 5 6 | 7 0 0 | 8 0 1295 9 0 10 | 11 0 0 | 12 0 1296 ------------------------------------- 1297 13 0 14 | 15 16 17 | 0 0 1298 Proc1 0 18 0 | 19 20 21 | 0 0 1299 0 0 0 | 22 23 0 | 24 0 1300 ------------------------------------- 1301 Proc2 25 26 27 | 0 0 28 | 29 0 1302 30 0 0 | 31 32 33 | 0 34 1303 .ve 1304 1305 This can be represented as a collection of submatrices as: 1306 1307 .vb 1308 A B C 1309 D E F 1310 G H I 1311 .ve 1312 1313 Where the submatrices A,B,C are owned by proc0, D,E,F are 1314 owned by proc1, G,H,I are owned by proc2. 1315 1316 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1317 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1318 The 'M','N' parameters are 8,8, and have the same values on all procs. 1319 1320 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1321 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1322 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1323 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1324 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1325 matrix, ans [DF] as another SeqSELL matrix. 1326 1327 When d_nz, o_nz parameters are specified, d_nz storage elements are 1328 allocated for every row of the local diagonal submatrix, and o_nz 1329 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1330 One way to choose d_nz and o_nz is to use the max nonzerors per local 1331 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1332 In this case, the values of d_nz,o_nz are: 1333 .vb 1334 proc0 : dnz = 2, o_nz = 2 1335 proc1 : dnz = 3, o_nz = 2 1336 proc2 : dnz = 1, o_nz = 4 1337 .ve 1338 We are allocating m*(d_nz+o_nz) storage locations for every proc. This 1339 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1340 for proc3. i.e we are using 12+15+10=37 storage locations to store 1341 34 values. 1342 1343 When d_nnz, o_nnz parameters are specified, the storage is specified 1344 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1345 In the above case the values for d_nnz,o_nnz are: 1346 .vb 1347 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1348 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1349 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1350 .ve 1351 Here the space allocated is according to nz (or maximum values in the nnz 1352 if nnz is provided) for DIAGONAL and OFF-DIAGONAL submatrices, i.e (2+2+3+2)*3+(1+4)*2=37 1353 1354 Level: intermediate 1355 1356 .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatCreatesell(), 1357 MATMPISELL, MatGetInfo(), PetscSplitOwnership() 1358 @*/ 1359 PetscErrorCode MatMPISELLSetPreallocation(Mat B,PetscInt d_nz,const PetscInt d_nnz[],PetscInt o_nz,const PetscInt o_nnz[]) 1360 { 1361 PetscFunctionBegin; 1362 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 1363 PetscValidType(B,1); 1364 PetscTryMethod(B,"MatMPISELLSetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(B,d_nz,d_nnz,o_nz,o_nnz)); 1365 PetscFunctionReturn(0); 1366 } 1367 1368 /*MC 1369 MATMPISELL - MATMPISELL = "mpisell" - A matrix type to be used for MPI sparse matrices, 1370 based on the sliced Ellpack format 1371 1372 Options Database Keys: 1373 . -mat_type sell - sets the matrix type to "seqsell" during a call to MatSetFromOptions() 1374 1375 Level: beginner 1376 1377 .seealso: MatCreateSell(), MATSEQSELL, MATSELL, MATSEQAIJ, MATAIJ, MATMPIAIJ 1378 M*/ 1379 1380 /*@C 1381 MatCreateSELL - Creates a sparse parallel matrix in SELL format. 1382 1383 Collective 1384 1385 Input Parameters: 1386 + comm - MPI communicator 1387 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 1388 This value should be the same as the local size used in creating the 1389 y vector for the matrix-vector product y = Ax. 1390 . n - This value should be the same as the local size used in creating the 1391 x vector for the matrix-vector product y = Ax. (or PETSC_DECIDE to have 1392 calculated if N is given) For square matrices n is almost always m. 1393 . M - number of global rows (or PETSC_DETERMINE to have calculated if m is given) 1394 . N - number of global columns (or PETSC_DETERMINE to have calculated if n is given) 1395 . d_rlenmax - max number of nonzeros per row in DIAGONAL portion of local submatrix 1396 (same value is used for all local rows) 1397 . d_rlen - array containing the number of nonzeros in the various rows of the 1398 DIAGONAL portion of the local submatrix (possibly different for each row) 1399 or NULL, if d_rlenmax is used to specify the nonzero structure. 1400 The size of this array is equal to the number of local rows, i.e 'm'. 1401 . o_rlenmax - max number of nonzeros per row in the OFF-DIAGONAL portion of local 1402 submatrix (same value is used for all local rows). 1403 - o_rlen - array containing the number of nonzeros in the various rows of the 1404 OFF-DIAGONAL portion of the local submatrix (possibly different for 1405 each row) or NULL, if o_rlenmax is used to specify the nonzero 1406 structure. The size of this array is equal to the number 1407 of local rows, i.e 'm'. 1408 1409 Output Parameter: 1410 . A - the matrix 1411 1412 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1413 MatXXXXSetPreallocation() paradigm instead of this routine directly. 1414 [MatXXXXSetPreallocation() is, for example, MatSeqSELLSetPreallocation] 1415 1416 Notes: 1417 If the *_rlen parameter is given then the *_rlenmax parameter is ignored 1418 1419 m,n,M,N parameters specify the size of the matrix, and its partitioning across 1420 processors, while d_rlenmax,d_rlen,o_rlenmax,o_rlen parameters specify the approximate 1421 storage requirements for this matrix. 1422 1423 If PETSC_DECIDE or PETSC_DETERMINE is used for a particular argument on one 1424 processor than it must be used on all processors that share the object for 1425 that argument. 1426 1427 The user MUST specify either the local or global matrix dimensions 1428 (possibly both). 1429 1430 The parallel matrix is partitioned across processors such that the 1431 first m0 rows belong to process 0, the next m1 rows belong to 1432 process 1, the next m2 rows belong to process 2 etc.. where 1433 m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores 1434 values corresponding to [m x N] submatrix. 1435 1436 The columns are logically partitioned with the n0 columns belonging 1437 to 0th partition, the next n1 columns belonging to the next 1438 partition etc.. where n0,n1,n2... are the input parameter 'n'. 1439 1440 The DIAGONAL portion of the local submatrix on any given processor 1441 is the submatrix corresponding to the rows and columns m,n 1442 corresponding to the given processor. i.e diagonal matrix on 1443 process 0 is [m0 x n0], diagonal matrix on process 1 is [m1 x n1] 1444 etc. The remaining portion of the local submatrix [m x (N-n)] 1445 constitute the OFF-DIAGONAL portion. The example below better 1446 illustrates this concept. 1447 1448 For a square global matrix we define each processor's diagonal portion 1449 to be its local rows and the corresponding columns (a square submatrix); 1450 each processor's off-diagonal portion encompasses the remainder of the 1451 local matrix (a rectangular submatrix). 1452 1453 If o_rlen, d_rlen are specified, then o_rlenmax, and d_rlenmax are ignored. 1454 1455 When calling this routine with a single process communicator, a matrix of 1456 type SEQSELL is returned. If a matrix of type MATMPISELL is desired for this 1457 type of communicator, use the construction mechanism: 1458 MatCreate(...,&A); MatSetType(A,MATMPISELL); MatSetSizes(A, m,n,M,N); MatMPISELLSetPreallocation(A,...); 1459 1460 Options Database Keys: 1461 - -mat_sell_oneindex - Internally use indexing starting at 1 1462 rather than 0. Note that when calling MatSetValues(), 1463 the user still MUST index entries starting at 0! 1464 1465 Example usage: 1466 1467 Consider the following 8x8 matrix with 34 non-zero values, that is 1468 assembled across 3 processors. Lets assume that proc0 owns 3 rows, 1469 proc1 owns 3 rows, proc2 owns 2 rows. This division can be shown 1470 as follows: 1471 1472 .vb 1473 1 2 0 | 0 3 0 | 0 4 1474 Proc0 0 5 6 | 7 0 0 | 8 0 1475 9 0 10 | 11 0 0 | 12 0 1476 ------------------------------------- 1477 13 0 14 | 15 16 17 | 0 0 1478 Proc1 0 18 0 | 19 20 21 | 0 0 1479 0 0 0 | 22 23 0 | 24 0 1480 ------------------------------------- 1481 Proc2 25 26 27 | 0 0 28 | 29 0 1482 30 0 0 | 31 32 33 | 0 34 1483 .ve 1484 1485 This can be represented as a collection of submatrices as: 1486 1487 .vb 1488 A B C 1489 D E F 1490 G H I 1491 .ve 1492 1493 Where the submatrices A,B,C are owned by proc0, D,E,F are 1494 owned by proc1, G,H,I are owned by proc2. 1495 1496 The 'm' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1497 The 'n' parameters for proc0,proc1,proc2 are 3,3,2 respectively. 1498 The 'M','N' parameters are 8,8, and have the same values on all procs. 1499 1500 The DIAGONAL submatrices corresponding to proc0,proc1,proc2 are 1501 submatrices [A], [E], [I] respectively. The OFF-DIAGONAL submatrices 1502 corresponding to proc0,proc1,proc2 are [BC], [DF], [GH] respectively. 1503 Internally, each processor stores the DIAGONAL part, and the OFF-DIAGONAL 1504 part as SeqSELL matrices. for eg: proc1 will store [E] as a SeqSELL 1505 matrix, ans [DF] as another SeqSELL matrix. 1506 1507 When d_rlenmax, o_rlenmax parameters are specified, d_rlenmax storage elements are 1508 allocated for every row of the local diagonal submatrix, and o_rlenmax 1509 storage locations are allocated for every row of the OFF-DIAGONAL submat. 1510 One way to choose d_rlenmax and o_rlenmax is to use the max nonzerors per local 1511 rows for each of the local DIAGONAL, and the OFF-DIAGONAL submatrices. 1512 In this case, the values of d_rlenmax,o_rlenmax are: 1513 .vb 1514 proc0 : d_rlenmax = 2, o_rlenmax = 2 1515 proc1 : d_rlenmax = 3, o_rlenmax = 2 1516 proc2 : d_rlenmax = 1, o_rlenmax = 4 1517 .ve 1518 We are allocating m*(d_rlenmax+o_rlenmax) storage locations for every proc. This 1519 translates to 3*(2+2)=12 for proc0, 3*(3+2)=15 for proc1, 2*(1+4)=10 1520 for proc3. i.e we are using 12+15+10=37 storage locations to store 1521 34 values. 1522 1523 When d_rlen, o_rlen parameters are specified, the storage is specified 1524 for every row, corresponding to both DIAGONAL and OFF-DIAGONAL submatrices. 1525 In the above case the values for d_nnz,o_nnz are: 1526 .vb 1527 proc0: d_nnz = [2,2,2] and o_nnz = [2,2,2] 1528 proc1: d_nnz = [3,3,2] and o_nnz = [2,1,1] 1529 proc2: d_nnz = [1,1] and o_nnz = [4,4] 1530 .ve 1531 Here the space allocated is still 37 though there are 34 nonzeros because 1532 the allocation is always done according to rlenmax. 1533 1534 Level: intermediate 1535 1536 .seealso: MatCreate(), MatCreateSeqSELL(), MatSetValues(), MatMPISELLSetPreallocation(), MatMPISELLSetPreallocationSELL(), 1537 MATMPISELL, MatCreateMPISELLWithArrays() 1538 @*/ 1539 PetscErrorCode MatCreateSELL(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt d_rlenmax,const PetscInt d_rlen[],PetscInt o_rlenmax,const PetscInt o_rlen[],Mat *A) 1540 { 1541 PetscMPIInt size; 1542 1543 PetscFunctionBegin; 1544 PetscCall(MatCreate(comm,A)); 1545 PetscCall(MatSetSizes(*A,m,n,M,N)); 1546 PetscCallMPI(MPI_Comm_size(comm,&size)); 1547 if (size > 1) { 1548 PetscCall(MatSetType(*A,MATMPISELL)); 1549 PetscCall(MatMPISELLSetPreallocation(*A,d_rlenmax,d_rlen,o_rlenmax,o_rlen)); 1550 } else { 1551 PetscCall(MatSetType(*A,MATSEQSELL)); 1552 PetscCall(MatSeqSELLSetPreallocation(*A,d_rlenmax,d_rlen)); 1553 } 1554 PetscFunctionReturn(0); 1555 } 1556 1557 PetscErrorCode MatMPISELLGetSeqSELL(Mat A,Mat *Ad,Mat *Ao,const PetscInt *colmap[]) 1558 { 1559 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1560 PetscBool flg; 1561 1562 PetscFunctionBegin; 1563 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&flg)); 1564 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"This function requires a MATMPISELL matrix as input"); 1565 if (Ad) *Ad = a->A; 1566 if (Ao) *Ao = a->B; 1567 if (colmap) *colmap = a->garray; 1568 PetscFunctionReturn(0); 1569 } 1570 1571 /*@C 1572 MatMPISELLGetLocalMatCondensed - Creates a SeqSELL matrix from an MATMPISELL matrix by taking all its local rows and NON-ZERO columns 1573 1574 Not Collective 1575 1576 Input Parameters: 1577 + A - the matrix 1578 . scall - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 1579 - row, col - index sets of rows and columns to extract (or NULL) 1580 1581 Output Parameter: 1582 . A_loc - the local sequential matrix generated 1583 1584 Level: developer 1585 1586 .seealso: MatGetOwnershipRange(), MatMPISELLGetLocalMat() 1587 1588 @*/ 1589 PetscErrorCode MatMPISELLGetLocalMatCondensed(Mat A,MatReuse scall,IS *row,IS *col,Mat *A_loc) 1590 { 1591 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1592 PetscInt i,start,end,ncols,nzA,nzB,*cmap,imark,*idx; 1593 IS isrowa,iscola; 1594 Mat *aloc; 1595 PetscBool match; 1596 1597 PetscFunctionBegin; 1598 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATMPISELL,&match)); 1599 PetscCheck(match,PetscObjectComm((PetscObject)A), PETSC_ERR_SUP,"Requires MATMPISELL matrix as input"); 1600 PetscCall(PetscLogEventBegin(MAT_Getlocalmatcondensed,A,0,0,0)); 1601 if (!row) { 1602 start = A->rmap->rstart; end = A->rmap->rend; 1603 PetscCall(ISCreateStride(PETSC_COMM_SELF,end-start,start,1,&isrowa)); 1604 } else { 1605 isrowa = *row; 1606 } 1607 if (!col) { 1608 start = A->cmap->rstart; 1609 cmap = a->garray; 1610 nzA = a->A->cmap->n; 1611 nzB = a->B->cmap->n; 1612 PetscCall(PetscMalloc1(nzA+nzB, &idx)); 1613 ncols = 0; 1614 for (i=0; i<nzB; i++) { 1615 if (cmap[i] < start) idx[ncols++] = cmap[i]; 1616 else break; 1617 } 1618 imark = i; 1619 for (i=0; i<nzA; i++) idx[ncols++] = start + i; 1620 for (i=imark; i<nzB; i++) idx[ncols++] = cmap[i]; 1621 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,ncols,idx,PETSC_OWN_POINTER,&iscola)); 1622 } else { 1623 iscola = *col; 1624 } 1625 if (scall != MAT_INITIAL_MATRIX) { 1626 PetscCall(PetscMalloc1(1,&aloc)); 1627 aloc[0] = *A_loc; 1628 } 1629 PetscCall(MatCreateSubMatrices(A,1,&isrowa,&iscola,scall,&aloc)); 1630 *A_loc = aloc[0]; 1631 PetscCall(PetscFree(aloc)); 1632 if (!row) { 1633 PetscCall(ISDestroy(&isrowa)); 1634 } 1635 if (!col) { 1636 PetscCall(ISDestroy(&iscola)); 1637 } 1638 PetscCall(PetscLogEventEnd(MAT_Getlocalmatcondensed,A,0,0,0)); 1639 PetscFunctionReturn(0); 1640 } 1641 1642 #include <../src/mat/impls/aij/mpi/mpiaij.h> 1643 1644 PetscErrorCode MatConvert_MPISELL_MPIAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1645 { 1646 Mat_MPISELL *a=(Mat_MPISELL*)A->data; 1647 Mat B; 1648 Mat_MPIAIJ *b; 1649 1650 PetscFunctionBegin; 1651 PetscCheck(A->assembled,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1652 1653 if (reuse == MAT_REUSE_MATRIX) { 1654 B = *newmat; 1655 } else { 1656 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 1657 PetscCall(MatSetType(B,MATMPIAIJ)); 1658 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 1659 PetscCall(MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs)); 1660 PetscCall(MatSeqAIJSetPreallocation(B,0,NULL)); 1661 PetscCall(MatMPIAIJSetPreallocation(B,0,NULL,0,NULL)); 1662 } 1663 b = (Mat_MPIAIJ*) B->data; 1664 1665 if (reuse == MAT_REUSE_MATRIX) { 1666 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_REUSE_MATRIX, &b->A)); 1667 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_REUSE_MATRIX, &b->B)); 1668 } else { 1669 PetscCall(MatDestroy(&b->A)); 1670 PetscCall(MatDestroy(&b->B)); 1671 PetscCall(MatDisAssemble_MPISELL(A)); 1672 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->A, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->A)); 1673 PetscCall(MatConvert_SeqSELL_SeqAIJ(a->B, MATSEQAIJ, MAT_INITIAL_MATRIX, &b->B)); 1674 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1675 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1676 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1677 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1678 } 1679 1680 if (reuse == MAT_INPLACE_MATRIX) { 1681 PetscCall(MatHeaderReplace(A,&B)); 1682 } else { 1683 *newmat = B; 1684 } 1685 PetscFunctionReturn(0); 1686 } 1687 1688 PetscErrorCode MatConvert_MPIAIJ_MPISELL(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 1689 { 1690 Mat_MPIAIJ *a=(Mat_MPIAIJ*)A->data; 1691 Mat B; 1692 Mat_MPISELL *b; 1693 1694 PetscFunctionBegin; 1695 PetscCheck(A->assembled,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Matrix must be assembled"); 1696 1697 if (reuse == MAT_REUSE_MATRIX) { 1698 B = *newmat; 1699 } else { 1700 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 1701 PetscCall(MatSetType(B,MATMPISELL)); 1702 PetscCall(MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N)); 1703 PetscCall(MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs)); 1704 PetscCall(MatSeqAIJSetPreallocation(B,0,NULL)); 1705 PetscCall(MatMPIAIJSetPreallocation(B,0,NULL,0,NULL)); 1706 } 1707 b = (Mat_MPISELL*) B->data; 1708 1709 if (reuse == MAT_REUSE_MATRIX) { 1710 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_REUSE_MATRIX, &b->A)); 1711 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_REUSE_MATRIX, &b->B)); 1712 } else { 1713 PetscCall(MatDestroy(&b->A)); 1714 PetscCall(MatDestroy(&b->B)); 1715 PetscCall(MatDisAssemble_MPIAIJ(A)); 1716 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->A, MATSEQSELL, MAT_INITIAL_MATRIX, &b->A)); 1717 PetscCall(MatConvert_SeqAIJ_SeqSELL(a->B, MATSEQSELL, MAT_INITIAL_MATRIX, &b->B)); 1718 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1719 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 1720 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1721 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1722 } 1723 1724 if (reuse == MAT_INPLACE_MATRIX) { 1725 PetscCall(MatHeaderReplace(A,&B)); 1726 } else { 1727 *newmat = B; 1728 } 1729 PetscFunctionReturn(0); 1730 } 1731 1732 PetscErrorCode MatSOR_MPISELL(Mat matin,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1733 { 1734 Mat_MPISELL *mat=(Mat_MPISELL*)matin->data; 1735 Vec bb1=NULL; 1736 1737 PetscFunctionBegin; 1738 if (flag == SOR_APPLY_UPPER) { 1739 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1740 PetscFunctionReturn(0); 1741 } 1742 1743 if (its > 1 || ~flag & SOR_ZERO_INITIAL_GUESS || flag & SOR_EISENSTAT) { 1744 PetscCall(VecDuplicate(bb,&bb1)); 1745 } 1746 1747 if ((flag & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) { 1748 if (flag & SOR_ZERO_INITIAL_GUESS) { 1749 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1750 its--; 1751 } 1752 1753 while (its--) { 1754 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1755 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1756 1757 /* update rhs: bb1 = bb - B*x */ 1758 PetscCall(VecScale(mat->lvec,-1.0)); 1759 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1760 1761 /* local sweep */ 1762 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_SYMMETRIC_SWEEP,fshift,lits,1,xx)); 1763 } 1764 } else if (flag & SOR_LOCAL_FORWARD_SWEEP) { 1765 if (flag & SOR_ZERO_INITIAL_GUESS) { 1766 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1767 its--; 1768 } 1769 while (its--) { 1770 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1771 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1772 1773 /* update rhs: bb1 = bb - B*x */ 1774 PetscCall(VecScale(mat->lvec,-1.0)); 1775 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1776 1777 /* local sweep */ 1778 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_FORWARD_SWEEP,fshift,lits,1,xx)); 1779 } 1780 } else if (flag & SOR_LOCAL_BACKWARD_SWEEP) { 1781 if (flag & SOR_ZERO_INITIAL_GUESS) { 1782 PetscCall((*mat->A->ops->sor)(mat->A,bb,omega,flag,fshift,lits,1,xx)); 1783 its--; 1784 } 1785 while (its--) { 1786 PetscCall(VecScatterBegin(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1787 PetscCall(VecScatterEnd(mat->Mvctx,xx,mat->lvec,INSERT_VALUES,SCATTER_FORWARD)); 1788 1789 /* update rhs: bb1 = bb - B*x */ 1790 PetscCall(VecScale(mat->lvec,-1.0)); 1791 PetscCall((*mat->B->ops->multadd)(mat->B,mat->lvec,bb,bb1)); 1792 1793 /* local sweep */ 1794 PetscCall((*mat->A->ops->sor)(mat->A,bb1,omega,SOR_BACKWARD_SWEEP,fshift,lits,1,xx)); 1795 } 1796 } else SETERRQ(PetscObjectComm((PetscObject)matin),PETSC_ERR_SUP,"Parallel SOR not supported"); 1797 1798 PetscCall(VecDestroy(&bb1)); 1799 1800 matin->factorerrortype = mat->A->factorerrortype; 1801 PetscFunctionReturn(0); 1802 } 1803 1804 /*MC 1805 MATMPISELL - MATMPISELL = "MPISELL" - A matrix type to be used for parallel sparse matrices. 1806 1807 Options Database Keys: 1808 . -mat_type MPISELL - sets the matrix type to "MPISELL" during a call to MatSetFromOptions() 1809 1810 Level: beginner 1811 1812 .seealso: MatCreateSELL() 1813 M*/ 1814 PETSC_EXTERN PetscErrorCode MatCreate_MPISELL(Mat B) 1815 { 1816 Mat_MPISELL *b; 1817 PetscMPIInt size; 1818 1819 PetscFunctionBegin; 1820 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 1821 PetscCall(PetscNewLog(B,&b)); 1822 B->data = (void*)b; 1823 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 1824 B->assembled = PETSC_FALSE; 1825 B->insertmode = NOT_SET_VALUES; 1826 b->size = size; 1827 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B),&b->rank)); 1828 /* build cache for off array entries formed */ 1829 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash)); 1830 1831 b->donotstash = PETSC_FALSE; 1832 b->colmap = NULL; 1833 b->garray = NULL; 1834 b->roworiented = PETSC_TRUE; 1835 1836 /* stuff used for matrix vector multiply */ 1837 b->lvec = NULL; 1838 b->Mvctx = NULL; 1839 1840 /* stuff for MatGetRow() */ 1841 b->rowindices = NULL; 1842 b->rowvalues = NULL; 1843 b->getrowactive = PETSC_FALSE; 1844 1845 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_MPISELL)); 1846 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_MPISELL)); 1847 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_MPISELL)); 1848 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatMPISELLSetPreallocation_C",MatMPISELLSetPreallocation_MPISELL)); 1849 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpisell_mpiaij_C",MatConvert_MPISELL_MPIAIJ)); 1850 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatDiagonalScaleLocal_C",MatDiagonalScaleLocal_MPISELL)); 1851 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATMPISELL)); 1852 PetscFunctionReturn(0); 1853 } 1854